Merge branch 'master' into composite_operator_datacontainer

21 files changed, 1889 insertions, 1294 deletions

diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index d1ad26b..dab2dd9 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -3,7 +3,7 @@
 #   Visual Analytics and Imaging System Group of the Science Technology
 #   Facilities Council, STFC
 
-#   Copyright 2018 Edoardo Pasca
+#   Copyright 2018-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.
@@ -26,6 +26,7 @@ import numpy
 import sys
 from datetime import timedelta, datetime
 import warnings
+from functools import reduce
 
 def find_key(dic, val):
     """return the key of dictionary dic given the value"""
@@ -729,6 +730,16 @@ class DataContainer(object):
     ## reductions
     def sum(self, out=None, *args, **kwargs):
         return self.as_array().sum(*args, **kwargs)
+    def squared_norm(self):
+        '''return the squared euclidean norm of the DataContainer viewed as a vector'''
+        shape = self.shape
+        size = reduce(lambda x,y:x*y, shape, 1)
+        y = numpy.reshape(self.as_array(), (size, ))
+        return numpy.dot(y, y.conjugate())
+    def norm(self):
+        '''return the euclidean norm of the DataContainer viewed as a vector'''
+        return numpy.sqrt(self.squared_norm())
+    
     
     
 class ImageData(DataContainer):
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/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py
index a736277..15638a9 100755
--- a/Wrappers/Python/ccpi/optimisation/algs.py
+++ b/Wrappers/Python/ccpi/optimisation/algs.py
@@ -158,9 +158,8 @@ def FBPD(x_init, operator=None, constraint=None, data_fidelity=None,\
     memopt = opt['memopts'] if 'memopts' in opt.keys() else False
     
     # step-sizes
-    tau   = 2 / (data_fidelity.L + 2);
-    sigma = (1/tau - data_fidelity.L/2) / regulariser.L;
-    
+    tau   = 2 / (data_fidelity.L + 2)
+    sigma = (1/tau - data_fidelity.L/2) / regulariser.L
     inv_sigma = 1/sigma
 
     # initialization
diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py
index c7a6143..9b9fc36 100755
--- a/Wrappers/Python/ccpi/optimisation/funcs.py
+++ b/Wrappers/Python/ccpi/optimisation/funcs.py
@@ -217,10 +217,10 @@ class ZeroFun(Function):
 class Norm1(Function):
     
     def __init__(self,gamma):
+        super(Norm1, self).__init__()
         self.gamma = gamma
         self.L = 1
         self.sign_x = None
-        super(Norm1, self).__init__()
     
     def __call__(self,x,out=None):
         if out is None:
diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py
index a0e1713..54cebcd 100755
--- a/Wrappers/Python/ccpi/optimisation/ops.py
+++ b/Wrappers/Python/ccpi/optimisation/ops.py
@@ -48,9 +48,15 @@ class Operator(object):
         raise NotImplementedError
     def allocate_adjoint(self):
         '''Allocates memory on the X space'''
+<<<<<<< HEAD
         raise NotImplementedError
     def range_dim(self):
         raise NotImplementedError
+=======
+        raise NotImplementedError
+    def range_dim(self):
+        raise NotImplementedError
+>>>>>>> master
     def domain_dim(self):
         raise NotImplementedError
     def __rmul__(self, other):
@@ -93,9 +99,10 @@ 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):
         
@@ -236,10 +243,11 @@ def PowerMethodNonsquare(op,numiters , x0=None):
     # Loop
     for it in numpy.arange(numiters):
         x1 = op.adjoint(op.direct(x0))
-        x1norm = numpy.sqrt((x1*x1).sum())
+        #x1norm = numpy.sqrt((x1*x1).sum())
+        x1norm = x1.norm()
         #print ("x0 **********" ,x0)
         #print ("x1 **********" ,x1)
