In [1]:
# -*- coding: utf-8 -*-
# Copyright 2019 - 2022 United Kingdom Research and Innovation
# Copyright 2019 - 2022 The University of Manchester
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
Primal Dual Hybrid Gradient Algorithm¶
We use the Primal Dual Hybrid Algorithm (PDHG) introduced by Chambolle & Pock for Tomography Reconstruction. We will solve the following minimisation problem under three different regularisation terms, i.e.,
- $\|\cdot\|_{1}$ or
- Tikhonov regularisation or
- with $L=\nabla$ and Total variation:
$$ u^{*} =\underset{u}{\operatorname{argmin}} \frac{1}{2} \| \mathcal{A} u - g\|^{2} + \alpha\,\mathrm{TV}(u) + \mathbb{I}_{\{u\geq 0\}}(u). \tag{all reg} $$
where,
$g$ is the Acquisition data obtained from the detector.
$\mathcal{A}$ is the projection operator ( Radon transform ) that maps from an image-space to an acquisition space, i.e., $\mathcal{A} : \mathbb{X} \rightarrow \mathbb{Y}, $ where $\mathbb{X}$ is an ImageGeometry and $\mathbb{Y}$ is an AcquisitionGeometry.
$\alpha$: regularising parameter that measures the trade-off between the fidelity and the regulariser terms.
The total variation (isotropic) is defined as $$\mathrm{TV}(u) = \|\nabla u \|_{2,1} = \sum \sqrt{ (\partial_{y}u)^{2} + (\partial_{x}u)^{2} }$$
$\mathbb{I}_{{u\geq 0}}(u) : =
$, $\quad$ a non-negativity constraint for the minimiser $u$.
In [2]:
# Import libraries
from cil.framework import BlockDataContainer
from cil.optimisation.functions import L2NormSquared, L1Norm, BlockFunction, MixedL21Norm, IndicatorBox, TotalVariation
from cil.optimisation.operators import GradientOperator, BlockOperator
from cil.optimisation.algorithms import PDHG, SIRT
from cil.plugins.astra.operators import ProjectionOperator
from cil.plugins.astra.processors import FBP
from cil.plugins.ccpi_regularisation.functions import FGP_TV
from cil.utilities.display import show2D, show1D, show_geometry
from cil.utilities.jupyter import islicer
from cil.io import ZEISSDataReader
from cil.processors import Binner, TransmissionAbsorptionConverter, Slicer
import matplotlib.pyplot as plt
import numpy as np
import os
Data information¶
We use the Walnut found in Jørgensen_et_all. In total, there are 6 individual micro Computed Tomography datasets in the native Zeiss TXRM/TXM format. The six datasets were acquired at the 3D Imaging Center at Technical University of Denmark in 2014 (HDTomo3D in 2016) as part of the ERC-funded project High-Definition Tomography (HDTomo) headed by Prof. Per Christian Hansen.
This example requires the dataset walnut.zip from https://zenodo.org/record/4822516 :
If running locally please download the data and update the path variable below.
Load walnut data¶
In [3]:
path = '../../../data/'
In [4]:
reader = ZEISSDataReader()
filename = os.path.join(path, "valnut_tomo-A.txrm")
data3D = ZEISSDataReader(file_name=filename).read()
In [5]:
# reorder data to match default order for Astra/Tigre operator
data3D.reorder('astra')
# Get Image and Acquisition geometries
ag3D = data3D.geometry
ig3D = ag3D.get_ImageGeometry()
Acquisition and Image geometry information¶
In [6]:
print(ag3D)
3D Cone-beam tomography System configuration: Source position: [ 0. , -105.05081177, 0. ] Rotation axis position: [0., 0., 0.] Rotation axis direction: [0., 0., 1.] Detector position: [ 0. , 45.08757401, 0. ] Detector direction x: [1., 0., 0.] Detector direction y: [0., 0., 1.] Panel configuration: Number of pixels: [1024 1024] Pixel size: [0.0658543 0.0658543] Pixel origin: bottom-left Channel configuration: Number of channels: 1 Acquisition description: Number of positions: 1601 Angles 0-20 in radians: [-3.1415665, -3.1377017, -3.1337626, -3.1298182, -3.125836 , -3.1219127, -3.1180956, -3.1140666, -3.1101887, -3.1062822, -3.1022923, -3.0984268, -3.0944946, -3.0905435, -3.0865552, -3.082691 , -3.0787866, -3.074828 , -3.0708766, -3.0669732] Distances in units: units distance
In [7]:
print(ig3D)
Number of channels: 1 channel_spacing: 1.0 voxel_num : x1024,y1024,z1024 voxel_size : x0.04607780456542968,y0.04607780456542968,z0.04607780456542968 center : x0,y0,z0
Show Acquisition geometry and full 3D sinogram.¶
In [8]:
show_geometry(ag3D)
Out[8]:
<cil.utilities.display.show_geometry at 0x7fdba39ddfc0>
In [9]:
show2D(data3D, slice_list = [('vertical',512), ('angle',800), ('horizontal',512)], cmap="inferno", num_cols=3, size=(15,15))
Out[9]:
<cil.utilities.display.show2D at 0x7fda77dacc10>
Slice through projections¶
In [10]:
islicer(data3D, direction=1, cmap="inferno")
HBox(children=(Output(), Box(children=(Play(value=800, interval=500, max=1600), VBox(children=(Label(value='Sl…
For demonstration purposes, we extract the central slice and select only 160 angles from the total 1601 angles.¶
- We use the
Slicerprocessor with step size of 10. - We use the
Binnerprocessor to crop and bin the acquisition data in order to reduce the field of view. - We use the
TransmissionAbsorptionConverterto convert from transmission measurements to absorption based on the Beer-Lambert law.
Note: To avoid circular artifacts in the reconstruction space, we subtract the mean value of a background Region of interest (ROI), i.e., ROI that does not contain the walnut.
In [11]:
# Extract vertical slice
data2D = data3D.get_slice(vertical='centre')
# Select every 10 angles
sliced_data = Slicer(roi={'angle':(0,1600,10)})(data2D)
# Reduce background regions
binned_data = Binner(roi={'horizontal':(120,-120,2)})(sliced_data)
# Create absorption data
absorption_data = TransmissionAbsorptionConverter()(binned_data)
# Remove circular artifacts
absorption_data -= np.mean(absorption_data.as_array()[80:100,0:30])
#Add some gaussian noise
absorption_data+=np.random.normal(0, 0.1*np.mean(absorption_data))
In [12]:
# Get Image and Acquisition geometries for one slice
ag2D = absorption_data.geometry
ag2D.set_angles(ag2D.angles, initial_angle=0.2, angle_unit='radian')
ig2D = ag2D.get_ImageGeometry()
In [13]:
print(" Acquisition Geometry 2D: {} with labels {}".format(ag2D.shape, ag2D.dimension_labels))
print(" Image Geometry 2D: {} with labels {}".format(ig2D.shape, ig2D.dimension_labels))
Acquisition Geometry 2D: (160, 392) with labels ('angle', 'horizontal')
Image Geometry 2D: (392, 392) with labels ('horizontal_y', 'horizontal_x')
Define Projection Operator¶
We can define our projection operator using our astra plugin that wraps the Astra-Toolbox library.
In [14]:
A = ProjectionOperator(ig2D, ag2D, device = "gpu")
PDHG - implicit TV (using CIL)¶
In [15]:
alpha_tv=0.03
F = 0.5 * L2NormSquared(b=absorption_data)
G = alpha_tv * TotalVariation(max_iteration=100, lower=0.)
K = A
# Setup and run PDHG
pdhg_tv_implicit_cil = PDHG(f = F, g = G, operator = K,
max_iteration = 200,
update_objective_interval = 5)
pdhg_tv_implicit_cil.run(verbose=1)
Iter Max Iter Time/Iter Objective
[s]
0 200 0.000 8.12271e+03
5 200 1.127 2.14863e+02
10 200 1.031 1.05014e+02
15 200 1.011 9.89716e+01
20 200 1.001 8.49490e+01
25 200 0.995 7.63064e+01
30 200 0.994 7.43717e+01
35 200 0.990 7.33507e+01
40 200 0.987 7.23722e+01
45 200 0.989 7.21300e+01
50 200 0.989 7.18422e+01
55 200 0.991 7.16938e+01
60 200 0.992 7.16845e+01
65 200 0.992 7.16506e+01
70 200 0.994 7.16081e+01
75 200 0.995 7.15891e+01
80 200 0.998 7.15806e+01
85 200 0.997 7.15811e+01
90 200 0.996 7.15837e+01
95 200 0.995 7.15837e+01
100 200 0.994 7.15835e+01
105 200 0.994 7.15814e+01
110 200 0.998 7.15780e+01
115 200 0.998 7.15776e+01
120 200 0.999 7.15789e+01
125 200 0.999 7.15804e+01
130 200 1.001 7.15808e+01
135 200 1.002 7.15801e+01
140 200 1.002 7.15802e+01
145 200 1.002 7.15799e+01
150 200 1.001 7.15800e+01
155 200 0.998 7.15797e+01
160 200 0.998 7.15798e+01
165 200 0.998 7.15794e+01
170 200 0.997 7.15794e+01
175 200 0.996 7.15791e+01
180 200 0.997 7.15795e+01
185 200 0.997 7.15794e+01
190 200 0.998 7.15799e+01
195 200 0.999 7.15798e+01
200 200 0.998 7.15801e+01
-------------------------------------------------------
200 200 0.998 7.15801e+01
Stop criterion has been reached.
PDHG - implicit TV (using CIL) with warm start¶
In [16]:
F = 0.5 * L2NormSquared(b=absorption_data)
G = alpha_tv * TotalVariation(max_iteration=5, lower=0., warmstart=True)
K = A
# Setup and run PDHG
pdhg_tv_implicit_cil_warm_start = PDHG(f = F, g = G, operator = K,
max_iteration = 200,
update_objective_interval = 5)
pdhg_tv_implicit_cil_warm_start.run(verbose=1)
Iter Max Iter Time/Iter Objective
[s]
0 200 0.000 8.12271e+03
5 200 0.064 1.58997e+02
10 200 0.065 1.05208e+02
15 200 0.063 9.13542e+01
20 200 0.063 8.63515e+01
25 200 0.062 8.49822e+01
30 200 0.062 7.81616e+01
35 200 0.062 7.55360e+01
40 200 0.062 7.43816e+01
45 200 0.063 7.35991e+01
50 200 0.063 7.35409e+01
55 200 0.063 7.29401e+01
60 200 0.063 7.24558e+01
65 200 0.063 7.20665e+01
70 200 0.063 7.18876e+01
75 200 0.064 7.18398e+01
80 200 0.064 7.17169e+01
85 200 0.064 7.15857e+01
90 200 0.064 7.15228e+01
95 200 0.064 7.14661e+01
100 200 0.064 7.14238e+01
105 200 0.064 7.14002e+01
110 200 0.064 7.13730e+01
115 200 0.064 7.13474e+01
120 200 0.064 7.13328e+01
125 200 0.065 7.13191e+01
130 200 0.065 7.13032e+01
135 200 0.064 7.12866e+01
140 200 0.064 7.12748e+01
145 200 0.064 7.12706e+01
150 200 0.065 7.12674e+01
155 200 0.065 7.12610e+01
160 200 0.065 7.12522e+01
165 200 0.065 7.12458e+01
170 200 0.064 7.12403e+01
175 200 0.064 7.12352e+01
180 200 0.064 7.12307e+01
185 200 0.064 7.12269e+01
190 200 0.064 7.12237e+01
195 200 0.064 7.12203e+01
200 200 0.064 7.12178e+01
-------------------------------------------------------
200 200 0.064 7.12178e+01
Stop criterion has been reached.
PDHG - implicit TV (using FGP_TV)¶
In [17]:
F = 0.5 * L2NormSquared(b=absorption_data)
G = (alpha_tv/ig2D.voxel_size_y) * FGP_TV(max_iteration=100, nonnegativity = True, device = 'gpu')
K = A
# Setup and run PDHG
pdhg_tv_implicit_regtk = PDHG(f = F, g = G, operator = K,
max_iteration = 200,
update_objective_interval = 5)
pdhg_tv_implicit_regtk.run(verbose=1)
Iter Max Iter Time/Iter Objective
[s]
0 200 0.000 8.12271e+03
5 200 0.025 2.14863e+02
10 200 0.024 1.05013e+02
15 200 0.023 9.89711e+01
20 200 0.023 8.49473e+01
25 200 0.023 7.63081e+01
30 200 0.023 7.43733e+01
35 200 0.023 7.33531e+01
40 200 0.023 7.23740e+01
45 200 0.023 7.21279e+01
50 200 0.023 7.18407e+01
55 200 0.023 7.16935e+01
60 200 0.022 7.16833e+01
65 200 0.022 7.16491e+01
70 200 0.022 7.16072e+01
75 200 0.022 7.15885e+01
80 200 0.023 7.15804e+01
85 200 0.022 7.15810e+01
90 200 0.022 7.15831e+01
95 200 0.022 7.15831e+01
100 200 0.022 7.15830e+01
105 200 0.022 7.15814e+01
110 200 0.022 7.15779e+01
115 200 0.022 7.15771e+01
120 200 0.022 7.15788e+01
125 200 0.022 7.15799e+01
130 200 0.022 7.15804e+01
135 200 0.022 7.15797e+01
140 200 0.022 7.15799e+01
145 200 0.022 7.15795e+01
150 200 0.022 7.15797e+01
155 200 0.022 7.15795e+01
160 200 0.022 7.15795e+01
165 200 0.022 7.15791e+01
170 200 0.022 7.15791e+01
175 200 0.022 7.15787e+01
180 200 0.022 7.15793e+01
185 200 0.022 7.15792e+01
190 200 0.022 7.15797e+01
195 200 0.022 7.15795e+01
200 200 0.022 7.15801e+01
-------------------------------------------------------
200 200 0.022 7.15801e+01
Stop criterion has been reached.
In [18]:
# Show reconstruction and ground truth
show2D([pdhg_tv_implicit_regtk.solution,
pdhg_tv_implicit_cil_warm_start.solution,
pdhg_tv_implicit_cil.solution,
(pdhg_tv_implicit_cil.solution-pdhg_tv_implicit_cil_warm_start.solution).abs()],
fix_range=[(0,0.055),(0,0.055),(0,0.055), None], num_cols=4,
title = ['TV (CIL)','TV (CIL - warm start)', 'TV (CCPI-RegTk)', 'CIL with/without warm start abs diff' ],
cmap = 'inferno')
print(' Absolute error between Regtk solution and warm start solution: ', np.linalg.norm(pdhg_tv_implicit_regtk.solution.as_array()-pdhg_tv_implicit_cil_warm_start.solution.as_array()))
# Plot middle line profile
show1D([pdhg_tv_implicit_regtk.solution,pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil_warm_start.solution], slice_list=[('horizontal_y',int(ig2D.voxel_num_y/2))],
label = ['TV (CCPi-RegTk)','TV (CIL)', 'TV (CIL-warm start)'], title='Middle Line Profiles')
print(pdhg_tv_implicit_regtk.objective)
print(pdhg_tv_implicit_cil.objective)
print(pdhg_tv_implicit_cil_warm_start.objective)
plt.figure()
iter_range = np.arange(0,201,5)
plt.semilogy(iter_range, pdhg_tv_implicit_regtk.objective, label='implicit PDHG (Regtk)')
plt.semilogy(iter_range, pdhg_tv_implicit_cil.objective, label='implicit PDHG (CIL)')
plt.semilogy(iter_range, pdhg_tv_implicit_cil_warm_start.objective, label='implicit PDHG (CIL Warm start)')
plt.xlabel('Iterations')
plt.ylabel('Objective function')
plt.ylim(71,73)
plt.legend()
plt.show()
Absolute error between Regtk solution and warm start solution: 0.03850115 [8122.7080078125, 214.86310958862305, 105.01270294189453, 98.97105407714844, 84.94731521606445, 76.30808639526367, 74.37328720092773, 73.35308074951172, 72.3740406036377, 72.12789344787598, 71.84066009521484, 71.69347953796387, 71.6832504272461, 71.64911651611328, 71.60719299316406, 71.58854675292969, 71.58039093017578, 71.58100509643555, 71.58312225341797, 71.58305931091309, 71.58297920227051, 71.58143043518066, 71.57787132263184, 71.57711410522461, 71.57875061035156, 71.57989501953125, 71.58036041259766, 71.57970428466797, 71.5799388885498, 71.57954788208008, 71.57974624633789, 71.57947158813477, 71.57951354980469, 71.57908248901367, 71.57912063598633, 71.57867240905762, 71.5793228149414, 71.57917213439941, 71.57971000671387, 71.57953643798828, 71.58005905151367] [8122.7080078125, 214.86308959960937, 105.01370727539063, 98.97164794921875, 84.94897247314452, 76.30642196655273, 74.37172058105469, 73.35065368652343, 72.37217788696289, 72.13004875183105, 71.8421728515625, 71.69379280090331, 71.68450500488281, 71.65055541992187, 71.60814971923827, 71.58905410766602, 71.58061897277832, 71.58109245300292, 71.58366699218749, 71.58365692138672, 71.58352951049804, 71.58144241333008, 71.5780233001709, 71.57756576538085, 71.57887283325195, 71.58042854309082, 71.58077728271485, 71.58009613037109, 71.58015953063965, 71.57985015869141, 71.58000396728515, 71.57969398498534, 71.57981239318848, 71.57939331054688, 71.57941101074218, 71.5791480255127, 71.57951141357421, 71.57944473266602, 71.57992416381836, 71.57977928161621, 71.58014846801757] [8122.7080078125, 158.9965350341797, 105.20823959350585, 91.35418487548827, 86.35148330688477, 84.98220001220703, 78.16156143188476, 75.53604095458985, 74.38160850524902, 73.59913505554199, 73.54087921142579, 72.940147857666, 72.4558348083496, 72.0665428161621, 71.88759712219237, 71.83979850769043, 71.71689781188965, 71.58573318481444, 71.52284255981445, 71.46614715576172, 71.42379371643067, 71.40017517089844, 71.37297790527344, 71.34740554809571, 71.33279983520508, 71.31906188964844, 71.3031842803955, 71.28661170959472, 71.27475189208984, 71.27059600830077, 71.26735610961913, 71.26102172851563, 71.25219146728516, 71.24583976745606, 71.24025482177734, 71.23521453857421, 71.23066017150879, 71.2269416809082, 71.22365844726562, 71.22034950256347, 71.21776794433593]
Trying different number of inner iterations for warm start¶
In [19]:
F = 0.5 * L2NormSquared(b=absorption_data)
G = (alpha_tv/ig2D.voxel_size_y) * FGP_TV(max_iteration=100, nonnegativity = True, device = 'gpu')
K = A
# Setup and run PDHG
pdhg_tv_implicit_regtk = PDHG(f = F, g = G, operator = K,
max_iteration = 400,
update_objective_interval = 5)
pdhg_tv_implicit_regtk.run(verbose=1)
F = 0.5 * L2NormSquared(b=absorption_data)
G = alpha_tv * TotalVariation(max_iteration=100, lower=0.)
K = A
# Setup and run PDHG
pdhg_tv_implicit_cil = PDHG(f = F, g = G, operator = K,
max_iteration = 400,
update_objective_interval = 5)
pdhg_tv_implicit_cil.run(verbose=1)
F = 0.5 * L2NormSquared(b=absorption_data)
K = A
plt.figure(figsize=(12,12))
iter_range = np.arange(0,401,5)
plt.semilogy(iter_range, pdhg_tv_implicit_regtk.objective, label='implicit PDHG (Regtk)')
plt.semilogy(iter_range, pdhg_tv_implicit_cil.objective, label='implicit PDHG (CIL)')
for i in range(0,30,5):
G = alpha_tv * TotalVariation(max_iteration=i, lower=0., warmstart=True)
# Setup and run PDHG
pdhg_tv_implicit_cil_warm_start = PDHG(f = F, g = G, operator = K,
max_iteration = 400,
update_objective_interval = 5)
pdhg_tv_implicit_cil_warm_start.run(verbose=1)
plt.semilogy(iter_range, pdhg_tv_implicit_cil_warm_start.objective, label='implicit PDHG (CIL Warm start '+str(i)+' iterations')
plt.xlabel('Iterations')
plt.ylabel('Objective function')
plt.ylim(70,76)
plt.legend()
plt.show()
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.022 2.14863e+02
10 400 0.022 1.05013e+02
15 400 0.022 9.89711e+01
20 400 0.022 8.49473e+01
25 400 0.022 7.63081e+01
30 400 0.022 7.43733e+01
35 400 0.022 7.33531e+01
40 400 0.022 7.23740e+01
45 400 0.022 7.21279e+01
50 400 0.022 7.18407e+01
55 400 0.022 7.16935e+01
60 400 0.022 7.16833e+01
65 400 0.022 7.16491e+01
70 400 0.022 7.16072e+01
75 400 0.022 7.15885e+01
80 400 0.022 7.15804e+01
85 400 0.022 7.15810e+01
90 400 0.022 7.15831e+01
95 400 0.022 7.15831e+01
100 400 0.022 7.15830e+01
105 400 0.022 7.15814e+01
110 400 0.022 7.15779e+01
115 400 0.022 7.15771e+01
120 400 0.022 7.15788e+01
125 400 0.022 7.15799e+01
130 400 0.022 7.15804e+01
135 400 0.022 7.15797e+01
140 400 0.022 7.15799e+01
145 400 0.022 7.15795e+01
150 400 0.022 7.15797e+01
155 400 0.022 7.15795e+01
160 400 0.022 7.15795e+01
165 400 0.022 7.15791e+01
170 400 0.022 7.15791e+01
175 400 0.022 7.15787e+01
180 400 0.022 7.15793e+01
185 400 0.022 7.15792e+01
190 400 0.022 7.15797e+01
195 400 0.022 7.15795e+01
200 400 0.022 7.15801e+01
205 400 0.022 7.15796e+01
210 400 0.022 7.15799e+01
215 400 0.022 7.15797e+01
220 400 0.022 7.15799e+01
225 400 0.022 7.15798e+01
230 400 0.022 7.15800e+01
235 400 0.022 7.15797e+01
240 400 0.022 7.15799e+01
245 400 0.022 7.15796e+01
250 400 0.022 7.15799e+01
255 400 0.023 7.15796e+01
260 400 0.023 7.15798e+01
265 400 0.023 7.15796e+01
270 400 0.023 7.15797e+01
275 400 0.023 7.15795e+01
280 400 0.023 7.15798e+01
285 400 0.023 7.15796e+01
290 400 0.023 7.15798e+01
295 400 0.023 7.15796e+01
300 400 0.023 7.15798e+01
305 400 0.023 7.15796e+01
310 400 0.023 7.15798e+01
315 400 0.023 7.15795e+01
320 400 0.023 7.15798e+01
325 400 0.023 7.15796e+01
330 400 0.023 7.15798e+01
335 400 0.023 7.15795e+01
340 400 0.023 7.15798e+01
345 400 0.023 7.15796e+01
350 400 0.023 7.15797e+01
355 400 0.022 7.15795e+01
360 400 0.023 7.15798e+01
365 400 0.023 7.15796e+01
370 400 0.022 7.15797e+01
375 400 0.022 7.15795e+01
380 400 0.022 7.15798e+01
385 400 0.023 7.15796e+01
390 400 0.023 7.15797e+01
395 400 0.023 7.15795e+01
400 400 0.023 7.15798e+01
-------------------------------------------------------
400 400 0.023 7.15798e+01
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.988 2.14863e+02
10 400 0.985 1.05014e+02
15 400 0.975 9.89716e+01
20 400 0.989 8.49490e+01
25 400 0.998 7.63064e+01
30 400 0.994 7.43717e+01
35 400 0.999 7.33507e+01
40 400 1.001 7.23722e+01
45 400 0.997 7.21300e+01
50 400 0.990 7.18422e+01
55 400 0.993 7.16938e+01
60 400 0.988 7.16845e+01
65 400 0.993 7.16506e+01
70 400 0.997 7.16081e+01
75 400 0.995 7.15891e+01
80 400 0.996 7.15806e+01
85 400 1.000 7.15811e+01
90 400 1.002 7.15837e+01
95 400 1.003 7.15837e+01
100 400 1.002 7.15835e+01
105 400 1.001 7.15814e+01
110 400 1.000 7.15780e+01
115 400 1.002 7.15776e+01
120 400 1.002 7.15789e+01
125 400 1.002 7.15804e+01
130 400 1.004 7.15808e+01
135 400 1.005 7.15801e+01
140 400 1.007 7.15802e+01
145 400 1.005 7.15799e+01
150 400 1.006 7.15800e+01
155 400 1.004 7.15797e+01
160 400 1.004 7.15798e+01
165 400 1.004 7.15794e+01
170 400 1.005 7.15794e+01
175 400 1.001 7.15791e+01
180 400 1.002 7.15795e+01
185 400 1.002 7.15794e+01
190 400 1.002 7.15799e+01
195 400 1.003 7.15798e+01
200 400 1.004 7.15801e+01
205 400 1.005 7.15799e+01
210 400 1.004 7.15800e+01
215 400 1.005 7.15797e+01
220 400 1.004 7.15801e+01
225 400 1.003 7.15798e+01
230 400 1.002 7.15800e+01
235 400 1.002 7.15798e+01
240 400 1.002 7.15800e+01
245 400 1.001 7.15798e+01
250 400 1.002 7.15800e+01
255 400 1.001 7.15797e+01
260 400 1.001 7.15800e+01
265 400 1.002 7.15797e+01
270 400 1.001 7.15799e+01
275 400 1.001 7.15796e+01
280 400 1.001 7.15799e+01
285 400 1.001 7.15797e+01
290 400 1.001 7.15799e+01
295 400 1.002 7.15797e+01
300 400 1.002 7.15799e+01
305 400 1.003 7.15797e+01
310 400 1.002 7.15799e+01
315 400 1.002 7.15796e+01
320 400 1.002 7.15799e+01
325 400 1.003 7.15797e+01
330 400 1.004 7.15799e+01
335 400 1.004 7.15796e+01
340 400 1.003 7.15800e+01
345 400 1.003 7.15797e+01
350 400 1.003 7.15799e+01
355 400 1.004 7.15796e+01
360 400 1.004 7.15799e+01
365 400 1.004 7.15797e+01
370 400 1.005 7.15799e+01
375 400 1.005 7.15796e+01
380 400 1.005 7.15799e+01
385 400 1.005 7.15797e+01
390 400 1.005 7.15799e+01
395 400 1.005 7.15796e+01
400 400 1.005 7.15800e+01
-------------------------------------------------------
400 400 1.005 7.15800e+01
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.013 8.12271e+03
10 400 0.013 8.12271e+03
15 400 0.013 8.12271e+03
20 400 0.013 8.12271e+03
25 400 0.013 8.12271e+03
30 400 0.013 8.12271e+03
35 400 0.013 8.12271e+03
40 400 0.013 8.12271e+03
45 400 0.013 8.12271e+03
50 400 0.013 8.12271e+03
55 400 0.013 8.12271e+03
60 400 0.013 8.12271e+03
65 400 0.013 8.12271e+03
70 400 0.013 8.12271e+03
75 400 0.013 8.12271e+03
80 400 0.013 8.12271e+03
85 400 0.013 8.12271e+03
90 400 0.013 8.12271e+03
95 400 0.013 8.12271e+03
100 400 0.013 8.12271e+03
105 400 0.013 8.12271e+03
110 400 0.013 8.12271e+03
115 400 0.013 8.12271e+03
120 400 0.013 8.12271e+03
125 400 0.013 8.12271e+03
130 400 0.013 8.12271e+03
135 400 0.013 8.12271e+03
140 400 0.013 8.12271e+03
145 400 0.013 8.12271e+03
150 400 0.013 8.12271e+03
155 400 0.013 8.12271e+03
160 400 0.013 8.12271e+03
165 400 0.013 8.12271e+03
170 400 0.013 8.12271e+03
175 400 0.013 8.12271e+03
180 400 0.013 8.12271e+03
185 400 0.013 8.12271e+03
190 400 0.013 8.12271e+03
195 400 0.013 8.12271e+03
200 400 0.013 8.12271e+03
205 400 0.013 8.12271e+03
210 400 0.013 8.12271e+03
215 400 0.013 8.12271e+03
220 400 0.013 8.12271e+03
225 400 0.013 8.12271e+03
230 400 0.013 8.12271e+03
235 400 0.013 8.12271e+03
240 400 0.013 8.12271e+03
245 400 0.013 8.12271e+03
250 400 0.013 8.12271e+03
255 400 0.013 8.12271e+03
260 400 0.013 8.12271e+03
265 400 0.013 8.12271e+03
270 400 0.013 8.12271e+03
275 400 0.013 8.12271e+03
280 400 0.013 8.12271e+03
285 400 0.013 8.12271e+03
290 400 0.013 8.12271e+03
295 400 0.013 8.12271e+03
300 400 0.013 8.12271e+03
305 400 0.013 8.12271e+03
310 400 0.013 8.12271e+03
315 400 0.013 8.12271e+03
320 400 0.013 8.12271e+03
325 400 0.013 8.12271e+03
330 400 0.013 8.12271e+03
335 400 0.013 8.12271e+03
340 400 0.013 8.12271e+03
345 400 0.013 8.12271e+03
350 400 0.013 8.12271e+03
355 400 0.013 8.12271e+03
360 400 0.013 8.12271e+03
365 400 0.013 8.12271e+03
370 400 0.013 8.12271e+03
375 400 0.013 8.12271e+03
380 400 0.013 8.12271e+03
385 400 0.013 8.12271e+03
390 400 0.013 8.12271e+03
395 400 0.013 8.12271e+03
400 400 0.013 8.12271e+03
-------------------------------------------------------
400 400 0.013 8.12271e+03
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.064 1.58997e+02
10 400 0.062 1.05208e+02
15 400 0.062 9.13542e+01
20 400 0.063 8.63515e+01
25 400 0.062 8.49822e+01
30 400 0.063 7.81616e+01
35 400 0.061 7.55360e+01
40 400 0.061 7.43816e+01
45 400 0.062 7.35991e+01
50 400 0.063 7.35409e+01
55 400 0.064 7.29401e+01
60 400 0.064 7.24558e+01
65 400 0.065 7.20665e+01
70 400 0.066 7.18876e+01
75 400 0.066 7.18398e+01
80 400 0.066 7.17169e+01
85 400 0.067 7.15857e+01
90 400 0.066 7.15228e+01
95 400 0.066 7.14661e+01
100 400 0.066 7.14238e+01
105 400 0.067 7.14002e+01
110 400 0.067 7.13730e+01
115 400 0.067 7.13474e+01
120 400 0.067 7.13328e+01
125 400 0.067 7.13191e+01
130 400 0.067 7.13032e+01
135 400 0.067 7.12866e+01
140 400 0.067 7.12748e+01
145 400 0.067 7.12706e+01
150 400 0.067 7.12674e+01
155 400 0.067 7.12610e+01
160 400 0.067 7.12522e+01
165 400 0.067 7.12458e+01
170 400 0.067 7.12403e+01
175 400 0.067 7.12352e+01
180 400 0.067 7.12307e+01
185 400 0.067 7.12269e+01
190 400 0.068 7.12237e+01
195 400 0.068 7.12203e+01
200 400 0.068 7.12178e+01
205 400 0.068 7.12156e+01
210 400 0.068 7.12127e+01
215 400 0.068 7.12095e+01
220 400 0.068 7.12063e+01
225 400 0.068 7.12036e+01
230 400 0.068 7.12013e+01
235 400 0.068 7.11993e+01
240 400 0.068 7.11975e+01
245 400 0.067 7.11953e+01
250 400 0.067 7.11929e+01
255 400 0.067 7.11905e+01
260 400 0.067 7.11883e+01
265 400 0.067 7.11861e+01
270 400 0.067 7.11837e+01
275 400 0.067 7.11818e+01
280 400 0.067 7.11802e+01
285 400 0.067 7.11784e+01
290 400 0.067 7.11766e+01
295 400 0.067 7.11746e+01
300 400 0.067 7.11726e+01
305 400 0.067 7.11708e+01
310 400 0.067 7.11692e+01
315 400 0.067 7.11679e+01
320 400 0.067 7.11668e+01
325 400 0.067 7.11655e+01
330 400 0.067 7.11640e+01
335 400 0.067 7.11627e+01
340 400 0.067 7.11614e+01
345 400 0.066 7.11601e+01
350 400 0.066 7.11589e+01
355 400 0.066 7.11577e+01
360 400 0.066 7.11565e+01
365 400 0.066 7.11553e+01
370 400 0.066 7.11540e+01
375 400 0.066 7.11527e+01
380 400 0.066 7.11515e+01
385 400 0.066 7.11504e+01
390 400 0.066 7.11494e+01
395 400 0.066 7.11484e+01
400 400 0.066 7.11474e+01
-------------------------------------------------------
400 400 0.066 7.11474e+01
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.113 1.71721e+02
10 400 0.112 1.00645e+02
15 400 0.112 9.14914e+01
20 400 0.110 8.53634e+01
25 400 0.113 7.99447e+01
30 400 0.113 7.51679e+01
35 400 0.113 7.36352e+01
40 400 0.114 7.30386e+01
45 400 0.114 7.24725e+01
50 400 0.114 7.22318e+01
55 400 0.114 7.19175e+01
60 400 0.114 7.16012e+01
65 400 0.115 7.14779e+01
70 400 0.115 7.13913e+01
75 400 0.116 7.13185e+01
80 400 0.116 7.12793e+01
85 400 0.116 7.12432e+01
90 400 0.116 7.12276e+01
95 400 0.117 7.12127e+01
100 400 0.116 7.11986e+01
105 400 0.116 7.11829e+01
110 400 0.116 7.11699e+01
115 400 0.117 7.11618e+01
120 400 0.117 7.11575e+01
125 400 0.117 7.11538e+01
130 400 0.117 7.11487e+01
135 400 0.117 7.11453e+01
140 400 0.117 7.11418e+01
145 400 0.117 7.11390e+01
150 400 0.117 7.11374e+01
155 400 0.117 7.11358e+01
160 400 0.117 7.11328e+01
165 400 0.117 7.11297e+01
170 400 0.117 7.11274e+01
175 400 0.117 7.11258e+01
180 400 0.117 7.11245e+01
185 400 0.116 7.11232e+01
190 400 0.116 7.11220e+01
195 400 0.115 7.11208e+01
200 400 0.115 7.11193e+01
205 400 0.115 7.11179e+01
210 400 0.115 7.11165e+01
215 400 0.115 7.11153e+01
220 400 0.115 7.11141e+01
225 400 0.115 7.11129e+01
230 400 0.115 7.11117e+01
235 400 0.115 7.11107e+01
240 400 0.115 7.11099e+01
245 400 0.115 7.11091e+01
250 400 0.115 7.11082e+01
255 400 0.115 7.11072e+01
260 400 0.115 7.11063e+01
265 400 0.115 7.11054e+01
270 400 0.115 7.11044e+01
275 400 0.114 7.11035e+01
280 400 0.114 7.11028e+01
285 400 0.114 7.11023e+01
290 400 0.114 7.11016e+01
295 400 0.114 7.11008e+01
300 400 0.114 7.11000e+01
305 400 0.114 7.10994e+01
310 400 0.113 7.10989e+01
315 400 0.113 7.10984e+01
320 400 0.113 7.10978e+01
325 400 0.113 7.10971e+01
330 400 0.113 7.10964e+01
335 400 0.113 7.10957e+01
340 400 0.113 7.10950e+01
345 400 0.113 7.10944e+01
350 400 0.113 7.10938e+01
355 400 0.113 7.10934e+01
360 400 0.113 7.10929e+01
365 400 0.113 7.10923e+01
370 400 0.113 7.10915e+01
375 400 0.113 7.10907e+01
380 400 0.113 7.10902e+01
385 400 0.113 7.10900e+01
390 400 0.113 7.10897e+01
395 400 0.114 7.10893e+01
400 400 0.114 7.10888e+01
-------------------------------------------------------
400 400 0.114 7.10888e+01
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.167 1.85117e+02
10 400 0.172 9.92438e+01
15 400 0.170 9.37167e+01
20 400 0.173 8.46188e+01
25 400 0.176 7.87773e+01
30 400 0.173 7.41055e+01
35 400 0.172 7.32907e+01
40 400 0.172 7.24861e+01
45 400 0.170 7.20799e+01
50 400 0.170 7.18621e+01
55 400 0.170 7.15708e+01
60 400 0.169 7.14134e+01
65 400 0.169 7.13155e+01
70 400 0.168 7.12594e+01
75 400 0.167 7.12148e+01
80 400 0.167 7.11902e+01
85 400 0.166 7.11680e+01
90 400 0.167 7.11507e+01
95 400 0.167 7.11420e+01
100 400 0.167 7.11350e+01
105 400 0.167 7.11272e+01
110 400 0.167 7.11229e+01
115 400 0.167 7.11169e+01
120 400 0.167 7.11125e+01
125 400 0.167 7.11102e+01
130 400 0.168 7.11088e+01
135 400 0.168 7.11069e+01
140 400 0.168 7.11049e+01
145 400 0.169 7.11029e+01
150 400 0.168 7.11009e+01
155 400 0.168 7.10991e+01
160 400 0.168 7.10977e+01
165 400 0.167 7.10966e+01
170 400 0.167 7.10953e+01
175 400 0.168 7.10940e+01
180 400 0.167 7.10925e+01
185 400 0.167 7.10909e+01
190 400 0.167 7.10897e+01
195 400 0.167 7.10890e+01
200 400 0.167 7.10885e+01
205 400 0.167 7.10879e+01
210 400 0.167 7.10870e+01
215 400 0.167 7.10861e+01
220 400 0.167 7.10853e+01
225 400 0.167 7.10845e+01
230 400 0.167 7.10837e+01
235 400 0.167 7.10830e+01
240 400 0.166 7.10824e+01
245 400 0.166 7.10818e+01
250 400 0.166 7.10811e+01
255 400 0.166 7.10807e+01
260 400 0.166 7.10803e+01
265 400 0.166 7.10798e+01
270 400 0.166 7.10794e+01
275 400 0.166 7.10789e+01
280 400 0.166 7.10786e+01
285 400 0.167 7.10782e+01
290 400 0.167 7.10779e+01
295 400 0.167 7.10774e+01
300 400 0.167 7.10770e+01
305 400 0.167 7.10766e+01
310 400 0.167 7.10762e+01
315 400 0.167 7.10759e+01
320 400 0.167 7.10757e+01
325 400 0.167 7.10754e+01
330 400 0.167 7.10750e+01
335 400 0.167 7.10745e+01
340 400 0.167 7.10742e+01
345 400 0.167 7.10739e+01
350 400 0.166 7.10738e+01
355 400 0.166 7.10737e+01
360 400 0.167 7.10735e+01
365 400 0.167 7.10731e+01
370 400 0.167 7.10728e+01
375 400 0.167 7.10726e+01
380 400 0.167 7.10724e+01
385 400 0.167 7.10723e+01
390 400 0.167 7.10721e+01
395 400 0.167 7.10719e+01
400 400 0.167 7.10718e+01
-------------------------------------------------------
400 400 0.167 7.10718e+01
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.212 1.93965e+02
10 400 0.217 9.95994e+01
15 400 0.223 9.49914e+01
20 400 0.221 8.44162e+01
25 400 0.220 7.79351e+01
30 400 0.219 7.37225e+01
35 400 0.217 7.31768e+01
40 400 0.216 7.22947e+01
45 400 0.216 7.19475e+01
50 400 0.216 7.16872e+01
55 400 0.216 7.14623e+01
60 400 0.215 7.13256e+01
65 400 0.214 7.12585e+01
70 400 0.214 7.12063e+01
75 400 0.214 7.11779e+01
80 400 0.214 7.11522e+01
85 400 0.214 7.11339e+01
90 400 0.214 7.11201e+01
95 400 0.214 7.11119e+01
100 400 0.213 7.11078e+01
105 400 0.212 7.11044e+01
110 400 0.212 7.10996e+01
115 400 0.212 7.10940e+01
120 400 0.212 7.10895e+01
125 400 0.212 7.10876e+01
130 400 0.212 7.10869e+01
135 400 0.211 7.10863e+01
140 400 0.211 7.10850e+01
145 400 0.210 7.10832e+01
150 400 0.210 7.10817e+01
155 400 0.210 7.10807e+01
160 400 0.210 7.10798e+01
165 400 0.210 7.10790e+01
170 400 0.210 7.10782e+01
175 400 0.210 7.10777e+01
180 400 0.210 7.10771e+01
185 400 0.210 7.10765e+01
190 400 0.211 7.10759e+01
195 400 0.210 7.10755e+01
200 400 0.210 7.10751e+01
205 400 0.210 7.10743e+01
210 400 0.210 7.10736e+01
215 400 0.210 7.10732e+01
220 400 0.210 7.10730e+01
225 400 0.210 7.10727e+01
230 400 0.210 7.10723e+01
235 400 0.210 7.10718e+01
240 400 0.211 7.10715e+01
245 400 0.211 7.10714e+01
250 400 0.211 7.10711e+01
255 400 0.212 7.10709e+01
260 400 0.212 7.10707e+01
265 400 0.212 7.10705e+01
270 400 0.212 7.10703e+01
275 400 0.212 7.10701e+01
280 400 0.212 7.10699e+01
285 400 0.212 7.10698e+01
290 400 0.212 7.10696e+01
295 400 0.212 7.10694e+01
300 400 0.212 7.10693e+01
305 400 0.212 7.10691e+01
310 400 0.212 7.10689e+01
315 400 0.212 7.10688e+01
320 400 0.212 7.10687e+01
325 400 0.212 7.10685e+01
330 400 0.212 7.10684e+01
335 400 0.212 7.10682e+01
340 400 0.212 7.10681e+01
345 400 0.212 7.10679e+01
350 400 0.212 7.10677e+01
355 400 0.212 7.10676e+01
360 400 0.212 7.10674e+01
365 400 0.212 7.10673e+01
370 400 0.212 7.10671e+01
375 400 0.212 7.10670e+01
380 400 0.212 7.10669e+01
385 400 0.212 7.10668e+01
390 400 0.212 7.10667e+01
395 400 0.212 7.10666e+01
400 400 0.212 7.10665e+01
-------------------------------------------------------
400 400 0.212 7.10665e+01
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 400 0.000 8.12271e+03
5 400 0.272 1.99302e+02
10 400 0.264 1.00195e+02
15 400 0.262 9.59857e+01
20 400 0.260 8.43527e+01
25 400 0.260 7.73410e+01
30 400 0.262 7.35644e+01
35 400 0.261 7.31222e+01
40 400 0.263 7.21785e+01
45 400 0.262 7.18266e+01
50 400 0.262 7.16061e+01
55 400 0.262 7.14000e+01
60 400 0.263 7.12714e+01
65 400 0.264 7.12256e+01
70 400 0.263 7.11695e+01
75 400 0.263 7.11531e+01
80 400 0.264 7.11269e+01
85 400 0.264 7.11078e+01
90 400 0.264 7.10990e+01
95 400 0.264 7.10959e+01
100 400 0.264 7.10917e+01
105 400 0.266 7.10872e+01
110 400 0.266 7.10839e+01
115 400 0.266 7.10811e+01
120 400 0.267 7.10788e+01
125 400 0.267 7.10775e+01
130 400 0.268 7.10766e+01
135 400 0.268 7.10759e+01
140 400 0.269 7.10750e+01
145 400 0.269 7.10739e+01
150 400 0.270 7.10731e+01
155 400 0.270 7.10725e+01
160 400 0.270 7.10717e+01
165 400 0.270 7.10712e+01
170 400 0.270 7.10709e+01
175 400 0.269 7.10705e+01
180 400 0.270 7.10702e+01
185 400 0.270 7.10699e+01
190 400 0.270 7.10697e+01
195 400 0.270 7.10696e+01
200 400 0.270 7.10693e+01
205 400 0.269 7.10690e+01
210 400 0.269 7.10688e+01
215 400 0.269 7.10685e+01
220 400 0.268 7.10683e+01
225 400 0.268 7.10682e+01
230 400 0.268 7.10679e+01
235 400 0.268 7.10677e+01
240 400 0.268 7.10674e+01
245 400 0.268 7.10672e+01
250 400 0.267 7.10670e+01
255 400 0.267 7.10668e+01
260 400 0.267 7.10667e+01
265 400 0.266 7.10665e+01
270 400 0.266 7.10664e+01
275 400 0.266 7.10663e+01
280 400 0.266 7.10662e+01
285 400 0.266 7.10661e+01
290 400 0.266 7.10659e+01
295 400 0.266 7.10658e+01
300 400 0.266 7.10658e+01
305 400 0.266 7.10657e+01
310 400 0.267 7.10656e+01
315 400 0.266 7.10655e+01
320 400 0.266 7.10654e+01
325 400 0.267 7.10653e+01
330 400 0.266 7.10652e+01
335 400 0.266 7.10652e+01
340 400 0.266 7.10651e+01
345 400 0.266 7.10650e+01
350 400 0.266 7.10650e+01
355 400 0.266 7.10649e+01
360 400 0.266 7.10649e+01
365 400 0.266 7.10648e+01
370 400 0.266 7.10648e+01
375 400 0.266 7.10647e+01
380 400 0.266 7.10647e+01
385 400 0.266 7.10646e+01
390 400 0.266 7.10646e+01
395 400 0.266 7.10645e+01
400 400 0.266 7.10645e+01
-------------------------------------------------------
400 400 0.266 7.10645e+01
Stop criterion has been reached.
Absolute error in the prximal calcultion printed for each iteration¶
In [22]:
F = 0.5 * L2NormSquared(b=absorption_data)
G_new = alpha_tv*TotalVariation(max_iteration=5, lower=0., warmstart=True)
G_FGP_TV=(alpha_tv/ig2D.voxel_size_y)*FGP_TV(max_iteration=500, nonnegativity = True, device = 'gpu')
K = A
# Setup and run PDHG
pdhg_tv_implicit_cil = PDHG(f = F, g = G_FGP_TV, operator = K,
max_iteration = 400,
update_objective_interval = 50)
hold=np.zeros(400)
for i in range(400):
pdhg_tv_implicit_cil.__next__()
#print(np.linalg.norm(pdhg_tv_implicit_cil.solution.as_array()))
#print(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, 0.5).as_array())
#print(G_new.proximal(pdhg_tv_implicit_cil.solution, 0.5).as_array())
hold[i]=np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array()-G_new.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())/np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())
print('{0:.15f}'.format(hold[i]))
plt.figure()
plt.semilogy(range(400), hold)
plt.title('Comparison of proximal values of warm start TV and CCPI TV reg at PDHG (Regtk) iteration values')
plt.ylabel('Absolute error')
plt.xlabel('Iteration number')
plt.show()
print('Next example ')
F = 0.5 * L2NormSquared(b=absorption_data)
G_new = alpha_tv*TotalVariation(max_iteration=5, warmstart=True)
G_FGP_TV=(alpha_tv/ig2D.voxel_size_y)*FGP_TV(max_iteration=500, device = 'gpu')
K = A
# Setup and run PDHG
pdhg_tv_implicit_cil = PDHG(f = F, g = G_new, operator = K,
max_iteration = 400,
update_objective_interval = 50)
for i in range(400):
pdhg_tv_implicit_cil.__next__()
hold[i]=np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array()-G_new.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())/np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())
print('{0:.15f}'.format(hold[i]))
plt.figure()
plt.semilogy(range(400), hold)
plt.title('Comparison of proximal values of warm start TV and CCPI TV reg at PDHG (CIL-warm start) iteration values')
plt.ylabel('Absolute error')
plt.xlabel('Iteration number')
plt.show()
/tmp/ipykernel_1298247/2984884879.py:16: RuntimeWarning: invalid value encountered in float_scalars hold[i]=np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array()-G_new.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())/np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())
nan 0.020799819380045 0.021384317427874 0.033312588930130 0.041921649128199 0.043050855398178 0.044340092688799 0.046354092657566 0.048461508005857 0.049331754446030 0.049025006592274 0.047629836946726 0.044929951429367 0.041443478316069 0.037805583328009 0.034432861953974 0.031737826764584 0.029794363304973 0.028351578861475 0.027389535680413 0.026745310053229 0.026083055883646 0.025350060313940 0.024762144312263 0.024176770821214 0.023660127073526 0.023142427206039 0.022607652470469 0.022292049601674 0.021852213889360 0.021534141153097 0.021271744742990 0.021099086850882 0.020900765433908 0.020724283531308 0.020611926913261 0.020454827696085 0.020272128283978 0.020074548199773 0.019858028739691 0.019637780264020 0.019415618851781 0.019192153587937 0.018973046913743 0.018754156306386 0.018522590398788 0.018290841951966 0.018084855750203 0.017882673069835 0.017698820680380 0.017526533454657 0.017383899539709 0.017282161861658 0.017190827056766 0.017112428322434 0.017024483531713 0.016940958797932 0.016864009201527 0.016790103167295 0.016718877479434 0.016651149839163 0.016568399965763 0.016478333622217 0.016384471207857 0.016295792534947 0.016205381602049 0.016105266287923 0.015996431931853 0.015885453671217 0.015774061903358 0.015666736289859 0.015563121996820 0.015463100746274 0.015366944484413 0.015273279510438 0.015183814801276 0.015096340328455 0.015014175325632 0.014936955645680 0.014866853132844 0.014802990481257 0.014743614010513 0.014688221737742 0.014634391292930 0.014582896605134 0.014532226137817 0.014480254612863 0.014426308684051 0.014372186735272 0.014317475259304 0.014264217577875 0.014211135916412 0.014157537370920 0.014103732071817 0.014049308374524 0.013994182460010 0.013938329182565 0.013881776481867 0.013824261724949 0.013766625896096 0.013709085993469 0.013652008026838 0.013595969416201 0.013540034182370 0.013484480790794 0.013429556973279 0.013375485315919 0.013322478160262 0.013270355761051 0.013219906017184 0.013170844875276 0.013122999109328 0.013076215982437 0.013030333444476 0.012985426932573 0.012941243126988 0.012897842563689 0.012855316512287 0.012813423760235 0.012772134505212 0.012731450609863 0.012691387906671 0.012651369906962 0.012611770071089 0.012572408653796 0.012533335015178 0.012494218535721 0.012455478310585 0.012416675686836 0.012377838604152 0.012338911183178 0.012299864552915 0.012260959483683 0.012222290039062 0.012183626182377 0.012145542539656 0.012107403017581 0.012069347314537 0.012031538411975 0.011994232423604 0.011957048438489 0.011920105665922 0.011883604340255 0.011847326532006 0.011811387725174 0.011775834485888 0.011740620248020 0.011705878190696 0.011671487241983 0.011637453921139 0.011603986844420 0.011571066454053 0.011538523249328 0.011506373994052 0.011474826373160 0.011443492956460 0.011412693187594 0.011382143944502 0.011351940222085 0.011322130449116 0.011292596347630 0.011263293214142 0.011234030127525 0.011205048300326 0.011176412925124 0.011147704906762 0.011119483038783 0.011091241613030 0.011063321493566 0.011035440489650 0.011007796041667 0.010980352759361 0.010953054763377 0.010925811715424 0.010898941196501 0.010872153565288 0.010845541954041 0.010819165967405 0.010792879387736 0.010766777209938 0.010741044767201 0.010715378448367 0.010689979419112 0.010664644651115 0.010639727115631 0.010614922270179 0.010590255260468 0.010565860196948 0.010541666299105 0.010517595335841 0.010493807494640 0.010470166802406 0.010446652770042 0.010423382744193 0.010400356724858 0.010377462953329 0.010354672558606 0.010332133620977 0.010309672914445 0.010287385433912 0.010265220887959 0.010243286378682 0.010221400298178 0.010199592448771 0.010177855379879 0.010156258940697 0.010134860873222 0.010113518685102 0.010092250071466 0.010071176104248 0.010050274431705 0.010029318742454 0.010008671320975 0.009988023899496 0.009967545978725 0.009947078302503 0.009926824830472 0.009906590916216 0.009886508807540 0.009866550564766 0.009846702218056 0.009826914407313 0.009807257913053 0.009787688963115 0.009768088348210 0.009748763404787 0.009729426354170 0.009710201993585 0.009691022336483 0.009672027081251 0.009653185494244 0.009634474292397 0.009615871123970 0.009597321972251 0.009578797966242 0.009560639970005 0.009542313404381 0.009524028748274 0.009505901485682 0.009487580507994 0.009469498880208 0.009451366961002 0.009433213621378 0.009415220469236 0.009397248737514 0.009379052557051 0.009361074306071 0.009342967532575 0.009324998594820 0.009307011961937 0.009289229288697 0.009271519258618 0.009253704920411 0.009235978126526 0.009218352846801 0.009200588800013 0.009183201007545 0.009165644645691 0.009148376993835 0.009130938909948 0.009113729931414 0.009096583351493 0.009079404175282 0.009062312543392 0.009045282378793 0.009028269909322 0.009011310525239 0.008994473144412 0.008977798745036 0.008960950188339 0.008944416418672 0.008927810937166 0.008911322802305 0.008894818834960 0.008878519758582 0.008862196467817 0.008846043609083 0.008829923346639 0.008813810534775 0.008797804825008 0.008781856857240 0.008766098879278 0.008750371634960 0.008734646253288 0.008719237521291 0.008703622967005 0.008688339963555 0.008672747761011 0.008657477796078 0.008642172440886 0.008626976050436 0.008611928671598 0.008596763946116 0.008581751026213 0.008566764183342 0.008551796898246 0.008536930195987 0.008522063493729 0.008507140912116 0.008492273278534 0.008477303199470 0.008462473750114 0.008447571657598 0.008432507514954 0.008417699486017 0.008402906358242 0.008388193324208 0.008373588323593 0.008358911611140 0.008344350382686 0.008329815231264 0.008315353654325 0.008301015943289 0.008286654017866 0.008272417820990 0.008258274756372 0.008244100958109 0.008229991421103 0.008215999230742 0.008202060125768 0.008188263513148 0.008174491114914 0.008160576224327 0.008147012442350 0.008133314549923 0.008119843900204 0.008106333203614 0.008092856034636 0.008079635910690 0.008066354319453 0.008053119294345 0.008039951324463 0.008026870898902 0.008013871498406 0.008000898174942 0.007988117635250 0.007975244894624 0.007962475530803 0.007949607446790 0.007937044836581 0.007924462668598 0.007911867462099 0.007899365387857 0.007886989042163 0.007874594070017 0.007862209342420 0.007850079797208 0.007837787270546 0.007825664244592 0.007813588716090 0.007801536004990 0.007789441850036 0.007777496241033 0.007765587884933 0.007753815967590 0.007741956505924 0.007730168756098 0.007718429900706 0.007706765551120 0.007695002947003 0.007683433592319 0.007671759463847 0.007660272996873 0.007648772560060 0.007637389935553 0.007625898811966 0.007614559959620 0.007603175006807 0.007591864559799 0.007580614183098 0.007569485343993 0.007558348122984 0.007547246292233 0.007536197546870 0.007525241933763 0.007514239754528 0.007503360044211 0.007492462638766 0.007481636013836 0.007470798213035 0.007460161112249 0.007449476979673 0.007438709493726 0.007428253535181 0.007417724002153 0.007407157216221 0.007396685425192 0.007386199664325 0.007375694345683 0.007365324068815 0.007355127949268 0.007344848942012 0.007334636058658 0.007324434816837 0.007314192596823 0.007304301951081 0.007294111885130 0.007284175604582 0.007274150848389 0.007264244370162 Next example nan
/tmp/ipykernel_1298247/2984884879.py:41: RuntimeWarning: invalid value encountered in float_scalars hold[i]=np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array()-G_new.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())/np.linalg.norm(G_FGP_TV.proximal(pdhg_tv_implicit_cil.solution, pdhg_tv_implicit_cil.tau).as_array())
0.029623007401824 0.051463138312101 0.060652390122414 0.064053997397423 0.080338470637798 0.091146059334278 0.086835384368896 0.075095549225807 0.066289030015469 0.065516695380211 0.065589822828770 0.063110716640949 0.058823939412832 0.057353351265192 0.061581049114466 0.068392552435398 0.074081249535084 0.076473690569401 0.075923353433609 0.072983093559742 0.067984402179718 0.062063418328762 0.055978287011385 0.049281783401966 0.042695127427578 0.039296623319387 0.040602579712868 0.041958440095186 0.041936922818422 0.041155569255352 0.040271371603012 0.039253354072571 0.037964228540659 0.036519519984722 0.035012576729059 0.033552922308445 0.032384194433689 0.031532708555460 0.031184166669846 0.031313393265009 0.031884945929050 0.032674264162779 0.033165600150824 0.033261474221945 0.032880887389183 0.032071515917778 0.030892191454768 0.029429819434881 0.027882147580385 0.026386935263872 0.025060936808586 0.023922393098474 0.023056661710143 0.022462941706181 0.022049317136407 0.021790260449052 0.021737795323133 0.021773239597678 0.021872024983168 0.022019866853952 0.022174360230565 0.022309415042400 0.022429853677750 0.022499836981297 0.022479727864265 0.022369915619493 0.022155307233334 0.021886080503464 0.021557629108429 0.021204160526395 0.020857909694314 0.020557988435030 0.020329663529992 0.020173605531454 0.020093526691198 0.020070081576705 0.020082853734493 0.020111048594117 0.020137432962656 0.020152274519205 0.020133748650551 0.020081480965018 0.019996410235763 0.019886409863830 0.019764542579651 0.019636582583189 0.019513620063663 0.019412502646446 0.019338848069310 0.019291685894132 0.019266514107585 0.019252140074968 0.019240213558078 0.019221086055040 0.019188165664673 0.019128002226353 0.019050288945436 0.018958097323775 0.018865045160055 0.018731493502855 0.018595099449158 0.018461331725121 0.018341148272157 0.018239950761199 0.018164556473494 0.018117444589734 0.018096150830388 0.018096923828125 0.018109502270818 0.018126402050257 0.018141582608223 0.018146086484194 0.018135767430067 0.018111884593964 0.018075661733747 0.018028963357210 0.017978558316827 0.017930492758751 0.017892396077514 0.017869470641017 0.017866028472781 0.017882086336613 0.017915511503816 0.017964184284210 0.018018390983343 0.018072009086609 0.018122663721442 0.018162231892347 0.018185926601291 0.018195742741227 0.018188903108239 0.018166121095419 0.018130375072360 0.018083171918988 0.018029093742371 0.017970109358430 0.017909364774823 0.017850009724498 0.017791517078876 0.017737455666065 0.017687913030386 0.017641842365265 0.017599614337087 0.017560971900821 0.017526511102915 0.017494561150670 0.017464602366090 0.017436491325498 0.017410986125469 0.017387712374330 0.017366759479046 0.017348365858197 0.017333002761006 0.017320757731795 0.017311546951532 0.017305411398411 0.017302019521594 0.017300834879279 0.017301265150309 0.017303397879004 0.017307009547949 0.017311235889792 0.017316164448857 0.017321202903986 0.017325842753053 0.017329482361674 0.017332492396235 0.017334889620543 0.017335424199700 0.017334438860416 0.017330650240183 0.017324367538095 0.017314461991191 0.017301702871919 0.017286710441113 0.017270233482122 0.017251797020435 0.017231697216630 0.017211109399796 0.017190055921674 0.017168415710330 0.017146700993180 0.017125427722931 0.017103662714362 0.017083011567593 0.017062736675143 0.017043512314558 0.017025090754032 0.017008036375046 0.016991944983602 0.016977254301310 0.016963746398687 0.016951307654381 0.016940213739872 0.016929993405938 0.016921287402511 0.016913559287786 0.016906872391701 0.016900980845094 0.016895970329642 0.016892237588763 0.016889268532395 0.016886493191123 0.016884177923203 0.016881842166185 0.016879713162780 0.016877787187696 0.016875529661775 0.016873074695468 0.016870390623808 0.016866931691766 0.016862638294697 0.016857827082276 0.016852591186762 0.016846535727382 0.016839319840074 0.016831945627928 0.016823654994369 0.016814803704619 0.016805408522487 0.016795495525002 0.016785522922873 0.016774974763393 0.016764082014561 0.016753070056438 0.016741814091802 0.016730703413486 0.016719888895750 0.016709184274077 0.016698973253369 0.016689024865627 0.016679516062140 0.016670590266585 0.016662230715156 0.016653902828693 0.016646092757583 0.016638554632664 0.016631569713354 0.016624981537461 0.016618479043245 0.016612309962511 0.016606179997325 0.016600281000137 0.016594719141722 0.016589369624853 0.016584536060691 0.016579698771238 0.016575030982494 0.016570374369621 0.016565563157201 0.016560869291425 0.016556125134230 0.016551408916712 0.016546580940485 0.016541842371225 0.016536811366677 0.016531681641936 0.016526222229004 0.016520688310266 0.016515063121915 0.016509322449565 0.016503706574440 0.016497986391187 0.016492303460836 0.016486490145326 0.016480604186654 0.016474889591336 0.016468843445182 0.016463097184896 0.016457319259644 0.016451491042972 0.016445523127913 0.016439592465758 0.016433419659734 0.016427069902420 0.016420839354396 0.016414821147919 0.016409033909440 0.016403453424573 0.016397939994931 0.016392543911934 0.016387309879065 0.016382185742259 0.016377067193389 0.016371967270970 0.016367010772228 0.016361815854907 0.016356831416488 0.016351763159037 0.016346842050552 0.016341922804713 0.016337119042873 0.016332427039742 0.016327898949385 0.016323205083609 0.016318948939443 0.016314605250955 0.016310581937432 0.016306815668941 0.016303012147546 0.016299473121762 0.016295989975333 0.016292395070195 0.016288688406348 0.016285117715597 0.016281481832266 0.016277670860291 0.016273882240057 0.016269866377115 0.016265725716949 0.016261536628008 0.016257083043456 0.016252728179097 0.016247963532805 0.016243131831288 0.016238098964095 0.016232717782259 0.016227303072810 0.016221595928073 0.016215557232499 0.016209872439504 0.016204100102186 0.016198545694351 0.016193026676774 0.016187347471714 0.016181878745556 0.016176613047719 0.016171377152205 0.016166225075722 0.016161266714334 0.016156293451786 0.016151366755366 0.016146607697010 0.016141951084137 0.016137462109327 0.016132876276970 0.016128450632095 0.016124004498124 0.016119720414281 0.016115464270115 0.016111385077238 0.016107317060232 0.016103472560644 0.016099685803056 0.016096178442240 0.016092766076326 0.016089441254735 0.016086123883724 0.016082951799035 0.016079643741250 0.016076220199466 0.016072759404778 0.016069319099188 0.016065686941147 0.016062036156654 0.016058241948485 0.016054656356573 0.016051093116403 0.016047500073910 0.016043934971094 0.016040218994021 0.016036149114370 0.016032090410590 0.016027878969908 0.016023453325033 0.016019077971578 0.016014453023672 0.016009762883186 0.016004838049412 0.015999918803573 0.015995012596250 0.015990206971765 0.015985697507858 0.015981394797564 0.015977207571268 0.015973214060068 0.015969092026353 0.015965240076184 0.015961563214660 0.015958040952682 0.015954585745931 0.015950968489051 0.015947245061398 0.015943422913551 0.015939682722092 0.015935907140374 0.015932293608785 0.015928735956550 0.015925448387861 0.015922183170915 0.015919141471386 0.015916200354695 0.015913339331746 0.015910709276795 0.015908038243651 0.015905493870378 0.015902912244201 0.015900351107121 0.015897706151009
Out[22]:
Text(0.5, 0, 'Iteration number')
In [ ]:
In [ ]: