Merge remote-tracking branch 'origin/demos' into update_fix_test

13 files changed, 15 insertions, 1713 deletions

diff --git a/Wrappers/Python/ccpi/framework/BlockGeometry.py b/Wrappers/Python/ccpi/framework/BlockGeometry.py
index ed44d99..e58035b 100755
--- a/Wrappers/Python/ccpi/framework/BlockGeometry.py
+++ b/Wrappers/Python/ccpi/framework/BlockGeometry.py
@@ -10,6 +10,12 @@ from ccpi.framework import BlockDataContainer
 #from ccpi.optimisation.operators import Operator, LinearOperator

  

 class BlockGeometry(object):

+    

+    RANDOM = 'random'

+    RANDOM_INT = 'random_int'

+    

+    

+    

     '''Class to hold Geometry as column vector'''

     #__array_priority__ = 1

     def __init__(self, *args, **kwargs):

diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 5f04363..cbdc420 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -104,7 +104,7 @@ class BlockOperator(Operator):
         index = row*self.shape[1]+col
         return self.operators[index]
     
-    def calculate_norm(self, **kwargs):
+    def norm(self, **kwargs):
         norm = [op.norm(**kwargs)**2 for op in self.operators]
         return numpy.sqrt(sum(norm))    
     
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
index d848b9f..6937fa0 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
@@ -63,6 +63,7 @@ from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
                           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:
@@ -211,7 +212,7 @@ else:
     plt.show()
 
 
-##%% Check with CVX solution
+#%% Check with CVX solution
 
 from ccpi.optimisation.operators import SparseFiniteDiff
 
@@ -231,8 +232,7 @@ if cvx_not_installable:
     DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
     
     # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = pnorm( u - noisy_data.as_array(),1)
+    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():
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
index f179e95..4f7639e 100644
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
@@ -52,12 +52,12 @@ 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
+#if int(numpy.version.version.split('.')[1]) > 12:
+from skimage.util import random_noise
+#else:
+#    from demoutil import random_noise
 
-import scipy.io
+#import scipy.io
 
 # user supplied input
 if len(sys.argv) > 1:
diff --git a/Wrappers/Python/demos/PDHG_examples/TGV_Denoising/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/TGV_Denoising/PDHG_TGV_Denoising_SaltPepper.py
deleted file mode 100644
index 15c0a05..0000000
--- a/Wrappers/Python/demos/PDHG_examples/TGV_Denoising/PDHG_TGV_Denoising_SaltPepper.py
+++ /dev/null
@@ -1,245 +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_{x} \alpha * ||\nabla x - w||_{2,1} +
-                     \beta *  || E w ||_{2,1} +
-                     \frac{1}{2} * || x - g ||_{2}^{2}
-
-             \alpha: Regularization parameter
-             \beta: Regularization parameter
-             
-             \nabla: Gradient operator 
-              E: Symmetrized Gradient operator
-             
-             g: Noisy Data with Salt & Pepper Noise
-                          
-             Method = 0 ( PDHG - split ) :  K = [ \nabla, - Identity
-                                                  ZeroOperator, E 
-                                                  Identity, ZeroOperator]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  K = [ \nabla, - Identity
-                                                  ZeroOperator, E ] 
-                                                                
-"""
-
-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, Identity, \
-                        Gradient, SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TGV SaltPepper denoising
-
-N = 100
-
-data = np.zeros((N,N))
-
-x1 = np.linspace(0, int(N/2), N)
-x2 = np.linspace(int(N/2), 0., N)
-xv, yv = np.meshgrid(x1, x2)
-
-xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T
-
-data = xv
-data = ImageData(data/data.max())
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
-noisy_data = ImageData(n1)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(15,15))
-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 Parameters
-alpha = 0.8
-beta = numpy.sqrt(2)* alpha
-
-method = '1'
-
-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() 
-    f3 = L1Norm(b=noisy_data)    
-    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(L1Norm(b=noisy_data), 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 = 50
-pdhg.run(2000, verbose = False)
-
-#%%
-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()[0].as_array())
-plt.title('TGV Reconstruction')
-plt.colorbar()
-plt.show()
-## 
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/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:    
-    
-    u = Variable(ig.shape)
-    w1 = Variable((N, N))
-    w2 = Variable((N, N))
-    
-    # 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 = []
-    fidelity = pnorm(u - noisy_data.as_array(),1)    
-    solver = MOSEK
-
-        # 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 = 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,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u.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 (PDHG) {} '.format(pdhg.objective[-1][0]))
-
-
-
-
-
diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian.py
deleted file mode 100644
index 9d00ee1..0000000
--- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Gaussian.py
+++ /dev/null
@@ -1,225 +0,0 @@
-#========================================================================
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
-
-#   Copyright 2017 UKRI-STFC
-#   Copyright 2017 University of Manchester
-
-#   Licensed under the Apache License, Version 2.0 (the "License");
-#   you may not use this file except in compliance with the License.
-#   You may obtain a copy of the License at
-
-#   http://www.apache.org/licenses/LICENSE-2.0
-
-#   Unless required by applicable law or agreed to in writing, software
-#   distributed under the License is distributed on an "AS IS" BASIS,
-#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-#   See the License for the specific language governing permissions and
-#   limitations under the License.
-#=========================================================================
-""" 
-
-Total Variation 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 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, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.framework import TestData
-import os, sys
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-
-# Load Data                      
-N = 200
-M = 300
-
-
-# user can change the size of the input data
-# you can choose between 
-# TestData.PEPPERS 2D + Channel
-# TestData.BOAT 2D 
-# TestData.CAMERA 2D
-# TestData.RESOLUTION_CHART 2D 
-# TestData.SIMPLE_PHANTOM_2D 2D
-data = loader.load(TestData.BOAT, size=(N,M), scale=(0,1))
-
-ig = data.geometry
-ag = ig
-                      
-# Create Noisy data. Add Gaussian noise
-np.random.seed(10)
-noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.1, size=data.shape) )
-
-print ("min {} max {}".format(data.as_array().min(), data.as_array().max()))
-
-# Show Ground Truth and Noisy Data
-plt.figure()
-plt.subplot(1,3,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(1,3,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(1,3,3)
-plt.imshow((data - noisy_data).as_array())
-plt.title('diff')
-plt.colorbar()
-
-plt.show()
-
-# Regularisation Parameter
-alpha = .1
-
-method = '0'
-
-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      
-    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)
-pdhg.max_iteration = 10000
-pdhg.update_objective_interval = 100
-pdhg.run(1000, verbose=True)
-
-# Show Results
-plt.figure()
-plt.subplot(1,3,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.clim(0,1)
-plt.subplot(1,3,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.clim(0,1)
-plt.subplot(1,3,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.clim(0,1)
-plt.colorbar()
-plt.show()
-
-plt.plot(np.linspace(0,N,M), noisy_data.as_array()[int(N/2),:], label = 'Noisy data')
-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 = '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([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = 0.5 * 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)
-    
-    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,N,M), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,M), u.value[int(N/2),:], label = 'CVX')
-    plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'Truth')
-    
-    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/TV_Denoising/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Poisson.py
deleted file mode 100644
index 1d887c1..0000000
--- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_Poisson.py
+++ /dev/null
@@ -1,212 +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} + \int x - g * log(x)
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator  
-             
-             g: Noisy Data with Poisson Noise
-             
-             
-             Method = 0 ( PDHG - split ) :  K = [ \nabla,
-                                                 Identity]
-                          
-                                                                    
-             Method = 1 (PDHG - explicit ):  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, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, IndicatorBox, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Poisson denoising
-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
-
-# Create noisy data. Apply 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 = 2
-
-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 = KullbackLeibler(noisy_data)    
-    f = BlockFunction(f1, f2)                                        
-    g = IndicatorBox(lower=0)
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  KullbackLeibler(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)
-pdhg.max_iteration = 3000
-pdhg.update_objective_interval = 200
-pdhg.run(3000, verbose=True)
-
-
-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,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/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    
-    u1 = Variable(ig.shape)
-    q = Variable()
-    w = Variable()
-    
-    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([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 = ECOS
-    obj =  Minimize( regulariser +  q)
-    prob = Problem(obj, constraints)
-    result = prob.solve(verbose = True, solver = solver)
-
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u1.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(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,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), 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 (PDHG) {} '.format(pdhg.objective[-1][0]))
-
-
-
-
-
diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_SaltPepper.py
deleted file mode 100644
index c5709c3..0000000
--- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/PDHG_TV_Denoising_SaltPepper.py
+++ /dev/null
@@ -1,213 +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    
-             
-             
-"""
-
-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, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Salt & Pepper denoising
-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
-
-# Create noisy data. Apply Salt & Pepper noise
-n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
-noisy_data = ImageData(n1)
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(15,15))
-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 = 2
-
-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 = L1Norm(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  L1Norm(b = noisy_data)
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
-
-# 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 = 50
-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,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/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([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = pnorm( u - noisy_data.as_array(),1)
-    
-    # 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 = 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,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u.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 (PDHG) {} '.format(pdhg.objective[-1][0]))
-
-
-
-
-
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
deleted file mode 100644
index e74e1c6..0000000
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
+++ /dev/null
@@ -1,194 +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) Tomography 2D using PDHG algorithm:
-
-
-Problem:     min_{x>0} \alpha * ||\nabla x - w||_{2,1} + \beta *  || E w ||_{2,1} +
-                       \frac{1}{2}||Au - g||^{2}
-                     
-            min_{u>0} \alpha * ||\nabla u - w||_{2,1} + \beta *  || E w ||_{2,1} +
-                      int A u - g log(Au + \eta)                     
-
-             \alpha: Regularization parameter
-             \beta: Regularization parameter
-             
-             \nabla: Gradient operator 
-              E: Symmetrized Gradient operator
-              A: System Matrix
-             
-             g: Noisy Sinogram 
-                          
-              K = [ \nabla, - Identity
-                   ZeroOperator, E 
-                   A, ZeroOperator]
-                                         
-"""
-
-from ccpi.framework import AcquisitionGeometry, AcquisitionData, 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, \
-                                    SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import IndicatorBox, KullbackLeibler, ZeroFunction,\
-                      MixedL21Norm, BlockFunction, L2NormSquared
-
-from ccpi.astra.ops import AstraProjectorSimple
-import os, sys
-
-
-import tomophantom
-from tomophantom import TomoP2D
-
-# user supplied input
-if len(sys.argv) > 1:
-    which_noise = int(sys.argv[1])
-else:
-    which_noise = 1 
-    
-# Load Piecewise smooth Shepp-Logan phantom 
-model = 2 # select a model number from the library
-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 phantom (2D)
-phantom_2D = TomoP2D.Model(model, N, path_library2D)
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-data = ImageData(phantom_2D)
-
-#Create Acquisition Data 
-detectors = N
-angles = np.linspace(0, np.pi, N)
-ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-
-#device = input('Available device: GPU==1 / CPU==0 ')
-device = '1'
-if device=='1':
-    dev = 'gpu'
-else:
-    dev = 'cpu'
-
-Aop = AstraProjectorSimple(ig, ag, 'cpu')
-sin = Aop.direct(data)
-
-# Create noisy sinogram.
-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
-op11 = Gradient(ig)
-op12 = Identity(op11.range_geometry())
-
-op22 = SymmetrizedGradient(op11.domain_geometry())    
-op21 = ZeroOperator(ig, op22.range_geometry())
-    
-op31 = Aop
-op32 = ZeroOperator(op22.domain_geometry(), ag)
-
-operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) 
-normK = operator.norm()
-
-# Create functions
-if noise == 'poisson':
-    alpha = 3
-    beta = 6
-    f3 = KullbackLeibler(noisy_data)    
-    g =  BlockFunction(IndicatorBox(lower=0), ZeroFunction()) 
-    
-    # Primal & dual stepsizes
-    sigma = 1
-    tau = 1/(sigma*normK**2)    
-    
-elif noise == 'gaussian':   
-    alpha = 20
-    beta = 50
-    f3 = 0.5 * L2NormSquared(b=noisy_data)                                         
-    g = BlockFunction(ZeroFunction(), ZeroFunction())
-    
-    # Primal & dual stepsizes
-    sigma = 10
-    tau = 1/(sigma*normK**2)     
-    
-f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm()  
-f = BlockFunction(f1, f2, f3)         
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 3000
-pdhg.update_objective_interval = 500
-pdhg.run(3000)
-#%%
-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()[0].as_array())
-plt.title('TGV 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()[0].as_array()[:, int(N/2)], label = 'TGV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py
deleted file mode 100644
index f179e95..0000000
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/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/Tomo/PDHG_TV_Tomo2D_gaussian.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_gaussian.py
deleted file mode 100644
index bd1bb69..0000000
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_gaussian.py
+++ /dev/null
@@ -1,212 +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_x, x>0  \alpha * ||\nabla x||_{2,1} + \frac{1}{2}||Ax - g||^{2}
-
-             \nabla: Gradient operator 
-             
-             A: Projection Matrix
-             g: Noisy sinogram corrupted with Gaussian 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
-
-from ccpi.astra.ops import AstraProjectorSimple
-
-
-
-# Create phantom for TV 2D tomography 
-N = 100
-x = np.zeros((N,N))
-x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-
-data = ImageData(x)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-
-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, dev)
-sin = Aop.direct(data)
-
-# Create noisy data. Apply Poisson noise
-n1 = np.random.normal(0, 3, size=ag.shape)
-noisy_data = AcquisitionData(n1 + sin.as_array(), 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 = 50
-
-# Create operators
-op1 = Gradient(ig)
-op2 = Aop
-
-# 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()
-
-# Compute operator Norm
-normK = operator.norm()
-    
-# Primal & dual stepsizes
-sigma = 10
-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 = 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,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-#%% Check with CVX solution
-
-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)
-    
-    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), 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('strip', proj_geom, vol_geom)
-
-    matrix_id = astra.projector.matrix(proj_id)
-
-    ProjMat = astra.matrix.get(matrix_id)
-    
-    tmp = noisy_data.as_array().ravel()
-    
-    fidelity = 0.5 * sum_squares(ProjMat * u - tmp)
-
-    solver = MOSEK
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)    
-         
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - np.reshape(u.value, (N,N) ))
-           
-    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(np.reshape(u.value, (N, N)))
-    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,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N,N) )[int(N/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/Tomo/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py
deleted file mode 100644
index e7e33cb..0000000
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py
+++ /dev/null
@@ -1,230 +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} + 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 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 KullbackLeibler, \
-                      MixedL21Norm, BlockFunction, IndicatorBox
-
-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 = 50
-M = 50
-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 = np.random.poisson(eta + sin.as_array())
-
-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 = 2
-
-# Create operators
-op1 = Gradient(ig)
-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 = IndicatorBox(lower=0)
-
-
-# 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 = 3000
-pdhg.update_objective_interval = 500
-pdhg.run(3000, verbose = True)
-
-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,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-#%% Check with CVX solution
-#
-#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 = alpha * sum(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('strip', proj_geom, vol_geom)
-#
-#    matrix_id = astra.projector.matrix(proj_id)
-#
-#    ProjMat = astra.matrix.get(matrix_id)
-#    
-#    tmp = noisy_data.as_array().ravel()
-#    
-#    fidelity = sum(kl_div(tmp, ProjMat * u))    
-#    
-#    constraints = [q>=fidelity, u>=0]   
-#    solver = SCS
-#    obj =  Minimize( regulariser +  q)
-#    prob = Problem(obj, constraints)
-#    result = prob.solve(verbose = True, solver = solver)    
-#             
-#    diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N)))
-#           
-#    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(np.reshape(u.value, (N, N)))
-#    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,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-#    plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N, N))[int(N/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/Tomo/phantom.mat b/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat
deleted file mode 100755
index c465bbe..0000000
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat
+++ /dev/null