Merge pull request #3 from vais-ral/cleandemos

Tidy up all astra demos except demoIP

4 files changed, 195 insertions, 162 deletions

diff --git a/Wrappers/Python/ccpi/astra/processors.py b/Wrappers/Python/ccpi/astra/processors.py
index 16c1f78..02c9070 100755..100644
--- a/Wrappers/Python/ccpi/astra/processors.py
+++ b/Wrappers/Python/ccpi/astra/processors.py
@@ -3,7 +3,7 @@ from ccpi.astra.utils import convert_geometry_to_astra
 import astra
 
 
-class AstraForwardProjector(DataSetProcessor):
+class AstraForwardProjector(DataProcessor):
     '''AstraForwardProjector
     
     Forward project ImageData to AcquisitionData using ASTRA proj_id.
@@ -25,7 +25,7 @@ class AstraForwardProjector(DataSetProcessor):
                   'device'  : device
                   }
         
-        #DataSetProcessor.__init__(self, **kwargs)
+        #DataProcessor.__init__(self, **kwargs)
         super(AstraForwardProjector, self).__init__(**kwargs)
         
         self.set_ImageGeometry(volume_geometry)
@@ -77,7 +77,7 @@ class AstraForwardProjector(DataSetProcessor):
         #return AcquisitionData(array=DATA, geometry=self.sinogram_geometry)
         return DATA
 
-class AstraBackProjector(DataSetProcessor):
+class AstraBackProjector(DataProcessor):
     '''AstraBackProjector
     
     Back project AcquisitionData to ImageData using ASTRA proj_id.
@@ -99,7 +99,7 @@ class AstraBackProjector(DataSetProcessor):
                   'device'  : device
                   }
         
-        #DataSetProcessor.__init__(self, **kwargs)
+        #DataProcessor.__init__(self, **kwargs)
         super(AstraBackProjector, self).__init__(**kwargs)
         
         self.set_ImageGeometry(volume_geometry)
@@ -202,4 +202,4 @@ class AstraBackProjectorMC(AstraBackProjector):
                     DATA.as_array()[k], 
                     self.proj_id)
             astra.data2d.delete(rec_id)
-        return IM
\ No newline at end of file
+        return IM
diff --git a/Wrappers/Python/wip/demo_sophiabeads.py b/Wrappers/Python/wip/demo_sophiabeads.py
index 33ce9ec..53a0bf8 100755
--- a/Wrappers/Python/wip/demo_sophiabeads.py
+++ b/Wrappers/Python/wip/demo_sophiabeads.py
@@ -1,4 +1,12 @@
 
+# This demo shows how to load a Nikon XTek micro-CT data set and reconstruct
+# the central slice using the CGLS method. The SophiaBeads dataset with 256 
+# projections is used as test data and can be obtained from here:
+# https://zenodo.org/record/16474
+# The filename with full path to the .xtekct file should be given as string 
+# input to XTEKReader to  load in the data.
+
+# Do all imports
 from ccpi.io.reader import XTEKReader
 import numpy as np
 import matplotlib.pyplot as plt
@@ -7,7 +15,7 @@ from ccpi.astra.ops import AstraProjectorSimple
 from ccpi.optimisation.algs import CGLS
 
 # Set up reader object and read the data
-datareader = XTEKReader("C:/Users/mbbssjj2/Documents/SophiaBeads_256_averaged/SophiaBeads_256_averaged.xtekct")
+datareader = XTEKReader("REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_256_averaged.xtekct")
 data = datareader.getAcquisitionData()
 
 # Extract central slice, scale and negative-log transform
@@ -56,7 +64,7 @@ Aop = AstraProjectorSimple(ig2d, ag2d,"gpu")
 # Set initial guess for CGLS reconstruction
 x_init = ImageData(np.zeros((N,N)),geometry=ig2d)
 
-# Run CGLS reconstruction
+# Run 50-iteration CGLS reconstruction
 num_iter = 50
 x, it, timing, criter = CGLS(Aop,data2d,num_iter,x_init)
 
