From 0f542eafe04c4fe7568f83e935859665ffc77e3b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 5 Jul 2019 14:25:10 +0100 Subject: Norm2Sq and FISTA to give hints of why they fail (#351) * Norm2Sq and FISTA to give hints of why they fail * added denoising tests from demos * L1Norm store instance of ShrinkageOperator * relax limit for python2 * added source of tests --- .../Python/ccpi/optimisation/algorithms/FISTA.py | 3 +- .../Python/ccpi/optimisation/functions/L1Norm.py | 11 +- .../Python/ccpi/optimisation/functions/Norm2Sq.py | 6 +- Wrappers/Python/test/test_algorithms.py | 229 ++++++++++++++++++++- 4 files changed, 239 insertions(+), 10 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py index 347dacd..e23116b 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -73,7 +73,8 @@ class FISTA(Algorithm): self.f = f self.g = g - + if f.L is None: + raise ValueError('Error: Fidelity Function\'s Lipschitz constant is set to None') self.invL = 1/f.L self.t = 1 self.update_objective() diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py index 97d03b9..cc4bef8 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -35,7 +35,8 @@ class L1Norm(Function): def __init__(self, **kwargs): super(L1Norm, self).__init__() - self.b = kwargs.get('b',None) + self.b = kwargs.get('b',None) + self.shinkage_operator = ShrinkageOperator() def __call__(self, x): @@ -69,14 +70,14 @@ class L1Norm(Function): if out is None: if self.b is not None: - return self.b + ShrinkageOperator.__call__(self, x - self.b, tau) + return self.b + self.shinkage_operator(x - self.b, tau) else: - return ShrinkageOperator.__call__(self, x, tau) + return self.shinkage_operator(x, tau) else: if self.b is not None: - out.fill(self.b + ShrinkageOperator.__call__(self, x - self.b, tau)) + out.fill(self.b + self.shinkage_operator(x - self.b, tau)) else: - out.fill(ShrinkageOperator.__call__(self, x, tau)) + out.fill(self.shinkage_operator(x, tau)) def proximal_conjugate(self, x, tau, out=None): diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py index c8f5e46..0da6e50 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py +++ b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py @@ -50,9 +50,11 @@ class Norm2Sq(Function): try: self.L = 2.0*self.c*(self.A.norm()**2) except AttributeError as ae: - pass + warnings.warn('{} could not calculate Lipschitz Constant. {}'.format( + self.__class__.__name__, ae)) except NotImplementedError as noe: - pass + warnings.warn('{} could not calculate Lipschitz Constant. {}'.format( + self.__class__.__name__, noe)) #def grad(self,x): # return self.gradient(x, out=None) diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index f00467c..1577fa6 100755 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -30,8 +30,12 @@ from ccpi.optimisation.algorithms import GradientDescent from ccpi.optimisation.algorithms import CGLS from ccpi.optimisation.algorithms import FISTA +from ccpi.optimisation.algorithms import PDHG - +from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff +from ccpi.optimisation.functions import MixedL21Norm, BlockFunction, L1Norm, KullbackLeibler +from ccpi.framework import TestData +import os ,sys class TestAlgorithms(unittest.TestCase): def setUp(self): @@ -108,8 +112,229 @@ class TestAlgorithms(unittest.TestCase): alg.run(20, verbose=True) self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + def test_FISTA_Norm2Sq(self): + print ("Test FISTA Norm2Sq") + ig = ImageGeometry(127,139,149) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + x_init = ig.allocate(ImageGeometry.RANDOM) + identity = Identity(ig) + + #### it seems FISTA does not work with Nowm2Sq + norm2sq = Norm2Sq(identity, b) + #norm2sq.L = 2 * norm2sq.c * identity.norm()**2 + #norm2sq = FunctionOperatorComposition(L2NormSquared(b=b), identity) + opt = {'tol': 1e-4, 'memopt':False} + print ("initial objective", norm2sq(x_init)) + alg = FISTA(x_init=x_init, f=norm2sq, g=ZeroFunction()) + alg.max_iteration = 2 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + def test_FISTA_catch_Lipschitz(self): + print ("Test FISTA catch Lipschitz") + ig = ImageGeometry(127,139,149) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + x_init = ig.allocate(ImageGeometry.RANDOM) + identity = Identity(ig) + + #### it seems FISTA does not work with Nowm2Sq + norm2sq = Norm2Sq(identity, b) + print ('Lipschitz', norm2sq.L) + norm2sq.L = None + #norm2sq.L = 2 * norm2sq.c * identity.norm()**2 + #norm2sq = FunctionOperatorComposition(L2NormSquared(b=b), identity) + opt = {'tol': 1e-4, 'memopt':False} + print ("initial objective", norm2sq(x_init)) + try: + alg = FISTA(x_init=x_init, f=norm2sq, g=ZeroFunction()) + self.assertTrue(False) + except ValueError as ve: + print (ve) + self.assertTrue(True) + def test_PDHG_Denoising(self): + print ("PDHG Denoising with 3 noises") + # adapted from demo PDHG_TV_Color_Denoising.py in CIL-Demos repository + + loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) + data = loader.load(TestData.PEPPERS, size=(256,256)) + ig = data.geometry + ag = ig + + which_noise = 0 + # Create noisy data. + noises = ['gaussian', 'poisson', 's&p'] + noise = noises[which_noise] + + def setup(data, noise): + if noise == 's&p': + n1 = TestData.random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2, seed=10) + elif noise == 'poisson': + scale = 5 + n1 = TestData.random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale + elif noise == 'gaussian': + n1 = TestData.random_noise(data.as_array(), mode = noise, seed = 10) + else: + raise ValueError('Unsupported Noise ', noise) + noisy_data = ig.allocate() + noisy_data.fill(n1) + + # Regularisation Parameter depending on the noise distribution + if noise == 's&p': + alpha = 0.8 + elif noise == 'poisson': + alpha = 1 + elif noise == 'gaussian': + alpha = .3 + # fidelity + if noise == 's&p': + g = L1Norm(b=noisy_data) + elif noise == 'poisson': + g = KullbackLeibler(noisy_data) + elif noise == 'gaussian': + g = 0.5 * L2NormSquared(b=noisy_data) + return noisy_data, alpha, g + + noisy_data, alpha, g = setup(data, noise) + operator = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) + + f1 = alpha * MixedL21Norm() + + + + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + + # Setup and run the PDHG algorithm + pdhg1 = PDHG(f=f1,g=g,operator=operator, tau=tau, sigma=sigma) + pdhg1.max_iteration = 2000 + pdhg1.update_objective_interval = 200 + pdhg1.run(1000) + + rmse = (pdhg1.get_output() - data).norm() / data.as_array().size + print ("RMSE", rmse) + self.assertLess(rmse, 2e-4) + + which_noise = 1 + noise = noises[which_noise] + noisy_data, alpha, g = setup(data, noise) + operator = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) + + f1 = alpha * MixedL21Norm() + + + + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + + # Setup and run the PDHG algorithm + pdhg1 = PDHG(f=f1,g=g,operator=operator, tau=tau, sigma=sigma) + pdhg1.max_iteration = 2000 + pdhg1.update_objective_interval = 200 + pdhg1.run(1000) + + rmse = (pdhg1.get_output() - data).norm() / data.as_array().size + print ("RMSE", rmse) + self.assertLess(rmse, 2e-4) + + + which_noise = 2 + noise = noises[which_noise] + noisy_data, alpha, g = setup(data, noise) + operator = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) + + f1 = alpha * MixedL21Norm() + + + + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + + # Setup and run the PDHG algorithm + pdhg1 = PDHG(f=f1,g=g,operator=operator, tau=tau, sigma=sigma) + pdhg1.max_iteration = 2000 + pdhg1.update_objective_interval = 200 + pdhg1.run(1000) + + rmse = (pdhg1.get_output() - data).norm() / data.as_array().size + print ("RMSE", rmse) + self.assertLess(rmse, 2e-4) + + def test_FISTA_Denoising(self): + print ("FISTA Denoising Poisson Noise Tikhonov") + # adapted from demo FISTA_Tikhonov_Poisson_Denoising.py in CIL-Demos repository + loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) + data = loader.load(TestData.SHAPES) + ig = data.geometry + ag = ig + N=300 + # Create Noisy data with Poisson noise + scale = 5 + n1 = TestData.random_noise( data.as_array()/scale, mode = 'poisson', seed = 10)*scale + noisy_data = ImageData(n1) + + # Regularisation Parameter + alpha = 10 + + # Setup and run the FISTA algorithm + operator = Gradient(ig) + fid = KullbackLeibler(noisy_data) + + def KL_Prox_PosCone(x, tau, out=None): + + if out is None: + tmp = 0.5 *( (x - fid.bnoise - tau) + ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b ) .sqrt() ) + return tmp.maximum(0) + else: + tmp = 0.5 *( (x - fid.bnoise - tau) + + ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b ) .sqrt() + ) + x.add(fid.bnoise, out=out) + out -= tau + out *= out + tmp = fid.b * (4 * tau) + out.add(tmp, out=out) + out.sqrt(out=out) + + x.subtract(fid.bnoise, out=tmp) + tmp -= tau + + out += tmp + + out *= 0.5 + + # ADD the constraint here + out.maximum(0, out=out) + + fid.proximal = KL_Prox_PosCone + + reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator) + + x_init = ig.allocate() + fista = FISTA(x_init=x_init , f=reg, g=fid) + fista.max_iteration = 3000 + fista.update_objective_interval = 500 + fista.run(3000, verbose=True) + rmse = (fista.get_output() - data).norm() / data.as_array().size + print ("RMSE", rmse) + self.assertLess(rmse, 4.2e-4) - def assertNumpyArrayEqual(self, first, second): res = True try: -- cgit v1.2.3