summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorEdoardo Pasca <edo.paskino@gmail.com>2019-02-28 18:58:33 +0000
committerGitHub <noreply@github.com>2019-02-28 18:58:33 +0000
commit13054c64cff2b506a7affa9d21445f368abbb15b (patch)
treea0101b07c8ceebb7ccce8884158c8573f48c7c34
parent0abb52cd5929d809b3ac1a651d85f465ada80137 (diff)
downloadframework-13054c64cff2b506a7affa9d21445f368abbb15b.tar.gz
framework-13054c64cff2b506a7affa9d21445f368abbb15b.tar.bz2
framework-13054c64cff2b506a7affa9d21445f368abbb15b.tar.xz
framework-13054c64cff2b506a7affa9d21445f368abbb15b.zip
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
Diffstat
-rw-r--r--.gitattributes1
-rw-r--r--Wrappers/Python/ccpi/framework.py12
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/Algorithm.py157
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/CGLS.py87
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/FBPD.py86
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/FISTA.py121
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py72
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/__init__.py29
-rwxr-xr-xWrappers/Python/ccpi/optimisation/ops.py48
-rw-r--r--Wrappers/Python/conda-recipe/meta.yaml9
-rwxr-xr-xWrappers/Python/conda-recipe/run_test.py985
-rw-r--r--Wrappers/Python/setup.py8
-rw-r--r--Wrappers/Python/test/__init__.py0
-rw-r--r--Wrappers/Python/test/skip_TestReader.py (renamed from Wrappers/Python/test/TestReader.py)266
-rwxr-xr-xWrappers/Python/test/test_DataContainer.py499
-rwxr-xr-xWrappers/Python/test/test_NexusReader.py96
-rwxr-xr-xWrappers/Python/test/test_algorithms.py123
-rwxr-xr-xWrappers/Python/test/test_run_test.py435
18 files changed, 1903 insertions, 1131 deletions
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..d9bd16b
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1 @@
+*.py text eol=lf
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/__init__.py b/Wrappers/Python/test/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/Wrappers/Python/test/__init__.py
diff --git a/Wrappers/Python/test/TestReader.py b/Wrappers/Python/test/skip_TestReader.py
index 51db052..ec7ed58 100644
--- a/Wrappers/Python/test/TestReader.py
+++ b/Wrappers/Python/test/skip_TestReader.py
@@ -1,134 +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']
+# -*- 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()
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-04 03:12:58 +0000