Merge remote-tracking branch 'origin/composite_operator_datacontainer' into demo_ccpi

8 files changed, 482 insertions, 82 deletions

diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index 7516447..7236e0e 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -707,6 +707,10 @@ class DataContainer(object):
     def maximum(self, x2, *args, **kwargs):
         return self.pixel_wise_binary(numpy.maximum, x2, *args, **kwargs)
     
+    def minimum(self,x2, out=None, *args, **kwargs):
+        return self.pixel_wise_binary(numpy.minimum, x2=x2, out=out, *args, **kwargs)
+
+    
     ## unary operations
     def pixel_wise_unary(self, pwop, *args,  **kwargs):
         out = kwargs.get('out', None)
@@ -760,12 +764,26 @@ class DataContainer(object):
         return numpy.sqrt(self.squared_norm())
     def dot(self, other, *args, **kwargs):
         '''return the inner product of 2 DataContainers viewed as vectors'''
+        method = kwargs.get('method', 'reduce')
+        if method not in ['numpy','reduce']:
+            raise ValueError('dot: specified method not valid. Expecting numpy or reduce got {} '.format(
+                    method))
         if self.shape == other.shape:
-            return numpy.dot(self.as_array().ravel(), other.as_array().ravel())
+            # return (self*other).sum()
+            if method == 'numpy':
+                return numpy.dot(self.as_array().ravel(), other.as_array())
+            elif method == 'reduce':
+                # see https://github.com/vais-ral/CCPi-Framework/pull/273
+                # notice that Python seems to be smart enough to use
+                # the appropriate type to hold the result of the reduction
+                sf = reduce(lambda x,y: x + y[0]*y[1],
+                            zip(self.as_array().ravel(),
+                                other.as_array().ravel()),
+                            0)
+                return sf
         else:
             raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape))
-    
-    
+   
     
     
     
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
new file mode 100644
index 0000000..30584d4
--- /dev/null
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2018 Edoardo Pasca
+
+#   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.
+
+from ccpi.optimisation.algorithms import Algorithm
+
+class SIRT(Algorithm):
+
+    '''Simultaneous Iterative Reconstruction Technique
+
+    Parameters:
+      x_init: initial guess
+      operator: operator for forward/backward projections
+      data: data to operate on
+      constraint: Function with prox-method, for example IndicatorBox to 
+                  enforce box constraints, default is None).
+    '''
+    def __init__(self, **kwargs):
+        super(SIRT, self).__init__()
+        self.x          = kwargs.get('x_init', None)
+        self.operator   = kwargs.get('operator', None)
+        self.data       = kwargs.get('data', None)
+        self.constraint = kwargs.get('constraint', None)
+        if self.x is not None and self.operator is not None and \
+           self.data is not None:
+            print ("Calling from creator")
+            self.set_up(x_init=kwargs['x_init'],
+                        operator=kwargs['operator'],
+                        data=kwargs['data'],
+                        constraint=kwargs['constraint'])
+
+    def set_up(self, x_init, operator , data, constraint=None ):
+        
+        self.x = x_init.copy()
+        self.operator = operator
+        self.data = data
+        self.constraint = constraint
+        
+        self.r = data.copy()
+        
+        self.relax_par = 1.0
+        
+        # Set up scaling matrices D and M.
+        self.M = 1/self.operator.direct(self.operator.domain_geometry().allocate(value=1.0))        
+        self.D = 1/self.operator.adjoint(self.operator.range_geometry().allocate(value=1.0))
+
+
+    def update(self):
+        
+        self.r = self.data - self.operator.direct(self.x)
+        
+        self.x += self.relax_par * (self.D*self.operator.adjoint(self.M*self.r))
+        
+        if self.constraint != None:
+            self.x = self.constraint.prox(self.x,None)
+
+    def update_objective(self):
+        self.loss.append(self.r.squared_norm())
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
index f562973..2dbacfc 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
@@ -24,6 +24,7 @@ Created on Thu Feb 21 11:03:13 2019
 
 from .Algorithm import Algorithm
 from .CGLS import CGLS
+from .SIRT import SIRT
 from .GradientDescent import GradientDescent
 from .FISTA import FISTA
 from .FBPD import FBPD
diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py
index 2f819d3..f5ba85e 100755
--- a/Wrappers/Python/ccpi/optimisation/algs.py
+++ b/Wrappers/Python/ccpi/optimisation/algs.py
@@ -20,13 +20,8 @@
 import numpy
 import time
 
-from ccpi.optimisation.functions import Function
-from ccpi.optimisation.functions import ZeroFunction
-from ccpi.framework import ImageData 
-from ccpi.framework import AcquisitionData
-from ccpi.optimisation.spdhg import spdhg 
-from ccpi.optimisation.spdhg import KullbackLeibler
-from ccpi.optimisation.spdhg import KullbackLeiblerConvexConjugate
+
+
 
 def FISTA(x_init, f=None, g=None, opt=None):
     '''Fast Iterative Shrinkage-Thresholding Algorithm
@@ -280,10 +275,6 @@ def SIRT(x_init, operator , data , opt=None, constraint=None):
     tol = opt['tol']
     max_iter = opt['iter']
     
-    # Set default constraint to unconstrained
-    if constraint==None:
-        constraint = Function()
-    
     x = x_init.clone()
     
     timing = numpy.zeros(max_iter)
@@ -293,21 +284,18 @@ def SIRT(x_init, operator , data , opt=None, constraint=None):
     relax_par = 1.0
     
     # Set up scaling matrices D and M.
-    im1 = ImageData(geometry=x_init.geometry)
-    im1.array[:] = 1.0
-    M = 1/operator.direct(im1)
-    del im1
-    aq1 = AcquisitionData(geometry=M.geometry)
-    aq1.array[:] = 1.0
-    D = 1/operator.adjoint(aq1)
-    del aq1
+    M = 1/operator.direct(operator.domain_geometry().allocate(value=1.0))
+    D = 1/operator.adjoint(operator.range_geometry().allocate(value=1.0))
     
     # algorithm loop
     for it in range(0, max_iter):
         t = time.time()
         r = data - operator.direct(x)
         
-        x = constraint.prox(x + relax_par * (D*operator.adjoint(M*r)),None)
+        x = x + relax_par * (D*operator.adjoint(M*r))
+
+        if constraint != None:
+            x = constraint.prox(x,None)
         
         timing[it] = time.time() - t
         if it > 0:
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 8e8ab87..e92d4c6 100755
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -455,6 +455,11 @@ class TestDataContainer(unittest.TestCase):
             self.assertTrue(False)
         except ValueError as ve:
             self.assertTrue(True)
+            
+        print ("test dot numpy")
+        n0 = (ds0 * ds1).sum()
+        n1 = ds0.as_array().ravel().dot(ds1.as_array().ravel())
+        self.assertEqual(n0, n1)
         
         
 
diff --git a/Wrappers/Python/wip/compare_CGLS_algos.py b/Wrappers/Python/wip/compare_CGLS_algos.py
new file mode 100644
index 0000000..119752c
--- /dev/null
+++ b/Wrappers/Python/wip/compare_CGLS_algos.py
@@ -0,0 +1,127 @@
+# This demo illustrates how to use the SIRT algorithm without and with 
+# nonnegativity and box constraints. The ASTRA 2D projectors are used.
+
+# First make all imports
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \
+    AcquisitionData
+from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT
+from ccpi.astra.operators import AstraProjectorSimple
+
+from ccpi.optimisation.algorithms import CGLS as CGLSalg
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
+test_case = 1
+
+# 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
+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.figure()
+#plt.imshow(x)
+#plt.title('Phantom image')
+#plt.show()
+
+# 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
+
+# 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, 'cpu')
+
+# 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)
+z = Aop.adjoint(b)
+
+#plt.figure()
+#plt.imshow(b.array)
+#plt.title('Simulated data')
+#plt.show()
+
+#plt.figure()
+#plt.imshow(z.array)
+#plt.title('Backprojected data')
+#plt.show()
+
+# 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': 7}
+
+# First a CGLS reconstruction using the function version of CGLS can be done:
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt)
+
+#plt.figure()
+#plt.imshow(x_CGLS.array)
+#plt.title('CGLS')
+#plt.colorbar()
+#plt.show()
+
+#plt.figure()
+#plt.semilogy(criter_CGLS)
+#plt.title('CGLS criterion')
+#plt.show()
+
+
+
+# Now CLGS using the algorithm class
+CGLS_alg = CGLSalg()
+CGLS_alg.set_up(x_init, Aop, b )
+CGLS_alg.max_iteration = 2000
+CGLS_alg.run(opt['iter'])
+x_CGLS_alg = CGLS_alg.get_output()
+
+#plt.figure()
+#plt.imshow(x_CGLS_alg.as_array())
+#plt.title('CGLS ALG')
+#plt.colorbar()
+#plt.show()
+
+#plt.figure()
+#plt.semilogy(CGLS_alg.objective)
+#plt.title('CGLS criterion')
+#plt.show()
+
+print(criter_CGLS)
+print(CGLS_alg.objective)
+
+print((x_CGLS - x_CGLS_alg).norm())
\ No newline at end of file
diff --git a/Wrappers/Python/wip/demo_test_sirt.py b/Wrappers/Python/wip/demo_SIRT.py
index 6f5a44d..5a85d41 100644
--- a/Wrappers/Python/wip/demo_test_sirt.py
+++ b/Wrappers/Python/wip/demo_SIRT.py
@@ -4,9 +4,9 @@
 # First make all imports
 from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \
     AcquisitionData
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT
-from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D, IndicatorBox
+from ccpi.optimisation.functions import IndicatorBox
 from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.optimisation.algorithms import SIRT, CGLS
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -25,6 +25,7 @@ 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.figure()
 plt.imshow(x)
 plt.title('Phantom image')
 plt.show()
@@ -69,11 +70,13 @@ Aop = AstraProjectorSimple(ig, ag, 'gpu')
 b = Aop.direct(Phantom)
 z = Aop.adjoint(b)
 
-plt.imshow(b.array)
+plt.figure()
+plt.imshow(b.as_array())
 plt.title('Simulated data')
 plt.show()
 
-plt.imshow(z.array)
+plt.figure()
+plt.imshow(z.as_array())
 plt.title('Backprojected data')
 plt.show()
 
@@ -81,96 +84,122 @@ plt.show()
 # 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}
+opt = {'tol': 1e-4, 'iter': 100}
 
-# 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')
+# First run a simple CGLS reconstruction:
+CGLS_alg = CGLS()
+CGLS_alg.set_up(x_init, Aop, b )
+CGLS_alg.max_iteration = 2000
+CGLS_alg.run(opt['iter'])
+x_CGLS_alg = CGLS_alg.get_output()
+
+plt.figure()
+plt.imshow(x_CGLS_alg.as_array())
+plt.title('CGLS ALG')
 plt.colorbar()
 plt.show()
 
-plt.semilogy(criter_CGLS)
+plt.figure()
+plt.semilogy(CGLS_alg.objective)
 plt.title('CGLS criterion')
 plt.show()
 
-# A SIRT unconstrained reconstruction can be done: similarly:
-x_SIRT, it_SIRT, timing_SIRT, criter_SIRT = SIRT(x_init, Aop, b, opt)
 
-plt.imshow(x_SIRT.array)
+# A SIRT reconstruction can be done simply by replacing CGLS by SIRT.
+# In the first instance, no constraints are enforced.
+SIRT_alg = SIRT()
+SIRT_alg.set_up(x_init, Aop, b )
+SIRT_alg.max_iteration = 2000
+SIRT_alg.run(opt['iter'])
+x_SIRT_alg = SIRT_alg.get_output()
+
+plt.figure()
+plt.imshow(x_SIRT_alg.as_array())
 plt.title('SIRT unconstrained')
 plt.colorbar()
 plt.show()
 
-plt.semilogy(criter_SIRT)
+plt.figure()
+plt.semilogy(SIRT_alg.objective)
 plt.title('SIRT unconstrained criterion')
 plt.show()
 
-# A SIRT nonnegativity constrained reconstruction can be done using the 
-# additional input "constraint" set to a box indicator function with 0 as the 
-# lower bound and the default upper bound of infinity:
-x_SIRT0, it_SIRT0, timing_SIRT0, criter_SIRT0 = SIRT(x_init, Aop, b, opt,
-                                                      constraint=IndicatorBox(lower=0))
+# The SIRT algorithm is stopped after the specified number of iterations has 
+# been run. It can be resumed by calling the run command again, which will run 
+# it for the specificed number of iterations
+SIRT_alg.run(opt['iter'])
+x_SIRT_alg2 = SIRT_alg.get_output()
 
-plt.imshow(x_SIRT0.array)
-plt.title('SIRT nonneg')
+plt.figure()
+plt.imshow(x_SIRT_alg2.as_array())
+plt.title('SIRT unconstrained, extra iterations')
 plt.colorbar()
 plt.show()
 
-plt.semilogy(criter_SIRT0)
-plt.title('SIRT nonneg criterion')
+plt.figure()
+plt.semilogy(SIRT_alg.objective)
+plt.title('SIRT unconstrained criterion, extra iterations')
 plt.show()
 
-# A SIRT reconstruction with box constraints on [0,1] can also be done:
-x_SIRT01, it_SIRT01, timing_SIRT01, criter_SIRT01 = SIRT(x_init, Aop, b, opt,
-         constraint=IndicatorBox(lower=0,upper=1))
 
-plt.imshow(x_SIRT01.array)
-plt.title('SIRT box(0,1)')
+# A SIRT nonnegativity constrained reconstruction can be done using the 
+# additional input "constraint" set to a box indicator function with 0 as the 
+# lower bound and the default upper bound of infinity. First setup a new 
+# instance of the SIRT algorithm.
+SIRT_alg0 = SIRT()
+SIRT_alg0.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0) )
+SIRT_alg0.max_iteration = 2000
+SIRT_alg0.run(opt['iter'])
+x_SIRT_alg0 = SIRT_alg0.get_output()
+
+plt.figure()
+plt.imshow(x_SIRT_alg0.as_array())
+plt.title('SIRT nonnegativity constrained')
 plt.colorbar()
 plt.show()
 
-plt.semilogy(criter_SIRT01)
-plt.title('SIRT box(0,1) criterion')
-plt.show()
-
-# The indicator function can also be used with the FISTA algorithm to do 
-# least squares with nonnegativity constraint.
-
-# 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 constraints
-x_fista, it, timing, criter = FISTA(x_init, f, None,opt)
-
-plt.imshow(x_fista.array)
-plt.title('FISTA Least squares')
+plt.figure()
+plt.semilogy(SIRT_alg0.objective)
+plt.title('SIRT nonnegativity criterion')
 plt.show()
 
-plt.semilogy(criter)
-plt.title('FISTA Least squares criterion')
-plt.show()
 
-# Run FISTA for least squares with nonnegativity constraint
-x_fista0, it0, timing0, criter0 = FISTA(x_init, f, IndicatorBox(lower=0),opt)
+# A SIRT reconstruction with box constraints on [0,1] can also be done. 
+SIRT_alg01 = SIRT()
+SIRT_alg01.set_up(x_init, Aop, b, constraint=IndicatorBox(lower=0,upper=1) )
+SIRT_alg01.max_iteration = 2000
+SIRT_alg01.run(opt['iter'])
+x_SIRT_alg01 = SIRT_alg01.get_output()
 
-plt.imshow(x_fista0.array)
-plt.title('FISTA Least squares nonneg')
+plt.figure()
+plt.imshow(x_SIRT_alg01.as_array())
+plt.title('SIRT boc(0,1)')
+plt.colorbar()
 plt.show()
 
-plt.semilogy(criter0)
-plt.title('FISTA Least squares nonneg criterion')
+plt.figure()
+plt.semilogy(SIRT_alg01.objective)
+plt.title('SIRT box(0,1) criterion')
 plt.show()
 
-# Run FISTA for least squares with box constraint [0,1]
-x_fista01, it01, timing01, criter01 = FISTA(x_init, f, IndicatorBox(lower=0,upper=1),opt)
+# The test image has values in the range [0,1], so enforcing values in the 
+# reconstruction to be within this interval improves a lot. Just for fun
+# we can also easily see what happens if we choose a narrower interval as 
+# constrint in the reconstruction, lower bound 0.2, upper bound 0.8. 
+SIRT_alg0208 = SIRT()
+SIRT_alg0208.set_up(x_init,Aop,b,constraint=IndicatorBox(lower=0.2,upper=0.8))
+SIRT_alg0208.max_iteration = 2000
+SIRT_alg0208.run(opt['iter'])
+x_SIRT_alg0208 = SIRT_alg0208.get_output()
 
-plt.imshow(x_fista01.array)
-plt.title('FISTA Least squares box(0,1)')
+plt.figure()
+plt.imshow(x_SIRT_alg0208.as_array())
+plt.title('SIRT boc(0.2,0.8)')
+plt.colorbar()
 plt.show()
 
-plt.semilogy(criter01)
-plt.title('FISTA Least squares box(0,1) criterion')
+plt.figure()
+plt.semilogy(SIRT_alg0208.objective)
+plt.title('SIRT box(0.2,0.8) criterion')
 plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/demo_box_constraints_FISTA.py b/Wrappers/Python/wip/demo_box_constraints_FISTA.py
new file mode 100644
index 0000000..2f9e0c6
--- /dev/null
+++ b/Wrappers/Python/wip/demo_box_constraints_FISTA.py
@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 17 14:46:21 2019
+
+@author: jakob
+
+Demonstrate the use of box constraints in FISTA
+"""
+
+# First make all imports
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \
+    AcquisitionData
+from ccpi.optimisation.algorithms import FISTA
+from ccpi.optimisation.functions import Norm2sq, IndicatorBox
+from ccpi.astra.ops import AstraProjectorSimple
+
+from ccpi.optimisation.operators import Identity
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# 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
+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.figure()
+plt.imshow(x)
+plt.title('Phantom image')
+plt.show()
+
+# 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
+
+test_case = 1
+
+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
+
+# 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')
+
+Aop = Identity(ig,ig)
+
+# 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)
+z = Aop.adjoint(b)
+
+plt.figure()
+plt.imshow(b.array)
+plt.title('Simulated data')
+plt.show()
+
+plt.figure()
+plt.imshow(z.array)
+plt.title('Backprojected data')
+plt.show()
+
+# 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': 100}
+
+
+
+# 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 constraints
+FISTA_alg = FISTA()
+FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
+FISTA_alg.max_iteration = 2000
+FISTA_alg.run(opt['iter'])
+x_FISTA = FISTA_alg.get_output()
+
+plt.figure()
+plt.imshow(x_FISTA.array)
+plt.title('FISTA unconstrained')
+plt.colorbar()
+plt.show()
+
+plt.figure()
+plt.semilogy(FISTA_alg.objective)
+plt.title('FISTA unconstrained criterion')
+plt.show()
+
+# Run FISTA for least squares with lower bound 0.1
+FISTA_alg0 = FISTA()
+FISTA_alg0.set_up(x_init=x_init, f=f, g=IndicatorBox(lower=0.1), opt=opt)
+FISTA_alg0.max_iteration = 2000
+FISTA_alg0.run(opt['iter'])
+x_FISTA0 = FISTA_alg0.get_output()
+
+plt.figure()
+plt.imshow(x_FISTA0.array)
+plt.title('FISTA lower bound 0.1')
+plt.colorbar()
+plt.show()
+
+plt.figure()
+plt.semilogy(FISTA_alg0.objective)
+plt.title('FISTA criterion, lower bound 0.1')
+plt.show()
+
+# Run FISTA for least squares with box constraint [0.1,0.8]
+FISTA_alg0 = FISTA()
+FISTA_alg0.set_up(x_init=x_init, f=f, g=IndicatorBox(lower=0.1,upper=0.8), opt=opt)
+FISTA_alg0.max_iteration = 2000
+FISTA_alg0.run(opt['iter'])
+x_FISTA0 = FISTA_alg0.get_output()
+
+plt.figure()
+plt.imshow(x_FISTA0.array)
+plt.title('FISTA box(0.1,0.8) constrained')
+plt.colorbar()
+plt.show()
+
+plt.figure()
+plt.semilogy(FISTA_alg0.objective)
+plt.title('FISTA criterion, box(0.1,0.8) constrained criterion')
+plt.show()
\ No newline at end of file