new dir old demos

8 files changed, 1703 insertions, 0 deletions

diff --git a/Wrappers/Python/wip/old_demos/demo_colourbay.py b/Wrappers/Python/wip/old_demos/demo_colourbay.py
new file mode 100644
index 0000000..5dbf2e1
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/demo_colourbay.py
@@ -0,0 +1,137 @@
+# This script demonstrates how to load a mat-file with UoM colour-bay data 
+# into the CIL optimisation framework and run (simple) multichannel 
+# reconstruction methods.
+
+# All third-party imports.
+import numpy
+from scipy.io import loadmat
+import matplotlib.pyplot as plt
+
+# All own imports.
+from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData
+from ccpi.astra.ops import AstraProjectorMC
+from ccpi.optimisation.algs import CGLS, FISTA
+from ccpi.optimisation.funcs import Norm2sq, Norm1
+
+# Load full data and permute to expected ordering. Change path as necessary.
+# The loaded X has dims 80x60x80x150, which is pix x angle x pix x channel.
+# Permute (numpy.transpose) puts into our default ordering which is 
+# (channel, angle, vertical, horizontal).
+
+pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/ColourBay/spectral_data_sets/CarbonPd/'
+filename = 'carbonPd_full_sinogram_stripes_removed.mat'
+
+X = loadmat(pathname + filename)
+X = numpy.transpose(X['SS'],(3,1,2,0))
+
+# Store geometric variables for reuse
+num_channels = X.shape[0]
+num_pixels_h = X.shape[3]
+num_pixels_v = X.shape[2]
+num_angles = X.shape[1]
+
+# Display a single projection in a single channel
+plt.imshow(X[100,5,:,:])
+plt.title('Example of a projection image in one channel' )
+plt.show()
+
+# Set angles to use
+angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False)
+
+# Define full 3D acquisition geometry and data container.
+# Geometric info is taken from the txt-file in the same dir as the mat-file
+ag = AcquisitionGeometry('cone',
+                         '3D',
+                         angles,
+                         pixel_num_h=num_pixels_h,
+                         pixel_size_h=0.25,
+                         pixel_num_v=num_pixels_v,
+                         pixel_size_v=0.25,                            
+                         dist_source_center=233.0, 
+                         dist_center_detector=245.0,
+                         channels=num_channels)
+data = AcquisitionData(X, geometry=ag)
+
+# Reduce to central slice by extracting relevant parameters from data and its
+# geometry. Perhaps create function to extract central slice automatically?
+data2d = data.subset(vertical=40)
+ag2d = AcquisitionGeometry('cone',
+                         '2D',
+                         ag.angles,
+                         pixel_num_h=ag.pixel_num_h,
+                         pixel_size_h=ag.pixel_size_h,
+                         pixel_num_v=1,
+                         pixel_size_v=ag.pixel_size_h,                            
+                         dist_source_center=ag.dist_source_center, 
+                         dist_center_detector=ag.dist_center_detector,
+                         channels=ag.channels)
+data2d.geometry = ag2d
+
+# Set up 2D Image Geometry.
+# First need the geometric magnification to scale the voxel size relative
+# to the detector pixel size.
+mag = (ag.dist_source_center + ag.dist_center_detector)/ag.dist_source_center
+ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, 
+                     voxel_num_y=ag2d.pixel_num_h,  
+                     voxel_size_x=ag2d.pixel_size_h/mag, 
+                     voxel_size_y=ag2d.pixel_size_h/mag, 
+                     channels=X.shape[0])
+
+# Create GPU multichannel projector/backprojector operator with ASTRA.
+Aall = AstraProjectorMC(ig2d,ag2d,'gpu')
+
+# Compute and simple backprojction and display one channel as image.
+Xbp = Aall.adjoint(data2d)
+plt.imshow(Xbp.subset(channel=100).array)
+plt.show()
+
+# Set initial guess ImageData with zeros for algorithms, and algorithm options.
+x_init = ImageData(numpy.zeros((num_channels,num_pixels_v,num_pixels_h)),
+                   geometry=ig2d,
+                   dimension_labels=['channel','horizontal_y','horizontal_x'])
+opt_CGLS = {'tol': 1e-4, 'iter': 5}
+
+# Run CGLS algorithm and display one channel.
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aall, data2d, opt_CGLS)
+
+plt.imshow(x_CGLS.subset(channel=100).array)
+plt.title('CGLS')
+plt.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS Criterion vs iterations')
+plt.show()
+
+# 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(Aall,data2d,c=0.5)
+
+# Options for FISTA algorithm.
+opt = {'tol': 1e-4, 'iter': 100}
+
+# Run FISTA for least squares without regularization and display one channel
+# reconstruction as image.
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
+
+plt.imshow(x_fista0.subset(channel=100).array)
+plt.title('FISTA LS')
+plt.show()
+
+plt.semilogy(criter0)
+plt.title('FISTA LS Criterion vs iterations')
+plt.show()
+
+# Set up 1-norm regularisation (over all channels), solve with FISTA, and 
+# display one channel of reconstruction.
+lam = 0.1
+g0 = Norm1(lam)
+
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
+
+plt.imshow(x_fista1.subset(channel=100).array)
+plt.title('FISTA LS+1')
+plt.show()
+
+plt.semilogy(criter1)
+plt.title('FISTA LS+1 Criterion vs iterations')
+plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/old_demos/demo_compare_cvx.py b/Wrappers/Python/wip/old_demos/demo_compare_cvx.py
new file mode 100644
index 0000000..27b1c97
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/demo_compare_cvx.py
@@ -0,0 +1,306 @@
+
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer
+from ccpi.optimisation.algs import FISTA, FBPD, CGLS
+from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D, Norm2
+
+from ccpi.optimisation.ops import LinearOperatorMatrix, TomoIdentity
+from ccpi.optimisation.ops import Identity
+from ccpi.optimisation.ops import FiniteDiff2D
+
+# Requires CVXPY, see http://www.cvxpy.org/
+# CVXPY can be installed in anaconda using
+# conda install -c cvxgrp cvxpy libgcc
+
+# Whether to use or omit CVXPY
+use_cvxpy = True
+if use_cvxpy:
+    from cvxpy import *
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Problem data.
+m = 30
+n = 20
+np.random.seed(1)
+Amat = np.random.randn(m, n)
+A = LinearOperatorMatrix(Amat)
+bmat = np.random.randn(m)
+bmat.shape = (bmat.shape[0],1)
+
+# A = Identity()
+# Change n to equal to m.
+
+b = DataContainer(bmat)
+
+# Regularization parameter
+lam = 10
+opt = {'memopt':True}
+# Create object instances with the test data A and b.
+f = Norm2sq(A,b,c=0.5, memopt=True)
+g0 = ZeroFun()
+
+# Initial guess
+x_init = DataContainer(np.zeros((n,1)))
+
+f.grad(x_init)
+
+# Run FISTA for least squares plus zero function.
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0 , opt=opt)
+
+# Print solution and final objective/criterion value for comparison
+print("FISTA least squares plus zero function solution and objective value:")
+print(x_fista0.array)
+print(criter0[-1])
+
+if use_cvxpy:
+    # Compare to CVXPY
+    
+    # Construct the problem.
+    x0 = Variable(n)
+    objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]) )
+    prob0 = Problem(objective0)
+    
+    # The optimal objective is returned by prob.solve().
+    result0 = prob0.solve(verbose=False,solver=SCS,eps=1e-9)
+    
+    # The optimal solution for x is stored in x.value and optimal objective value 
+    # is in result as well as in objective.value
+    print("CVXPY least squares plus zero function solution and objective value:")
+    print(x0.value)
+    print(objective0.value)
+
+# Plot criterion curve to see FISTA converge to same value as CVX.
+iternum = np.arange(1,1001)
+plt.figure()
+plt.loglog(iternum[[0,-1]],[objective0.value, objective0.value], label='CVX LS')
+plt.loglog(iternum,criter0,label='FISTA LS')
+plt.legend()
+plt.show()
+
+# Create 1-norm object instance
+g1 = Norm1(lam)
+
+g1(x_init)
+x_rand = DataContainer(np.reshape(np.random.rand(n),(n,1)))
+x_rand2 = DataContainer(np.reshape(np.random.rand(n-1),(n-1,1)))
+v = g1.prox(x_rand,0.02)
+#vv = g1.prox(x_rand2,0.02)
+vv = v.copy() 
+vv *= 0
+print (">>>>>>>>>>vv" , vv.as_array())
+vv.fill(v)
+print (">>>>>>>>>>fill" , vv.as_array())
+g1.proximal(x_rand, 0.02, out=vv)
+print (">>>>>>>>>>v" , v.as_array())
+print (">>>>>>>>>>gradient" , vv.as_array())
+
+print (">>>>>>>>>>" , (v-vv).as_array())
+import sys
+#sys.exit(0)
+# Combine with least squares and solve using generic FISTA implementation
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1,opt=opt)
+
+# Print for comparison
+print("FISTA least squares plus 1-norm solution and objective value:")
+print(x_fista1)
+print(criter1[-1])
+
+if use_cvxpy:
+    # Compare to CVXPY
+    
+    # Construct the problem.
+    x1 = Variable(n)
+    objective1 = Minimize(0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1,1) )
+    prob1 = Problem(objective1)
+    
+    # The optimal objective is returned by prob.solve().
+    result1 = prob1.solve(verbose=False,solver=SCS,eps=1e-9)
+    
+    # The optimal solution for x is stored in x.value and optimal objective value 
+    # is in result as well as in objective.value
+    print("CVXPY least squares plus 1-norm solution and objective value:")
+    print(x1.value)
+    print(objective1.value)
+    
+# Now try another algorithm FBPD for same problem:
+x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init,Identity(), None, f, g1)
+print(x_fbpd1)
+print(criterfbpd1[-1])
+
+# Plot criterion curve to see both FISTA and FBPD converge to same value.
+# Note that FISTA is very efficient for 1-norm minimization so it beats
+# FBPD in this test by a lot. But FBPD can handle a larger class of problems 
+# than FISTA can.
+plt.figure()
+plt.loglog(iternum[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1')
+plt.loglog(iternum,criter1,label='FISTA LS+1')
+plt.legend()
+plt.show()
+
+plt.figure()
+plt.loglog(iternum[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1')
+plt.loglog(iternum,criter1,label='FISTA LS+1')
+plt.loglog(iternum,criterfbpd1,label='FBPD LS+1')
+plt.legend()
+plt.show()
+
+# Now try 1-norm and TV denoising with FBPD, first 1-norm.
+
+# 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 = 64
+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()
+
+# Identity operator for denoising
+I = TomoIdentity(ig)
+
+# Data and add noise
+y = I.direct(Phantom)
+y.array = y.array + 0.1*np.random.randn(N, N)
+
+plt.imshow(y.array)
+plt.title('Noisy image')
+plt.show()
+
+
+###################
+# Data fidelity term
+f_denoise = Norm2sq(I,y,c=0.5,memopt=True)
+
+# 1-norm regulariser
+lam1_denoise = 1.0
+g1_denoise = Norm1(lam1_denoise)
+
+# Initial guess
+x_init_denoise = ImageData(np.zeros((N,N)))
+
+# Combine with least squares and solve using generic FISTA implementation
+x_fista1_denoise, it1_denoise, timing1_denoise, criter1_denoise = FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt)
+
+print(x_fista1_denoise)
+print(criter1_denoise[-1])
+
+#plt.imshow(x_fista1_denoise.as_array())
+#plt.title('FISTA LS+1')
+#plt.show()
+
+# Now denoise LS + 1-norm with FBPD
+x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise, \
+  criterfbpd1_denoise = FBPD(x_init_denoise, I, None, f_denoise, g1_denoise)
+print(x_fbpd1_denoise)
+print(criterfbpd1_denoise[-1])
+
+#plt.imshow(x_fbpd1_denoise.as_array())
+#plt.title('FBPD LS+1')
+#plt.show()
+
+if use_cvxpy:
+    # Compare to CVXPY
+    
+    # Construct the problem.
+    x1_denoise = Variable(N**2,1)
+    objective1_denoise = Minimize(0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise,1) )
+    prob1_denoise = Problem(objective1_denoise)
+    
+    # The optimal objective is returned by prob.solve().
+    result1_denoise = prob1_denoise.solve(verbose=False,solver=SCS,eps=1e-12)
+    
+    # The optimal solution for x is stored in x.value and optimal objective value 
+    # is in result as well as in objective.value
+    print("CVXPY least squares plus 1-norm solution and objective value:")
+    print(x1_denoise.value)
+    print(objective1_denoise.value)
+
+x1_cvx = x1_denoise.value
+x1_cvx.shape = (N,N)
+
+
+
+#plt.imshow(x1_cvx)
+#plt.title('CVX LS+1')
+#plt.show()
+
+fig = plt.figure()
+plt.subplot(1,4,1)
+plt.imshow(y.array)
+plt.title("LS+1")
+plt.subplot(1,4,2)
+plt.imshow(x_fista1_denoise.as_array())
+plt.title("fista")
+plt.subplot(1,4,3)
+plt.imshow(x_fbpd1_denoise.as_array())
+plt.title("fbpd")
+plt.subplot(1,4,4)
+plt.imshow(x1_cvx)
+plt.title("cvx")
+plt.show()
+
+##############################################################
+# Now TV with FBPD and Norm2
+lam_tv = 0.1
+gtv = TV2D(lam_tv)
+norm2 = Norm2(lam_tv)
+op = FiniteDiff2D()
+#gtv(gtv.op.direct(x_init_denoise))
+
+opt_tv = {'tol': 1e-4, 'iter': 10000}
+
+x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, \
+ criterfbpdtv_denoise = FBPD(x_init_denoise, op, None, \
+                             f_denoise, norm2 ,opt=opt_tv)
+print(x_fbpdtv_denoise)
+print(criterfbpdtv_denoise[-1])
+
+plt.imshow(x_fbpdtv_denoise.as_array())
+plt.title('FBPD TV')
+#plt.show()
+
+if use_cvxpy:
+    # Compare to CVXPY
+    
+    # Construct the problem.
+    xtv_denoise = Variable((N,N))
+    #print (xtv_denoise.value.shape)
+    objectivetv_denoise = Minimize(0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise) )
+    probtv_denoise = Problem(objectivetv_denoise)
+    
+    # The optimal objective is returned by prob.solve().
+    resulttv_denoise = probtv_denoise.solve(verbose=False,solver=SCS,eps=1e-12)
+    
+    # The optimal solution for x is stored in x.value and optimal objective value 
+    # is in result as well as in objective.value
+    print("CVXPY least squares plus 1-norm solution and objective value:")
+    print(xtv_denoise.value)
+    print(objectivetv_denoise.value)
+    
+plt.imshow(xtv_denoise.value)
+plt.title('CVX TV')
+#plt.show()
+
+fig = plt.figure()
+plt.subplot(1,3,1)
+plt.imshow(y.array)
+plt.title("TV2D")
+plt.subplot(1,3,2)
+plt.imshow(x_fbpdtv_denoise.as_array())
+plt.title("fbpd tv denoise")
+plt.subplot(1,3,3)
+plt.imshow(xtv_denoise.value)
+plt.title("CVX tv")
+plt.show()
+
+
+
+plt.loglog([0,opt_tv['iter']], [objectivetv_denoise.value,objectivetv_denoise.value], label='CVX TV')
+plt.loglog(criterfbpdtv_denoise, label='FBPD TV')
diff --git a/Wrappers/Python/wip/old_demos/demo_gradient_descent.py b/Wrappers/Python/wip/old_demos/demo_gradient_descent.py
new file mode 100755
index 0000000..4d6647e
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/demo_gradient_descent.py
@@ -0,0 +1,295 @@
+
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer
+from ccpi.optimisation.algs import FISTA, FBPD, CGLS
+from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D, Norm2
+
+from ccpi.optimisation.ops import LinearOperatorMatrix, TomoIdentity
+from ccpi.optimisation.ops import Identity
+from ccpi.optimisation.ops import FiniteDiff2D
+
+# Requires CVXPY, see http://www.cvxpy.org/
+# CVXPY can be installed in anaconda using
+# conda install -c cvxgrp cvxpy libgcc
+
+# Whether to use or omit CVXPY
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+class Algorithm(object):
+    def __init__(self, *args, **kwargs):
+        pass
+    def set_up(self, *args, **kwargs):
+        raise NotImplementedError()
+    def update(self):
+        raise NotImplementedError()
+    
+    def should_stop(self):
+        raise NotImplementedError()
+    
+    def __iter__(self):
+        return self
+    
+    def __next__(self):
+        if self.should_stop():
+            raise StopIteration()
+        else:
+            self.update()
+        
+class GradientDescent(Algorithm):
+    x = None
+    rate = 0
+    objective_function = None
+    regulariser = None
+    iteration = 0
+    stop_cryterion = 'max_iter'
+    __max_iteration = 0
+    __loss = []
+    def __init__(self, **kwargs):
+        args = ['x_init', 'objective_function', 'rate']
+        present = True
+        for k,v in kwargs.items():
+            if k in args:
+                args.pop(args.index(k))
+        if len(args) == 0:
+            return self.set_up(x_init=kwargs['x_init'],
+                               objective_function=kwargs['objective_function'],
+                               rate=kwargs['rate'])
+    
+    def should_stop(self):
+        return self.iteration >= self.max_iteration
+    
+    def set_up(self, x_init, objective_function, rate):
+        self.x = x_init.copy()
+        self.x_update = x_init.copy()
+        self.objective_function = objective_function
+        self.rate = rate
+        self.__loss.append(objective_function(x_init))
+        
+    def update(self):
+        
+        self.objective_function.gradient(self.x, out=self.x_update)
+        self.x_update *= -self.rate
+        self.x += self.x_update
+        self.__loss.append(self.objective_function(self.x))
+        self.iteration += 1
+        
+    def get_output(self):
+        return self.x
+    def get_current_loss(self):
+        return self.__loss[-1]
+    @property
+    def loss(self):
+        return self.__loss
+    @property
+    def max_iteration(self):
+        return self.__max_iteration
+    @max_iteration.setter
+    def max_iteration(self, value):
+        assert isinstance(value, int)
+        self.__max_iteration = value
+        
+    
+
+
+
+# Problem data.
+m = 30
+n = 20
+np.random.seed(1)
+Amat = np.random.randn(m, n)
+A = LinearOperatorMatrix(Amat)
+bmat = np.random.randn(m)
+bmat.shape = (bmat.shape[0],1)
+
+# A = Identity()
+# Change n to equal to m.
+
+b = DataContainer(bmat)
+
+# Regularization parameter
+lam = 10
+opt = {'memopt':True}
+# Create object instances with the test data A and b.
+f = Norm2sq(A,b,c=0.5, memopt=True)
+g0 = ZeroFun()
+
+# Initial guess
+x_init = DataContainer(np.zeros((n,1)))
+
+f.grad(x_init)
+
+# Run FISTA for least squares plus zero function.
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0 , opt=opt)
+
+# Print solution and final objective/criterion value for comparison
+print("FISTA least squares plus zero function solution and objective value:")
+print(x_fista0.array)
+print(criter0[-1])
+
+gd = GradientDescent(x_init=x_init, objective_function=f, rate=0.001)
+gd.max_iteration = 5000
+
+for i,el in enumerate(gd):
+    if i%100 == 0:
+        print ("\rIteration {} Loss: {}".format(gd.iteration, 
+               gd.get_current_loss()))
+
+
+#%%
+
+
+#
+#if use_cvxpy:
+#    # Compare to CVXPY
+#    
+#    # Construct the problem.
+#    x0 = Variable(n)
+#    objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]) )
+#    prob0 = Problem(objective0)
+#    
+#    # The optimal objective is returned by prob.solve().
+#    result0 = prob0.solve(verbose=False,solver=SCS,eps=1e-9)
+#    
+#    # The optimal solution for x is stored in x.value and optimal objective value 
+#    # is in result as well as in objective.value
+#    print("CVXPY least squares plus zero function solution and objective value:")
+#    print(x0.value)
+#    print(objective0.value)
+#
+## Plot criterion curve to see FISTA converge to same value as CVX.
+#iternum = np.arange(1,1001)
+#plt.figure()
+#plt.loglog(iternum[[0,-1]],[objective0.value, objective0.value], label='CVX LS')
+#plt.loglog(iternum,criter0,label='FISTA LS')
+#plt.legend()
+#plt.show()
+#
+## Create 1-norm object instance
+#g1 = Norm1(lam)
+#
+#g1(x_init)
+#x_rand = DataContainer(np.reshape(np.random.rand(n),(n,1)))
+#x_rand2 = DataContainer(np.reshape(np.random.rand(n-1),(n-1,1)))
+#v = g1.prox(x_rand,0.02)
+##vv = g1.prox(x_rand2,0.02)
+#vv = v.copy() 
+#vv *= 0
+#print (">>>>>>>>>>vv" , vv.as_array())
+#vv.fill(v)
+#print (">>>>>>>>>>fill" , vv.as_array())
+#g1.proximal(x_rand, 0.02, out=vv)
+#print (">>>>>>>>>>v" , v.as_array())
+#print (">>>>>>>>>>gradient" , vv.as_array())
+#
+#print (">>>>>>>>>>" , (v-vv).as_array())
+#import sys
+##sys.exit(0)
+## Combine with least squares and solve using generic FISTA implementation
+#x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1,opt=opt)
+#
+## Print for comparison
+#print("FISTA least squares plus 1-norm solution and objective value:")
+#print(x_fista1)
+#print(criter1[-1])
+#
+#if use_cvxpy:
+#    # Compare to CVXPY
+#    
+#    # Construct the problem.
+#    x1 = Variable(n)
+#    objective1 = Minimize(0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1,1) )
+#    prob1 = Problem(objective1)
+#    
+#    # The optimal objective is returned by prob.solve().
+#    result1 = prob1.solve(verbose=False,solver=SCS,eps=1e-9)
+#    
+#    # The optimal solution for x is stored in x.value and optimal objective value 
+#    # is in result as well as in objective.value
+#    print("CVXPY least squares plus 1-norm solution and objective value:")
+#    print(x1.value)
+#    print(objective1.value)
+#    
+## Now try another algorithm FBPD for same problem:
+#x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init,Identity(), None, f, g1)
+#print(x_fbpd1)
+#print(criterfbpd1[-1])
+#
+## Plot criterion curve to see both FISTA and FBPD converge to same value.
+## Note that FISTA is very efficient for 1-norm minimization so it beats
+## FBPD in this test by a lot. But FBPD can handle a larger class of problems 
+## than FISTA can.
+#plt.figure()
+#plt.loglog(iternum[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1')
+#plt.loglog(iternum,criter1,label='FISTA LS+1')
+#plt.legend()
+#plt.show()
+#
+#plt.figure()
+#plt.loglog(iternum[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1')
+#plt.loglog(iternum,criter1,label='FISTA LS+1')
+#plt.loglog(iternum,criterfbpd1,label='FBPD LS+1')
+#plt.legend()
+#plt.show()
+
+# Now try 1-norm and TV denoising with FBPD, first 1-norm.
+
+# 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 = 64
+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()
+
+# Identity operator for denoising
+I = TomoIdentity(ig)
+
+# Data and add noise
+y = I.direct(Phantom)
+y.array = y.array + 0.1*np.random.randn(N, N)
+
+plt.imshow(y.array)
+plt.title('Noisy image')
+plt.show()
+
+
+###################
+# Data fidelity term
+f_denoise = Norm2sq(I,y,c=0.5,memopt=True)
+
+# 1-norm regulariser
+lam1_denoise = 1.0
+g1_denoise = Norm1(lam1_denoise)
+
+# Initial guess
+x_init_denoise = ImageData(np.zeros((N,N)))
+
+# Combine with least squares and solve using generic FISTA implementation
+x_fista1_denoise, it1_denoise, timing1_denoise, criter1_denoise = \
+   FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt)
+
+print(x_fista1_denoise)
+print(criter1_denoise[-1])
+
+f_2 = 
+gd = GradientDescent(x_init=x_init_denoise, 
+                     objective_function=f, rate=0.001)
+gd.max_iteration = 5000
+
+for i,el in enumerate(gd):
+    if i%100 == 0:
+        print ("\rIteration {} Loss: {}".format(gd.iteration, 
+               gd.get_current_loss()))
+
+plt.imshow(gd.get_output().as_array())
+plt.title('GD image')
+plt.show()
+    
diff --git a/Wrappers/Python/wip/old_demos/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/old_demos/demo_imat_multichan_RGLTK.py
new file mode 100644
index 0000000..8370c78
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/demo_imat_multichan_RGLTK.py
@@ -0,0 +1,151 @@
+# This script demonstrates how to load IMAT fits data
+# into the CIL optimisation framework and run reconstruction methods.
+#
+# Demo to reconstruct energy-discretized channels of IMAT data
+
+# needs dxchange: conda install -c conda-forge dxchange
+# needs astropy: conda install astropy
+
+
+# All third-party imports.
+import numpy as np
+import matplotlib.pyplot as plt
+from dxchange.reader import read_fits
+from astropy.io import fits
+
+# All own imports.
+from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData, DataContainer
+from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple
+from ccpi.optimisation.algs import CGLS, FISTA
+from ccpi.optimisation.funcs import Norm2sq, Norm1
+from ccpi.plugins.regularisers import FGP_TV
+
+# set main parameters here
+n = 512 
+totalAngles = 250 # total number of projection angles
+# spectral discretization parameter
+num_average = 145 # channel discretization frequency (total number of averaged channels)
+numChannels = 2970 # 2970
+totChannels = round(numChannels/num_average) # the resulting number of channels
+Projections_stack = np.zeros((num_average,n,n),dtype='uint16')
+ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32')
+
+#########################################################################
+print ("Loading the data...")
+MainPath = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/' # path to data
+pathname0 = '{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/Sample/')
+counterFileName = 4675
+# A main loop over all available angles 
+for ll in range(0,totalAngles,1):
+    pathnameData = '{!s}{!s}{!s}/'.format(pathname0,'angle',str(ll))
+    filenameCurr = '{!s}{!s}{!s}'.format('IMAT0000',str(counterFileName),'_Tomo_test_000_')
+    counterT = 0
+    # loop over reduced channels (discretized)
+    for i in range(0,totChannels,1):
+        sumCount = 0
+        # loop over actual channels to obtain averaged one
+        for j in range(0,num_average,1):
+            if counterT < 10:
+                outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'0000',str(counterT))
+            if ((counterT  >= 10) & (counterT < 100)):
+                outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'000',str(counterT))
+            if ((counterT  >= 100) & (counterT < 1000)):
+                outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'00',str(counterT))
+            if ((counterT  >= 1000) & (counterT < 10000)):
+                outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'0',str(counterT))
+            try:
+                Projections_stack[j,:,:] = read_fits(outfile)
+            except:
+                print("Fits is corrupted, skipping no.", counterT)
+                sumCount -= 1
+            counterT += 1
+            sumCount += 1
+        AverageProj=np.sum(Projections_stack,axis=0)/sumCount # averaged projection over "num_average" channels
+        ProjAngleChannels[ll,i,:,:] = AverageProj
+    print("Angle is processed", ll)
+    counterFileName += 1
+#########################################################################
+
+flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits'))
+nonzero = flat1 > 0
+# Apply flat field and take negative log
+for ll in range(0,totalAngles,1):
+    for i in range(0,totChannels,1):
+        ProjAngleChannels[ll,i,nonzero] = ProjAngleChannels[ll,i,nonzero]/flat1[nonzero]
+
+eqzero = ProjAngleChannels == 0
+ProjAngleChannels[eqzero] = 1
+ProjAngleChannels_NormLog = -np.log(ProjAngleChannels) # normalised and neg-log data
+
+# extact sinogram over energy channels
+selectedVertical_slice = 256
+sino_all_channels = ProjAngleChannels_NormLog[:,:,:,selectedVertical_slice]
+# Set angles to use
+angles = np.linspace(-np.pi,np.pi,totalAngles,endpoint=False)
+
+# set the geometry
+ig = ImageGeometry(n,n)
+ag = AcquisitionGeometry('parallel',
+                             '2D',
+                             angles,
+                             n,1)
+Aop = AstraProjectorSimple(ig, ag, 'gpu')
+
+
+# loop to reconstruct energy channels
+REC_chan = np.zeros((totChannels, n, n), 'float32')
+for i in range(0,totChannels,1):
+    sino_channel = sino_all_channels[:,i,:] # extract a sinogram for i-th channel
+
+    print ("Initial guess")
+    x_init = ImageData(geometry=ig)
+    
+    # Create least squares object instance with projector and data.
+    print ("Create least squares object instance with projector and data.")
+    f = Norm2sq(Aop,DataContainer(sino_channel),c=0.5)
+    
+    print ("Run FISTA-TV for least squares")
+    lamtv = 5
+    opt = {'tol': 1e-4, 'iter': 200}
+    g_fgp = FGP_TV(lambdaReg = lamtv,
+                 iterationsTV=50,
+                 tolerance=1e-6,
+                 methodTV=0,
+                 nonnegativity=0,
+                 printing=0,
+                 device='gpu')
+    
+    x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp, opt)
+    REC_chan[i,:,:] = x_fista_fgp.array
+    """
+    plt.figure()
+    plt.subplot(121)
+    plt.imshow(x_fista_fgp.array, vmin=0, vmax=0.05)
+    plt.title('FISTA FGP TV')
+    plt.subplot(122)
+    plt.semilogy(criter_fgp)
+    plt.show()
+    """
+    """
+    print ("Run CGLS for least squares")
+    opt = {'tol': 1e-4, 'iter': 20}
+    x_init = ImageData(geometry=ig)
+    x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, DataContainer(sino_channel), opt=opt)
+    
+    plt.figure()
+    plt.imshow(x_CGLS.array,vmin=0, vmax=0.05)
+    plt.title('CGLS')
+    plt.show()
+    """
+# Saving images into fits using astrapy if required
+counter = 0
+filename = 'FISTA_TV_imat_slice'
+for i in range(totChannels):
+    im = REC_chan[i,:,:] 
+    add_val = np.min(im[:])
+    im += abs(add_val)
+    im = im/np.max(im[:])*65535
+    outfile = '{!s}_{!s}_{!s}.fits'.format(filename,str(selectedVertical_slice),str(counter))
+    hdu = fits.PrimaryHDU(np.uint16(im))
+    hdu.writeto(outfile, overwrite=True)
+    counter += 1
\ No newline at end of file
diff --git a/Wrappers/Python/wip/old_demos/demo_imat_whitebeam.py b/Wrappers/Python/wip/old_demos/demo_imat_whitebeam.py
new file mode 100644
index 0000000..e0d213e
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/demo_imat_whitebeam.py
@@ -0,0 +1,138 @@
+# This script demonstrates how to load IMAT fits data
+# into the CIL optimisation framework and run reconstruction methods.
+#
+# This demo loads the summedImg files which are the non-spectral images 
+# resulting from summing projections over all spectral channels.
+
+# needs dxchange: conda install -c conda-forge dxchange
+# needs astropy: conda install astropy
+
+
+# All third-party imports.
+import numpy
+from scipy.io import loadmat
+import matplotlib.pyplot as plt
+from dxchange.reader import read_fits
+
+# All own imports.
+from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData
+from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple
+from ccpi.optimisation.algs import CGLS, FISTA
+from ccpi.optimisation.funcs import Norm2sq, Norm1
+
+# Load and display a couple of summed projection as examples
+pathname0 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/'
+filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits'
+
+data0 = read_fits(pathname0 + filename0)
+
+pathname10 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/'
+filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits'
+
+data10 = read_fits(pathname10 + filename10)
+
+# Load a flat field (more are available, should we average over them?)
+flat1 = read_fits('/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')
+
+# Apply flat field and display after flat-field correction and negative log
+data0_rel = numpy.zeros(numpy.shape(flat1), dtype = float)
+nonzero = flat1 > 0
+data0_rel[nonzero] = data0[nonzero] / flat1[nonzero]
+data10_rel =  numpy.zeros(numpy.shape(flat1), dtype = float)
+data10_rel[nonzero] = data10[nonzero] / flat1[nonzero]
+
+plt.imshow(data0_rel)
+plt.colorbar()
+plt.show()
+
+plt.imshow(-numpy.log(data0_rel))
+plt.colorbar()
+plt.show()
+
+plt.imshow(data10_rel)
+plt.colorbar()
+plt.show()
+
+plt.imshow(-numpy.log(data10_rel))
+plt.colorbar()
+plt.show()
+
+# Set up for loading all summed images at 250 angles.
+pathname = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/'
+filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits'
+
+# Dimensions
+num_angles = 250
+imsize = 512
+
+# Initialise array
+data = numpy.zeros((num_angles,imsize,imsize))
+
+# Load only 0-249, as 250 is at repetition of zero degrees just like image 0
+for i in range(0,250):
+    curimfile = (pathname + filename).format(i, i+4675)
+    data[i,:,:] = read_fits(curimfile)
+
+# Apply flat field and take negative log
+nonzero = flat1 > 0
+for i in range(0,250):
+    data[i,nonzero] = data[i,nonzero]/flat1[nonzero]
+
+eqzero = data == 0
+data[eqzero] = 1
+
+data_rel = -numpy.log(data)
+
+# Permute order to get: angles, vertical, horizontal, as default in framework.
+data_rel = numpy.transpose(data_rel,(0,2,1))
+
+# Set angles to use
+angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False)
+
+# Create 3D acquisition geometry and acquisition data
+ag = AcquisitionGeometry('parallel',
+                         '3D',
+                         angles,
+                         pixel_num_h=imsize,
+                         pixel_num_v=imsize)
+b = AcquisitionData(data_rel, geometry=ag)
+
+# Reduce to single (noncentral) slice by extracting relevant parameters from data and its
+# geometry. Perhaps create function to extract central slice automatically?
+b2d = b.subset(vertical=128)
+ag2d = AcquisitionGeometry('parallel',
+                         '2D',
+                         ag.angles,
+                         pixel_num_h=ag.pixel_num_h)
+b2d.geometry = ag2d
+
+# Create 2D image geometry
+ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, 
+                     voxel_num_y=ag2d.pixel_num_h)
+
+# Create GPU projector/backprojector operator with ASTRA.
+Aop = AstraProjectorSimple(ig2d,ag2d,'gpu')
+
+# Demonstrate operator is working by applying simple backprojection.
+z = Aop.adjoint(b2d)
+plt.imshow(z.array)
+plt.title('Simple backprojection')
+plt.colorbar()
+plt.show()
+
+# Set initial guess ImageData with zeros for algorithms, and algorithm options.
+x_init = ImageData(numpy.zeros((imsize,imsize)),
+                   geometry=ig2d)
+opt_CGLS = {'tol': 1e-4, 'iter': 20}
+
+# Run CGLS algorithm and display reconstruction.
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b2d, opt_CGLS)
+
+plt.imshow(x_CGLS.array)
+plt.title('CGLS')
+plt.colorbar()
+plt.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS Criterion vs iterations')
+plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/old_demos/demo_memhandle.py b/Wrappers/Python/wip/old_demos/demo_memhandle.py
new file mode 100755
index 0000000..db48d73
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/demo_memhandle.py
@@ -0,0 +1,193 @@
+
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer
+from ccpi.optimisation.algs import FISTA, FBPD, CGLS
+from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D
+
+from ccpi.optimisation.ops import LinearOperatorMatrix, Identity
+from ccpi.optimisation.ops import TomoIdentity
+
+# Requires CVXPY, see http://www.cvxpy.org/
+# CVXPY can be installed in anaconda using
+# conda install -c cvxgrp cvxpy libgcc
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Problem data.
+m = 30
+n = 20
+np.random.seed(1)
+Amat = np.random.randn(m, n)
+A = LinearOperatorMatrix(Amat)
+bmat = np.random.randn(m)
+bmat.shape = (bmat.shape[0],1)
+
+
+
+# A = Identity()
+# Change n to equal to m.
+
+b = DataContainer(bmat)
+
+# Regularization parameter
+lam = 10
+
+# Create object instances with the test data A and b.
+f = Norm2sq(A,b,c=0.5, memopt=True)
+g0 = ZeroFun()
+
+# Initial guess
+x_init = DataContainer(np.zeros((n,1)))
+
+f.grad(x_init)
+opt = {'memopt': True}
+# Run FISTA for least squares plus zero function.
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0)
+x_fista0_m, it0_m, timing0_m, criter0_m = FISTA(x_init, f, g0, opt=opt)
+
+iternum = [i for i in range(len(criter0))]
+# Print solution and final objective/criterion value for comparison
+print("FISTA least squares plus zero function solution and objective value:")
+print(x_fista0.array)
+print(criter0[-1])
+
+
+# Plot criterion curve to see FISTA converge to same value as CVX.
+#iternum = np.arange(1,1001)
+plt.figure()
+plt.loglog(iternum,criter0,label='FISTA LS')
+plt.loglog(iternum,criter0_m,label='FISTA LS memopt')
+plt.legend()
+plt.show()
+#%%
+# Create 1-norm object instance
+g1 = Norm1(lam)
+
+g1(x_init)
+g1.prox(x_init,0.02)
+
+# Combine with least squares and solve using generic FISTA implementation
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1)
+x_fista1_m, it1_m, timing1_m, criter1_m = FISTA(x_init, f, g1, opt=opt)
+iternum = [i for i in range(len(criter1))]
+# Print for comparison
+print("FISTA least squares plus 1-norm solution and objective value:")
+print(x_fista1)
+print(criter1[-1])
+
+
+# Now try another algorithm FBPD for same problem:
+x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, None, f, g1)
+iternum = [i for i in range(len(criterfbpd1))]
+print(x_fbpd1)
+print(criterfbpd1[-1])
+
+# Plot criterion curve to see both FISTA and FBPD converge to same value.
+# Note that FISTA is very efficient for 1-norm minimization so it beats
+# FBPD in this test by a lot. But FBPD can handle a larger class of problems 
+# than FISTA can.
+plt.figure()
+plt.loglog(iternum,criter1,label='FISTA LS+1')
+plt.loglog(iternum,criter1_m,label='FISTA LS+1 memopt')
+plt.legend()
+plt.show()
+
+plt.figure()
+plt.loglog(iternum,criter1,label='FISTA LS+1')
+plt.loglog(iternum,criterfbpd1,label='FBPD LS+1')
+plt.legend()
+plt.show()
+#%%
+# Now try 1-norm and TV denoising with FBPD, first 1-norm.
+
+# 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 = 1000
+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()
+
+# Identity operator for denoising
+I = TomoIdentity(ig)
+
+# Data and add noise
+y = I.direct(Phantom)
+y.array +=  0.1*np.random.randn(N, N)
+
+plt.figure()
+plt.imshow(y.array)
+plt.title('Noisy image')
+plt.show()
+
+# Data fidelity term
+f_denoise = Norm2sq(I,y,c=0.5, memopt=True)
+
+# 1-norm regulariser
+lam1_denoise = 1.0
+g1_denoise = Norm1(lam1_denoise)
+
+# Initial guess
+x_init_denoise = ImageData(np.zeros((N,N)))
+opt = {'memopt': False, 'iter' : 50}
+# Combine with least squares and solve using generic FISTA implementation
+print ("no memopt")
+x_fista1_denoise, it1_denoise, timing1_denoise, \
+        criter1_denoise = FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt)
+opt = {'memopt': True, 'iter' : 50}        
+print ("yes memopt")
+x_fista1_denoise_m, it1_denoise_m, timing1_denoise_m, \
+      criter1_denoise_m = FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt)
+
+print(x_fista1_denoise)
+print(criter1_denoise[-1])
+
+plt.figure()
+plt.imshow(x_fista1_denoise.as_array())
+plt.title('FISTA LS+1')
+plt.show()
+
+plt.figure()
+plt.imshow(x_fista1_denoise_m.as_array())
+plt.title('FISTA LS+1 memopt')
+plt.show()
+
+plt.figure()
+plt.loglog(iternum,criter1_denoise,label='FISTA LS+1')
+plt.loglog(iternum,criter1_denoise_m,label='FISTA LS+1 memopt')
+plt.legend()
+plt.show()
+#%%
+# Now denoise LS + 1-norm with FBPD
+x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise, criterfbpd1_denoise = FBPD(x_init_denoise, None, f_denoise, g1_denoise)
+print(x_fbpd1_denoise)
+print(criterfbpd1_denoise[-1])
+
+plt.figure()
+plt.imshow(x_fbpd1_denoise.as_array())
+plt.title('FBPD LS+1')
+plt.show()
+
+
+# Now TV with FBPD
+lam_tv = 0.1
+gtv = TV2D(lam_tv)
+gtv(gtv.op.direct(x_init_denoise))
+
+opt_tv = {'tol': 1e-4, 'iter': 10000}
+
+x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, criterfbpdtv_denoise = FBPD(x_init_denoise, None, f_denoise, gtv,opt=opt_tv)
+print(x_fbpdtv_denoise)
+print(criterfbpdtv_denoise[-1])
+
+plt.imshow(x_fbpdtv_denoise.as_array())
+plt.title('FBPD TV')
+plt.show()
diff --git a/Wrappers/Python/wip/old_demos/demo_test_sirt.py b/Wrappers/Python/wip/old_demos/demo_test_sirt.py
new file mode 100644
index 0000000..6f5a44d
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/demo_test_sirt.py
@@ -0,0 +1,176 @@
+# 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.optimisation.funcs import Norm2sq, Norm1, TV2D, IndicatorBox
+from ccpi.astra.ops import AstraProjectorSimple
+
+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.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, 'gpu')
+
+# 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.imshow(b.array)
+plt.title('Simulated data')
+plt.show()
+
+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': 1000}
+
+# 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.colorbar()
+plt.show()
+
+plt.semilogy(criter_CGLS)
+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)
+plt.title('SIRT unconstrained')
+plt.colorbar()
+plt.show()
+
+plt.semilogy(criter_SIRT)
+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))
+
+plt.imshow(x_SIRT0.array)
+plt.title('SIRT nonneg')
+plt.colorbar()
+plt.show()
+
+plt.semilogy(criter_SIRT0)
+plt.title('SIRT nonneg criterion')
+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)')
+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.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)
+
+plt.imshow(x_fista0.array)
+plt.title('FISTA Least squares nonneg')
+plt.show()
+
+plt.semilogy(criter0)
+plt.title('FISTA Least squares nonneg 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)
+
+plt.imshow(x_fista01.array)
+plt.title('FISTA Least squares box(0,1)')
+plt.show()
+
+plt.semilogy(criter01)
+plt.title('FISTA Least squares box(0,1) criterion')
+plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/old_demos/multifile_nexus.py b/Wrappers/Python/wip/old_demos/multifile_nexus.py
new file mode 100755
index 0000000..d1ad438
--- /dev/null
+++ b/Wrappers/Python/wip/old_demos/multifile_nexus.py
@@ -0,0 +1,307 @@
+# -*- coding: utf-8 -*-

+"""

+Created on Wed Aug 15 16:00:53 2018

+

+@author: ofn77899

+"""

+

+import os

+from ccpi.io.reader import NexusReader

+

+from sys import getsizeof

+

+import matplotlib.pyplot as plt

+

+from ccpi.framework import DataProcessor, DataContainer

+from ccpi.processors import Normalizer

+from ccpi.processors import CenterOfRotationFinder

+import numpy

+

+

+class averager(object):

+    def __init__(self):

+        self.reset()

+        

+    def reset(self):

+        self.N = 0

+        self.avg = 0

+        self.min = 0

+        self.max = 0

+        self.var = 0

+        self.ske = 0

+        self.kur = 0

+    

+    def add_reading(self, val):

+        

+        if (self.N == 0):

+           self.avg = val;

+           self.min = val;

+           self.max = val;

+        elif (self.N == 1):

+        #//set min/max

+            self.max = val if val > self.max else self.max

+            self.min = val if val < self.min else self.min

+            

+      

+            thisavg = (self.avg + val)/2

+            #// initial value is in avg

+            self.var = (self.avg - thisavg)*(self.avg-thisavg) + (val - thisavg) * (val-thisavg)

+            self.ske = self.var * (self.avg - thisavg)

+            self.kur = self.var * self.var

+            self.avg = thisavg

+        else:

+            self.max = val if val > self.max else self.max

+            self.min = val if val < self.min else self.min

+            

+            M = self.N

+

+            #// b-factor =(<v>_N + v_(N+1)) / (N+1)

+            #float b = (val -avg) / (M+1);

+            b = (val -self.avg) / (M+1)

+            

+            self.kur = self.kur + (M *(b*b*b*b) * (1+M*M*M))- (4* b * self.ske) + (6 * b *b * self.var * (M-1))

+    

+            self.ske = self.ske + (M * b*b*b *(M-1)*(M+1)) - (3 * b * self.var * (M-1))

+    

+            #//var = var * ((M-1)/M) + ((val - avg)*(val - avg)/(M+1)) ;

+            self.var = self.var * ((M-1)/M) + (b * b * (M+1))

+            

+            self.avg = self.avg * (M/(M+1)) + val / (M+1)

+    

+        self.N += 1

+    

+    def stats(self, vector):

+        i = 0

+        while i < vector.size:

+            self.add_reading(vector[i])

+            i+=1

+

+avg = averager()

+a = numpy.linspace(0,39,40)

+avg.stats(a)

+print ("average" , avg.avg, a.mean())

+print ("variance" , avg.var, a.var())

+b = a - a.mean()

+b *= b

+b = numpy.sqrt(sum(b)/(a.size-1))

+print ("std" , numpy.sqrt(avg.var), b)

+#%%         

+

+class DataStatMoments(DataProcessor):

+    '''Normalization based on flat and dark

+    

+    This processor read in a AcquisitionData and normalises it based on 

+    the instrument reading with and without incident photons or neutrons.

+    

+    Input: AcquisitionData

+    Parameter: 2D projection with flat field (or stack)

+               2D projection with dark field (or stack)

+    Output: AcquisitionDataSetn

+    '''

+    

+    def __init__(self, axis, skewness=False, kurtosis=False, offset=0):

+        kwargs = {

+                  'axis'     : axis,

+                  'skewness' : skewness,

+                  'kurtosis' : kurtosis,

+                  'offset'   : offset,

+                  }

+        #DataProcessor.__init__(self, **kwargs)

+        super(DataStatMoments, self).__init__(**kwargs)

+        

+    

+    def check_input(self, dataset):

+        #self.N = dataset.get_dimension_size(self.axis)

+        return True

+    

+    @staticmethod

+    def add_sample(dataset, N, axis, stats=None, offset=0):

+        #dataset = self.get_input()

+        if (N == 0):

+            # get the axis number along to calculate the stats

+            

+            

+            #axs = dataset.get_dimension_size(self.axis)

+            # create a placeholder for the output

+            if stats is None:

+                ax = dataset.get_dimension_axis(axis)

+                shape = [dataset.shape[i] for i in range(len(dataset.shape)) if i != ax]

+                # output has 4 components avg, std, skewness and kurtosis + last avg+ (val-thisavg)

+                shape.insert(0, 4+2)

+                stats = numpy.zeros(shape)

+                

+            stats[0] = dataset.subset(**{axis:N-offset}).array[:]

+            

+            #avg = val

+        elif (N == 1):

+            # val

+            stats[5] = dataset.subset(**{axis:N-offset}).array

+            stats[4] = stats[0] + stats[5]

+            stats[4] /= 2         # thisavg

+            stats[5] -= stats[4]  # (val - thisavg) 

+            

+            #float thisavg = (avg + val)/2;

+            

+            #// initial value is in avg

+            #var = (avg - thisavg)*(avg-thisavg) + (val - thisavg) * (val-thisavg);

+            stats[1] = stats[5] * stats[5] + stats[5] * stats[5]

+            #skewness = var * (avg - thisavg);

+            stats[2] = stats[1] * stats[5]

+            #kurtosis = var * var;

+            stats[3] = stats[1] * stats[1]

+            #avg = thisavg;

+            stats[0] = stats[4]

+        else:

+        

+            #float M = (float)N;

+            M = N

+            #val

+            stats[4] = dataset.subset(**{axis:N-offset}).array

+            #// b-factor =(<v>_N + v_(N+1)) / (N+1)

+            stats[5] = stats[4] - stats[0]

+            stats[5] /= (M+1)    # b factor

+            #float b = (val -avg) / (M+1);

+    

+            #kurtosis = kurtosis + (M *(b*b*b*b) * (1+M*M*M))- (4* b * skewness) + (6 * b *b * var * (M-1));

+            #if self.kurtosis:

+            #    stats[3] += (M * stats[5] * stats[5] * stats[5] * stats[5]) - \

+            #                (4 * stats[5] * stats[2]) + \

+            #                (6 * stats[5] * stats[5] * stats[1] * (M-1))

+                        

+            #skewness = skewness + (M * b*b*b *(M-1)*(M+1)) - (3 * b * var * (M-1));

+            #if self.skewness:

+            #    stats[2] = stats[2] +  (M * stats[5]* stats[5] * stats[5] * (M-1)*(M-1) ) -\

+            #                3 * stats[5] * stats[1] * (M-1) 

+            #//var = var * ((M-1)/M) + ((val - avg)*(val - avg)/(M+1)) ;

+            #var = var * ((M-1)/M) + (b * b * (M+1));

+            stats[1] = ((M-1)/M) * stats[1] + (stats[5] * stats[5] * (M+1))

+           

+            #avg = avg * (M/(M+1)) + val / (M+1)

+            stats[0] = stats[0] * (M/(M+1)) + stats[4] / (M+1)

+    

+        N += 1

+        return stats , N

+    

+           

+    def process(self):

+        

+        data = self.get_input()

+        

+        #stats, i = add_sample(0)

+        N = data.get_dimension_size(self.axis)

+        ax = data.get_dimension_axis(self.axis)

+        stats = None

+        i = 0

+        while i < N:

+            stats , i = DataStatMoments.add_sample(data, i, self.axis, stats, offset=self.offset)

+        self.offset += N

+        labels = ['StatisticalMoments']

+        

+        labels += [data.dimension_labels[i] \

+                   for i in range(len(data.dimension_labels)) if i != ax]

+        y = DataContainer( stats[:4] , False, 

+                    dimension_labels=labels)

+        return y

+

+directory = r'E:\Documents\Dataset\CCPi\Nexus_test'

+data_path="entry1/instrument/pco1_hw_hdf_nochunking/data"

+

+reader = NexusReader(os.path.join( os.path.abspath(directory) , '74331.nxs'))

+

+print ("read flat")

+read_flat = NexusReader(os.path.join( os.path.abspath(directory) , '74240.nxs'))

+read_flat.data_path = data_path

+flatsslice = read_flat.get_acquisition_data_whole()

+avg = DataStatMoments('angle')

+

+avg.set_input(flatsslice)

+flats = avg.get_output()

+

+ave = averager()

+ave.stats(flatsslice.array[:,0,0])

+

+print ("avg" , ave.avg, flats.array[0][0][0])

+print ("var" , ave.var, flats.array[1][0][0])

+

+print ("read dark")

+read_dark = NexusReader(os.path.join( os.path.abspath(directory) , '74243.nxs'))

+read_dark.data_path = data_path

+

+## darks are very many so we proceed in batches

+total_size = read_dark.get_projection_dimensions()[0]

+print ("total_size", total_size)

+

+batchsize = 40

+if batchsize > total_size:

+    batchlimits = [batchsize * (i+1) for i in range(int(total_size/batchsize))] + [total_size-1]

+else:

+    batchlimits = [total_size]

+#avg.N = 0

+avg.offset = 0

+N = 0

+for batch in range(len(batchlimits)):  

+    print ("running batch " , batch)

+    bmax = batchlimits[batch]

+    bmin = bmax-batchsize

+    

+    darksslice = read_dark.get_acquisition_data_batch(bmin,bmax)

+    if batch == 0:

+        #create stats

+        ax = darksslice.get_dimension_axis('angle')

+        shape = [darksslice.shape[i] for i in range(len(darksslice.shape)) if i != ax]

+        # output has 4 components avg, std, skewness and kurtosis + last avg+ (val-thisavg)

+        shape.insert(0, 4+2)

+        print ("create stats shape ", shape)

+        stats = numpy.zeros(shape)

+    print ("N" , N)        

+    #avg.set_input(darksslice)

+    i = bmin

+    while i < bmax:

+        stats , i = DataStatMoments.add_sample(darksslice, i, 'angle', stats, bmin)

+        print ("{0}-{1}-{2}".format(bmin, i , bmax ) )

+    

+darks = stats

+#%%

+

+fig = plt.subplot(2,2,1)

+fig.imshow(flats.subset(StatisticalMoments=0).array)

+fig = plt.subplot(2,2,2)

+fig.imshow(numpy.sqrt(flats.subset(StatisticalMoments=1).array))

+fig = plt.subplot(2,2,3)

+fig.imshow(darks[0])

+fig = plt.subplot(2,2,4)

+fig.imshow(numpy.sqrt(darks[1]))

+plt.show()

+

+

+#%%

+norm = Normalizer(flat_field=flats.array[0,200,:], dark_field=darks[0,200,:])

+#norm.set_flat_field(flats.array[0,200,:])

+#norm.set_dark_field(darks.array[0,200,:])

+norm.set_input(reader.get_acquisition_data_slice(200))

+

+n = Normalizer.normalize_projection(norm.get_input().as_array(), flats.array[0,200,:], darks[0,200,:], 1e-5)

+#dn_n= Normalizer.estimate_normalised_error(norm.get_input().as_array(), flats.array[0,200,:], darks[0,200,:],

+#                                           numpy.sqrt(flats.array[1,200,:]), numpy.sqrt(darks[1,200,:]))

+#%%

+

+

+#%%

+fig = plt.subplot(2,1,1)

+

+

+fig.imshow(norm.get_input().as_array())

+fig = plt.subplot(2,1,2)

+fig.imshow(n)

+

+#fig = plt.subplot(3,1,3)

+#fig.imshow(dn_n)

+

+

+plt.show()

+

+

+

+

+

+