-        s[it] = (x1*x0).sum() / (x0*x0).sum()
+        s[it] = (x1*x0).sum() / (x0.squared_norm())
         x0 = (1.0/x1norm)*x1
     return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
 
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 57afd57..0000000
--- a/Wrappers/Python/conda-recipe/run_test.py
+++ /dev/null
@@ -1,1145 +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

-

-

-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)

-        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])

-

-

-

-

-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 test_FBPD_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)

-

-            # 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 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)

-

-            # 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, 1)

-            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]))

-

-class TestCompositeDataContainer(unittest.TestCase):

-

-    def test_one(self):

-        from ccpi.optimisation.ops import PowerMethodNonsquare

-        from ccpi.optimisation.ops import TomoIdentity

-        from ccpi.optimisation.funcs import Norm2sq, Norm1

-        from ccpi.framework import ImageGeometry, AcquisitionGeometry

-        from ccpi.optimisation.operators.CompositeOperator import CompositeDataContainer, CompositeDataContainer

-        import matplotlib.pyplot as plt

-    

-        ig0 = ImageGeometry(2,3,4)

-        ig1 = ImageGeometry(12,42,55,32)

-        

-        data0 = ImageData(geometry=ig0)

-        data1 = ImageData(geometry=ig1) + 1

-        

-        data2 = ImageData(geometry=ig0) + 2

-        data3 = ImageData(geometry=ig1) + 3

-        

-        cp0 = CompositeDataContainer(data0,data1)

-        cp1 = CompositeDataContainer(data2,data3)

-    #    

-        a = [ (el, ot) for el,ot in zip(cp0.containers,cp1.containers)]

-        print  (a[0][0].shape)

-        #cp2 = CompositeDataContainer(*a)

-        cp2 = cp0.add(cp1)

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 2.)

-        assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.)

-        

-        cp2 = cp0 + cp1 

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 2.)

-        assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.)

-        cp2 = cp0 + 1 

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)

-        cp2 = cp0 + [1 ,2]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 3., decimal = 5)

-        cp2 += cp1

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , +3. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)

-        

-        cp2 += 1

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , +4. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +7., decimal = 5)

-        

-        cp2 += [-2,-1]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 2. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 6., decimal = 5)

-        

-        

-        cp2 = cp0.subtract(cp1)

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == -2.)

-        assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.)

-        cp2 = cp0 - cp1

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == -2.)

-        assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.)

-        

-        cp2 = cp0 - 1 

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -1. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0, decimal = 5)

-        cp2 = cp0 - [1 ,2]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -1. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -1., decimal = 5)

-        

-        cp2 -= cp1

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -3. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5)

-        

-        cp2 -= 1

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -4. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -5., decimal = 5)

-        

-        cp2 -= [-2,-1]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -2. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5)

-        

-        

-        cp2 = cp0.multiply(cp1)

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)

-        assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.)

-        cp2 = cp0 * cp1

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)

-        assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.)

-        

-        cp2 = cp0 * 2 

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2, decimal = 5)

-        cp2 = cp0 * [3 ,2]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5)

-        

-        cp2 *= cp1

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0 , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)

-        

-        cp2 *= 1

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5)

-        

-        cp2 *= [-2,-1]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -6., decimal = 5)

-        

-        

-        cp2 = cp0.divide(cp1)

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4)

-        cp2 = cp0/cp1

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4)

-        

-        cp2 = cp0 / 2 

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)

-        cp2 = cp0 / [3 ,2]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)

-        

-        cp2 += 1

-        cp2 /= cp1

-        # TODO fix inplace division

-         

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1./2 , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1.5/3., decimal = 5)

-        

-        cp2 /= 1

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0.5 , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5)

-        

-        cp2 /= [-2,-1]

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -0.5/2. , decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -0.5, decimal = 5)

-        ####

-        

-        cp2 = cp0.power(cp1)

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)

-        cp2 = cp0**cp1

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)

-        

-        cp2 = cp0 ** 2 

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0., decimal=5)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1., decimal = 5)

-        

-        cp2 = cp0.maximum(cp1)

