diff options
8 files changed, 67 insertions, 429 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py index bc9b62d..d6c98c4 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -7,7 +7,8 @@ Created on Fri Mar 8 10:01:31 2019 """ import numpy as np -from ccpi.optimisation.funcs import Function +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function from ccpi.framework import BlockDataContainer class BlockFunction(Function): diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py index a67b884..0f3defe 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -7,7 +7,8 @@ Created on Fri Mar 8 09:55:36 2019 """ import numpy as np -from ccpi.optimisation.funcs import Function +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function class FunctionOperatorComposition(Function): diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py index ec7aa5b..f83de6f 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -7,12 +7,13 @@ Created on Wed Mar 6 19:42:34 2019 """ import numpy as np -from ccpi.optimisation.funcs import Function +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function from ccpi.framework import DataContainer, ImageData, ImageGeometry ############################ L1NORM FUNCTIONS ############################# -class SimpleL1NormEdo(Function): +class SimpleL1Norm(Function): def __init__(self, alpha=1): @@ -35,7 +36,7 @@ class SimpleL1NormEdo(Function): def proximal_conjugate(self, x, tau): return x.divide((x.abs()/self.alpha).maximum(1.0)) -class L1Norm(SimpleL1NormEdo): +class L1Norm(SimpleL1Norm): def __init__(self, alpha=1, **kwargs): diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 8b73be6..5817317 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -9,7 +9,8 @@ Created on Thu Feb 7 13:10:56 2019 """ import numpy as np -from ccpi.optimisation.funcs import Function +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function from ccpi.framework import DataContainer, ImageData, ImageGeometry diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py index b41cc26..9def741 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py @@ -7,7 +7,8 @@ Created on Wed Mar 6 19:44:10 2019 """ import numpy as np -from ccpi.optimisation.funcs import Function +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function from ccpi.framework import DataContainer, ImageData from ccpi.framework import BlockDataContainer @@ -33,8 +34,11 @@ class ZeroFun(Function): else: return x.maximum(0).sum() + x.maximum(0).sum() - def proximal(self,x,tau): - return x.copy() + def proximal(self,x,tau, out=None): + if out is None: + return x.copy() + else: + out.fill(x) def proximal_conjugate(self, x, tau): return 0
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/__init__.py b/Wrappers/Python/ccpi/optimisation/functions/__init__.py index 7ce617a..c4ba0a6 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- - +from .Function import Function from .ZeroFun import ZeroFun from .L1Norm import * from .L2NormSquared import * diff --git a/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py index 0ce1d28..8fe8620 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py @@ -7,9 +7,10 @@ Created on Wed Mar 6 19:43:12 2019 """ import numpy as np -from ccpi.optimisation.funcs import Function +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function from ccpi.framework import DataContainer, ImageData, ImageGeometry - +from ccpi.optimisation.functions.FunctionOperatorComposition import FunctionOperatorComposition ############################ mixed_L1,2NORM FUNCTIONS ############################# @@ -59,7 +60,7 @@ class mixed_L12Norm(Function): def composition_with(self, operator): if self.b is None: - return FunctionComposition(mixed_L12Norm(self.alpha), operator) + return FunctionOperatorComposition(operator, mixed_L12Norm(self.alpha)) else: - return FunctionComposition(mixed_L12Norm(self.alpha, b=self.b), operator) + return FunctionOperatorComposition(operator, mixed_L12Norm(self.alpha, b=self.b)) diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index 54a952a..0741d1c 100644 --- a/Wrappers/Python/test/test_functions.py +++ b/Wrappers/Python/test/test_functions.py @@ -6,438 +6,67 @@ Created on Sat Mar 2 19:24:37 2019 @author: evangelos """ -# -*- coding: utf-8 -*- - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Feb 7 13:10:56 2019 - -@author: evangelos -""" import numpy as np -from ccpi.optimisation.funcs import Function +#from ccpi.optimisation.funcs import Function +from ccpi.optimisation.functions import Function from ccpi.framework import DataContainer, ImageData, ImageGeometry -from operators import CompositeDataContainer, Identity, CompositeOperator +from ccpi.optimisation.operators import Identity +from ccpi.optimisation.operators import BlockOperator +from ccpi.framework import BlockDataContainer from numbers import Number -from GradientOperator import Gradient - - -class SimpleL2NormSq(Function): - - def __init__(self, alpha=1): - - super(SimpleL2NormSq, self).__init__() - self.alpha = alpha - - def __call__(self, x): - return self.alpha * x.power(2).sum() - - def gradient(self,x): - return 2 * self.alpha * x - - def convex_conjugate(self,x): - return (1/(4*self.alpha)) * x.power(2).sum() - - def proximal(self, x, tau): - return x.divide(1+2*tau*self.alpha) - - def proximal_conjugate(self, x, tau): - return x.divide(1 + tau/(2*self.alpha) ) - - -############################ L2NORM FUNCTIONS ############################# -class L2NormSq(SimpleL2NormSq): - - def __init__(self, alpha, **kwargs): - - super(L2NormSq, self).__init__(alpha) - self.alpha = alpha - self.b = kwargs.get('b',None) - self.L = 1 - - def __call__(self, x): - - if self.b is None: - return SimpleL2NormSq.__call__(self, x) - else: - return SimpleL2NormSq.__call__(self, x - self.b) - - def gradient(self, x): - - if self.b is None: - return 2*self.alpha * x - else: - return 2*self.alpha * (x - self.b) - - def composition_with(self, operator): - - if self.b is None: - return FunctionComposition(L2NormSq(self.alpha), operator) - else: - return FunctionComposition(L2NormSq(self.alpha, b=self.b), operator) - - def convex_conjugate(self, x): - - ''' The convex conjugate corresponds to the simple functional i.e., - f(x) = alpha * ||x - b||_{2}^{2} - ''' - - if self.b is None: - return SimpleL2NormSq.convex_conjugate(self, x) - else: - return SimpleL2NormSq.convex_conjugate(self, x) + (self.b * x).sum() - - def proximal(self, x, tau): - - ''' The proximal operator corresponds to the simple functional i.e., - f(x) = alpha * ||x - b||_{2}^{2} - - argmin_x { 0.5||x - u||^{2} + tau f(x) } - ''' - - if self.b is None: - return SimpleL2NormSq.proximal(self, x, tau) - else: - return self.b + SimpleL2NormSq.proximal(self, x - self.b , tau) - - - def proximal_conjugate(self, x, tau): - - ''' The proximal operator corresponds to the simple convex conjugate - functional i.e., f^{*}(x^{) - argmin_x { 0.5||x - u||^{2} + tau f(x) } - ''' - if self.b is None: - return SimpleL2NormSq.proximal_conjugate(self, x, tau) - else: - return SimpleL2NormSq.proximal_conjugate(self, x - tau * self.b, tau) - - -############################ L1NORM FUNCTIONS ############################# -class SimpleL1Norm(Function): - - def __init__(self, alpha=1): - - super(SimpleL1Norm, self).__init__() - self.alpha = alpha - - def __call__(self, x): - return self.alpha * x.abs().sum() - - def gradient(self,x): - return ValueError('Not Differentiable') - - def convex_conjugate(self,x): - return 0 - - def proximal(self, x, tau): - ''' Soft Threshold''' - return x.sign() * (x.abs() - tau * self.alpha).maximum(1.0) - - def proximal_conjugate(self, x, tau): - return x.divide((x.abs()/self.alpha).maximum(1.0)) - -class L1Norm(SimpleL1Norm): - - def __init__(self, alpha=1, **kwargs): - - super(L1Norm, self).__init__() - self.alpha = alpha - - self.A = kwargs.get('A',None) - self.b = kwargs.get('b',None) - - def __call__(self, x): - - if self.b is None: - return SimpleL1Norm.__call__(self, self.A.direct(x)) - else: - return SimpleL1Norm.__call__(self, self.A.direct(x) - self.b) - - def eval_norm(self, x): - - return SimpleL1Norm.__call__(self, x) - - def gradient(self, x): - return ValueError('Not Differentiable') - - def convex_conjugate(self,x): - if self.b is None: - return SimpleL1Norm.convex_conjugate(self, x) - else: - return SimpleL1Norm.convex_conjugate(self, x) + (self.b * x).sum() - - def proximal(self, x, tau): - - if self.b is None: - return SimpleL1Norm.proximal(self, x, tau) - else: - return self.b + SimpleL1Norm.proximal(self, x + self.b , tau) - - def proximal_conjugate(self, x, tau): - - if self.b is None: - return SimpleL1Norm.proximal_conjugate(self, x, tau) - else: - return SimpleL1Norm.proximal_conjugate(self, x - tau*self.b, tau) - - -############################ mixed_L1,2NORM FUNCTIONS ############################# -class mixed_L12Norm(Function): - - def __init__(self, alpha, **kwargs): +from ccpi.optimisation.operators import Gradient - super(mixed_L12Norm, self).__init__() - - self.alpha = alpha - self.b = kwargs.get('b',None) - self.sym_grad = kwargs.get('sym_grad',False) - - def __call__(self,x): - - if self.b is None: - tmp1 = x - else: - tmp1 = x - self.b -# - if self.sym_grad: - tmp = np.sqrt(tmp1.as_array()[0]**2 + tmp1.as_array()[1]**2 + 2*tmp1.as_array()[2]**2) - else: - tmp = ImageData(tmp1.power(2).sum(axis=0)).sqrt() - - return self.alpha*tmp.sum() - - def gradient(self,x): - return ValueError('Not Differentiable') - - def convex_conjugate(self,x): - return 0 - - def proximal(self, x, tau): - pass - - def proximal_conjugate(self, x, tau): - - if self.sym_grad: - tmp2 = np.sqrt(x.as_array()[0]**2 + x.as_array()[1]**2 + 2*x.as_array()[2]**2)/self.alpha - res = x.divide(ImageData(tmp2).maximum(1.0)) - else: - res = x.divide((ImageData(x.power(2).sum(axis=0)).sqrt()/self.alpha).maximum(1.0)) - - return res - - def composition_with(self, operator): - - if self.b is None: - return FunctionComposition(mixed_L12Norm(self.alpha), operator) - else: - return FunctionComposition(mixed_L12Norm(self.alpha, b=self.b), operator) - +from ccpi.optimisation.functions import SimpleL2NormSq +from ccpi.optimisation.functions import L2NormSq +from ccpi.optimisation.functions import SimpleL1Norm +from ccpi.optimisation.functions import L1Norm +# from ccpi.optimisation.functions.L2NormSquared import SimpleL2NormSq, L2NormSq +# from ccpi.optimisation.functions.L1Norm import SimpleL1Norm, L1Norm +from ccpi.optimisation.functions import mixed_L12Norm +from ccpi.optimisation.functions import ZeroFun -#%% - -class ZeroFun(Function): - - def __init__(self): - super(ZeroFun, self).__init__() - - def __call__(self,x): - return 0 - - def convex_conjugate(self, x): - ''' This is the support function sup <x, x^{*}> which in fact is the - indicator function for the set = {0} - So 0 if x=0, or inf if x neq 0 - ''' - - if x.shape[0]==1: - return x.maximum(0).sum() - else: - return x.get_item(0).maximum(0).sum() + x.get_item(1).maximum(0).sum() - - def proximal(self,x,tau): - return x.copy() - - def proximal_conjugate(self, x, tau): - return 0 - - -class CompositeFunction(Function): - - def __init__(self, *functions, blockMatrix): - - self.blockMatrix = blockMatrix - self.functions = functions - self.length = len(self.functions) - - def get_item(self, ind): - return self.functions[ind] - - def __call__(self,x): - - t = 0 - for i in range(x.shape[0]): - t += self.functions[i](x.get_item(i)) - return t - - - def convex_conjugate(self, x): - - z = 0 - t = 0 - for i in range(x.shape[0]): - t += self.functions[z].convex_conjugate(x.get_item(i)) - z += 1 - - return t - - def proximal(self, x, tau, out = None): - - if isinstance(tau, Number): - tau = CompositeDataContainer(tau) - out = [None]*self.length - for i in range(self.length): - out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(0)) - return CompositeDataContainer(*out) +from ccpi.optimisation.functions import FunctionOperatorComposition +import unittest - - def proximal_conjugate(self, x, tau, out = None): - - if isinstance(tau, Number): - tau = CompositeDataContainer(tau) - if isinstance(x, ImageData): - x = CompositeDataContainer(x) - out = [None]*self.length - for i in range(self.length): - out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(0)) - return CompositeDataContainer(*out) - # -class FunctionComposition(Function): - - def __init__(self, function, operator): - - self.function = function - self.alpha = self.function.alpha - self.b = self.function.b - self.operator = operator - - - super(FunctionComposition, self).__init__() - - ''' overide call and gradient ''' - def __call__(self, x): - return self.function(x) - - def gradient(self,x): - return self.operator.adjoint(self.function.gradient(self.operator.direct(x))) - - ''' Same as in the parent class''' - def proximal(self,x, tau): - return self.function.proximal(x, tau) - - def proximal_conjugate(self,x, tau): - return self.function.proximal_conjugate(x, tau) - def convex_conjugate(self,x): - return self.function.convex_conjugate(x) - -class FunctionComposition_new(Function): +class TestFunction(unittest.TestCase): + def test_Function(self): + - def __init__(self, operator, *functions): + N = 3 + ig = (N,N) + ag = ig + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Form Composite Operator + operator = BlockOperator((2,1), op1, op2 ) - self.functions = functions - self.operator = operator - self.length = len(self.functions) + # Create functions + noisy_data = ImageData(np.random.randint(10, size=ag)) -# if self.length==1: -# self.L = self.functions[0].alpha*(self.operator.norm()**2) - - # length == to operator.shape[0]# - super(FunctionComposition_new, self).__init__() - - def __call__(self, x): - - if isinstance(x, ImageData): - x = CompositeDataContainer(x) + d = ImageData(np.random.randint(10, size=ag)) - t = 0 - for i in range(x.shape[0]): - t += self.functions[i](x.get_item(i)) - return t - - def convex_conjugate(self, x): - - if isinstance(x, ImageData): - x = CompositeDataContainer(x) + f = mixed_L12Norm(alpha = 1).composition_with(op1) + g = L2NormSq(alpha=0.5, b=noisy_data) - t = 0 - for i in range(x.shape[0]): - t += self.functions[i].convex_conjugate(x.get_item(i)) - return t - - def proximal_conjugate(self, x, tau, out = None): - - if isinstance(x, ImageData): - x = CompositeDataContainer(x) - - out = [None]*self.length - for i in range(self.length): - out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau) + # Compare call of f + a1 = ImageData(op1.direct(d).power(2).sum(axis=0)).sqrt().sum() + #print(a1, f(d)) + self.assertEqual (a1, f(d)) - if self.length==1: - return ImageData(*out) - else: - return CompositeDataContainer(*out) - - def proximal(self, x, tau, out = None): - - if isinstance(x, ImageData): - x = CompositeDataContainer(x) - - out = [None]*self.length - for i in range(self.length): - out[i] = self.functions[i].proximal(x.get_item(i), tau) + # Compare call of g + a2 = g.alpha*(d - noisy_data).power(2).sum() + #print(a2, g(d)) + self.assertEqual(a2, g(d)) - if self.length==1: - return ImageData(*out) - else: - return CompositeDataContainer(*out) - - -if __name__ == '__main__': - - N = 3 - ig = (N,N) - ag = ig - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Form Composite Operator - operator = CompositeOperator((2,1), op1, op2 ) - - # Create functions - noisy_data = ImageData(np.random.randint(10, size=ag)) - - d = ImageData(np.random.randint(10, size=ag)) - - f = mixed_L12Norm(alpha = 1).composition_with(op1) - g = L2NormSq(alpha=0.5, b=noisy_data) - - # Compare call of f - a1 = ImageData(op1.direct(d).power(2).sum(axis=0)).sqrt().sum() - print(a1, f(d)) - - # Compare call of g - a2 = g.alpha*(d - noisy_data).power(2).sum() - print(a2, g(d)) - - # Compare convex conjugate of g - a3 = 0.5 * d.power(2).sum() + (d*noisy_data).sum() - print( a3, g.convex_conjugate(d)) + # Compare convex conjugate of g + a3 = 0.5 * d.power(2).sum() + (d*noisy_data).sum() + self.assertEqual(a3, g.convex_conjugate(d)) + #print( a3, g.convex_conjugate(d)) |