diff options
Diffstat (limited to 'Wrappers/Python')
5 files changed, 82 insertions, 77 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py index ace0ba7..f585c9b 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py +++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py @@ -60,7 +60,7 @@ class IndicatorBox(Function): if out is None: return (x.maximum(self.lower)).minimum(self.upper) - else: + else: x.maximum(self.lower, out=out) out.minimum(self.upper, out=out) diff --git a/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py index 49d4db6..19cf684 100644 --- a/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py +++ b/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py @@ -27,7 +27,7 @@ from ccpi.optimisation.algorithms import PDHG from ccpi.optimisation.operators import BlockOperator, Gradient, Identity, \ SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ +from ccpi.optimisation.functions import ZeroFunction, IndicatorBox, KullbackLeibler, \ MixedL21Norm, BlockFunction from ccpi.astra.ops import AstraProjectorSimple @@ -84,13 +84,16 @@ f1 = alpha * MixedL21Norm() f2 = beta * MixedL21Norm() f3 = KullbackLeibler(noisy_data) f = BlockFunction(f1, f2, f3) -g = ZeroFunction() + +g = BlockFunction(-1 * IndicatorBox(lower=0), ZeroFunction()) +#g = IndicatorBox(lower=0) +#g = ZeroFunction() # Compute operator Norm normK = operator.norm() # Primal & dual stepsizes -sigma = 1 +sigma = 10 tau = 1/(sigma*normK**2) diff --git a/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py index 63f8955..22972da 100644 --- a/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py +++ b/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py @@ -50,6 +50,8 @@ from ccpi.optimisation.operators import BlockOperator, Gradient from ccpi.optimisation.functions import IndicatorBox, L2NormSquared, BlockFunction from skimage.util import random_noise 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')) @@ -97,9 +99,9 @@ plt.title('Noisy Data') plt.colorbar() plt.show() -#%% + # Regularisation Parameter -alpha = 100 +alpha = 500 # Create operators op1 = Gradient(ig) diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py index 57f6fcd..15c0a05 100644 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py @@ -28,7 +28,7 @@ Problem: min_{x} \alpha * ||\nabla x - w||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} \alpha: Regularization parameter - \alpha: Regularization parameter + \beta: Regularization parameter \nabla: Gradient operator E: Symmetrized Gradient operator diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py index 8f39f86..c44f393 100644 --- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py @@ -153,74 +153,74 @@ plt.title('Middle Line Profiles') plt.show() -##%% Check with CVX solution +#%% 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 +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 |