Merge pull request #352 from vais-ral/Remove_demos

Removed demos

23 files changed, 0 insertions, 4526 deletions

diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
deleted file mode 100644
index 2ac050f..0000000
--- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
+++ /dev/null
@@ -1,211 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-Conjugate Gradient for (Regularized) Least Squares for Tomography
-
-
-Problem:     min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2}
-
-             A: Projection operator
-             g: Sinogram
-             L: Identity or Gradient Operator
-
-"""
-
-
-from ccpi.framework import ImageGeometry, ImageData, \
-                            AcquisitionGeometry, BlockDataContainer, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import CGLS
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
-       
-import tomophantom
-from tomophantom import TomoP2D
-from ccpi.astra.operators import AstraProjectorSimple 
-import os
-
-
-# Load  Shepp-Logan phantom 
-model = 1 # select a model number from the library
-N = 64 # set dimension of the phantom
-path = os.path.dirname(tomophantom.__file__)
-path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
-phantom_2D = TomoP2D.Model(model, N, path_library2D)
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-data = ImageData(phantom_2D)
-
-detectors =  N
-angles = np.linspace(0, np.pi, 180, dtype=np.float32)
-
-ag = AcquisitionGeometry('parallel','2D', angles, detectors)
-
-device = input('Available device: GPU==1 / CPU==0 ')
-
-if device =='1':
-    dev = 'gpu'
-else:
-    dev = 'cpu'
-
-Aop = AstraProjectorSimple(ig, ag, dev)    
-sin = Aop.direct(data)
-
-np.random.seed(10)
-noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,1,ag.shape))
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-plt.subplot(2,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(2,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Setup and run the simple CGLS algorithm  
-x_init = ig.allocate()  
-
-cgls1 = CGLS(x_init = x_init, operator = Aop, data = noisy_data)
-cgls1.max_iteration = 20
-cgls1.update_objective_interval = 5
-cgls1.run(20, verbose = True)
-
-# Setup and run the regularised CGLS algorithm  (Tikhonov with Identity)
-
-x_init = ig.allocate()  
-alpha1 = 50
-op1 = Identity(ig)
-
-block_op1 = BlockOperator( Aop, alpha1 * op1, shape=(2,1))
-block_data1 = BlockDataContainer(noisy_data, op1.range_geometry().allocate())
-   
-cgls2 = CGLS(x_init = x_init, operator = block_op1, data = block_data1)
-cgls2.max_iteration = 200
-cgls2.update_objective_interval = 10
-cgls2.run(200, verbose=True)
-
-# Setup and run the regularised CGLS algorithm  (Tikhonov with Gradient)
-
-x_init = ig.allocate() 
-alpha2 = 25
-op2 = Gradient(ig)
-
-block_op2 = BlockOperator( Aop, alpha2 * op2, shape=(2,1))
-block_data2 = BlockDataContainer(noisy_data, op2.range_geometry().allocate())
-   
-cgls3 = CGLS(x_init=x_init, operator = block_op2, data = block_data2)
-cgls3.max_iteration = 200
-cgls3.update_objective_interval = 5
-cgls3.run(200, verbose = True)
-
-# Show results
-plt.figure(figsize=(8,8))
-
-plt.subplot(2,2,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-
-plt.subplot(2,2,2)
-plt.imshow(cgls1.get_output().as_array())
-plt.title('GGLS')
-
-plt.subplot(2,2,3)
-plt.imshow(cgls2.get_output().as_array())
-plt.title('Regularised GGLS L = {} * Identity'.format(alpha1))
-
-plt.subplot(2,2,4)
-plt.imshow(cgls3.get_output().as_array())
-plt.title('Regularised GGLS L =  {} * Gradient'.format(alpha2))
-
-plt.show()
-
-
-
-print('Compare CVX vs Regularised CG with L = Gradient')
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-import astra
-import numpy
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-    
-    ##Construct problem    
-    u = Variable(N*N)
-    #q = Variable()
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-
-    regulariser = alpha2**2 * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    
-    # create matrix representation for Astra operator
-
-    vol_geom = astra.create_vol_geom(N, N)
-    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-
-    proj_id = astra.create_projector('line', proj_geom, vol_geom)
-
-    matrix_id = astra.projector.matrix(proj_id)
-
-    ProjMat = astra.matrix.get(matrix_id)
-    
-    fidelity = sum_squares( ProjMat * u - noisy_data.as_array().ravel()) 
-        
-    solver = MOSEK
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)    
-
-
-plt.figure(figsize=(10,20))
-
-plt.subplot(1,3,1)
-plt.imshow(np.reshape(u.value, (N, N)))    
-plt.colorbar()
-
-plt.subplot(1,3,2)
-plt.imshow(cgls3.get_output().as_array())    
-plt.colorbar()
-
-plt.subplot(1,3,3)
-plt.imshow(np.abs(cgls3.get_output().as_array() - np.reshape(u.value, (N, N)) ))    
-plt.colorbar()
-
-plt.show()
-
-print('Primal Objective (CVX) {} '.format(obj.value))
-print('Primal Objective (CGLS) {} '.format(cgls3.objective[-1]))
-
diff --git a/Wrappers/Python/demos/CGLS_examples/LinearSystem.py b/Wrappers/Python/demos/CGLS_examples/LinearSystem.py
deleted file mode 100644
index cc398cb..0000000
--- a/Wrappers/Python/demos/CGLS_examples/LinearSystem.py
+++ /dev/null
@@ -1,82 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-import numpy
-from ccpi.optimisation.operators import *
-from ccpi.optimisation.algorithms import *
-from ccpi.optimisation.functions import *
-from ccpi.framework import *
-
-# Problem dimension.
-m = 200
-n = 200
-
-numpy.random.seed(10)
-
-# Create matrix A and data b
-Amat = numpy.asarray( numpy.random.randn(m, n), dtype = numpy.float32)
-bmat = numpy.asarray( numpy.random.randn(m), dtype=numpy.float32)
-
-# Geometries for data and solution
-vgb = VectorGeometry(m)
-vgx = VectorGeometry(n)
-
-b = vgb.allocate(0, dtype=numpy.float32)
-b.fill(bmat)
-A = LinearOperatorMatrix(Amat)
-
-# Setup and run CGLS
-x_init = vgx.allocate()
-
-cgls = CGLS(x_init = x_init, operator = A, data = b)
-cgls.max_iteration = 2000
-cgls.update_objective_interval = 200
-cgls.run(2000, verbose = True)
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-# Compare with CVX
-x = Variable(n)
-obj = Minimize(sum_squares(A.A*x - b.as_array()))
-prob = Problem(obj)
-
-# choose solver
-if 'MOSEK' in installed_solvers():
-    solver = MOSEK
-else:
-    solver = SCS 
-
-result = prob.solve(solver = MOSEK)
-
-
-diff_sol = x.value - cgls.get_output().as_array()
- 
-print('Error |CVX - CGLS| = {}'.format(numpy.sum(numpy.abs(diff_sol))))
-print('CVX objective = {}'.format(obj.value))
-print('CGLS objective = {}'.format(cgls.objective[-1]))
-
-
-
-
diff --git a/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py b/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py
deleted file mode 100644
index c124fa1..0000000
--- a/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py
+++ /dev/null
@@ -1,146 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-Conjugate Gradient for (Regularized) Least Squares for Tomography
-
-
-Problem:     min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2}
-
-             A: Identity operator
-             g: Sinogram
-             L: Identity or Gradient Operator
-
-"""
-
-
-from ccpi.framework import ImageGeometry, ImageData, \
-                            AcquisitionGeometry, BlockDataContainer, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import CGLS
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
-from ccpi.framework import TestData
-       
-import os, sys
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-data = loader.load(TestData.SHAPES)
-ig = data.geometry
-ag = ig
-
-noisy_data = ImageData(TestData.random_noise(data.as_array(), mode = 'gaussian', seed = 1))
-#noisy_data = ImageData(data.as_array())
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-plt.subplot(2,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(2,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Setup and run the regularised CGLS algorithm  (Tikhonov with Gradient)
-x_init = ig.allocate() 
-alpha = 2
-op = Gradient(ig)
-
-block_op = BlockOperator( Identity(ig), alpha * op, shape=(2,1))
-block_data = BlockDataContainer(noisy_data, op.range_geometry().allocate())
-   
-cgls = CGLS(x_init=x_init, operator = block_op, data = block_data)
-cgls.max_iteration = 200
-cgls.update_objective_interval = 5
-cgls.run(200, verbose = True)
-
-# Show results
-plt.figure(figsize=(20,10))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy')
-plt.subplot(3,1,3)
-plt.imshow(cgls.get_output().as_array())
-plt.title('Regularised GGLS with Gradient')
-plt.show()
-
-#%%
-
-print('Compare CVX vs Regularised CG with L = Gradient')
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-    
-    ##Construct problem    
-    u = Variable(ig.shape)
-    #q = Variable()
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-
-    regulariser = alpha**2 * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-        
-    fidelity = sum_squares( u - noisy_data.as_array()) 
-        
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS 
-        
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = MOSEK)    
-
-
-plt.figure(figsize=(20,20))
-plt.subplot(3,1,1)
-plt.imshow(np.reshape(u.value, ig.shape))    
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(cgls.get_output().as_array())    
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(np.abs(cgls.get_output().as_array() - np.reshape(u.value, ig.shape) ))    
-plt.colorbar()
-plt.show()
-
-print('Primal Objective (CVX) {} '.format(obj.value))
-print('Primal Objective (CGLS) {} '.format(cgls.objective[-1]))
-
diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
deleted file mode 100644
index 09e350b..0000000
--- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
+++ /dev/null
@@ -1,199 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-Compare solutions of FISTA & PDHG & CGLS  
-                    & Astra Built-in algorithms for Least Squares
-
-
-Problem:     min_x || A x - g ||_{2}^{2}
-
-             A: Projection operator
-             g: Sinogram
-
-"""
-
-
-from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA
-
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition
-from ccpi.astra.ops import AstraProjectorSimple
-import astra   
-
-import tomophantom
-from tomophantom import TomoP2D
-import os
-
-# Load  Shepp-Logan phantom 
-model = 1 # select a model number from the library
-N = 128 # set dimension of the phantom
-path = os.path.dirname(tomophantom.__file__)
-path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
-phantom_2D = TomoP2D.Model(model, N, path_library2D)
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-data = ImageData(phantom_2D)
-
-detectors =  N
-angles = np.linspace(0, np.pi, 180, dtype=np.float32)
-
-ag = AcquisitionGeometry('parallel','2D', angles, detectors)
-
-device = input('Available device: GPU==1 / CPU==0 ')
-
-if device =='1':
-    dev = 'gpu'
-else:
-    dev = 'cpu'
-
-Aop = AstraProjectorSimple(ig, ag, dev)    
-sinogram = Aop.direct(data)
-
-
-
-###############################################################################
-# Setup and run Astra CGLS algorithm
-vol_geom = astra.create_vol_geom(N, N)
-proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-proj_id = astra.create_projector('line', proj_geom, vol_geom)
-
-# Create a sinogram id
-sinogram_id = astra.data2d.create('-sino', proj_geom, sinogram.as_array())
-
-# Create a data id
-rec_id = astra.data2d.create('-vol', vol_geom)
-
-cgls_astra = astra.astra_dict('CGLS')
-cgls_astra['ReconstructionDataId'] = rec_id
-cgls_astra['ProjectionDataId'] = sinogram_id
-cgls_astra['ProjectorId'] = proj_id
-
-# Create the algorithm object from the configuration structure
-alg_id = astra.algorithm.create(cgls_astra)
-
-astra.algorithm.run(alg_id, 500)
-recon_cgls_astra = ImageData(astra.data2d.get(rec_id))
-
-
-#%%
-
-# Setup and run the simple CGLS algorithm  
-x_init = ig.allocate()  
-
-cgls = CGLS(x_init = x_init, operator = Aop, data = sinogram)
-cgls.max_iteration = 500
-cgls.update_objective_interval = 100
-cgls.run(500, verbose = True)
-
-#%%
-
-###############################################################################
-# Setup and run the PDHG algorithm 
-operator = Aop
-f = L2NormSquared(b = sinogram)
-g = ZeroFunction()
-
-## Compute operator Norm
-normK = operator.norm()
-
-## Primal & dual stepsizes
-sigma = 0.02
-tau = 1/(sigma*normK**2)
-
-
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 1000
-pdhg.update_objective_interval = 100
-pdhg.run(1000, verbose=True)
-
-#%%
-###############################################################################
-# Setup and run the FISTA algorithm 
-fidelity = FunctionOperatorComposition(L2NormSquared(b=sinogram), Aop)
-regularizer = ZeroFunction()
-
-fista = FISTA(x_init=x_init , f=fidelity, g=regularizer)
-fista.max_iteration = 500
-fista.update_objective_interval = 100
-fista.run(500, verbose = True)
-
-#%% Show results
-
-plt.figure(figsize=(10,10))
-plt.suptitle('Reconstructions ', fontsize=16)
-
-plt.subplot(2,2,1)
-plt.imshow(cgls.get_output().as_array())
-plt.colorbar()
-plt.title('CGLS reconstruction')
-
-plt.subplot(2,2,2)
-plt.imshow(fista.get_output().as_array())
-plt.colorbar()
-plt.title('FISTA reconstruction')
-
-plt.subplot(2,2,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.colorbar()
-plt.title('PDHG reconstruction')
-
-plt.subplot(2,2,4)
-plt.imshow(recon_cgls_astra.as_array())
-plt.colorbar()
-plt.title('CGLS astra')
-
-diff1 = pdhg.get_output() - recon_cgls_astra
-diff2 = fista.get_output() - recon_cgls_astra
-diff3 = cgls.get_output() - recon_cgls_astra
-
-plt.figure(figsize=(15,15))
-
-plt.subplot(3,1,1)
-plt.imshow(diff1.abs().as_array())
-plt.title('Diff PDHG vs CGLS astra')
-plt.colorbar()
-
-plt.subplot(3,1,2)
-plt.imshow(diff2.abs().as_array())
-plt.title('Diff FISTA vs CGLS astra')
-plt.colorbar()
-
-plt.subplot(3,1,3)
-plt.imshow(diff3.abs().as_array())
-plt.title('Diff CLGS vs CGLS astra')
-plt.colorbar()
-
-
-cgls_astra_obj = fidelity(ImageData(recon_cgls_astra))
-
-print('Primal Objective (FISTA) {} '.format(fista.objective[-1]))
-print('Primal Objective (CGLS) {} '.format(cgls.objective[-1]))
-print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
-print('Primal Objective (CGLS_astra) {} '.format(cgls_astra_obj))
-
-
-
diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py
deleted file mode 100644
index 9b6d10f..0000000
--- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py
+++ /dev/null
@@ -1,133 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-Compare solutions of PDHG & "Block CGLS" algorithms for 
-
-
-Problem:     min_x alpha * ||\grad x ||^{2}_{2} + || A x - g ||_{2}^{2}
-
-
-             A: Projection operator
-             g: Sinogram
-
-"""
-
-
-from ccpi.framework import AcquisitionGeometry, BlockDataContainer, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, CGLS
-from ccpi.optimisation.operators import BlockOperator, Gradient
-
-from ccpi.optimisation.functions import ZeroFunction, BlockFunction, L2NormSquared  
-from ccpi.astra.ops import AstraProjectorSimple    
-from ccpi.framework import TestData
-import os, sys
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-
-# Create Ground truth phantom and Sinogram                 
-N = 150
-M = 150
-data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
-ig = data.geometry
-
-detectors = N
-angles = np.linspace(0, np.pi, N, dtype=np.float32)
-ag = AcquisitionGeometry('parallel','2D', angles, detectors)
-
-device = input('Available device: GPU==1 / CPU==0 ')
-if device=='1':
-    dev = 'gpu'
-else:
-    dev = 'cpu'
-    
-Aop = AstraProjectorSimple(ig, ag, dev)    
-sin = Aop.direct(data)
-
-noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,3,ig.shape))
-
-# Setup and run the CGLS algorithm  
-alpha = 50
-Grad = Gradient(ig)
-
-# Form Tikhonov as a Block CGLS structure
-op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1))
-block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate())
-
-x_init = ig.allocate()      
-cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data)
-cgls.max_iteration = 1000
-cgls.update_objective_interval = 200
-cgls.run(1000,verbose=False)
-
-
-#Setup and run the PDHG algorithm 
-
-# Create BlockOperator
-op_PDHG = BlockOperator(Grad, Aop, shape=(2,1) ) 
-# Create functions     
-f1 = 0.5 * alpha**2 * L2NormSquared()
-f2 = 0.5 * L2NormSquared(b = noisy_data)    
-f = BlockFunction(f1, f2)                                       
-g = ZeroFunction()
-
-## Compute operator Norm
-normK = op_PDHG.norm()
-
-## Primal & dual stepsizes
-sigma = 10
-tau = 1/(sigma*normK**2)
-
-pdhg = PDHG(f=f,g=g,operator=op_PDHG, tau=tau, sigma=sigma)
-pdhg.max_iteration = 1000
-pdhg.update_objective_interval = 200
-pdhg.run(1000, verbose=False)
-
-# Show results
-plt.figure(figsize=(10,10))
-
-plt.subplot(2,1,1)
-plt.imshow(cgls.get_output().as_array())
-plt.title('CGLS reconstruction')
-
-plt.subplot(2,1,2)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('PDHG reconstruction')
-
-plt.show()
-
-diff1 = pdhg.get_output() - cgls.get_output()
-
-plt.imshow(diff1.abs().as_array())
-plt.title('Diff PDHG vs CGLS')
-plt.colorbar()
-plt.show()
-
-plt.plot(np.linspace(0,N,M), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-plt.plot(np.linspace(0,N,M), cgls.get_output().as_array()[int(N/2),:], label = 'CGLS')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
diff --git a/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py b/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py
deleted file mode 100644
index c24ebac..0000000
--- a/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py
+++ /dev/null
@@ -1,124 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-
-""" 
-Compare solutions of FISTA & PDHG algorithms
-
-
-Problem:     min_x alpha * ||\grad x ||^{2}_{2} + || x - g ||_{1}
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator
-             
-             g: Noisy Data with Salt & Pepper Noise
-
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import FISTA, PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
-from ccpi.optimisation.functions import L2NormSquared, L1Norm, \
-                                         FunctionOperatorComposition, BlockFunction, ZeroFunction
-                                                
-from skimage.util import random_noise
-
-
-# Create Ground truth and Noisy data
-N = 100
-
-data = np.zeros((N,N))
-data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-data = ImageData(data)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 5
-
-###############################################################################
-# Setup and run the FISTA algorithm
-operator = Gradient(ig)
-fidelity = L1Norm(b=noisy_data)
-regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
-
-x_init = ig.allocate()
-opt = {'memopt':True}
-fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt)
-fista.max_iteration = 2000
-fista.update_objective_interval = 50
-fista.run(2000, verbose=False)
-###############################################################################
-
-
-###############################################################################
-# Setup and run the PDHG algorithm
-op1 = Gradient(ig)
-op2 = Identity(ig, ag)
-
-operator = BlockOperator(op1, op2, shape=(2,1) )   
-f = BlockFunction(alpha * L2NormSquared(), fidelity)                                        
-g = ZeroFunction()
-  
-normK = operator.norm()
-
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000, verbose=False)
-###############################################################################
-
-# Show results
-
-plt.figure(figsize=(10,10))
-
-plt.subplot(2,1,1)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('PDHG reconstruction')
-
-plt.subplot(2,1,2)
-plt.imshow(fista.get_output().as_array())
-plt.title('FISTA reconstruction')
-
-plt.show()
-
-diff1 = pdhg.get_output() - fista.get_output()
-
-plt.imshow(diff1.abs().as_array())
-plt.title('Diff PDHG vs FISTA')
-plt.colorbar()
-plt.show()
-
-
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
deleted file mode 100644
index 1a96a2d..0000000
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
+++ /dev/null
@@ -1,215 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.framework import AcquisitionGeometry
-from ccpi.optimisation.algorithms import FISTA
-from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, \
-                         L2NormSquared, FunctionOperatorComposition
-from ccpi.astra.operators import AstraProjectorSimple
-
-import numpy as np
-import matplotlib.pyplot as plt
-from ccpi.framework import TestData
-import os, sys
-from ccpi.optimisation.funcs import Norm2sq
-
-# Load Data 
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-                     
-N = 50
-M = 50
-data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
-
-ig = data.geometry
-
-# Show Ground Truth 
-plt.figure(figsize=(5,5))
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-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).
-
-#device = input('Available device: GPU==1 / CPU==0 ')
-
-#if device=='1':
-#    dev = 'gpu'
-#else:
-#    dev = 'cpu'
-
-test_case = 1
-dev = 'cpu'
-
-if test_case==1:
-    
-    detectors = N
-    angles_num = N    
-    det_w = 1.0
-    
-    angles = np.linspace(0, np.pi, angles_num, endpoint=False)
-    ag = AcquisitionGeometry('parallel',
-                             '2D',
-                             angles,
-                             detectors,det_w)
-  
-elif test_case==2:
-    
-    SourceOrig = 200
-    OrigDetec = 0
-    angles = np.linspace(0,2*np.pi,angles_num)
-    ag = AcquisitionGeometry('cone',
-                             '2D',
-                             angles,
-                             detectors,
-                             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, dev)
-    
-# Forward and backprojection are available as methods direct and adjoint. Here 
-# generate test data b and do simple backprojection to obtain z.
-sin = Aop.direct(data)
-back_proj = Aop.adjoint(sin)
-
-plt.figure(figsize=(5,5))
-plt.imshow(sin.array)
-plt.title('Simulated data')
-plt.show()
-
-plt.figure(figsize=(5,5))
-plt.imshow(back_proj.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:
-
-f = FunctionOperatorComposition(L2NormSquared(b=sin), Aop)
-#f = Norm2sq(Aop, sin, c=0.5, memopt=True)
-g = ZeroFunction()
-
-x_init = ig.allocate()
-fista = FISTA(x_init=x_init, f=f, g=g)
-fista.max_iteration = 1000
-fista.update_objective_interval = 100
-fista.run(1000, verbose=True)
-
-plt.figure()
-plt.imshow(fista.get_output().as_array())
-plt.title('FISTA unconstrained')
-plt.colorbar()
-plt.show()
-
-
-# Run FISTA for least squares with lower/upper bound 
-fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1))
-fista0.max_iteration = 1000
-fista0.update_objective_interval = 100
-fista0.run(1000, verbose=True)
-
-plt.figure()
-plt.imshow(fista0.get_output().as_array())
-plt.title('FISTA constrained in [0,1]')
-plt.colorbar()
-plt.show()
-
-#plt.figure()
-#plt.semilogy(fista0.objective)
-#plt.title('FISTA constrained in [0,1]')
-#plt.show()
-
-#%% Check with CVX solution
-
-import astra
-import numpy
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-    
-if cvx_not_installable:
-    
-    ##Construct problem    
-    u = Variable(N*N)
-
-    # create matrix representation for Astra operator
-    vol_geom = astra.create_vol_geom(N, N)
-    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-
-    proj_id = astra.create_projector('linear', proj_geom, vol_geom)
-
-    matrix_id = astra.projector.matrix(proj_id)
-    
-    ProjMat = astra.matrix.get(matrix_id)
-    
-    tmp = sin.as_array().ravel()
-    
-    fidelity = 0.5 * sum_squares(ProjMat * u - tmp)
-
-    solver = MOSEK
-    obj =  Minimize(fidelity)
-    constraints = [u>=0, u<=1]
-    prob = Problem(obj, constraints=constraints)
-    result = prob.solve(verbose = True, solver = solver)   
-         
-    diff_cvx = numpy.abs( fista0.get_output().as_array() - np.reshape(u.value, (N,M) ))
-           
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(fista0.get_output().as_array())
-    plt.title('FISTA solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(np.reshape(u.value, (N, M)))
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,N,M), fista0.get_output().as_array()[int(N/2),:], label = 'FISTA')
-    plt.plot(np.linspace(0,N,M), np.reshape(u.value, (N,M) )[int(N/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (FISTA) {} '.format(fista0.loss[1]))    
\ No newline at end of file
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
deleted file mode 100644
index 6007990..0000000
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
+++ /dev/null
@@ -1,201 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-
-"Tikhonov regularization" for Poisson denoising using FISTA algorithm:
-
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{2}^{2} + \int x - g * log(x) 
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator  
-             
-             g: Noisy Data with Poisson Noise
-               
-                       
-"""
-
-from ccpi.framework import ImageData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import FISTA
-
-from ccpi.optimisation.operators import Gradient
-from ccpi.optimisation.functions import KullbackLeibler, L2NormSquared, FunctionOperatorComposition
-
-from ccpi.framework import TestData
-import os, sys
-from skimage.util import random_noise
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-
-# Load Data                      
-N = 100
-M = 100
-data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
-
-ig = data.geometry
-ag = ig
-
-# Create Noisy data with Poisson noise
-n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
-noisy_data = ImageData(n1)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-plt.subplot(2,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(2,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Regularisation Parameter
-alpha = 10
-
-# Setup and run the FISTA algorithm
-operator = Gradient(ig)
-fid = KullbackLeibler(noisy_data)
-
-def KL_Prox_PosCone(x, tau, out=None):
-        
-    if out is None: 
-        tmp = 0.5 *( (x - fid.bnoise - tau) + ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b   ) .sqrt() )
-        return tmp.maximum(0)
-    else:            
-        tmp =  0.5 *( (x - fid.bnoise - tau) + 
-                    ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b   ) .sqrt()
-                    )
-        x.add(fid.bnoise, out=out)
-        out -= tau
-        out *= out
-        tmp = fid.b * (4 * tau)
-        out.add(tmp, out=out)
-        out.sqrt(out=out)
-        
-        x.subtract(fid.bnoise, out=tmp)
-        tmp -= tau
-        
-        out += tmp
-        
-        out *= 0.5
-        
-        # ADD the constraint here
-        out.maximum(0, out=out)
-        
-fid.proximal = KL_Prox_PosCone
-
-reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
-
-x_init = ig.allocate()
-fista = FISTA(x_init=x_init , f=reg, g=fid)
-fista.max_iteration = 2000
-fista.update_objective_interval = 500
-fista.run(2000, verbose=True)
-
-# Show results
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(fista.get_output().as_array())
-plt.title('Reconstruction')
-plt.colorbar()
-plt.show()
-
-plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,M), fista.get_output().as_array()[int(N/2),:], label = 'Reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-#%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u1 = Variable(ig.shape)
-    q = Variable()
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0))    
-    fidelity = sum(kl_div(noisy_data.as_array(), u1))  
-    
-    constraints = [q>=fidelity, u1>=0]    
-    
-    solver = SCS
-    obj =  Minimize( regulariser +  q)
-    prob = Problem(obj, constraints)
-    result = prob.solve(verbose = True, solver = solver)
-    
-    diff_cvx = numpy.abs( fista.get_output().as_array() - u1.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(fista.get_output().as_array())
-    plt.title('FISTA solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u1.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,N,M), fista.get_output().as_array()[int(N/2),:], label = 'FISTA')
-    plt.plot(np.linspace(0,N,M), u1.value[int(N/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (FISTA) {} '.format(fista.loss[1]))
-
-
-
diff --git a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
deleted file mode 100644
index e69060f..0000000
--- a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
+++ /dev/null
@@ -1,273 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.framework import ImageGeometry, ImageData, AcquisitionGeometry, AcquisitionData, BlockDataContainer
-
-import numpy as numpy                 
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, CGLS
-from ccpi.optimisation.algs import CGLS as CGLS_old
-
-from ccpi.optimisation.operators import BlockOperator, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.operators import AstraProjectorMC
-from scipy.io import loadmat
-import h5py
-
-#%%
-
-phantom = 'powder'
-
-if phantom == 'carbon':
-    pathname = '/media/newhd/shared/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))
-    X = X[80:100] # delete this to take all channels
-    
-elif phantom == 'powder':
-    pathname = '/media/newhd/shared/DataProcessed/' 
-    filename = 'S_180.mat'
-    path = pathname + filename
-    arrays = {}
-    f = h5py.File(path)
-    for k, v in f.items():
-        arrays[k] = numpy.array(v)
-    XX = arrays['S']    
-    X = numpy.transpose(XX,(0,2,1,3))
-    X = X[100:120]
-    
-    
-#%% Setup Geometry of Colorbay
-
-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[5,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=5).array)
-plt.show()
-
-#%% CGLS
-
-def callback(iteration, objective, x):
-    plt.imshow(x.as_array()[5])
-    plt.colorbar()
-    plt.show()
-
-x_init = ig2d.allocate()      
-cgls1 = CGLS(x_init=x_init, operator=Aall, data=data2d)
-cgls1.max_iteration = 100
-cgls1.update_objective_interval = 1
-cgls1.run(5,verbose=True, callback = callback)
-
-plt.imshow(cgls1.get_output().subset(channel=5).array)
-plt.title('CGLS')
-plt.show()
-
-#%% Tikhonov 
-
-alpha = 2.5
-Grad = Gradient(ig2d, correlation=Gradient.CORRELATION_SPACE) # use also CORRELATION_SPACECHANNEL
-
-# Form Tikhonov as a Block CGLS structure
-op_CGLS = BlockOperator( Aall, alpha * Grad, shape=(2,1))
-block_data = BlockDataContainer(data2d, Grad.range_geometry().allocate())
-    
-cgls2 = CGLS(x_init=x_init, operator=op_CGLS, data=block_data)
-cgls2.max_iteration = 100
-cgls2.update_objective_interval = 1
-
-cgls2.run(10,verbose=True, callback=callback)
-
-plt.imshow(cgls2.get_output().subset(channel=5).array)
-plt.title('Tikhonov')
-plt.show()
-
-#%% Total Variation
-
-# Regularisation Parameter
-#alpha_TV = 0.08 # for carbon
-alpha_TV = 0.08 # for powder
-
-# Create operators
-op1 = Gradient(ig2d, correlation=Gradient.CORRELATION_SPACE)
-op2 = Aall
-
-# Create BlockOperator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-# Create functions      
-f1 = alpha_TV * MixedL21Norm()
-f2 = 0.5 * L2NormSquared(b=data2d)    
-f = BlockFunction(f1, f2)  
-g = ZeroFunction()
-
-# Compute operator Norm
-#normK = 8.70320267279591 # For powder Run one time no need to compute again takes time
-normK = 14.60320657253632 # for carbon
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 100
-pdhg.run(1000, verbose =True, callback=callback)
-
-
-#%% Show sinograms
-channel_ind = [10,15,19]
-
-plt.figure(figsize=(15,15))
-
-plt.subplot(4,3,1)
-plt.imshow(data2d.subset(channel = channel_ind[0]).as_array())
-plt.title('Channel {}'.format(channel_ind[0]))
-plt.colorbar()
-
-plt.subplot(4,3,2)
-plt.imshow(data2d.subset(channel = channel_ind[1]).as_array())
-plt.title('Channel {}'.format(channel_ind[1]))
-plt.colorbar()
-
-plt.subplot(4,3,3)
-plt.imshow(data2d.subset(channel = channel_ind[2]).as_array())
-plt.title('Channel {}'.format(channel_ind[2]))
-plt.colorbar()
-
-###############################################################################
-# Show CGLS
-plt.subplot(4,3,4)
-plt.imshow(cgls1.get_output().subset(channel = channel_ind[0]).as_array())
-plt.colorbar()
-
-plt.subplot(4,3,5)
-plt.imshow(cgls1.get_output().subset(channel = channel_ind[1]).as_array())
-plt.colorbar()
-
-plt.subplot(4,3,6)
-plt.imshow(cgls1.get_output().subset(channel = channel_ind[2]).as_array())
-plt.colorbar()
-
-###############################################################################
-# Show Tikhonov
-
-plt.subplot(4,3,7)
-plt.imshow(cgls2.get_output().subset(channel = channel_ind[0]).as_array())
-plt.colorbar()
-
-plt.subplot(4,3,8)
-plt.imshow(cgls2.get_output().subset(channel = channel_ind[1]).as_array())
-plt.colorbar()
-
-plt.subplot(4,3,9)
-plt.imshow(cgls2.get_output().subset(channel = channel_ind[2]).as_array())
-plt.colorbar()
-
-
-###############################################################################
-# Show Total variation
-
-plt.subplot(4,3,10)
-plt.imshow(pdhg.get_output().subset(channel = channel_ind[0]).as_array())
-plt.colorbar()
-
-plt.subplot(4,3,11)
-plt.imshow(pdhg.get_output().subset(channel = channel_ind[1]).as_array())
-plt.colorbar()
-
-plt.subplot(4,3,12)
-plt.imshow(pdhg.get_output().subset(channel = channel_ind[2]).as_array())
-plt.colorbar()
-
-
-###############################################################################
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
deleted file mode 100755
index 9dbcf3e..0000000
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
+++ /dev/null
@@ -1,282 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-""" 
-
-Total Generalised Variation (TGV) Denoising using PDHG algorithm:
-
-
-Problem:     min_{u} \alpha * ||\nabla u - w||_{2,1} +
-                     \beta *  || E u ||_{2,1} +
-                     Fidelity(u, g)
-                     
-             \nabla: Gradient operator 
-             E: Symmetrized Gradient operator             
-             \alpha: Regularization parameter
-             \beta:  Regularization parameter             
-             
-             g: Noisy Data 
-                          
-             Fidelity =  1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian
-                         2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper
-                         3) Kullback Leibler (\int u - g * log(u) + Id_{u>0})  if Noise is Poisson
-                                
-                          
-             Method = 0 ( PDHG - split ) :  K = [ \nabla, - Identity
-                                                  ZeroOperator, E 
-                                                  Identity, ZeroOperator]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  K = [ \nabla, - Identity
-                                                  ZeroOperator, E ] 
-                          
-             Default: TGV denoising
-             noise = Gaussian
-             Fidelity = L2NormSquarred 
-             method = 0             
-                                                                
-"""
-
-from ccpi.framework import ImageData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, \
-                        Gradient, SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
-                      MixedL21Norm, BlockFunction, KullbackLeibler, L2NormSquared
-                      
-from ccpi.framework import TestData
-import os, sys
-if int(numpy.version.version.split('.')[1]) > 12:
-    from skimage.util import random_noise
-else:
-    from demoutil import random_noise
-
-# user supplied input
-if len(sys.argv) > 1:
-    which_noise = int(sys.argv[1])
-else:
-    which_noise = 0
-print ("Applying {} noise")
-
-if len(sys.argv) > 2:
-    method = sys.argv[2]
-else:
-    method = '0'
-print ("method ", method)
-
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-data = loader.load(TestData.SHAPES)
-ig = data.geometry
-ag = ig
-
-# Create noisy data. 
-noises = ['gaussian', 'poisson', 's&p']
-noise = noises[which_noise]
-if noise == 's&p':
-    n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2, seed=10)
-elif noise == 'poisson':
-    scale = 5
-    n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale
-elif noise == 'gaussian':
-    n1 = random_noise(data.as_array(), mode = noise, seed = 10)
-else:
-    raise ValueError('Unsupported Noise ', noise)
-noisy_data = ImageData(n1)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,5))
-plt.subplot(1,2,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,2,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Regularisation Parameter depending on the noise distribution
-if noise == 's&p':
-    alpha = 0.8
-elif noise == 'poisson':
-    alpha = .3
-elif noise == 'gaussian':
-    alpha = .2
-
-# TODO add ref why this choice
-beta = 2 * alpha
-
-# Fidelity
-if noise == 's&p':
-    f3 = L1Norm(b=noisy_data)
-elif noise == 'poisson':
-    f3 = KullbackLeibler(noisy_data)
-elif noise == 'gaussian':
-    f3 = 0.5 * L2NormSquared(b=noisy_data)
-
-if method == '0':
-    
-    # Create operators
-    op11 = Gradient(ig)
-    op12 = Identity(op11.range_geometry())
-    
-    op22 = SymmetrizedGradient(op11.domain_geometry())    
-    op21 = ZeroOperator(ig, op22.range_geometry())
-        
-    op31 = Identity(ig, ag)
-    op32 = ZeroOperator(op22.domain_geometry(), ag)
-    
-    operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) 
-        
-    f1 = alpha * MixedL21Norm()
-    f2 = beta * MixedL21Norm() 
-    
-    f = BlockFunction(f1, f2, f3)         
-    g = ZeroFunction()
-        
-else:
-    
-    # Create operators
-    op11 = Gradient(ig)
-    op12 = Identity(op11.range_geometry())
-    op22 = SymmetrizedGradient(op11.domain_geometry())    
-    op21 = ZeroOperator(ig, op22.range_geometry())    
-    
-    operator = BlockOperator(op11, -1*op12, op21, op22, shape=(2,2) )      
-    
-    f1 = alpha * MixedL21Norm()
-    f2 = beta * MixedL21Norm()     
-    
-    f = BlockFunction(f1, f2)         
-    g = BlockFunction(f3, ZeroFunction())
-     
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 100
-pdhg.run(2000)
-
-# Show results
-plt.figure(figsize=(20,5))
-plt.subplot(1,4,1)
-plt.imshow(data.subset(channel=0).as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,4,2)
-plt.imshow(noisy_data.subset(channel=0).as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(1,4,3)
-plt.imshow(pdhg.get_output()[0].as_array())
-plt.title('TGV Reconstruction')
-plt.colorbar()
-plt.subplot(1,4,4)
-plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'TGV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-#%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-if cvx_not_installable:    
-    
-    u = Variable(ig.shape)
-    w1 = Variable(ig.shape)
-    w2 = Variable(ig.shape)
-    
-    # create TGV regulariser
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \
-                                           DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \
-                  beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \
-                                      0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \
-                                      0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0  ) )  
-    
-    constraints = []
-    
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS      
-        
-    # fidelity
-    if noise == 's&p':
-        fidelity = pnorm( u - noisy_data.as_array(),1)
-    elif noise == 'poisson':
-        fidelity = sum(kl_div(noisy_data.as_array(), u)) 
-        solver = SCS
-    elif noise == 'gaussian':
-        fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
-        
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)
-    
-    diff_cvx = numpy.abs( pdhg.get_output()[0].as_array() - u.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output()[0].as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
deleted file mode 100755
index 6937fa0..0000000
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
+++ /dev/null
@@ -1,280 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-
-Problem:     min_{u},   \alpha * ||\nabla u||_{2,1} + Fidelity(u, g)
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: Noisy Data 
-                          
-             Fidelity =  1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian
-                         2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper
-                         3) Kullback Leibler (\int u - g * log(u) + Id_{u>0})  if Noise is Poisson
-                                                       
-             Method = 0 ( PDHG - split ) :  K = [ \nabla,
-                                                 Identity]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  K = \nabla    
-             
-             
-             Default: ROF denoising
-             noise = Gaussian
-             Fidelity = L2NormSquarred 
-             method = 0
-             
-             
-"""
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
-                      MixedL21Norm, BlockFunction, L2NormSquared,\
-                          KullbackLeibler
-from ccpi.framework import TestData
-import os, sys
-#import scipy.io
-if int(numpy.version.version.split('.')[1]) > 12:
-    from skimage.util import random_noise
-else:
-    from demoutil import random_noise
-
-# user supplied input
-if len(sys.argv) > 1:
-    which_noise = int(sys.argv[1])
-else:
-    which_noise = 0
-print ("Applying {} noise")
-
-if len(sys.argv) > 2:
-    method = sys.argv[2]
-else:
-    method = '0'
-print ("method ", method)
-
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-data = loader.load(TestData.SHAPES, size=(50,50))
-ig = data.geometry
-ag = ig
-
-# Create noisy data. 
-noises = ['gaussian', 'poisson', 's&p']
-noise = noises[which_noise]
-if noise == 's&p':
-    n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2)
-elif noise == 'poisson':
-    scale = 5
-    n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale
-elif noise == 'gaussian':
-    n1 = random_noise(data.as_array(), mode = noise, seed = 10)
-else:
-    raise ValueError('Unsupported Noise ', noise)
-noisy_data = ig.allocate()
-noisy_data.fill(n1)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,5))
-plt.subplot(1,2,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,2,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-
-# Regularisation Parameter depending on the noise distribution
-if noise == 's&p':
-    alpha = 0.8
-elif noise == 'poisson':
-    alpha = 1
-elif noise == 'gaussian':
-    alpha = .3
-
-# fidelity
-if noise == 's&p':
-    f2 = L1Norm(b=noisy_data)
-elif noise == 'poisson':
-    f2 = KullbackLeibler(noisy_data)
-elif noise == 'gaussian':
-    f2 = 0.5 * L2NormSquared(b=noisy_data)
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACE)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions      
-    f = BlockFunction(alpha * MixedL21Norm(), f2) 
-    g = ZeroFunction()
-    
-else:
-    
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g = f2
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 100
-pdhg.run(2000)
-
-
-if data.geometry.channels > 1:
-    plt.figure(figsize=(20,15))
-    for row in range(data.geometry.channels):
-        
-        plt.subplot(3,4,1+row*4)
-        plt.imshow(data.subset(channel=row).as_array())
-        plt.title('Ground Truth')
-        plt.colorbar()
-        plt.subplot(3,4,2+row*4)
-        plt.imshow(noisy_data.subset(channel=row).as_array())
-        plt.title('Noisy Data')
-        plt.colorbar()
-        plt.subplot(3,4,3+row*4)
-        plt.imshow(pdhg.get_output().subset(channel=row).as_array())
-        plt.title('TV Reconstruction')
-        plt.colorbar()
-        plt.subplot(3,4,4+row*4)
-        plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.subset(channel=row).as_array()[int(N/2),:], label = 'GTruth')
-        plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().subset(channel=row).as_array()[int(N/2),:], label = 'TV reconstruction')
-        plt.legend()
-        plt.title('Middle Line Profiles')
-    plt.show()
-    
-else:
-    plt.figure(figsize=(20,5))
-    plt.subplot(1,4,1)
-    plt.imshow(data.subset(channel=0).as_array())
-    plt.title('Ground Truth')
-    plt.colorbar()
-    plt.subplot(1,4,2)
-    plt.imshow(noisy_data.subset(channel=0).as_array())
-    plt.title('Noisy Data')
-    plt.colorbar()
-    plt.subplot(1,4,3)
-    plt.imshow(pdhg.get_output().subset(channel=0).as_array())
-    plt.title('TV Reconstruction')
-    plt.colorbar()
-    plt.subplot(1,4,4)
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth')
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'TV reconstruction')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-
-
-#%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u = Variable(ig.shape)
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([Constant(DX.matrix()) * vec(u), Constant(DY.matrix()) * vec(u)]), 2, axis = 0))
-    
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS      
-
-    # fidelity
-    if noise == 's&p':
-        fidelity = pnorm( u - noisy_data.as_array(),1)
-    elif noise == 'poisson':
-        fidelity = sum(kl_div(noisy_data.as_array(), u)) 
-        solver = SCS
-    elif noise == 'gaussian':
-        fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
-                
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py
deleted file mode 100644
index febe76d..0000000
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py
+++ /dev/null
@@ -1,243 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-""" 
-
-Total Variation (Dynamic) Denoising using PDHG algorithm and Tomophantom:
-
-
-Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: 2D Dynamic noisy data with Gaussian Noise
-                          
-             K = \nabla  
-                                                                
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared,  \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorMC
-
-import os
-import tomophantom
-from tomophantom import TomoP2D
-
-# Create phantom for TV 2D dynamic tomography 
-
-model = 102  # note that the selected model is temporal (2D + time)
-N = 128 # set dimension of the phantom
-
-path = os.path.dirname(tomophantom.__file__)
-path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
-phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D)
-
-plt.close('all')
-plt.figure(1)
-plt.rcParams.update({'font.size': 21})
-plt.title('{}''{}'.format('2D+t phantom using model no.',model))
-for sl in range(0,np.shape(phantom_2Dt)[0]):
-    im = phantom_2Dt[sl,:,:]
-    plt.imshow(im, vmin=0, vmax=1)
-#    plt.pause(.1)
-#    plt.draw
-
-# Setup geometries    
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
-data = ImageData(phantom_2Dt, geometry=ig)
-ag = ig
-                      
-# Create noisy data. Apply Gaussian noise
-np.random.seed(10)
-noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.25, size=ig.shape) )
-
-# time-frames index
-tindex = [8, 16, 24]
-
-fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10))
-plt.subplot(1,3,1)
-plt.imshow(noisy_data.as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-plt.subplot(1,3,2)
-plt.imshow(noisy_data.as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-plt.subplot(1,3,3)
-plt.imshow(noisy_data.as_array()[tindex[2],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-plt.show()   
-
-# Regularisation Parameter
-alpha = 0.3
-
-# Create operators
-op1 = Gradient(ig, correlation='Space')
-op2 = Gradient(ig, correlation='SpaceChannels')
-op3 = Identity(ig, ag)
-
-# Create BlockOperator
-operator1 = BlockOperator(op1, op3, shape=(2,1) ) 
-operator2 = BlockOperator(op2, op3, shape=(2,1) )
-
-# Create functions
-      
-f1 = alpha * MixedL21Norm()
-f2 = 0.5 * L2NormSquared(b = noisy_data)    
-f = BlockFunction(f1, f2)  
-                                                                        
-g = ZeroFunction()
-    
-# Compute operator Norm
-normK1 = operator1.norm()
-normK2 = operator2.norm()
-
-#%%
-# Primal & dual stepsizes
-sigma1 = 1
-tau1 = 1/(sigma1*normK1**2)
-
-sigma2 = 1
-tau2 = 1/(sigma2*normK2**2)
-
-# Setup and run the PDHG algorithm
-pdhg1 = PDHG(f=f,g=g,operator=operator1, tau=tau1, sigma=sigma1)
-pdhg1.max_iteration = 2000
-pdhg1.update_objective_interval = 200
-pdhg1.run(2000)
-
-# Setup and run the PDHG algorithm
-pdhg2 = PDHG(f=f,g=g,operator=operator2, tau=tau2, sigma=sigma2)
-pdhg2.max_iteration = 2000
-pdhg2.update_objective_interval = 200
-pdhg2.run(2000)
-
-
-#%%
-
-tindex = [8, 16, 24]
-fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
-
-plt.subplot(3,3,1)
-plt.imshow(phantom_2Dt[tindex[0],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-
-plt.subplot(3,3,2)
-plt.imshow(phantom_2Dt[tindex[1],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-
-plt.subplot(3,3,3)
-plt.imshow(phantom_2Dt[tindex[2],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-plt.subplot(3,3,4)
-plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.subplot(3,3,5)
-plt.imshow(pdhg1.get_output().as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.subplot(3,3,6)
-plt.imshow(pdhg1.get_output().as_array()[tindex[2],:,:])
-plt.axis('off')
-
-
-plt.subplot(3,3,7)
-plt.imshow(pdhg2.get_output().as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.subplot(3,3,8)
-plt.imshow(pdhg2.get_output().as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.subplot(3,3,9)
-plt.imshow(pdhg2.get_output().as_array()[tindex[2],:,:])
-plt.axis('off')
-
-
-im = plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:])
-
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
-cbar = fig.colorbar(im, cax=cb_ax)
-plt.show()
-
-#%%
-import matplotlib.animation as animation
-fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 30))
-ims1 = []
-ims2 = []
-ims3 = []
-for sl in range(0,np.shape(phantom_2Dt)[0]):
-    
-    plt.subplot(1,3,1)
-    im1 = plt.imshow(phantom_2Dt[sl,:,:], animated=True)
-    
-    plt.subplot(1,3,2)
-    im2 = plt.imshow(pdhg1.get_output().as_array()[sl,:,:])
-    
-    plt.subplot(1,3,3)
-    im3 = plt.imshow(pdhg2.get_output().as_array()[sl,:,:])
-    
-    ims1.append([im1])
-    ims2.append([im2])
-    ims3.append([im3])
-    
-    
-ani1 = animation.ArtistAnimation(fig, ims1, interval=500,
-                                repeat_delay=10)
-
-ani2 = animation.ArtistAnimation(fig, ims2, interval=500,
-                                repeat_delay=10)
-
-ani3 = animation.ArtistAnimation(fig, ims3, interval=500,
-                                repeat_delay=10)
-plt.show()        
-#    plt.pause(0.25)
-#    plt.show()
-
-
-
-
-
-
-
-
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py
deleted file mode 100644
index 15709cd..0000000
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py
+++ /dev/null
@@ -1,147 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-""" 
-
-Total Variation (3D) Denoising using PDHG algorithm and Tomophantom:
-
-
-Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: 3D Noisy Data with Gaussian Noise
-                          
-             Method = 0 ( PDHG - split ) :  K = [ \nabla,
-                                                 Identity]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  K = \nabla  
-                                                                
-"""
-
-from ccpi.framework import ImageData, ImageGeometry                       
-import matplotlib.pyplot as plt
-from ccpi.optimisation.algorithms import PDHG
-from ccpi.optimisation.operators import Gradient
-from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm
-
-from skimage.util import random_noise
-
-import timeit
-import os
-from tomophantom import TomoP3D
-import tomophantom
-
-# Create a phantom from Tomophantom
-print ("Building 3D phantom using TomoPhantom software")
-tic=timeit.default_timer()
-model = 13 # select a model number from the library
-N = 64 # Define phantom dimensions using a scalar value (cubic phantom)
-path = os.path.dirname(tomophantom.__file__)
-path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
-
-#This will generate a N x N x N phantom (3D)
-phantom_tm = TomoP3D.Model(model, N, path_library3D)
-
-# Create noisy data. Apply Gaussian noise
-ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N)
-ag = ig
-n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10)
-noisy_data = ImageData(n1)
-
-# Show results
-sliceSel = int(0.5*N)
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1)
-plt.title('Axial View')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1)
-plt.title('Coronal View')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1)
-plt.title('Sagittal View')
-plt.colorbar()
-plt.show()   
-
-# Regularisation Parameter
-alpha = 0.05
-    
-# Setup and run the PDHG algorithm
-operator = Gradient(ig)
-f =  alpha * MixedL21Norm()
-g =  0.5 * L2NormSquared(b = noisy_data)
-            
-normK = operator.norm()
-
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 1000
-pdhg.update_objective_interval = 200
-pdhg.run(1000, verbose = True)
-
-# Show results
-fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
-fig.suptitle('TV Reconstruction',fontsize=20)
-
-plt.subplot(2,3,1)
-plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Axial View')
-
-plt.subplot(2,3,2)
-plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Coronal View')
-
-plt.subplot(2,3,3)
-plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Sagittal View')
-
-
-plt.subplot(2,3,4)
-plt.imshow(pdhg.get_output().as_array()[sliceSel,:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.subplot(2,3,5)
-plt.imshow(pdhg.get_output().as_array()[:,sliceSel,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.subplot(2,3,6)
-plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1)
-plt.axis('off')
-im = plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1)
-
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
-cbar = fig.colorbar(im, cax=cb_ax)
-
-
-plt.show()
-
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
deleted file mode 100644
index 4f7639e..0000000
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
+++ /dev/null
@@ -1,173 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-
-Total Variation 2D Tomography Reconstruction using PDHG algorithm:
-
-
-Problem:     min_u  \alpha * ||\nabla u||_{2,1} + \frac{1}{2}||Au - g||^{2}
-             min_u, u>0  \alpha * ||\nabla u||_{2,1} + \int A u  - g log (Au + \eta)
-
-             \nabla: Gradient operator              
-             A: System Matrix
-             g: Noisy sinogram 
-             \eta: Background noise
-             
-             \alpha: Regularization parameter
- 
-"""
-
-from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction, KullbackLeibler, IndicatorBox
-                      
-from ccpi.astra.ops import AstraProjectorSimple
-from ccpi.framework import TestData
-from PIL import Image
-import os, sys
-#if int(numpy.version.version.split('.')[1]) > 12:
-from skimage.util import random_noise
-#else:
-#    from demoutil import random_noise
-
-#import scipy.io
-
-# user supplied input
-if len(sys.argv) > 1:
-    which_noise = int(sys.argv[1])
-else:
-    which_noise = 1
-
-# Load 256 shepp-logan
-data256 = scipy.io.loadmat('phantom.mat')['phantom256']
-data = ImageData(numpy.array(Image.fromarray(data256).resize((256,256))))
-N, M = data.shape
-ig = ImageGeometry(voxel_num_x=N, voxel_num_y=M)
-
-# Add it to testdata or use tomophantom
-#loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(50, 50))
-#ig = data.geometry
-
-# Create acquisition data and geometry
-detectors = N
-angles = np.linspace(0, np.pi, 180)
-ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-
-# Select device
-device = '0'
-#device = input('Available device: GPU==1 / CPU==0 ')
-if device=='1':
-    dev = 'gpu'
-else:
-    dev = 'cpu'
-    
-Aop = AstraProjectorSimple(ig, ag, dev)
-sin = Aop.direct(data)
-
-# Create noisy data. Apply Gaussian noise
-noises = ['gaussian', 'poisson']
-noise = noises[which_noise]
-
-if noise == 'poisson':
-    scale = 5
-    eta = 0
-    noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
-elif noise == 'gaussian':
-    n1 = np.random.normal(0, 1, size = ag.shape)
-    noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
-else:
-    raise ValueError('Unsupported Noise ', noise)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-plt.subplot(1,2,2)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,2,1)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Create operators
-op1 = Gradient(ig)
-op2 = Aop
-
-# Create BlockOperator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-# Compute operator Norm
-normK = operator.norm()
-
-# Create functions
-if noise == 'poisson':
-    alpha = 3
-    f2 = KullbackLeibler(noisy_data)  
-    g =  IndicatorBox(lower=0) 
-    sigma = 1
-    tau = 1/(sigma*normK**2)    
-    
-elif noise == 'gaussian':   
-    alpha = 20
-    f2 = 0.5 * L2NormSquared(b=noisy_data)                                         
-    g = ZeroFunction()
-    sigma = 10
-    tau = 1/(sigma*normK**2)     
-    
-f1 = alpha * MixedL21Norm() 
-f = BlockFunction(f1, f2)    
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000)
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
deleted file mode 100644
index 78d4980..0000000
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
+++ /dev/null
@@ -1,258 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-
-``Tikhonov`` Regularization Denoising using PDHG algorithm:
-
-
-Problem:     min_{u},   \alpha * ||\nabla u||_{2,2} + Fidelity(u, g)
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: Noisy Data 
-                          
-             Fidelity =  1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian
-                         2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper
-                         3) Kullback Leibler (\int u - g * log(u) + Id_{u>0})  if Noise is Poisson
-                                                       
-             Method = 0 ( PDHG - split ) :  K = [ \nabla,
-                                                 Identity]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  K = \nabla    
-             
-             
-             Default: Tikhonov denoising
-             noise = Gaussian
-             Fidelity = L2NormSquarred 
-             method = 0
-             
-"""      
-
-from ccpi.framework import ImageData, TestData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared,\
-      BlockFunction, KullbackLeibler, L1Norm
-
-import sys, os
-if int(numpy.version.version.split('.')[1]) > 12:
-    from skimage.util import random_noise
-else:
-    from demoutil import random_noise
-
-
-# user supplied input
-if len(sys.argv) > 1:
-    which_noise = int(sys.argv[1])
-else:
-    which_noise = 2
-print ("Applying {} noise")
-
-if len(sys.argv) > 2:
-    method = sys.argv[2]
-else:
-    method = '0'
-print ("method ", method)
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-data = loader.load(TestData.SHAPES)
-ig = data.geometry
-ag = ig
-
-# Create noisy data. 
-noises = ['gaussian', 'poisson', 's&p']
-noise = noises[which_noise]
-if noise == 's&p':
-    n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2)
-elif noise == 'poisson':
-    scale = 5
-    n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale
-elif noise == 'gaussian':
-    n1 = random_noise(data.as_array(), mode = noise, seed = 10)
-else:
-    raise ValueError('Unsupported Noise ', noise)
-noisy_data = ImageData(n1)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,5))
-plt.subplot(1,2,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,2,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Regularisation Parameter depending on the noise distribution
-if noise == 's&p':
-    alpha = 20
-elif noise == 'poisson':
-    alpha = 10
-elif noise == 'gaussian':
-    alpha = 5
-    
-# fidelity
-if noise == 's&p':
-    f2 = L1Norm(b=noisy_data)
-elif noise == 'poisson':
-    f2 = KullbackLeibler(noisy_data)
-elif noise == 'gaussian':
-    f2 = 0.5 * L2NormSquared(b=noisy_data)    
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions      
-    f1 = alpha * L2NormSquared()
-    f = BlockFunction(f1, f2)                                       
-    g = ZeroFunction()
-    
-else:
-    
-    operator = Gradient(ig)
-    g = f2
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 100
-pdhg.run(2000)
-
-
-plt.figure(figsize=(20,5))
-plt.subplot(1,4,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,4,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(1,4,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.subplot(1,4,4)
-plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'Tikhonov reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-
-##%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u = Variable(ig.shape)
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    # Define Total Variation as a regulariser
-    
-    regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS      
-    
-    # fidelity
-    if noise == 's&p':
-        fidelity = pnorm( u - noisy_data.as_array(),1)
-    elif noise == 'poisson':
-        fidelity = sum(kl_div(noisy_data.as_array(), u)) 
-        solver = SCS
-    elif noise == 'gaussian':
-        fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
-        
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
-#
-#
-#
-#
-#
diff --git a/Wrappers/Python/demos/PDHG_examples/IMATDemo.py b/Wrappers/Python/demos/PDHG_examples/IMATDemo.py
deleted file mode 100644
index 2051860..0000000
--- a/Wrappers/Python/demos/PDHG_examples/IMATDemo.py
+++ /dev/null
@@ -1,339 +0,0 @@
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Mar 25 12:50:27 2019
-
-@author: vaggelis
-"""
-
-from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData
-from astropy.io import fits
-import numpy
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import CGLS, PDHG
-from ccpi.optimisation.functions import MixedL21Norm, L2NormSquared, BlockFunction, ZeroFunction, KullbackLeibler, IndicatorBox
-from ccpi.optimisation.operators import Gradient, BlockOperator
-
-from ccpi.astra.operators import  AstraProjectorMC, AstraProjectorSimple
-
-import pickle
-
-
-# load file
-
-#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_282.fits' 
-#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_564.fits'
-#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_141.fits'
-filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_80_channels.fits'
-
-sino_handler = fits.open(filename_sino)
-sino = numpy.array(sino_handler[0].data, dtype=float)
-
-# change axis order: channels, angles, detectors
-sino_new = numpy.rollaxis(sino, 2)
-sino_handler.close()
-
-
-sino_shape = sino_new.shape
-
-num_channels = sino_shape[0] # channelss
-num_pixels_h = sino_shape[2] # detectors
-num_pixels_v = sino_shape[2] # detectors
-num_angles = sino_shape[1] # angles
-
-
-ig = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v, channels = num_channels)
-
-with open("/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/IMAT_data/golden_angles_new.txt") as f:
-    angles_string = [line.rstrip() for line in f]
-    angles = numpy.array(angles_string).astype(float)
-                                                                                 
-
-ag = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h, channels = num_channels)
-op_MC = AstraProjectorMC(ig, ag, 'gpu')
-
-sino_aqdata = AcquisitionData(sino_new, ag)
-result_bp = op_MC.adjoint(sino_aqdata)
-
-#%%
-
-channel = [40, 60]
-for j in range(2):        
-    z4 = sino_aqdata.as_array()[channel[j]] 
-    plt.figure(figsize=(10,6))
-    plt.imshow(z4, cmap='viridis')
-    plt.axis('off')
-    plt.savefig('Sino_141/Sinogram_ch_{}_.png'.format(channel[j]), bbox_inches='tight', transparent=True)
-    plt.show() 
-
-#%%
-
-def callback(iteration, objective, x):
-    plt.imshow(x.as_array()[40])
-    plt.colorbar()
-    plt.show()
-
-#%%
-# CGLS 
-
-x_init = ig.allocate()      
-cgls1 = CGLS(x_init=x_init, operator=op_MC, data=sino_aqdata)
-cgls1.max_iteration = 100
-cgls1.update_objective_interval = 2
-cgls1.run(20,verbose=True, callback=callback)
-
-plt.imshow(cgls1.get_output().subset(channel=20).array)
-plt.title('CGLS')
-plt.colorbar()
-plt.show()
-
-#%%
-with open('Sino_141/CGLS/CGLS_{}_iter.pkl'.format(20), 'wb') as f:
-    z = cgls1.get_output()
-    pickle.dump(z, f)  
-    
-#%% 
-#% Tikhonov Space
-    
-x_init = ig.allocate() 
-alpha = [1,3,5,10,20,50]
-
-for a in alpha:
-    
-    Grad = Gradient(ig, correlation = Gradient.CORRELATION_SPACE)
-    operator = BlockOperator(op_MC, a * Grad, shape=(2,1))
-    blockData = BlockDataContainer(sino_aqdata, \
-                                   Grad.range_geometry().allocate())
-    cgls2 = CGLS()
-    cgls2.max_iteration = 500
-    cgls2.set_up(x_init, operator, blockData)
-    cgls2.update_objective_interval = 50
-    cgls2.run(100,verbose=True)
-        
-    with open('Sino_141/CGLS_Space/CGLS_a_{}.pkl'.format(a), 'wb') as f:
-        z = cgls2.get_output()
-        pickle.dump(z, f)
-   
-#% Tikhonov SpaceChannels
-
-for a1 in alpha:
-    
-    Grad1 = Gradient(ig, correlation = Gradient.CORRELATION_SPACECHANNEL)
-    operator1 = BlockOperator(op_MC, a1 * Grad1, shape=(2,1))
-    blockData1 = BlockDataContainer(sino_aqdata, \
-                                   Grad1.range_geometry().allocate())
-    cgls3 = CGLS()
-    cgls3.max_iteration = 500
-    cgls3.set_up(x_init, operator1, blockData1)
-    cgls3.update_objective_interval = 10
-    cgls3.run(100, verbose=True)
-    
-    with open('Sino_141/CGLS_SpaceChannels/CGLS_a_{}.pkl'.format(a1), 'wb') as f1:
-        z1 = cgls3.get_output()
-        pickle.dump(z1, f1)
-    
-  
-    
-#%%
-#
-ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v)
-ag_tmp = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h)
-op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu')
-normK1 = op_tmp.norm()
-
-alpha_TV = [2, 5, 10] # for powder
-
-# Create operators
-op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL)
-op2 = op_MC
-
-# Create BlockOperator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-
-for alpha in alpha_TV:
-# Create functions      
-    f1 = alpha * MixedL21Norm()
-    
-    f2 = KullbackLeibler(sino_aqdata)    
-    f = BlockFunction(f1, f2)  
-    g = IndicatorBox(lower=0)
-    
-    # Compute operator Norm
-    normK = numpy.sqrt(8 + normK1**2) 
-    
-    # Primal & dual stepsizes
-    sigma = 1
-    tau = 1/(sigma*normK**2)
-    
-    # Setup and run the PDHG algorithm
-    pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-    pdhg.max_iteration = 5000
-    pdhg.update_objective_interval = 500
-#    pdhg.run(2000, verbose=True, callback=callback)   
-    pdhg.run(5000, verbose=True, callback=callback)  
-#    
-    with open('Sino_141/TV_SpaceChannels/TV_a = {}.pkl'.format(alpha), 'wb') as f3:
-        z3 = pdhg.get_output()
-        pickle.dump(z3, f3)    
-#    
-#
-#
-#
-##%%
-#        
-#ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v)
-#ag_tmp = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h)
-#op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu')
-#normK1 = op_tmp.norm()
-#
-#alpha_TV = 10 # for powder
-#
-## Create operators
-#op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL)
-#op2 = op_MC
-#
-## Create BlockOperator
-#operator = BlockOperator(op1, op2, shape=(2,1) ) 
-#
-#
-## Create functions      
-#f1 = alpha_TV * MixedL21Norm()
-#f2 = 0.5 * L2NormSquared(b=sino_aqdata)    
-#f = BlockFunction(f1, f2)  
-#g = ZeroFunction()
-#
-## Compute operator Norm
-##normK = 8.70320267279591 # For powder Run one time no need to compute again takes time
-#normK = numpy.sqrt(8 + normK1**2) # for carbon
-#
-## Primal & dual stepsizes
-#sigma = 0.1
-#tau = 1/(sigma*normK**2)
-#
-#def callback(iteration, objective, x):
-#    plt.imshow(x.as_array()[100])
-#    plt.colorbar()
-#    plt.show()
-#
-## Setup and run the PDHG algorithm
-#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-#pdhg.max_iteration = 2000
-#pdhg.update_objective_interval = 100
-#pdhg.run(2000, verbose=True)       
-#        
-#        
-#        
-#        
-#        
-#        
-        
-        
-        
-        
-        
-        
-        
-        
-        
-        
-#%%
-    
-#with open('/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/CGLS_Tikhonov/CGLS_Space/CGLS_Space_a = 50.pkl', 'wb') as f:
-#    z = cgls2.get_output()
-#    pickle.dump(z, f)
-#      
-        
-        #%%
-with open('Sino_141/CGLS_Space/CGLS_Space_a_20.pkl', 'rb') as f1:
-    x = pickle.load(f1)
-    
-with open('Sino_141/CGLS_SpaceChannels/CGLS_SpaceChannels_a_20.pkl', 'rb') as f1:
-    x1 = pickle.load(f1)    
-    
-   
-    
-#    
-plt.imshow(x.as_array()[40]*mask)
-plt.colorbar()
-plt.show()
-
-plt.imshow(x1.as_array()[40]*mask)
-plt.colorbar()
-plt.show()
-
-plt.plot(x.as_array()[40,100,:])
-plt.plot(x1.as_array()[40,100,:])
-plt.show()
-
-#%%
-
-# Show results
-
-def circ_mask(h, w, center=None, radius=None):
-
-    if center is None: # use the middle of the image
-        center = [int(w/2), int(h/2)]
-    if radius is None: # use the smallest distance between the center and image walls
-        radius = min(center[0], center[1], w-center[0], h-center[1])
-
-    Y, X = numpy.ogrid[:h, :w]
-    dist_from_center = numpy.sqrt((X - center[0])**2 + (Y-center[1])**2)
-
-    mask = dist_from_center <= radius
-    return mask
-
-mask = circ_mask(141, 141, center=None, radius = 55)
-plt.imshow(numpy.multiply(x.as_array()[40],mask))
-plt.show()
-#%%
-#channel = [100, 200, 300]
-#
-#for i in range(3):     
-#    tmp = cgls1.get_output().as_array()[channel[i]] 
-#    
-#    z =  tmp * mask
-#    plt.figure(figsize=(10,6))
-#    plt.imshow(z, vmin=0, cmap='viridis')
-#    plt.axis('off')
-##    plt.clim(0, 0.02)
-##    plt.colorbar()
-##    del z
-#    plt.savefig('CGLS_282/CGLS_Chan_{}.png'.format(channel[i]), bbox_inches='tight', transparent=True)
-#    plt.show()
-#    
-#    
-##%% Line Profiles
-#    
-#n1, n2, n3 =  cgs.get_output().as_array().shape
-#mask = circ_mask(564, 564, center=None, radius = 220)
-#material = ['Cu', 'Fe', 'Ni']
-#ycoords = [200, 300, 380]
-#    
-#for i in range(3):
-#    z = cgs.get_output().as_array()[channel[i]]  * mask
-#    
-#    for k1 in range(len(ycoords)):
-#        plt.plot(numpy.arange(0,n2), z[ycoords[k1],:])
-#        plt.title('Channel {}: {}'.format(channel[i], material[k1]))
-#        plt.savefig('CGLS/line_profile_chan_{}_material_{}.png'.\
-#                        format(channel[i], material[k1]), bbox_inches='tight')
-#        plt.show()
-#    
-#    
-#
-#
-#
-##%%
-#
-#%%
-        
-
-
-#%%
-
-#plt.imshow(pdhg.get_output().subset(channel=100).as_array())
-#plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_2D_time_denoising.py b/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_2D_time_denoising.py
deleted file mode 100644
index 045458a..0000000
--- a/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_2D_time_denoising.py
+++ /dev/null
@@ -1,169 +0,0 @@
-# -*- 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-2019 Evangelos Papoutsellis and 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.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorMC
-
-import os
-import tomophantom
-from tomophantom import TomoP2D
-
-# Create phantom for TV 2D dynamic tomography 
-
-model = 102  # note that the selected model is temporal (2D + time)
-N = 50 # set dimension of the phantom
-# one can specify an exact path to the parameters file
-# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
-path = os.path.dirname(tomophantom.__file__)
-path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
-#This will generate a N_size x N_size x Time frames phantom (2D + time)
-phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D)
-
-plt.close('all')
-plt.figure(1)
-plt.rcParams.update({'font.size': 21})
-plt.title('{}''{}'.format('2D+t phantom using model no.',model))
-for sl in range(0,np.shape(phantom_2Dt)[0]):
-    im = phantom_2Dt[sl,:,:]
-    plt.imshow(im, vmin=0, vmax=1)
-    plt.pause(.1)
-    plt.draw
-
-    
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
-data = ImageData(phantom_2Dt, geometry=ig)
-
-detectors = N
-angles = np.linspace(0,np.pi,N)
-
-ag = AcquisitionGeometry('parallel','2D', angles, detectors, channels = np.shape(phantom_2Dt)[0])
-Aop = AstraProjectorMC(ig, ag, 'gpu')
-sin = Aop.direct(data)
-
-scale = 2
-n1 = scale * np.random.poisson(sin.as_array()/scale)
-noisy_data = AcquisitionData(n1, ag)
-
-tindex = [3, 6, 10]
-
-fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10))
-plt.subplot(1,3,1)
-plt.imshow(noisy_data.as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-plt.subplot(1,3,2)
-plt.imshow(noisy_data.as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-plt.subplot(1,3,3)
-plt.imshow(noisy_data.as_array()[tindex[2],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-plt.show()   
-
-#%%
-# Regularisation Parameter
-alpha = 5
-
-# Create operators
-#op1 = Gradient(ig)
-op1 = Gradient(ig, correlation='SpaceChannels')
-op2 = Aop
-
-# Create BlockOperator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-# Create functions
-      
-f1 = alpha * MixedL21Norm()
-f2 = KullbackLeibler(noisy_data)    
-f = BlockFunction(f1, f2)  
-                                      
-g = ZeroFunction()
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000)
-
-
-#%%
-fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
-
-plt.subplot(2,3,1)
-plt.imshow(phantom_2Dt[tindex[0],:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-
-plt.subplot(2,3,2)
-plt.imshow(phantom_2Dt[tindex[1],:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-
-plt.subplot(2,3,3)
-plt.imshow(phantom_2Dt[tindex[2],:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-
-plt.subplot(2,3,4)
-plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.subplot(2,3,5)
-plt.imshow(pdhg.get_output().as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.subplot(2,3,6)
-plt.imshow(pdhg.get_output().as_array()[tindex[2],:,:])
-plt.axis('off')
-im = plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:])
-
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
-cbar = fig.colorbar(im, cax=cb_ax)
-
-
-plt.show()
-
diff --git a/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_TV_Tomo2D_time.py b/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_TV_Tomo2D_time.py
deleted file mode 100644
index 045458a..0000000
--- a/Wrappers/Python/demos/PDHG_examples/MultiChannel/PDHG_TV_Tomo2D_time.py
+++ /dev/null
@@ -1,169 +0,0 @@
-# -*- 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-2019 Evangelos Papoutsellis and 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.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorMC
-
-import os
-import tomophantom
-from tomophantom import TomoP2D
-
-# Create phantom for TV 2D dynamic tomography 
-
-model = 102  # note that the selected model is temporal (2D + time)
-N = 50 # set dimension of the phantom
-# one can specify an exact path to the parameters file
-# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
-path = os.path.dirname(tomophantom.__file__)
-path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
-#This will generate a N_size x N_size x Time frames phantom (2D + time)
-phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D)
-
-plt.close('all')
-plt.figure(1)
-plt.rcParams.update({'font.size': 21})
-plt.title('{}''{}'.format('2D+t phantom using model no.',model))
-for sl in range(0,np.shape(phantom_2Dt)[0]):
-    im = phantom_2Dt[sl,:,:]
-    plt.imshow(im, vmin=0, vmax=1)
-    plt.pause(.1)
-    plt.draw
-
-    
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
-data = ImageData(phantom_2Dt, geometry=ig)
-
-detectors = N
-angles = np.linspace(0,np.pi,N)
-
-ag = AcquisitionGeometry('parallel','2D', angles, detectors, channels = np.shape(phantom_2Dt)[0])
-Aop = AstraProjectorMC(ig, ag, 'gpu')
-sin = Aop.direct(data)
-
-scale = 2
-n1 = scale * np.random.poisson(sin.as_array()/scale)
-noisy_data = AcquisitionData(n1, ag)
-
-tindex = [3, 6, 10]
-
-fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10))
-plt.subplot(1,3,1)
-plt.imshow(noisy_data.as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-plt.subplot(1,3,2)
-plt.imshow(noisy_data.as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-plt.subplot(1,3,3)
-plt.imshow(noisy_data.as_array()[tindex[2],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-plt.show()   
-
-#%%
-# Regularisation Parameter
-alpha = 5
-
-# Create operators
-#op1 = Gradient(ig)
-op1 = Gradient(ig, correlation='SpaceChannels')
-op2 = Aop
-
-# Create BlockOperator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-# Create functions
-      
-f1 = alpha * MixedL21Norm()
-f2 = KullbackLeibler(noisy_data)    
-f = BlockFunction(f1, f2)  
-                                      
-g = ZeroFunction()
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000)
-
-
-#%%
-fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
-
-plt.subplot(2,3,1)
-plt.imshow(phantom_2Dt[tindex[0],:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-
-plt.subplot(2,3,2)
-plt.imshow(phantom_2Dt[tindex[1],:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-
-plt.subplot(2,3,3)
-plt.imshow(phantom_2Dt[tindex[2],:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-
-plt.subplot(2,3,4)
-plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.subplot(2,3,5)
-plt.imshow(pdhg.get_output().as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.subplot(2,3,6)
-plt.imshow(pdhg.get_output().as_array()[tindex[2],:,:])
-plt.axis('off')
-im = plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:])
-
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
-cbar = fig.colorbar(im, cax=cb_ax)
-
-
-plt.show()
-
diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py
deleted file mode 100644
index ddf5ace..0000000
--- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py
+++ /dev/null
@@ -1,115 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-Total Variation Denoising using PDHG algorithm:
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{2,1} + ||x-g||_{1}
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: Noisy Data with Salt & Pepper Noise
-             
-             
-             Method = 0 ( PDHG - split ) :  K = [ \nabla,
-                                                 Identity]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  K = \nabla    
-             
-             
-"""
-
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff
-from ccpi.optimisation.functions import MixedL21Norm, MixedL11Norm, L2NormSquared, BlockFunction, L1Norm                      
-from ccpi.framework import TestData, ImageGeometry
-import os, sys
-if int(numpy.version.version.split('.')[1]) > 12:
-    from skimage.util import random_noise
-else:
-    from demoutil import random_noise
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-data = loader.load(TestData.PEPPERS, size=(256,256))
-ig = data.geometry
-ag = ig
-
-# Create noisy data. 
-n1 = random_noise(data.as_array(), mode = 'gaussian', var = 0.15, seed = 50)
-noisy_data = ig.allocate()
-noisy_data.fill(n1)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,5))
-plt.subplot(1,2,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,2,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Regularisation Parameter
-operator = Gradient(ig, correlation=Gradient.CORRELATION_SPACE)
-f1 =  5 * MixedL21Norm()
-g = 0.5 * L2NormSquared(b=noisy_data)
-            
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-# Setup and run the PDHG algorithm
-pdhg1 = PDHG(f=f1,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg1.max_iteration = 2000
-pdhg1.update_objective_interval = 200
-pdhg1.run(1000)
-
-
-# Show results
-plt.figure(figsize=(10,10))
-plt.subplot(2,2,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(2,2,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(2,2,3)
-plt.imshow(pdhg1.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.subplot(2,2,4)
-plt.imshow(pdhg2.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-
diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_2D_time.py
deleted file mode 100644
index 14608db..0000000
--- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_2D_time.py
+++ /dev/null
@@ -1,192 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorMC
-
-import os
-import tomophantom
-from tomophantom import TomoP2D
-
-# Create phantom for TV 2D dynamic tomography 
-
-model = 102  # note that the selected model is temporal (2D + time)
-N = 128 # set dimension of the phantom
-# one can specify an exact path to the parameters file
-# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
-path = os.path.dirname(tomophantom.__file__)
-path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
-#This will generate a N_size x N_size x Time frames phantom (2D + time)
-phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D)
-
-plt.close('all')
-plt.figure(1)
-plt.rcParams.update({'font.size': 21})
-plt.title('{}''{}'.format('2D+t phantom using model no.',model))
-for sl in range(0,np.shape(phantom_2Dt)[0]):
-    im = phantom_2Dt[sl,:,:]
-    plt.imshow(im, vmin=0, vmax=1)
-#    plt.pause(.1)
-#    plt.draw
-
-    
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
-data = ImageData(phantom_2Dt, geometry=ig)
-ag = ig
-                      
-# Create Noisy data. Add Gaussian noise
-np.random.seed(10)
-noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.25, size=ig.shape) )
-
-tindex = [3, 6, 10]
-
-fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10))
-plt.subplot(1,3,1)
-plt.imshow(noisy_data.as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-plt.subplot(1,3,2)
-plt.imshow(noisy_data.as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-plt.subplot(1,3,3)
-plt.imshow(noisy_data.as_array()[tindex[2],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-plt.show()   
-
-#%%
-# Regularisation Parameter
-alpha = 0.3
-
-# Create operators
-#op1 = Gradient(ig)
-op1 = Gradient(ig, correlation='Space')
-op2 = Gradient(ig, correlation='SpaceChannels')
-
-op3 = Identity(ig, ag)
-
-# Create BlockOperator
-operator1 = BlockOperator(op1, op3, shape=(2,1) ) 
-operator2 = BlockOperator(op2, op3, shape=(2,1) )
-
-# Create functions
-      
-f1 = alpha * MixedL21Norm()
-f2 = 0.5 * L2NormSquared(b = noisy_data)    
-f = BlockFunction(f1, f2)  
-                                                                        
-g = ZeroFunction()
-    
-# Compute operator Norm
-normK1 = operator1.norm()
-normK2 = operator2.norm()
-
-#%%
-# Primal & dual stepsizes
-sigma1 = 1
-tau1 = 1/(sigma1*normK1**2)
-
-sigma2 = 1
-tau2 = 1/(sigma2*normK2**2)
-
-# Setup and run the PDHG algorithm
-pdhg1 = PDHG(f=f,g=g,operator=operator1, tau=tau1, sigma=sigma1)
-pdhg1.max_iteration = 2000
-pdhg1.update_objective_interval = 200
-pdhg1.run(2000)
-
-# Setup and run the PDHG algorithm
-pdhg2 = PDHG(f=f,g=g,operator=operator2, tau=tau2, sigma=sigma2)
-pdhg2.max_iteration = 2000
-pdhg2.update_objective_interval = 200
-pdhg2.run(2000)
-
-
-#%%
-
-tindex = [3, 6, 10]
-fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
-
-plt.subplot(3,3,1)
-plt.imshow(phantom_2Dt[tindex[0],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[0]))
-
-plt.subplot(3,3,2)
-plt.imshow(phantom_2Dt[tindex[1],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[1]))
-
-plt.subplot(3,3,3)
-plt.imshow(phantom_2Dt[tindex[2],:,:])
-plt.axis('off')
-plt.title('Time {}'.format(tindex[2]))
-
-plt.subplot(3,3,4)
-plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.subplot(3,3,5)
-plt.imshow(pdhg1.get_output().as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.subplot(3,3,6)
-plt.imshow(pdhg1.get_output().as_array()[tindex[2],:,:])
-plt.axis('off')
-
-
-plt.subplot(3,3,7)
-plt.imshow(pdhg2.get_output().as_array()[tindex[0],:,:])
-plt.axis('off')
-plt.subplot(3,3,8)
-plt.imshow(pdhg2.get_output().as_array()[tindex[1],:,:])
-plt.axis('off')
-plt.subplot(3,3,9)
-plt.imshow(pdhg2.get_output().as_array()[tindex[2],:,:])
-plt.axis('off')
-
-#%%
-im = plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:])
-
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
-cbar = fig.colorbar(im, cax=cb_ax)
-
-
-plt.show()
-
diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py
deleted file mode 100644
index 03dc2ef..0000000
--- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian_3D.py
+++ /dev/null
@@ -1,181 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-""" 
-
-Total Variation (3D) Denoising using PDHG algorithm:
-
-
-Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: Noisy Data with Gaussian Noise
-                          
-             Method = 0 ( PDHG - split ) :  K = [ \nabla,
-                                                 Identity]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  K = \nabla  
-                                                                
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-                         
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Gaussian denoising
-import timeit
-import os
-from tomophantom import TomoP3D
-import tomophantom
-
-print ("Building 3D phantom using TomoPhantom software")
-tic=timeit.default_timer()
-model = 13 # select a model number from the library
-N = 64 # Define phantom dimensions using a scalar value (cubic phantom)
-path = os.path.dirname(tomophantom.__file__)
-path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
-
-#This will generate a N x N x N phantom (3D)
-phantom_tm = TomoP3D.Model(model, N, path_library3D)
-
-#%%
-
-# Create noisy data. Add Gaussian noise
-ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N)
-ag = ig
-n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10)
-noisy_data = ImageData(n1)
-
-sliceSel = int(0.5*N)
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1)
-plt.title('Axial View')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1)
-plt.title('Coronal View')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1)
-plt.title('Sagittal View')
-plt.colorbar()
-plt.show()   
-
-#%%
-
-# Regularisation Parameter
-alpha = 0.05
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * MixedL21Norm()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000, verbose = True)
-
-
-#%%
-fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
-fig.suptitle('TV Reconstruction',fontsize=20)
-
-
-plt.subplot(2,3,1)
-plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Axial View')
-
-plt.subplot(2,3,2)
-plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Coronal View')
-
-plt.subplot(2,3,3)
-plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1)
-plt.axis('off')
-plt.title('Sagittal View')
-
-
-plt.subplot(2,3,4)
-plt.imshow(pdhg.get_output().as_array()[sliceSel,:,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.subplot(2,3,5)
-plt.imshow(pdhg.get_output().as_array()[:,sliceSel,:],vmin=0, vmax=1)
-plt.axis('off')
-plt.subplot(2,3,6)
-plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1)
-plt.axis('off')
-im = plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1)
-
-
-fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-
-cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
-cbar = fig.colorbar(im, cax=cb_ax)
-
-
-plt.show()
-
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_Tikhonov_Tomo2D.py
deleted file mode 100644
index 02cd053..0000000
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_Tikhonov_Tomo2D.py
+++ /dev/null
@@ -1,156 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#         http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# 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.
-#
-#=========================================================================
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{2}^{2} + int A x -g log(Ax + \eta)
-
-             \nabla: Gradient operator 
-             
-             A: Projection Matrix
-             g: Noisy sinogram corrupted with Poisson Noise
-             
-             \eta: Background Noise
-             \alpha: Regularization parameter
-             
-                
-                       
-"""
-
-
-from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorSimple
-from ccpi.framework import TestData
-import os, sys
-
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-
-# Load Data                      
-N = 100
-M = 100
-data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
-
-ig = data.geometry
-ag = ig
-
-#Create Acquisition Data and apply poisson noise
-
-detectors = N
-angles = np.linspace(0, np.pi, N)
-
-ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-
-device = input('Available device: GPU==1 / CPU==0 ')
-
-if device=='1':
-    dev = 'gpu'
-else:
-    dev = 'cpu'
-
-Aop = AstraProjectorSimple(ig, ag, 'cpu')
-sin = Aop.direct(data)
-
-# Create noisy data. Apply Poisson noise
-scale = 0.5
-eta = 0 
-n1 = scale * np.random.poisson(eta + sin.as_array()/scale)
-
-noisy_data = AcquisitionData(n1, ag)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-plt.subplot(2,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(2,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-
-# Regularisation Parameter
-alpha = 1000
-
-# Create operators
-op1 = Gradient(ig)
-op2 = Aop
-
-# Create BlockOperator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-# Create functions
-      
-f1 = alpha * L2NormSquared()
-f2 = 0.5 * L2NormSquared(b=noisy_data)    
-f = BlockFunction(f1, f2)  
-                     
-g = ZeroFunction()
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 500
-pdhg.run(2000)
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('Tikhonov Reconstruction')
-plt.colorbar()
-plt.show()
-## 
-plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,M), pdhg.get_output().as_array()[int(N/2),:], label = 'Tikhonov reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py
deleted file mode 100644
index 854f645..0000000
--- a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py
+++ /dev/null
@@ -1,238 +0,0 @@
-# -*- coding: utf-8 -*-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 22 14:53:03 2019
-
-@author: evangelos
-"""
-
-from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, \
-   AcquisitionGeometry, AcquisitionData
-
-import numpy as np                           
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, PDHG_old
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction, ScaledFunction
-
-from ccpi.plugins.operators import CCPiProjectorSimple
-from timeit import default_timer as timer
-from ccpi.reconstruction.parallelbeam import alg as pbalg
-import os
-
-try:
-    import tomophantom
-    from tomophantom import TomoP3D
-    no_tomophantom = False
-except ImportError as ie:
-    no_tomophantom = True
-
-#%%
-
-#%%###############################################################################
-# Create phantom for TV tomography
-
-#import os
-#import tomophantom
-#from tomophantom import TomoP2D
-#from tomophantom.supp.qualitymetrics import QualityTools
-
-#model = 1 # select a model number from the library
-#N = 150 # set dimension of the phantom
-## one can specify an exact path to the parameters file
-## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
-#path = os.path.dirname(tomophantom.__file__)
-#path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
-##This will generate a N_size x N_size phantom (2D)
-#phantom_2D = TomoP2D.Model(model, N, path_library2D)
-#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-#data = ImageData(phantom_2D, geometry=ig)
-
-N = 75
-#x = np.zeros((N,N))
-
-vert = 4
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert)
-
-angles_num = 100
-det_w = 1.0
-det_num = N
-
-angles = np.linspace(-90.,90.,N, dtype=np.float32) 
-# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, 
-#         horz detector pixel size, vert detector pixel count, 
-#         vert detector pixel size.
-ag = AcquisitionGeometry('parallel',
-                         '3D',
-                         angles,
-                         N, 
-                         det_w,
-                         vert,
-                         det_w)
-
-#no_tomophantom = True
-if no_tomophantom:
-    data = ig.allocate()
-    Phantom = data
-    # Populate image data by looping over and filling slices
-    i = 0
-    while i < vert:
-        if vert > 1:
-            x = Phantom.subset(vertical=i).array
-        else:
-            x = Phantom.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)] = 0.98
-        if vert > 1 :
-            Phantom.fill(x, vertical=i)
-        i += 1
-    
-    Aop = CCPiProjectorSimple(ig, ag, 'cpu')
-    sin = Aop.direct(data)
-else:
-        
-    model = 13 # select a model number from the library
-    N_size = N # Define phantom dimensions using a scalar value (cubic phantom)
-    path = os.path.dirname(tomophantom.__file__)
-    path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
-    #This will generate a N_size x N_size x N_size phantom (3D)
-    phantom_tm = TomoP3D.Model(model, N_size, path_library3D)
-    
-    #%%
-    Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal)
-    Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
-    #angles_num = int(0.5*np.pi*N_size); # angles number
-    #angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees
-    
-    print ("Building 3D analytical projection data with TomoPhantom")
-    projData3D_analyt = TomoP3D.ModelSino(model, 
-                                                      N_size, 
-                                                      Horiz_det, 
-                                                      Vert_det, 
-                                                      angles, 
-                                                      path_library3D)
-    
-    # tomophantom outputs in [vert,angles,horiz]
-    # we want [angle,vert,horiz]
-    data = np.transpose(projData3D_analyt, [1,0,2])
-    ag.pixel_num_h = Horiz_det
-    ag.pixel_num_v = Vert_det
-    sin = ag.allocate()
-    sin.fill(data)
-    ig.voxel_num_y = Vert_det
-    
-    Aop = CCPiProjectorSimple(ig, ag, 'cpu')
-    
-
-plt.imshow(sin.subset(vertical=0).as_array())
-plt.title('Sinogram')
-plt.colorbar()
-plt.show()
-
-
-#%%
-# Add Gaussian noise to the sinogram data
-np.random.seed(10)
-n1 = np.random.random(sin.shape)
-
-noisy_data = sin + ImageData(5*n1)
-
-plt.imshow(noisy_data.subset(vertical=0).as_array())
-plt.title('Noisy Sinogram')
-plt.colorbar()
-plt.show()
-
-
-#%% Works only with Composite Operator Structure of PDHG
-
-#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-
-# Create operators
-op1 = Gradient(ig)
-op2 = Aop
-
-# Form Composite Operator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-alpha = 50
-f = BlockFunction( alpha * MixedL21Norm(), \
-                   0.5 * L2NormSquared(b = noisy_data) )
-g = ZeroFunction()
-
-normK = Aop.norm()
-
-# Compute operator Norm
-normK = operator.norm()
-
-## Primal & dual stepsizes
-diag_precon = False 
-
-if diag_precon:
-    
-    def tau_sigma_precond(operator):
-        
-        tau = 1/operator.sum_abs_row()
-        sigma = 1/ operator.sum_abs_col()
-    
-        return tau, sigma
-
-    tau, sigma = tau_sigma_precond(operator)
-    
-else:    
-    sigma = 1
-    tau = 1/(sigma*normK**2)
-
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-niter = 50
-opt = {'niter':niter}
-opt1 = {'niter':niter, 'memopt': True}
-
-
-
-pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, max_iteration=niter)
-#pdhg1.max_iteration = 2000
-pdhg1.update_objective_interval = 100
-
-t1_old = timer()
-resold, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-t2_old = timer()
-
-pdhg1.run(niter)
-print (sum(pdhg1.timing))
-res = pdhg1.get_output().subset(vertical=0)
-
-#%%
-plt.figure()
-plt.subplot(1,4,1)
-plt.imshow(res.as_array())
-plt.title('Algorithm')
-plt.colorbar()
-plt.subplot(1,4,2)
-plt.imshow(resold.subset(vertical=0).as_array())
-plt.title('function')
-plt.colorbar()
-plt.subplot(1,4,3)
-plt.imshow((res - resold.subset(vertical=0)).abs().as_array())
-plt.title('diff')
-plt.colorbar()
-plt.subplot(1,4,4)
-plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'Algorithm')
-plt.plot(np.linspace(0,N,N), resold.subset(vertical=0).as_array()[int(N/2),:], label = 'function')
-plt.legend()
-plt.show()
-#
-print ("Time: No memopt in {}s, \n Time: Memopt in  {}s ".format(sum(pdhg1.timing), t2_old -t1_old))
-diff = (res - resold.subset(vertical=0)).abs().as_array().max()
-#
-print(" Max of abs difference is {}".format(diff))
-