From 180ceb3c6b480de12ca480fd0a3d06e326a68db5 Mon Sep 17 00:00:00 2001
From: epapoutsellis <epapoutsellis@gmail.com>
Date: Wed, 26 Jun 2019 17:30:02 +0100
Subject: fix cgls imple
---
.../ccpi/optimisation/algorithms/Algorithm.py | 2 +
.../Python/ccpi/optimisation/algorithms/CGLS.py | 53 +++--
.../Python/demos/CGLS_examples/CGLS_Tikhonov.py | 153 ---------------
.../Python/demos/CGLS_examples/CGLS_Tomography.py | 214 +++++++++++++++++++++
.../Python/demos/CGLS_examples/LinearSystem.py | 82 ++++++++
.../CGLS_examples/Regularised_CGLS_Denoising.py | 146 ++++++++++++++
.../demos/PDHG_examples/GatherAll/phantom.mat | Bin 5583 -> 0 bytes
7 files changed, 486 insertions(+), 164 deletions(-)
delete mode 100644 Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
create mode 100644 Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
create mode 100644 Wrappers/Python/demos/CGLS_examples/LinearSystem.py
create mode 100644 Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py
delete mode 100755 Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat
(limited to 'Wrappers/Python')
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
index c62d0ea..bf851d7 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
@@ -84,6 +84,8 @@ class Algorithm(object):
calling this method triggers update and update_objective
'''
if self.should_stop():
+ self.update_objective()
+ print(self.verbose_output())
raise StopIteration()
else:
time0 = time.time()
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
index 28b19f6..6cf4f06 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
@@ -34,19 +34,31 @@ class CGLS(Algorithm):
x_init: initial guess
operator: operator for forward/backward projections
data: data to operate on
+ tolerance: tolerance to stop algorithm
+
+ Reference:
+ https://web.stanford.edu/group/SOL/software/cgls/
+
'''
+
+
+
def __init__(self, **kwargs):
+
super(CGLS, self).__init__()
self.x = kwargs.get('x_init', None)
self.operator = kwargs.get('operator', None)
self.data = kwargs.get('data', None)
+ self.tolerance = kwargs.get('tolerance', None)
if self.x is not None and self.operator is not None and \
self.data is not None:
- print ("Calling from creator")
self.set_up(x_init =kwargs['x_init'],
operator=kwargs['operator'],
data =kwargs['data'])
+ if self.tolerance is None:
+ self.tolerance = 1e-6
+
def set_up(self, x_init, operator , data ):
self.x = x_init
@@ -55,10 +67,19 @@ class CGLS(Algorithm):
self.p = self.s
self.norms0 = self.s.norm()
+
+ ##
+ self.norms = self.s.norm()
+ ##
+
+
self.gamma = self.norms0**2
self.normx = self.x.norm()
self.xmax = self.normx
- self.configured = True
+
+ n = Norm2Sq(self.operator, self.data)
+ self.loss.append(n(self.x))
+ self.configured = True
# def set_up(self, x_init, operator , data ):
#
@@ -80,15 +101,15 @@ class CGLS(Algorithm):
# self.loss.append(n(x_init))
# self.configured = True
- def update(self):
- self.update_new()
+ #def update(self):
+ # self.update_new()
- def update_new(self):
+ def update(self):
self.q = self.operator.direct(self.p)
delta = self.q.squared_norm()
alpha = self.gamma/delta
-
+
self.x += alpha * self.p
self.r -= alpha * self.q
@@ -102,10 +123,7 @@ class CGLS(Algorithm):
self.normx = self.x.norm()
self.xmax = numpy.maximum(self.xmax, self.normx)
-
-
- if self.gamma<=1e-6:
- raise StopIteration()
+
# def update_new(self):
#
# Ad = self.operator.direct(self.d)
@@ -141,7 +159,20 @@ class CGLS(Algorithm):
raise StopIteration()
self.loss.append(a)
-# def should_stop(self):
+ def should_stop(self):
+
+ flag = (self.norms <= self.norms0 * self.tolerance) or (self.normx * self.tolerance >= 1);
+
+ #if self.gamma<=self.tolerance:
+ if flag == 1 or self.max_iteration_stop_cryterion():
+ print('Tolerance is reached: Iter: {}'.format(self.iteration))
+ self.update_objective()
+ return True
+
+
+ #raise StopIteration()
+
+
# if self.iteration > 0:
# x = self.get_last_objective()
# a = x > 0
diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
deleted file mode 100644
index 63a2254..0000000
--- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
+++ /dev/null
@@ -1,153 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0.txt
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-#
-#=========================================================================
-
-"""
-Compare solutions of PDHG & "Block CGLS" algorithms for
-
-
-Problem: min_x alpha * ||\grad x ||^{2}_{2} + || A x - g ||_{2}^{2}
-
-
- A: Projection operator
- g: Sinogram
-
-"""
-
-
-from ccpi.framework import AcquisitionGeometry, BlockDataContainer, AcquisitionData
-
-import numpy as np
-import numpy
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import CGLS
-from ccpi.optimisation.operators import BlockOperator, Gradient
-
-from ccpi.framework import TestData
-import os, sys
-from ccpi.astra.ops import AstraProjectorSimple
-
-# Load Data
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-N = 64
-M = 64
-data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
-
-ig = data.geometry
-
-detectors = N
-angles = np.linspace(0, np.pi, N, dtype=np.float32)
-
-ag = AcquisitionGeometry('parallel','2D', angles, detectors)
-
-device = input('Available device: GPU==1 / CPU==0 ')
-if device=='1':
- dev = 'gpu'
-else:
- dev = 'cpu'
-
-Aop = AstraProjectorSimple(ig, ag, dev)
-sin = Aop.direct(data)
-
-noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,3,ig.shape))
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(10,10))
-plt.subplot(2,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(2,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Setup and run the CGLS algorithm
-alpha = 5
-Grad = Gradient(ig)
-#
-## Form Tikhonov as a Block CGLS structure
-op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1))
-block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate())
-#
-x_init = ig.allocate()
-cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data)
-cgls.max_iteration = 1000
-cgls.update_objective_interval = 200
-cgls.run(1000,verbose=True)
-
-# Show results
-plt.figure(figsize=(5,5))
-plt.imshow(cgls.get_output().as_array())
-plt.title('CGLS reconstruction')
-plt.colorbar()
-plt.show()
-
-#%%
-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*M)
- #q = Variable()
-
- DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
- DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-
- regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(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 = sum_squares( ProjMat * u - noisy_data.as_array().ravel())
-
- solver = SCS
- obj = Minimize( regulariser + fidelity)
- prob = Problem(obj, constraints)
- result = prob.solve(verbose = True, solver = solver)
-
-
-
-
-
-
-
-
diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
new file mode 100644
index 0000000..abcb26c
--- /dev/null
+++ b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
@@ -0,0 +1,214 @@
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+#=========================================================================
+
+"""
+Conjugate Gradient for (Regularized) Least Squares for Tomography
+
+
+Problem: min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2}
+
+ A: Projection operator
+ g: Sinogram
+ L: Identity or Gradient Operator
+
+"""
+
+
+from ccpi.framework import ImageGeometry, ImageData, \
+ AcquisitionGeometry, BlockDataContainer, AcquisitionData
+
+import numpy as np
+import numpy
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import CGLS
+from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
+
+import tomophantom
+from tomophantom import TomoP2D
+from ccpi.astra.operators import AstraProjectorSimple
+import os
+
+
+# Load Shepp-Logan phantom
+model = 1 # select a model number from the library
+N = 64 # set dimension of the phantom
+path = os.path.dirname(tomophantom.__file__)
+path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
+phantom_2D = TomoP2D.Model(model, N, path_library2D)
+
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+data = ImageData(phantom_2D)
+
+detectors = N
+angles = np.linspace(0, np.pi, 90, dtype=np.float32)
+
+ag = AcquisitionGeometry('parallel','2D', angles, detectors)
+
+device = input('Available device: GPU==1 / CPU==0 ')
+
+if device =='1':
+ dev = 'gpu'
+else:
+ dev = 'cpu'
+
+Aop = AstraProjectorSimple(ig, ag, dev)
+sin = Aop.direct(data)
+
+np.random.seed(10)
+noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,1,ag.shape))
+#noisy_data = AcquisitionData( sin.as_array() )
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,10))
+plt.subplot(2,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(2,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+
+# Setup and run the simple CGLS algorithm
+x_init = ig.allocate()
+
+cgls1 = CGLS(x_init = x_init, operator = Aop, data = noisy_data)
+cgls1.max_iteration = 20
+cgls1.update_objective_interval = 5
+cgls1.run(20, verbose = True)
+
+
+# Setup and run the regularised CGLS algorithm (Tikhonov with Identity)
+
+x_init = ig.allocate()
+alpha1 = 50
+op1 = Identity(ig)
+
+block_op1 = BlockOperator( Aop, alpha1 * op1, shape=(2,1))
+block_data1 = BlockDataContainer(noisy_data, op1.range_geometry().allocate())
+
+cgls2 = CGLS(x_init = x_init, operator = block_op1, data = block_data1)
+cgls2.max_iteration = 200
+cgls2.update_objective_interval = 10
+cgls2.run(200, verbose=True)
+
+# Setup and run the regularised CGLS algorithm (Tikhonov with Gradient)
+
+x_init = ig.allocate()
+alpha2 = 25
+op2 = Gradient(ig)
+
+block_op2 = BlockOperator( Aop, alpha2 * op2, shape=(2,1))
+block_data2 = BlockDataContainer(noisy_data, op2.range_geometry().allocate())
+
+cgls3 = CGLS(x_init=x_init, operator = block_op2, data = block_data2)
+cgls3.max_iteration = 200
+cgls3.update_objective_interval = 5
+cgls3.run(200, verbose = True)
+
+# Show results
+plt.figure(figsize=(8,8))
+
+plt.subplot(2,2,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+
+plt.subplot(2,2,2)
+plt.imshow(cgls1.get_output().as_array())
+plt.title('GGLS')
+
+plt.subplot(2,2,3)
+plt.imshow(cgls2.get_output().as_array())
+plt.title('Regularised GGLS L = {} * Identity'.format(alpha1))
+
+plt.subplot(2,2,4)
+plt.imshow(cgls3.get_output().as_array())
+plt.title('Regularised GGLS L = {} * Gradient'.format(alpha2))
+
+plt.show()
+
+
+
+print('Compare CVX vs Regularised CG with L = Gradient')
+
+from ccpi.optimisation.operators import SparseFiniteDiff
+import astra
+import numpy
+
+try:
+ from cvxpy import *
+ cvx_not_installable = True
+except ImportError:
+ cvx_not_installable = False
+
+
+if cvx_not_installable:
+
+ ##Construct problem
+ u = Variable(N*N)
+ #q = Variable()
+
+ DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+ DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+
+ regulariser = alpha2**2 * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+
+ # create matrix representation for Astra operator
+
+ vol_geom = astra.create_vol_geom(N, N)
+ proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
+
+ proj_id = astra.create_projector('line', proj_geom, vol_geom)
+
+ matrix_id = astra.projector.matrix(proj_id)
+
+ ProjMat = astra.matrix.get(matrix_id)
+
+ fidelity = sum_squares( ProjMat * u - noisy_data.as_array().ravel())
+
+ solver = MOSEK
+ obj = Minimize( regulariser + fidelity)
+ prob = Problem(obj)
+ result = prob.solve(verbose = True, solver = solver)
+
+
+plt.figure(figsize=(10,20))
+
+plt.subplot(1,3,1)
+plt.imshow(np.reshape(u.value, (N, N)))
+plt.colorbar()
+
+plt.subplot(1,3,2)
+plt.imshow(cgls3.get_output().as_array())
+plt.colorbar()
+
+plt.subplot(1,3,3)
+plt.imshow(np.abs(cgls3.get_output().as_array() - np.reshape(u.value, (N, N)) ))
+plt.colorbar()
+
+plt.show()
+
+print('Primal Objective (CVX) {} '.format(obj.value))
+print('Primal Objective (CGLS) {} '.format(cgls3.objective[-1]))
+
diff --git a/Wrappers/Python/demos/CGLS_examples/LinearSystem.py b/Wrappers/Python/demos/CGLS_examples/LinearSystem.py
new file mode 100644
index 0000000..cc398cb
--- /dev/null
+++ b/Wrappers/Python/demos/CGLS_examples/LinearSystem.py
@@ -0,0 +1,82 @@
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+#=========================================================================
+
+import numpy
+from ccpi.optimisation.operators import *
+from ccpi.optimisation.algorithms import *
+from ccpi.optimisation.functions import *
+from ccpi.framework import *
+
+# Problem dimension.
+m = 200
+n = 200
+
+numpy.random.seed(10)
+
+# Create matrix A and data b
+Amat = numpy.asarray( numpy.random.randn(m, n), dtype = numpy.float32)
+bmat = numpy.asarray( numpy.random.randn(m), dtype=numpy.float32)
+
+# Geometries for data and solution
+vgb = VectorGeometry(m)
+vgx = VectorGeometry(n)
+
+b = vgb.allocate(0, dtype=numpy.float32)
+b.fill(bmat)
+A = LinearOperatorMatrix(Amat)
+
+# Setup and run CGLS
+x_init = vgx.allocate()
+
+cgls = CGLS(x_init = x_init, operator = A, data = b)
+cgls.max_iteration = 2000
+cgls.update_objective_interval = 200
+cgls.run(2000, verbose = True)
+
+try:
+ from cvxpy import *
+ cvx_not_installable = True
+except ImportError:
+ cvx_not_installable = False
+
+# Compare with CVX
+x = Variable(n)
+obj = Minimize(sum_squares(A.A*x - b.as_array()))
+prob = Problem(obj)
+
+# choose solver
+if 'MOSEK' in installed_solvers():
+ solver = MOSEK
+else:
+ solver = SCS
+
+result = prob.solve(solver = MOSEK)
+
+
+diff_sol = x.value - cgls.get_output().as_array()
+
+print('Error |CVX - CGLS| = {}'.format(numpy.sum(numpy.abs(diff_sol))))
+print('CVX objective = {}'.format(obj.value))
+print('CGLS objective = {}'.format(cgls.objective[-1]))
+
+
+
+
diff --git a/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py b/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py
new file mode 100644
index 0000000..c124fa1
--- /dev/null
+++ b/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py
@@ -0,0 +1,146 @@
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+#=========================================================================
+
+"""
+Conjugate Gradient for (Regularized) Least Squares for Tomography
+
+
+Problem: min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2}
+
+ A: Identity operator
+ g: Sinogram
+ L: Identity or Gradient Operator
+
+"""
+
+
+from ccpi.framework import ImageGeometry, ImageData, \
+ AcquisitionGeometry, BlockDataContainer, AcquisitionData
+
+import numpy as np
+import numpy
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import CGLS
+from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
+from ccpi.framework import TestData
+
+import os, sys
+
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+data = loader.load(TestData.SHAPES)
+ig = data.geometry
+ag = ig
+
+noisy_data = ImageData(TestData.random_noise(data.as_array(), mode = 'gaussian', seed = 1))
+#noisy_data = ImageData(data.as_array())
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,10))
+plt.subplot(2,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(2,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+# Setup and run the regularised CGLS algorithm (Tikhonov with Gradient)
+x_init = ig.allocate()
+alpha = 2
+op = Gradient(ig)
+
+block_op = BlockOperator( Identity(ig), alpha * op, shape=(2,1))
+block_data = BlockDataContainer(noisy_data, op.range_geometry().allocate())
+
+cgls = CGLS(x_init=x_init, operator = block_op, data = block_data)
+cgls.max_iteration = 200
+cgls.update_objective_interval = 5
+cgls.run(200, verbose = True)
+
+# Show results
+plt.figure(figsize=(20,10))
+plt.subplot(3,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.subplot(3,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy')
+plt.subplot(3,1,3)
+plt.imshow(cgls.get_output().as_array())
+plt.title('Regularised GGLS with Gradient')
+plt.show()
+
+#%%
+
+print('Compare CVX vs Regularised CG with L = Gradient')
+
+from ccpi.optimisation.operators import SparseFiniteDiff
+
+
+try:
+ from cvxpy import *
+ cvx_not_installable = True
+except ImportError:
+ cvx_not_installable = False
+
+
+if cvx_not_installable:
+
+ ##Construct problem
+ u = Variable(ig.shape)
+ #q = Variable()
+
+ DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+ DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+
+ regulariser = alpha**2 * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+
+ fidelity = sum_squares( u - noisy_data.as_array())
+
+ # choose solver
+ if 'MOSEK' in installed_solvers():
+ solver = MOSEK
+ else:
+ solver = SCS
+
+ obj = Minimize( regulariser + fidelity)
+ prob = Problem(obj)
+ result = prob.solve(verbose = True, solver = MOSEK)
+
+
+plt.figure(figsize=(20,20))
+plt.subplot(3,1,1)
+plt.imshow(np.reshape(u.value, ig.shape))
+plt.colorbar()
+plt.subplot(3,1,2)
+plt.imshow(cgls.get_output().as_array())
+plt.colorbar()
+plt.subplot(3,1,3)
+plt.imshow(np.abs(cgls.get_output().as_array() - np.reshape(u.value, ig.shape) ))
+plt.colorbar()
+plt.show()
+
+print('Primal Objective (CVX) {} '.format(obj.value))
+print('Primal Objective (CGLS) {} '.format(cgls.objective[-1]))
+
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat b/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat
deleted file mode 100755
index c465bbe..0000000
Binary files a/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat and /dev/null differ
--
cgit v1.2.3