From 13054c64cff2b506a7affa9d21445f368abbb15b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 28 Feb 2019 18:58:33 +0000 Subject: Proposal of Algorithm class (#179) * initial revision * Removed class members of Algorithm class added update_objective * initial version. Fix inline __idiv__ * First implementation of CompositeOperator/DataContainer * removed __getitem__ added get_item added shape * added CGLS * working unit test, initial tomography test * added reverse multiplication of operator with number * added operators directory * fixed typo * added unittest for CompositeDataContainer * fix TomoIdentity with scalar * check numerical types from numpy * add default stop criterion and run method * add run method * first working implementation of CGLS with CompositeOperator/DataContainer notice problem with _rmul_ and _mul_ methods precedence with numpy. * new Algorithm class and algorithms in separate files Added new Algorithm class and derivatives in different files for GradientDescent, CGLS, FBPD, FISTA * added algorithms and restored CIL_VERSION env variable * removed Algorithms.py * modified run and renamed a few members/methods * uses squared_norm * renamed get_current_objective to get_last_objective update_objective can be issued every N iteration, default 1. fixed run method to run N iterations within the stop criterion. * load class as module files * force py line endings to LF * updates * call super __init__ as first thing * unit tests are now to be found in test directory unit tests are now split in several files in the directory test * install algorithms module * Implementation with Algorithm * skip Reader tests * unittest for linux * commented not needed import Iterable * removed explicit return from __init__ * remove composite operator file --- Wrappers/Python/ccpi/framework.py | 12 +- .../ccpi/optimisation/algorithms/Algorithm.py | 157 ++++ .../Python/ccpi/optimisation/algorithms/CGLS.py | 87 ++ .../Python/ccpi/optimisation/algorithms/FBPD.py | 86 ++ .../Python/ccpi/optimisation/algorithms/FISTA.py | 121 +++ .../optimisation/algorithms/GradientDescent.py | 72 ++ .../ccpi/optimisation/algorithms/__init__.py | 29 + Wrappers/Python/ccpi/optimisation/ops.py | 48 +- Wrappers/Python/conda-recipe/meta.yaml | 9 + Wrappers/Python/conda-recipe/run_test.py | 985 --------------------- Wrappers/Python/setup.py | 8 +- Wrappers/Python/test/TestReader.py | 134 --- Wrappers/Python/test/__init__.py | 0 Wrappers/Python/test/skip_TestReader.py | 134 +++ Wrappers/Python/test/test_DataContainer.py | 499 +++++++++++ Wrappers/Python/test/test_NexusReader.py | 96 ++ Wrappers/Python/test/test_algorithms.py | 123 +++ Wrappers/Python/test/test_run_test.py | 435 +++++++++ 18 files changed, 1903 insertions(+), 1132 deletions(-) create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/__init__.py delete mode 100755 Wrappers/Python/conda-recipe/run_test.py delete mode 100644 Wrappers/Python/test/TestReader.py create mode 100644 Wrappers/Python/test/__init__.py create mode 100644 Wrappers/Python/test/skip_TestReader.py create mode 100755 Wrappers/Python/test/test_DataContainer.py create mode 100755 Wrappers/Python/test/test_NexusReader.py create mode 100755 Wrappers/Python/test/test_algorithms.py create mode 100755 Wrappers/Python/test/test_run_test.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index e1a2dff..0c23628 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -473,7 +473,10 @@ class DataContainer(object): else: raise ValueError('*:Wrong shape: {0} and {1}'.format(self.shape, other.shape)) - elif isinstance(other, (int, float, complex)): + elif isinstance(other, (int, float, complex,\ + numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\ + numpy.float, numpy.float16, numpy.float32, numpy.float64, \ + numpy.complex)): return type(self)(self.as_array() * other, deep_copy=True, dimension_labels=self.dimension_labels, @@ -559,6 +562,8 @@ class DataContainer(object): # __isub__ def __idiv__(self, other): + return self.__itruediv__(other) + def __itruediv__(self, other): if isinstance(other, (int, float)) : numpy.divide(self.array, other, out=self.array) elif issubclass(type(other), DataContainer): @@ -632,6 +637,10 @@ class DataContainer(object): if out is None: if isinstance(x2, (int, float, complex)): out = pwop(self.as_array() , x2 , *args, **kwargs ) + elif isinstance(x2, (numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\ + numpy.float, numpy.float16, numpy.float32, numpy.float64, \ + numpy.complex)): + out = pwop(self.as_array() , x2 , *args, **kwargs ) elif issubclass(type(x2) , DataContainer): out = pwop(self.as_array() , x2.as_array() , *args, **kwargs ) return type(self)(out, @@ -731,6 +740,7 @@ class DataContainer(object): '''return the euclidean norm of the DataContainer viewed as a vector''' return numpy.sqrt(self.squared_norm()) + class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' def __init__(self, diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py new file mode 100755 index 0000000..680b268 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -0,0 +1,157 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import time +from numbers import Integral + +class Algorithm(object): + '''Base class for iterative algorithms + + provides the minimal infrastructure. + Algorithms are iterables so can be easily run in a for loop. They will + stop as soon as the stop cryterion is met. + The user is required to implement the set_up, __init__, update and + and update_objective methods + + A courtesy method run is available to run n iterations. The method accepts + a callback function that receives the current iteration number and the actual objective + value and can be used to trigger print to screens and other user interactions. The run + method will stop when the stopping cryterion is met. + ''' + + def __init__(self): + '''Constructor + + Set the minimal number of parameters: + iteration: current iteration number + max_iteration: maximum number of iterations + memopt: whether to use memory optimisation () + timing: list to hold the times it took to run each iteration + update_objectice_interval: the interval every which we would save the current + objective. 1 means every iteration, 2 every 2 iteration + and so forth. This is by default 1 and should be increased + when evaluating the objective is computationally expensive. + ''' + self.iteration = 0 + self.__max_iteration = 0 + self.__loss = [] + self.memopt = False + self.timing = [] + self.update_objective_interval = 1 + def set_up(self, *args, **kwargs): + '''Set up the algorithm''' + raise NotImplementedError() + def update(self): + '''A single iteration of the algorithm''' + raise NotImplementedError() + + def should_stop(self): + '''default stopping cryterion: number of iterations + + The user can change this in concrete implementatition of iterative algorithms.''' + return self.max_iteration_stop_cryterion() + + def max_iteration_stop_cryterion(self): + '''default stop cryterion for iterative algorithm: max_iteration reached''' + return self.iteration >= self.max_iteration + def __iter__(self): + '''Algorithm is an iterable''' + return self + def next(self): + '''Algorithm is an iterable + + python2 backwards compatibility''' + return self.__next__() + def __next__(self): + '''Algorithm is an iterable + + calling this method triggers update and update_objective + ''' + if self.should_stop(): + raise StopIteration() + else: + time0 = time.time() + self.update() + self.timing.append( time.time() - time0 ) + if self.iteration % self.update_objective_interval == 0: + self.update_objective() + self.iteration += 1 + def get_output(self): + '''Returns the solution found''' + return self.x + def get_last_loss(self): + '''Returns the last stored value of the loss function + + if update_objective_interval is 1 it is the value of the objective at the current + iteration. If update_objective_interval > 1 it is the last stored value. + ''' + return self.__loss[-1] + def get_last_objective(self): + '''alias to get_last_loss''' + return self.get_last_loss() + def update_objective(self): + '''calculates the objective with the current solution''' + raise NotImplementedError() + @property + def loss(self): + '''returns the list of the values of the objective during the iteration + + The length of this list may be shorter than the number of iterations run when + the update_objective_interval > 1 + ''' + return self.__loss + @property + def objective(self): + '''alias of loss''' + return self.loss + @property + def max_iteration(self): + '''gets the maximum number of iterations''' + return self.__max_iteration + @max_iteration.setter + def max_iteration(self, value): + '''sets the maximum number of iterations''' + assert isinstance(value, int) + self.__max_iteration = value + @property + def update_objective_interval(self): + return self.__update_objective_interval + @update_objective_interval.setter + def update_objective_interval(self, value): + if isinstance(value, Integral): + if value >= 1: + self.__update_objective_interval = value + else: + raise ValueError('Update objective interval must be an integer >= 1') + else: + raise ValueError('Update objective interval must be an integer >= 1') + def run(self, iterations, verbose=False, callback=None): + '''run n iterations and update the user with the callback if specified''' + if self.should_stop(): + print ("Stop cryterion has been reached.") + i = 0 + for _ in self: + if verbose: + print ("Iteration {}/{}, objective {}".format(self.iteration, + self.max_iteration, self.get_last_objective()) ) + if callback is not None: + callback(self.iteration, self.get_last_objective()) + i += 1 + if i == iterations: + break + diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py new file mode 100755 index 0000000..7194eb8 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -0,0 +1,87 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:11:23 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +#from collections.abc import Iterable +class CGLS(Algorithm): + + '''Conjugate Gradient Least Squares algorithm + + Parameters: + x_init: initial guess + operator: operator for forward/backward projections + data: data to operate on + ''' + def __init__(self, **kwargs): + super(CGLS, self).__init__() + self.x = kwargs.get('x_init', None) + self.operator = kwargs.get('operator', None) + self.data = kwargs.get('data', None) + if self.x is not None and self.operator is not None and \ + self.data is not None: + print ("Calling from creator") + self.set_up(x_init =kwargs['x_init'], + operator=kwargs['operator'], + data =kwargs['data']) + + def set_up(self, x_init, operator , data ): + + self.r = data.copy() + self.x = x_init.copy() + + self.operator = operator + self.d = operator.adjoint(self.r) + + + self.normr2 = self.d.squared_norm() + #if isinstance(self.normr2, Iterable): + # self.normr2 = sum(self.normr2) + #self.normr2 = numpy.sqrt(self.normr2) + #print ("set_up" , self.normr2) + + def update(self): + + Ad = self.operator.direct(self.d) + #norm = (Ad*Ad).sum() + #if isinstance(norm, Iterable): + # norm = sum(norm) + norm = Ad.squared_norm() + + alpha = self.normr2/norm + self.x += (self.d * alpha) + self.r -= (Ad * alpha) + s = self.operator.adjoint(self.r) + + normr2_new = s.squared_norm() + #if isinstance(normr2_new, Iterable): + # normr2_new = sum(normr2_new) + #normr2_new = numpy.sqrt(normr2_new) + #print (normr2_new) + + beta = normr2_new/self.normr2 + self.normr2 = normr2_new + self.d = s + beta*self.d + + def update_objective(self): + self.loss.append(self.r.squared_norm()) \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py new file mode 100755 index 0000000..322e9eb --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py @@ -0,0 +1,86 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:09:03 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.funcs import ZeroFun + +class FBPD(Algorithm): + '''FBPD Algorithm + + Parameters: + x_init: initial guess + f: constraint + g: data fidelity + h: regularizer + opt: additional algorithm + ''' + constraint = None + data_fidelity = None + regulariser = None + def __init__(self, **kwargs): + pass + def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\ + regulariser=None, opt=None): + + # default inputs + if constraint is None: + self.constraint = ZeroFun() + else: + self.constraint = constraint + if data_fidelity is None: + data_fidelity = ZeroFun() + else: + self.data_fidelity = data_fidelity + if regulariser is None: + self.regulariser = ZeroFun() + else: + self.regulariser = regulariser + + # algorithmic parameters + + + # step-sizes + self.tau = 2 / (self.data_fidelity.L + 2) + self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L + + self.inv_sigma = 1/self.sigma + + # initialization + self.x = x_init + self.y = operator.direct(self.x) + + + def update(self): + + # primal forward-backward step + x_old = self.x + self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) ) + self.x = self.constraint.prox(self.x, self.tau); + + # dual forward-backward step + self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old); + self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma); + + # time and criterion + self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x)) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py new file mode 100755 index 0000000..bc4489e --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 21 11:07:30 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.funcs import ZeroFun +import numpy + +class FISTA(Algorithm): + '''Fast Iterative Shrinkage-Thresholding Algorithm + + Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding + algorithm for linear inverse problems. + SIAM journal on imaging sciences,2(1), pp.183-202. + + Parameters: + x_init: initial guess + f: data fidelity + g: regularizer + h: + opt: additional algorithm + ''' + + def __init__(self, **kwargs): + '''initialisation can be done at creation time if all + proper variables are passed or later with set_up''' + super(FISTA, self).__init__() + self.f = None + self.g = None + self.invL = None + self.t_old = 1 + args = ['x_init', 'f', 'g', 'opt'] + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(kwargs['x_init'], + f=kwargs['f'], + g=kwargs['g'], + opt=kwargs['opt']) + + def set_up(self, x_init, f=None, g=None, opt=None): + + # default inputs + if f is None: + self.f = ZeroFun() + else: + self.f = f + if g is None: + g = ZeroFun() + self.g = g + else: + self.g = g + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-4, 'memopt':False} + + self.tol = opt['tol'] if 'tol' in opt.keys() else 1e-4 + memopt = opt['memopt'] if 'memopt' in opt.keys() else False + self.memopt = memopt + + # initialization + if memopt: + self.y = x_init.clone() + self.x_old = x_init.clone() + self.x = x_init.clone() + self.u = x_init.clone() + else: + self.x_old = x_init.copy() + self.y = x_init.copy() + + #timing = numpy.zeros(max_iter) + #criter = numpy.zeros(max_iter) + + + self.invL = 1/f.L + + self.t_old = 1 + + def update(self): + # algorithm loop + #for it in range(0, max_iter): + + if self.memopt: + # u = y - invL*f.grad(y) + # store the result in x_old + self.f.gradient(self.y, out=self.u) + self.u.__imul__( -self.invL ) + self.u.__iadd__( self.y ) + # x = g.prox(u,invL) + self.g.proximal(self.u, self.invL, out=self.x) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + # y = x + (t_old-1)/t*(x-x_old) + self.x.subtract(self.x_old, out=self.y) + self.y.__imul__ ((self.t_old-1)/self.t) + self.y.__iadd__( self.x ) + + self.x_old.fill(self.x) + self.t_old = self.t + + + else: + u = self.y - self.invL*self.f.grad(self.y) + + self.x = self.g.prox(u,self.invL) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + self.y = self.x + (self.t_old-1)/self.t*(self.x-self.x_old) + + self.x_old = self.x.copy() + self.t_old = self.t + + def update_objective(self): + self.loss.append( self.f(self.x) + self.g(self.x) ) \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py new file mode 100755 index 0000000..7794b4d --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py @@ -0,0 +1,72 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:05:09 2019 + +@author: ofn77899 +""" +from ccpi.optimisation.algorithms import Algorithm + +class GradientDescent(Algorithm): + '''Implementation of Gradient Descent algorithm + ''' + + def __init__(self, **kwargs): + '''initialisation can be done at creation time if all + proper variables are passed or later with set_up''' + super(GradientDescent, self).__init__() + self.x = None + self.rate = 0 + self.objective_function = None + self.regulariser = None + args = ['x_init', 'objective_function', 'rate'] + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(x_init=kwargs['x_init'], + objective_function=kwargs['objective_function'], + rate=kwargs['rate']) + + def should_stop(self): + '''stopping cryterion, currently only based on number of iterations''' + return self.iteration >= self.max_iteration + + def set_up(self, x_init, objective_function, rate): + '''initialisation of the algorithm''' + self.x = x_init.copy() + if self.memopt: + self.x_update = x_init.copy() + self.objective_function = objective_function + self.rate = rate + self.loss.append(objective_function(x_init)) + self.iteration = 0 + + def update(self): + '''Single iteration''' + if self.memopt: + self.objective_function.gradient(self.x, out=self.x_update) + self.x_update *= -self.rate + self.x += self.x_update + else: + self.x += -self.rate * self.objective_function.grad(self.x) + + def update_objective(self): + self.loss.append(self.objective_function(self.x)) + \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py new file mode 100755 index 0000000..52fe6d7 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -0,0 +1,29 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:03:13 2019 + +@author: ofn77899 +""" + +from .Algorithm import Algorithm +from .CGLS import CGLS +from .GradientDescent import GradientDescent +from .FISTA import FISTA +from .FBPD import FBPD \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index b2f996d..e9e7f44 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -24,26 +24,49 @@ from ccpi.framework import AcquisitionData from ccpi.framework import ImageData from ccpi.framework import ImageGeometry from ccpi.framework import AcquisitionGeometry - +from numbers import Number # Maybe operators need to know what types they take as inputs/outputs # to not just use generic DataContainer class Operator(object): + '''Operator that maps from a space X -> Y''' + def __init__(self, **kwargs): + self.scalar = 1 + def is_linear(self): + '''Returns if the operator is linear''' + return False def direct(self,x, out=None): - return x - def adjoint(self,x, out=None): - return x + raise NotImplementedError def size(self): # To be defined for specific class raise NotImplementedError - def get_max_sing_val(self): + def norm(self): raise NotImplementedError def allocate_direct(self): + '''Allocates memory on the Y space''' raise NotImplementedError def allocate_adjoint(self): + '''Allocates memory on the X space''' + raise NotImplementedError + def range_dim(self): raise NotImplementedError + def domain_dim(self): + raise NotImplementedError + def __rmul__(self, other): + '''reverse multiplication of Operator with number sets the variable scalar in the Operator''' + assert isinstance(other, Number) + self.scalar = other + return self +class LinearOperator(Operator): + '''Operator that maps from a space X -> Y''' + def is_linear(self): + '''Returns if the operator is linear''' + return True + def adjoint(self,x, out=None): + raise NotImplementedError + class Identity(Operator): def __init__(self): self.s1 = 1.0 @@ -70,21 +93,26 @@ class Identity(Operator): class TomoIdentity(Operator): def __init__(self, geometry, **kwargs): + super(TomoIdentity, self).__init__() self.s1 = 1.0 self.geometry = geometry - super(TomoIdentity, self).__init__() + def direct(self,x,out=None): + if out is None: + if self.scalar != 1: + return x * self.scalar return x.copy() else: + if self.scalar != 1: + out.fill(x * self.scalar) + return out.fill(x) + return def adjoint(self,x, out=None): - if out is None: - return x.copy() - else: - out.fill(x) + return self.direct(x, out) def size(self): return NotImplemented diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index 1b7cae6..8ded429 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -12,6 +12,15 @@ test: requires: - python-wget - cvxpy # [not win] + + source_files: + - ./test # [win] + - ./ccpi/Wrappers/Python/test # [not win] + + commands: + - python -c "import os; print (os.getcwd())" + - python -m unittest discover # [win] + - python -m unittest discover -s ccpi/Wrappers/Python/test # [not win] requirements: build: diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py deleted file mode 100755 index b52af55..0000000 --- a/Wrappers/Python/conda-recipe/run_test.py +++ /dev/null @@ -1,985 +0,0 @@ -import unittest -import numpy -import numpy as np -from ccpi.framework import DataContainer -from ccpi.framework import ImageData -from ccpi.framework import AcquisitionData -from ccpi.framework import ImageGeometry -from ccpi.framework import AcquisitionGeometry -import sys -from timeit import default_timer as timer -from ccpi.optimisation.algs import FISTA -from ccpi.optimisation.algs import FBPD -from ccpi.optimisation.funcs import Norm2sq -from ccpi.optimisation.funcs import ZeroFun -from ccpi.optimisation.funcs import Norm1 -from ccpi.optimisation.funcs import TV2D -from ccpi.optimisation.funcs import Norm2 - -from ccpi.optimisation.ops import LinearOperatorMatrix -from ccpi.optimisation.ops import TomoIdentity -from ccpi.optimisation.ops import Identity -from ccpi.optimisation.ops import PowerMethodNonsquare - - -import numpy.testing -import wget -import os -from ccpi.io.reader import NexusReader - - -try: - from cvxpy import * - cvx_not_installable = False -except ImportError: - cvx_not_installable = True - - -def aid(x): - # This function returns the memory - # block address of an array. - return x.__array_interface__['data'][0] - - -def dt(steps): - return steps[-1] - steps[-2] - - -class TestDataContainer(unittest.TestCase): - def create_simple_ImageData(self): - N = 64 - ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) - Phantom = ImageData(geometry=ig) - - x = Phantom.as_array() - - x[int(round(N/4)):int(round(3*N/4)), - int(round(N/4)):int(round(3*N/4))] = 0.5 - x[int(round(N/8)):int(round(7*N/8)), - int(round(3*N/8)):int(round(5*N/8))] = 1 - - return (ig, Phantom) - - def create_DataContainer(self, X,Y,Z, value=1): - steps = [timer()] - a = value * numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - return ds - - def test_creation_nocopy(self): - shape = (2, 3, 4, 5) - size = shape[0] - for i in range(1, len(shape)): - size = size * shape[i] - #print("a refcount " , sys.getrefcount(a)) - a = numpy.asarray([i for i in range(size)]) - #print("a refcount " , sys.getrefcount(a)) - a = numpy.reshape(a, shape) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W']) - #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a), 3) - self.assertEqual(ds.dimension_labels, {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'}) - - def testGb_creation_nocopy(self): - X, Y, Z = 512, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - print("test clone") - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a), 3) - ds1 = ds.copy() - self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) - ds1 = ds.clone() - self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) - - def testInlineAlgebra(self): - print("Test Inline Algebra") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - print(t0) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - #ds.__iadd__( 2 ) - ds += 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 3.) - #ds.__isub__( 2 ) - ds -= 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - #ds.__imul__( 2 ) - ds *= 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 2.) - #ds.__idiv__( 2 ) - ds /= 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - ds1 = ds.copy() - #ds1.__iadd__( 1 ) - ds1 += 1 - #ds.__iadd__( ds1 ) - ds += ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 3.) - #ds.__isub__( ds1 ) - ds -= ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - #ds.__imul__( ds1 ) - ds *= ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 2.) - #ds.__idiv__( ds1 ) - ds /= ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - def test_unary_operations(self): - print("Test unary operations") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = -numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - print(t0) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - - ds.sign(out=ds) - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], -1.) - - ds.abs(out=ds) - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - ds.__imul__(2) - ds.sqrt(out=ds) - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], - numpy.sqrt(2., dtype='float32')) - - def test_binary_operations(self): - self.binary_add() - self.binary_subtract() - self.binary_multiply() - self.binary_divide() - - def binary_add(self): - print("Test binary add") - X, Y, Z = 512, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.add(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.add(ds1, out=ds)", dt(steps)) - steps.append(timer()) - ds2 = ds.add(ds1) - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.add(ds1)", dt(steps)) - - self.assertLess(t1, t2) - self.assertEqual(ds.as_array()[0][0][0], 2.) - - ds0 = ds - ds0.add(2, out=ds0) - steps.append(timer()) - print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0]) - self.assertEqual(4., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.add(2) - steps.append(timer()) - print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - - def binary_subtract(self): - print("Test binary subtract") - X, Y, Z = 512, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.subtract(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.subtract(ds1, out=ds)", dt(steps)) - self.assertEqual(0., ds.as_array()[0][0][0]) - - steps.append(timer()) - ds2 = ds.subtract(ds1) - self.assertEqual(-1., ds2.as_array()[0][0][0]) - - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.subtract(ds1)", dt(steps)) - - self.assertLess(t1, t2) - - del ds1 - ds0 = ds.copy() - steps.append(timer()) - ds0.subtract(2, out=ds0) - #ds0.__isub__( 2 ) - steps.append(timer()) - print("ds0.subtract(2,out=ds0)", dt( - steps), -2., ds0.as_array()[0][0][0]) - self.assertEqual(-2., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.subtract(2) - steps.append(timer()) - print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - self.assertEqual(-2., ds0.as_array()[0][0][0]) - self.assertEqual(-4., ds3.as_array()[0][0][0]) - - def binary_multiply(self): - print("Test binary multiply") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.multiply(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.multiply(ds1, out=ds)", dt(steps)) - steps.append(timer()) - ds2 = ds.multiply(ds1) - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.multiply(ds1)", dt(steps)) - - self.assertLess(t1, t2) - - ds0 = ds - ds0.multiply(2, out=ds0) - steps.append(timer()) - print("ds0.multiply(2,out=ds0)", dt( - steps), 2., ds0.as_array()[0][0][0]) - self.assertEqual(2., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.multiply(2) - steps.append(timer()) - print("ds3 = ds0.multiply(2)", dt(steps), 4., ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - self.assertEqual(4., ds3.as_array()[0][0][0]) - self.assertEqual(2., ds.as_array()[0][0][0]) - - def binary_divide(self): - print("Test binary divide") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.divide(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.divide(ds1, out=ds)", dt(steps)) - steps.append(timer()) - ds2 = ds.divide(ds1) - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.divide(ds1)", dt(steps)) - - self.assertLess(t1, t2) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - ds0 = ds - ds0.divide(2, out=ds0) - steps.append(timer()) - print("ds0.divide(2,out=ds0)", dt(steps), 0.5, ds0.as_array()[0][0][0]) - self.assertEqual(0.5, ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.divide(2) - steps.append(timer()) - print("ds3 = ds0.divide(2)", dt(steps), 0.25, ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - self.assertEqual(.25, ds3.as_array()[0][0][0]) - self.assertEqual(.5, ds.as_array()[0][0][0]) - - def test_creation_copy(self): - shape = (2, 3, 4, 5) - size = shape[0] - for i in range(1, len(shape)): - size = size * shape[i] - #print("a refcount " , sys.getrefcount(a)) - a = numpy.asarray([i for i in range(size)]) - #print("a refcount " , sys.getrefcount(a)) - a = numpy.reshape(a, shape) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, True, ['X', 'Y', 'Z', 'W']) - #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a), 2) - - def test_subset(self): - shape = (4, 3, 2) - a = [i for i in range(2*3*4)] - a = numpy.asarray(a) - a = numpy.reshape(a, shape) - ds = DataContainer(a, True, ['X', 'Y', 'Z']) - sub = ds.subset(['X']) - res = True - try: - numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([0, 6, 12, 18])) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - sub = ds.subset(['X'], Y=2, Z=0) - res = True - try: - numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([4, 10, 16, 22])) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - sub = ds.subset(['Y']) - try: - numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([0, 2, 4])) - res = True - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - sub = ds.subset(['Z']) - try: - numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([0, 1])) - res = True - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - sub = ds.subset(['Z'], X=1, Y=2) - try: - numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([10, 11])) - res = True - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - print(a) - sub = ds.subset(['X', 'Y'], Z=1) - res = True - try: - numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([[1, 3, 5], - [7, 9, 11], - [13, 15, 17], - [19, 21, 23]])) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def test_ImageData(self): - # create ImageData from geometry - vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2) - vol = ImageData(geometry=vgeometry) - self.assertEqual(vol.shape, (2, 3, 4)) - - vol1 = vol + 1 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape)) - - vol1 = vol - 1 - self.assertNumpyArrayEqual(vol1.as_array(), -numpy.ones(vol.shape)) - - vol1 = 2 * (vol + 1) - self.assertNumpyArrayEqual(vol1.as_array(), 2 * numpy.ones(vol.shape)) - - vol1 = (vol + 1) / 2 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2) - - vol1 = vol + 1 - self.assertEqual(vol1.sum(), 2*3*4) - vol1 = (vol + 2) ** 2 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4) - - self.assertEqual(vol.number_of_dimensions, 3) - - def test_AcquisitionData(self): - sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), - geom_type='parallel', pixel_num_v=3, - pixel_num_h=5, channels=2) - sino = AcquisitionData(geometry=sgeometry) - self.assertEqual(sino.shape, (2, 10, 3, 5)) - - def assertNumpyArrayEqual(self, first, second): - res = True - try: - numpy.testing.assert_array_equal(first, second) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): - res = True - try: - numpy.testing.assert_array_almost_equal(first, second, decimal) - except AssertionError as err: - res = False - print(err) - print("expected " , second) - print("actual " , first) - - self.assertTrue(res) - def test_DataContainerChaining(self): - dc = self.create_DataContainer(256,256,256,1) - - dc.add(9,out=dc)\ - .subtract(1,out=dc) - self.assertEqual(1+9-1,dc.as_array().flatten()[0]) - def test_reduction(self): - print ("test reductions") - dc = self.create_DataContainer(2,2,2,value=1) - sqnorm = dc.squared_norm() - norm = dc.norm() - self.assertEqual(sqnorm, 8.0) - numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7) - - - -class TestAlgorithms(unittest.TestCase): - def assertNumpyArrayEqual(self, first, second): - res = True - try: - numpy.testing.assert_array_equal(first, second) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): - res = True - try: - numpy.testing.assert_array_almost_equal(first, second, decimal) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def test_FISTA_cvx(self): - if not cvx_not_installable: - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0], 1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt': True} - # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) - g0 = ZeroFun() - - # Initial guess - x_init = DataContainer(np.zeros((n, 1))) - - f.grad(x_init) - - # Run FISTA for least squares plus zero function. - x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt) - - # Print solution and final objective/criterion value for comparison - print("FISTA least squares plus zero function solution and objective value:") - print(x_fista0.array) - print(criter0[-1]) - - # Compare to CVXPY - - # Construct the problem. - x0 = Variable(n) - objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0])) - prob0 = Problem(objective0) - - # The optimal objective is returned by prob.solve(). - result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus zero function solution and objective value:") - print(x0.value) - print(objective0.value) - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fista0.array), x0.value, 6) - else: - self.assertTrue(cvx_not_installable) - - def test_FISTA_Norm1_cvx(self): - if not cvx_not_installable: - opt = {'memopt': True} - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0], 1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt': True} - # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) - g0 = ZeroFun() - - # Initial guess - x_init = DataContainer(np.zeros((n, 1))) - - # Create 1-norm object instance - g1 = Norm1(lam) - - g1(x_init) - g1.prox(x_init, 0.02) - - # Combine with least squares and solve using generic FISTA implementation - x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt) - - # Print for comparison - print("FISTA least squares plus 1-norm solution and objective value:") - print(x_fista1.as_array().squeeze()) - print(criter1[-1]) - - # Compare to CVXPY - - # Construct the problem. - x1 = Variable(n) - objective1 = Minimize( - 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) - prob1 = Problem(objective1) - - # The optimal objective is returned by prob.solve(). - result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1.value) - print(objective1.value) - - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fista1.array), x1.value, 6) - else: - self.assertTrue(cvx_not_installable) - - def skip_test_FBPD_Norm1_cvx(self): - print ("test_FBPD_Norm1_cvx") - if not cvx_not_installable: - opt = {'memopt': True} - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0], 1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt': True} - # Initial guess - x_init = DataContainer(np.random.randn(n, 1)) - - # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) - f.L = PowerMethodNonsquare(A, 25, x_init)[0] - print ("Lipschitz", f.L) - g0 = ZeroFun() - - - # Create 1-norm object instance - g1 = Norm1(lam) - - # Compare to CVXPY - - # Construct the problem. - x1 = Variable(n) - objective1 = Minimize( - 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) - prob1 = Problem(objective1) - - # The optimal objective is returned by prob.solve(). - result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1.value) - print(objective1.value) - - # Now try another algorithm FBPD for same problem: - x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, - Identity(), None, f, g1) - print(x_fbpd1) - print(criterfbpd1[-1]) - - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fbpd1.array), x1.value, 6) - else: - self.assertTrue(cvx_not_installable) - # Plot criterion curve to see both FISTA and FBPD converge to same value. - # Note that FISTA is very efficient for 1-norm minimization so it beats - # FBPD in this test by a lot. But FBPD can handle a larger class of problems - # than FISTA can. - - # Now try 1-norm and TV denoising with FBPD, first 1-norm. - - # Set up phantom size NxN by creating ImageGeometry, initialising the - # ImageData object with this geometry and empty array and finally put some - # data into its array, and display as image. - def skip_test_FISTA_denoise_cvx(self): - if not cvx_not_installable: - opt = {'memopt': True} - N = 64 - ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) - Phantom = ImageData(geometry=ig) - - x = Phantom.as_array() - - x[int(round(N/4)):int(round(3*N/4)), - int(round(N/4)):int(round(3*N/4))] = 0.5 - x[int(round(N/8)):int(round(7*N/8)), - int(round(3*N/8)):int(round(5*N/8))] = 1 - - # Identity operator for denoising - I = TomoIdentity(ig) - - # Data and add noise - y = I.direct(Phantom) - y.array = y.array + 0.1*np.random.randn(N, N) - - # Data fidelity term - f_denoise = Norm2sq(I, y, c=0.5, memopt=True) - x_init = ImageData(geometry=ig) - f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0] - - # 1-norm regulariser - lam1_denoise = 1.0 - g1_denoise = Norm1(lam1_denoise) - - # Initial guess - x_init_denoise = ImageData(np.zeros((N, N))) - - # Combine with least squares and solve using generic FISTA implementation - x_fista1_denoise, it1_denoise, timing1_denoise, \ - criter1_denoise = \ - FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) - - print(x_fista1_denoise) - print(criter1_denoise[-1]) - - # Now denoise LS + 1-norm with FBPD - x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\ - criterfbpd1_denoise = \ - FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) - print(x_fbpd1_denoise) - print(criterfbpd1_denoise[-1]) - - # Compare to CVXPY - - # Construct the problem. - x1_denoise = Variable(N**2) - objective1_denoise = Minimize( - 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1)) - prob1_denoise = Problem(objective1_denoise) - - # The optimal objective is returned by prob.solve(). - result1_denoise = prob1_denoise.solve( - verbose=False, solver=SCS, eps=1e-12) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1_denoise.value) - print(objective1_denoise.value) - self.assertNumpyArrayAlmostEqual( - x_fista1_denoise.array.flatten(), x1_denoise.value, 5) - - self.assertNumpyArrayAlmostEqual( - x_fbpd1_denoise.array.flatten(), x1_denoise.value, 5) - x1_cvx = x1_denoise.value - x1_cvx.shape = (N, N) - - # Now TV with FBPD - lam_tv = 0.1 - gtv = TV2D(lam_tv) - gtv(gtv.op.direct(x_init_denoise)) - - opt_tv = {'tol': 1e-4, 'iter': 10000} - - x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\ - criterfbpdtv_denoise = \ - FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv, opt=opt_tv) - print(x_fbpdtv_denoise) - print(criterfbpdtv_denoise[-1]) - - # Compare to CVXPY - - # Construct the problem. - xtv_denoise = Variable((N, N)) - objectivetv_denoise = Minimize( - 0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise)) - probtv_denoise = Problem(objectivetv_denoise) - - # The optimal objective is returned by prob.solve(). - resulttv_denoise = probtv_denoise.solve( - verbose=False, solver=SCS, eps=1e-12) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(xtv_denoise.value) - print(objectivetv_denoise.value) - - self.assertNumpyArrayAlmostEqual( - x_fbpdtv_denoise.as_array(), xtv_denoise.value, 1) - - else: - self.assertTrue(cvx_not_installable) - - -class TestFunction(unittest.TestCase): - def assertNumpyArrayEqual(self, first, second): - res = True - try: - numpy.testing.assert_array_equal(first, second) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def create_simple_ImageData(self): - N = 64 - ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) - Phantom = ImageData(geometry=ig) - - x = Phantom.as_array() - - x[int(round(N/4)):int(round(3*N/4)), - int(round(N/4)):int(round(3*N/4))] = 0.5 - x[int(round(N/8)):int(round(7*N/8)), - int(round(3*N/8)):int(round(5*N/8))] = 1 - - return (ig, Phantom) - - def _test_Norm2(self): - print("test Norm2") - opt = {'memopt': True} - ig, Phantom = self.create_simple_ImageData() - x = Phantom.as_array() - print(Phantom) - print(Phantom.as_array()) - - norm2 = Norm2() - v1 = norm2(x) - v2 = norm2(Phantom) - self.assertEqual(v1, v2) - - p1 = norm2.prox(Phantom, 1) - print(p1) - p2 = norm2.prox(x, 1) - self.assertNumpyArrayEqual(p1.as_array(), p2) - - p3 = norm2.proximal(Phantom, 1) - p4 = norm2.proximal(x, 1) - self.assertNumpyArrayEqual(p3.as_array(), p2) - self.assertNumpyArrayEqual(p3.as_array(), p4) - - out = Phantom.copy() - p5 = norm2.proximal(Phantom, 1, out=out) - self.assertEqual(id(p5), id(out)) - self.assertNumpyArrayEqual(p5.as_array(), p3.as_array()) -# ============================================================================= -# def test_upper(self): -# self.assertEqual('foo'.upper(), 'FOO') -# -# def test_isupper(self): -# self.assertTrue('FOO'.isupper()) -# self.assertFalse('Foo'.isupper()) -# -# def test_split(self): -# s = 'hello world' -# self.assertEqual(s.split(), ['hello', 'world']) -# # check that s.split fails when the separator is not a string -# with self.assertRaises(TypeError): -# s.split(2) -# ============================================================================= - -class TestNexusReader(unittest.TestCase): - - def setUp(self): - wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') - self.filename = '24737_fd.nxs' - - def tearDown(self): - os.remove(self.filename) - - - def testGetDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct") - - def testGetProjectionDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct") - - def testLoadProjectionWithoutDimensions(self): - nr = NexusReader(self.filename) - projections = nr.load_projection() - self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionWithDimensions(self): - nr = NexusReader(self.filename) - projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) - self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionCompareSingle(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) - - def testLoadProjectionCompareMulti(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) - - def testLoadProjectionCompareRandom(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) - numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) - - def testLoadProjectionCompareFull(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) - - def testLoadFlatCompareFull(self): - nr = NexusReader(self.filename) - flats_full = nr.load_flat() - flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) - - def testLoadDarkCompareFull(self): - nr = NexusReader(self.filename) - darks_full = nr.load_dark() - darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) - - def testProjectionAngles(self): - nr = NexusReader(self.filename) - angles = nr.get_projection_angles() - self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") - - def test_get_acquisition_data_subset(self): - nr = NexusReader(self.filename) - key = nr.get_image_keys() - sl = nr.get_acquisition_data_subset(0,10) - data = nr.get_acquisition_data().subset(['vertical','horizontal']) - - self.assertTrue(sl.shape , (10,data.shape[1])) - - -if __name__ == '__main__': - unittest.main() - diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index b584344..eaf124b 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -23,12 +23,16 @@ import os import sys -cil_version='0.11.3' +cil_version=os.environ['CIL_VERSION'] +if cil_version == '': + print("Please set the environmental variable CIL_VERSION") + sys.exit(1) setup( name="ccpi-framework", version=cil_version, - packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation'], + packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation', + 'ccpi.optimisation.algorithms'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine diff --git a/Wrappers/Python/test/TestReader.py b/Wrappers/Python/test/TestReader.py deleted file mode 100644 index 51db052..0000000 --- a/Wrappers/Python/test/TestReader.py +++ /dev/null @@ -1,134 +0,0 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -''' -Unit tests for Readers - -@author: Mr. Srikanth Nagella -''' -import unittest - -import numpy.testing -import wget -import os -from ccpi.io.reader import NexusReader -from ccpi.io.reader import XTEKReader -#@unittest.skip -class TestNexusReader(unittest.TestCase): - - def setUp(self): - wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') - self.filename = '24737_fd.nxs' - - def tearDown(self): - os.remove(self.filename) - - - def testGetDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.getSinogramDimensions(), (135, 91, 160), "Sinogram dimensions are not correct") - - def testGetProjectionDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.getProjectionDimensions(), (91,135,160), "Projection dimensions are not correct") - - def testLoadProjectionWithoutDimensions(self): - nr = NexusReader(self.filename) - projections = nr.loadProjection() - self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionWithDimensions(self): - nr = NexusReader(self.filename) - projections = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) - self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionCompareSingle(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) - - def testLoadProjectionCompareMulti(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(0,3), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) - - def testLoadProjectionCompareRandom(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(1,8), slice(5,10), slice(8,20))) - numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) - - def testLoadProjectionCompareFull(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) - - def testLoadFlatCompareFull(self): - nr = NexusReader(self.filename) - flats_full = nr.loadFlat() - flats_part = nr.loadFlat((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) - - def testLoadDarkCompareFull(self): - nr = NexusReader(self.filename) - darks_full = nr.loadDark() - darks_part = nr.loadDark((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) - - def testProjectionAngles(self): - nr = NexusReader(self.filename) - angles = nr.getProjectionAngles() - self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") - -class TestXTEKReader(unittest.TestCase): - - def setUp(self): - testpath, filename = os.path.split(os.path.realpath(__file__)) - testpath = os.path.normpath(os.path.join(testpath, '..','..','..')) - self.filename = os.path.join(testpath,'data','SophiaBeads','SophiaBeads_64_averaged.xtekct') - - def tearDown(self): - pass - - def testLoad(self): - xtekReader = XTEKReader(self.filename) - self.assertEqual(xtekReader.geometry.pixel_num_h, 500, "Detector pixel X size is not correct") - self.assertEqual(xtekReader.geometry.pixel_num_v, 500, "Detector pixel Y size is not correct") - self.assertEqual(xtekReader.geometry.dist_source_center, -80.6392412185669, "Distance from source to center is not correct") - self.assertEqual(xtekReader.geometry.dist_center_detector, (1007.006 - 80.6392412185669), "Distance from center to detector is not correct") - - def testReadAngles(self): - xtekReader = XTEKReader(self.filename) - angles = xtekReader.readAngles() - self.assertEqual(angles.shape, (63,), "Angles doesn't match") - self.assertAlmostEqual(angles[46], -085.717, 3, "46th Angles doesn't match") - - def testLoadProjection(self): - xtekReader = XTEKReader(self.filename) - pixels = xtekReader.loadProjection() - self.assertEqual(pixels.shape, (63, 500, 500), "projections shape doesn't match") - - - -if __name__ == "__main__": - #import sys;sys.argv = ['', 'TestNexusReader.testLoad'] - unittest.main() \ No newline at end of file diff --git a/Wrappers/Python/test/__init__.py b/Wrappers/Python/test/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/Wrappers/Python/test/skip_TestReader.py b/Wrappers/Python/test/skip_TestReader.py new file mode 100644 index 0000000..ec7ed58 --- /dev/null +++ b/Wrappers/Python/test/skip_TestReader.py @@ -0,0 +1,134 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +''' +Unit tests for Readers + +@author: Mr. Srikanth Nagella +''' +import unittest + +import numpy.testing +import wget +import os +from ccpi.io.reader import NexusReader +from ccpi.io.reader import XTEKReader +#@unittest.skip +class TestNexusReader(unittest.TestCase): + + def setUp(self): + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' + + def tearDown(self): + os.remove(self.filename) + + + def testGetDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.getSinogramDimensions(), (135, 91, 160), "Sinogram dimensions are not correct") + + def testGetProjectionDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.getProjectionDimensions(), (91,135,160), "Projection dimensions are not correct") + + def testLoadProjectionWithoutDimensions(self): + nr = NexusReader(self.filename) + projections = nr.loadProjection() + self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionWithDimensions(self): + nr = NexusReader(self.filename) + projections = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) + self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionCompareSingle(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) + + def testLoadProjectionCompareMulti(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(0,3), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) + + def testLoadProjectionCompareRandom(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(1,8), slice(5,10), slice(8,20))) + numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) + + def testLoadProjectionCompareFull(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) + + def testLoadFlatCompareFull(self): + nr = NexusReader(self.filename) + flats_full = nr.loadFlat() + flats_part = nr.loadFlat((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) + + def testLoadDarkCompareFull(self): + nr = NexusReader(self.filename) + darks_full = nr.loadDark() + darks_part = nr.loadDark((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) + + def testProjectionAngles(self): + nr = NexusReader(self.filename) + angles = nr.getProjectionAngles() + self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") + +class TestXTEKReader(unittest.TestCase): + + def setUp(self): + testpath, filename = os.path.split(os.path.realpath(__file__)) + testpath = os.path.normpath(os.path.join(testpath, '..','..','..')) + self.filename = os.path.join(testpath,'data','SophiaBeads','SophiaBeads_64_averaged.xtekct') + + def tearDown(self): + pass + + def testLoad(self): + xtekReader = XTEKReader(self.filename) + self.assertEqual(xtekReader.geometry.pixel_num_h, 500, "Detector pixel X size is not correct") + self.assertEqual(xtekReader.geometry.pixel_num_v, 500, "Detector pixel Y size is not correct") + self.assertEqual(xtekReader.geometry.dist_source_center, -80.6392412185669, "Distance from source to center is not correct") + self.assertEqual(xtekReader.geometry.dist_center_detector, (1007.006 - 80.6392412185669), "Distance from center to detector is not correct") + + def testReadAngles(self): + xtekReader = XTEKReader(self.filename) + angles = xtekReader.readAngles() + self.assertEqual(angles.shape, (63,), "Angles doesn't match") + self.assertAlmostEqual(angles[46], -085.717, 3, "46th Angles doesn't match") + + def testLoadProjection(self): + xtekReader = XTEKReader(self.filename) + pixels = xtekReader.loadProjection() + self.assertEqual(pixels.shape, (63, 500, 500), "projections shape doesn't match") + + + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'TestNexusReader.testLoad'] + unittest.main() \ No newline at end of file diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py new file mode 100755 index 0000000..05f3fe8 --- /dev/null +++ b/Wrappers/Python/test/test_DataContainer.py @@ -0,0 +1,499 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:08:21 2019 + +@author: ofn77899 +""" +import sys +import unittest +import numpy +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from timeit import default_timer as timer + + +def dt(steps): + return steps[-1] - steps[-2] +def aid(x): + # This function returns the memory + # block address of an array. + return x.__array_interface__['data'][0] + + + +class TestDataContainer(unittest.TestCase): + def create_simple_ImageData(self): + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + return (ig, Phantom) + + def create_DataContainer(self, X,Y,Z, value=1): + steps = [timer()] + a = value * numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + return ds + + def test_creation_nocopy(self): + shape = (2, 3, 4, 5) + size = shape[0] + for i in range(1, len(shape)): + size = size * shape[i] + #print("a refcount " , sys.getrefcount(a)) + a = numpy.asarray([i for i in range(size)]) + #print("a refcount " , sys.getrefcount(a)) + a = numpy.reshape(a, shape) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 3) + self.assertEqual(ds.dimension_labels, {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'}) + + def testGb_creation_nocopy(self): + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print("test clone") + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 3) + ds1 = ds.copy() + self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) + ds1 = ds.clone() + self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) + + def testInlineAlgebra(self): + print("Test Inline Algebra") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print(t0) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + #ds.__iadd__( 2 ) + ds += 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 3.) + #ds.__isub__( 2 ) + ds -= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + #ds.__imul__( 2 ) + ds *= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 2.) + #ds.__idiv__( 2 ) + ds /= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds1 = ds.copy() + #ds1.__iadd__( 1 ) + ds1 += 1 + #ds.__iadd__( ds1 ) + ds += ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 3.) + #ds.__isub__( ds1 ) + ds -= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + #ds.__imul__( ds1 ) + ds *= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 2.) + #ds.__idiv__( ds1 ) + ds /= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + def test_unary_operations(self): + print("Test unary operations") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = -numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print(t0) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + + ds.sign(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], -1.) + + ds.abs(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds.__imul__(2) + ds.sqrt(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], + numpy.sqrt(2., dtype='float32')) + + def test_binary_operations(self): + self.binary_add() + self.binary_subtract() + self.binary_multiply() + self.binary_divide() + + def binary_add(self): + print("Test binary add") + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.add(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.add(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.add(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.add(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 2.) + + ds0 = ds + ds0.add(2, out=ds0) + steps.append(timer()) + print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0]) + self.assertEqual(4., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.add(2) + steps.append(timer()) + print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + + def binary_subtract(self): + print("Test binary subtract") + X, Y, Z = 512, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.subtract(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.subtract(ds1, out=ds)", dt(steps)) + self.assertEqual(0., ds.as_array()[0][0][0]) + + steps.append(timer()) + ds2 = ds.subtract(ds1) + self.assertEqual(-1., ds2.as_array()[0][0][0]) + + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.subtract(ds1)", dt(steps)) + + self.assertLess(t1, t2) + + del ds1 + ds0 = ds.copy() + steps.append(timer()) + ds0.subtract(2, out=ds0) + #ds0.__isub__( 2 ) + steps.append(timer()) + print("ds0.subtract(2,out=ds0)", dt( + steps), -2., ds0.as_array()[0][0][0]) + self.assertEqual(-2., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.subtract(2) + steps.append(timer()) + print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(-2., ds0.as_array()[0][0][0]) + self.assertEqual(-4., ds3.as_array()[0][0][0]) + + def binary_multiply(self): + print("Test binary multiply") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.multiply(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.multiply(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.multiply(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.multiply(ds1)", dt(steps)) + + self.assertLess(t1, t2) + + ds0 = ds + ds0.multiply(2, out=ds0) + steps.append(timer()) + print("ds0.multiply(2,out=ds0)", dt( + steps), 2., ds0.as_array()[0][0][0]) + self.assertEqual(2., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.multiply(2) + steps.append(timer()) + print("ds3 = ds0.multiply(2)", dt(steps), 4., ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(4., ds3.as_array()[0][0][0]) + self.assertEqual(2., ds.as_array()[0][0][0]) + + def binary_divide(self): + print("Test binary divide") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.divide(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.divide(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.divide(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.divide(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds0 = ds + ds0.divide(2, out=ds0) + steps.append(timer()) + print("ds0.divide(2,out=ds0)", dt(steps), 0.5, ds0.as_array()[0][0][0]) + self.assertEqual(0.5, ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.divide(2) + steps.append(timer()) + print("ds3 = ds0.divide(2)", dt(steps), 0.25, ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(.25, ds3.as_array()[0][0][0]) + self.assertEqual(.5, ds.as_array()[0][0][0]) + + def test_creation_copy(self): + shape = (2, 3, 4, 5) + size = shape[0] + for i in range(1, len(shape)): + size = size * shape[i] + #print("a refcount " , sys.getrefcount(a)) + a = numpy.asarray([i for i in range(size)]) + #print("a refcount " , sys.getrefcount(a)) + a = numpy.reshape(a, shape) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, True, ['X', 'Y', 'Z', 'W']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 2) + + def test_subset(self): + shape = (4, 3, 2) + a = [i for i in range(2*3*4)] + a = numpy.asarray(a) + a = numpy.reshape(a, shape) + ds = DataContainer(a, True, ['X', 'Y', 'Z']) + sub = ds.subset(['X']) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([0, 6, 12, 18])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['X'], Y=2, Z=0) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([4, 10, 16, 22])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['Y']) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([0, 2, 4])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['Z']) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([0, 1])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + sub = ds.subset(['Z'], X=1, Y=2) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([10, 11])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + print(a) + sub = ds.subset(['X', 'Y'], Z=1) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([[1, 3, 5], + [7, 9, 11], + [13, 15, 17], + [19, 21, 23]])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def test_ImageData(self): + # create ImageData from geometry + vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2) + vol = ImageData(geometry=vgeometry) + self.assertEqual(vol.shape, (2, 3, 4)) + + vol1 = vol + 1 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape)) + + vol1 = vol - 1 + self.assertNumpyArrayEqual(vol1.as_array(), -numpy.ones(vol.shape)) + + vol1 = 2 * (vol + 1) + self.assertNumpyArrayEqual(vol1.as_array(), 2 * numpy.ones(vol.shape)) + + vol1 = (vol + 1) / 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2) + + vol1 = vol + 1 + self.assertEqual(vol1.sum(), 2*3*4) + vol1 = (vol + 2) ** 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4) + + self.assertEqual(vol.number_of_dimensions, 3) + + def test_AcquisitionData(self): + sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5, channels=2) + sino = AcquisitionData(geometry=sgeometry) + self.assertEqual(sino.shape, (2, 10, 3, 5)) + + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + print("expected " , second) + print("actual " , first) + + self.assertTrue(res) + def test_DataContainerChaining(self): + dc = self.create_DataContainer(256,256,256,1) + + dc.add(9,out=dc)\ + .subtract(1,out=dc) + self.assertEqual(1+9-1,dc.as_array().flatten()[0]) + def test_reduction(self): + print ("test reductions") + dc = self.create_DataContainer(2,2,2,value=1) + sqnorm = dc.squared_norm() + norm = dc.norm() + self.assertEqual(sqnorm, 8.0) + numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7) + + + +if __name__ == '__main__': + unittest.main() + \ No newline at end of file diff --git a/Wrappers/Python/test/test_NexusReader.py b/Wrappers/Python/test/test_NexusReader.py new file mode 100755 index 0000000..55543ba --- /dev/null +++ b/Wrappers/Python/test/test_NexusReader.py @@ -0,0 +1,96 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:05:00 2019 + +@author: ofn77899 +""" + +import unittest +import wget +import os +from ccpi.io.reader import NexusReader +import numpy + + +class TestNexusReader(unittest.TestCase): + + def setUp(self): + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' + + def tearDown(self): + os.remove(self.filename) + + + def testGetDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct") + + def testGetProjectionDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct") + + def testLoadProjectionWithoutDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection() + self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionWithDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionCompareSingle(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) + + def testLoadProjectionCompareMulti(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) + + def testLoadProjectionCompareRandom(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) + numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) + + def testLoadProjectionCompareFull(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) + + def testLoadFlatCompareFull(self): + nr = NexusReader(self.filename) + flats_full = nr.load_flat() + flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) + + def testLoadDarkCompareFull(self): + nr = NexusReader(self.filename) + darks_full = nr.load_dark() + darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) + + def testProjectionAngles(self): + nr = NexusReader(self.filename) + angles = nr.get_projection_angles() + self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") + + def test_get_acquisition_data_subset(self): + nr = NexusReader(self.filename) + key = nr.get_image_keys() + sl = nr.get_acquisition_data_subset(0,10) + data = nr.get_acquisition_data().subset(['vertical','horizontal']) + + self.assertTrue(sl.shape , (10,data.shape[1])) + + + +if __name__ == '__main__': + unittest.main() + diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py new file mode 100755 index 0000000..b5959b5 --- /dev/null +++ b/Wrappers/Python/test/test_algorithms.py @@ -0,0 +1,123 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:11:36 2019 + +@author: ofn77899 +""" + +import unittest +import numpy +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from ccpi.optimisation.ops import TomoIdentity +from ccpi.optimisation.funcs import Norm2sq +from ccpi.optimisation.algorithms import GradientDescent +from ccpi.optimisation.algorithms import CGLS +from ccpi.optimisation.algorithms import FISTA +from ccpi.optimisation.algorithms import FBPD + + + + +class TestAlgorithms(unittest.TestCase): + def setUp(self): + #wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + #self.filename = '24737_fd.nxs' + # we use TomoIdentity as the operator and solve the simple least squares + # problem for a random-valued ImageData or AcquisitionData b? + # Then we know the minimiser is b itself + + # || I x -b ||^2 + + # create an ImageGeometry + ig = ImageGeometry(12,13,14) + pass + + def tearDown(self): + #os.remove(self.filename) + pass + + def test_GradientDescent(self): + print ("Test GradientDescent") + ig = ImageGeometry(12,13,14) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + norm2sq = Norm2sq(identity, b) + + alg = GradientDescent(x_init=x_init, + objective_function=norm2sq, + rate=0.3) + alg.max_iteration = 20 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + def test_CGLS(self): + print ("Test CGLS") + ig = ImageGeometry(124,153,154) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + alg = CGLS(x_init=x_init, operator=identity, data=b) + alg.max_iteration = 1 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + + def test_FISTA(self): + print ("Test FISTA") + ig = ImageGeometry(127,139,149) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + x_init = ImageData(geometry=ig) + x_init.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + norm2sq = Norm2sq(identity, b) + opt = {'tol': 1e-4, 'memopt':False} + alg = FISTA(x_init=x_init, f=norm2sq, g=None, opt=opt) + alg.max_iteration = 2 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + + + + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + + + + +if __name__ == '__main__': + unittest.main() + \ No newline at end of file diff --git a/Wrappers/Python/test/test_run_test.py b/Wrappers/Python/test/test_run_test.py new file mode 100755 index 0000000..d0b87f5 --- /dev/null +++ b/Wrappers/Python/test/test_run_test.py @@ -0,0 +1,435 @@ +import unittest +import numpy +import numpy as np +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from ccpi.optimisation.algs import FISTA +from ccpi.optimisation.algs import FBPD +from ccpi.optimisation.funcs import Norm2sq +from ccpi.optimisation.funcs import ZeroFun +from ccpi.optimisation.funcs import Norm1 +from ccpi.optimisation.funcs import TV2D +from ccpi.optimisation.funcs import Norm2 + +from ccpi.optimisation.ops import LinearOperatorMatrix +from ccpi.optimisation.ops import TomoIdentity +from ccpi.optimisation.ops import Identity +from ccpi.optimisation.ops import PowerMethodNonsquare + + +import numpy.testing + +try: + from cvxpy import * + cvx_not_installable = False +except ImportError: + cvx_not_installable = True + + +def aid(x): + # This function returns the memory + # block address of an array. + return x.__array_interface__['data'][0] + + +def dt(steps): + return steps[-1] - steps[-2] + + + + +class TestAlgorithms(unittest.TestCase): + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def test_FISTA_cvx(self): + if not cvx_not_installable: + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + f.grad(x_init) + + # Run FISTA for least squares plus zero function. + x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt) + + # Print solution and final objective/criterion value for comparison + print("FISTA least squares plus zero function solution and objective value:") + print(x_fista0.array) + print(criter0[-1]) + + # Compare to CVXPY + + # Construct the problem. + x0 = Variable(n) + objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0])) + prob0 = Problem(objective0) + + # The optimal objective is returned by prob.solve(). + result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus zero function solution and objective value:") + print(x0.value) + print(objective0.value) + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista0.array), x0.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def test_FISTA_Norm1_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + # Create 1-norm object instance + g1 = Norm1(lam) + + g1(x_init) + g1.prox(x_init, 0.02) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt) + + # Print for comparison + print("FISTA least squares plus 1-norm solution and objective value:") + print(x_fista1.as_array().squeeze()) + print(criter1[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def skip_test_FBPD_Norm1_cvx(self): + print ("test_FBPD_Norm1_cvx") + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Initial guess + x_init = DataContainer(np.random.randn(n, 1)) + + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + f.L = PowerMethodNonsquare(A, 25, x_init)[0] + print ("Lipschitz", f.L) + g0 = ZeroFun() + + + # Create 1-norm object instance + g1 = Norm1(lam) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + # Now try another algorithm FBPD for same problem: + x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, + Identity(), None, f, g1) + print(x_fbpd1) + print(criterfbpd1[-1]) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fbpd1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) + # Plot criterion curve to see both FISTA and FBPD converge to same value. + # Note that FISTA is very efficient for 1-norm minimization so it beats + # FBPD in this test by a lot. But FBPD can handle a larger class of problems + # than FISTA can. + + # Now try 1-norm and TV denoising with FBPD, first 1-norm. + + # Set up phantom size NxN by creating ImageGeometry, initialising the + # ImageData object with this geometry and empty array and finally put some + # data into its array, and display as image. + def skip_test_FISTA_denoise_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + # Identity operator for denoising + I = TomoIdentity(ig) + + # Data and add noise + y = I.direct(Phantom) + y.array = y.array + 0.1*np.random.randn(N, N) + + # Data fidelity term + f_denoise = Norm2sq(I, y, c=0.5, memopt=True) + x_init = ImageData(geometry=ig) + f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0] + + # 1-norm regulariser + lam1_denoise = 1.0 + g1_denoise = Norm1(lam1_denoise) + + # Initial guess + x_init_denoise = ImageData(np.zeros((N, N))) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1_denoise, it1_denoise, timing1_denoise, \ + criter1_denoise = \ + FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) + + print(x_fista1_denoise) + print(criter1_denoise[-1]) + + # Now denoise LS + 1-norm with FBPD + x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\ + criterfbpd1_denoise = \ + FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) + print(x_fbpd1_denoise) + print(criterfbpd1_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1_denoise = Variable(N**2) + objective1_denoise = Minimize( + 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1)) + prob1_denoise = Problem(objective1_denoise) + + # The optimal objective is returned by prob.solve(). + result1_denoise = prob1_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1_denoise.value) + print(objective1_denoise.value) + self.assertNumpyArrayAlmostEqual( + x_fista1_denoise.array.flatten(), x1_denoise.value, 5) + + self.assertNumpyArrayAlmostEqual( + x_fbpd1_denoise.array.flatten(), x1_denoise.value, 5) + x1_cvx = x1_denoise.value + x1_cvx.shape = (N, N) + + # Now TV with FBPD + lam_tv = 0.1 + gtv = TV2D(lam_tv) + gtv(gtv.op.direct(x_init_denoise)) + + opt_tv = {'tol': 1e-4, 'iter': 10000} + + x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\ + criterfbpdtv_denoise = \ + FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv, opt=opt_tv) + print(x_fbpdtv_denoise) + print(criterfbpdtv_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + xtv_denoise = Variable((N, N)) + objectivetv_denoise = Minimize( + 0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise)) + probtv_denoise = Problem(objectivetv_denoise) + + # The optimal objective is returned by prob.solve(). + resulttv_denoise = probtv_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(xtv_denoise.value) + print(objectivetv_denoise.value) + + self.assertNumpyArrayAlmostEqual( + x_fbpdtv_denoise.as_array(), xtv_denoise.value, 1) + + else: + self.assertTrue(cvx_not_installable) + + +class TestFunction(unittest.TestCase): + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def create_simple_ImageData(self): + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + return (ig, Phantom) + + def _test_Norm2(self): + print("test Norm2") + opt = {'memopt': True} + ig, Phantom = self.create_simple_ImageData() + x = Phantom.as_array() + print(Phantom) + print(Phantom.as_array()) + + norm2 = Norm2() + v1 = norm2(x) + v2 = norm2(Phantom) + self.assertEqual(v1, v2) + + p1 = norm2.prox(Phantom, 1) + print(p1) + p2 = norm2.prox(x, 1) + self.assertNumpyArrayEqual(p1.as_array(), p2) + + p3 = norm2.proximal(Phantom, 1) + p4 = norm2.proximal(x, 1) + self.assertNumpyArrayEqual(p3.as_array(), p2) + self.assertNumpyArrayEqual(p3.as_array(), p4) + + out = Phantom.copy() + p5 = norm2.proximal(Phantom, 1, out=out) + self.assertEqual(id(p5), id(out)) + self.assertNumpyArrayEqual(p5.as_array(), p3.as_array()) +# ============================================================================= +# def test_upper(self): +# self.assertEqual('foo'.upper(), 'FOO') +# +# def test_isupper(self): +# self.assertTrue('FOO'.isupper()) +# self.assertFalse('Foo'.isupper()) +# +# def test_split(self): +# s = 'hello world' +# self.assertEqual(s.split(), ['hello', 'world']) +# # check that s.split fails when the separator is not a string +# with self.assertRaises(TypeError): +# s.split(2) +# ============================================================================= + + + +if __name__ == '__main__': + unittest.main() + -- cgit v1.2.3