From d6b245ab0e0a7acc8380e3da9b2f262581da352e Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Thu, 28 Sep 2023 13:35:41 +0100
Subject: [PATCH 1/9] add symmetric flag to PowerMethod (#1470)

* add method flag to PowerMethod

* force LinearOperator::calculate_norm to use 'compose_with_adjoint' method

* Unit tests and case of nil-potent matrix

* Debugging and tests

* Changed warnings to logging.warning

closes #1468

---------

Signed-off-by: Edoardo Pasca <edo.paskino@gmail.com>
Signed-off-by: Margaret Duff <43645617+MargaretDuff@users.noreply.github.com>
Signed-off-by: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
Co-authored-by: Margaret Duff <margaret.duff@stfc.ac.uk>
Co-authored-by: Margaret Duff <43645617+MargaretDuff@users.noreply.github.com>
Co-authored-by: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
---
 CHANGELOG.md                                  |   2 +-
 .../cil/optimisation/operators/Operator.py    | 322 ++++++++++--------
 Wrappers/Python/test/test_BlockOperator.py    |   1 -
 Wrappers/Python/test/test_Operator.py         |  53 ++-
 4 files changed, 233 insertions(+), 145 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3b8ce3c690..06d756794a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,9 +9,9 @@
   - Fix for show_geometry bug for 2D data
   - Added warmstart capability to proximal evaluation of the CIL TotalVariation function.
   - FBP split processing bug fix - now respects panel origin
+  - Bug fix in the LinearOperator norm with an additional flag for the algorithm linearOperator.PowerMethod
   - Tidied up documentation in the framework folder
 
-
 * 23.0.1
   - Fix bug with NikonReader requiring ROI to be set in constructor.
 
diff --git a/Wrappers/Python/cil/optimisation/operators/Operator.py b/Wrappers/Python/cil/optimisation/operators/Operator.py
index cc2eb44bb4..adf89b9fba 100644
--- a/Wrappers/Python/cil/optimisation/operators/Operator.py
+++ b/Wrappers/Python/cil/optimisation/operators/Operator.py
@@ -20,7 +20,8 @@
 from numbers import Number
 import numpy
 import functools
-import warnings
+import logging
+
 
 class Operator(object):
     """
@@ -45,11 +46,11 @@ def __init__(self, domain_geometry, **kwargs):
     def is_linear(self):
         '''Returns if the operator is linear'''
         return False
-    def direct(self,x, out=None):
+
+    def direct(self, x, out=None):
         '''Returns the application of the Operator on x'''
         raise NotImplementedError
 
-
     def norm(self, **kwargs):
         '''Returns the norm of the Operator. On first call the norm will be calculated using the operator's calculate_norm
         method. Subsequent calls will return the cached norm.
@@ -60,7 +61,7 @@ def norm(self, **kwargs):
         '''
 
         if len(kwargs) != 0:
-            warnings.warn('norm: the norm method does not use any parameters.\n\
+            logging.warning('norm: the norm method does not use any parameters.\n\
                 For LinearOperators you can use PowerMethod to calculate the norm with non-default parameters and use set_norm to set it')
 
         if self._norm is None:
@@ -68,7 +69,7 @@ def norm(self, **kwargs):
 
         return self._norm
 
-    def set_norm(self,norm=None):
+    def set_norm(self, norm=None):
         '''Sets the norm of the operator to a custom value.
         '''
         self._norm = norm
@@ -76,44 +77,49 @@ def set_norm(self,norm=None):
     def calculate_norm(self):
         '''Calculates the norm of the Operator'''
         raise NotImplementedError
+
     def range_geometry(self):
         '''Returns the range of the Operator: Y space'''
         return self._range_geometry
+
     def domain_geometry(self):
         '''Returns the domain of the Operator: X space'''
         return self._domain_geometry
+
     @property
     def domain(self):
         return self.domain_geometry()
+
     @property
     def range(self):
         return self.range_geometry()
+
     def __rmul__(self, scalar):
         '''Defines the multiplication by a scalar on the left
 
         returns a ScaledOperator'''
         return ScaledOperator(self, scalar)
-    
+
     def compose(self, *other, **kwargs):
-        # TODO: check equality of domain and range of operators        
-        #if self.operator2.range_geometry != self.operator1.domain_geometry:
-        #    raise ValueError('Cannot compose operators, check domain geometry of {} and range geometry of {}'.format(self.operato1,self.operator2))    
-        
-        return CompositionOperator(self, *other, **kwargs) 
+        # TODO: check equality of domain and range of operators
+        # if self.operator2.range_geometry != self.operator1.domain_geometry:
+        #    raise ValueError('Cannot compose operators, check domain geometry of {} and range geometry of {}'.format(self.operato1,self.operator2))
+
+        return CompositionOperator(self, *other, **kwargs)
 
     def __add__(self, other):
         return SumOperator(self, other)
 
     def __mul__(self, scalar):
-        return self.__rmul__(scalar)    
-    
+        return self.__rmul__(scalar)
+
     def __neg__(self):
         """ Return -self """
-        return -1 * self    
-        
+        return -1 * self
+
     def __sub__(self, other):
         """ Returns the subtraction of the operators."""
-        return self + (-1) * other   
+        return self + (-1) * other
 
 
 class LinearOperator(Operator):
@@ -129,25 +135,26 @@ class LinearOperator(Operator):
     range_geometry : ImageGeometry or AcquisitionGeometry, optional, default None
         range of the operator 
     """
+
     def __init__(self, domain_geometry, **kwargs):
         super(LinearOperator, self).__init__(domain_geometry, **kwargs)
+
     def is_linear(self):
         '''Returns if the operator is linear'''
         return True
-    def adjoint(self,x, out=None):
+
+    def adjoint(self, x, out=None):
         '''returns the adjoint/inverse operation
-        
+
         only available to linear operators'''
         raise NotImplementedError
-    
-    @staticmethod
-    def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5,  return_all=False):
 
+    @staticmethod
+    def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5,  return_all=False, method='auto'):
         r"""Power method or Power iteration algorithm 
-        
-        The Power method computes the largest (dominant) eigenvalue of a square matrix in magnitude, e.g.,
-        absolute value in the real case and module in the complex case.        
-        For the non-square case, the power method is applied for the matrix :math: A^{T}*A.
+
+        The Power method computes the largest (dominant) eigenvalue of a matrix in magnitude, e.g.,
+        absolute value in the real case and modulus in the complex case.        
 
         Parameters
         ----------
@@ -159,9 +166,11 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5,  ret
             Starting point for the Power method.
         tolerance: positive:`float`, default = 1e-5
             Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance.                    
-        return_all: boolean, default = False
+        return_all: `boolean`, default = False
             Toggles the verbosity of the return
-    
+        method: `string` one of `"auto"`, `"composed_with_adjoint"` and `"direct_only"`, default = `"auto"` 
+            The default `auto` lets the code choose the method, this can be specified with `"direct_only"` or `"composed_with_adjoint"`
+
 
         Returns
         -------
@@ -172,6 +181,23 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5,  ret
             Corresponding eigenvector of the dominant eigenvalue. Only returned if return_all is True.
         list of eigenvalues: :obj:`list`
             List of eigenvalues. Only returned if return_all is True.
+        convergence: `boolean`
+            Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True.
+
+
+        Note
+        -----
+        The power method contains two different algorithms chosen by the `method` flag. 
+
+        In the case `method="direct_only"`, for operator, :math:`A`, the power method computes the iterations 
+        :math:`x_{k+1} = A (x_k/\|x_{k}\|)` initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. 
+
+        In the case `method="composed_with_adjoint"`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` 
+        returning the square root of this value, i.e. the iterations:
+        :math:`x_{k+1} = A^TA (x_k/\|x_{k}\|)` and returning  :math:`\sqrt{\|x_k\|}`.
+
+        The default flag is `method="auto"`, the algorithm checks to see if the `operator.domain_geometry() == operator.range_geometry()` and if so
+        uses the method "direct_only" and if not the method "composed_with_adjoint".
 
         Examples
         --------    
@@ -189,14 +215,25 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5,  ret
 
         """
 
-        # Default case: non-symmetric
-        symmetric = False
-        try:
-            if operator.domain_geometry()==operator.range_geometry():
-                symmetric = True
-        except AssertionError:
-            # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972
-            pass
+        allowed_methods = ["auto", "direct_only", "composed_with_adjoint"]
+
+        if method not in allowed_methods:
+            raise ValueError("The argument 'method' can be set to one of {0} got {1}".format(
+                allowed_methods, method))
+
+        apply_adjoint = True
+        if method == "direct_only":
+            apply_adjoint = False
+        if method == "auto":
+            try:
+                geometries_match = operator.domain_geometry() == operator.range_geometry()
+
+            except AssertionError:
+                # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972
+                pass
+            else:
+                if geometries_match:
+                    apply_adjoint = False
 
         if initial is None:
             x0 = operator.domain_geometry().allocate('random')
@@ -211,58 +248,59 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance = 1e-5,  ret
 
         # initial guess for dominant eigenvalue
         eig_old = 1.