diff --git a/Wrappers/Python/wip/simple_demo_astra.py b/Wrappers/Python/wip/simple_demo_astra.py
index 26f3ff4..b4f0c61 100755
--- a/Wrappers/Python/wip/simple_demo_astra.py
+++ b/Wrappers/Python/wip/simple_demo_astra.py
@@ -1,189 +1,219 @@
-#import sys
-#sys.path.append("..")
 
+# This demo illustrates how ASTRA 2D projectors can be used with
+# the modular optimisation framework. The demo sets up a 2D test case and 
+# demonstrates reconstruction using CGLS, as well as FISTA for least squares 
+# and 1-norm regularisation and FBPD for 1-norm and TV regularisation.
+
+# First make all imports
 from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry
 from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1 , TV2D
+from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D
 from ccpi.astra.ops import AstraProjectorSimple
 
 import numpy as np
 import matplotlib.pyplot as plt
 
-test_case = 1   # 1=parallel2D, 2=cone2D
+# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
+test_case = 1
 
-# Set up phantom
+# Set up phantom size NxN by creating ImageGeometry, initialising the 
+# ImageData object with this geometry and empty array and finally put some
+# data into its array, and display as image.
 N = 128
-
-
-vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
-Phantom = ImageData(geometry=vg)
+ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+Phantom = ImageData(geometry=ig)
 
 x = Phantom.as_array()
 x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
 x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
 
 plt.imshow(x)
+plt.title('Phantom image')
 plt.show()
 
-# Set up measurement geometry
-angles_num = 20; # angles number
+# Set up AcquisitionGeometry object to hold the parameters of the measurement
+# setup geometry: # Number of angles, the actual angles from 0 to 
+# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
+# pixel relative to an object pixel, the number of detector pixels, and the 
+# source-origin and origin-detector distance (here the origin-detector distance 
+# set to 0 to simulate a "virtual detector" with same detector pixel size as
+# object pixel size).
+angles_num = 20
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
 
 if test_case==1:
     angles = np.linspace(0,np.pi,angles_num,endpoint=False)
+    ag = AcquisitionGeometry('parallel',
+                             '2D',
+                             angles,
+                             det_num,det_w)
 elif test_case==2:
     angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
+    ag = AcquisitionGeometry('cone',
+                             '2D',
+                             angles,
+                             det_num,
+                             det_w,
+                             dist_source_center=SourceOrig, 
+                             dist_center_detector=OrigDetec)
 else:
     NotImplemented
 
-det_w = 1.0
-det_num = N
-SourceOrig = 200
-OrigDetec = 0
+# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
+# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
+Aop = AstraProjectorSimple(ig, ag, 'gpu')
 
-# Parallelbeam geometry test
-if test_case==1:
-    pg = AcquisitionGeometry('parallel',
-                          '2D',
-                          angles,
-                          det_num,det_w)
-elif test_case==2:
-    pg = AcquisitionGeometry('cone',
-                          '2D',
-                          angles,
-                          det_num,
-                          det_w,
-                          dist_source_center=SourceOrig, 
-                          dist_center_detector=OrigDetec)
-
-# ASTRA operator using volume and sinogram geometries
-Aop = AstraProjectorSimple(vg, pg, 'cpu')
-
-# Unused old astra projector without geometry
-# Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
-#                      OrigDetec, angles, 
-#                      N,'fanbeam','gpu') 
-
-# Try forward and backprojection
+# Forward and backprojection are available as methods direct and adjoint. Here 
+# generate test data b and do simple backprojection to obtain z.
 b = Aop.direct(Phantom)
-out2 = Aop.adjoint(b)
+z = Aop.adjoint(b)
 
 plt.imshow(b.array)
+plt.title('Simulated data')
 plt.show()
 
-plt.imshow(out2.array)
+plt.imshow(z.array)
+plt.title('Backprojected data')
 plt.show()
 
