diff options
26 files changed, 2123 insertions, 1 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py new file mode 100644 index 0000000..7488310 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Feb 4 16:18:06 2019 + +@author: evangelos +""" + +from ccpi.framework import ImageData +import numpy as np +import matplotlib.pyplot as plt +import time +from Operators.CompositeOperator import CompositeOperator +from Operators.CompositeDataContainer import CompositeDataContainer + +def PDHG(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-6, 'niter': 500, 'show_iter': 100, \ + 'memopt': False} + + if sigma is None and tau is None: + raise ValueError('Need sigma*tau||K||^2<1') + + niter = opt['niter'] if 'niter' in opt.keys() else 1000 + tol = opt['tol'] if 'tol' in opt.keys() else 1e-4 + memopt = opt['memopt'] if 'memopt' in opt.keys() else False + show_iter = opt['show_iter'] if 'show_iter' in opt.keys() else False + stop_crit = opt['stop_crit'] if 'stop_crit' in opt.keys() else False + + if isinstance(operator, CompositeOperator): +# if isinstance(operator, CompositeOperator_DataContainer): + x_old = operator.alloc_domain_dim() + y_old = operator.alloc_range_dim() + else: + x_old = ImageData(np.zeros(operator.domain_dim())) + y_old = ImageData(np.zeros(operator.range_dim())) + + + xbar = x_old + x_tmp = x_old + x = x_old + + y_tmp = y_old + y = y_tmp + + # relaxation parameter + theta = 1 + + t = time.time() + + objective = [] + + for i in range(niter): + + # Gradient descent, Dual problem solution + y_tmp = y_old + sigma * operator.direct(xbar) + y = f.proximal_conjugate(y_tmp, sigma) + + # Gradient ascent, Primal problem solution + x_tmp = x_old - tau * operator.adjoint(y) + x = g.proximal(x_tmp, tau) + + #Update + xbar = x + theta * (x - x_old) + + x_old = x + y_old = y + +# pdgap + print(f(x) + g(x) + f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) ) + + + + + +# # TV denoising, pdgap with composite +# +# primal_obj = f.get_item(0).alpha * ImageData(operator.compMat[0][0].direct(x.get_item(0)).power(2).sum(axis=0)).sqrt().sum() +\ +# 0.5*( (operator.compMat[1][0].direct(x.get_item(0)) - f.get_item(1).b).power(2).sum()) +# dual_obj = 0.5 * ((y.get_item(1).power(2)).sum()) + ( y.get_item(1)*f.get_item(1).b ).sum() + + # TV denoising, pdgap with no composite + + + +# primal_obj = f.get_item(0).alpha * ImageData(operator.compMat[0][0].direct(x.get_item(0)).power(2).sum(axis=0)).sqrt().sum() +\ +# 0.5*( (operator.compMat[1][0].direct(x.get_item(0)) - f.get_item(1).b).power(2).sum()) +# dual_obj = 0.5 * ((y.get_item(1).power(2)).sum()) + ( y.get_item(1)*f.get_item(1).b ).sum() + + +# print(primal_obj) +# objective = primal_obj +# + + + + t_end = time.time() + + return x, t_end - t, objective + diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py new file mode 100644 index 0000000..bc9b62d --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -0,0 +1,69 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 10:01:31 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.optimisation.funcs import Function +from ccpi.framework import BlockDataContainer + +class BlockFunction(Function): + + def __init__(self, operator, *functions): + + self.functions = functions + self.operator = operator + self.length = len(self.functions) + + super(BlockFunction, self).__init__() + + def __call__(self, x): + + tmp = self.operator.direct(x) + + t = 0 + for i in range(tmp.shape[0]): + t += self.functions[i](tmp.get_item(i)) + return t + + def call_adjoint(self, x): + + tmp = operator.adjoint(x) + + t = 0 + for i in range(tmp.shape[0]): + t += self.functions[i](tmp.get_item(i)) + return t + + def convex_conjugate(self, x): + + ''' Convex_conjugate does not take into account the BlockOperator''' + 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): + + ''' proximal_conjugate does not take into account the BlockOperator''' + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau) + + return BlockDataContainer(*out) + + def proximal(self, x, tau, out = None): + + ''' proximal does not take into account the BlockOperator''' + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau) + + return BlockDataContainer(*out) + + def gradient(self,x, out=None): + pass
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionComposition.py new file mode 100644 index 0000000..5b6defc --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionComposition.py @@ -0,0 +1,175 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:45:06 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.optimisation.funcs import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry +from ccpi.framework import BlockDataContainer + +class FunctionOperatorComposition(Function): + + def __init__(self, function, operator): + + self.function = functions + self.operator = operator + self.grad_Lipschitz_cnst = 2*self.function.alpha*operator.norm()**2 + super(FunctionOperatorComposition, self).__init__() + + def __call__(self, x): + + return self.function(operator.direct(x)) + + def call_adjoint(self, x): + + return self.function(operator.adjoint(x)) + + def convex_conjugate(self, x): + + return self.function.convex_conjugate(x) + + def proximal(self, x, tau): + + ''' proximal does not take into account the Operator''' + + return self.function.proximal(x, tau, out=None) + + def proximal_conjugate(self, x, tau): + + ''' proximal conjugate does not take into account the Operator''' + + return self.function.proximal_conjugate(x, tau, out=None) + + def gradient(self, x): + + ''' Gradient takes into account the Operator''' + + return self.adjoint(self.function.gradient(self.operator.direct(x))) + + +class BlockFunction(Function): + + def __init__(self, operator, *functions): + + self.functions = functions + self.operator = operator + self.length = len(self.functions) + + super(BlockFunction, self).__init__() + + def __call__(self, x): + + tmp = operator.direct(x) + + t = 0 + for i in range(tmp.shape[0]): + t += self.functions[i](tmp.get_item(i)) + return t + + def call_adjoint(self, x): + + tmp = operator.adjoint(x) + + t = 0 + for i in range(tmp.shape[0]): + t += self.functions[i](tmp.get_item(i)) + return t + + def convex_conjugate(self, x): + + ''' Convex_conjugate does not take into account the BlockOperator''' + 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): + + ''' proximal_conjugate does not take into account the BlockOperator''' + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau) + + return CompositeDataContainer(*out) + + def proximal(self, x, tau, out = None): + + ''' proximal does not take into account the BlockOperator''' + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau) + + return CompositeDataContainer(*out) + + def gradient(self,x, out=None): + pass + + +class FunctionComposition_new(Function): + + def __init__(self, operator, *functions): + + self.functions = functions + self.operator = operator + self.length = len(self.functions) + + super(FunctionComposition_new, self).__init__() + + def __call__(self, x): + + if isinstance(x, ImageData): + x = CompositeDataContainer(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): + + if isinstance(x, ImageData): + x = CompositeDataContainer(x) + + 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) + + 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) + + if self.length==1: + return ImageData(*out) + else: + return CompositeDataContainer(*out) + + +if __name__ == '__main__': + + from operators import Operator + from IdentityOperator import Identity
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py new file mode 100644 index 0000000..a67b884 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 09:55:36 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.optimisation.funcs import Function + + +class FunctionOperatorComposition(Function): + + def __init__(self, operator, function): + + self.function = function + self.operator = operator + self.L = 2*self.function.alpha*operator.norm()**2 + super(FunctionOperatorComposition, self).__init__() + + def __call__(self, x): + + return self.function(self.operator.direct(x)) + + def call_adjoint(self, x): + + return self.function(self.operator.adjoint(x)) + + def convex_conjugate(self, x): + + ''' convex_conjugate does not take into account the Operator''' + return self.function.convex_conjugate(x) + + def proximal(self, x, tau): + + ''' proximal does not take into account the Operator''' + + return self.function.proximal(x, tau, out=None) + + def proximal_conjugate(self, x, tau, out=None): + + ''' proximal conjugate does not take into account the Operator''' + + return self.function.proximal_conjugate(x, tau) + + def gradient(self, x): + + ''' Gradient takes into account the Operator''' + + return self.operator.adjoint(self.function.gradient(self.operator.direct(x))) + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py new file mode 100644 index 0000000..ec7aa5b --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:42:34 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.optimisation.funcs import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry + + +############################ L1NORM FUNCTIONS ############################# +class SimpleL1NormEdo(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(0) + + def proximal_conjugate(self, x, tau): + return x.divide((x.abs()/self.alpha).maximum(1.0)) + +class L1Norm(SimpleL1NormEdo): + + def __init__(self, alpha=1, **kwargs): + + super(L1Norm, self).__init__() + self.alpha = alpha + self.b = kwargs.get('b',None) + + def __call__(self, x): + + if self.b is None: + return SimpleL1Norm.__call__(self, x) + else: + return SimpleL1Norm.__call__(self, x - self.b) + + 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) +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py new file mode 100644 index 0000000..8b73be6 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -0,0 +1,101 @@ +# -*- 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.framework import DataContainer, ImageData, ImageGeometry + + +class SimpleL2NormSq(Function): + + def __init__(self, alpha=1): + + super(SimpleL2NormSq, self).__init__() + self.alpha = alpha + + # Lispchitz constant of gradient + self.L = 2*self.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) + + 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 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) + + diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py new file mode 100644 index 0000000..b41cc26 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:44:10 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.optimisation.funcs import Function +from ccpi.framework import DataContainer, ImageData +from ccpi.framework import BlockDataContainer + +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: + if isinstance(x, CompositeDataContainer): + return x.get_item(0).maximum(0).sum() + x.get_item(1).maximum(0).sum() + else: + return x.maximum(0).sum() + x.maximum(0).sum() + + def proximal(self,x,tau): + return x.copy() + + 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 new file mode 100644 index 0000000..7ce617a --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py @@ -0,0 +1,10 @@ +# -*- coding: utf-8 -*- + + +from .ZeroFun import ZeroFun +from .L1Norm import * +from .L2NormSquared import * +from .mixed_L12Norm import * +from .FunctionOperatorComposition import FunctionOperatorComposition +from .BlockFunction import BlockFunction + diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/BlockFunction.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/BlockFunction.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..660532e --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/BlockFunction.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/FunctionComposition.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/FunctionComposition.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..075cdfb --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/FunctionComposition.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/FunctionOperatorComposition.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/FunctionOperatorComposition.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..f564eff --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/FunctionOperatorComposition.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/L1Norm.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/L1Norm.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..4ef959d --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/L1Norm.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/L2NormSquared.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/L2NormSquared.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..ad1b296 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/L2NormSquared.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/ZeroFun.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/ZeroFun.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..f2bf70f --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/ZeroFun.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/__init__.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/__init__.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..1321257 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/__init__.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/__pycache__/mixed_L12Norm.cpython-36.pyc b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/mixed_L12Norm.cpython-36.pyc Binary files differnew file mode 100644 index 0000000..d43e3ad --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/__pycache__/mixed_L12Norm.cpython-36.pyc diff --git a/Wrappers/Python/ccpi/optimisation/functions/functions.py b/Wrappers/Python/ccpi/optimisation/functions/functions.py new file mode 100644 index 0000000..f40abb9 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/functions.py @@ -0,0 +1,311 @@ +# -*- 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.framework import DataContainer, ImageData, ImageGeometry +from operators import CompositeDataContainer, Identity, CompositeOperator +from numbers import Number + + +############################ L2NORM FUNCTIONS ############################# +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 ) + + +class L2NormSq(SimpleL2NormSq): + + def __init__(self, A, b = None, alpha=1, **kwargs): + + super(L2NormSq, self).__init__(alpha=alpha) + self.alpha = alpha + self.A = A + self.b = b + + def __call__(self, x): + + if self.b is None: + return SimpleL2NormSq.__call__(self, self.A.direct(x)) + else: + return SimpleL2NormSq.__call__(self, self.A.direct(x) - self.b) + + 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 gradient(self, x): + + if self.b is None: + return 2*self.alpha * self.A.adjoint(self.A.direct(x)) + else: + return 2*self.alpha * self.A.adjoint(self.A.direct(x) - self.b) + + 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, A, b = None, alpha=1, **kwargs): + + super(L1Norm, self).__init__() + self.alpha = alpha + self.A = A + self.b = b + + 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 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, A, b=None, alpha=1, **kwargs): + + super(mixed_L12Norm, self).__init__() + self.alpha = alpha + self.A = A + self.b = b + + self.sym_grad = kwargs.get('sym_grad',False) + + + + def gradient(self,x): + return ValueError('Not Differentiable') + + + def __call__(self,x): + + y = self.A.direct(x) + eucl_norm = ImageData(y.power(2).sum(axis=0)).sqrt() + eucl_norm.__isub__(self.b) + return eucl_norm.sum() * self.alpha + + def convex_conjugate(self,x): + return 0 + + def proximal_conjugate(self, x, tau): + + if self.b is None: + + 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)) + + else: + res = (x - tau*self.b)/ ((x - tau*self.b)).abs().maximum(1.0) + + return res + + +#%% + +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 + ''' + return x.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, *args): + self.functions = args + self.length = len(self.functions) + + def get_item(self, ind): + return self.functions[ind] + + def __call__(self,x): + + t = 0 + for i in range(self.length): + for j in range(x.shape[0]): + t +=self.functions[i](x.get_item(j)) + 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_conjugate(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) + + + def proximal_conjugate(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_conjugate(x.get_item(i), tau.get_item(0)) + 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 + alpha = 1 + noisy_data = ImageData(np.random.randint(10, size=ag)) + f = CompositeFunction(L1Norm(op1,alpha), \ + L2NormSq(op2, noisy_data, c = 0.5, memopt = False) ) + + u = ImageData(np.random.randint(10, size=ig)) + uComp = CompositeDataContainer(u) + + print(f(uComp)) # This is f(Kx) = f1(K1*u) + f2(K2*u) + + f1 = L1Norm(op1,alpha) + f2 = L2NormSq(op2, noisy_data, c = 0.5, memopt = False) + + print(f1(u) + f2(u)) + + + diff --git a/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py new file mode 100644 index 0000000..0ce1d28 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:43:12 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.optimisation.funcs import Function +from ccpi.framework import DataContainer, ImageData, ImageGeometry + + + +############################ mixed_L1,2NORM FUNCTIONS ############################# +class mixed_L12Norm(Function): + + def __init__(self, alpha, **kwargs): + + 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) + diff --git a/Wrappers/Python/ccpi/optimisation/functions/test_functions.py b/Wrappers/Python/ccpi/optimisation/functions/test_functions.py new file mode 100644 index 0000000..54a952a --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/test_functions.py @@ -0,0 +1,474 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +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.framework import DataContainer, ImageData, ImageGeometry +from operators import CompositeDataContainer, Identity, CompositeOperator +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): + + 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) + + +#%% + +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) + + + 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): + + def __init__(self, operator, *functions): + + self.functions = functions + self.operator = operator + self.length = len(self.functions) + +# 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) + + 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) + + 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) + + 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) + + 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)) + + + + + +# +# f1 = L2NormSq(alpha=1, b=noisy_data) +# print(f1(noisy_data)) +# +# f2 = L2NormSq(alpha=5, b=noisy_data).composition_with(op2) +# print(f2(noisy_data)) +# +# print(f1.gradient(noisy_data).as_array()) +# print(f2.gradient(noisy_data).as_array()) +## +# print(f1.proximal(noisy_data,1).as_array()) +# print(f2.proximal(noisy_data,1).as_array()) +# +# +# f3 = mixed_L12Norm(alpha = 1).composition_with(op1) +# print(f3(noisy_data)) +# +# print(ImageData(op1.direct(noisy_data).power(2).sum(axis=0)).sqrt().sum()) +# +# print( 5*(op2.direct(d) - noisy_data).power(2).sum(), f2(d)) +# +# from functions import mixed_L12Norm as mixed_L12Norm_old +# +# print(mixed_L12Norm_old(op1,None,alpha)(noisy_data)) + + + # + + diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py new file mode 100644 index 0000000..16cd215 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -0,0 +1,314 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:51:17 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator +from ccpi.optimisation.ops import PowerMethodNonsquare +from ccpi.framework import ImageData, BlockDataContainer +import numpy as np + +class FiniteDiff(Operator): + + # Works for Neum/Symmetric & periodic boundary conditions + # TODO add central differences??? + # TODO not very well optimised, too many conditions + # TODO add discretisation step, should get that from imageGeometry + + # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x'] + # Grad_order = ['channels', 'direction_y', 'direction_x'] + # Grad_order = ['direction_z', 'direction_y', 'direction_x'] + # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x'] + + def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'): + + super(FiniteDiff, self).__init__() + + self.gm_domain = gm_domain + self.gm_range = gm_range + self.direction = direction + self.bnd_cond = bnd_cond + + if self.gm_range is None: + self.gm_range = self.gm_domain + + if self.direction + 1 > len(gm_domain): + raise ValueError('Gradient directions more than geometry domain') + + def direct(self, x, out=None): + +# x_asarr = x.as_array() + x_asarr = x + x_sz = len(x.shape) + + if out is None: + out = np.zeros(x.shape) + + fd_arr = out + + ######################## Direct for 2D ############################### + if x_sz == 2: + + if self.direction == 1: + + np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = fd_arr[:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0], x_asarr[:,-1], out = fd_arr[:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 0: + + np.subtract( x_asarr[1:], x_asarr[0:-1], out = fd_arr[0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:], x_asarr[-1,:], out = fd_arr[-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Direct for 3D ############################### + elif x_sz == 3: + + if self.direction == 0: + + np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = fd_arr[0:-1,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = fd_arr[-1,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + + np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = fd_arr[:,0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = fd_arr[:,-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + + if self.direction == 2: + + np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = fd_arr[:,:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = fd_arr[:,:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Direct for 4D ############################### + elif x_sz == 4: + + if self.direction == 0: + np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = fd_arr[0:-1,:,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = fd_arr[-1,:,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = fd_arr[:,0:-1,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = fd_arr[:,-1,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = fd_arr[:,:,0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = fd_arr[:,:,-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 3: + np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = fd_arr[:,:,:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = fd_arr[:,:,:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + else: + raise NotImplementedError + + res = out + return res + + def adjoint(self, x, out=None): + +# x_asarr = x.as_array() + x_asarr = x + x_sz = len(x.shape) + + if out is None: + out = np.zeros(x.shape) + + fd_arr = out + + ######################## Adjoint for 2D ############################### + if x_sz == 2: + + if self.direction == 1: + + np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = fd_arr[:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0], 0, out = fd_arr[:,0] ) + np.subtract( -x_asarr[:,-2], 0, out = fd_arr[:,-1] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0], x_asarr[:,-1], out = fd_arr[:,0] ) + + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 0: + + np.subtract( x_asarr[1:,:], x_asarr[0:-1,:], out = fd_arr[1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:], 0, out = fd_arr[0,:] ) + np.subtract( -x_asarr[-2,:], 0, out = fd_arr[-1,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:], x_asarr[-1,:], out = fd_arr[0,:] ) + + else: + raise ValueError('No valid boundary conditions') + + ######################## Adjoint for 3D ############################### + elif x_sz == 3: + + if self.direction == 0: + + np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = fd_arr[1:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:,:], 0, out = fd_arr[0,:,:] ) + np.subtract( -x_asarr[-2,:,:], 0, out = fd_arr[-1,:,:] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = fd_arr[0,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = fd_arr[:,1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0,:], 0, out = fd_arr[:,0,:] ) + np.subtract( -x_asarr[:,-2,:], 0, out = fd_arr[:,-1,:] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = fd_arr[:,0,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = fd_arr[:,:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,0], 0, out = fd_arr[:,:,0] ) + np.subtract( -x_asarr[:,:,-2], 0, out = fd_arr[:,:,-1] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = fd_arr[:,:,0] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Adjoint for 4D ############################### + elif x_sz == 4: + + if self.direction == 0: + np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = fd_arr[1:,:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:,:,:], 0, out = fd_arr[0,:,:,:] ) + np.subtract( -x_asarr[-2,:,:,:], 0, out = fd_arr[-1,:,:,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = fd_arr[0,:,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = fd_arr[:,1:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0,:,:], 0, out = fd_arr[:,0,:,:] ) + np.subtract( -x_asarr[:,-2,:,:], 0, out = fd_arr[:,-1,:,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = fd_arr[:,0,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = fd_arr[:,:,1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,0,:], 0, out = fd_arr[:,:,0,:] ) + np.subtract( -x_asarr[:,:,-2,:], 0, out = fd_arr[:,:,-1,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = fd_arr[:,:,0,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 3: + np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = fd_arr[:,:,:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,:,0], 0, out = fd_arr[:,:,:,0] ) + np.subtract( -x_asarr[:,:,:,-2], 0, out = fd_arr[:,:,:,-1] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = fd_arr[:,:,:,0] ) + else: + raise ValueError('No valid boundary conditions') + + else: + raise NotImplementedError + + res = out + return res + + def range_dim(self): + return self.gm_range + + def domain_dim(self): + return self.gm_domain + + def norm(self): + x0 = ImageData(np.random.random_sample(self.domain_dim())) + self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0) + return self.s1 + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py new file mode 100644 index 0000000..3dcc1bd --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:50:04 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator +from ccpi.optimisation.ops import PowerMethodNonsquare +from ccpi.framework import ImageData, BlockDataContainer +import numpy as np +from ccpi.optimisation.operators import FiniteDiff + +#%% + +class Gradient(Operator): + + def __init__(self, gm_domain, gm_range=None, bnd_cond = 'Neumann', **kwargs): + + super(Gradient, self).__init__() + + self.gm_domain = gm_domain # Domain of Grad Operator + self.gm_range = gm_range # Range of Grad Operator + self.bnd_cond = bnd_cond # Boundary conditions of Finite Differences + + + if self.gm_range is None: + self.gm_range = ((len(self.gm_domain),)+self.gm_domain) + + # Kwargs Default options + self.memopt = kwargs.get('memopt',False) + self.correlation = kwargs.get('correlation','Space') + + #TODO not tested yet, operator norm??? + self.voxel_size = kwargs.get('voxel_size',[1]*len(gm_domain)) + + + def direct(self, x, out=None): + + tmp = np.zeros(self.gm_range) + for i in range(len(self.gm_domain)): + tmp[i] = FiniteDiff(self.gm_domain, direction = i, bnd_cond = self.bnd_cond).direct(x.as_array())/self.voxel_size[i] +# return type(x)(tmp) + return type(x)(tmp) + + def adjoint(self, x, out=None): + + tmp = np.zeros(self.gm_domain) + for i in range(len(self.gm_domain)): + tmp+=FiniteDiff(self.gm_domain, direction = i, bnd_cond = self.bnd_cond).adjoint(x.as_array()[i])/self.voxel_size[i] + return type(x)(-tmp) + + def alloc_domain_dim(self): + return ImageData(np.zeros(self.gm_domain)) + + def alloc_range_dim(self): + return ImageData(np.zeros(self.range_dim)) + + def domain_dim(self): + return self.gm_domain + + def range_dim(self): + return self.gm_range + + def norm(self): +# return np.sqrt(4*len(self.domainDim())) + #TODO this takes time for big ImageData + # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12) + x0 = ImageData(np.random.random_sample(self.domain_dim())) + self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0) + return self.s1 + + +if __name__ == '__main__': + + N, M = (200,300) + ig = (N,M) + G = Gradient(ig) + u = DataContainer(np.random.randint(10, size=G.domain_dim())) + w = DataContainer(np.random.randint(10, size=G.range_dim())) +# w = [DataContainer(np.random.randint(10, size=G.domain_dim())),\ +# DataContainer(np.random.randint(10, size=G.domain_dim()))] + + # domain_dim + print('Domain {}'.format(G.domain_dim())) + + # range_dim + print('Range {}'.format(G.range_dim())) + + # Direct + z = G.direct(u) + + # Adjoint + z1 = G.adjoint(w) + + print(z) + print(z1) + + LHS = (G.direct(u)*w).sum() + RHS = (u * G.adjoint(w)).sum() +# + print(LHS,RHS) + print(G.norm()) + +# print(G.adjoint(G.direct(u))) + + + + + + + + + + + + + + + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py new file mode 100644 index 0000000..d49cb30 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:30:51 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator + + +class Identity(Operator): + + def __init__(self, gm_domain, gm_range=None): + + self.gm_domain = gm_domain + self.gm_range = gm_range + if self.gm_range is None: + self.gm_range = self.gm_domain + + super(Identity, self).__init__() + + def direct(self,x,out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def adjoint(self,x, out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def norm(self): + return 1.0 + + def domain_dim(self): + return self.gm_domain + + def range_dim(self): + return self.gm_range
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py new file mode 100644 index 0000000..d908e49 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:53:55 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator +from ccpi.optimisation.operators import FiniteDiff +from ccpi.optimisation.ops import PowerMethodNonsquare +from ccpi.framework import ImageData, DataContainer +import numpy as np + + +class SymmetrizedGradient(Operator): + + def __init__(self, gm_domain, gm_range, bnd_cond = 'Neumann', **kwargs): + + super(SymmetrizedGradient, self).__init__() + + self.gm_domain = gm_domain # Domain of Grad Operator + self.gm_range = gm_range # Range of Grad Operator + self.bnd_cond = bnd_cond # Boundary conditions of Finite Differences + + # Kwargs Default options + self.memopt = kwargs.get('memopt',False) + self.correlation = kwargs.get('correlation','Space') + + #TODO not tested yet, operator norm??? + self.voxel_size = kwargs.get('voxel_size',[1]*len(gm_domain)) + + + def direct(self, x, out=None): + + tmp = np.zeros(self.gm_range) + tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) + tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) + tmp[2] = 0.5 * (FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) + + FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) ) + + return type(x)(tmp) + + + def adjoint(self, x, out=None): + + tmp = np.zeros(self.gm_domain) + + tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) + \ + FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) + + tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) + \ + FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1]) + + return type(x)(tmp) + + def alloc_domain_dim(self): + return ImageData(np.zeros(self.gm_domain)) + + def alloc_range_dim(self): + return ImageData(np.zeros(self.range_dim)) + + def domain_dim(self): + return self.gm_domain + + def range_dim(self): + return self.gm_range + + def norm(self): +# return np.sqrt(4*len(self.domainDim())) + #TODO this takes time for big ImageData + # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12) + x0 = ImageData(np.random.random_sample(self.domain_dim())) + self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0) + return self.s1 + + + +if __name__ == '__main__': + + ########################################################################### + ## Symmetrized Gradient + + N, M = 2, 3 + ig = (N,M) + ig2 = (2,) + ig + ig3 = (3,) + ig + u1 = DataContainer(np.random.randint(10, size=ig2)) + w1 = DataContainer(np.random.randint(10, size=ig3)) + + E = SymmetrizedGradient(ig2,ig3) + + d1 = E.direct(u1) + d2 = E.adjoint(w1) + + LHS = (d1.as_array()[0]*w1.as_array()[0] + \ + d1.as_array()[1]*w1.as_array()[1] + \ + 2*d1.as_array()[2]*w1.as_array()[2]).sum() + + RHS = (u1.as_array()[0]*d2.as_array()[0] + \ + u1.as_array()[1]*d2.as_array()[1]).sum() + + + print(LHS, RHS, E.norm()) + + +# + + + + + + + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py new file mode 100644 index 0000000..a7c5f09 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:25:53 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.framework import ImageData +from ccpi.optimisation.operators import Operator + +class ZeroOp(Operator): + + def __init__(self, gm_domain, gm_range): + self.gm_domain = gm_domain + self.gm_range = gm_range + super(ZeroOp, self).__init__() + + def direct(self,x,out=None): + if out is None: + return ImageData(np.zeros(self.gm_range)) + else: + return ImageData(np.zeros(self.gm_range)) + + def adjoint(self,x, out=None): + if out is None: + return ImageData(np.zeros(self.gm_domain)) + else: + return ImageData(np.zeros(self.gm_domain)) + + def norm(self): + return 0 + + def domain_dim(self): + return self.gm_domain + + def range_dim(self): + return self.gm_range
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py index cc307e0..1e86efc 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py @@ -10,3 +10,11 @@ from .LinearOperator import LinearOperator from .ScaledOperator import ScaledOperator
from .BlockOperator import BlockOperator
from .BlockScaledOperator import BlockScaledOperator
+
+
+from .FiniteDifferenceOperator import FiniteDiff
+from .GradientOperator import Gradient
+from .SymmetrizedGradientOperator import SymmetrizedGradient
+from .IdentityOperator import Identity
+from .ZeroOperator import ZeroOp
+
diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index 630e33e..87930b5 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -34,7 +34,8 @@ setup( packages=['ccpi' , 'ccpi.io', 'ccpi.framework', 'ccpi.optimisation', 'ccpi.optimisation.operators', - 'ccpi.optimisation.algorithms'], + 'ccpi.optimisation.algorithms', + 'ccpi.optimisation.functions'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine |