From 743d40f9b98c3b15475fafdee26c9290833f3388 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 9 Oct 2019 12:30:51 +0100 Subject: Finite diff for sirf (#367) * python2 compatibility import future * add staticmethod dot to test LinearOperators * a little more efficient code * skips all tests if module wget is not present * removed sirf import and simplified code --- .../Python/ccpi/framework/BlockDataContainer.py | 42 +++++- Wrappers/Python/ccpi/framework/TestData.py | 7 +- .../Python/ccpi/optimisation/algorithms/CGLS.py | 5 + .../Python/ccpi/optimisation/algorithms/FISTA.py | 5 + .../optimisation/algorithms/GradientDescent.py | 5 + .../Python/ccpi/optimisation/algorithms/PDHG.py | 5 + .../Python/ccpi/optimisation/algorithms/SIRT.py | 5 + .../ccpi/optimisation/functions/BlockFunction.py | 7 +- .../Python/ccpi/optimisation/functions/Function.py | 5 + .../functions/FunctionOperatorComposition.py | 18 ++- .../ccpi/optimisation/functions/IndicatorBox.py | 7 +- .../ccpi/optimisation/functions/KullbackLeibler.py | 5 + .../Python/ccpi/optimisation/functions/L1Norm.py | 7 +- .../ccpi/optimisation/functions/L2NormSquared.py | 5 + .../ccpi/optimisation/functions/MixedL21Norm.py | 45 +++++-- .../Python/ccpi/optimisation/functions/Norm2Sq.py | 28 +++- .../ccpi/optimisation/functions/ScaledFunction.py | 5 + .../ccpi/optimisation/functions/ZeroFunction.py | 5 + .../ccpi/optimisation/operators/BlockOperator.py | 42 +++++- .../optimisation/operators/BlockScaledOperator.py | 7 +- .../operators/FiniteDifferenceOperator.py | 25 +++- .../optimisation/operators/GradientOperator.py | 5 + .../optimisation/operators/IdentityOperator.py | 7 +- .../ccpi/optimisation/operators/LinearOperator.py | 52 ++++++- .../optimisation/operators/LinearOperatorMatrix.py | 6 + .../Python/ccpi/optimisation/operators/Operator.py | 6 + .../optimisation/operators/ShrinkageOperator.py | 7 +- .../optimisation/operators/SparseFiniteDiff.py | 7 +- .../operators/SymmetrizedGradientOperator.py | 11 +- .../ccpi/optimisation/operators/ZeroOperator.py | 7 +- Wrappers/Python/test/test_NexusReader.py | 149 +++++++++++---------- 31 files changed, 417 insertions(+), 125 deletions(-) diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index b5116e5..8247f24 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -26,7 +26,8 @@ import functools from ccpi.framework import DataContainer #from ccpi.framework import AcquisitionData, ImageData #from ccpi.optimisation.operators import Operator, LinearOperator - + + class BlockDataContainer(object): '''Class to hold DataContainers as column vector @@ -102,7 +103,10 @@ class BlockDataContainer(object): raise ValueError('List/ numpy array can only contain numbers {}'\ .format(type(ot))) return len(self.containers) == len(other) - elif issubclass(other.__class__, DataContainer): + elif isinstance(other, BlockDataContainer): + return len(self.containers) == len(other.containers) + else: + # this should work for other as DataContainers and children ret = True for i, el in enumerate(self.containers): if isinstance(el, BlockDataContainer): @@ -110,9 +114,9 @@ class BlockDataContainer(object): else: a = el.shape == other.shape ret = ret and a + # probably will raise return ret - #return self.get_item(0).shape == other.shape - return len(self.containers) == len(other.containers) + def get_item(self, row): if row > self.shape[0]: @@ -180,7 +184,7 @@ class BlockDataContainer(object): if not self.is_compatible(other): raise ValueError('Incompatible for divide') out = kwargs.get('out', None) - if isinstance(other, Number) or issubclass(other.__class__, DataContainer): + if isinstance(other, Number): # try to do algebra with one DataContainer. Will raise error if not compatible kw = kwargs.copy() res = [] @@ -240,8 +244,32 @@ class BlockDataContainer(object): return type(self)(*res, shape=self.shape) return type(self)(*[ operation(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape) else: - raise ValueError('Incompatible type {}'.format(type(other))) - + # try to do algebra with one DataContainer. Will raise error if not compatible + kw = kwargs.copy() + res = [] + for i,el in enumerate(self.containers): + if operation == BlockDataContainer.ADD: + op = el.add + elif operation == BlockDataContainer.SUBTRACT: + op = el.subtract + elif operation == BlockDataContainer.MULTIPLY: + op = el.multiply + elif operation == BlockDataContainer.DIVIDE: + op = el.divide + elif operation == BlockDataContainer.POWER: + op = el.power + else: + raise ValueError('Unsupported operation', operation) + if out is not None: + kw['out'] = out.get_item(i) + op(other, *args, **kw) + else: + res.append(op(other, *args, **kw)) + if out is not None: + return + else: + return type(self)(*res, shape=self.shape) + def power(self, other, *args, **kwargs): if not self.is_compatible(other): diff --git a/Wrappers/Python/ccpi/framework/TestData.py b/Wrappers/Python/ccpi/framework/TestData.py index 2f4c685..74d37be 100755 --- a/Wrappers/Python/ccpi/framework/TestData.py +++ b/Wrappers/Python/ccpi/framework/TestData.py @@ -15,6 +15,11 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.framework import ImageData, ImageGeometry, DataContainer import numpy import numpy as np @@ -334,4 +339,4 @@ class TestData(object): if clip: out = np.clip(out, low_clip, 1.0) - return out \ No newline at end of file + return out diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py index c6c1d4c..d2e5b29 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -17,6 +17,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.algorithms import Algorithm import numpy diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py index e23116b..5d79b67 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -17,6 +17,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.algorithms import Algorithm from ccpi.optimisation.functions import ZeroFunction import numpy diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py index d060690..f79651a 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py @@ -20,6 +20,11 @@ # #========================================================================= +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.algorithms import Algorithm class GradientDescent(Algorithm): diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 0968872..7bc4e11 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -17,6 +17,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.algorithms import Algorithm diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index 2b49ab0..8feef87 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -20,6 +20,11 @@ # #========================================================================= +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.algorithms import Algorithm class SIRT(Algorithm): diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py index a6ac66c..ee3ad78 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -20,6 +20,11 @@ # #========================================================================= +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.functions import Function from ccpi.framework import BlockDataContainer from numbers import Number @@ -232,4 +237,4 @@ if __name__ == '__main__': - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/functions/Function.py b/Wrappers/Python/ccpi/optimisation/functions/Function.py index 7156995..48c6d30 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/Function.py +++ b/Wrappers/Python/ccpi/optimisation/functions/Function.py @@ -20,6 +20,11 @@ # #========================================================================= +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + 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 58d4f27..4162134 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -39,8 +39,11 @@ class FunctionOperatorComposition(Function): self.function = function self.operator = operator - self.L = function.L * operator.norm()**2 - + try: + self.L = function.L * operator.norm()**2 + except Error as er: + self.L = None + warnings.warn("Lipschitz constant was not calculated") def __call__(self, x): @@ -56,12 +59,13 @@ class FunctionOperatorComposition(Function): ''' + tmp = self.operator.range_geometry().allocate() + self.operator.direct(x, out=tmp) + self.function.gradient(tmp, out=tmp) if out is None: - return self.operator.adjoint(self.function.gradient(self.operator.direct(x))) + #return self.operator.adjoint(self.function.gradient(self.operator.direct(x))) + return self.operator.adjoint(tmp) else: - tmp = self.operator.range_geometry().allocate() - self.operator.direct(x, out=tmp) - self.function.gradient(tmp, out=tmp) self.operator.adjoint(tmp, out=out) @@ -122,4 +126,4 @@ if __name__ == '__main__': - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py index 51d08d1..fd34d96 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py +++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py @@ -16,6 +16,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.functions import Function import numpy @@ -145,4 +150,4 @@ if __name__ == '__main__': - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index f88c339..d71f597 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -17,6 +17,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + import numpy from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py index cc4bef8..09e550e 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -17,6 +17,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction from ccpi.optimisation.operators import ShrinkageOperator @@ -172,4 +177,4 @@ if __name__ == '__main__': - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index a625f07..92e0116 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -18,6 +18,11 @@ # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py index 8cbed67..378cbda 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py @@ -18,10 +18,16 @@ # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.functions import Function, ScaledFunction from ccpi.framework import BlockDataContainer import functools +import numpy class MixedL21Norm(Function): @@ -45,11 +51,11 @@ class MixedL21Norm(Function): ''' if not isinstance(x, BlockDataContainer): raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) - - tmp = [ el**2 for el in x.containers ] - res = sum(tmp).sqrt().sum() + tmp = x.get_item(0) * 0 + for el in x.containers: + tmp += el.power(2.) + return tmp.sqrt().sum() - return res def gradient(self, x, out=None): return ValueError('Not Differentiable') @@ -84,16 +90,29 @@ class MixedL21Norm(Function): if out is None: - tmp = [ el*el for el in x.containers] - res = sum(tmp).sqrt().maximum(1.0) - frac = [el/res for el in x.containers] - return BlockDataContainer(*frac) + # tmp = [ el*el for el in x.containers] + # res = sum(tmp).sqrt().maximum(1.0) + # frac = [el/res for el in x.containers] + # return BlockDataContainer(*frac) + tmp = x.get_item(0) * 0 + for el in x.containers: + tmp += el.power(2.) + tmp.sqrt(out=tmp) + tmp.maximum(1.0, out=tmp) + frac = [ el.divide(tmp) for el in x.containers ] + return BlockDataContainer(*frac) + 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) + if False: + res = res1.sqrt().maximum(1.0) + x.divide(res, out=out) + else: + res1.sqrt(out=res1) + res1.maximum(1.0, out=res1) + x.divide(res1, out=out) def __rmul__(self, scalar): @@ -106,6 +125,12 @@ class MixedL21Norm(Function): return ScaledFunction(self, scalar) +def sqrt_maximum(x, a): + y = numpy.sqrt(x) + if y >= a: + return y + else: + return a # if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py index 0da6e50..01c4f38 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py +++ b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py @@ -17,7 +17,12 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals +from ccpi.optimisation.operators import LinearOperator from ccpi.optimisation.functions import Function import warnings @@ -50,8 +55,13 @@ class Norm2Sq(Function): try: self.L = 2.0*self.c*(self.A.norm()**2) except AttributeError as ae: - warnings.warn('{} could not calculate Lipschitz Constant. {}'.format( + if self.A.is_linear(): + Anorm = LinearOperator.PowerMethod(self.A, 10)[0] + self.L = 2.0 * self.c * (Anorm*Anorm) + else: + warnings.warn('{} could not calculate Lipschitz Constant. {}'.format( self.__class__.__name__, ae)) + except NotImplementedError as noe: warnings.warn('{} could not calculate Lipschitz Constant. {}'.format( self.__class__.__name__, noe)) @@ -65,22 +75,28 @@ class Norm2Sq(Function): # return self.c*( ( (self.A.direct(x)-self.b)**2).sum() ) #else: y = self.A.direct(x) - y.__isub__(self.b) + y.subtract(self.b, out=y) #y.__imul__(y) #return y.sum() * self.c try: + if self.c == 1: + return y.squared_norm() return y.squared_norm() * self.c except AttributeError as ae: - # added for compatibility with SIRF - return (y.norm()**2) * self.c + # added for compatibility with SIRF + warnings.warn('squared_norm method not found! Proceeding with norm.') + yn = y.norm() + if self.c == 1: + return yn * yn + return (yn * yn) * self.c def gradient(self, x, out=None): if out is not None: #return 2.0*self.c*self.A.adjoint( self.A.direct(x) - self.b ) self.A.direct(x, out=self.range_tmp) - self.range_tmp -= self.b + self.range_tmp.subtract(self.b , out=self.range_tmp) self.A.adjoint(self.range_tmp, out=out) #self.direct_placehold.multiply(2.0*self.c, out=out) - out *= (self.c * 2.0) + out.multiply (self.c * 2.0, out=out) else: return (2.0*self.c)*self.A.adjoint(self.A.direct(x) - self.b) diff --git a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py index 3e689e6..a123e8d 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py @@ -17,6 +17,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + 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 ca52f31..19db668 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py @@ -18,6 +18,11 @@ # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.functions import Function class ZeroFunction(Function): diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index e3a02ec..23cb799 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -17,11 +17,22 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + import numpy import functools from ccpi.framework import ImageData, BlockDataContainer, DataContainer -from ccpi.optimisation.operators import Operator +from ccpi.optimisation.operators import Operator, LinearOperator from ccpi.framework import BlockGeometry +try: + from sirf import SIRF + from sirf.SIRF import DataContainer as SIRFDataContainer + has_sirf = True +except ImportError as ie: + has_sirf = False class BlockOperator(Operator): r'''A Block matrix containing Operators @@ -115,7 +126,23 @@ class BlockOperator(Operator): return self.operators[index] def norm(self, **kwargs): - norm = [op.norm(**kwargs)**2 for op in self.operators] + '''Returns the norm of the BlockOperator + + if the operator in the block do not have method norm defined, i.e. they are SIRF + AcquisitionModel's we use PowerMethod if applicable, otherwise we raise an Error + ''' + norm = [] + for op in self.operators: + if hasattr(op, 'norm'): + norm.append(op.norm(**kwargs) ** 2.) + else: + # use Power method + if op.is_linear(): + norm.append( + LinearOperator.PowerMethod(op, 20)[0] + ) + else: + raise TypeError('Operator {} does not have a norm method and is not linear'.format(op)) return numpy.sqrt(sum(norm)) def direct(self, x, out=None): @@ -188,7 +215,8 @@ class BlockOperator(Operator): prod += self.get_item(row, col).adjoint(x_b.get_item(row)) res.append(prod) if self.shape[1]==1: - return ImageData(*res) + # the output is a single DataContainer, so we can take it out + return res[0] else: return BlockDataContainer(*res, shape=shape) else: @@ -196,7 +224,8 @@ class BlockOperator(Operator): for col in range(self.shape[1]): for row in range(self.shape[0]): if row == 0: - if issubclass(out.__class__, DataContainer): + if issubclass(out.__class__, DataContainer) or \ + ( has_sirf and issubclass(out.__class__, SIRFDataContainer) ): self.get_item(row, col).adjoint( x_b.get_item(row), out=out) @@ -206,7 +235,8 @@ class BlockOperator(Operator): x_b.get_item(row), out=out.get_item(col)) else: - if issubclass(out.__class__, DataContainer): + if issubclass(out.__class__, DataContainer) or \ + ( has_sirf and issubclass(out.__class__, SIRFDataContainer) ): out += self.get_item(row,col).adjoint( x_b.get_item(row)) else: @@ -423,4 +453,4 @@ if __name__ == '__main__': - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py index c23c23a..eeecee9 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py @@ -18,6 +18,11 @@ # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from numbers import Number import numpy from ccpi.optimisation.operators import ScaledOperator @@ -82,4 +87,4 @@ class BlockScaledOperator(ScaledOperator): @property def T(self): '''Return the transposed of self''' - return type(self)(self.operator.T, self.scalar) \ No newline at end of file + return type(self)(self.operator.T, self.scalar) diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py index 9b5ae24..3cc4309 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -15,6 +15,11 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.operators import LinearOperator import numpy as np @@ -66,8 +71,11 @@ class FiniteDiff(LinearOperator): x_sz = len(x.shape) if out is None: - out = np.zeros_like(x_asarr) + res = x * 0 + #out = np.zeros_like(x_asarr) + out = res.as_array() else: + res = out out = out.as_array() out[:]=0 @@ -180,7 +188,9 @@ class FiniteDiff(LinearOperator): raise NotImplementedError # res = out #/self.voxel_size - return type(x)(out) + #return type(x)(out) + res.fill(out) + return res def adjoint(self, x, out=None): @@ -189,8 +199,11 @@ class FiniteDiff(LinearOperator): x_sz = len(x.shape) if out is None: - out = np.zeros_like(x_asarr) + #out = np.zeros_like(x_asarr) + res = x * 0 + out = res.as_array() else: + res = out out = out.as_array() out[:]=0 @@ -319,7 +332,9 @@ class FiniteDiff(LinearOperator): raise NotImplementedError out *= -1 #/self.voxel_size - return type(x)(out) + #return type(x)(out) + res.fill(out) + return res def range_geometry(self): @@ -389,4 +404,4 @@ if __name__ == '__main__': - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index baebc61..3c32a93 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -15,6 +15,11 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer import numpy diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py index d8f86a4..e95234b 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py @@ -16,6 +16,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.operators import LinearOperator import scipy.sparse as sp import numpy as np @@ -113,4 +118,4 @@ if __name__ == '__main__': - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py index 8514699..f4d97b8 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py @@ -16,6 +16,11 @@ # See the License for the specific language governing permissions and # limitations under the License. +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.optimisation.operators import Operator import numpy @@ -39,7 +44,7 @@ class LinearOperator(Operator): # Initialise random if x_init is None: - x0 = operator.domain_geometry().allocate(type(operator.domain_geometry()).RANDOM_INT) + x0 = operator.domain_geometry().allocate('random') else: x0 = x_init.copy() @@ -51,7 +56,11 @@ class LinearOperator(Operator): operator.direct(x0,out=y_tmp) operator.adjoint(y_tmp,out=x1) x1norm = x1.norm() - s[it] = x1.dot(x0) / x0.squared_norm() + if hasattr(x0, 'squared_norm'): + s[it] = x1.dot(x0) / x0.squared_norm() + else: + x0norm = x0.norm() + s[it] = x1.dot(x0) / (x0norm * x0norm) x1.multiply((1.0/x1norm), out=x0) return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 @@ -62,4 +71,41 @@ class LinearOperator(Operator): s1, sall, svec = LinearOperator.PowerMethod(self, iterations, x_init=x0) return s1 - + @staticmethod + def dot_test(operator, domain_init=None, range_init=None, verbose=False): + '''Does a dot linearity test on the operator + + Evaluates if the following equivalence holds + + :math: .. + + Ax\times y = y \times A^Tx + + :param operator: operator to test + :param range_init: optional initialisation container in the operator range + :param domain_init: optional initialisation container in the operator domain + :returns: boolean, True if the test is passed. + ''' + if range_init is None: + y = operator.range_geometry().allocate('random_int') + else: + y = range_init + if domain_init is None: + x = operator.domain_geometry().allocate('random_int') + else: + x = domain_init + + fx = operator.direct(x) + by = operator.adjoint(y) + a = fx.dot(y) + b = by.dot(x) + if verbose: + print ('Left hand side {}, \nRight hand side {}'.format(a, b)) + try: + numpy.testing.assert_almost_equal(abs((a-b)/a), 0, decimal=4) + return True + except AssertionError as ae: + print (ae) + return False + + diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py index bc3312d..7d18ea1 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py @@ -15,6 +15,12 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + import numpy from scipy.sparse.linalg import svds from ccpi.framework import VectorGeometry diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py index 2678bf2..87059e6 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -15,6 +15,12 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + + from ccpi.optimisation.operators.ScaledOperator import ScaledOperator class Operator(object): diff --git a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py index c1f7ca4..9239d90 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py @@ -15,6 +15,11 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + from ccpi.framework import DataContainer @@ -31,4 +36,4 @@ class ShrinkageOperator(): def __call__(self, x, tau, out=None): return x.sign() * (x.abs() - tau).maximum(0) - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py index 91d5ca9..698a993 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py @@ -15,6 +15,11 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + import scipy.sparse as sp import numpy as np @@ -155,4 +160,4 @@ if __name__ == '__main__': u_per_sp_adjoint3D = sFD_per3D.adjoint(arr3D) np.testing.assert_array_almost_equal(u_per_adjoint3D.as_array(), u_per_sp_adjoint3D.as_array(), decimal=4) - \ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index 92f8f90..8d14cf8 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -15,9 +15,16 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals -from ccpi.optimisation.operators import Gradient, Operator, LinearOperator, ScaledOperator -from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer + +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 diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py index 5f1de30..c37e15e 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -15,6 +15,11 @@ # 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 __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + import numpy as np from ccpi.framework import ImageData @@ -84,4 +89,4 @@ class ZeroOperator(LinearOperator): '''Returns domain_geometry of ZeroOperator''' - return self.gm_range \ No newline at end of file + return self.gm_range diff --git a/Wrappers/Python/test/test_NexusReader.py b/Wrappers/Python/test/test_NexusReader.py index 71a05c2..992ce4f 100755 --- a/Wrappers/Python/test/test_NexusReader.py +++ b/Wrappers/Python/test/test_NexusReader.py @@ -17,7 +17,11 @@ # limitations under the License. import unittest -import wget +has_wget = True +try: + import wget +except ImportError as ie: + has_wget = False import os from ccpi.io.reader import NexusReader import numpy @@ -26,80 +30,85 @@ import numpy class TestNexusReader(unittest.TestCase): def setUp(self): - wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') - self.filename = '24737_fd.nxs' + if has_wget: + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' def tearDown(self): - os.remove(self.filename) + if has_wget: + os.remove(self.filename) def testAll(self): - # def testGetDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct") - - # def testGetProjectionDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct") - - # def testLoadProjectionWithoutDimensions(self): - nr = NexusReader(self.filename) - projections = nr.load_projection() - self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") - - # def testLoadProjectionWithDimensions(self): - nr = NexusReader(self.filename) - projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) - self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") - - # def testLoadProjectionCompareSingle(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) - - # def testLoadProjectionCompareMulti(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) - - # def testLoadProjectionCompareRandom(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) - numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) - - # def testLoadProjectionCompareFull(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) - - # def testLoadFlatCompareFull(self): - nr = NexusReader(self.filename) - flats_full = nr.load_flat() - flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) - - # def testLoadDarkCompareFull(self): - nr = NexusReader(self.filename) - darks_full = nr.load_dark() - darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) - - # def testProjectionAngles(self): - nr = NexusReader(self.filename) - angles = nr.get_projection_angles() - self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") - - # def test_get_acquisition_data_subset(self): - nr = NexusReader(self.filename) - key = nr.get_image_keys() - sl = nr.get_acquisition_data_subset(0,10) - data = nr.get_acquisition_data().subset(['vertical','horizontal']) - - self.assertTrue(sl.shape , (10,data.shape[1])) + if has_wget: + # def testGetDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct") + # def testGetProjectionDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct") + + # def testLoadProjectionWithoutDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection() + self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") + + # def testLoadProjectionWithDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") + + # def testLoadProjectionCompareSingle(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) + + # def testLoadProjectionCompareMulti(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) + + # def testLoadProjectionCompareRandom(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) + numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) + + # def testLoadProjectionCompareFull(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) + + # def testLoadFlatCompareFull(self): + nr = NexusReader(self.filename) + flats_full = nr.load_flat() + flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) + + # def testLoadDarkCompareFull(self): + nr = NexusReader(self.filename) + darks_full = nr.load_dark() + darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) + + # def testProjectionAngles(self): + nr = NexusReader(self.filename) + angles = nr.get_projection_angles() + self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") + + # def test_get_acquisition_data_subset(self): + nr = NexusReader(self.filename) + key = nr.get_image_keys() + sl = nr.get_acquisition_data_subset(0,10) + data = nr.get_acquisition_data().subset(['vertical','horizontal']) + + self.assertTrue(sl.shape , (10,data.shape[1])) + else: + # skips all tests if module wget is not present + self.assertFalse(has_wget) if __name__ == '__main__': -- cgit v1.2.3