From ab405e3d4714611ac3964f5a563d809977d8daa6 Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com> Date: Fri, 5 Jul 2019 10:46:21 +0100 Subject: Fix documentation (#343) * fix documentation algs * fix docum functions dir * fix docum operator dir * docstring title * updated math in docstring --- .../ccpi/optimisation/algorithms/Algorithm.py | 34 +++++---- .../Python/ccpi/optimisation/algorithms/CGLS.py | 59 ++++++++------- .../Python/ccpi/optimisation/algorithms/FISTA.py | 72 +++++++++++-------- .../optimisation/algorithms/GradientDescent.py | 43 ++++++----- .../Python/ccpi/optimisation/algorithms/PDHG.py | 81 ++++++++++++++------- .../Python/ccpi/optimisation/algorithms/SIRT.py | 57 +++++++++------ .../ccpi/optimisation/functions/BlockFunction.py | 78 ++++++++++---------- .../Python/ccpi/optimisation/functions/Function.py | 36 +++++----- .../functions/FunctionOperatorComposition.py | 56 ++++++++------- .../ccpi/optimisation/functions/IndicatorBox.py | 38 +++++++--- .../ccpi/optimisation/functions/KullbackLeibler.py | 83 ++++++++++++---------- .../Python/ccpi/optimisation/functions/L1Norm.py | 76 ++++++++++++-------- .../ccpi/optimisation/functions/L2NormSquared.py | 83 ++++++++++------------ .../ccpi/optimisation/functions/MixedL21Norm.py | 83 +++++++++++----------- .../Python/ccpi/optimisation/functions/Norm2Sq.py | 39 +++++----- .../ccpi/optimisation/functions/ScaledFunction.py | 32 +++++---- .../ccpi/optimisation/functions/ZeroFunction.py | 72 +++++++++++++------ .../ccpi/optimisation/operators/BlockOperator.py | 45 ++++++------ .../optimisation/operators/BlockScaledOperator.py | 71 +++++++++--------- .../operators/FiniteDifferenceOperator.py | 48 ++++++++----- .../optimisation/operators/GradientOperator.py | 40 +++++++++-- .../optimisation/operators/IdentityOperator.py | 39 ++++++++-- .../ccpi/optimisation/operators/LinearOperator.py | 1 - .../optimisation/operators/LinearOperatorMatrix.py | 8 +-- .../ccpi/optimisation/operators/ScaledOperator.py | 26 ++++--- .../optimisation/operators/ShrinkageOperator.py | 5 ++ .../optimisation/operators/SparseFiniteDiff.py | 4 ++ .../operators/SymmetrizedGradientOperator.py | 39 +++++----- .../ccpi/optimisation/operators/ZeroOperator.py | 37 +++++++++- docs/source/optimisation.rst | 4 ++ 30 files changed, 837 insertions(+), 552 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 1168db3..fec37c5 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -1,20 +1,26 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. import time, functools from numbers import Integral diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py index 82fbb96..c6c1d4c 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -1,40 +1,47 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. from ccpi.optimisation.algorithms import Algorithm -from ccpi.optimisation.functions import Norm2Sq import numpy class CGLS(Algorithm): - '''Conjugate Gradient Least Squares algorithm + r'''Conjugate Gradient Least Squares algorithm + + Problem: + + .. math:: - Parameters: - x_init: initial guess - operator: operator for forward/backward projections - data: data to operate on - tolerance: tolerance to stop algorithm + \min || A x - b ||^2_2 + + | + + Parameters : + + :parameter operator : Linear operator for the inverse problem + :parameter x_init : Initial guess ( Default x_init = 0) + :parameter data : Acquired data to reconstruct + :parameter tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm Reference: https://web.stanford.edu/group/SOL/software/cgls/ - ''' - def __init__(self, **kwargs): super(CGLS, self).__init__() @@ -59,9 +66,7 @@ class CGLS(Algorithm): self.p = self.s self.norms0 = self.s.norm() - ## self.norms = self.s.norm() - ## self.gamma = self.norms0**2 self.normx = self.x.norm() diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py index fa1d8d8..347dacd 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -1,43 +1,60 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. from ccpi.optimisation.algorithms import Algorithm from ccpi.optimisation.functions import ZeroFunction import numpy class FISTA(Algorithm): - '''Fast Iterative Shrinkage-Thresholding Algorithm - Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding - algorithm for linear inverse problems. - SIAM journal on imaging sciences,2(1), pp.183-202. + r'''Fast Iterative Shrinkage-Thresholding Algorithm + + Problem : + + .. math:: + + \min_{x} f(x) + g(x) + + | - Parameters: - x_init: initial guess - f: data fidelity - g: regularizer - opt: additional options + Parameters : + + :parameter x_init : Initial guess ( Default x_init = 0) + :parameter f : Differentiable function + :parameter g : Convex function with " simple " proximal operator + + + Reference: + + Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding + algorithm for linear inverse problems. + SIAM journal on imaging sciences,2(1), pp.183-202. ''' def __init__(self, **kwargs): - '''initialisation can be done at creation time if all + + '''creator + + initialisation can be done at creation time if all proper variables are passed or later with set_up''' + super(FISTA, self).__init__() f = kwargs.get('f', None) g = kwargs.get('g', ZeroFunction()) @@ -67,12 +84,11 @@ class FISTA(Algorithm): self.f.gradient(self.y, out=self.u) self.u.__imul__( -self.invL ) self.u.__iadd__( self.y ) - # x = g.prox(u,invL) + self.g.proximal(self.u, self.invL, out=self.x) self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) -# self.x.subtract(self.x_old, out=self.y) self.y = self.x - self.x_old self.y.__imul__ ((self.t_old-1)/self.t) self.y.__iadd__( self.x ) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py index 8c2b693..d060690 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py @@ -1,24 +1,33 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from ccpi.optimisation.algorithms import Algorithm class GradientDescent(Algorithm): - '''Implementation of Gradient Descent algorithm + ''' + + Gradient Descent algorithm + + ''' def __init__(self, **kwargs): @@ -35,7 +44,7 @@ class GradientDescent(Algorithm): self.set_up(x_init=x_init, objective_function=objective_function, rate=rate) def should_stop(self): - '''stopping cryterion, currently only based on number of iterations''' + '''stopping criterion, currently only based on number of iterations''' return self.iteration >= self.max_iteration def set_up(self, x_init, objective_function, rate): diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 8f37765..0968872 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -1,32 +1,60 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from ccpi.optimisation.algorithms import Algorithm -from ccpi.framework import ImageData, DataContainer -import numpy as np -import numpy -import time -from ccpi.optimisation.operators import BlockOperator -from ccpi.framework import BlockDataContainer -from ccpi.optimisation.functions import FunctionOperatorComposition + class PDHG(Algorithm): - '''Primal Dual Hybrid Gradient''' + r'''Primal Dual Hybrid Gradient + + Problem: + + .. math:: + + \min_{x} f(Kx) + g(x) + | + + Parameters : + + :parameter operator : Linear Operator = K + :parameter f : Convex function with "simple" proximal of its conjugate. + :parameter g : Convex function with "simple" proximal + :parameter sigma : Step size parameter for Primal problem + :parameter tau : Step size parameter for Dual problem + + Remark: Convergence is guaranted provided that + + .. math:: \tau \sigma \|K\|^{2} <1 + + + Reference : + + + (a) A. Chambolle and T. Pock (2011), "A first-order primal–dual algorithm for convex + problems with applications to imaging", J. Math. Imaging Vision 40, 120–145. + + + (b) E. Esser, X. Zhang and T. F. Chan (2010), "A general framework for a class of first + order primal–dual algorithms for convex optimization in imaging science", + SIAM J. Imaging Sci. 3, 1015–1046. + ''' def __init__(self, **kwargs): super(PDHG, self).__init__(max_iteration=kwargs.get('max_iteration',0)) @@ -78,7 +106,7 @@ class PDHG(Algorithm): def update(self): - # Gradient descent, Dual problem solution + # Gradient ascent for the dual variable self.operator.direct(self.xbar, out=self.y_tmp) self.y_tmp *= self.sigma self.y_tmp += self.y_old @@ -86,7 +114,7 @@ class PDHG(Algorithm): # self.y = self.f.proximal_conjugate(self.y_old, self.sigma) self.f.proximal_conjugate(self.y_tmp, self.sigma, out=self.y) - # Gradient ascent, Primal problem solution + # Gradient descent for the primal variable self.operator.adjoint(self.y, out=self.x_tmp) self.x_tmp *= -1*self.tau self.x_tmp += self.x_old @@ -98,6 +126,7 @@ class PDHG(Algorithm): self.xbar *= self.theta self.xbar += self.x + self.x_old.fill(self.x) self.y_old.fill(self.y) @@ -106,4 +135,4 @@ class PDHG(Algorithm): p1 = self.f(self.operator.direct(self.x)) + self.g(self.x) d1 = -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(-1*self.operator.adjoint(self.y))) - self.loss.append([p1, d1, p1-d1]) \ No newline at end of file + self.loss.append([p1, d1, p1-d1]) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index ca5b084..2b49ab0 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -1,33 +1,46 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= from ccpi.optimisation.algorithms import Algorithm class SIRT(Algorithm): - '''Simultaneous Iterative Reconstruction Technique + r'''Simultaneous Iterative Reconstruction Technique + + Problem: + + .. math:: + + A x = b + | Parameters: - x_init: initial guess - operator: operator for forward/backward projections - data: data to operate on - constraint: Function with prox-method, for example IndicatorBox to - enforce box constraints, default is None). + + :parameter operator : Linear operator for the inverse problem + :parameter x_init : Initial guess + :parameter data : Acquired data to reconstruct + :parameter constraint : Function proximal method + e.g. x\in[0, 1], IndicatorBox to enforce box constraints + Default is None). ''' def __init__(self, **kwargs): super(SIRT, self).__init__() diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py index 6d4b4e1..a6ac66c 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -1,20 +1,24 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= from ccpi.optimisation.functions import Function from ccpi.framework import BlockDataContainer @@ -22,29 +26,29 @@ from numbers import Number class BlockFunction(Function): - '''BlockFunction acts as a separable sum function, i.e., + r'''BlockFunction acts as a separable sum function: f = [f_1,...,f_n] - f = [f_1,...,f_n] - - f([x_1,...,x_n]) = f_1(x_1) + .... + f_n(x_n) + .. math:: + + f([x_1,...,x_n]) = f_1(x_1) + .... + f_n(x_n) + | ''' + def __init__(self, *functions): super(BlockFunction, self).__init__() self.functions = functions self.length = len(self.functions) - - - - + def __call__(self, x): - '''Evaluates the BlockFunction at a BlockDataContainer x + r'''Evaluates the BlockFunction at a BlockDataContainer x - :param: x (BlockDataContainer): must have as many rows as self.length + :param: x (BlockDataContainer): must have as many rows as self.length - returns sum(f_i(x_i)) + returns ..math:: sum(f_i(x_i)) + ''' if self.length != x.shape[0]: @@ -56,9 +60,9 @@ class BlockFunction(Function): def convex_conjugate(self, x): - ''' Evaluate convex conjugate of BlockFunction at x + r'''Convex conjugate of BlockFunction at x - returns sum(f_i^{*}(x_i)) + .. math:: returns sum(f_i^{*}(x_i)) ''' t = 0 @@ -69,9 +73,9 @@ class BlockFunction(Function): def proximal_conjugate(self, x, tau, out = None): - ''' Evaluate Proximal Operator of tau * f(\cdot) at x - - prox_{tau*f}(x) = sum_{i} prox_{tau*f_{i}}(x_{i}) + r'''Proximal operator of BlockFunction at x: + + .. math:: prox_{tau*f}(x) = sum_{i} prox_{tau*f_{i}}(x_{i}) ''' @@ -99,11 +103,9 @@ class BlockFunction(Function): def proximal(self, x, tau, out = None): - ''' Evaluate Proximal Operator of tau * f^{*}(\cdot) at x - - prox_{tau*f^{*}}(x) = sum_{i} prox_{tau*f^{*}_{i}}(x_{i}) - + r'''Proximal operator of the convex conjugate of BlockFunction at x: + .. math:: prox_{tau*f^{*}}(x) = sum_{i} prox_{tau*f^{*}_{i}}(x_{i}) ''' if out is None: @@ -130,9 +132,9 @@ class BlockFunction(Function): def gradient(self,x, out=None): - ''' Evaluate gradient of f at x: f'(x) + r'''Evaluates gradient of BlockFunction at x - returns: BlockDataContainer [f_{1}'(x_{1}), ... , f_{n}'(x_{n})] + returns: BlockDataContainer .. math:: [f_{1}'(x_{1}), ... , f_{n}'(x_{n})] ''' diff --git a/Wrappers/Python/ccpi/optimisation/functions/Function.py b/Wrappers/Python/ccpi/optimisation/functions/Function.py index 9a6009f..7156995 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/Function.py +++ b/Wrappers/Python/ccpi/optimisation/functions/Function.py @@ -1,20 +1,24 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= import warnings from ccpi.optimisation.functions.ScaledFunction import ScaledFunction diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py index 24ebbf0..58d4f27 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -1,30 +1,35 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from ccpi.optimisation.functions import Function from ccpi.optimisation.functions import ScaledFunction class FunctionOperatorComposition(Function): - ''' Function composition with Operator, i.e., f(Ax) + '''Function composition with Operator: (f o A)(x) = f(Ax) - A: operator - f: function + : parameter A: operator + : parameter f: function ''' @@ -39,17 +44,18 @@ class FunctionOperatorComposition(Function): def __call__(self, x): - ''' Evaluate FunctionOperatorComposition at x - - returns f(Ax) - - ''' + '''Evaluates f(Ax)''' return self.function(self.operator.direct(x)) def gradient(self, x, out=None): -# - ''' Gradient takes into account the Operator''' + + '''Evaluates gradient of f(Ax): + + ..math :: A^{T}f'(Ax) + + ''' + if out is None: return self.operator.adjoint(self.function.gradient(self.operator.direct(x))) else: diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py index 2001a56..51d08d1 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py +++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py @@ -18,24 +18,31 @@ from ccpi.optimisation.functions import Function import numpy -from ccpi.framework import ImageData class IndicatorBox(Function): - '''Box constraints indicator function. - Calling returns 0 if argument is within the box. The prox operator is projection onto the box. - Only implements one scalar lower and one upper as constraint on all elements. Should generalise - to vectors to allow different constraints one elements. -''' + + r'''Indicator function for box constraint + + .. math:: + f(x) = \mathbb{I}_{[a, b]} = \begin{cases} + + 0, if x\in[a, b] + \infty, otherwise + \end{cases} + + ''' def __init__(self,lower=-numpy.inf,upper=numpy.inf): - # Do nothing + super(IndicatorBox, self).__init__() self.lower = lower self.upper = upper - + def __call__(self,x): + + '''Evaluates IndicatorBox at x''' if (numpy.all(x.array>=self.lower) and numpy.all(x.array <= self.upper) ): @@ -48,12 +55,18 @@ class IndicatorBox(Function): return ValueError('Not Differentiable') def convex_conjugate(self,x): - # support function sup , z \in [lower, upper] - # ???? + + '''Convex conjugate of IndicatorBox at x''' + return x.maximum(0).sum() def proximal(self, x, tau, out=None): + r'''Proximal operator of IndicatorBox at x + + .. math:: prox_{\tau * f}(x) + ''' + if out is None: return (x.maximum(self.lower)).minimum(self.upper) else: @@ -61,6 +74,11 @@ class IndicatorBox(Function): out.minimum(self.upper, out=out) def proximal_conjugate(self, x, tau, out=None): + + r'''Proximal operator of the convex conjugate of IndicatorBox at x: + + ..math:: prox_{\tau * f^{*}}(x) + ''' if out is None: diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 2a87205..f88c339 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -1,20 +1,22 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. import numpy from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction @@ -23,13 +25,14 @@ import scipy.special class KullbackLeibler(Function): - ''' - - KL_div(x, y + back) = int x * log(x/(y+back)) - x + (y+back) + r'''Kullback-Leibler divergence function - Assumption: y>=0 - back>=0 - + .. math:: + f(x, y) = \begin{cases} x \log(x / y) - x + y & x > 0, y > 0 \\ + y & x = 0, y \ge 0 \\ + \infty & \text{otherwise} + \end{cases} + ''' def __init__(self, data, **kwargs): @@ -42,29 +45,24 @@ class KullbackLeibler(Function): def __call__(self, x): - - ''' - - x - y * log( x + bnoise) + y * log(y) - y + bnoise - - - ''' - + + '''Evaluates KullbackLeibler at x''' + ind = x.as_array()>0 tmp = scipy.special.kl_div(self.b.as_array()[ind], x.as_array()[ind]) return numpy.sum(tmp) - def log(self, datacontainer): '''calculates the in-place log of the datacontainer''' - if not functools.reduce(lambda x,y: x and y>0, - datacontainer.as_array().ravel(), True): + if not functools.reduce(lambda x,y: x and y>0, datacontainer.as_array().ravel(), True): raise ValueError('KullbackLeibler. Cannot calculate log of negative number') datacontainer.fill( numpy.log(datacontainer.as_array()) ) def gradient(self, x, out=None): + '''Evaluates gradient of KullbackLeibler at x''' + if out is None: return 1 - self.b/(x + self.bnoise) else: @@ -76,11 +74,19 @@ class KullbackLeibler(Function): def convex_conjugate(self, x): + '''Convex conjugate of KullbackLeibler at x''' + xlogy = - scipy.special.xlogy(self.b.as_array(), 1 - x.as_array()) return numpy.sum(xlogy) def proximal(self, x, tau, out=None): + r'''Proximal operator of KullbackLeibler at x + + .. math:: prox_{\tau * f}(x) + + ''' + if out is None: return 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b ) .sqrt() ) else: @@ -103,6 +109,11 @@ class KullbackLeibler(Function): out *= 0.5 def proximal_conjugate(self, x, tau, out=None): + + r'''Proximal operator of the convex conjugate of KullbackLeibler at x: + + .. math:: prox_{\tau * f^{*}}(x) + ''' if out is None: @@ -110,7 +121,6 @@ class KullbackLeibler(Function): return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) else: - #tmp = x + tau * self.bnoise tmp = tau * self.bnoise tmp += x tmp -= 1 @@ -126,10 +136,9 @@ class KullbackLeibler(Function): def __rmul__(self, scalar): - ''' Multiplication of L2NormSquared with a scalar - - Returns: ScaledFunction - + '''Multiplication of KullbackLeibler with a scalar + + Returns: ScaledFunction ''' return ScaledFunction(self, scalar) diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py index 08c3d17..97d03b9 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -1,20 +1,21 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction @@ -23,16 +24,14 @@ from ccpi.optimisation.operators import ShrinkageOperator class L1Norm(Function): - ''' - - Class: L1Norm - - Cases: a) f(x) = ||x||_{1} - - b) f(x) = ||x - b||_{1} - - ''' - + r'''L1Norm function: + + Cases considered (with/without data): + a) .. math:: f(x) = ||x||_{1} + b) .. math:: f(x) = ||x - b||_{1} + + ''' + def __init__(self, **kwargs): super(L1Norm, self).__init__() @@ -40,7 +39,7 @@ class L1Norm(Function): def __call__(self, x): - ''' Evaluate L1Norm at x: f(x) ''' + '''Evaluates L1Norm at x''' y = x if self.b is not None: @@ -48,11 +47,12 @@ class L1Norm(Function): return y.abs().sum() def gradient(self,x): - #TODO implement subgradient??? + return ValueError('Not Differentiable') def convex_conjugate(self,x): - #TODO implement Indicator infty??? + + '''Convex conjugate of L1Norm at x''' y = 0 if self.b is not None: @@ -61,7 +61,11 @@ class L1Norm(Function): def proximal(self, x, tau, out=None): - # TODO implement shrinkage operator, we will need it later e.g SplitBregman + r'''Proximal operator of L1Norm at x + + ..math:: prox_{\tau * f}(x) + + ''' if out is None: if self.b is not None: @@ -76,6 +80,12 @@ class L1Norm(Function): def proximal_conjugate(self, x, tau, out=None): + r'''Proximal operator of the convex conjugate of L1Norm at x: + + .. math:: prox_{\tau * f^{*}}(x) + + ''' + if out is None: if self.b is not None: return (x - tau*self.b).divide((x - tau*self.b).abs().maximum(1.0)) @@ -88,6 +98,12 @@ class L1Norm(Function): out.fill(x.divide(x.abs().maximum(1.0)) ) def __rmul__(self, scalar): + + '''Multiplication of L2NormSquared with a scalar + + Returns: ScaledFunction + ''' + return ScaledFunction(self, scalar) diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 2afae25..a625f07 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -1,44 +1,47 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction -from ccpi.optimisation.functions import FunctionOperatorComposition class L2NormSquared(Function): - ''' + r'''L2NormSquared function: - Cases: a) f(x) = \|x\|^{2}_{2} - - b) f(x) = ||x - b||^{2}_{2} - + Cases considered (with/without data): + a) .. math:: f(x) = \|x\|^{2}_{2} + b) .. math:: f(x) = \|\|x - b\|\|^{2}_{2} + ''' def __init__(self, **kwargs): super(L2NormSquared, self).__init__() - self.b = kwargs.get('b',None) + self.b = kwargs.get('b',None) + + # Lipschitz constant self.L = 2 def __call__(self, x): - ''' Evaluate L2NormSquared at x: f(x) ''' + '''Evaluates L2NormSquared at x''' y = x if self.b is not None: @@ -51,7 +54,7 @@ class L2NormSquared(Function): def gradient(self, x, out=None): - ''' Evaluate gradient of L2NormSquared at x: f'(x) ''' + '''Evaluates gradient of L2NormSquared at x''' if out is not None: @@ -70,7 +73,7 @@ class L2NormSquared(Function): def convex_conjugate(self, x): - ''' Evaluate convex conjugate of L2NormSquared at x: f^{*}(x)''' + '''Convex conjugate of L2NormSquared at x''' tmp = 0 @@ -82,11 +85,10 @@ class L2NormSquared(Function): def proximal(self, x, tau, out = None): - ''' Evaluate Proximal Operator of tau * f(\cdot) at x: - - prox_{tau*f(\cdot)}(x) = \argmin_{z} \frac{1}{2}|| z - x ||^{2}_{2} + tau * f(z) - - ''' + r'''Proximal operator of L2NormSquared at x + + .. math:: prox_{\tau * f}(x) + ''' if out is None: @@ -109,11 +111,9 @@ class L2NormSquared(Function): def proximal_conjugate(self, x, tau, out=None): - ''' Evaluate Proximal Operator of tau * f^{*}(\cdot) at x (i.e., the convex conjugate of f) : - - prox_{tau*f(\cdot)}(x) = \argmin_{z} \frac{1}{2}|| z - x ||^{2}_{2} + tau * f^{*}(z) - - ''' + r'''Proximal operator of the convex conjugate of L2NormSquared at x: + + .. math:: prox_{\tau * f^{*}}(x)''' if out is None: if self.b is not None: @@ -129,20 +129,13 @@ class L2NormSquared(Function): def __rmul__(self, scalar): - ''' Multiplication of L2NormSquared with a scalar - - Returns: ScaledFunction - - ''' + '''Multiplication of L2NormSquared with a scalar + + Returns: ScaledFunction''' return ScaledFunction(self, scalar) - def composition(self, operator): - - return FunctionOperatorComposition(operator) - - if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py index d7d23a5..8cbed67 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py @@ -1,20 +1,22 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from ccpi.optimisation.functions import Function, ScaledFunction from ccpi.framework import BlockDataContainer @@ -24,8 +26,10 @@ import functools class MixedL21Norm(Function): - ''' - f(x) = ||x||_{2,1} = \sum |x|_{2} + r'''MixedL21Norm: .. math:: f(x) = ||x||_{2,1} = \int \|x\|_{2} dx + + where x is a vector/tensor vield + ''' def __init__(self, **kwargs): @@ -35,10 +39,9 @@ class MixedL21Norm(Function): def __call__(self, x): - ''' Evaluates L2,1Norm at point x + '''Evaluates MixedL21Norm at point x - :param: x is a BlockDataContainer - + :param: x: is a BlockDataContainer ''' if not isinstance(x, BlockDataContainer): raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) @@ -53,26 +56,31 @@ class MixedL21Norm(Function): def convex_conjugate(self,x): - ''' This is the Indicator function of ||\cdot||_{2, \infty} - which is either 0 if ||x||_{2, \infty} or \infty + r'''Convex conjugate of of MixedL21Norm: + + Indicator function of .. math:: ||\cdot||_{2, \infty} + which is either 0 if .. math:: ||x||_{2, \infty}<1 or \infty + ''' return 0.0 - #tmp = [ el**2 for el in x.containers ] - #print(sum(tmp).sqrt().as_array().max()) - #return sum(tmp).sqrt().as_array().max() def proximal(self, x, tau, out=None): - ''' - For this we need to define a MixedL2,2 norm acting on BDC, - different form L2NormSquared which acts on DC - + r'''Proximal operator of MixedL21Norm at x: + + .. math:: prox_{\tau * f(x) ''' pass def proximal_conjugate(self, x, tau, out=None): + + r'''Proximal operator of the convex conjugate of MixedL21Norm at x: + + .. math:: prox_{\tau * f^{*}}(x) + + ''' if out is None: @@ -80,26 +88,19 @@ class MixedL21Norm(Function): res = sum(tmp).sqrt().maximum(1.0) frac = [el/res for el in x.containers] return BlockDataContainer(*frac) - - - #TODO this is slow, why??? -# return x.divide(x.pnorm().maximum(1.0)) + else: res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ) res = res1.sqrt().maximum(1.0) x.divide(res, out=out) - -# x.divide(sum([el*el for el in x.containers]).sqrt().maximum(1.0), out=out) - #TODO this is slow, why ??? -# x.divide(x.pnorm().maximum(1.0), out=out) - + def __rmul__(self, scalar): - ''' Multiplication of L2NormSquared with a scalar - - Returns: ScaledFunction + '''Multiplication of MixedL21Norm with a scalar + + Returns: ScaledFunction ''' return ScaledFunction(self, scalar) diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py index eea300d..c8f5e46 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py +++ b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py @@ -1,36 +1,37 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from ccpi.optimisation.functions import Function -import numpy import warnings # Define a class for squared 2-norm class Norm2Sq(Function): - ''' - f(x) = c*||A*x-b||_2^2 + r'''.. math:: f(x) = c*||A*x-b||_2^2 which has - grad[f](x) = 2*c*A^T*(A*x-b) + .. math:: grad[f](x) = 2*c*A^T*(A*x-b) and Lipschitz constant - L = 2*c*||A||_2^2 = 2*s1(A)^2 + .. math:: L = 2*c*||A||_2^2 = 2*s1(A)^2 where s1(A) is the largest singular value of A. diff --git a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py index df008e0..3e689e6 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py @@ -1,20 +1,22 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from numbers import Number import numpy import warnings diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py index fb990f8..ca52f31 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py @@ -1,53 +1,81 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from ccpi.optimisation.functions import Function -from ccpi.framework import BlockDataContainer class ZeroFunction(Function): - ''' ZeroFunction: f(x) = 0 - - + r'''ZeroFunction: .. math:: f(x) = 0, + + Maps evely element x\in X to zero ''' def __init__(self): super(ZeroFunction, self).__init__() def __call__(self,x): + + '''Evaluates ZeroFunction at x''' return 0 + + def gradient(self, x, out=None): + + '''Evaluates gradient of ZeroFunction at x''' + + if out is None: + return 0 + else: + out *= 0 + def convex_conjugate(self, x): - ''' This is the support function sup which in fact is the - indicator function for the set = {0} - So 0 if x=0, or inf if x neq 0 + r''' Convex conjugate of ZeroFunction: support function .. math:: sup + + 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, out=None): + + r'''Proximal operator of ZeroFunction at x + + .. math:: prox_{\tau * f}(x) + ''' + if out is None: return x.copy() else: out.fill(x) def proximal_conjugate(self, x, tau, out = None): + + r'''Proximal operator of the convex conjugate of ZeroFunction at x: + + .. math:: prox_{\tau * f^{*}}(x) + + ''' + if out is None: return 0 else: diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index 5dfd6ba..e3a02ec 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -1,37 +1,37 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). - -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. import numpy -from numbers import Number import functools -from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer, DataContainer -from ccpi.optimisation.operators import Operator, LinearOperator -from ccpi.optimisation.operators.BlockScaledOperator import BlockScaledOperator +from ccpi.framework import ImageData, BlockDataContainer, DataContainer +from ccpi.optimisation.operators import Operator from ccpi.framework import BlockGeometry class BlockOperator(Operator): - '''A Block matrix containing Operators + r'''A Block matrix containing Operators The Block Framework is a generic strategy to treat variational problems in the following form: .. math:: - min Regulariser + Fidelity + \min Regulariser + Fidelity BlockOperators have a generic shape M x N, and when applied on an @@ -192,7 +192,6 @@ class BlockOperator(Operator): else: return BlockDataContainer(*res, shape=shape) else: - #tmp = self.domain_geometry().allocate() for col in range(self.shape[1]): for row in range(self.shape[0]): @@ -240,7 +239,7 @@ class BlockOperator(Operator): def __rmul__(self, scalar): '''Defines the left multiplication with a scalar - Args: scalar (number or iterable containing numbers): + :paramer scalar: (number or iterable containing numbers): Returns: a block operator with Scaled Operators inside''' if isinstance (scalar, list) or isinstance(scalar, tuple) or \ diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py index 30b750d..c23c23a 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py @@ -1,44 +1,51 @@ # -*- coding: utf-8 -*- -# CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL). +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. -# Copyright 2017 UKRI-STFC -# Copyright 2017 University of Manchester -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. from numbers import Number import numpy from ccpi.optimisation.operators import ScaledOperator import functools class BlockScaledOperator(ScaledOperator): + '''ScaledOperator - A class to represent the scalar multiplication of an Operator with a scalar. - It holds an operator and a scalar. Basically it returns the multiplication - of the result of direct and adjoint of the operator with the scalar. - For the rest it behaves like the operator it holds. + A class to represent the scalar multiplication of an Operator with a scalar. + It holds an operator and a scalar. Basically it returns the multiplication + of the result of direct and adjoint of the operator with the scalar. + For the rest it behaves like the operator it holds. + + Args: + :param operator (Operator): a Operator or LinearOperator + :param scalar (Number): a scalar multiplier + Example: + The scaled operator behaves like the following: + .. code-block:: python - Args: - operator (Operator): a Operator or LinearOperator - scalar (Number): a scalar multiplier - Example: - The scaled operator behaves like the following: - sop = ScaledOperator(operator, scalar) - sop.direct(x) = scalar * operator.direct(x) - sop.adjoint(x) = scalar * operator.adjoint(x) - sop.norm() = operator.norm() - sop.range_geometry() = operator.range_geometry() - sop.domain_geometry() = operator.domain_geometry() + sop = ScaledOperator(operator, scalar) + sop.direct(x) = scalar * operator.direct(x) + sop.adjoint(x) = scalar * operator.adjoint(x) + sop.norm() = operator.norm() + sop.range_geometry() = operator.range_geometry() + sop.domain_geometry() = operator.domain_geometry() + ''' def __init__(self, operator, scalar, shape=None): if shape is None: @@ -75,10 +82,4 @@ class BlockScaledOperator(ScaledOperator): @property def T(self): '''Return the transposed of self''' - #print ("transpose before" , self.shape) - #shape = (self.shape[1], self.shape[0]) - ##self.shape = shape - ##self.operator.shape = shape - #print ("transpose" , shape) - #return self return type(self)(self.operator.T, self.scalar) \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py index cb8cb65..9b5ae24 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -16,25 +16,31 @@ # See the License for the specific language governing permissions and # limitations under the License. from ccpi.optimisation.operators import LinearOperator -from ccpi.framework import ImageData, BlockDataContainer import numpy as np class FiniteDiff(LinearOperator): - # 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'] + '''Finite Difference Operator: + + Computes first-order forward/backward differences + on 2D, 3D, 4D ImageData + under Neumann/Periodic boundary conditions + + Order of the Gradient ( ImageGeometry may contain channels ): + + 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__() - '''FIXME: domain and range should be geometries''' + self.gm_domain = gm_domain self.gm_range = gm_range @@ -44,6 +50,7 @@ class FiniteDiff(LinearOperator): # Domain Geometry = Range Geometry if not stated if self.gm_range is None: self.gm_range = self.gm_domain + # check direction and "length" of geometry if self.direction + 1 > len(self.gm_domain.shape): raise ValueError('Gradient directions more than geometry domain') @@ -63,9 +70,7 @@ class FiniteDiff(LinearOperator): else: out = out.as_array() out[:]=0 - - ######################## Direct for 2D ############################### if x_sz == 2: @@ -317,11 +322,22 @@ class FiniteDiff(LinearOperator): return type(x)(out) def range_geometry(self): - '''Returns the range geometry''' + + ''' + + Returns the range_geometry of FiniteDiff + + ''' + return self.gm_range def domain_geometry(self): - '''Returns the domain geometry''' + + ''' + + Returns the domain_geometry of FiniteDiff + + ''' return self.gm_domain diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index b3b23bf..baebc61 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -23,12 +23,29 @@ from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff #%% class Gradient(LinearOperator): + + + r'''Gradient Operator: .. math:: \nabla : X -> Y + + Computes first-order forward/backward differences + on 2D, 3D, 4D ImageData + under Neumann/Periodic boundary conditions + + Example (2D): u\in X, \nabla(u) = [\partial_{y} u, \partial_{x} u] + u^{*}\in Y, \nabla^{*}(u^{*}) = \partial_{y} v1 + \partial_{x} v2 + + 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'] + + + ''' + + CORRELATION_SPACE = "Space" CORRELATION_SPACECHANNEL = "SpaceChannels" - # 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, bnd_cond = 'Neumann', **kwargs): super(Gradient, self).__init__() @@ -72,8 +89,7 @@ class Gradient(LinearOperator): self.bnd_cond = bnd_cond - # Call FiniteDiff operator - + # Call FiniteDiff operator self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond) @@ -114,12 +130,24 @@ class Gradient(LinearOperator): def domain_geometry(self): + + '''Returns domain_geometry of Gradient''' + return self.gm_domain def range_geometry(self): + + '''Returns range_geometry of Gradient''' + return self.gm_range def __rmul__(self, scalar): + + '''Multiplication of Gradient with a scalar + + Returns: ScaledOperator + ''' + return ScaledOperator(self, scalar) ########################################################################### diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py index 54dce08..d8f86a4 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py @@ -19,11 +19,18 @@ from ccpi.optimisation.operators import LinearOperator import scipy.sparse as sp import numpy as np -from ccpi.framework import ImageData class Identity(LinearOperator): + '''Identity: Id: X -> Y, Id(x) = x\in Y + + X : gm_domain + Y : gm_range ( Default: Y = X ) + + ''' + + def __init__(self, gm_domain, gm_range=None): self.gm_domain = gm_domain @@ -34,38 +41,58 @@ class Identity(LinearOperator): super(Identity, self).__init__() def direct(self,x,out=None): + + '''Returns Id(x)''' + if out is None: return x.copy() else: out.fill(x) def adjoint(self,x, out=None): + + '''Returns Id(x)''' + + if out is None: return x.copy() else: out.fill(x) def calculate_norm(self, **kwargs): + + '''Evaluates operator norm of Identity''' + return 1.0 - def domain_geometry(self): + def domain_geometry(self): + + '''Returns domain_geometry of Identity''' + return self.gm_domain def range_geometry(self): + + '''Returns range_geometry of Identity''' + return self.gm_range + + ########################################################################### + ############### For preconditioning ###################################### + ########################################################################### def matrix(self): return sp.eye(np.prod(self.gm_domain.shape)) def sum_abs_row(self): - return self.gm_range.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))) - + return self.gm_range.allocate(1) + def sum_abs_col(self): - return self.gm_domain.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F'))) - + return self.gm_domain.allocate(1) + if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py index b254d08..8514699 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py @@ -17,7 +17,6 @@ # limitations under the License. from ccpi.optimisation.operators import Operator -from ccpi.framework import ImageGeometry import numpy diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py index fbc3c83..bc3312d 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py @@ -17,14 +17,12 @@ # limitations under the License. import numpy from scipy.sparse.linalg import svds -from ccpi.framework import DataContainer -from ccpi.framework import AcquisitionData -from ccpi.framework import VectorData from ccpi.framework import VectorGeometry -from ccpi.framework import AcquisitionGeometry -from numbers import Number from ccpi.optimisation.operators import LinearOperator + class LinearOperatorMatrix(LinearOperator): + '''Matrix wrapped into a LinearOperator''' + def __init__(self,A): self.A = A M_A, N_A = self.A.shape diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py index b9ebef1..c5db47d 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py @@ -19,23 +19,33 @@ from numbers import Number import numpy class ScaledOperator(object): + + '''ScaledOperator + A class to represent the scalar multiplication of an Operator with a scalar. It holds an operator and a scalar. Basically it returns the multiplication of the result of direct and adjoint of the operator with the scalar. For the rest it behaves like the operator it holds. + Args: - operator (Operator): a Operator or LinearOperator - scalar (Number): a scalar multiplier + :param operator (Operator): a Operator or LinearOperator + :param scalar (Number): a scalar multiplier + Example: The scaled operator behaves like the following: - sop = ScaledOperator(operator, scalar) - sop.direct(x) = scalar * operator.direct(x) - sop.adjoint(x) = scalar * operator.adjoint(x) - sop.norm() = operator.norm() - sop.range_geometry() = operator.range_geometry() - sop.domain_geometry() = operator.domain_geometry() + + .. code-block:: python + + sop = ScaledOperator(operator, scalar) + sop.direct(x) = scalar * operator.direct(x) + sop.adjoint(x) = scalar * operator.adjoint(x) + sop.norm() = operator.norm() + sop.range_geometry() = operator.range_geometry() + sop.domain_geometry() = operator.domain_geometry() + ''' + def __init__(self, operator, scalar): super(ScaledOperator, self).__init__() if not isinstance (scalar, Number): diff --git a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py index 1a9f1d8..c1f7ca4 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py @@ -20,6 +20,11 @@ from ccpi.framework import DataContainer class ShrinkageOperator(): + r'''Proximal Operator for .. math:: f(x) = \|\| x \|\|_{1} + + prox_{\tau * f}(x) = x.sign() * \max( |x| - \tau, 0 ) + ''' + def __init__(self): pass diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py index 6b4d822..91d5ca9 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py @@ -22,6 +22,10 @@ from ccpi.framework import ImageData class SparseFiniteDiff(object): + + '''Create Sparse Matrices for the Finite Difference Operator''' + + def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'): super(SparseFiniteDiff, self).__init__() diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index bed56b6..92f8f90 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -19,25 +19,34 @@ from ccpi.optimisation.operators import Gradient, Operator, LinearOperator, ScaledOperator from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer import numpy -from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff +from ccpi.optimisation.operators import FiniteDiff class SymmetrizedGradient(Gradient): - ''' Symmetrized Gradient, denoted by E: V --> W - where V is the Range of the Gradient Operator - and W is the Range of the Symmetrized Gradient. + r'''Symmetrized Gradient Operator: E: V -> W + + V : range of the Gradient Operator + W : range of the Symmetrized Gradient + + Example (2D): + .. math:: + v = (v1, v2), + + Ev = 0.5 * ( \nabla\cdot v + (\nabla\cdot c)^{T} ) + + \begin{matrix} + \partial_{y} v1 & 0.5 * (\partial_{x} v1 + \partial_{y} v2) \\ + 0.5 * (\partial_{x} v1 + \partial_{y} v2) & \partial_{x} v2 + \end{matrix} + | ''' def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) - - ''' - Domain of SymGrad is the Range of Gradient - ''' - + self.gm_domain = self.gm_range self.bnd_cond = bnd_cond @@ -57,6 +66,8 @@ class SymmetrizedGradient(Gradient): def direct(self, x, out=None): + '''Returns E(v)''' + if out is None: tmp = [] @@ -232,12 +243,4 @@ if __name__ == '__main__': print(E1.norm()) print(E2.norm()) - - - - - -# -# -# -# + diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py index b9345da..5f1de30 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -22,6 +22,22 @@ from ccpi.optimisation.operators import LinearOperator class ZeroOperator(LinearOperator): + r'''ZeroOperator: O: X -> Y, maps any element of x\in X into the zero element in Y + O(x) = O_{Y} + + X : gm_domain + Y : gm_range ( Default: Y = X ) + + + Note: + .. math:: + + O^{*}: Y^{*} -> X^{*} (Adjoint) + + < O(x), y > = < x, O^{*}(y) > + + ''' + def __init__(self, gm_domain, gm_range=None): super(ZeroOperator, self).__init__() @@ -33,22 +49,39 @@ class ZeroOperator(LinearOperator): def direct(self,x,out=None): + + '''Returns O(x)''' + + if out is None: return self.gm_range.allocate() else: out.fill(self.gm_range.allocate()) def adjoint(self,x, out=None): + + '''Returns O^{*}(y)''' + if out is None: return self.gm_domain.allocate() else: out.fill(self.gm_domain.allocate()) def calculate_norm(self, **kwargs): - return 0. + + '''Evaluates operator norm of ZeroOperator''' + + return 0 - def domain_geometry(self): + def domain_geometry(self): + + '''Returns domain_geometry of ZeroOperator''' + + return self.gm_domain def range_geometry(self): + + '''Returns domain_geometry of ZeroOperator''' + return self.gm_range \ No newline at end of file diff --git a/docs/source/optimisation.rst b/docs/source/optimisation.rst index 33bf376..eec54e1 100644 --- a/docs/source/optimisation.rst +++ b/docs/source/optimisation.rst @@ -55,6 +55,10 @@ The :code:`Algorithm` provides the infrastructure to continue iteration, to acce :members: .. autoclass:: ccpi.optimisation.algorithms.FISTA :members: +.. autoclass:: ccpi.optimisation.algorithms.PDHG + :members: +.. autoclass:: ccpi.optimisation.algorithms.SIRT + :members: Operator ======== -- cgit v1.2.3