-# Create least squares object instance with projector and data.
-f = Norm2sq(Aop,b,c=0.5)
+# Using the test data b, different reconstruction methods can now be set up as
+# demonstrated in the rest of this file. In general all methods need an initial 
+# guess and some algorithm options to be set:
+x_init = ImageData(np.zeros(x.shape),geometry=ig)
+opt = {'tol': 1e-4, 'iter': 1000}
 
-# Initial guess
-x_init = ImageData(np.zeros(x.shape),geometry=vg)
+# First a CGLS reconstruction can be done:
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt)
+
+plt.imshow(x_CGLS.array)
+plt.title('CGLS')
+plt.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS criterion')
+plt.show()
+
+# CGLS solves the simple least-squares problem. The same problem can be solved 
+# by FISTA by setting up explicitly a least squares function object and using 
+# no regularisation:
+
+# Create least squares object instance with projector, test data and a constant 
+# coefficient of 0.5:
+f = Norm2sq(Aop,b,c=0.5)
 
 # Run FISTA for least squares without regularization
-x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None)
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None,opt)
 
 plt.imshow(x_fista0.array)
-plt.title('FISTA0')
+plt.title('FISTA Least squares')
 plt.show()
 
-# Now least squares plus 1-norm regularization
+plt.semilogy(criter0)
+plt.title('FISTA Least squares criterion')
+plt.show()
+
+# FISTA can also solve regularised forms by specifying a second function object
+# such as 1-norm regularisation with choice of regularisation parameter lam:
+
+# Create 1-norm function object
 lam = 0.1
 g0 = Norm1(lam)
 
 # Run FISTA for least squares plus 1-norm function.
-x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0)
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
 
 plt.imshow(x_fista1.array)
-plt.title('FISTA1')
+plt.title('FISTA Least squares plus 1-norm regularisation')
 plt.show()
 
 plt.semilogy(criter1)
+plt.title('FISTA Least squares plus 1-norm regularisation criterion')
 plt.show()
 
-# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
-opt = {'tol': 1e-4, 'iter': 100}
-x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
+# The least squares plus 1-norm regularisation problem can also be solved by 
+# other algorithms such as the Forward Backward Primal Dual algorithm. This
+# algorithm minimises the sum of three functions and the least squares and 
+# 1-norm functions should be given as the second and third function inputs. 
+# In this test case, this algorithm requires more iterations to converge, so
+# new options are specified.
+opt_FBPD = {'tol': 1e-4, 'iter': 7000}
+x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt_FBPD)
 
 plt.imshow(x_fbpd1.array)
-plt.title('FBPD1')
+plt.title('FBPD for least squares plus 1-norm regularisation')
 plt.show()
 
 plt.semilogy(criter_fbpd1)
+plt.title('FBPD for least squares plus 1-norm regularisation criterion')
 plt.show()
 
-# Now FBPD for least squares plus TV
+# The FBPD algorithm can also be used conveniently for TV regularisation:
+
+# Specify TV function object
 #lamtv = 1
 #gtv = TV2D(lamtv)
-
-#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt)
-
+#
+#x_fbpdtv,it_fbpdtv,timing_fbpdtv,criter_fbpdtv=FBPD(x_init,None,f,gtv,opt_FBPD)
+#
 #plt.imshow(x_fbpdtv.array)
 #plt.show()
-
+#
 #plt.semilogy(criter_fbpdtv)
 #plt.show()
 
 
-# Run CGLS, which should agree with the FISTA0
-x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt )
-
-plt.imshow(x_CGLS.array)
-plt.title('CGLS')
-#plt.title('CGLS recon, compare FISTA0')
-plt.show()
-
-plt.semilogy(criter_CGLS)
-plt.title('CGLS criterion')
-plt.show()
-
-
-#%%
-
+# Compare all reconstruction and criteria
 clims = (0,1)
 cols = 3
 rows = 2
 current = 1
+
 fig = plt.figure()
-# projections row
 a=fig.add_subplot(rows,cols,current)
 a.set_title('phantom {0}'.format(np.shape(Phantom.as_array())))
-
 imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1])
+plt.axis('off')
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA0')
-imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1])
+a.set_title('CGLS')
+imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1])
+plt.axis('off')
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA1')
-imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
+a.set_title('FISTA LS')
+imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1])
+plt.axis('off')
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FBPD1')
-imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1])
+a.set_title('FISTA LS+1')
+imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
+plt.axis('off')
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('CGLS')
-imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1])
+a.set_title('FBPD LS+1')
+imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1])
+plt.axis('off')
 
 #current = current + 1
 #a=fig.add_subplot(rows,cols,current)
 #a.set_title('FBPD TV')
 #imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1])
+#plt.axis('off')
 
 fig = plt.figure()
-# projections row
 b=fig.add_subplot(1,1,1)
 b.set_title('criteria')
-imgplot = plt.loglog(criter0 , label='FISTA0')
-imgplot = plt.loglog(criter1 , label='FISTA1')
-imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
 imgplot = plt.loglog(criter_CGLS, label='CGLS')
+imgplot = plt.loglog(criter0 , label='FISTA LS')
+imgplot = plt.loglog(criter1 , label='FISTA LS+1')
+imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1')
 #imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV')
-b.legend(loc='right')
+b.legend(loc='lower left')
 plt.show()
-#%%
\ No newline at end of file
diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py
index 7369d8f..3b78eb3 100755
--- a/Wrappers/Python/wip/simple_mc_demo.py
+++ b/Wrappers/Python/wip/simple_mc_demo.py
@@ -1,6 +1,10 @@
-#import sys
-#sys.path.append("..")
 
+# This demo demonstrates a simple multichannel reconstruction case. A 
+# synthetic 3-channel phantom image is set up, data is simulated and the FISTA 
+# algorithm is used to compute least squares and least squares with 1-norm 
+# regularisation reconstructions.
+
+# Do all imports
 from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry
 from ccpi.optimisation.algs import FISTA
 from ccpi.optimisation.funcs import Norm2sq, Norm1
@@ -9,15 +13,16 @@ from ccpi.astra.ops import AstraProjectorMC
 import numpy
 import matplotlib.pyplot as plt
 
-test_case = 1   # 1=parallel2D, 2=cone2D
+# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
+test_case = 1
 
-# Set up phantom
+# Set up phantom NxN pixels and 3 channels. Set up the ImageGeometry and fill 
+# some test data in each of the channels. Display each channel as image.
 N = 128
-M = 100
 numchannels = 3
 
-vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels)
-Phantom = ImageData(geometry=vg)
+ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels)
+Phantom = ImageData(geometry=ig)
 
 x = Phantom.as_array()
 x[0 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.0
@@ -29,65 +34,52 @@ x[1 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 1.2
 x[2 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.5
 x[2 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.2
 
-#x = numpy.zeros((N,N,1,numchannels))
-
-#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0
-#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0
-
-#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7
-#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2
-
-#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5
-#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2
-
 f, axarr = plt.subplots(1,numchannels)
 for k in numpy.arange(3):
     axarr[k].imshow(x[k],vmin=0,vmax=2.5)
 plt.show()
 
-#vg = ImageGeometry(N,N,None, 1,1,None,channels=numchannels)
-
-
-#Phantom = ImageData(x,geometry=vg)
-
-# Set up measurement geometry
-angles_num = 20; # angles number
-
-if test_case==1:
-    angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False)
-elif test_case==2:
-    angles = numpy.linspace(0,2*numpy.pi,angles_num,endpoint=False)
-else:
-    NotImplemented
-
+# Set up AcquisitionGeometry object to hold the parameters of the measurement
+# setup geometry: # Number of angles, the actual angles from 0 to 
+# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
+# pixel relative to an object pixel, the number of detector pixels, and the 
+# source-origin and origin-detector distance (here the origin-detector distance 
+# set to 0 to simulate a "virtual detector" with same detector pixel size as
+# object pixel size).
+angles_num = 20
 det_w = 1.0
 det_num = N
 SourceOrig = 200
 OrigDetec = 0
 
-# Parallelbeam geometry test
 if test_case==1:
-    pg = AcquisitionGeometry('parallel',
-                          '2D',
-                          angles,
-                          det_num,
-                          det_w,
-                          channels=numchannels)
+    angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False)
+    ag = AcquisitionGeometry('parallel',
+                             '2D',
+                             angles,
+                             det_num,
+                             det_w,
+                             channels=numchannels)
 elif test_case==2:
-    pg = AcquisitionGeometry('cone',
-                          '2D',
-                          angles,
-                          det_num,
-                          det_w,
-                          dist_source_center=SourceOrig, 
-                          dist_center_detector=OrigDetec,
-                          channels=numchannels)
-
-# ASTRA operator using volume and sinogram geometries
-Aop = AstraProjectorMC(vg, pg, 'cpu')
+    angles = numpy.linspace(0,2*numpy.pi,angles_num,endpoint=False)
+    ag = AcquisitionGeometry('cone',
+                            '2D',
+                            angles,
+                            det_num,
+                            det_w,
+                            dist_source_center=SourceOrig, 
+                            dist_center_detector=OrigDetec,
+                            channels=numchannels)
+else:
+    NotImplemented
 
+# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
+# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
+Aop = AstraProjectorMC(ig, ag, 'gpu')
 
-# Try forward and backprojection
+# Forward and backprojection are available as methods direct and adjoint. Here 
+# generate test data b and do simple backprojection to obtain z. Applies 
+# channel by channel
 b = Aop.direct(Phantom)
 
 fb, axarrb = plt.subplots(1,numchannels)
@@ -95,26 +87,27 @@ for k in numpy.arange(3):
     axarrb[k].imshow(b.as_array()[k],vmin=0,vmax=250)
 plt.show()
 
-out2 = Aop.adjoint(b)
+z = Aop.adjoint(b)
 
 fo, axarro = plt.subplots(1,numchannels)
 for k in range(3):
-    axarro[k].imshow(out2.as_array()[k],vmin=0,vmax=3500)
+    axarro[k].imshow(z.as_array()[k],vmin=0,vmax=3500)
 plt.show()
 
-# Create least squares object instance with projector and data.
-f = Norm2sq(Aop,b,c=0.5)
-
-# Initial guess
-x_init = ImageData(numpy.zeros(x.shape),geometry=vg)
-
-# FISTA options
+# Using the test data b, different reconstruction methods can now be set up as
+# demonstrated in the rest of this file. In general all methods need an initial 
+# guess and some algorithm options to be set:
+x_init = ImageData(np.zeros(x.shape),geometry=ig)
 opt = {'tol': 1e-4, 'iter': 200}
 
+# Create least squares object instance with projector, test data and a constant 
+# coefficient of 0.5. Note it is least squares over all channels:
+f = Norm2sq(Aop,b,c=0.5)
+
 # Run FISTA for least squares without regularization
 x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
 
-
+# Display reconstruction and criteration
 ff0, axarrf0 = plt.subplots(1,numchannels)
 for k in numpy.arange(3):
     axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5)
@@ -124,14 +117,16 @@ plt.semilogy(criter0)
 plt.title('Criterion vs iterations, least squares')
 plt.show()
 
-# Now least squares plus 1-norm regularization
+# FISTA can also solve regularised forms by specifying a second function object
+# such as 1-norm regularisation with choice of regularisation parameter lam. 
+# Again the regulariser is over all channels:
 lam = 0.1
 g0 = Norm1(lam)
 
-
 # Run FISTA for least squares plus 1-norm function.
 x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
 
+# Display reconstruction and criteration
 ff1, axarrf1 = plt.subplots(1,numchannels)
 for k in numpy.arange(3):
     axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5)
@@ -139,4 +134,4 @@ plt.show()
 
 plt.semilogy(criter1)
 plt.title('Criterion vs iterations, least squares plus 1-norm regu')
-plt.show()
\ No newline at end of file
+plt.show()