-        assert (cp2.get_item(0,0).as_array()[0][0][0] == cp1.get_item(0,0).as_array()[0][0][0])

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], cp2.get_item(1,0).as_array()[0][0][0], decimal=4)

-        

-        

-        cp2 = cp0.abs()

-        numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0], 0., decimal=4)

-        numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4)

-        

-        cp2 = cp0.subtract(cp1)

-        s = cp2.sign()

-        numpy.testing.assert_almost_equal(s.get_item(0,0).as_array()[0][0][0], -1., decimal=4)

-        numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], -1., decimal=4)

-        

-        cp2 = cp0.add(cp1)

-        s = cp2.sqrt()

-        numpy.testing.assert_almost_equal(s.get_item(0,0).as_array()[0][0][0], numpy.sqrt(2), decimal=4)

-        numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], numpy.sqrt(4), decimal=4)

-        

-        s = cp0.sum()

-        numpy.testing.assert_almost_equal(s[0], 0, decimal=4)

-        s0 = 1

-        s1 = 1

-        for i in cp0.get_item(0,0).shape:

-            s0 *= i

-        for i in cp0.get_item(1,0).shape:

-            s1 *= i

-            

-        numpy.testing.assert_almost_equal(s[1], cp0.get_item(0,0).as_array()[0][0][0]*s0 +cp0.get_item(1,0).as_array()[0][0][0]*s1, decimal=4)

-        

-

-if __name__ == '__main__':

-    unittest.main()

-    

diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py
index 00eeed0..8ebb857 100644
--- a/Wrappers/Python/setup.py
+++ b/Wrappers/Python/setup.py
@@ -23,12 +23,20 @@ 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,
+<<<<<<< HEAD
     packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation', 'ccpi.optimisation.operators'],
+=======
+    packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation', 
+              'ccpi.optimisation.algorithms'],
+>>>>>>> master
 
     # 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()
+    
diff --git a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py
index 0ec116f..8370c78 100644
--- a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py
+++ b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py
@@ -32,8 +32,8 @@ ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32')
 
 #########################################################################
 print ("Loading the data...")
-MainPath = '/media/algol/336F96987817D4B4/DATA/IMAT_DATA/' # path to data
-pathname0 = '{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/Sample/')
+MainPath = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/' # path to data
+pathname0 = '{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/Sample/')
 counterFileName = 4675
 # A main loop over all available angles 
 for ll in range(0,totalAngles,1):
@@ -66,7 +66,7 @@ for ll in range(0,totalAngles,1):
     counterFileName += 1
 #########################################################################
 
-flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits'))
+flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits'))
 nonzero = flat1 > 0
 # Apply flat field and take negative log
 for ll in range(0,totalAngles,1):
diff --git a/Wrappers/Python/wip/demo_imat_whitebeam.py b/Wrappers/Python/wip/demo_imat_whitebeam.py
index 482c1ae..e0d213e 100644
--- a/Wrappers/Python/wip/demo_imat_whitebeam.py
+++ b/Wrappers/Python/wip/demo_imat_whitebeam.py
@@ -21,18 +21,18 @@ from ccpi.optimisation.algs import CGLS, FISTA
 from ccpi.optimisation.funcs import Norm2sq, Norm1
 
 # Load and display a couple of summed projection as examples
-pathname0 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle0/'
+pathname0 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/'
 filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits'
 
 data0 = read_fits(pathname0 + filename0)
 
-pathname10 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle10/'
+pathname10 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/'
 filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits'
 
 data10 = read_fits(pathname10 + filename10)
 
 # Load a flat field (more are available, should we average over them?)
-flat1 = read_fits('/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')
+flat1 = read_fits('/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')
 
 # Apply flat field and display after flat-field correction and negative log
 data0_rel = numpy.zeros(numpy.shape(flat1), dtype = float)
@@ -58,7 +58,7 @@ plt.colorbar()
 plt.show()
 
 # Set up for loading all summed images at 250 angles.
-pathname = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle{}/'
+pathname = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/'
 filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits'
 
 # Dimensions