summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorVaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com>2019-07-05 10:46:21 +0100
committerEdoardo Pasca <edo.paskino@gmail.com>2019-07-05 10:46:21 +0100
commitab405e3d4714611ac3964f5a563d809977d8daa6 (patch)
tree45815071f1a8200315ebec5606872bad643d1d75
parent5fa800621a28380d1731f7d4fa8cf426e0c04b70 (diff)
downloadframework-ab405e3d4714611ac3964f5a563d809977d8daa6.tar.gz
framework-ab405e3d4714611ac3964f5a563d809977d8daa6.tar.bz2
framework-ab405e3d4714611ac3964f5a563d809977d8daa6.tar.xz
framework-ab405e3d4714611ac3964f5a563d809977d8daa6.zip
Fix documentation (#343)
* fix documentation algs * fix docum functions dir * fix docum operator dir * docstring title * updated math in docstring
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/Algorithm.py34
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/CGLS.py59
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/FISTA.py72
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py43
-rw-r--r--Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py81
-rw-r--r--Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py57
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py78
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/Function.py36
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py56
-rwxr-xr-xWrappers/Python/ccpi/optimisation/functions/IndicatorBox.py38
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py83
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/L1Norm.py76
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py83
-rwxr-xr-xWrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py83
-rwxr-xr-xWrappers/Python/ccpi/optimisation/functions/Norm2Sq.py39
-rwxr-xr-xWrappers/Python/ccpi/optimisation/functions/ScaledFunction.py32
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py72
-rwxr-xr-xWrappers/Python/ccpi/optimisation/operators/BlockOperator.py45
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py71
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py48
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py40
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py39
-rwxr-xr-xWrappers/Python/ccpi/optimisation/operators/LinearOperator.py1
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py8
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py26
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py5
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py4
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py39
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py37
-rw-r--r--docs/source/optimisation.rst4
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 <x, z>, 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 <x, x^{*}> 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 <x, x^{*}>
+
+ 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
========