diff options
14 files changed, 932 insertions, 186 deletions
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()) + + + |