-
-        eig_list = []
+        if return_all:
+            eig_list = []
+            convergence_check = True
         diff = numpy.finfo('d').max
         i = 0
         while (i < max_iteration and diff > tolerance):
-            i+=1
-            
-            operator.direct(x0, out = y_tmp)
+            operator.direct(x0, out=y_tmp)
 
-            if symmetric:                
-                #swap datacontainer references
+            if not apply_adjoint:
+                # swap datacontainer references
                 tmp = x0
                 x0 = y_tmp
                 y_tmp = tmp
             else:
-                operator.adjoint(y_tmp,out=x0)
-            
+                operator.adjoint(y_tmp, out=x0)
+
             # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization
             x0_norm = x0.norm()
+            if x0_norm < tolerance:
+                logging.warning(
+                    'The operator has at least one zero eigenvector and is likely to be nilpotent')
+                eig_new = 0.
+                break
             x0 /= x0_norm
 
-            eig_new =  numpy.abs(x0_norm)
-            if not symmetric:
+            eig_new = numpy.abs(x0_norm)
+            if apply_adjoint:
                 eig_new = numpy.sqrt(eig_new)
-
             diff = numpy.abs(eig_new - eig_old)
-            eig_list.append(eig_new)   
-            eig_old = eig_new       
+            if return_all:
+                eig_list.append(eig_new)
+            eig_old = eig_new
+            i += 1
+
+        if return_all and i == max_iteration:
+            convergence_check = False
 
         if return_all:
-            return eig_new, i, x0, eig_list
+            return eig_new, i, x0, eig_list, convergence_check
         else:
             return eig_new
-            
 
     def calculate_norm(self):
-        
         r""" Returns the norm of the LinearOperator calculated by the PowerMethod with default values.
                 """
-        return LinearOperator.PowerMethod(self)
-
+        return LinearOperator.PowerMethod(self, method="composed_with_adjoint")
 
     @staticmethod
     def dot_test(operator, domain_init=None, range_init=None, tolerance=1e-6, **kwargs):
         r'''Does a dot linearity test on the operator
-        
         Evaluates if the following equivalence holds
-        
         .. math::
-        
           Ax\times y = y \times A^Tx
-        
         :param operator: operator to test the dot_test
         :param range_init: optional initialisation container in the operator range 
         :param domain_init: optional initialisation container in the operator domain 
@@ -270,53 +308,52 @@ def dot_test(operator, domain_init=None, range_init=None, tolerance=1e-6, **kwar
         :type : int, default = 1
         :param tolerance: Check if the following expression is below the tolerance
         .. math:: 
-        
+
             |Ax\times y - y \times A^Tx|/(\|A\|\|x\|\|y\| + 1e-12) < tolerance
-        
+
         :type : float, default 1e-6
         :returns: boolean, True if the test is passed.       
         '''
 
         seed = kwargs.get('seed', 1)
-    
+
         if range_init is None:
-            y = operator.range_geometry().allocate('random', seed = seed + 10)
+            y = operator.range_geometry().allocate('random', seed=seed + 10)
         else:
             y = range_init
         if domain_init is None:
-            x = operator.domain_geometry().allocate('random', seed = seed)
+            x = operator.domain_geometry().allocate('random', seed=seed)
         else:
             x = domain_init
-        
+
         fx = operator.direct(x)
         by = operator.adjoint(y)
         a = fx.dot(y)
         b = by.dot(x).conjugate()
 
-        # Check relative tolerance but normalised with respect to 
+        # Check relative tolerance but normalised with respect to
         # operator, x and y norms and avoid zero division
-        error = numpy.abs( a - b )/ (operator.norm()*x.norm()*y.norm() + 1e-12)
-            
+        error = numpy.abs(a - b) / (operator.norm()*x.norm()*y.norm() + 1e-12)
+
         if error < tolerance:
             return True
         else:
-            print ('Left hand side  {}, \nRight hand side {}'.format(a, b))
-            return False    
-        
-        
+            print('Left hand side  {}, \nRight hand side {}'.format(a, b))
+            return False
+
+
 class ScaledOperator(Operator):
-    
-    
+
     '''ScaledOperator
 
     A class to represent the scalar multiplication of an Operator with a scalar.
     It holds an operator and a scalar. Basically it returns the multiplication
     of the result of direct and adjoint of the operator with the scalar.
     For the rest it behaves like the operator it holds.
-    
+
     :param operator: a Operator or LinearOperator
     :param scalar: a scalar multiplier
-    
+
     Example:
        The scaled operator behaves like the following:
 
@@ -330,7 +367,7 @@ class ScaledOperator(Operator):
       sop.domain_geometry() = operator.domain_geometry()
 
     '''
-    
+
     def __init__(self, operator, scalar, **kwargs):
         '''creator
 
@@ -338,12 +375,13 @@ def __init__(self, operator, scalar, **kwargs):
         :param scalar: a scalar multiplier
         :type scalar: Number'''
 
-        super(ScaledOperator, self).__init__(domain_geometry=operator.domain_geometry(), 
+        super(ScaledOperator, self).__init__(domain_geometry=operator.domain_geometry(),
                                              range_geometry=operator.range_geometry())
-        if not isinstance (scalar, Number):
+        if not isinstance(scalar, Number):
             raise TypeError('expected scalar: got {}'.format(type(scalar)))
         self.scalar = scalar
         self.operator = operator
+
     def direct(self, x, out=None):
         '''direct method'''
         if out is None:
@@ -353,6 +391,7 @@ def direct(self, x, out=None):
         else:
             self.operator.direct(x, out=out)
             out *= self.scalar
+
     def adjoint(self, x, out=None):
         '''adjoint method'''
         if self.operator.is_linear():
@@ -365,62 +404,62 @@ def adjoint(self, x, out=None):
                 out *= self.scalar
         else:
             raise TypeError('Operator is not linear')
+
     def norm(self, **kwargs):
         '''norm of the operator'''
         return numpy.abs(self.scalar) * self.operator.norm(**kwargs)
 
     def is_linear(self):
         '''returns whether the operator is linear
-        
+
         :returns: boolean '''
         return self.operator.is_linear()
 
 
 ###############################################################################
 ################   SumOperator  ###########################################
-###############################################################################      
-    
+###############################################################################
+
 class SumOperator(Operator):
-    
-    
+
     def __init__(self, operator1, operator2):
-                
+
         self.operator1 = operator1
         self.operator2 = operator2
-        
+
         # if self.operator1.domain_geometry() != self.operator2.domain_geometry():
-        #     raise ValueError('Domain geometry of {} is not equal with domain geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))    
-                
+        #     raise ValueError('Domain geometry of {} is not equal with domain geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))
+
         # if self.operator1.range_geometry() != self.operator2.range_geometry():
-        #     raise ValueError('Range geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))    
-            
-        self.linear_flag = self.operator1.is_linear() and self.operator2.is_linear()            
-        
+        #     raise ValueError('Range geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))
+
+        self.linear_flag = self.operator1.is_linear() and self.operator2.is_linear()
+
         super(SumOperator, self).__init__(domain_geometry=self.operator1.domain_geometry(),
-                                          range_geometry=self.operator1.range_geometry()) 
-                                  
+                                          range_geometry=self.operator1.range_geometry())
+
     def direct(self, x, out=None):
-        
+
         if out is None:
             return self.operator1.direct(x) + self.operator2.direct(x)
         else:
             self.operator1.direct(x, out=out)
-            out.add(self.operator2.direct(x), out=out)     
+            out.add(self.operator2.direct(x), out=out)
 
     def adjoint(self, x, out=None):
-        
-        if self.linear_flag:        
+
+        if self.linear_flag:
             if out is None:
                 return self.operator1.adjoint(x) + self.operator2.adjoint(x)
             else:
                 self.operator1.adjoint(x, out=out)
-                out.add(self.operator2.adjoint(x), out=out) 
+                out.add(self.operator2.adjoint(x), out=out)
         else:
             raise ValueError('No adjoint operation with non-linear operators')
-                                        
+
     def is_linear(self):
-        return self.linear_flag 
-    
+        return self.linear_flag
+
     def calculate_norm(self):
         if self.is_linear():
             return LinearOperator.calculate_norm(self)
@@ -429,43 +468,46 @@ def calculate_norm(self):
 
 ###############################################################################
 ################   Composition  ###########################################
-###############################################################################             
+###############################################################################
+
 
 class CompositionOperator(Operator):
-    
+
     def __init__(self, *operators, **kwargs):
-        
+
         # get a reference to the operators
         self.operators = operators
