From adf4163c145e6ddc16899a92a06c3282f144d88c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 21 Feb 2020 16:41:57 +0000 Subject: Operator composition (#493) * CompositionOperator and some refactoring * Added SumOperator and CompositionOperator added domain_geometry and optional range_geometry as parameter of Operator. These are saved in _domain_geometry and _range_geometry. Updated all operator and tests to take notice of this change. * fighting with Composition * fix direct and adjoint for CompositeOperator * fix composition Operator * remove target OUTPUT to trigger build * added unit tests * add test for ZeroOperator * fixes * add numba * removed hard coded path * removed comments * fix direct/adjoint with out parameter * removed __rmul__ * use calculate_norm inherited from LinearOperator * removed hard coded version tag * removed cached string version * use add_custom_target instead add_custom_command --- Wrappers/Python/CMakeLists.txt | 54 +-- .../Python/ccpi/framework/BlockDataContainer.py | 2 +- .../ccpi/optimisation/functions/KullbackLeibler.py | 140 +++++-- .../operators/FiniteDifferenceOperator.py | 25 +- .../optimisation/operators/GradientOperator.py | 59 +-- .../optimisation/operators/IdentityOperator.py | 24 +- .../optimisation/operators/LinearOperatorMatrix.py | 17 +- .../Python/ccpi/optimisation/operators/Operator.py | 433 ++++++++++++++++++++- .../optimisation/operators/SparseFiniteDiff.py | 12 +- .../operators/SymmetrizedGradientOperator.py | 42 +- .../ccpi/optimisation/operators/ZeroOperator.py | 33 +- .../Python/ccpi/optimisation/operators/__init__.py | 9 +- Wrappers/Python/test/test_BlockOperator.py | 4 +- Wrappers/Python/test/test_Operator.py | 264 +++++++++++++ 14 files changed, 932 insertions(+), 186 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt index 79bc912..47c6406 100644 --- a/Wrappers/Python/CMakeLists.txt +++ b/Wrappers/Python/CMakeLists.txt @@ -3,7 +3,7 @@ option (BUILD_PYTHON_WRAPPER "Build Python Wrapper" ON) if (BUILD_PYTHON_WRAPPER) find_package(PythonInterp REQUIRED) - set(PYTHON_DEST_DIR "" CACHE PATH "Directory of the Python wrappers") + #set(PYTHON_DEST_DIR "" CACHE PATH "Directory of the Python wrappers") if (PYTHON_DEST_DIR) set(PYTHON_DEST "${PYTHON_DEST_DIR}") else() @@ -34,60 +34,68 @@ if (BUILD_PYTHON_WRAPPER) if (PYTHONINTERP_FOUND) message("Python found " ${PYTHON_EXECUTABLE}) set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup-cil.py.in") - set(SETUP_PY "${CMAKE_CURRENT_BINARY_DIR}/setup.py") + set(SETUP_PY "${CMAKE_CURRENT_SOURCE_DIR}/setup.py") #set(DEPS "${CMAKE_CURRENT_SOURCE_DIR}/module/__init__.py") - set (DEPS "${CMAKE_BINARY_DIR}/") - set(OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp") - - #configure_file(${SETUP_PY_IN} ${SETUP_PY}) - + # set (DEPS "${CMAKE_CURRENT_SOURCE_DIR}/ccpi") + set(OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/timestamp") + file(GLOB_RECURSE DEPS ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/*.py ) + + # add to add_custom_command DEPENDS the list of python files of the project. + # as a hack I remove ${OUTPUT}. This should trigger the new build. + file( REMOVE ${OUTPUT} ) + #file( REMOVE_RECURSE ${CMAKE_CURRENT_BINARY_DIR}/build/ ) + #message(STATUS "We should have removed the build directory now") + #file( COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + #file( COPY ${CMAKE_CURRENT_SOURCE_DIR}/test DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + #file( COPY ${CMAKE_CURRENT_SOURCE_DIR}/setup.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + if (CONDA_BUILD) - add_custom_command(OUTPUT ${OUTPUT} - #COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} + add_custom_target(pythonsetup ALL COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION} ${PYTHON_EXECUTABLE} ${SETUP_PY} -vv install - #echo "EDO" - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT} DEPENDS cilacc) else() if (WIN32) - add_custom_command(OUTPUT ${OUTPUT} - COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} + add_custom_target(pythonsetup ALL COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION} PREFIX=${CMAKE_SOURCE_DIR}/src/ LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/include LIBRARY_LIB=${CMAKE_BINARY_DIR}/${CMAKE_BUILD_TYPE} ${PYTHON_EXECUTABLE} ${SETUP_PY} build_py + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT} DEPENDS cilacc) else() - add_custom_command(OUTPUT ${OUTPUT} - COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env PREFIX=${CMAKE_SOURCE_DIR}/src/ + add_custom_target(pythonsetup ALL + COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION} + PREFIX=${CMAKE_SOURCE_DIR}/src/ LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/include LIBRARY_LIB=${CMAKE_BINARY_DIR}/ - ${PYTHON_EXECUTABLE} ${SETUP_PY} build_py --build-lib ${CMAKE_CURRENT_BINARY_DIR}/build/lib/ + ${PYTHON_EXECUTABLE} ${SETUP_PY} build_py --verbose --build-lib=${CMAKE_CURRENT_BINARY_DIR}/build/lib + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT} - DEPENDS cilacc + DEPENDS cilacc ) endif() #set (PYTHON_DEST ${CMAKE_INSTALL_PREFIX}/python/) install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/build/lib/ccpi - - DESTINATION ${PYTHON_DEST}) + DESTINATION ${PYTHON_DEST} ) install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/data/ DESTINATION ${CMAKE_INSTALL_PREFIX}/share/ccpi) #file(TOUCH ${PYTHON_DEST}/edo/__init__.py) endif() - add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT}) - + #add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT}) + add_custom_target(PythonWrapper ALL DEPENDS pythonsetup) + endif() endif() diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index a0d139b..310132a 100644 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -261,7 +261,7 @@ class BlockDataContainer(object): kw['out'] = out.get_item(i) if operation == BlockDataContainer.AXPBY: kw['y'] = ot - el.axpby(kw['a'], kw['b'], kw['y'], kw['out'], kw['dtype'], kw['num_threads']) + el.axpby(kw['a'], kw['b'], kw['y'], kw['out'], dtype=kw['dtype'], num_threads=kw['num_threads']) else: op(ot, *args, **kw) else: diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 39bd983..8fb7ab7 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -26,6 +26,66 @@ from ccpi.optimisation.functions import Function import functools import scipy.special +try: + import numba + from numba import jit, prange + has_numba = True + '''Some parallelisation of KL calls''' + @jit(nopython=True) + def kl_proximal(x,b, bnoise, tau, out, eta): + for i in prange(x.size): + out.flat[i] = 0.5 * ( + ( x.flat[i] + eta.flat[i] - bnoise.flat[i] - tau ) +\ + numpy.sqrt( (x.flat[i] + eta.flat[i] + bnoise.flat[i] - tau)**2. + \ + (4. * tau * b.flat[i]) + ) + ) + @jit(nopython=True) + def kl_proximal_conjugate(x, b, bnoise, tau, out): + #z = x + tau * self.bnoise + #return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) + + for i in prange(x.size): + z = x.flat[i] + ( tau * bnoise.flat[i] ) + out.flat[i] = 0.5 * ( + (z + 1) - numpy.sqrt((z-1)*(z-1) + 4 * tau * b.flat[i]) + ) + @jit(nopython=True) + def kl_gradient(x, b, bnoise, out, eta): + for i in prange(x.size): + out.flat[i] = 1 - b.flat[i]/(x.flat[i] + bnoise.flat[i] + eta.flat[i]) + + @jit(nopython=True) + def kl_div(x, y, eta): + accumulator = 0. + for i in prange(x.size): + X = x.flat[i] + Y = y.flat[i] + eta.flat[i] + if X > 0 and Y > 0: + # out.flat[i] = X * numpy.log(X/Y) - X + Y + accumulator += X * numpy.log(X/Y) - X + Y + elif X == 0 and Y >= 0: + # out.flat[i] = Y + accumulator += Y + else: + # out.flat[i] = numpy.inf + return numpy.inf + return accumulator + + # force a jit + x = numpy.asarray(numpy.random.random((10,10)), dtype=numpy.float32) + b = numpy.asarray(numpy.random.random((10,10)), dtype=numpy.float32) + bnoise = numpy.zeros_like(x) + eta = numpy.zeros_like(x) + out = numpy.empty_like(x) + tau = 1. + kl_div(b,x,eta) + kl_gradient(x,b,bnoise,out, eta) + kl_proximal(x,b, bnoise, tau, out, eta) + kl_proximal_conjugate(x,b, bnoise, tau, out) + +except ImportError as ie: + has_numba = False class KullbackLeibler(Function): @@ -84,11 +144,14 @@ class KullbackLeibler(Function): r"""Returns the value of the KullbackLeibler function at :math:`(b, x + \eta)`. To avoid infinity values, we consider only pixels/voxels for :math:`x+\eta\geq0`. """ - - tmp_sum = (x + self.eta).as_array() - ind = tmp_sum >= 0 - tmp = scipy.special.kl_div(self.b.as_array()[ind], tmp_sum[ind]) - return numpy.sum(tmp) + if has_numba: + # tmp = numpy.empty_like(x.as_array()) + return kl_div(self.b.as_array(), x.as_array(), self.eta.as_array()) + else: + tmp_sum = (x + self.eta).as_array() + ind = tmp_sum >= 0 + tmp = scipy.special.kl_div(self.b.as_array()[ind], tmp_sum[ind]) + return numpy.sum(tmp) def log(self, datacontainer): '''calculates the in-place log of the datacontainer''' @@ -106,15 +169,26 @@ class KullbackLeibler(Function): We require the :math:`x+\eta>0` otherwise we have inf values. """ - - tmp_sum_array = (x + self.eta).as_array() - if out is None: - tmp_out = x.geometry.allocate() - tmp_out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0] - return tmp_out - else: - x.add(self.eta, out=out) - out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0] + if has_numba: + if out is None: + out = (x * 0.) + out_np = out.as_array() + kl_gradient(x.as_array(), self.b.as_array(), self.bnoise.as_array(), out_np, self.eta.as_array()) + # out.fill(out_np) + return out + else: + out_np = out.as_array() + kl_gradient(x.as_array(), self.b.as_array(), self.bnoise.as_array(), out_np, self.eta.as_array()) + # out.fill(out_np) + else: + tmp_sum_array = (x + self.eta).as_array() + if out is None: + tmp_out = x.geometry.allocate() + tmp_out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0] + return tmp_out + else: + x.add(self.eta, out=out) + out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0] def convex_conjugate(self, x): @@ -142,18 +216,30 @@ class KullbackLeibler(Function): where :math:`z = x + \tau \eta` """ - - if out is None: - return 0.5 *( (x - self.eta - tau) + ( (x + self.eta - tau)**2 + 4*tau*self.b ) .sqrt() ) - else: - x.add(self.eta, out=out) - out -= tau - out *= out - out.add(self.b * (4 * tau), out=out) - out.sqrt(out=out) - out.subtract(tau, out = out) - out.add(x, out=out) - out *= 0.5 + if has_numba: + if out is None: + out = (x * 0.) + # out_np = numpy.empty_like(out.as_array(), dtype=numpy.float64) + out_np = out.as_array() + kl_proximal(x.as_array(), self.b.as_array(), self.bnoise.as_array(), tau, out_np, self.eta.as_array()) + out.fill(out_np) + return out + else: + out_np = out.as_array() + kl_proximal(x.as_array(), self.b.as_array(), self.bnoise.as_array(), tau, out_np, self.eta.as_array()) + # out.fill(out_np) + else: + if out is None: + return 0.5 *( (x - self.eta - tau) + ( (x + self.eta - tau)**2 + 4*tau*self.b ) .sqrt() ) + else: + x.add(self.eta, out=out) + out -= tau + out *= out + out.add(self.b * (4 * tau), out=out) + out.sqrt(out=out) + out.subtract(tau, out = out) + out.add(x, out=out) + out *= 0.5 def proximal_conjugate(self, x, tau, out=None): @@ -174,7 +260,7 @@ class KullbackLeibler(Function): self.b.multiply(4*tau, out=out) - out.add((tmp)**2, out=out) + out.add(tmp.power(2), out=out) out.sqrt(out=out) out *= -1 tmp += 2 diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py index 1df7745..f8aef7d 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -56,8 +56,7 @@ class FiniteDiff(LinearOperator): :type bnd_cond: str, default :code:`Neumann` ''' - super(FiniteDiff, self).__init__() - + self.gm_domain = gm_domain self.gm_range = gm_range @@ -75,6 +74,8 @@ class FiniteDiff(LinearOperator): #self.voxel_size = kwargs.get('voxel_size',1) # this wrongly assumes a homogeneous voxel size # self.voxel_size = self.gm_domain.voxel_size_x + super(FiniteDiff, self).__init__(domain_geometry=gm_domain, + range_geometry=self.gm_range) def direct(self, x, out=None): @@ -350,24 +351,24 @@ class FiniteDiff(LinearOperator): #else: # out.fill(outa) - def range_geometry(self): + # def range_geometry(self): - ''' + # ''' - Returns the range_geometry of FiniteDiff + # Returns the range_geometry of FiniteDiff - ''' + # ''' - return self.gm_range + # return self.gm_range - def domain_geometry(self): + # def domain_geometry(self): - ''' + # ''' - Returns the domain_geometry of FiniteDiff + # Returns the domain_geometry of FiniteDiff - ''' - return self.gm_domain + # ''' + # return self.gm_domain if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index a5feca3..fd51873 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -75,24 +75,29 @@ class Gradient(LinearOperator): CORRELATION_SPACE = CORRELATION_SPACE CORRELATION_SPACECHANNEL = CORRELATION_SPACECHANNEL - def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): + def __init__(self, domain_geometry, bnd_cond = 'Neumann', **kwargs): """Constructor method """ - super(Gradient, self).__init__() - + backend = kwargs.get('backend',C) correlation = kwargs.get('correlation',CORRELATION_SPACE) - if correlation == CORRELATION_SPACE and gm_domain.channels > 1: + if correlation == CORRELATION_SPACE and domain_geometry.channels > 1: #numpy implementation only for now backend = NUMPY warnings.warn("Warning: correlation='Space' on multi-channel dataset will use `numpy` backend") if backend == NUMPY: - self.operator = Gradient_numpy(gm_domain, bnd_cond=bnd_cond, **kwargs) + self.operator = Gradient_numpy(domain_geometry, bnd_cond=bnd_cond, **kwargs) else: - self.operator = Gradient_C(gm_domain, bnd_cond=bnd_cond, **kwargs) + self.operator = Gradient_C(domain_geometry, bnd_cond=bnd_cond, **kwargs) + + super(Gradient, self).__init__(domain_geometry=domain_geometry, + range_geometry=self.operator.range_geometry()) + + self.gm_range = self.range_geometry() + self.gm_domain = self.domain_geometry() def direct(self, x, out=None): @@ -120,16 +125,16 @@ class Gradient(LinearOperator): """ return self.operator.adjoint(x, out=out) - def domain_geometry(self): - '''Returns domain_geometry of Gradient''' + # def domain_geometry(self): + # '''Returns domain_geometry of Gradient''' - return self.operator.gm_domain + # return self.operator.gm_domain - def range_geometry(self): + # def range_geometry(self): - '''Returns range_geometry of Gradient''' + # '''Returns range_geometry of Gradient''' - return self.operator.gm_range + # return self.operator.gm_range class Gradient_numpy(LinearOperator): @@ -143,7 +148,6 @@ class Gradient_numpy(LinearOperator): :param correlation: optional, :code:`SpaceChannel` or :code:`Space` :type correlation: str, optional, default :code:`Space` ''' - super(Gradient_numpy, self).__init__() self.gm_domain = gm_domain # Domain of Grad Operator @@ -190,6 +194,10 @@ class Gradient_numpy(LinearOperator): # Call FiniteDiff operator self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond) + + super(Gradient_numpy, self).__init__(domain_geometry=self.gm_domain, + range_geometry=self.gm_range) + print("Initialised GradientOperator with numpy backend") def direct(self, x, out=None): @@ -240,14 +248,6 @@ class Gradient_numpy(LinearOperator): return self.gm_range - def __rmul__(self, scalar): - - '''Multiplication of Gradient with a scalar - - Returns: ScaledOperator - ''' - - return ScaledOperator(self, scalar) ########################################################################### ############### For preconditioning ###################################### @@ -348,8 +348,6 @@ class Gradient_C(LinearOperator): def __init__(self, gm_domain, gm_range=None, bnd_cond = NEUMANN, **kwargs): - super(Gradient_C, self).__init__() - self.num_threads = kwargs.get('num_threads',NUM_THREADS) self.gm_domain = gm_domain @@ -375,6 +373,9 @@ class Gradient_C(LinearOperator): else: raise ValueError('Number of dimensions not supported, expected 2, 3 or 4, got {}'.format(len(gm_domain.shape))) #self.num_threads + # super(Gradient_C, self).__init__() + super(Gradient_C, self).__init__(domain_geometry=self.gm_domain, + range_geometry=self.gm_range) print("Initialised GradientOperator with C backend running with ", cilacc.openMPtest(self.num_threads)," threads") @staticmethod @@ -418,17 +419,17 @@ class Gradient_C(LinearOperator): if return_val is True: return out - def domain_geometry(self): + # def domain_geometry(self): - '''Returns domain_geometry of Gradient''' + # '''Returns domain_geometry of Gradient''' - return self.gm_domain + # return self.gm_domain - def range_geometry(self): + # def range_geometry(self): - '''Returns range_geometry of Gradient''' + # '''Returns range_geometry of Gradient''' - return self.gm_range + # return self.gm_range if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py index a4c6ae3..2c14c5d 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py @@ -35,14 +35,14 @@ class Identity(LinearOperator): ''' - def __init__(self, gm_domain, gm_range=None): + def __init__(self, domain_geometry, range_geometry=None): - self.gm_domain = gm_domain - self.gm_range = gm_range - if self.gm_range is None: - self.gm_range = self.gm_domain - super(Identity, self).__init__() + if range_geometry is None: + range_geometry = domain_geometry + + super(Identity, self).__init__(domain_geometry=domain_geometry, + range_geometry=range_geometry) def direct(self,x,out=None): @@ -68,18 +68,6 @@ class Identity(LinearOperator): '''Evaluates operator norm of Identity''' return 1.0 - - 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 ########################################################################### diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py index 6f13532..95bcc92 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py @@ -39,15 +39,16 @@ class LinearOperatorMatrix(LinearOperator): ''' self.A = A M_A, N_A = self.A.shape - self.gm_domain = VectorGeometry(N_A) - self.gm_range = VectorGeometry(M_A) + domain_geometry = VectorGeometry(N_A) + range_geometry = VectorGeometry(M_A) self.s1 = None # Largest singular value, initially unknown - super(LinearOperatorMatrix, self).__init__() + super(LinearOperatorMatrix, self).__init__(domain_geometry=domain_geometry, + range_geometry=range_geometry) def direct(self,x, out=None): if out is None: - tmp = self.gm_range.allocate() + tmp = self.range_geometry().allocate() tmp.fill(numpy.dot(self.A,x.as_array())) return tmp else: @@ -58,7 +59,7 @@ class LinearOperatorMatrix(LinearOperator): def adjoint(self,x, out=None): if out is None: - tmp = self.gm_domain.allocate() + tmp = self.domain_geometry().allocate() tmp.fill(numpy.dot(self.A.transpose(),x.as_array())) return tmp else: @@ -70,8 +71,4 @@ class LinearOperatorMatrix(LinearOperator): def calculate_norm(self, **kwargs): # If unknown, compute and store. If known, simply return it. return svds(self.A,1,return_singular_vectors=False)[0] - - def domain_geometry(self): - return self.gm_domain - def range_geometry(self): - return self.gm_range + diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py index 14d53e8..d49bc1a 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -18,13 +18,23 @@ from __future__ import absolute_import from __future__ import division from __future__ import print_function - -from ccpi.optimisation.operators.ScaledOperator import ScaledOperator +from numbers import Number +import numpy +import functools class Operator(object): '''Operator that maps from a space X -> Y''' - def __init__(self, **kwargs): - self.__norm = None + def __init__(self, domain_geometry, **kwargs): + r''' + Creator + + :param domain_geometry: domain of the operator + :param range_geometry: range of the operator + :type range_geometry: optional, default None + ''' + self._norm = None + self._domain_geometry = domain_geometry + self._range_geometry = kwargs.get('range_geometry', None) def is_linear(self): '''Returns if the operator is linear''' @@ -51,20 +61,425 @@ class Operator(object): :parameter force: forces the recalculation of the norm :type force: boolean, default :code:`False` ''' - if self.__norm is None or kwargs.get('force', False): - self.__norm = self.calculate_norm(**kwargs) - return self.__norm + if self._norm is None or kwargs.get('force', False): + self._norm = self.calculate_norm(**kwargs) + return self._norm def calculate_norm(self, **kwargs): '''Calculates the norm of the Operator''' raise NotImplementedError def range_geometry(self): '''Returns the range of the Operator: Y space''' - raise NotImplementedError + return self._range_geometry def domain_geometry(self): '''Returns the domain of the Operator: X space''' - raise NotImplementedError + return self._domain_geometry def __rmul__(self, scalar): '''Defines the multiplication by a scalar on the left returns a ScaledOperator''' return ScaledOperator(self, scalar) + + def compose(self, *other, **kwargs): + # TODO: check equality of domain and range of operators + #if self.operator2.range_geometry != self.operator1.domain_geometry: + # raise ValueError('Cannot compose operators, check domain geometry of {} and range geometry of {}'.format(self.operato1,self.operator2)) + + return CompositionOperator(self, *other, **kwargs) + + def __add__(self, other): + return SumOperator(self, other) + + def __mul__(self, scalar): + return self.__rmul__(scalar) + + def __neg__(self): + """ Return -self """ + return -1 * self + + def __sub__(self, other): + """ Returns the subtraction of the operators.""" + return self + (-1) * other + + +class LinearOperator(Operator): + '''A Linear Operator that maps from a space X <-> Y''' + def __init__(self, domain_geometry, **kwargs): + super(LinearOperator, self).__init__(domain_geometry, **kwargs) + def is_linear(self): + '''Returns if the operator is linear''' + return True + def adjoint(self,x, out=None): + '''returns the adjoint/inverse operation + + only available to linear operators''' + raise NotImplementedError + + @staticmethod + def PowerMethod(operator, iterations, x_init=None): + '''Power method to calculate iteratively the Lipschitz constant + + :param operator: input operator + :type operator: :code:`LinearOperator` + :param iterations: number of iterations to run + :type iteration: int + :param x_init: starting point for the iteration in the operator domain + :returns: tuple with: L, list of L at each iteration, the data the iteration worked on. + ''' + + # Initialise random + if x_init is None: + x0 = operator.domain_geometry().allocate('random') + else: + x0 = x_init.copy() + + x1 = operator.domain_geometry().allocate() + y_tmp = operator.range_geometry().allocate() + s = numpy.zeros(iterations) + # Loop + for it in numpy.arange(iterations): + operator.direct(x0,out=y_tmp) + operator.adjoint(y_tmp,out=x1) + x1norm = x1.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 + + def calculate_norm(self, **kwargs): + '''Returns the norm of the LinearOperator as calculated by the PowerMethod + + :param iterations: number of iterations to run + :type iterations: int, optional, default = 25 + :param x_init: starting point for the iteration in the operator domain + :type x_init: same type as domain, a subclass of :code:`DataContainer`, optional, None + :parameter force: forces the recalculation of the norm + :type force: boolean, default :code:`False` + ''' + x0 = kwargs.get('x_init', None) + iterations = kwargs.get('iterations', 25) + 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, **kwargs): + r'''Does a dot linearity test on the operator + + Evaluates if the following equivalence holds + + .. math:: + + Ax\times y = y \times A^Tx + + The equivalence is tested within a user specified precision + + .. code:: + + abs(desired-actual) < 1.5 * 10**(-decimal) + + :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. + :param decimal: desired precision + :type decimal: int, optional, default 4 + ''' + 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)) + print ("decimal ", kwargs.get('decimal', 4)) + try: + numpy.testing.assert_almost_equal(abs((a-b)/a), 0., decimal=kwargs.get('decimal',4)) + return True + except AssertionError as ae: + print (ae) + return False + + +class ScaledOperator(Operator): + + + '''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. + + :param operator: a Operator or LinearOperator + :param scalar: a scalar multiplier + + Example: + The scaled operator behaves like the following: + + .. 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, **kwargs): + '''creator + + :param operator: a Operator or LinearOperator + :param scalar: a scalar multiplier + :type scalar: float''' + + super(ScaledOperator, self).__init__(domain_geometry=operator.domain_geometry(), + range_geometry=operator.range_geometry()) + if not isinstance (scalar, Number): + raise TypeError('expected scalar: got {}'.format(type(scalar))) + self.scalar = scalar + self.operator = operator + def direct(self, x, out=None): + '''direct method''' + if out is None: + return self.scalar * self.operator.direct(x, out=out) + else: + self.operator.direct(x, out=out) + out *= self.scalar + def adjoint(self, x, out=None): + '''adjoint method''' + if self.operator.is_linear(): + if out is None: + return self.scalar * self.operator.adjoint(x, out=out) + else: + self.operator.adjoint(x, out=out) + out *= self.scalar + else: + raise TypeError('Operator is not linear') + def norm(self, **kwargs): + '''norm of the operator''' + return numpy.abs(self.scalar) * self.operator.norm(**kwargs) + # def range_geometry(self): + # '''range of the operator''' + # return self.operator.range_geometry() + # def domain_geometry(self): + # '''domain of the operator''' + # return self.operator.domain_geometry() + def is_linear(self): + '''returns whether the operator is linear + + :returns: boolean ''' + return self.operator.is_linear() + + +############################################################################### +################ SumOperator ########################################### +############################################################################### + +class SumOperator(Operator): + + + def __init__(self, operator1, operator2): + + self.operator1 = operator1 + self.operator2 = operator2 + + # if self.operator1.domain_geometry() != self.operator2.domain_geometry(): + # raise ValueError('Domain geometry of {} is not equal with domain geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__)) + + # if self.operator1.range_geometry() != self.operator2.range_geometry(): + # raise ValueError('Range geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__)) + + self.linear_flag = self.operator1.is_linear() and self.operator2.is_linear() + + super(SumOperator, self).__init__(domain_geometry=self.operator1.domain_geometry(), + range_geometry=self.operator1.range_geometry()) + + def direct(self, x, out=None): + + if out is None: + return self.operator1.direct(x) + self.operator2.direct(x) + else: + #TODO check if allcating with tmp is faster + self.operator1.direct(x, out=out) + out.add(self.operator2.direct(x), out=out) + + def adjoint(self, x, out=None): + + if self.linear_flag: + if out is None: + return self.operator1.adjoint(x) + self.operator2.adjoint(x) + else: + #TODO check if allcating with tmp is faster + self.operator1.adjoint(x, out=out) + out.add(self.operator2.adjoint(x), out=out) + else: + raise ValueError('No adjoint operation with non-linear operators') + + def is_linear(self): + return self.linear_flag + + def calculate_norm(self, **kwargs): + if self.is_linear(): + return LinearOperator.calculate_norm(self, **kwargs) + +############################################################################### +################ Composition ########################################### +############################################################################### + +class Composition2Operator(Operator): + + def __init__(self, operator1, operator2): + + self.operator1 = operator1 + self.operator2 = operator2 + + self.linear_flag = self.operator1.is_linear() and self.operator2.is_linear() + + if self.operator2.range_geometry() != self.operator1.domain_geometry(): + raise ValueError('Domain geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__)) + + super(Composition2Operator, self).__init__(self.operator1.domain_geometry(), + self.operator2.range_geometry()) + + def direct(self, x, out = None): + + if out is None: + return self.operator1.direct(self.operator2.direct(x)) + else: + tmp = self.operator2.range_geometry().allocate() + self.operator2.direct(x, out = tmp) + self.operator1.direct(tmp, out = out) + + def adjoint(self, x, out = None): + + if self.linear_flag: + + if out is None: + return self.operator2.adjoint(self.operator1.adjoint(x)) + else: + + tmp = self.operator1.domain_geometry().allocate() + self.operator1.adjoint(x, out=tmp) + self.operator2.adjoint(tmp, out=out) + else: + raise ValueError('No adjoint operation with non-linear operators') + + + def is_linear(self): + return self.linear_flag + + def calculate_norm(self, **kwargs): + if self.is_linear(): + return LinearOperator.calculate_norm(self, **kwargs) + +class CompositionOperator(Operator): + + def __init__(self, *operators, **kwargs): + + # get a reference to the operators + self.operators = operators + + self.linear_flag = functools.reduce(lambda x,y: x and y.is_linear(), + self.operators, True) + self.preallocate = kwargs.get('preallocate', False) + if self.preallocate: + self.tmp_domain = [op.domain_geometry().allocate() for op in self.operators[:-1]] + self.tmp_range = [op.range_geometry().allocate() for op in self.operators[1:]] + # pass + + # TODO address the equality of geometries + # if self.operator2.range_geometry() != self.operator1.domain_geometry(): + # raise ValueError('Domain geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__)) + + super(CompositionOperator, self).__init__( + domain_geometry=self.operators[-1].domain_geometry(), + range_geometry=self.operators[0].range_geometry()) + + def direct(self, x, out = None): + + if out is None: + #return self.operator1.direct(self.operator2.direct(x)) + # return functools.reduce(lambda X,operator: operator.direct(X), + # self.operators[::-1][1:], + # self.operators[-1].direct(x)) + for i,operator in enumerate(self.operators[::-1]): + if i == 0: + step = operator.direct(x) + else: + step = operator.direct(step) + return step + + else: + # tmp = self.operator2.range_geometry().allocate() + # self.operator2.direct(x, out = tmp) + # self.operator1.direct(tmp, out = out) + + # out.fill ( + # functools.reduce(lambda X,operator: operator.direct(X), + # self.operators[::-1][1:], + # self.operators[-1].direct(x)) + # ) + + # TODO this is a bit silly but will handle the pre allocation later + + for i,operator in enumerate(self.operators[::-1]): + if i == 0: + operator.direct(x, out=self.tmp_range[i]) + elif i == len(self.operators) - 1: + operator.direct(self.tmp_range[i-1], out=out) + else: + operator.direct(self.tmp_range[i-1], out=self.tmp_range[i]) + + + def adjoint(self, x, out = None): + + if self.linear_flag: + + if out is not None: + # return self.operator2.adjoint(self.operator1.adjoint(x)) + # return functools.reduce(lambda X,operator: operator.adjoint(X), + # self.operators[1:], + # self.operators[0].adjoint(x)) + + for i,operator in enumerate(self.operators): + if i == 0: + operator.adjoint(x, out=self.tmp_domain[i]) + elif i == len(self.operators) - 1: + step = operator.adjoint(self.tmp_domain[i-1], out=out) + else: + operator.adjoint(self.tmp_domain[i-1], out=self.tmp_domain[i]) + return + + else: + for i,operator in enumerate(self.operators): + if i == 0: + step = operator.adjoint(x) + else: + step = operator.adjoint(step) + + return step + else: + raise ValueError('No adjoint operation with non-linear operators') + + + def is_linear(self): + return self.linear_flag + + def calculate_norm(self, **kwargs): + if self.is_linear(): + return LinearOperator.calculate_norm(self, **kwargs) + + + diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py index 747db65..310d4e8 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py @@ -29,16 +29,16 @@ class SparseFiniteDiff(object): '''Create Sparse Matrices for the Finite Difference Operator''' - def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'): + def __init__(self, domain_geometry, range_geometry=None, + direction=0, bnd_cond = 'Neumann'): - super(SparseFiniteDiff, self).__init__() - self.gm_domain = gm_domain - self.gm_range = gm_range + super(SparseFiniteDiff, self).__init__(domain_geometry=domain_geometry, + range_geometry=range_geometry) self.direction = direction self.bnd_cond = bnd_cond - if self.gm_range is None: - self.gm_range = self.gm_domain + if self.range_geometry is None: + self.range_geometry = self.domain_geometry self.get_dims = [i for i in gm_domain.shape] diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index 6fd2e86..9c2300f 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -48,33 +48,36 @@ class SymmetrizedGradient(LinearOperator): CORRELATION_SPACE = "Space" CORRELATION_SPACECHANNEL = "SpaceChannels" - def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): + def __init__(self, domain_geometry, bnd_cond = 'Neumann', **kwargs): '''creator - :param gm_domain: domain of the operator + :param domain_geometry: domain of the operator :param bnd_cond: boundary condition, either :code:`Neumann` or :code:`Periodic`. :type bnd_cond: str, optional, default :code:`Neumann` :param correlation: :code:`SpaceChannel` or :code:`Channel` :type correlation: str, optional, default :code:`Channel` ''' - super(SymmetrizedGradient, self).__init__() + # super(SymmetrizedGradient, self).__init__(domain_geometry=domain_geometry) - self.gm_domain = gm_domain self.bnd_cond = bnd_cond self.correlation = kwargs.get('correlation',SymmetrizedGradient.CORRELATION_SPACE) - tmp_gm = len(self.gm_domain.geometries)*self.gm_domain.geometries + tmp_gm = len(domain_geometry.geometries)*domain_geometry.geometries - self.gm_range = BlockGeometry(*tmp_gm) # Define FD operator. We need one geometry from the BlockGeometry of the domain - self.FD = FiniteDiff(self.gm_domain.get_item(0), direction = 0, bnd_cond = self.bnd_cond) + self.FD = FiniteDiff(domain_geometry.get_item(0), direction = 0, + bnd_cond = self.bnd_cond) - if self.gm_domain.shape[0]==2: + if domain_geometry.shape[0]==2: self.order_ind = [0,2,1,3] else: self.order_ind = [0,3,6,1,4,7,2,5,8] - + + super(SymmetrizedGradient, self).__init__( + domain_geometry=domain_geometry, + range_geometry=BlockGeometry(*tmp_gm)) + def direct(self, x, out=None): @@ -83,7 +86,7 @@ class SymmetrizedGradient(LinearOperator): if out is None: tmp = [] - for i in range(self.gm_domain.shape[0]): + for i in range(self.domain_geometry().shape[0]): for j in range(x.shape[0]): self.FD.direction = i tmp.append(self.FD.adjoint(x.get_item(j))) @@ -97,7 +100,7 @@ class SymmetrizedGradient(LinearOperator): else: ind = 0 - for i in range(self.gm_domain.shape[0]): + for i in range(self.domain_geometry().shape[0]): for j in range(x.shape[0]): self.FD.direction = i self.FD.adjoint(x.get_item(j), out=out[ind]) @@ -110,12 +113,12 @@ class SymmetrizedGradient(LinearOperator): if out is None: - tmp = [None]*self.gm_domain.shape[0] + tmp = [None]*self.domain_geometry().shape[0] i = 0 - for k in range(self.gm_domain.shape[0]): + for k in range(self.domain_geometry().shape[0]): tmp1 = 0 - for j in range(self.gm_domain.shape[0]): + for j in range(self.domain_geometry().shape[0]): self.FD.direction = j tmp1 += self.FD.direct(x[i]) i+=1 @@ -125,23 +128,18 @@ class SymmetrizedGradient(LinearOperator): else: - tmp = self.gm_domain.allocate() + tmp = self.domain_geometry().allocate() i = 0 - for k in range(self.gm_domain.shape[0]): + for k in range(self.domain_geometry().shape[0]): tmp1 = 0 - for j in range(self.gm_domain.shape[0]): + for j in range(self.domain_geometry().shape[0]): self.FD.direction = j self.FD.direct(x[i], out=tmp[j]) i+=1 tmp1+=tmp[j] out[k].fill(tmp1) - - def domain_geometry(self): - return self.gm_domain - def range_geometry(self): - return self.gm_range if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py index 0feafd8..bc4f08e 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -41,14 +41,12 @@ class ZeroOperator(LinearOperator): ''' - def __init__(self, gm_domain, gm_range=None): - - super(ZeroOperator, self).__init__() + def __init__(self, domain_geometry, range_geometry=None): + if range_geometry is None: + range_geometry = domain_geometry.clone() + super(ZeroOperator, self).__init__(domain_geometry=domain_geometry, + range_geometry=range_geometry) - self.gm_domain = gm_domain - self.gm_range = gm_range - if self.gm_range is None: - self.gm_range = self.gm_domain def direct(self,x,out=None): @@ -57,18 +55,18 @@ class ZeroOperator(LinearOperator): if out is None: - return self.gm_range.allocate() + return self.range_geometry().allocate() else: - out.fill(self.gm_range.allocate()) + out.fill(self.range_geometry.allocate()) def adjoint(self,x, out=None): '''Returns O^{*}(y)''' if out is None: - return self.gm_domain.allocate() + return self.domain_geometry().allocate() else: - out.fill(self.gm_domain.allocate()) + out.fill(self.domain_geometry().allocate()) def calculate_norm(self, **kwargs): @@ -76,15 +74,4 @@ class ZeroOperator(LinearOperator): return 0 - 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/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py index 15c8932..cc25c92 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py @@ -16,11 +16,12 @@ # See the License for the specific language governing permissions and # limitations under the License. -from .Operator import Operator -from .LinearOperator import LinearOperator -from .ScaledOperator import ScaledOperator +from .Operator import Operator, LinearOperator, ScaledOperator, SumOperator,\ + CompositionOperator, Composition2Operator +#from .LinearOperator import LinearOperator +#from .ScaledOperator import ScaledOperator from .BlockOperator import BlockOperator -from .BlockScaledOperator import BlockScaledOperator +# from .BlockScaledOperator import BlockScaledOperator from .SparseFiniteDiff import SparseFiniteDiff from .ShrinkageOperator import ShrinkageOperator diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py index d0f568b..2f33f08 100644 --- a/Wrappers/Python/test/test_BlockOperator.py +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -81,8 +81,8 @@ class TestBlockOperator(unittest.TestCase): ImageGeometry(10,22,31) , \ ImageGeometry(10,20,31) ] x = [ g.allocate() for g in ig ] - ops = [ Identity(g, gm_range=r) for g,r in zip(ig, rg0) ] - ops += [ Identity(g, gm_range=r) for g,r in zip(ig, rg1) ] + ops = [ Identity(g, range_geometry=r) for g,r in zip(ig, rg0) ] + ops += [ Identity(g, range_geometry=r) for g,r in zip(ig, rg1) ] K = BlockOperator(*ops, shape=(2,3)) print ("K col comp? " , K.column_wise_compatible()) diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 57dd41f..8312c6d 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -23,7 +23,12 @@ import numpy from timeit import default_timer as timer from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff from ccpi.optimisation.operators import LinearOperator, LinearOperatorMatrix +import numpy +from ccpi.optimisation.operators import SumOperator, Gradient,\ + ZeroOperator, SymmetrizedGradient, CompositionOperator +from ccpi.framework import TestData +import os def dt(steps): return steps[-1] - steps[-2] @@ -221,6 +226,10 @@ class TestOperator(CCPiTestClass): for n in [norm, norm2, norm3, norm4, norm5]: print ("norm {}", format(n)) + + + + class TestGradients(CCPiTestClass): def setUp(self): @@ -339,6 +348,7 @@ class TestGradients(CCPiTestClass): + class TestBlockOperator(unittest.TestCase): def assertBlockDataContainerEqual(self, container1, container2): print ("assert Block Data Container Equal") @@ -645,6 +655,260 @@ class TestBlockOperator(unittest.TestCase): u1 = B.adjoint(w) self.assertEqual((w * w1).sum() , (u1*u).sum()) +class TestOperatorCompositionSum(unittest.TestCase): + def setUp(self): + + self.data = TestData().load(TestData.BOAT, size=(128,128)) + self.ig = self.data.geometry + + def test_SumOperator(self): + + # numpy.random.seed(1) + ig = self.ig + data = self.data + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + # z = ZeroOperator(domain_geometry=ig) + # sym = SymmetrizedGradient(domain_geometry=ig) + + c = SumOperator(Id1,Id2) + out = c.direct(data) + + numpy.testing.assert_array_almost_equal(out.as_array(),3 * data.as_array()) + + + def test_CompositionOperator_direct1(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = CompositionOperator(G, Id2) + + out1 = G.direct(data) + out2 = d.direct(data) + + + numpy.testing.assert_array_almost_equal(out2.get_item(0).as_array(), out1.get_item(0).as_array()) + numpy.testing.assert_array_almost_equal(out2.get_item(1).as_array(), out1.get_item(1).as_array()) + + def test_CompositionOperator_direct2(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = CompositionOperator(G, Id2) + + out1 = G.direct(data) + + d_out = d.direct(data) + + d1 = Id2.direct(data) + d2 = G.direct(d1) + + numpy.testing.assert_array_almost_equal(d2.get_item(0).as_array(), + d_out.get_item(0).as_array()) + numpy.testing.assert_array_almost_equal(d2.get_item(1).as_array(), + d_out.get_item(1).as_array()) + + def test_CompositionOperator_direct3(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = CompositionOperator(G, Id2) + + out1 = G.direct(data) + + d_out = d.direct(data) + + d1 = Id2.direct(data) + d2 = G.direct(d1) + + numpy.testing.assert_array_almost_equal(d2.get_item(0).as_array(), + d_out.get_item(0).as_array()) + numpy.testing.assert_array_almost_equal(d2.get_item(1).as_array(), + d_out.get_item(1).as_array()) + + G2Id = G.compose(2*Id2) + d2g = G2Id.direct(data) + + numpy.testing.assert_array_almost_equal(d2g.get_item(0).as_array(), + 2 * d_out.get_item(0).as_array()) + numpy.testing.assert_array_almost_equal(d2g.get_item(1).as_array(), + 2 * d_out.get_item(1).as_array()) + + def test_CompositionOperator_direct4(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = CompositionOperator(G, Id1, Id2) + + out1 = G.direct(data) + + d_out = d.direct(data) + + numpy.testing.assert_array_almost_equal(d_out.get_item(0).as_array(), + 2 * out1.get_item(0).as_array()) + numpy.testing.assert_array_almost_equal(d_out.get_item(1).as_array(), + 2 * out1.get_item(1).as_array()) + + def test_CompositionOperator_adjoint1(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = CompositionOperator(G, Id2) + da = d.direct(data) + + out1 = G.adjoint(da) + out2 = d.adjoint(da) + + numpy.testing.assert_array_almost_equal(out2.as_array(), out1.as_array()) + + def test_CompositionOperator_adjoint2(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = CompositionOperator(G, Id1) + da = d.direct(data) + + out1 = G.adjoint(da) + out2 = d.adjoint(da) + + numpy.testing.assert_array_almost_equal(out2.as_array(), 2 * out1.as_array()) + def test_CompositionOperator_adjoint3(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = G.compose(Id1) + da = d.direct(data) + + out1 = G.adjoint(da) + out2 = d.adjoint(da) + + numpy.testing.assert_array_almost_equal(out2.as_array(), 2 * out1.as_array()) + + + def test_CompositionOperator_adjoint4(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + + d = G.compose(-Id1) + da = d.direct(data) + + out1 = G.adjoint(da) + out2 = d.adjoint(da) + + numpy.testing.assert_array_almost_equal(out2.as_array(), - 2 * out1.as_array()) + def test_CompositionOperator_adjoint5(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 3 * Identity(ig) + Id = Id1 - Identity(ig) + d = G.compose(Id) + da = d.direct(data) + + out1 = G.adjoint(da) + out2 = d.adjoint(da) + + numpy.testing.assert_array_almost_equal(out2.as_array(), 2 * out1.as_array()) + + def test_CompositionOperator_adjoint6(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 3 * Identity(ig) + Id = ZeroOperator(ig) + d = G.compose(Id) + da = d.direct(data) + + out1 = G.adjoint(da) + out2 = d.adjoint(da) + + numpy.testing.assert_array_almost_equal(out2.as_array(), 0 * out1.as_array()) + + def stest_CompositionOperator_direct4(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + sym = SymmetrizedGradient(domain_geometry=ig) + Id2 = Identity(ig) + + d = CompositionOperator(sym, Id2) + + out1 = G.direct(data) + out2 = d.direct(data) + + + numpy.testing.assert_array_almost_equal(out2.get_item(0).as_array(), out1.get_item(0).as_array()) + numpy.testing.assert_array_almost_equal(out2.get_item(1).as_array(), out1.get_item(1).as_array()) + + def test_CompositionOperator_adjoint7(self): + ig = self.ig + data = self.data + G = Gradient(domain_geometry=ig) + + + Id1 = 2 * Identity(ig) + Id2 = Identity(ig) + + d = CompositionOperator(G, Id1, Id2) + + out1 = G.direct(data) + out2 = G.adjoint(out1) + + d_out = d.adjoint(out1) + + numpy.testing.assert_array_almost_equal(d_out.as_array(), + 2 * out2.as_array()) + numpy.testing.assert_array_almost_equal(d_out.as_array(), + 2 * out2.as_array()) + + + -- cgit v1.2.3