From 3e0b5a62f3264c3c043c0b219a86b56aadd0d61e Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 17 Apr 2019 11:53:35 +0100 Subject: wip tests and demos --- .../Python/ccpi/optimisation/algorithms/PDHG.py | 35 +++++------ .../ccpi/optimisation/functions/KullbackLeibler.py | 59 +++++++++++------- Wrappers/Python/wip/pdhg_TV_denoising.py | 69 ++++++++++++++++++++-- Wrappers/Python/wip/pdhg_tv_denoising_poisson.py | 66 ++++++++++++++++++--- 4 files changed, 178 insertions(+), 51 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index a165e55..bc080f8 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -166,15 +166,15 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): y_old.fill(y) - if i%100==0: - - p1 = f(operator.direct(x)) + g(x) - d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) ) - primal.append(p1) - dual.append(d1) - pdgap.append(p1-d1) - print(p1, d1, p1-d1) - +# if i%10==0: +## +# p1 = f(operator.direct(x)) + g(x) +## print(p1) +# d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) ) +# primal.append(p1) +# dual.append(d1) +# pdgap.append(p1-d1) +# print(p1, d1, p1-d1) else: @@ -196,14 +196,15 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): x_old.fill(x) y_old.fill(y) - if i%100==0: - - p1 = f(operator.direct(x)) + g(x) - d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) ) - primal.append(p1) - dual.append(d1) - pdgap.append(p1-d1) - print(p1, d1, p1-d1) +# if i%10==0: +## +# p1 = f(operator.direct(x)) + g(x) +## print(p1) +# d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) ) +# primal.append(p1) +# dual.append(d1) +# pdgap.append(p1-d1) +# print(p1, d1, p1-d1) t_end = time.time() diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 40dddd7..7889cad 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -39,19 +39,14 @@ class KullbackLeibler(Function): def __call__(self, x): # TODO check - - self.sum_value = x + self.bnoise - if (self.sum_value.as_array()<0).any(): - self.sum_value = numpy.inf - - if self.sum_value==numpy.inf: - return numpy.inf - else: - tmp = self.sum_value.as_array() - return (x - self.b * ImageData( numpy.log(tmp))).sum() - -# return numpy.sum( x.as_array() - self.b.as_array() * numpy.log(self.sum_value.as_array())) + + tmp = x + self.bnoise + ind = tmp.as_array()>0 + res = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind]) + + return sum(res) + def gradient(self, x, out=None): @@ -64,10 +59,12 @@ class KullbackLeibler(Function): def convex_conjugate(self, x): - tmp = self.b.as_array()/( 1 - x.as_array() ) + tmp = self.b/( 1 - x ) + ind = tmp.as_array()>0 + + sel return (self.b * ( ImageData( numpy.log(tmp) ) - 1 ) - self.bnoise * (x - 1)).sum() -# return self.b * ( ImageData(numpy.log(self.b/(1-x)) - 1 )) - self.bnoise * (x - 1) def proximal(self, x, tau, out=None): @@ -105,19 +102,39 @@ class KullbackLeibler(Function): if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + import numpy N, M = 2,3 ig = ImageGeometry(N, M) - data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - - bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) + data = ImageData(numpy.random.randint(-10, 10, size=(M, N))) + x = ImageData(numpy.random.randint(-10, 10, size=(M, N))) - f = KullbackLeibler(data, bnoise=bnoise) - print(f.sum_value) + bnoise = ImageData(numpy.random.randint(-10, 10, size=(M, N))) - print(f(x)) + f = KullbackLeibler(data) + print(f(x)) + +# numpy.random.seed(10) +# +# +# x = numpy.random.randint(-10, 10, size = (2,3)) +# b = numpy.random.randint(1, 10, size = (2,3)) +# +# ind1 = x>0 +# +# res = x[ind1] - b * numpy.log(x[ind1]) +# +## ind = x>0 +# +## y = x[ind] +# +# +# +# +# \ No newline at end of file diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py index e6d38e9..0e7f628 100755 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising.py @@ -8,7 +8,8 @@ Created on Fri Feb 22 14:53:03 2019 from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer -import numpy as np +import numpy as np +import numpy import matplotlib.pyplot as plt from ccpi.optimisation.algorithms import PDHG, PDHG_old @@ -89,13 +90,12 @@ t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) t2 = timer() -print(" Run memopt") 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()) @@ -115,10 +115,67 @@ 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() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(primal[-1])) + + + + diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py index fbe0d9b..16e6b43 100644 --- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py +++ b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py @@ -6,7 +6,8 @@ Created on Fri Feb 22 14:53:03 2019 @author: evangelos """ -import numpy as np +import numpy as np +import numpy import matplotlib.pyplot as plt from ccpi.framework import ImageData, ImageGeometry @@ -46,6 +47,7 @@ plt.show() alpha = 2 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") + method = '0' if method == '0': @@ -57,7 +59,7 @@ if method == '0': operator = BlockOperator(op1, op2, shape=(2,1) ) f1 = alpha * MixedL21Norm() - f2 = KullbackLeibler(noisy_data) + f2 = 0.5 * KullbackLeibler(noisy_data) f = BlockFunction(f1, f2 ) g = ZeroFunction() @@ -69,9 +71,8 @@ else: ########################################################################### operator = Gradient(ig) f = alpha * MixedL21Norm() - g = KullbackLeibler(noisy_data) + g = 0.5 * KullbackLeibler(noisy_data) ########################################################################### -#%% # Compute operator Norm normK = operator.norm() @@ -87,13 +88,11 @@ t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) t2 = timer() -print(" Run memopt") - 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()) @@ -120,6 +119,59 @@ 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()) + + + 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() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(primal[-1])) + + + #pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) #pdhg.max_iteration = 2000 #pdhg.update_objective_interval = 10 -- cgit v1.2.3