-        
-        self.linear_flag = functools.reduce(lambda x,y: x and y.is_linear(),
+
+        self.linear_flag = functools.reduce(lambda x, y: x and y.is_linear(),
                                             self.operators, True)
         # self.preallocate = kwargs.get('preallocate', False)
         self.preallocate = False
         if self.preallocate:
-            self.tmp_domain = [op.domain_geometry().allocate() for op in self.operators[:-1]]
-            self.tmp_range = [op.range_geometry().allocate() for op in self.operators[1:]]
+            self.tmp_domain = [op.domain_geometry().allocate()
+                               for op in self.operators[:-1]]
+            self.tmp_range = [op.range_geometry().allocate()
+                              for op in self.operators[1:]]
             # pass
-        
+
         # TODO address the equality of geometries
         # if self.operator2.range_geometry() != self.operator1.domain_geometry():
-        #     raise ValueError('Domain geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))    
-                
+        #     raise ValueError('Domain geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))
+
         super(CompositionOperator, self).__init__(
             domain_geometry=self.operators[-1].domain_geometry(),
-            range_geometry=self.operators[0].range_geometry()) 
-        
-    def direct(self, x, out = None):
+            range_geometry=self.operators[0].range_geometry())
+
+    def direct(self, x, out=None):
 
         if out is None:
-            #return self.operator1.direct(self.operator2.direct(x))
-            # return functools.reduce(lambda X,operator: operator.direct(X), 
+            # return self.operator1.direct(self.operator2.direct(x))
+            # return functools.reduce(lambda X,operator: operator.direct(X),
             #                        self.operators[::-1][1:],
             #                        self.operators[-1].direct(x))
             if self.preallocate:
                 pass
             else:
-                for i,operator in enumerate(self.operators[::-1]):
+                for i, operator in enumerate(self.operators[::-1]):
                     if i == 0:
                         step = operator.direct(x)
                     else:
@@ -476,81 +518,81 @@ def direct(self, x, out = None):
             # tmp = self.operator2.range_geometry().allocate()
             # self.operator2.direct(x, out = tmp)
             # self.operator1.direct(tmp, out = out)
-            
+
             # out.fill (
-            #     functools.reduce(lambda X,operator: operator.direct(X), 
+            #     functools.reduce(lambda X,operator: operator.direct(X),
             #                        self.operators[::-1][1:],
             #                        self.operators[-1].direct(x))
             # )
-            
+
             # TODO this is a bit silly but will handle the pre allocation later
             if self.preallocate:
 
-                for i,operator in enumerate(self.operators[::-1]):
+                for i, operator in enumerate(self.operators[::-1]):
                     if i == 0:
                         operator.direct(x, out=self.tmp_range[i])
                     elif i == len(self.operators) - 1:
                         operator.direct(self.tmp_range[i-1], out=out)
                     else:
-                        operator.direct(self.tmp_range[i-1], out=self.tmp_range[i])
+                        operator.direct(
+                            self.tmp_range[i-1], out=self.tmp_range[i])
             else:
-                for i,operator in enumerate(self.operators[::-1]):
+                for i, operator in enumerate(self.operators[::-1]):
                     if i == 0:
                         step = operator.direct(x)
                     else:
                         step = operator.direct(step)
                 out.fill(step)
 
-            
-    def adjoint(self, x, out = None):
-        
-        if self.linear_flag: 
-            
+    def adjoint(self, x, out=None):
+
+        if self.linear_flag:
+
             if out is not None:
                 # return self.operator2.adjoint(self.operator1.adjoint(x))
-                # return functools.reduce(lambda X,operator: operator.adjoint(X), 
+                # return functools.reduce(lambda X,operator: operator.adjoint(X),
                 #                    self.operators[1:],
                 #                    self.operators[0].adjoint(x))
                 if self.preallocate:
-                    for i,operator in enumerate(self.operators):
+                    for i, operator in enumerate(self.operators):
                         if i == 0:
                             operator.adjoint(x, out=self.tmp_domain[i])
                         elif i == len(self.operators) - 1:
-                            step = operator.adjoint(self.tmp_domain[i-1], out=out)
+                            step = operator.adjoint(
+                                self.tmp_domain[i-1], out=out)
                         else:
-                            operator.adjoint(self.tmp_domain[i-1], out=self.tmp_domain[i])
+                            operator.adjoint(
+                                self.tmp_domain[i-1], out=self.tmp_domain[i])
                     return
                 else:
-                    for i,operator in enumerate(self.operators):
+                    for i, operator in enumerate(self.operators):
                         if i == 0:
                             step = operator.adjoint(x)
                         else:
                             step = operator.adjoint(step)
                     out.fill(step)
-                
+
             else:
                 if self.preallocate:
                     pass
                 else:
-                    for i,operator in enumerate(self.operators):
+                    for i, operator in enumerate(self.operators):
                         if i == 0:
                             step = operator.adjoint(x)
                         else:
                             step = operator.adjoint(step)
-                    
+
                     return step
         else:
             raise ValueError('No adjoint operation with non-linear operators')
-            
 
     def is_linear(self):
-        return self.linear_flag             
-            
+        return self.linear_flag
 
     def calculate_norm(self):
         '''Returns the norm of the CompositionOperator, that is the product of the norms
         of its operators.'''
         norm = 1.
         for operator in self.operators:
-                norm *= operator.norm()
+            norm *= operator.norm()
         return norm
diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py
index 3e81ab4cad..93e5bbd3e2 100644
--- a/Wrappers/Python/test/test_BlockOperator.py
+++ b/Wrappers/Python/test/test_BlockOperator.py
@@ -187,7 +187,6 @@ def test_FiniteDiffOperator(self):
         G = FiniteDifferenceOperator(ig, direction=0, bnd_cond = 'Neumann')
         logging.info("{} {}".format(type(u), str(u.as_array())))    
         logging.info(str(G.direct(u).as_array()))
-
         # Gradient Operator norm, for one direction should be close to 2
         numpy.testing.assert_allclose(G.norm(), numpy.sqrt(4), atol=0.1)
 
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 78d7e30c8c..6a96610fca 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -35,6 +35,7 @@
 import logging
 from testclass import CCPiTestClass
 
+
 initialise_tests()
 
 def dt(steps):
@@ -284,6 +285,7 @@ def test_PowerMethod(self):
 
         # Test with the norm       
         res2 = M1op.norm()
+        res1 = M1op.PowerMethod(M1op,100, method="composed_with_adjoint")
         numpy.testing.assert_almost_equal(res1,res2, decimal=4)
 
 
@@ -302,8 +304,36 @@ def test_PowerMethod(self):
         # 3x3 complex matrix, (real+complex eigenvalue), dominant eigenvalue = 3.1624439599276974
         M1 = numpy.array([[2,0,0],[1,2j,1j],[3, 3-1j,3]])
         M1op = MatrixOperator(M1)
-        res1 = M1op.PowerMethod(M1op,100)
-        numpy.testing.assert_almost_equal(res1,3.1624439599276974, decimal=4)                
+        res1 = M1op.PowerMethod(M1op,120)
+        numpy.testing.assert_almost_equal(res1,3.1624439599276974, decimal=3) 
+
+        # 2x2 non-diagonalisable nilpotent matrix
+        M1=numpy.array([[0.,1.], [0.,0.]])
+        M1op = MatrixOperator(M1)
+        res1 = M1op.PowerMethod(M1op,5)
+        numpy.testing.assert_almost_equal(res1,0, decimal=4) 
+
+        # 2x2 non-diagonalisable nilpotent matrix where method="composed_with_adjoint"
+        M1=numpy.array([[0.,1.], [0.,0.]])
+        M1op = MatrixOperator(M1)
+        res1 = M1op.PowerMethod(M1op,5, method="composed_with_adjoint")
+        numpy.testing.assert_almost_equal(res1,1, decimal=4) 
+
+
+        # 2x2 matrix, max absolute eigenvalue is not unique and initial vector chosen for non-convergence
+        
+        M1=numpy.array([[2.,1.], [0.,-2.]])
+        M1op = MatrixOperator(M1)
+        _,_,_,_,convergence = M1op.PowerMethod(M1op,100, initial=DataContainer(numpy.array([1.,1.])), return_all=True)
+        numpy.testing.assert_equal(convergence,False) 
+
+        # 2x2 matrix, max absolute eigenvalue is not unique and initial vector chosen for convergence
+        
+        M1=numpy.array([[2.,1.,0.],[0.,1.,1.], [0.,0.,1.]])
+        M1op = MatrixOperator(M1)
+        res1,_,_,_,convergence = M1op.PowerMethod(M1op,100, return_all=True)
+        numpy.testing.assert_almost_equal(res1,2., decimal=4) 
+        numpy.testing.assert_equal(convergence,True)     
 
         # Gradient Operator (float)
         ig = ImageGeometry(30,30)
@@ -320,7 +350,13 @@ def test_PowerMethod(self):
         # Identity Operator
         Id = IdentityOperator(ig)
         res1 = Id.PowerMethod(Id,100)
-        numpy.testing.assert_almost_equal(res1,1.0, decimal=4)                
+        numpy.testing.assert_almost_equal(res1,1.0, decimal=4)        
+
+        # Test errors produced if not a valid method
+        try:
+            res1 = Id.PowerMethod(Id,100, method='gobledigook')
+        except ValueError:
+            pass      
 
 
     def test_Norm(self):
@@ -347,6 +383,17 @@ def test_Norm(self):
         #recalculates norm
         self.assertAlmostEqual(G.norm(), numpy.sqrt(8), 2)
 
