delete old demos from wip

9 files changed, 0 insertions, 1904 deletions

diff --git a/Wrappers/Python/wip/fista_TV_denoising.py b/Wrappers/Python/wip/fista_TV_denoising.py
deleted file mode 100644
index a9712fc..0000000
--- a/Wrappers/Python/wip/fista_TV_denoising.py
+++ /dev/null
@@ -1,72 +0,0 @@
-#!/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
-
-import numpy as np                           
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, PDHG_old, FISTA
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction
-                      
-from ccpi.optimisation.algs import FISTA                      
-
-from skimage.util import random_noise
-
-from timeit import default_timer as timer
-def dt(steps):
-    return steps[-1] - steps[-2]
-
-# Create phantom for TV 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
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.05, seed=10)
-noisy_data = ImageData(n1)
-
-
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy data')
-plt.show()
-
-# Regularisation Parameter
-alpha = 2
-
-operator = Gradient(ig)
-g = alpha * MixedL21Norm()
-f = 0.5 * L2NormSquared(b = noisy_data)
-    
-x_init = ig.allocate()
-opt = {'niter':2000}
-
-
-x = FISTA(x_init, f, g, opt)
-
-#fista = FISTA()
-#fista.set_up(x_init, f, g, opt )
-#fista.max_iteration = 10
-#
-#fista.run(2000)
-#plt.figure(figsize=(15,15))
-#plt.subplot(3,1,1)
-#plt.imshow(fista.get_output().as_array())
-#plt.title('no memopt class')
-
-
-
diff --git a/Wrappers/Python/wip/pdhg_TGV_denoising.py b/Wrappers/Python/wip/pdhg_TGV_denoising.py
deleted file mode 100644
index 1c570cb..0000000
--- a/Wrappers/Python/wip/pdhg_TGV_denoising.py
+++ /dev/null
@@ -1,230 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 22 14:53:03 2019
-
-@author: evangelos
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, PDHG_old
-
-from ccpi.optimisation.operators import BlockOperator, Identity, \
-                        Gradient, SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-from timeit import default_timer as timer
-#def dt(steps):
-#    return steps[-1] - steps[-2]
-
-# Create phantom for TGV Gaussian 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 = data/data.max()
-
-plt.imshow(data)
-plt.show()
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.005, seed=10)
-noisy_data = ImageData(n1)
-
-
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy data')
-plt.show()
-
-alpha, beta = 0.1, 0.5
-
-#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-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 = 0.5 * L2NormSquared(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(0.5 * L2NormSquared(b = noisy_data), ZeroFunction())
-     
-## Compute operator Norm
-normK = operator.norm()
-#
-## Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-##
-opt = {'niter':500}
-opt1 = {'niter':500, 'memopt': True}
-#
-t1 = timer()
-res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-t2 = timer()
-#
-t3 = timer()
-res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
-t4 = timer()
-#
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(res[0].as_array())
-plt.title('no memopt')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(res1[0].as_array())
-plt.title('memopt')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow((res1[0] - res[0]).abs().as_array())
-plt.title('diff')
-plt.colorbar()
-plt.show()
-
-print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3))
-    
-
-######
-
-#%%
-
-#def update_plot(it_update, objective, x):
-#    
-##    sol = pdhg.get_output()
-#    plt.figure()
-#    plt.imshow(x[0].as_array())
-#    plt.show()
-#    
-#    
-##def stop_criterion(x,)    
-#
-#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-#pdhg.max_iteration = 2000
-#pdhg.update_objective_interval = 100
-#
-#pdhg.run(4000, verbose=False, callback=update_plot)
-
-
-#%%
-
-
-
-
-
-
-
-
-#%% 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 = 0.5 * sum_squares(u - noisy_data.as_array())    
-#    solver = MOSEK
-#
-#    obj =  Minimize( regulariser +  fidelity)
-#    prob = Problem(obj)
-#    result = prob.solve(verbose = True, solver = solver)
-#    
-#    diff_cvx = numpy.abs( res[0].as_array() - u.value )   
-#    
-#    # Show result
-#    plt.figure(figsize=(15,15))
-#    plt.subplot(3,1,1)
-#    plt.imshow(res[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), res[0].as_array()[int(N/2),:], label = 'PDHG')
-#    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-#    plt.legend()   
-#    
-#    print('Primal Objective (CVX) {} '.format(obj.value))
-#    print('Primal Objective (PDHG) {} '.format(primal[-1])) 
-#    print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max()))
-#    
-    
-   
diff --git a/Wrappers/Python/wip/pdhg_TGV_tomography2D.py b/Wrappers/Python/wip/pdhg_TGV_tomography2D.py
deleted file mode 100644
index ee3b089..0000000
--- a/Wrappers/Python/wip/pdhg_TGV_tomography2D.py
+++ /dev/null
@@ -1,200 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 22 14:53:03 2019
-
-@author: evangelos
-"""
-
-from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, PDHG_old
-
-from ccpi.optimisation.operators import BlockOperator, Identity, \
-                        Gradient, SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-from timeit import default_timer as timer
-from ccpi.astra.ops import AstraProjectorSimple
-
-#def dt(steps):
-#    return steps[-1] - steps[-2]
-
-# Create phantom for TGV Gaussian 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())
-
-plt.imshow(data.as_array())
-plt.show()
-
-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)
-Aop = AstraProjectorSimple(ig, ag, 'cpu')
-sin = Aop.direct(data)
-
-plt.imshow(sin.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.as_array())
-plt.title('Noisy Sinogram')
-plt.colorbar()
-plt.show()
-
-#%%
-
-alpha, beta = 20, 50
-
-
-# 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) ) 
-
-
-f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm() 
-f3 = 0.5 * L2NormSquared(b = noisy_data)    
-f = BlockFunction(f1, f2, f3)         
-g = ZeroFunction()
-             
-## Compute operator Norm
-normK = operator.norm()
-#
-## Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-##
-opt = {'niter':5000}
-opt1 = {'niter':5000, 'memopt': True}
-#
-t1 = timer()
-res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-t2 = timer()
-#
-plt.imshow(res[0].as_array())
-plt.show()
-
-
-#t3 = timer()
-#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
-#t4 = timer()
-#
-#plt.figure(figsize=(15,15))
-#plt.subplot(3,1,1)
-#plt.imshow(res[0].as_array())
-#plt.title('no memopt')
-#plt.colorbar()
-#plt.subplot(3,1,2)
-#plt.imshow(res1[0].as_array())
-#plt.title('memopt')
-#plt.colorbar()
-#plt.subplot(3,1,3)
-#plt.imshow((res1[0] - res[0]).abs().as_array())
-#plt.title('diff')
-#plt.colorbar()
-#plt.show()
-#
-#print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3))
-    
-
-#%% 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 = 0.5 * sum_squares(u - noisy_data.as_array())    
-#    solver = MOSEK
-#
-#    obj =  Minimize( regulariser +  fidelity)
-#    prob = Problem(obj)
-#    result = prob.solve(verbose = True, solver = solver)
-#    
-#    diff_cvx = numpy.abs( res[0].as_array() - u.value )   
-#    
-#    # Show result
-#    plt.figure(figsize=(15,15))
-#    plt.subplot(3,1,1)
-#    plt.imshow(res[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), res[0].as_array()[int(N/2),:], label = 'PDHG')
-#    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-#    plt.legend()   
-#    
-#    print('Primal Objective (CVX) {} '.format(obj.value))
-#    print('Primal Objective (PDHG) {} '.format(primal[-1])) 
-#    print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max()))
-    
-    
-   
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py
deleted file mode 100755
index b16e8f9..0000000
--- a/Wrappers/Python/wip/pdhg_TV_denoising.py
+++ /dev/null
@@ -1,204 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 22 14:53:03 2019
-
-@author: evangelos
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-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
-
-from skimage.util import random_noise
-
-from timeit import default_timer as timer
-#def dt(steps):
-#    return steps[-1] - steps[-2]
-
-# Create phantom for TV Gaussian 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
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.05, seed=10)
-noisy_data = ImageData(n1)
-
-
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy data')
-plt.show()
-
-# Regularisation Parameter
-alpha = 0.5
-
-#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-method = '1'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Form Composite Operator
-    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:
-    
-    ###########################################################################
-    #         No Composite #
-    ###########################################################################
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-    
-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:
-    # Compute operator Norm
-    normK = operator.norm()
-    
-    # Primal & dual stepsizes
-    sigma = 1
-    tau = 1/(sigma*normK**2)
-
-
-opt = {'niter':5000}
-opt1 = {'niter':5000, 'memopt': True}
-
-t1 = timer()
-res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-t2 = timer()
-
-t3 = timer()
-res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
-t4 = timer()
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(res.as_array())
-plt.title('no memopt')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(res1.as_array())
-plt.title('memopt')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow((res1 - res).abs().as_array())
-plt.title('diff')
-plt.colorbar()
-plt.show()
-# 
-plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt')
-plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt')
-plt.legend()
-plt.show()
-#
-print ("Time: No memopt in {}s, \n Time: Memopt in  {}s ".format(t2-t1, t4 -t3))
-diff = (res1 - res).abs().as_array().max()
-#
-print(" Max of abs difference is {}".format(diff))
-
-
-#%% 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 = solver)
-    
-    diff_cvx = numpy.abs( res.as_array() - u.value )
-    
-    # Show result
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(res.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), res1.as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-    plt.legend()
-    
-
-    
-    
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(primal[-1]))
-
-
-
-
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising3D.py b/Wrappers/Python/wip/pdhg_TV_denoising3D.py
deleted file mode 100644
index 06ecfa2..0000000
--- a/Wrappers/Python/wip/pdhg_TV_denoising3D.py
+++ /dev/null
@@ -1,360 +0,0 @@
-#!/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
-
-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, FunctionOperatorComposition, BlockFunction
-
-from skimage.util import random_noise
-
-from timeit import default_timer as timer
-def dt(steps):
-    return steps[-1] - steps[-2]
-
-#%%
-
-# Create phantom for TV 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_size = 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_size x N_size x N_size phantom (3D)
-phantom_tm = TomoP3D.Model(model, N_size, path_library3D)
-#toc=timeit.default_timer()
-#Run_time = toc - tic
-#print("Phantom has been built in {} seconds".format(Run_time))
-#
-#sliceSel = int(0.5*N_size)
-##plt.gray()
-#plt.figure() 
-#plt.subplot(131)
-#plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1)
-#plt.title('3D Phantom, axial view')
-#
-#plt.subplot(132)
-#plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1)
-#plt.title('3D Phantom, coronal view')
-#
-#plt.subplot(133)
-#plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1)
-#plt.title('3D Phantom, sagittal view')
-#plt.show()
-
-#%%
-
-N = N_size
-ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N)
-
-n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10)
-noisy_data = ImageData(n1)
-#plt.imshow(noisy_data.as_array()[:,:,32])
-
-#%%
-
-# Regularisation Parameter
-alpha = 0.02
-
-#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig)
-
-    # Form Composite Operator
-    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:
-    
-    ###########################################################################
-    #         No Composite #
-    ###########################################################################
-    operator = Gradient(ig)
-    f = alpha * FunctionOperatorComposition(operator, MixedL21Norm())
-    g = L2NormSquared(b = noisy_data)
-    
-    ###########################################################################
-#%%
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-opt = {'niter':2000}
-opt1 = {'niter':2000, 'memopt': True}
-
-#t1 = timer()
-#res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-#t2 = timer()
-
-
-t3 = timer()
-res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
-t4 = timer()
-
-#import cProfile
-#cProfile.run('res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) ')
-###
-print ("No memopt in {}s, memopt in  {}s ".format(t2-t1, t4 -t3))
-#
-##    
-##%%
-#
-#plt.figure(figsize=(10,10)) 
-#plt.subplot(311)
-#plt.imshow(res1.as_array()[sliceSel,:,:])
-#plt.colorbar()
-#plt.title('3D Phantom, axial view')
-#
-#plt.subplot(312)
-#plt.imshow(res1.as_array()[:,sliceSel,:])
-#plt.colorbar()
-#plt.title('3D Phantom, coronal view')
-#
-#plt.subplot(313)
-#plt.imshow(res1.as_array()[:,:,sliceSel])
-#plt.colorbar()
-#plt.title('3D Phantom, sagittal view')
-#plt.show()
-#
-#plt.figure(figsize=(10,10)) 
-#plt.subplot(311)
-#plt.imshow(res.as_array()[sliceSel,:,:])
-#plt.colorbar()
-#plt.title('3D Phantom, axial view')
-#
-#plt.subplot(312)
-#plt.imshow(res.as_array()[:,sliceSel,:])
-#plt.colorbar()
-#plt.title('3D Phantom, coronal view')
-#
-#plt.subplot(313)
-#plt.imshow(res.as_array()[:,:,sliceSel])
-#plt.colorbar()
-#plt.title('3D Phantom, sagittal view')
-#plt.show()
-#
-#diff  =  (res1 - res).abs()
-#
-#plt.figure(figsize=(10,10)) 
-#plt.subplot(311)
-#plt.imshow(diff.as_array()[sliceSel,:,:])
-#plt.colorbar()
-#plt.title('3D Phantom, axial view')
-#
-#plt.subplot(312)
-#plt.imshow(diff.as_array()[:,sliceSel,:])
-#plt.colorbar()
-#plt.title('3D Phantom, coronal view')
-#
-#plt.subplot(313)
-#plt.imshow(diff.as_array()[:,:,sliceSel])
-#plt.colorbar()
-#plt.title('3D Phantom, sagittal view')
-#plt.show()
-#
-#
-#
-#
-##%%
-#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-#pdhg.max_iteration = 2000
-#pdhg.update_objective_interval = 100
-####
-#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-#pdhgo.max_iteration = 2000
-#pdhgo.update_objective_interval = 100
-####
-#steps = [timer()]
-#pdhgo.run(2000)
-#steps.append(timer())
-#t1 = dt(steps)
-##
-#pdhg.run(2000)
-#steps.append(timer())
-#t2 = dt(steps)
-#
-#print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1))
-#res = pdhg.get_output()
-#res1 = pdhgo.get_output()
-
-#%%
-#plt.figure(figsize=(15,15))
-#plt.subplot(3,1,1)
-#plt.imshow(res.as_array())
-#plt.title('no memopt')
-#plt.colorbar()
-#plt.subplot(3,1,2)
-#plt.imshow(res1.as_array())
-#plt.title('memopt')
-#plt.colorbar()
-#plt.subplot(3,1,3)
-#plt.imshow((res1 - res).abs().as_array())
-#plt.title('diff')
-#plt.colorbar()
-#plt.show()
-
-
-#plt.figure(figsize=(15,15))
-#plt.subplot(3,1,1)
-#plt.imshow(pdhg.get_output().as_array())
-#plt.title('no memopt class')
-#plt.colorbar()
-#plt.subplot(3,1,2)
-#plt.imshow(res.as_array())
-#plt.title('no memopt')
-#plt.colorbar()
-#plt.subplot(3,1,3)
-#plt.imshow((pdhg.get_output() - res).abs().as_array())
-#plt.title('diff')
-#plt.colorbar()
-#plt.show()
-#    
-#    
-#
-#plt.figure(figsize=(15,15))
-#plt.subplot(3,1,1)
-#plt.imshow(pdhgo.get_output().as_array())
-#plt.title('memopt class')
-#plt.colorbar()
-#plt.subplot(3,1,2)
-#plt.imshow(res1.as_array())
-#plt.title('no memopt')
-#plt.colorbar()
-#plt.subplot(3,1,3)
-#plt.imshow((pdhgo.get_output() - res1).abs().as_array())
-#plt.title('diff')
-#plt.colorbar()
-#plt.show()
-    
-
-    
-
-    
-#    print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1))
-#    res = pdhg.get_output()
-#    res1 = pdhgo.get_output()
-#    
-#    diff = (res-res1)
-#    print ("diff norm {} max {}".format(diff.norm(), diff.abs().as_array().max()))
-#    print ("Sum ( abs(diff) )  {}".format(diff.abs().sum()))
-#    
-#    
-#    plt.figure(figsize=(5,5))
-#    plt.subplot(1,3,1)
-#    plt.imshow(res.as_array())
-#    plt.colorbar()
-#    #plt.show()
-#     
-#    #plt.figure(figsize=(5,5))
-#    plt.subplot(1,3,2)
-#    plt.imshow(res1.as_array())
-#    plt.colorbar()
-    
-#plt.show()
-
-
-
-#=======
-## opt = {'niter':2000, 'memopt': True}
-#
-## res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-# 
-#>>>>>>> origin/pdhg_fix
-#
-#
-## opt = {'niter':2000, 'memopt': False}
-## res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-#
-## plt.figure(figsize=(5,5))
-## plt.subplot(1,3,1)
-## plt.imshow(res.as_array())
-## plt.title('memopt')
-## plt.colorbar()
-## plt.subplot(1,3,2)
-## plt.imshow(res1.as_array())
-## plt.title('no memopt')
-## plt.colorbar()
-## plt.subplot(1,3,3)
-## plt.imshow((res1 - res).abs().as_array())
-## plt.title('diff')
-## plt.colorbar()
-## plt.show()
-#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-#pdhg.max_iteration = 2000
-#pdhg.update_objective_interval = 100
-#
-#
-#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-#pdhgo.max_iteration = 2000
-#pdhgo.update_objective_interval = 100
-#
-#steps = [timer()]
-#pdhgo.run(200)
-#steps.append(timer())
-#t1 = dt(steps)
-#
-#pdhg.run(200)
-#steps.append(timer())
-#t2 = dt(steps)
-#
-#print ("Time difference {} {} {}".format(t1,t2,t2-t1))
-#sol = pdhg.get_output().as_array()
-##sol = result.as_array()
-##
-#fig = plt.figure()
-#plt.subplot(1,3,1)
-#plt.imshow(noisy_data.as_array())
-#plt.colorbar()
-#plt.subplot(1,3,2)
-#plt.imshow(sol)
-#plt.colorbar()
-#plt.subplot(1,3,3)
-#plt.imshow(pdhgo.get_output().as_array())
-#plt.colorbar()
-#
-#plt.show()
-###
-##
-####
-##plt.plot(np.linspace(0,N,N), data[int(N/2),:], label = 'GTruth')
-##plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
-##plt.legend()
-##plt.show()
-#
-#
-##%%
-##
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py
deleted file mode 100644
index bd330fc..0000000
--- a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py
+++ /dev/null
@@ -1,268 +0,0 @@
-#!/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
-
-import numpy as np 
-import numpy                          
-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, L1Norm, \
-                      MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction
-
-from skimage.util import random_noise
-
-from timeit import default_timer as timer
-
-
-
-# ############################################################################
-# Create phantom for TV 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
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data, mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
-noisy_data = ImageData(n1)
-
-plt.imshow(noisy_data.as_array())
-plt.colorbar()
-plt.show()
-
-#%%
-
-# Regularisation Parameter
-alpha = 2
-
-#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-method = '0'
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    f1 = alpha * MixedL21Norm()
-    f2 = L1Norm(b = noisy_data)
-    
-    f = BlockFunction(f1, f2 )                                        
-    g = ZeroFunction()
-    
-else:
-    
-    ###########################################################################
-    #         No Composite #
-    ###########################################################################
-    operator = Gradient(ig)
-    f = alpha * MixedL21Norm()
-    g = L1Norm(b = noisy_data)
-    ###########################################################################
-#%%
-    
-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:
-    # Compute operator Norm
-    normK = operator.norm()
-    
-    # Primal & dual stepsizes
-    sigma = 1
-    tau = 1/(sigma*normK**2)
-
-
-opt = {'niter':5000}
-opt1 = {'niter':5000, 'memopt': True}
-
-t1 = timer()
-res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-t2 = timer()
-
-
-t3 = timer()
-res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
-t4 = timer()
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(res.as_array())
-plt.title('no memopt')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(res1.as_array())
-plt.title('memopt')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow((res1 - res).abs().as_array())
-plt.title('diff')
-plt.colorbar()
-plt.show()
-# 
-plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt')
-plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt')
-plt.legend()
-plt.show()
-#
-print ("Time: No memopt in {}s, \n Time: Memopt in  {}s ".format(t2-t1, t4 -t3))
-diff = (res1 - res).abs().as_array().max()
-#
-print(" Max of abs difference is {}".format(diff))
-
-#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-#pdhg.max_iteration = 2000
-#pdhg.update_objective_interval = 10
-#
-#pdhg.run(2000)
-
-    
-
-#sol = pdhg.get_output().as_array()
-##sol = result.as_array()
-##
-#fig = plt.figure()
-#plt.subplot(1,2,1)
-#plt.imshow(noisy_data.as_array())
-##plt.colorbar()
-#plt.subplot(1,2,2)
-#plt.imshow(sol)
-##plt.colorbar()
-#plt.show()
-##
-
-##
-#plt.plot(np.linspace(0,N,N), data[int(N/2),:], label = 'GTruth')
-#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
-#plt.legend()
-#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)
-    
-    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( res.as_array() - u.value )
-    
-# Show result
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(res.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), res1.as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-    plt.legend()
-    
-
-    
-    
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(primal[-1]))    
-    
-    
-
-#try_cvx = input("Do you want CVX comparison (0/1)")
-#
-#if try_cvx=='0':
-#
-#    from cvxpy import *
-#    import sys
-#    sys.path.insert(0,'/Users/evangelos/Desktop/Projects/CCPi/CCPi-Framework/Wrappers/Python/ccpi/optimisation/cvx_scripts')
-#    from cvx_functions import TV_cvx
-#
-#    u = Variable((N, N))
-#    fidelity = pnorm( u - noisy_data.as_array(),1)
-#    regulariser = alpha * TV_cvx(u)
-#    solver = MOSEK
-#    obj =  Minimize( regulariser +  fidelity)
-#    constraints = []
-#    prob = Problem(obj, constraints)
-#
-#    # Choose solver (SCS is fast but less accurate than MOSEK)
-#    result = prob.solve(verbose = True, solver = solver)
-#
-#    print('Objective value is {} '.format(obj.value))
-#
-#    diff_pdhg_cvx = np.abs(u.value - res.as_array())
-#    plt.imshow(diff_pdhg_cvx)
-#    plt.colorbar()
-#    plt.title('|CVX-PDHG|')        
-#    plt.show()
-#
-#    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-#    plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'PDHG')
-#    plt.legend()
-#    plt.show()
-#
-#else:
-#    print('No CVX solution available')
-
-
-
-
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
deleted file mode 100644
index e123739..0000000
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py
+++ /dev/null
@@ -1,231 +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, 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, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.operators import AstraProjectorSimple
-from timeit import default_timer as timer
-
-
-#%%###############################################################################
-# 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))
-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)
-Aop = AstraProjectorSimple(ig, ag, 'cpu')
-sin = Aop.direct(data)
-
-plt.imshow(sin.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.as_array())
-plt.title('Noisy Sinogram')
-plt.colorbar()
-plt.show()
-
-# Create operators
-op1 = Gradient(ig)
-op2 = Aop
-
-# Form Composite Operator
-operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-alpha = 10
-f = BlockFunction( alpha * MixedL21Norm(), \
-                   0.5 * L2NormSquared(b = noisy_data) )
-g = ZeroFunction()
-
-# Compute operator Norm
-normK = operator.norm()
-
-## Primal & dual stepsizes
-diag_precon = True 
-
-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:  
-    normK = operator.norm()
-    sigma = 1
-    tau = 1/(sigma*normK**2)
-
-# Compute operator Norm
-
-
-# Primal & dual stepsizes
-
-opt = {'niter':200}
-opt1 = {'niter':200, 'memopt': True}
-
-t1 = timer()
-res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-t2 = timer()
-
-
-t3 = timer()
-res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
-t4 = timer()
-#
-#
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(res.as_array())
-plt.title('no memopt')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(res1.as_array())
-plt.title('memopt')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow((res1 - res).abs().as_array())
-plt.title('diff')
-plt.colorbar()
-plt.show()
-# 
-plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt')
-plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt')
-plt.legend()
-plt.show()
-#
-print ("Time: No memopt in {}s, \n Time: Memopt in  {}s ".format(t2-t1, t4 -t3))
-diff = (res1 - res).abs().as_array().max()
-#
-print(" Max of abs difference is {}".format(diff))
-
-
-#%% 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:
-#    
-#
-#    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)
-#    
-#    fidelity = 0.5 * sum_squares(ProjMat * u - noisy_data.as_array().ravel())
-#    
-#    # 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)    
-#
-##%%
-#    
-#    u_rs = np.reshape(u.value, (N,N))
-#        
-#    diff_cvx = numpy.abs( res.as_array() - u_rs )
-#    
-#    # Show result
-#    plt.figure(figsize=(15,15))
-#    plt.subplot(3,1,1)
-#    plt.imshow(res.as_array())
-#    plt.title('PDHG solution')
-#    plt.colorbar()
-#    
-#    plt.subplot(3,1,2)
-#    plt.imshow(u_rs)
-#    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), res1.as_array()[int(N/2),:], label = 'PDHG')
-#    plt.plot(np.linspace(0,N,N), u_rs[int(N/2),:], label = 'CVX')
-#    plt.legend()
-#    
-#    
-#    print('Primal Objective (CVX) {} '.format(obj.value))
-#    print('Primal Objective (PDHG) {} '.format(primal[-1]))
\ No newline at end of file
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
deleted file mode 100644
index 5423b22..0000000
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
+++ /dev/null
@@ -1,152 +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.astra.ops import AstraProjectorSimple, AstraProjectorMC
-from skimage.util import random_noise
-
-
-#%%###############################################################################
-# Create phantom for TV tomography
-
-import numpy as np
-import matplotlib.pyplot as plt
-import os
-import tomophantom
-from tomophantom import TomoP2D
-
-model = 102  # note that the selected model is temporal (2D + time)
-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 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
-
-#N = 150
-#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
-
-#%%
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
-data = ImageData(phantom_2Dt, geometry=ig)
-
-
-
-detectors = 150
-angles = np.linspace(0,np.pi,100)
-
-ag = AcquisitionGeometry('parallel','2D',angles, detectors, channels = np.shape(phantom_2Dt)[0])
-Aop = AstraProjectorMC(ig, ag, 'gpu')
-sin = Aop.direct(data)
-
-plt.imshow(sin.as_array()[10])
-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.as_array()[10])
-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()
-
-# Compute operator Norm
-normK = operator.norm()
-
-## Primal & dual stepsizes
-
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-#sigma = 1/normK
-#tau = 1/normK
-
-opt = {'niter':2000}
-
-res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
- 
-plt.figure(figsize=(5,5))
-plt.imshow(res.as_array())
-plt.colorbar()
-plt.show()
-
-#sigma = 10
-#tau = 1/(sigma*normK**2)
-#
-#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-#pdhg.max_iteration = 5000
-#pdhg.update_objective_interval = 20
-#
-#pdhg.run(5000)
-#
-##%%
-#sol = pdhg.get_output().as_array()
-#fig = plt.figure()
-#plt.subplot(1,2,1)
-#plt.imshow(noisy_data.as_array())
-##plt.colorbar()
-#plt.subplot(1,2,2)
-#plt.imshow(sol)
-##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), sol[int(N/2),:], label = 'Recon')
-plt.legend()
-plt.show()
-
-
diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py
deleted file mode 100644
index cb37598..0000000
--- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py
+++ /dev/null
@@ -1,187 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 22 14:53:03 2019
-
-@author: evangelos
-"""
-
-import numpy as np 
-import numpy                        
-import matplotlib.pyplot as plt
-
-from ccpi.framework import ImageData, ImageGeometry
-
-from ccpi.optimisation.algorithms import PDHG, PDHG_old
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-                 
-
-from skimage.util import random_noise
-from timeit import default_timer as timer
-
-
-
-# ############################################################################
-# Create phantom for TV Poisson denoising
-
-N = 200
-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
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data, mode = 'poisson', seed = 10)
-noisy_data = ImageData(n1)
-
-plt.imshow(noisy_data.as_array())
-plt.colorbar()
-plt.show()
-
-# Regularisation Parameter
-alpha = 2
-
-#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-
-method = '0'
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Form Composite Operator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    f1 = alpha * MixedL21Norm()
-    f2 = KullbackLeibler(noisy_data)
-    
-    f = BlockFunction(f1, f2 )                                        
-    g = ZeroFunction()
-    
-else:
-    
-    ###########################################################################
-    #         No Composite #
-    ###########################################################################
-    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)
-
-opt = {'niter':2000}
-opt1 = {'niter':2000, 'memopt': True}
-
-t1 = timer()
-res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
-t2 = timer()
-
-t3 = timer()
-res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
-t4 = timer()
-
-print(pdgap[-1])
-
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(res.as_array())
-plt.title('no memopt')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(res1.as_array())
-plt.title('memopt')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow((res1 - res).abs().as_array())
-plt.title('diff')
-plt.colorbar()
-plt.show()
-# 
-plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt')
-plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt')
-plt.legend()
-plt.show()
-
-print ("Time: No memopt in {}s, \n Time: Memopt in  {}s ".format(t2-t1, t4 -t3))
-diff = (res1 - res).abs().as_array().max()
-
-print(" Max of abs difference is {}".format(diff))
-
-
-#%% 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')
-    
-    # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0))
-    
-    fidelity = sum( u1 - multiply(noisy_data.as_array(), log(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( res.as_array() - u1.value )
-    
-    # Show result
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(res.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), res1.as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX')
-    plt.legend()  
-    
-    
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(primal[-1]))
-    
-  
-
-
-