+         # 2x2 real matrix, dominant eigenvalue = 2. Check norm uses the right flag for power method 
+        M1 = numpy.array([[1,0],[1,2]], dtype=float)
+        M1op = MatrixOperator(M1)
+        res1 = M1op.norm()
+        res2 = M1op.PowerMethod(M1op,100)
+        res3 = M1op.PowerMethod(M1op,100, method="composed_with_adjoint")
+        res4 = M1op.PowerMethod(M1op,100, method="direct_only")
+        numpy.testing.assert_almost_equal(res1,res3, decimal=4)
+        self.assertNotEqual(res1, res2)
+        self.assertNotEqual(res1,res4)
+
 
     def test_ProjectionMap(self):
         # Check if direct is correct

From 9d6b4d5bfd0d3adf558a23125d93fdb1599aa728 Mon Sep 17 00:00:00 2001
From: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
Date: Thu, 28 Sep 2023 13:41:13 +0100
Subject: [PATCH 2/9] Padder, Slicer, Binner to respect panel origin (#1461)

* handle geometry correctly for other data origins

* unit and functionality tests w. all panel origins

* move detector shift to framework geometry

---------

Signed-off-by: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
---
 CHANGELOG.md                                  |   3 +-
 Wrappers/Python/cil/framework/framework.py    |  44 +++-
 Wrappers/Python/cil/processors/Padder.py      |   9 +-
 Wrappers/Python/cil/processors/Slicer.py      |   6 +-
 .../Python/test/test_AcquisitionGeometry.py   |  78 ++++++
 Wrappers/Python/test/test_DataProcessor.py    | 226 +++++++++++++++++-
 6 files changed, 349 insertions(+), 17 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 06d756794a..2ea88e4ce1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,8 +7,9 @@
   - Optimisation in L2NormSquared
   - Added support for partitioner, when partitions have size 1
   - Fix for show_geometry bug for 2D data
+  - FBP split processing bug fix - now respects panel origin set in geometry
+  - Binner/Padder/Slicer bug fix - now respects panel origin set in geometry
   - Added warmstart capability to proximal evaluation of the CIL TotalVariation function.
-  - FBP split processing bug fix - now respects panel origin
   - Bug fix in the LinearOperator norm with an additional flag for the algorithm linearOperator.PowerMethod
   - Tidied up documentation in the framework folder
 
diff --git a/Wrappers/Python/cil/framework/framework.py b/Wrappers/Python/cil/framework/framework.py
index 02558e00bd..163e071286 100644
--- a/Wrappers/Python/cil/framework/framework.py
+++ b/Wrappers/Python/cil/framework/framework.py
@@ -521,7 +521,7 @@ def __init__ (self, dof):
     @staticmethod  
     def create_vector(val):
         try:
-            vec = numpy.asarray(val, dtype=numpy.float64).reshape(len(val))
+            vec = numpy.array(val, dtype=numpy.float64).reshape(len(val))
         except:
             raise ValueError("Can't convert to numpy array")
    
@@ -1984,6 +1984,48 @@ def configured(self):
             configured = False
         return configured
 
+    def shift_detector_in_plane(self,
+                                          pixel_offset,
+                                          direction='horizontal'):
+        """
+        Adjusts the position of the detector in a specified direction within the imaging plane.
+
+        Parameters:
+        -----------
+        pixel_offset : float
+            The number of pixels to adjust the detector's position by.
+        direction : {'horizontal', 'vertical'}, optional
+            The direction in which to adjust the detector's position. Defaults to 'horizontal'.
+
+        Notes:
+        ------
+        - If `direction` is 'horizontal':
+            - If the panel's origin is 'left', positive offsets translate the detector to the right.
+            - If the panel's origin is 'right', positive offsets translate the detector to the left.
+
+        - If `direction` is 'vertical':
+            - If the panel's origin is 'bottom', positive offsets translate the detector upward.
+            - If the panel's origin is 'top', positive offsets translate the detector downward.
+
+        Returns:
+        --------
+        None
+        """
+
+        if direction == 'horizontal':
+            pixel_size = self.panel.pixel_size[0]
+            pixel_direction = self.system.detector.direction_x
+
+        elif direction == 'vertical':
+            pixel_size = self.panel.pixel_size[1]
+            pixel_direction = self.system.detector.direction_y
+
+        if 'bottom' in self.panel.origin or 'left' in self.panel.origin:
+            self.system.detector.position -= pixel_offset * pixel_direction * pixel_size
+        else:
+            self.system.detector.position += pixel_offset * pixel_direction * pixel_size
+
+
     def __str__(self):
         repres = ""
         if self.configured:
diff --git a/Wrappers/Python/cil/processors/Padder.py b/Wrappers/Python/cil/processors/Padder.py
index 879771db83..73f71420b6 100644
--- a/Wrappers/Python/cil/processors/Padder.py
+++ b/Wrappers/Python/cil/processors/Padder.py
@@ -555,7 +555,6 @@ def _process_acquisition_geometry(self):
                 continue
 
             offset = (self._pad_width_param[i][0] -self._pad_width_param[i][1])*0.5
-            system_detector = geometry.config.system.detector
 
             if dim == 'channel':
                 geometry.set_channels(num_channels= geometry.config.channels.num_channels + \
@@ -563,7 +562,6 @@ def _process_acquisition_geometry(self):
             elif dim == 'angle':
                 # extrapolate pre-values from a[1]-a[0]
                 # extrapolate post-values from a[-1]-a[-2]
-
                 a = self._geometry.angles
                 end_values = (
                     a[0]-(a[1]-a[0] )* self._pad_width_param[i][0],
@@ -574,12 +572,13 @@ def _process_acquisition_geometry(self):
             elif dim == 'vertical':
                 geometry.config.panel.num_pixels[1] += self._pad_width_param[i][0] 
                 geometry.config.panel.num_pixels[1] += self._pad_width_param[i][1]
-                system_detector.position = system_detector.position - offset * system_detector.direction_y * geometry.config.panel.pixel_size[1]
+                geometry.config.shift_detector_in_plane(offset, dim)
+                    
             elif dim == 'horizontal':
                 geometry.config.panel.num_pixels[0] += self._pad_width_param[i][0]
                 geometry.config.panel.num_pixels[0] += self._pad_width_param[i][1]
-                system_detector.position = system_detector.position - offset * system_detector.direction_x * geometry.config.panel.pixel_size[0]
-        
+                geometry.config.shift_detector_in_plane(offset, dim)
+
         return geometry
 
     def _process_image_geometry(self):
diff --git a/Wrappers/Python/cil/processors/Slicer.py b/Wrappers/Python/cil/processors/Slicer.py
index 4f8584b8c1..8aceca3d17 100644
--- a/Wrappers/Python/cil/processors/Slicer.py
+++ b/Wrappers/Python/cil/processors/Slicer.py
@@ -269,7 +269,6 @@ def _process_acquisition_geometry(self):
         Creates the new acquisition geometry
         """
         geometry_new = self._geometry.copy()
-        system_detector = geometry_new.config.system.detector
 
         processed_dims = self._processed_dims.copy()
 
@@ -283,7 +282,7 @@ def _process_acquisition_geometry(self):
                 if n_elements > 1:
                     # difference in end indices, minus differences in start indices, divided by 2
                     pixel_offset = ((self._shape_in[vert_ind] -1 - self._pixel_indices[vert_ind][1]) - self._pixel_indices[vert_ind][0])*0.5
-                    system_detector.position = system_detector.position - pixel_offset * system_detector.direction_y * geometry_new.config.panel.pixel_size[1]
+                    geometry_new.config.shift_detector_in_plane(pixel_offset, 'vertical')
                     geometry_new.config.panel.num_pixels[1] = n_elements
                 else:
                     try:
@@ -314,7 +313,8 @@ def _process_acquisition_geometry(self):
 
             elif axis == 'horizontal':
                 pixel_offset = ((self._shape_in[i] -1 - self._pixel_indices[i][1]) - self._pixel_indices[i][0])*0.5
-                system_detector.position = system_detector.position - pixel_offset * system_detector.direction_x * geometry_new.config.panel.pixel_size[0]
+
+                geometry_new.config.shift_detector_in_plane(pixel_offset, axis)
                 geometry_new.config.panel.num_pixels[0] = n_elements
                 geometry_new.config.panel.pixel_size[0] *= roi.step
 
diff --git a/Wrappers/Python/test/test_AcquisitionGeometry.py b/Wrappers/Python/test/test_AcquisitionGeometry.py
index 8b48dca847..07b4f4e052 100644
--- a/Wrappers/Python/test/test_AcquisitionGeometry.py
+++ b/Wrappers/Python/test/test_AcquisitionGeometry.py
@@ -157,6 +157,84 @@ def test_create_Cone3D(self):
         np.testing.assert_allclose(AG.config.system.rotation_axis.position, rotation_axis_position, rtol=1E-6)
         np.testing.assert_allclose(AG.config.system.rotation_axis.direction, rotation_axis_direction, rtol=1E-6)
 
+
+    def test_shift_detector_origin_bottom_left(self):
+        initial_position = np.array([2.5, -1.3, 10.2])
+        pixel_size = np.array([0.5, 0.7])
+        detector_direction_x=np.array([1/math.sqrt(2),0,1/math.sqrt(2)])
+        detector_direction_y=np.array([-1/math.sqrt(2),0,1/math.sqrt(2)])
+
+        geometry = AcquisitionGeometry.create_Parallel3D(detector_position=initial_position, detector_direction_x=detector_direction_x, detector_direction_y=detector_direction_y)\
+                                      .set_panel([10, 5], [0.5, 0.7], origin='bottom-left')\
+                                      .set_angles([0])
+        # Test horizontal shift to the left
+        shift = -1.5
+        geometry.config.shift_detector_in_plane(shift, 'horizontal')
+        updated_position = geometry.config.system.detector.position
+        expected_position = initial_position - detector_direction_x * pixel_size[0] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+        # Test horizontal shift to the right
+        shift = 3.0
+        geometry.config.shift_detector_in_plane(shift, 'horizontal')
+        updated_position = geometry.config.system.detector.position
+        expected_position = expected_position - detector_direction_x * pixel_size[0] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+        # Test vertical shift down
+        shift = -1.5
+        geometry.config.shift_detector_in_plane(shift, 'vertical')
+        updated_position = geometry.config.system.detector.position
+        expected_position = expected_position - detector_direction_y * pixel_size[1] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+        # Test vertical shift up
+        shift = 3.0
+        geometry.config.shift_detector_in_plane(shift, 'vertical')
+        updated_position = geometry.config.system.detector.position
+        expected_position = expected_position - detector_direction_y * pixel_size[1] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+
+    def test_shift_detector_origin_top_right(self):
+        initial_position = np.array([2.5, -1.3, 10.2])
+        detector_direction_x=np.array([1/math.sqrt(2),0,1/math.sqrt(2)])
+        detector_direction_y=np.array([-1/math.sqrt(2),0,1/math.sqrt(2)])
+
+        pixel_size = np.array([0.5, 0.7])
+        geometry = AcquisitionGeometry.create_Parallel3D(detector_position=initial_position, detector_direction_x=detector_direction_x, detector_direction_y=detector_direction_y)\
+                                      .set_panel([10, 5], [0.5, 0.7], origin='top-right')\
+                                      .set_angles([0])
+
+        # Test horizontal shift to the right
+        shift = -1.5
+        geometry.config.shift_detector_in_plane(shift, 'horizontal')
+        updated_position = geometry.config.system.detector.position
+        expected_position = initial_position + detector_direction_x * pixel_size[0] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+        # Test horizontal shift to the left
+        shift = 3.0
+        geometry.config.shift_detector_in_plane(shift, 'horizontal')
+        updated_position = geometry.config.system.detector.position
+        expected_position = expected_position + detector_direction_x * pixel_size[0] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+        # Test vertical shift up
+        shift = -1.5
+        geometry.config.shift_detector_in_plane(shift, 'vertical')
+        updated_position = geometry.config.system.detector.position
+        expected_position = expected_position + detector_direction_y * pixel_size[1] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+        # Test vertical shift down
+        shift = 3.0
+        geometry.config.shift_detector_in_plane(shift, 'vertical')
+        updated_position = geometry.config.system.detector.position
+        expected_position = expected_position +  detector_direction_y * pixel_size[1] * shift
+        np.testing.assert_array_almost_equal(updated_position, expected_position)
+
+
     def test_SystemConfiguration(self):
         
         #SystemConfiguration error handeling
diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py
index eaae2d2bea..8462c7219f 100644
--- a/Wrappers/Python/test/test_DataProcessor.py
+++ b/Wrappers/Python/test/test_DataProcessor.py
@@ -237,7 +237,8 @@ def test_process_acquisition_geometry_parallel2D(self):
 
     def test_process_acquisition_geometry_parallel3D(self):
 
-        ag = AcquisitionGeometry.create_Parallel3D().set_angles(numpy.linspace(0,360,360,endpoint=False)).set_panel([128,64],[0.1,0.2]).set_channels(4)
+        angles_full = numpy.linspace(0,360,360,endpoint=False)
+        ag = AcquisitionGeometry.create_Parallel3D().set_angles(angles_full).set_panel([128,64],[0.1,0.2], origin='bottom-left').set_channels(4)
 
         rois = [
                 # same as input
@@ -248,13 +249,21 @@ def test_process_acquisition_geometry_parallel3D(self):
                 
                 # bin to single dimension
                 {'vertical':(31,33,2)},
+
+        
+                # crop  asymmetrically
+                {'vertical':(10,None,None),'horizontal':(None,-20,None)}                
         ]
 
+        #calculate offsets for geometry4
+        ag4_offset_h = -ag.pixel_size_h*20/2
+        ag4_offset_v = ag.pixel_size_v*10/2
 
         ag_gold = [
                 ag.copy(),
-                AcquisitionGeometry.create_Parallel3D().set_angles(numpy.linspace(0.5,360.5,180,endpoint=False)).set_panel([8,8],[1.6,1.6]).set_channels(1),
-                AcquisitionGeometry.create_Parallel2D().set_angles(numpy.linspace(0,360,360,endpoint=False)).set_panel(128,[0.1,0.4]).set_channels(4),        
+                AcquisitionGeometry.create_Parallel3D().set_angles(numpy.linspace(0.5,360.5,180,endpoint=False)).set_panel([8,8],[1.6,1.6], origin='bottom-left').set_channels(1),
+                AcquisitionGeometry.create_Parallel2D().set_angles(angles_full).set_panel(128,[0.1,0.4], origin='bottom-left').set_channels(4),
+                AcquisitionGeometry.create_Parallel3D(detector_position=[ag4_offset_h, 0, ag4_offset_v]).set_angles(angles_full).set_panel([108,54],[0.1,0.2], origin='bottom-left').set_channels(4)
         ]
 
         for i, roi in enumerate(rois):
@@ -265,6 +274,48 @@ def test_process_acquisition_geometry_parallel3D(self):
             self.assertEqual(ag_gold[i], ag_out, msg="Binning acquisition geometry with roi {}".format(i))
 
 
+    def test_process_acquisition_geometry_parallel3D_origin(self):
+
+        #tests the geometry output with a non-default origin choice
+
+        angles_full = numpy.linspace(0,360,360,endpoint=False)
+        ag = AcquisitionGeometry.create_Parallel3D().set_angles(angles_full).set_panel([128,64],[0.1,0.2], origin='top-right').set_channels(4)
+
+        rois = [
+                # same as input
+                {'channel':(None,None,None),'angle':(None,None,None),'vertical':(None,None,None),'horizontal':(None,None,None)},
+
+                # bin all
+                {'channel':(None,None,4),'angle':(None,None,2),'vertical':(None,None,8),'horizontal':(None,None,16)},
+                
+                # bin to single dimension
+                {'vertical':(31,33,2)},
+
+        
+                # crop  asymmetrically
+                {'vertical':(10,None,None),'horizontal':(None,-20,None)}                
+        ]
+
+        #calculate offsets for geometry4
+        ag4_offset_h = ag.pixel_size_h*20/2
+        ag4_offset_v = -ag.pixel_size_v*10/2
+
+        ag_gold = [
+                ag.copy(),
+                AcquisitionGeometry.create_Parallel3D().set_angles(numpy.linspace(0.5,360.5,180,endpoint=False)).set_panel([8,8],[1.6,1.6], origin='top-right').set_channels(1),
+                AcquisitionGeometry.create_Parallel2D().set_angles(angles_full).set_panel(128,[0.1,0.4], origin='top-right').set_channels(4),
+                AcquisitionGeometry.create_Parallel3D(detector_position=[ag4_offset_h, 0, ag4_offset_v]).set_angles(angles_full).set_panel([108,54],[0.1,0.2], origin='top-right').set_channels(4)
+        ]
+
+        for i, roi in enumerate(rois):
+            proc = Binner(roi=roi)
+            proc.set_input(ag)
+            ag_out = proc._process_acquisition_geometry()
+
+            self.assertEqual(ag_gold[i], ag_out, msg="Binning acquisition geometry with roi {0}. \nExpected:\n{1}\nGot\n{2}".format(i,ag_gold[i], ag_out))
+
+
+
     def test_process_acquisition_geometry_cone2D(self):
 
         ag = AcquisitionGeometry.create_Cone2D([0,-50],[0,50]).set_angles(numpy.linspace(0,360,360,endpoint=False)).set_panel(128,0.1).set_channels(4)
@@ -680,6 +731,41 @@ def test_aqdata_full(self):
         # not a very tight tolerance as binning and fp at lower res are not identical operations.
         numpy.testing.assert_allclose(fp_roi.array, fp_binned.array, atol=0.06)
 
+
+    @unittest.skipUnless(has_astra and has_nvidia, "ASTRA GPU not installed")
+    def test_aqdata_full_origin(self):
+        """
+        This test bins a sinogram. It then uses that geometry for the forward projection.
+
+        This ensures the offsets are correctly set and the same window of data is output in both cases.
+        """
+        ag = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get().geometry
+        ag.set_labels(['vertical','angle','horizontal'])
+        ag.config.panel.origin='top-right'
+    
+        phantom = dataexample.SIMULATED_SPHERE_VOLUME.get()
+
+        PO = AstraProjectionOperator(phantom.geometry, ag)
+        fp_full = PO.direct(phantom)
+
+        roi = {'angle':(25,26),'vertical':(5,62,2),'horizontal':(-75,0,2)}
+        binner = Binner(roi)
+
+        binner.set_input(ag)
+        ag_roi = binner.get_output()
+
+        binner.set_input(fp_full)
+        fp_binned = binner.get_output()
+
+        self.assertEqual(ag_roi, fp_binned.geometry, msg="Binned geometries not equal")
+
+        PO = AstraProjectionOperator(phantom.geometry, ag_roi)
+        fp_roi = PO.direct(phantom)
+
+        # not a very tight tolerance as binning and fp at lower res are not identical operations.
+        numpy.testing.assert_allclose(fp_roi.array, fp_binned.array, atol=0.08)
+
+
     @unittest.skip
     @unittest.skipUnless(has_tigre and has_nvidia, "TIGRE GPU not installed")
     def test_aqdata_full_tigre(self):
@@ -784,7 +870,8 @@ def test_process_acquisition_geometry_parallel2D(self):
 
     def test_process_acquisition_geometry_parallel3D(self):
 
-        ag = AcquisitionGeometry.create_Parallel3D().set_angles(numpy.linspace(0,360,360,endpoint=False)).set_panel([128,64],[0.1,0.2]).set_channels(4)
+        angles_full = numpy.linspace(0,360,360,endpoint=False)
+        ag = AcquisitionGeometry.create_Parallel3D().set_angles(angles_full).set_panel([128,64],[0.1,0.2], origin='bottom-left').set_channels(4)
 
         rois = [
                 # same as input
@@ -795,22 +882,81 @@ def test_process_acquisition_geometry_parallel3D(self):
                 
                 # slice to single dimension
                 {'vertical':(31,33,2)},
+        
+                # slice  asymmetrically
+                {'vertical':(10,None,None),'horizontal':(None,-20,None)}
         ]
 
+        #calculate offsets for geometry2
         pix_end_h1 = ag.pixel_num_h -1
         pix_start_h2 = 0
         pix_end_h2 = 7 * 16 # last pixel index sliced multiplied by step size
-        offset_h = ag.pixel_size_h*((pix_start_h2)-(pix_end_h1-pix_end_h2))/2
+        ag2_offset_h = ag.pixel_size_h*((pix_start_h2)-(pix_end_h1-pix_end_h2))/2
 
         pix_end_v1 = ag.pixel_num_v -1
         pix_start_v2 = 0
         pix_end_v2 = 7 * 8 # last pixel index sliced multiplied by step size
-        offset_v = ag.pixel_size_v*((pix_start_v2)-(pix_end_v1- pix_end_v2))/2
+        ag2_offset_v = ag.pixel_size_v*((pix_start_v2)-(pix_end_v1- pix_end_v2))/2
+
+        #calculate offsets for geometry4
+        ag4_offset_h = -ag.pixel_size_h*20/2
+        ag4_offset_v = ag.pixel_size_v*10/2
 
         ag_gold = [
                 ag.copy(),
-                AcquisitionGeometry.create_Parallel3D(detector_position=[offset_h, 0, offset_v]).set_angles(numpy.linspace(0,360,180,endpoint=False)).set_panel([8,8],[1.6,1.6]).set_channels(1),
+                AcquisitionGeometry.create_Parallel3D(detector_position=[ag2_offset_h, 0, ag2_offset_v]).set_angles(numpy.linspace(0,360,180,endpoint=False)).set_panel([8,8],[1.6,1.6]).set_channels(1),
                 AcquisitionGeometry.create_Parallel2D().set_angles(numpy.linspace(0,360,360,endpoint=False)).set_panel(128,[0.1,0.4]).set_channels(4),        
+                AcquisitionGeometry.create_Parallel3D(detector_position=[ag4_offset_h, 0, ag4_offset_v]).set_angles(angles_full).set_panel([108,54],[0.1,0.2]).set_channels(4)
+        ]
+
+        for i, roi in enumerate(rois):
+            proc = Slicer(roi=roi)
+            proc.set_input(ag)
+            ag_out = proc._process_acquisition_geometry()
+
+            self.assertEqual(ag_gold[i], ag_out, msg="Slicing acquisition geometry with roi {0}. \nExpected:\n{1}\nGot\n{2}".format(i,ag_gold[i], ag_out))
+
+    def test_process_acquisition_geometry_parallel3D_origin(self):
+
+        #tests the geometry output with a non-default origin choice
+
+        angles_full = numpy.linspace(0,360,360,endpoint=False)
+        ag = AcquisitionGeometry.create_Parallel3D().set_angles(angles_full).set_panel([128,64],[0.1,0.2], origin='top-right').set_channels(4)
+
+        rois = [
+                # same as input
+                {'channel':(None,None,None),'angle':(None,None,None),'vertical':(None,None,None),'horizontal':(None,None,None)},
+
+                # slice all
+                {'channel':(None,None,4),'angle':(None,None,2),'vertical':(None,None,8),'horizontal':(None,None,16)},
+                
+                # slice to single dimension
+                {'vertical':(31,33,2)},
+        
+                # slice  asymmetrically
+                {'vertical':(10,None,None),'horizontal':(None,-20,None)}
+        ]
+
+        #calculate offsets for geometry2
+        pix_end_h1 = ag.pixel_num_h -1
+        pix_start_h2 = 0
+        pix_end_h2 = 7 * 16 # last pixel index sliced multiplied by step size
+        ag2_offset_h = -ag.pixel_size_h*((pix_start_h2)-(pix_end_h1-pix_end_h2))/2
+
+        pix_end_v1 = ag.pixel_num_v -1
+        pix_start_v2 = 0
+        pix_end_v2 = 7 * 8 # last pixel index sliced multiplied by step size
+        ag2_offset_v = -ag.pixel_size_v*((pix_start_v2)-(pix_end_v1- pix_end_v2))/2
+
+        #calculate offsets for geometry4
+        ag4_offset_h = ag.pixel_size_h*20/2
+        ag4_offset_v = -ag.pixel_size_v*10/2
+
+        ag_gold = [
+                ag.copy(),
+                AcquisitionGeometry.create_Parallel3D(detector_position=[ag2_offset_h, 0, ag2_offset_v]).set_angles(numpy.linspace(0,360,180,endpoint=False)).set_panel([8,8],[1.6,1.6], origin='top-right').set_channels(1),
+                AcquisitionGeometry.create_Parallel2D().set_angles(numpy.linspace(0,360,360,endpoint=False)).set_panel(128,[0.1,0.4], origin='top-right').set_channels(4),        
+                AcquisitionGeometry.create_Parallel3D(detector_position=[ag4_offset_h, 0, ag4_offset_v]).set_angles(angles_full).set_panel([108,54],[0.1,0.2], origin='top-right').set_channels(4)
         ]
 
         for i, roi in enumerate(rois):
@@ -1119,6 +1265,42 @@ def test_aqdata_full(self):
 
         numpy.testing.assert_allclose(fp_roi.array, fp_sliced.array, 1e-4)
 
+
+    @unittest.skipUnless(has_astra and has_nvidia, "ASTRA GPU not installed")
+    def test_aqdata_full_origin(self):
+        """
+        This test slices a sinogram. It then uses that geometry for the forward projection.
+
+        This ensures the offsets are correctly set and the same window of data is output in both cases.
+        """
+
+        ag = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get().geometry
+        ag.set_labels(['vertical','angle','horizontal'])
+        ag.config.panel.origin='top-right'
+
+        phantom = dataexample.SIMULATED_SPHERE_VOLUME.get()
+
+        PO = AstraProjectionOperator(phantom.geometry, ag)
+        fp_full = PO.direct(phantom)
+
+        roi = {'angle':(25,30,2),'vertical':(5,50,2),'horizontal':(-50,0,2)}
+        slicer = Slicer(roi)
+
+        slicer.set_input(ag)
+        ag_roi = slicer.get_output()
+
+        slicer.set_input(fp_full)
+        fp_sliced = slicer.get_output()
+
+        numpy.testing.assert_allclose(fp_sliced.array, fp_full.array[5:50:2,25:30:2,-50::2])
+
+        self.assertEqual(ag_roi, fp_sliced.geometry, msg="Sliced geometries not equal")
+
+        PO = AstraProjectionOperator(phantom.geometry, ag_roi)
+        fp_roi = PO.direct(phantom)
+        numpy.testing.assert_allclose(fp_roi.array, fp_sliced.array, 1e-4)
+
+
     @unittest.skip
     @unittest.skipUnless(has_tigre and has_nvidia, "TIGRE GPU not installed")
     def test_aqdata_full_tigre(self):
@@ -1359,6 +1541,9 @@ def setUp(self):
         self.ag_pad_width = {'channel':(1,2),'vertical':(3,4),'horizontal':(5,6)}
         self.ag_padded = AcquisitionGeometry.create_Parallel3D(detector_position=[-0.05, 0., -0.15]).set_angles([0,90,180,270]).set_panel([27,23],[0.1,0.1]).set_channels(7)
 
+        self.ag2 = AcquisitionGeometry.create_Parallel3D(detector_position=[-0.1, 0.,-0.2]).set_angles([0,90,180,270]).set_panel([16,16],[0.1,0.1],origin='top-right').set_channels(4)
+        self.ag2_padded = AcquisitionGeometry.create_Parallel3D(detector_position=[-0.15, 0., -0.25]).set_angles([0,90,180,270]).set_panel([27,23],[0.1,0.1],origin='top-right').set_channels(7)
+
         self.ig = ImageGeometry(5,4,3,center_x=0.5,center_y=1,center_z=-0.5,channels=2)
         self.ig_pad_width = {'channel':(1,2),'vertical':(3,2),'horizontal_x':(2,1), 'horizontal_y':(2,3)}
         self.ig_padded = ImageGeometry(8,9,8,center_x=0,center_y=1.5,center_z=-1, channels=5)
@@ -1487,6 +1672,17 @@ def test_process_acquisition_geometry(self):
         msg="Padder failed with geometry mismatch. Got:\n{0}\nExpected:\n{1}".format(geometry_padded, geometry_gold))
 
 
+    def test_process_acquisition_geometry_origin(self):
+        geometry = self.ag2
+
+        proc = Padder('constant', pad_width=self.ag_pad_width, pad_values=0.0)
+        proc.set_input(geometry)
+        geometry_padded = proc._process_acquisition_geometry()
+
+        self.assertEquals(geometry_padded, self.ag2_padded,
+        msg="Padder failed with geometry mismatch. Got:\n{0}\nExpected:\n{1}".format(geometry_padded, self.ag2_padded))
+
+    
     def test_process_image_geometry(self):
 
         geometry = self.ig
@@ -1830,6 +2026,22 @@ def test_pad_ad_full(self):
         numpy.testing.assert_allclose(recon_orig.array, recon_new.array, atol=1e-4)
 
 
+        # switch panel origin
+        data.geometry.config.panel.origin='top-right'
+
+        recon_orig = FBP(data).run(verbose=0)
+        recon_orig.apply_circular_mask()
+
+        proc = Padder('constant',pad_width=(5,40),pad_values=0.0)
+        proc.set_input(data)
+        data_padded = proc.get_output()
+
+        recon_new = FBP(data_padded, recon_orig.geometry).run(verbose=0)
+        recon_new.apply_circular_mask()
+
+        numpy.testing.assert_allclose(recon_orig.array, recon_new.array, atol=1e-4)
+
+
     @unittest.skipUnless(has_astra and has_nvidia, "ASTRA GPU not installed")
     def test_pad_id_full(self):
         """

From fb95d17772aea6d6e1f9b0e3f5642afafebc0d1e Mon Sep 17 00:00:00 2001
From: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
Date: Fri, 29 Sep 2023 17:09:59 +0100
Subject: [PATCH 3/9] Pre-release update for v21.1.0 (#1511)

* Update changelog

* Update readme

* build variants
---
 CHANGELOG.md                   |  2 +-
 README.md                      |  4 ++--
 recipe/conda_build_config.yaml | 22 +++++++++-------------
 3 files changed, 12 insertions(+), 16 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2ea88e4ce1..349e90a9c7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,5 @@
 
-* x.x.x
+* 23.1.0
   - Fix bug in IndicatorBox proximal_conjugate
   - Allow CCPi Regulariser functions for non CIL object
   - Add norm for CompositionOperator
diff --git a/README.md b/README.md
index a3ed88f8c5..dc8d5d1ca4 100644
--- a/README.md
+++ b/README.md
@@ -19,13 +19,13 @@ The documentation for CIL can be accessed [here](https://tomographicimaging.gith
 Binary installation of CIL can be done with `conda`. Install a new environment using:
 
 ```bash
-conda create --name cil -c conda-forge -c intel -c ccpi cil=23.0.1
+conda create --name cil -c conda-forge -c intel -c ccpi cil=23.1.0
 ```
 
 To install CIL and the additional packages and plugins needed to run the [CIL demos](https://github.com/TomographicImaging/CIL-Demos) install the environment with:
 
 ```bash
-conda create --name cil -c conda-forge -c intel -c ccpi cil=23.0.1 astra-toolbox tigre ccpi-regulariser tomophantom "ipywidgets<8"
+conda create --name cil -c conda-forge -c intel -c ccpi cil=23.1.0 astra-toolbox tigre ccpi-regulariser tomophantom "ipywidgets<8"
 ```
 where,
 
diff --git a/recipe/conda_build_config.yaml b/recipe/conda_build_config.yaml
index cd7d480cef..2ba2a05f82 100644
--- a/recipe/conda_build_config.yaml
+++ b/recipe/conda_build_config.yaml
@@ -22,18 +22,16 @@
 python:
   - 3.8
   - 3.9
-  - 3.10
+  - 3.10 
   - 3.8
   - 3.9
-  - 3.10
+  - 3.10 
   - 3.8
   - 3.9
-  - 3.10
-  # - 3.11
-  # - 3.8 
-  # - 3.9
-  # - 3.10 
-  # - 3.11
+  - 3.10 
+  - 3.8
+  - 3.9
+  - 3.10 
 numpy:
   - 1.21
   - 1.21
@@ -44,11 +42,9 @@ numpy:
   - 1.23
   - 1.23
   - 1.23
-  # - 1.23
-  # - 1.24
-  # - 1.24
-  # - 1.24
-  # - 1.24
+  - 1.24
+  - 1.24
+  - 1.24
 zip_keys:
   - python
   - numpy

From 224a8ea87cfcf4bc0691325c2614bf110bcd48e0 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <casper.dcl@physics.org>
Date: Tue, 3 Oct 2023 10:39:36 +0100
Subject: [PATCH 4/9] ci: update workflow

---
 .github/workflows/conda_and_docs_build.yml | 72 ++++++++--------------
 1 file changed, 27 insertions(+), 45 deletions(-)

diff --git a/.github/workflows/conda_and_docs_build.yml b/.github/workflows/conda_and_docs_build.yml
index 1fe806f25e..c1bf4ba240 100644
--- a/.github/workflows/conda_and_docs_build.yml
+++ b/.github/workflows/conda_and_docs_build.yml
@@ -1,29 +1,27 @@
 #  -*- coding: utf-8 -*-
 # Copyright 2021 United Kingdom Research and Innovation
 # Copyright 2021 The University of Manchester
-    
+
 # 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.
-    
+
 # Authors:
 # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
-
 name: conda_and_docs_build
-
 on:
   release:
     types: [published]
   push:
-    branches: [ master ]
+    branches: [master]
     tags:
       - '**'
     paths-ignore:
@@ -33,9 +31,8 @@ on:
     - 'scripts/**'
     - 'NOTICE.txt'
     - 'README.md'
-
   pull_request:
-    branches: [ master ]
+    branches: [master]
     paths-ignore:
     - 'CHANGELOG.md'
     - 'CITATION.cff'
@@ -43,25 +40,24 @@ on:
     - 'scripts/**'
     - 'NOTICE.txt'
     - 'README.md'
-
 jobs:
   conda_build:
     runs-on: ubuntu-20.04
     steps:
-    - uses: actions/checkout@v3.1.0
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 0
     - name: conda-build
       uses: paskino/conda-package-publish-action@v1.4.4
       with:
-        subDir: 'recipe'
-        channels: '-c conda-forge -c intel -c astra-toolbox -c ccpi'
+        subDir: recipe
+        channels: -c conda-forge -c intel -c astra-toolbox -c ccpi
         convert_win: false
         convert_osx: false
         test_pyver: 3.9
         test_npver: 1.22
     - name: Upload artifact of the conda package.
-      uses: actions/upload-artifact@v3.1.1
+      uses: actions/upload-artifact@v3
       with:
         name: cil-package
         path: recipe/linux-64/cil*
@@ -69,56 +65,42 @@ jobs:
     needs: conda_build
     runs-on: ubuntu-20.04
     steps:
-    - uses: actions/checkout@v3.1.0
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 0
-    - name: change directory
-      run: |
-        ls
+    - run: ls
     - name: Download artifact of the conda package.
-      uses: actions/download-artifact@v3.0.1
+      uses: actions/download-artifact@v3
       with:
-        name: 'cil-package'
+        name: cil-package
         path: 'conda_package'
     - uses: conda-incubator/setup-miniconda@v2
       with:
         python-version: 3.9
     - uses: lauramurgatroyd/build-sphinx-action@v0.1.4
       with:
-        DOCS_PATH: 'docs'
-        CONDA_BUILD_ENV_FILEPATH: 'docs/docs_environment.yml'
-        ARTIFACT_NAME: 'DocumentationHTML'
-        PACKAGE_FOLDER_PATH: 'conda_package'
-        PACKAGE_NAME: 'cil'
-        PACKAGE_CONDA_CHANNELS: 'conda-forge -c intel -c astra-toolbox -c ccpi'
-        BUILD_SUBDIR_NAME: 'nightly'
+        DOCS_PATH: docs
+        CONDA_BUILD_ENV_FILEPATH: docs/docs_environment.yml
+        ARTIFACT_NAME: DocumentationHTML
+        PACKAGE_FOLDER_PATH: conda_package
+        PACKAGE_NAME: cil
+        PACKAGE_CONDA_CHANNELS: conda-forge -c intel -c astra-toolbox -c ccpi
+        BUILD_SUBDIR_NAME: nightly
         python_version: 3.9
   docs_publish:
     needs: docs_build
-    runs-on: ubuntu-latest
     if: github.ref == 'refs/heads/master'
+    runs-on: ubuntu-latest
     steps:
     - name: Download artifact of the html output.
-      uses: actions/download-artifact@v3.0.1
+      uses: actions/download-artifact@v3
       with:
         name: DocumentationHTML
         path: docs/build
-    - name: Commit documentation changes
-      run: |
-        git clone https://github.com/TomographicImaging/CIL.git --branch gh-pages --single-branch gh-pages
-        cp -r docs/build/* gh-pages/
-        cd gh-pages
-        touch .nojekyll
-        git config --local user.email "action@github.com"
-        git config --local user.name "GitHub Action"
-        git add .
-        git commit -m "Update documentation" -a || true
-        # The above command will fail if no changes were present, so we ignore
-        # that.
     - name: Push changes
-      uses: ad-m/github-push-action@v0.6.0
+      uses: casperdcl/push-dir@v1
       with:
+        message: Update documentation
         branch: gh-pages
-        directory: gh-pages
-        github_token: ${{ secrets.GITHUB_TOKEN }}
-
+        dir: docs/build
+        nojekyll: true

From 5ab3f16549ca56ac6869e2a0a4ab7d9675df4386 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <casper.dcl@physics.org>
Date: Tue, 3 Oct 2023 14:10:47 +0100
Subject: [PATCH 5/9] drop -c astra-toolbox

Co-authored-by: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
Signed-off-by: Casper da Costa-Luis <casper.dcl@physics.org>
---
 .github/workflows/conda_and_docs_build.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/conda_and_docs_build.yml b/.github/workflows/conda_and_docs_build.yml
index c1bf4ba240..df63bd32c3 100644
--- a/.github/workflows/conda_and_docs_build.yml
+++ b/.github/workflows/conda_and_docs_build.yml
@@ -51,7 +51,7 @@ jobs:
       uses: paskino/conda-package-publish-action@v1.4.4
       with:
         subDir: recipe
-        channels: -c conda-forge -c intel -c astra-toolbox -c ccpi
+        channels: -c conda-forge -c intel -c ccpi
         convert_win: false
         convert_osx: false
         test_pyver: 3.9
@@ -84,7 +84,7 @@ jobs:
         ARTIFACT_NAME: DocumentationHTML
         PACKAGE_FOLDER_PATH: conda_package
         PACKAGE_NAME: cil
-        PACKAGE_CONDA_CHANNELS: conda-forge -c intel -c astra-toolbox -c ccpi
+        PACKAGE_CONDA_CHANNELS: conda-forge -c intel -c ccpi
         BUILD_SUBDIR_NAME: nightly
         python_version: 3.9
   docs_publish:

From 0df6e837c1d821e3b88e8f823d7881dbf17ff481 Mon Sep 17 00:00:00 2001
From: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
Date: Tue, 3 Oct 2023 14:15:13 +0100
Subject: [PATCH 6/9] Update README.md with libmamba instructions (#1519)

* Update README.md

Signed-off-by: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
---
 README.md | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index dc8d5d1ca4..f3140200c0 100644
--- a/README.md
+++ b/README.md
@@ -16,7 +16,14 @@ The documentation for CIL can be accessed [here](https://tomographicimaging.gith
 
 # Installation of CIL
 
-Binary installation of CIL can be done with `conda`. Install a new environment using:
+Binary installation of CIL can be achieved with `conda` or `mamba`. `mamba`'s environment solver (`libmamba`) provides faster environment resolution so we recommend this route, although in many cases the default `conda` will still work but may be very slow.
+
+`miniconda` is a minimal installer for `conda`. Installation instructions can be found [here](https://docs.conda.io/projects/miniconda/en/latest/).
+You can then force your `conda` installation to use `libmamba` instructions are [here](https://www.anaconda.com/blog/a-faster-conda-for-a-growing-community).
+
+Alternatively, `mamba` can be installed via `https://github.com/conda-forge/miniforge` which is another minimal installer for `conda` with optional support for `mamba`. In this case replace `conda` with `mamba` in the below commands.
+
+Install a new environment using:
 
 ```bash
 conda create --name cil -c conda-forge -c intel -c ccpi cil=23.1.0

From fc49977571a07a2bebfabe4d165350bee636f24d Mon Sep 17 00:00:00 2001
From: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
Date: Tue, 3 Oct 2023 14:16:43 +0100
Subject: [PATCH 7/9] Update README.md

Signed-off-by: Gemma Fardell <47746591+gfardell@users.noreply.github.com>
---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index f3140200c0..15a7d816df 100644
--- a/README.md
+++ b/README.md
@@ -21,7 +21,7 @@ Binary installation of CIL can be achieved with `conda` or `mamba`. `mamba`'s en
 `miniconda` is a minimal installer for `conda`. Installation instructions can be found [here](https://docs.conda.io/projects/miniconda/en/latest/).
 You can then force your `conda` installation to use `libmamba` instructions are [here](https://www.anaconda.com/blog/a-faster-conda-for-a-growing-community).
 
-Alternatively, `mamba` can be installed via `https://github.com/conda-forge/miniforge` which is another minimal installer for `conda` with optional support for `mamba`. In this case replace `conda` with `mamba` in the below commands.
+Alternatively, `mamba` can be installed via [`miniforge`](https://github.com/conda-forge/miniforge) which is another minimal installer for `conda` with optional support for `mamba`. In this case replace `conda` with `mamba` in the below commands.
 
 Install a new environment using:
 

From c3dd1b6cbd5ee95571a52bbbe34b7ad593ce15f9 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <casper.dcl@physics.org>
Date: Thu, 5 Oct 2023 08:50:05 +0100
Subject: [PATCH 8/9] ci: fix docs push (#1521)

---
 .github/workflows/conda_and_docs_build.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/conda_and_docs_build.yml b/.github/workflows/conda_and_docs_build.yml
index df63bd32c3..ae2f9558c3 100644
--- a/.github/workflows/conda_and_docs_build.yml
+++ b/.github/workflows/conda_and_docs_build.yml
@@ -92,6 +92,7 @@ jobs:
     if: github.ref == 'refs/heads/master'
     runs-on: ubuntu-latest
     steps:
+    - uses: actions/checkout@v4
     - name: Download artifact of the html output.
       uses: actions/download-artifact@v3
       with:

From 2b3e2cd393c7ae15e27a7539c144c5197afbacef Mon Sep 17 00:00:00 2001
From: Laura Murgatroyd <60604372+lauramurgatroyd@users.noreply.github.com>
Date: Thu, 5 Oct 2023 17:14:57 +0100
Subject: [PATCH 9/9] Update script for creating CIL test environment (#1524)

* update order of command and versions

* Update scripts/create_local_env_for_cil_development_tests.sh

Signed-off-by: Laura Murgatroyd <60604372+lauramurgatroyd@users.noreply.github.com>

---------

Signed-off-by: Laura Murgatroyd <60604372+lauramurgatroyd@users.noreply.github.com>
---
 scripts/create_local_env_for_cil_development_tests.sh | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/scripts/create_local_env_for_cil_development_tests.sh b/scripts/create_local_env_for_cil_development_tests.sh
index 6086d342c1..71519b7e21 100644
--- a/scripts/create_local_env_for_cil_development_tests.sh
+++ b/scripts/create_local_env_for_cil_development_tests.sh
@@ -5,7 +5,7 @@
 
 name=cil
 python=3.8
-numpy=1.20
+numpy=1.21
 
 while getopts h:n:p:e:v: option
  do
@@ -41,12 +41,12 @@ fi
 
 set -x 
 
-${conda_cmd} -c conda-forge -c intel  -c ccpi/label/dev -c ccpi -c astra-toolbox -c astra-toolbox/label/dev \
-        python=$python numpy=$numpy \
-        cil-data tigre=2.2 ccpi-regulariser=21.0.0 tomophantom=2.0.0  astra-toolbox'>=1.9.9.dev5,<2.1' \
+${conda_cmd} python=$python numpy=$numpy \
+        cil-data tigre=2.4 ccpi-regulariser=22.0.0 tomophantom=2.0.0  astra-toolbox'>=1.9.9.dev5,<2.1' \
         cvxpy python-wget scikit-image packaging \
         cmake'>=3.16' setuptools  \
         ipp-include ipp-devel ipp \
         ipywidgets scipy matplotlib \
         h5py pillow libgcc-ng dxchange olefile pywavelets numba tqdm \
+        -c conda-forge -c intel  -c ccpi/label/dev -c ccpi -c astra-toolbox -c astra-toolbox/label/dev \
         --override-channels