From ed8462641703ec7b9862b1994d1b14cc919fb2b2 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 8 Apr 2025 13:48:10 +0100
Subject: [PATCH 01/25] initial changes on datacontainer and image data

Co-authored-by: MikeSullivan7 <MikeSullivan7@users.noreply.github.com>
---
 .../Python/cil/framework/data_container.py    |  5 +-
 Wrappers/Python/cil/framework/image_data.py   | 51 +++++++++++++------
 Wrappers/Python/test/torch_test.py            |  9 ++++
 3 files changed, 47 insertions(+), 18 deletions(-)
 create mode 100644 Wrappers/Python/test/torch_test.py

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 4febd7a449..b2997eadd5 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -93,13 +93,14 @@ def size(self):
     __container_priority__ = 1
     def __init__ (self, array, deep_copy=True, dimension_labels=None,
                   **kwargs):
-        if type(array) == numpy.ndarray:
+        #if type(array) == numpy.ndarray:
+        if hasattr(array, "__array_interface__"):
             if deep_copy:
                 self.array = array.copy()
             else:
                 self.array = array
         else:
-            raise TypeError('Array must be NumpyArray, passed {0}'\
+            raise TypeError('Array must be adhere to the array interface protocol {0}'\
                             .format(type(array)))
 
         #Don't set for derived classes
diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index 65cac02214..c8619a8729 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -15,11 +15,14 @@
 #
 # Authors:
 # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
-import numpy
+import numpy as np
 
 from .data_container import DataContainer
 from .labels import ImageDimension, Backend
 
+from array_api_compat import array_namespace # https://data-apis.org/array-api-compat/
+# import array interface as xp
+
 class ImageData(DataContainer):
     '''DataContainer for holding 2D or 3D DataContainer'''
     __container_priority__ = 1
@@ -46,27 +49,32 @@ def __init__(self,
                  deep_copy=False,
                  geometry=None,
                  **kwargs):
+        '''Default to NumPy array if array is not passed'''
 
-        dtype = kwargs.get('dtype', numpy.float32)
-
-
-        if geometry is None:
-            raise AttributeError("ImageData requires a geometry")
-
-
-        labels = kwargs.get('dimension_labels', None)
-        if labels is not None and labels != geometry.dimension_labels:
-                raise ValueError("Deprecated: 'dimension_labels' cannot be set with 'allocate()'. Use 'geometry.set_labels()' to modify the geometry before using allocate.")
-
+        # xp = get_namespace()
+        
+        # TODO: change
         if array is None:
-            array = numpy.empty(geometry.shape, dtype=dtype)
+            # defaults to a numpy array
+            xp = np
+            array = xp.empty(geometry.shape, dtype=dtype)
         elif issubclass(type(array) , DataContainer):
+            # this is a reference so we might make a mess modifying both ends
             array = array.as_array()
-        elif issubclass(type(array) , numpy.ndarray):
+            xp = array_namespace(array)
+        # elif issubclass(type(array) , numpy.ndarray):
             # remove singleton dimensions
-            array = numpy.squeeze(array)
+            # array = numpy.squeeze(array)
         else:
-            raise TypeError('array must be a CIL type DataContainer or numpy.ndarray got {}'.format(type(array)))
+            # Consider array as an object is compliant to the array API
+            # https://docs.scipy.org/doc/scipy-1.15.2/dev/api-dev/array_api.html 
+            # this might raise an exception but that's fine
+            xp = array_namespace(array)
+            # remove singleton dimensions
+            array = xp.squeeze(array)
+            
+        # else:
+        #     raise TypeError('array must be a CIL type DataContainer or numpy.ndarray got {}'.format(type(array)))
 
         if array.shape != geometry.shape:
             raise ValueError('Shape mismatch {} {}'.format(array.shape, geometry.shape))
@@ -74,6 +82,17 @@ def __init__(self,
         if array.ndim not in [2,3,4]:
             raise ValueError('Number of dimensions are not 2 or 3 or 4 : {0}'.format(array.ndim))
 
+
+        dtype = kwargs.get('dtype', xp.float32)
+
+        if geometry is None:
+            raise AttributeError("ImageData requires a geometry")
+
+
+        labels = kwargs.get('dimension_labels', None)
+        if labels is not None and labels != geometry.dimension_labels:
+                raise ValueError("Deprecated: 'dimension_labels' cannot be set with 'allocate()'. Use 'geometry.set_labels()' to modify the geometry before using allocate.")
+
         super(ImageData, self).__init__(array, deep_copy, geometry=geometry, **kwargs)
 
     def __eq__(self, other):
diff --git a/Wrappers/Python/test/torch_test.py b/Wrappers/Python/test/torch_test.py
new file mode 100644
index 0000000000..91be29215c
--- /dev/null
+++ b/Wrappers/Python/test/torch_test.py
@@ -0,0 +1,9 @@
+import torch
+from cil.framework import DataContainer
+
+a = torch.zeros([2, 4], dtype=torch.float32)
+
+print (a[0], a[1])
+
+cil_a = DataContainer
+

From c090068559ac5c8a90ab7f6696053928b37b63a8 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 8 Apr 2025 14:07:35 +0100
Subject: [PATCH 02/25] torch test working

---
 Wrappers/Python/cil/framework/data_container.py | 17 ++++++++++-------
 Wrappers/Python/test/torch_test.py              | 17 +++++++++++++++--
 2 files changed, 25 insertions(+), 9 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index b2997eadd5..f45ef43298 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -22,6 +22,7 @@
 from numbers import Number
 
 import numpy
+from array_api_compat import array_namespace
 
 from .cilacc import cilacc
 from cil.utilities.multiprocessing import NUM_THREADS
@@ -94,14 +95,16 @@ def size(self):
     def __init__ (self, array, deep_copy=True, dimension_labels=None,
                   **kwargs):
         #if type(array) == numpy.ndarray:
-        if hasattr(array, "__array_interface__"):
-            if deep_copy:
-                self.array = array.copy()
-            else:
-                self.array = array
+        # if hasattr(array, "__array_interface__"):
+        if deep_copy:
+            # self.array = array.copy()
+            xp = array_namespace(array)
+            self.array = xp.asarray(array, copy=True)
         else:
-            raise TypeError('Array must be adhere to the array interface protocol {0}'\
-                            .format(type(array)))
+            self.array = array
+        # else:
+        #     raise TypeError('Array must be adhere to the array interface protocol {0}'\
+        #                     .format(type(array)))
 
         #Don't set for derived classes
         if type(self) is DataContainer:
diff --git a/Wrappers/Python/test/torch_test.py b/Wrappers/Python/test/torch_test.py
index 91be29215c..d9a0179272 100644
--- a/Wrappers/Python/test/torch_test.py
+++ b/Wrappers/Python/test/torch_test.py
@@ -1,9 +1,22 @@
 import torch
+import numpy 
 from cil.framework import DataContainer
 
+
+b = numpy.ones([2,4], dtype=numpy.float32)
+cil_np0 = DataContainer(b, deep_copy=True)
+cil_np1 = DataContainer(b, deep_copy=False)
+
+numpy.testing.assert_allclose(numpy.zeros([2,4], dtype=numpy.float32), 
+                              (cil_np0-cil_np1).as_array())
 a = torch.zeros([2, 4], dtype=torch.float32)
 
-print (a[0], a[1])
+print (f"torch array has __array_interface__? {hasattr(a, '__array_interface__')}")
+# https://github.com/pytorch/pytorch/issues/54138
 
-cil_a = DataContainer
+# print (a[0], a[1])
+cil_tc0 = DataContainer(a, deep_copy=True)
+cil_tc1 = DataContainer(a, deep_copy=False)
 
+numpy.testing.assert_allclose(numpy.zeros([2,4], dtype=numpy.float32), 
+                              (cil_tc0-cil_tc1).as_array().numpy())

From dc7c6f0c798eff886339cc73bf01fa5d0c7d0b7d Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 8 Apr 2025 15:22:43 +0100
Subject: [PATCH 03/25] updates to the pixelwise_binary

Co-authored-by: MikeSullivan7 <MikeSullivan7@users.noreply.github.com>
Co-authored-by: Gemma Fardell <gfardell@users.noreply.github.com>
---
 .../Python/cil/framework/data_container.py    | 65 ++++++++++---------
 Wrappers/Python/test/torch_test.py            | 58 ++++++++++++++++-
 2 files changed, 91 insertions(+), 32 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index f45ef43298..f1e9fddf8c 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -414,10 +414,8 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
                 out = pwop(self.as_array() , x2 , *args, **kwargs )
             elif issubclass(x2.__class__ , DataContainer):
                 out = pwop(self.as_array() , x2.as_array() , *args, **kwargs )
-            elif isinstance(x2, numpy.ndarray):
-                out = pwop(self.as_array() , x2 , *args, **kwargs )
             else:
-                raise TypeError('Expected x2 type as number or DataContainer, got {}'.format(type(x2)))
+                out = pwop(self.as_array() , x2 , *args, **kwargs )
             geom = self.geometry
             if geom is not None:
                 geom = self.geometry.copy()
@@ -438,17 +436,12 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
                 #       geometry=self.geometry)
                 return out
             raise ValueError(f"Wrong size for data memory: out {out.shape} x2 {x2.shape} expected {self.shape}")
-        elif issubclass(type(out), DataContainer) and \
-             isinstance(x2, (Number, numpy.ndarray)):
-            if self.check_dimensions(out):
-                if isinstance(x2, numpy.ndarray) and\
-                    not (x2.shape == self.shape and x2.dtype == self.dtype):
-                    raise ValueError(f"Wrong size for data memory: out {out.shape} x2 {x2.shape} expected {self.shape}")
+        elif issubclass(type(out), DataContainer) and self.check_dimensions(out):
                 kwargs['out']=out.as_array()
                 pwop(self.as_array(), x2, *args, **kwargs )
                 return out
-            raise ValueError(f"Wrong size for data memory: {out.shape} {self.shape}")
-        elif issubclass(type(out), numpy.ndarray):
+        # elif issubclass(type(out), numpy.ndarray):
+        else:
             if self.array.shape == out.shape and self.array.dtype == out.dtype:
                 kwargs['out'] = out
                 pwop(self.as_array(), x2, *args, **kwargs)
@@ -456,41 +449,48 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
                 #       deep_copy=False,
                 #       dimension_labels=self.dimension_labels,
                 #       geometry=self.geometry)
-        else:
-            raise ValueError(f"incompatible class: {pwop.__name__} {type(out)}")
+        # else:
+        #     raise ValueError(f"incompatible class: {pwop.__name__} {type(out)}")
 
     def add(self, other, *args, **kwargs):
         if hasattr(other, '__container_priority__') and \
            self.__class__.__container_priority__ < other.__class__.__container_priority__:
             return other.add(self, *args, **kwargs)
-        return self.pixel_wise_binary(numpy.add, other, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_binary(xp.add, other, *args, **kwargs)
 
     def subtract(self, other, *args, **kwargs):
         if hasattr(other, '__container_priority__') and \
            self.__class__.__container_priority__ < other.__class__.__container_priority__:
             return other.subtract(self, *args, **kwargs)
-        return self.pixel_wise_binary(numpy.subtract, other, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_binary(xp.subtract, other, *args, **kwargs)
 
     def multiply(self, other, *args, **kwargs):
         if hasattr(other, '__container_priority__') and \
            self.__class__.__container_priority__ < other.__class__.__container_priority__:
             return other.multiply(self, *args, **kwargs)
-        return self.pixel_wise_binary(numpy.multiply, other, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_binary(xp.multiply, other, *args, **kwargs)
 
     def divide(self, other, *args, **kwargs):
         if hasattr(other, '__container_priority__') and \
            self.__class__.__container_priority__ < other.__class__.__container_priority__:
             return other.divide(self, *args, **kwargs)
-        return self.pixel_wise_binary(numpy.divide, other, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_binary(xp.divide, other, *args, **kwargs)
 
     def power(self, other, *args, **kwargs):
-        return self.pixel_wise_binary(numpy.power, other, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_binary(xp.power, other, *args, **kwargs)
 
     def maximum(self, x2, *args, **kwargs):
-        return self.pixel_wise_binary(numpy.maximum, x2, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_binary(xp.maximum, x2, *args, **kwargs)
 
     def minimum(self,x2, out=None, *args, **kwargs):
-        return self.pixel_wise_binary(numpy.minimum, x2=x2, out=out, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_binary(xp.minimum, x2=x2, out=out, *args, **kwargs)
 
 
     def sapyb(self, a, y, b, out=None, num_threads=NUM_THREADS):
@@ -658,32 +658,39 @@ def pixel_wise_unary(self, pwop, *args,  **kwargs):
                 pwop(self.as_array(), *args, **kwargs )
             else:
                 raise ValueError(f"Wrong size for data memory: {out.shape} {self.shape}")
-        elif issubclass(type(out), numpy.ndarray):
+        # elif issubclass(type(out), numpy.ndarray):
+        else:
             if self.array.shape == out.shape and self.array.dtype == out.dtype:
                 kwargs['out'] = out
                 pwop(self.as_array(), *args, **kwargs)
-        else:
-            raise ValueError("incompatible class: {pwop.__name__} {type(out)}")
+        # else:
+        #     raise ValueError("incompatible class: {pwop.__name__} {type(out)}")
 
     def abs(self, *args,  **kwargs):
-        return self.pixel_wise_unary(numpy.abs, *args,  **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_unary(xp.abs, *args,  **kwargs)
 
     def sign(self, *args,  **kwargs):
-        return self.pixel_wise_unary(numpy.sign, *args,  **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_unary(xp.sign, *args,  **kwargs)
 
     def sqrt(self, *args,  **kwargs):
-        return self.pixel_wise_unary(numpy.sqrt, *args,  **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_unary(xp.sqrt, *args,  **kwargs)
 
     def conjugate(self, *args,  **kwargs):
-        return self.pixel_wise_unary(numpy.conjugate, *args,  **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_unary(xp.conjugate, *args,  **kwargs)
 
     def exp(self, *args, **kwargs):
         '''Applies exp pixel-wise to the DataContainer'''
-        return self.pixel_wise_unary(numpy.exp, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_unary(xp.exp, *args, **kwargs)
 
     def log(self, *args, **kwargs):
         '''Applies log pixel-wise to the DataContainer'''
-        return self.pixel_wise_unary(numpy.log, *args, **kwargs)
+        xp = array_namespace(self.array)
+        return self.pixel_wise_unary(xp.log, *args, **kwargs)
 
     ## reductions
     def squared_norm(self, **kwargs):
diff --git a/Wrappers/Python/test/torch_test.py b/Wrappers/Python/test/torch_test.py
index d9a0179272..4fac0cc634 100644
--- a/Wrappers/Python/test/torch_test.py
+++ b/Wrappers/Python/test/torch_test.py
@@ -9,14 +9,66 @@
 
 numpy.testing.assert_allclose(numpy.zeros([2,4], dtype=numpy.float32), 
                               (cil_np0-cil_np1).as_array())
-a = torch.zeros([2, 4], dtype=torch.float32)
+
+
+
+a = [ torch.ones([2, 4], dtype=torch.float32),
+      torch.ones([2, 4], dtype=torch.float32),
+      torch.ones([2, 4], dtype=torch.float32)]
+
 
 print (f"torch array has __array_interface__? {hasattr(a, '__array_interface__')}")
 # https://github.com/pytorch/pytorch/issues/54138
 
 # print (a[0], a[1])
-cil_tc0 = DataContainer(a, deep_copy=True)
-cil_tc1 = DataContainer(a, deep_copy=False)
+cil_tc0 = DataContainer(a[0], deep_copy=True)
+cil_tc1 = DataContainer(a[0], deep_copy=False)
 
 numpy.testing.assert_allclose(numpy.zeros([2,4], dtype=numpy.float32), 
                               (cil_tc0-cil_tc1).as_array().numpy())
+
+# Torch - NumPy = Torch
+c = a[0] - b
+print(c)
+
+# NumPy - Torch = TypeError: unsupported operand type(s) for -: 'numpy.ndarray' and 'Tensor'
+# d = b - a = TypeError: unsupported operand type(s) for -: 'numpy.ndarray' and 'Tensor'
+# in CIL we do subtraction as below and it works
+d = numpy.subtract(b,a[0])
+print("d", d)
+
+# CIL DataContainer and mixed stuff
+mixed_diff = cil_np0.subtract(cil_tc0)
+print (mixed_diff.array)
+print(f"Type of mixed_diff cil_np0.subtract(cil_tc0) {type(mixed_diff.array)}")
+numpy.testing.assert_allclose(numpy.zeros([2,4], dtype=numpy.float32), 
+                              (mixed_diff).as_array().numpy())
+
+cil_tc2 = DataContainer(a[1], deep_copy=False)
+print(f"cil_tc2.array ype of mixed_diff {type(cil_tc2.array)}")
+
+# cil_tc0.subtract(cil_tc1, out=cil_tc2.array)
+from array_api_compat import array_namespace
+# array_namespace will raise a TypeError if multiple namespaces for array input
+xp = array_namespace(*a)
+print(xp)
+xp.subtract(a[1], a[0], out=a[2])
+print ("after subtract", a[2])
+
+cil_tc2 = DataContainer(a[2], deep_copy=False)
+cil_tc2 += 1
+
+print("################")
+print (cil_tc2.array)
+cil_tc0.subtract(cil_tc1, out=cil_tc2)
+
+print (cil_tc2.array)
+# from array_api_compat import array_namespace
+# array_namespace([1,2])
+cil_tc2+=1
+# mixed types with output
+# output must match the type of the first array
+print("################")
+print (cil_np1.array)
+mixed_diff = cil_np0.subtract(cil_tc0, out=cil_np1)
+print (cil_np1.array)
\ No newline at end of file

From 07b00cc38f208209f13bcd27427afd55fd69830b Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 8 Apr 2025 16:06:44 +0100
Subject: [PATCH 04/25] sapyb

---
 .../Python/cil/framework/data_container.py    | 28 +++++++++------
 Wrappers/Python/test/torch_test.py            | 34 ++++++++++++++++++-
 2 files changed, 50 insertions(+), 12 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index f1e9fddf8c..88642a119c 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -527,9 +527,11 @@ def sapyb(self, a, y, b, out=None, num_threads=NUM_THREADS):
                 self._axpby(a, b, y, out, out.dtype, num_threads)
                 return out
             except RuntimeError as rte:
-                warnings.warn("sapyb defaulting to Python due to: {}".format(rte))
+                warnings.warn("sapyb defaulting to NumPy due to: {}".format(rte))
             except TypeError as te:
-                warnings.warn("sapyb defaulting to Python due to: {}".format(te))
+                warnings.warn("sapyb defaulting to NumPy due to: {}".format(te))
+            except ValueError as ve:
+                warnings.warn("sapyb defaulting to NumPy due to: {}".format(ve))
             finally:
                 pass
 
@@ -561,15 +563,19 @@ def _axpby(self, a, b, y, out, dtype=numpy.float32, num_threads=NUM_THREADS):
         c_double_p = ctypes.POINTER(ctypes.c_double)
 
         #convert a and b to numpy arrays and get the reference to the data (length = 1 or ndx.size)
+        # CAREFUL
+        # this will fail if the numpy array cannot be created without copying to RAM
         try:
-            nda = a.as_array()
-        except:
-            nda = numpy.asarray(a)
+            a = a.as_array()
+        except AttributeError:
+            pass
+        nda = numpy.asarray(a, copy=False)
 
         try:
-            ndb = b.as_array()
-        except:
-            ndb = numpy.asarray(b)
+            b = b.as_array()
+        except AttributeError:
+            pass
+        ndb = numpy.asarray(b, copy=False)
 
         a_vec = 0
         if nda.size > 1:
@@ -580,9 +586,9 @@ def _axpby(self, a, b, y, out, dtype=numpy.float32, num_threads=NUM_THREADS):
             b_vec = 1
 
         # get the reference to the data
-        ndx = self.as_array()
-        ndy = y.as_array()
-        ndout = out.as_array()
+        ndx   = numpy.asarray(self.as_array(), copy=False)
+        ndy   = numpy.asarray(y.as_array(),    copy=False)
+        ndout = numpy.asarray(out.as_array(),  copy=False)
 
         if ndout.dtype != dtype:
             raise Warning("out array of type {0} does not match requested dtype {1}. Using {0}".format(ndout.dtype, dtype))
diff --git a/Wrappers/Python/test/torch_test.py b/Wrappers/Python/test/torch_test.py
index 4fac0cc634..828e31bad6 100644
--- a/Wrappers/Python/test/torch_test.py
+++ b/Wrappers/Python/test/torch_test.py
@@ -71,4 +71,36 @@
 print("################")
 print (cil_np1.array)
 mixed_diff = cil_np0.subtract(cil_tc0, out=cil_np1)
-print (cil_np1.array)
\ No newline at end of file
+print (cil_np1.array)
+
+print("########## test sapyb with pytorch ########## ")
+a = [ torch.ones([2, 4], dtype=torch.float32),
+      torch.ones([2, 4], dtype=torch.float32),
+      torch.ones([2, 4], dtype=torch.float32)]
+
+cil_tc = [ DataContainer(el, deep_copy=False) for el in a]
+
+# 2 * 1 + 3 * 1
+res = cil_tc[1].sapyb(2, cil_tc[1], 3)
+cil_tc[1].sapyb(2, cil_tc[1], 3, out=cil_tc[2])
+
+numpy.testing.assert_allclose(cil_tc[2].as_array().numpy(), 
+                              (res).as_array().numpy())
+
+print(res.array)
+
+print("########## test sapyb with numpy ########## ")
+b = [ numpy.ones([2, 4], dtype=numpy.float32),
+      numpy.ones([2, 4], dtype=numpy.float32),
+      numpy.ones([2, 4], dtype=numpy.float32)]
+
+cil_np = [ DataContainer(el, deep_copy=False) for el in b]
+
+# 2 * 1 + 3 * 1
+res = cil_np[1].sapyb(2, cil_np[1], 3)
+cil_np[1].sapyb(2, cil_np[1], 3, out=cil_np[2])
+
+numpy.testing.assert_allclose(cil_np[2].as_array(), 
+                              (res).as_array())
+
+print(res.array)
\ No newline at end of file

From 1f1fbd06fabe1e9490547b7783088cf8d6684162 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 8 Apr 2025 17:27:00 +0100
Subject: [PATCH 05/25] fix ImageData

Co-authored-by: Gemma Fardell <gfardell@users.noreply.github.com>
---
 Wrappers/Python/cil/framework/image_data.py | 39 ++++++++++-----------
 1 file changed, 19 insertions(+), 20 deletions(-)

diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index c8619a8729..330be43cc0 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -49,22 +49,22 @@ def __init__(self,
                  deep_copy=False,
                  geometry=None,
                  **kwargs):
-        '''Default to NumPy array if array is not passed'''
+        '''Default to np array if array is not passed'''
 
-        # xp = get_namespace()
         
         # TODO: change
         if array is None:
-            # defaults to a numpy array
+            # defaults to a np array
             xp = np
+            dtype = kwargs.get('dtype', xp.float32)
             array = xp.empty(geometry.shape, dtype=dtype)
         elif issubclass(type(array) , DataContainer):
             # this is a reference so we might make a mess modifying both ends
             array = array.as_array()
             xp = array_namespace(array)
-        # elif issubclass(type(array) , numpy.ndarray):
+        # elif issubclass(type(array) , np.ndarray):
             # remove singleton dimensions
-            # array = numpy.squeeze(array)
+            # array = np.squeeze(array)
         else:
             # Consider array as an object is compliant to the array API
             # https://docs.scipy.org/doc/scipy-1.15.2/dev/api-dev/array_api.html 
@@ -74,7 +74,7 @@ def __init__(self,
             array = xp.squeeze(array)
             
         # else:
-        #     raise TypeError('array must be a CIL type DataContainer or numpy.ndarray got {}'.format(type(array)))
+        #     raise TypeError('array must be a CIL type DataContainer or np.ndarray got {}'.format(type(array)))
 
         if array.shape != geometry.shape:
             raise ValueError('Shape mismatch {} {}'.format(array.shape, geometry.shape))
@@ -83,8 +83,7 @@ def __init__(self,
             raise ValueError('Number of dimensions are not 2 or 3 or 4 : {0}'.format(array.ndim))
 
 
-        dtype = kwargs.get('dtype', xp.float32)
-
+        
         if geometry is None:
             raise AttributeError("ImageData requires a geometry")
 
@@ -98,11 +97,11 @@ def __init__(self,
     def __eq__(self, other):
         '''
         Check if two ImageData objects are equal. This is done by checking if the geometry, data and dtype are equal.
-        Also, if the other object is a numpy.ndarray, it will check if the data and dtype are equal.
+        Also, if the other object is a np.ndarray, it will check if the data and dtype are equal.
         
         Parameters
         ----------
-        other: ImageData or numpy.ndarray
+        other: ImageData or np.ndarray
             The object to compare with.
         
         Returns
@@ -112,11 +111,11 @@ def __eq__(self, other):
         '''
 
         if isinstance(other, ImageData):
-            if numpy.array_equal(self.as_array(), other.as_array()) \
+            if np.array_equal(self.as_array(), other.as_array()) \
                 and self.geometry == other.geometry \
                 and self.dtype == other.dtype:
                 return True 
-        elif numpy.array_equal(self.as_array(), other) and self.dtype==other.dtype:
+        elif np.array_equal(self.as_array(), other) and self.dtype==other.dtype:
             return True
         else:
             return False
@@ -137,7 +136,7 @@ def get_slice(self,channel=None, vertical=None, horizontal_x=None, horizontal_y=
         if vertical == 'centre':
             dim = self.geometry.dimension_labels.index('vertical')
             centre_slice_pos = (self.geometry.shape[dim]-1) / 2.
-            ind0 = int(numpy.floor(centre_slice_pos))
+            ind0 = int(np.floor(centre_slice_pos))
 
             w2 = centre_slice_pos - ind0
             out = DataContainer.get_slice(self, channel=channel, vertical=ind0, horizontal_x=horizontal_x, horizontal_y=horizontal_y)
@@ -182,10 +181,10 @@ def apply_circular_mask(self, radius=0.99, in_place=True):
         y_range = (ig.voxel_num_y-1)/2
         x_range = (ig.voxel_num_x-1)/2
 
-        Y, X = numpy.ogrid[-y_range:y_range+1,-x_range:x_range+1]
+        Y, X = np.ogrid[-y_range:y_range+1,-x_range:x_range+1]
 
         # use centre from geometry in units distance to account for aspect ratio of pixels
-        dist_from_center = numpy.sqrt((X*ig.voxel_size_x+ ig.center_x)**2 + (Y*ig.voxel_size_y+ig.center_y)**2)
+        dist_from_center = np.sqrt((X*ig.voxel_size_x+ ig.center_x)**2 + (Y*ig.voxel_size_y+ig.center_y)**2)
 
         size_x = ig.voxel_num_x * ig.voxel_size_x
         size_y = ig.voxel_num_y * ig.voxel_size_y
@@ -197,17 +196,17 @@ def apply_circular_mask(self, radius=0.99, in_place=True):
 
         # approximate the voxel as a circle and get the radius
         # ie voxel area = 1, circle of area=1 has r = 0.56
-        r=((ig.voxel_size_x * ig.voxel_size_y )/numpy.pi)**(1/2)
+        r=((ig.voxel_size_x * ig.voxel_size_y )/np.pi)**(1/2)
 
         # we have the voxel centre distance to mask. voxels with distance greater than |r| are fully inside or outside.
         # values on the border region between -r and r are preserved
         mask =(radius_applied-dist_from_center).clip(-r,r)
 
         #  rescale to -pi/2->+pi/2
-        mask *= (0.5*numpy.pi)/r
+        mask *= (0.5*np.pi)/r
 
         # the sin of the linear distance gives us an approximation of area of the circle to include in the mask
-        numpy.sin(mask, out = mask)
+        np.sin(mask, out = mask)
 
         # rescale the data 0 - 1
         mask = 0.5 + mask * 0.5
@@ -224,13 +223,13 @@ def apply_circular_mask(self, radius=0.99, in_place=True):
 
         if in_place == True:
             self.reorder(labels)
-            numpy.multiply(self.array, mask, out=self.array)
+            np.multiply(self.array, mask, out=self.array)
             self.reorder(labels_orig)
 
         else:
             image_data_out = self.copy()
             image_data_out.reorder(labels)
-            numpy.multiply(image_data_out.array, mask, out=image_data_out.array)
+            np.multiply(image_data_out.array, mask, out=image_data_out.array)
             image_data_out.reorder(labels_orig)
 
             return image_data_out

From e8b7a9237650c66aaf857abdb7e2d80d35863934 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 8 Apr 2025 22:09:45 +0100
Subject: [PATCH 06/25] add array_api_compat to recipe

---
 recipe/meta.yaml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/recipe/meta.yaml b/recipe/meta.yaml
index 3305bcc3ca..60b22a49d9 100644
--- a/recipe/meta.yaml
+++ b/recipe/meta.yaml
@@ -72,6 +72,7 @@ requirements:
     - tqdm
     - numba
     - zenodo_get >=1.6
+    - array-api-compat
 
   #optional packages with version dependancies
   run_constrained:

From 0bba4b0eeb5094ff894e3a975e2b4791a4183419 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 8 Apr 2025 22:33:45 +0100
Subject: [PATCH 07/25] fix GHA

---
 scripts/requirements-test.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/scripts/requirements-test.yml b/scripts/requirements-test.yml
index 9abef82307..5f27111299 100644
--- a/scripts/requirements-test.yml
+++ b/scripts/requirements-test.yml
@@ -49,5 +49,6 @@ dependencies:
   - tqdm
   - zenodo_get >=1.6
   - pip
+  - array-api-compat
   - pip:
     - unittest-parametrize

From 6827b4ab323854e449c76897f0e399f111ed0d3e Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 9 Apr 2025 11:16:41 +0100
Subject: [PATCH 08/25] Add parametrize to unittest

Co-authored-by: Hannah Robarts <hrobarts@users.noreply.github.com>
Co-authored-by: MikeSullivan7 <MikeSullivan7@users.noreply.github.com>
---
 Wrappers/Python/test/test_DataContainer.py | 23 ++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index b1003dddd2..9f029d6650 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -19,6 +19,7 @@
 import logging
 import sys
 
+import array_api_compat.numpy
 import numpy as np
 
 from cil.framework import (DataContainer, ImageGeometry, ImageData, VectorGeometry, AcquisitionData,
@@ -28,6 +29,12 @@
 from testclass import CCPiTestClass
 from utils import initialise_tests
 
+import array_api_compat
+from unittest_parametrize import param, parametrize, ParametrizedTestCase
+import numpy
+import torch
+
+
 log = logging.getLogger(__name__)
 initialise_tests()
 
@@ -37,19 +44,23 @@ def aid(x):
     return x.as_array().__array_interface__['data'][0]
 
 
-class TestDataContainer(CCPiTestClass):
-    def create_DataContainer(self, X,Y,Z, value=1):
+class TestDataContainer(ParametrizedTestCase, CCPiTestClass):
+    def create_DataContainer(self, X,Y,Z, value=1, backend='numpy'):
         a = value * np.ones((X, Y, Z), dtype='float32')
         #print("a refcount " , sys.getrefcount(a))
         ds = DataContainer(a, False, ['X', 'Y', 'Z'])
         return ds
 
-    def test_creation_nocopy(self):
+    @parametrize("xp, raise_error, err_type", 
+        [param(numpy, None, None, id="numpy"), 
+         param(torch, None, None, id="torch"),
+         ]) 
+    def test_creation_nocopy(self, xp, raise_error, err_type):
         shape = 2, 3, 4, 5
-        size = np.prod(shape)
-        a = np.arange(size)
+        size = xp.prod(xp.asarray(shape))
+        a = xp.arange(size)
         #print("a refcount " , sys.getrefcount(a))
-        a = np.reshape(a, shape)
+        a = xp.reshape(a, shape)
         #print("a refcount " , sys.getrefcount(a))
         ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W'])
         #print("a refcount " , sys.getrefcount(a))

From b0b32d86962c1c258660e776052755bca1a152b5 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 9 Apr 2025 22:54:22 +0100
Subject: [PATCH 09/25] added tests adn bugfixes

---
 .../Python/cil/framework/acquisition_data.py  | 20 +++---
 .../Python/cil/framework/data_container.py    | 29 ++++++---
 Wrappers/Python/test/test_DataContainer.py    | 63 +++++++++++++------
 3 files changed, 78 insertions(+), 34 deletions(-)

diff --git a/Wrappers/Python/cil/framework/acquisition_data.py b/Wrappers/Python/cil/framework/acquisition_data.py
index 207f5b947e..bd3cace537 100644
--- a/Wrappers/Python/cil/framework/acquisition_data.py
+++ b/Wrappers/Python/cil/framework/acquisition_data.py
@@ -20,7 +20,8 @@
 from .labels import AcquisitionDimension, Backend
 from .data_container import DataContainer
 from .partitioner import Partitioner
-
+import array_api_compat
+from array_api_compat import array_namespace # https://data-apis.org/array-api-compat/
 
 class AcquisitionData(DataContainer, Partitioner):
     '''DataContainer for holding 2D or 3D sinogram'''
@@ -49,7 +50,7 @@ def __init__(self,
                  geometry = None,
                  **kwargs):
 
-        dtype = kwargs.get('dtype', numpy.float32)
+        
 
         if geometry is None:
             raise AttributeError("AcquisitionData requires a geometry")
@@ -59,15 +60,18 @@ def __init__(self,
                 raise ValueError("Deprecated: 'dimension_labels' cannot be set with 'allocate()'. Use 'geometry.set_labels()' to modify the geometry before using allocate.")
 
         if array is None:
-            array = numpy.empty(geometry.shape, dtype=dtype)
+            xp = array_api_compat.numpy
+            dtype = kwargs.get('dtype', xp.float32)
+            array = xp.empty(geometry.shape, dtype=dtype)
         elif issubclass(type(array) , DataContainer):
             array = array.as_array()
-        elif issubclass(type(array) , numpy.ndarray):
-            # remove singleton dimensions
-            array = numpy.squeeze(array)
         else:
-            raise TypeError('array must be a CIL type DataContainer or numpy.ndarray got {}'.format(type(array)))
-
+            # xp = array_namespace(array)
+            # idx = numpy.where(numpy.asarray(array.shape) == 1)
+            # # remove singleton dimensions
+            # for axis in idx:    
+            #     array = xp.squeeze(array, axis=axis)
+            array = array.squeeze()
         if array.shape != geometry.shape:
             raise ValueError('Shape mismatch got {} expected {}'.format(array.shape, geometry.shape))
 
diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 88642a119c..3952d55552 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -237,15 +237,13 @@ def fill(self, array, **dimension):
 
         dc.fill(some_data, vertical=1, horizontal_x=32)
         will copy the data in some_data into the data container.
+        https://data-apis.org/array-api/latest/design_topics/copies_views_and_mutation.html
+        https://data-apis.org/array-api/latest/API_specification/generated/array_api.array.__setitem__.html#array_api.array.__setitem__
         '''
         if id(array) == id(self.array):
             return
         if dimension == {}:
-            if isinstance(array, numpy.ndarray):
-                numpy.copyto(self.array, array)
-            elif isinstance(array, Number):
-                self.array.fill(array)
-            elif issubclass(array.__class__ , DataContainer):
+            if issubclass(array.__class__ , DataContainer):
 
                 try:
                     if self.dimension_labels != array.dimension_labels:
@@ -260,7 +258,9 @@ def fill(self, array, **dimension):
                                      'Expecting shape {0} got {1}'.format(
                                      self.shape,array.shape))
             else:
-                raise TypeError('Can fill only with number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
+                # raise TypeError('Can fill only with number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
+                xp = array_namespace(self.as_array())
+                self.array.__setitem__(slice(None, None, None), array)
         else:
 
             axis = [':']* self.number_of_dimensions
@@ -276,7 +276,7 @@ def fill(self, array, **dimension):
                     command += ','
                 command += str(el)
                 i+=1
-
+            
             if isinstance(array, numpy.ndarray):
                 command = command + "] = array[:]"
             elif issubclass(array.__class__, DataContainer):
@@ -661,7 +661,20 @@ def pixel_wise_unary(self, pwop, *args,  **kwargs):
         elif issubclass(type(out), DataContainer):
             if self.check_dimensions(out):
                 kwargs['out'] = out.as_array()
-                pwop(self.as_array(), *args, **kwargs )
+                try:
+                    pwop(self.as_array(), *args, **kwargs )
+                except TypeError as te:
+                    import inspect
+                    txt = inspect.getmembers(pwop)
+                    msg = "Error out parameter not supported: "
+                    for el in txt:
+                        if el[0] in ['__name__', '__module__']:
+                            msg += f"{el[0]}: {el[1]} "
+                    import warnings
+                    warnings.warn(msg)
+                    kwargs.pop('out')
+                    out.fill(pwop(self.as_array(), *args, **kwargs ))
+                    return out
             else:
                 raise ValueError(f"Wrong size for data memory: {out.shape} {self.shape}")
         # elif issubclass(type(out), numpy.ndarray):
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 9f029d6650..87e3d71148 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -79,14 +79,22 @@ def testGb_creation_nocopy(self):
         ds1 = ds.clone()
         self.assertNotEqual(aid(ds), aid(ds1))
 
-    def test_ndim(self):
-        x_np = np.arange(0, 60).reshape(3,4,5)
+    @parametrize("xp, raise_error, err_type", 
+        [param(numpy, None, None, id="numpy"), 
+         param(torch, None, None, id="torch"),
+         ]) 
+    def test_ndim(self, xp, raise_error, err_type):
+        x_np = xp.arange(0, 60).reshape(3,4,5)
         x_cil = DataContainer(x_np)
         self.assertEqual(x_np.ndim, x_cil.ndim)
         self.assertEqual(3, x_cil.ndim)
 
-    def test_DataContainer_equal(self):
-        array = np.linspace(-1, 1, 32, dtype=np.float32).reshape(4, 8)
+    @parametrize("xp, raise_error, err_type", 
+        [param(numpy, None, None, id="numpy"), 
+         param(torch, None, None, id="torch"),
+         ]) 
+    def test_DataContainer_equal(self, xp, raise_error, err_type):
+        array = xp.linspace(-1, 1, 32, dtype=xp.float32).reshape(4, 8)
         data = DataContainer(array)
         data1 = data.copy()
 
@@ -99,8 +107,12 @@ def test_DataContainer_equal(self):
         data1 += 1
         self.assertFalse((data ==  data1).all())
 
-    def test_AcquisitionData_equal(self):
-        array = np.linspace(-1, 1, 32, dtype=np.float32).reshape(4, 8)
+    @parametrize("xp, raise_error, err_type", 
+        [param(numpy, None, None, id="numpy"), 
+         param(torch, None, None, id="torch"),
+         ]) 
+    def test_AcquisitionData_equal(self, xp, raise_error, err_type):
+        array = xp.linspace(-1, 1, 32, dtype=xp.float32).reshape(4, 8)
         geom = AcquisitionGeometry.create_Parallel3D().set_panel((8, 4)).set_channels(1).set_angles([1])
         data = AcquisitionData(array, geometry=geom)
 
@@ -123,8 +135,9 @@ def test_AcquisitionData_equal(self):
 
         # Check the equality of two AcquisitionData with different dtypes
         data_different_dtype = data.copy()
-        data_different_dtype.array = data_different_dtype.array.astype(np.float64)
-        self.assertFalse(data == data_different_dtype)
+        # data_different_dtype.array = data_different_dtype.array.astype(np.float64)
+        arr = xp.asarray(data_different_dtype.array, dtype=xp.float64)
+        self.assertFalse(data == arr)
 
 
         # Check the equality of two AcquisitionData with different labels
@@ -133,8 +146,12 @@ def test_AcquisitionData_equal(self):
         data_different_labels.geometry.set_labels([AcquisitionDimension("ANGLE"), AcquisitionDimension("CHANNEL") ])
         self.assertFalse(data == data_different_labels)
 
-    def test_ImageData_equal(self):
-        array = np.linspace(-1, 1, 32, dtype=np.float32).reshape(4, 8)
+    @parametrize("xp, device, raise_error, err_type", 
+        [param(numpy, None, None, None, id="numpy"), 
+         param(torch, None, None, None, id="torch_cpu"),
+         ]) 
+    def test_ImageData_equal(self, xp, device, raise_error, err_type):
+        array = xp.linspace(-1, 1, 32, dtype=xp.float32).reshape(4, 8)
         geom = ImageGeometry(voxel_num_x=8, voxel_num_y=4)
         data = ImageData(array, geometry=geom)
 
@@ -166,8 +183,11 @@ def test_ImageData_equal(self):
         data_different_labels.geometry.set_labels([ImageDimension("VERTICAL"), ImageDimension("HORIZONTAL_X")])
         self.assertFalse(data == data_different_labels)
 
-
-    def testInlineAlgebra(self):
+    @parametrize("xp, device, raise_error, err_type", 
+        [param(numpy, None, None, None, id="numpy"), 
+         param(torch, None, None, None, id="torch_cpu"),
+         ]) 
+    def testInlineAlgebra(self, xp, device, raise_error, err_type):
         X, Y, Z = 8, 16, 32
         a = np.ones((X, Y, Z), dtype='float32')
         b = np.ones((X, Y, Z), dtype='float32')
@@ -202,10 +222,14 @@ def testInlineAlgebra(self):
         np.testing.assert_array_almost_equal(ds.as_array(), b)
         # self.assertEqual(ds.as_array()[0][0][0], 1.)
 
-    def test_unary_operations(self):
+    @parametrize("xp, device, raise_error, err_type", 
+        [param(numpy, None, None, None, id="numpy"), 
+         param(torch, None, None, None, id="torch_cpu"),
+         ])
+    def test_unary_operations(self, xp, device, raise_error, err_type):
         X, Y, Z = 8, 16, 32
-        a = -np.ones((X, Y, Z), dtype='float32')
-        b = np.ones((X, Y, Z), dtype='float32')
+        a = -xp.ones((X, Y, Z), dtype=xp.float32)
+        b = xp.ones((X, Y, Z), dtype=xp.float32)
         ds = DataContainer(a, False, ['X', 'Y', 'Z'])
 
         ds.sign(out=ds)
@@ -1319,11 +1343,14 @@ def test_negation(self):
 
         np.testing.assert_array_equal(ds.as_array(), -a)
 
-
-    def test_fill_dimension_ImageData(self):
+    @parametrize("xp, device, raise_error, err_type", 
+        [param(numpy, None, None, None, id="numpy"), 
+         param(torch, None, None, None, id="torch_cpu"),
+         ])
+    def test_fill_dimension_ImageData(self, xp, device, raise_error, err_type):
         ig = ImageGeometry(2,3,4)
         u = ig.allocate(0)
-        a = np.ones((4,2))
+        a = xp.ones((4,2))
         # default_labels = [ImageDimension["VERTICAL"], ImageDimension["HORIZONTAL_Y"], ImageDimension["HORIZONTAL_X"]]
 
         data = u.as_array()

From 4d54f3b88e9803985cdf7bc47f6c123c32f7092c Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Fri, 11 Apr 2025 22:28:08 +0100
Subject: [PATCH 10/25] work in progress

---
 .../Python/cil/framework/data_container.py    | 55 ++++++++++---------
 1 file changed, 29 insertions(+), 26 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 3952d55552..bd12e9d0de 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -174,14 +174,19 @@ def get_slice(self, **kw):
         dimension_labels_list = list(self.dimension_labels)
 
         #remove axes from array and labels
+        squeeze_axes = []
         for key, value in kw.items():
             if value is not None:
                 axis = dimension_labels_list.index(key)
+                squeeze_axes.append(axis)
                 dimension_labels_list.remove(key)
                 if new_array is None:
                     new_array = self.as_array()
                 new_array = new_array.take(indices=value, axis=axis)
 
+        print ("new_array.shape", new_array.shape)
+        # xp = array_namespace(new_array)
+        # new_array = xp.squeeze(new_array, axis=tuple(squeeze_axes))
         if new_array.ndim > 1:
             return DataContainer(new_array, False, dimension_labels_list, suppress_warning=True)
         from .vector_data import VectorData
@@ -252,7 +257,9 @@ def fill(self, array, **dimension):
                     pass
 
                 if self.array.shape == array.shape:
-                    numpy.copyto(self.array, array.array)
+                    # numpy.copyto(self.array, array.array)
+                    xp = array_namespace(self.as_array())
+                    self.array.__setitem__(slice(None, None, None), array.array)
                 else:
                     raise ValueError('Cannot fill with the provided array.' + \
                                      'Expecting shape {0} got {1}'.format(
@@ -262,31 +269,27 @@ def fill(self, array, **dimension):
                 xp = array_namespace(self.as_array())
                 self.array.__setitem__(slice(None, None, None), array)
         else:
-
-            axis = [':']* self.number_of_dimensions
-            dimension_labels = tuple(self.dimension_labels)
-            for k,v in dimension.items():
-                i = dimension_labels.index(k)
-                axis[i] = v
-
-            command = 'self.array['
-            i = 0
-            for el in axis:
-                if i > 0:
-                    command += ','
-                command += str(el)
-                i+=1
-            
-            if isinstance(array, numpy.ndarray):
-                command = command + "] = array[:]"
-            elif issubclass(array.__class__, DataContainer):
-                command = command + "] = array.as_array()[:]"
-            elif isinstance (array, Number):
-                command = command + "] = array"
-            else:
-                raise TypeError('Can fill only with number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
-            exec(command)
-
+            slices = []
+            where = 0
+            for i,el in enumerate(self.dimension_labels):
+                for k,v in dimension.items():
+                    if el == k:
+                        slices.append(slice(v,v+1,None))
+                        where = i
+                    else:
+                        slices.append(slice(None, None, None))
+            # give new shape to the array to insert so that it fits
+            # the shape of the target data
+            try:
+                 
+                xp = array_namespace(array)
+                xp.expand_dims(array, axis=where)
+                self.array.__setitem__(tuple(slices), array)
+                xp.squeeze(array, axis=where)
+
+            except AttributeError as ae:
+                self.array.__setitem__(tuple(slices), array)
+                
 
     def check_dimensions(self, other):
         return self.shape == other.shape

From 5db9a1fe5a8871ca61f4d259029bac4a312141c3 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Mon, 2 Jun 2025 11:34:14 +0100
Subject: [PATCH 11/25] modified test

---
 Wrappers/Python/test/torch_test.py | 20 +++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/Wrappers/Python/test/torch_test.py b/Wrappers/Python/test/torch_test.py
index 828e31bad6..622bd98ce5 100644
--- a/Wrappers/Python/test/torch_test.py
+++ b/Wrappers/Python/test/torch_test.py
@@ -103,4 +103,22 @@
 numpy.testing.assert_allclose(cil_np[2].as_array(), 
                               (res).as_array())
 
-print(res.array)
\ No newline at end of file
+print(res.array)
+
+torch.squeeze(a[0])
+
+shape = [3,5,1,5]
+
+print (numpy.where(numpy.asarray(shape) == 1))
+
+tt = torch.asarray([1,2,-1,10])
+print(torch.sign(tt))
+ts = tt*0
+torch.sign(tt, out=ts)
+print (ts)
+ts = tt*0
+xp = array_namespace(tt)
+# xp.sign(tt, out=ts)
+print (ts)
+tt.sign(out=ts)
+print (ts)
\ No newline at end of file

From 101175ebdfda50b085c6d1c73771f00d173eb623 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Mon, 2 Jun 2025 13:48:24 +0100
Subject: [PATCH 12/25] convert numpy ComplexWarning to Exception

slightly relax tolerance of tests
---
 Wrappers/Python/cil/framework/data_container.py | 3 +++
 Wrappers/Python/test/test_DataContainer.py      | 6 ++++--
 Wrappers/Python/test/test_functions.py          | 9 ++++++---
 3 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index bd12e9d0de..8fe193547d 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -266,8 +266,11 @@ def fill(self, array, **dimension):
                                      self.shape,array.shape))
             else:
                 # raise TypeError('Can fill only with number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
+                import warnings
+                warnings.filterwarnings('error')
                 xp = array_namespace(self.as_array())
                 self.array.__setitem__(slice(None, None, None), array)
+                warnings.resetwarnings()
         else:
             slices = []
             where = 0
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 87e3d71148..7caffa4089 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -719,7 +719,9 @@ def test_VectorGeometry_allocate_dtype(self):
 
     def complex_allocate_geometry_test(self, geometry):
         data = geometry.allocate(dtype=np.complex64)
-        r = (1 + 1j*1)* np.ones(data.shape, dtype=data.dtype)
+        from array_api_compat import array_namespace
+        xp = array_namespace(data.array)
+        r = (1 + 1j*1)* xp.ones(data.shape, dtype=data.dtype)
         data.fill(r)
         self.assertAlmostEqual(data.squared_norm(), data.size * 2)
         np.testing.assert_almost_equal(data.abs().array, np.abs(r))
@@ -728,7 +730,7 @@ def complex_allocate_geometry_test(self, geometry):
         try:
             data1.fill(r)
             self.assertTrue(False)
-        except TypeError as err:
+        except numpy.exceptions.ComplexWarning as err:
             log.info(str(err))
             self.assertTrue(True)
 
diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py
index 8fa9c1558a..d6c96db8ad 100644
--- a/Wrappers/Python/test/test_functions.py
+++ b/Wrappers/Python/test/test_functions.py
@@ -1139,7 +1139,8 @@ def soft_shrinkage_test(self, x):
         tau = -1.
         with self.assertWarns(UserWarning):
             ret = soft_shrinkage(-0.5 *x, tau)
-            np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()))
+            np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()), 
+                                       atol=3e-7, rtol=2e-7)
         tau = -3j
         with self.assertRaises(ValueError):
             ret = soft_shrinkage(-0.5 *x, tau)
@@ -1162,7 +1163,8 @@ def soft_shrinkage_test(self, x):
         np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()))
         tau = -1.* np.ones_like(x.as_array())
         ret = soft_shrinkage(-0.5 *x, tau)
-        np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()))
+        np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()),
+                                   rtol=3e-7, atol=2e-7)
 
         # tau DataContainer
         tau = 1. * x
@@ -1182,7 +1184,8 @@ def soft_shrinkage_test(self, x):
         np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()))
         tau = -1. * x
         ret = soft_shrinkage(-0.5 *x, tau)
-        np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()))
+        np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()),
+                                   rtol=3e-7, atol=2e-7)
 
         np.testing.assert_allclose(ret.as_array().imag, np.zeros_like(ret.as_array().imag), atol=1e-6, rtol=1e-6)
 

From c53115ba622d52c8ec0c2842532b298cfa6c82c2 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 3 Jun 2025 09:27:14 +0100
Subject: [PATCH 13/25] add function to squeeze an array on multiple dimensions
 with array_api_compat

and fix numpy specific unit test
---
 Wrappers/Python/cil/framework/data_container.py | 17 +++++++++++++++++
 Wrappers/Python/cil/framework/image_data.py     |  7 ++-----
 Wrappers/Python/test/test_DataContainer.py      |  6 +-----
 3 files changed, 20 insertions(+), 10 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 8fe193547d..a535cc711f 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -27,6 +27,23 @@
 from .cilacc import cilacc
 from cil.utilities.multiprocessing import NUM_THREADS
 
+def squeeze_array(array):
+    '''squeezes the array, removing all singleton dimensions recursively
+    
+    Parameters
+    ----------
+    array : array-like
+        The array to squeeze
+    '''
+    xp = array_namespace(array)
+    # find and remove singleton dimensions
+    s = xp.nonzero(xp.asarray(array.shape) == 1)[0]
+    singletons = s.tolist()
+    if len(singletons) > 1:
+        # remove singleton dimension recursively
+        return squeeze_array(xp.squeeze(array, axis=singletons[0]))
+    return array
+
 
 class DataContainer(object):
     '''Generic class to hold data
diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index 330be43cc0..79a6a91109 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -17,7 +17,7 @@
 # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
 import numpy as np
 
-from .data_container import DataContainer
+from .data_container import DataContainer, squeeze_array
 from .labels import ImageDimension, Backend
 
 from array_api_compat import array_namespace # https://data-apis.org/array-api-compat/
@@ -69,10 +69,7 @@ def __init__(self,
             # Consider array as an object is compliant to the array API
             # https://docs.scipy.org/doc/scipy-1.15.2/dev/api-dev/array_api.html 
             # this might raise an exception but that's fine
-            xp = array_namespace(array)
-            # remove singleton dimensions
-            array = xp.squeeze(array)
-            
+            array = squeeze_array(array)
         # else:
         #     raise TypeError('array must be a CIL type DataContainer or np.ndarray got {}'.format(type(array)))
 
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 7caffa4089..d509d70860 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -173,8 +173,7 @@ def test_ImageData_equal(self, xp, device, raise_error, err_type):
         self.assertFalse(data == data_different_shape)
 
         # Check the equality of two ImageData with different dtypes
-        data_different_dtype = data.copy()
-        data_different_dtype.array = data_different_dtype.array.astype(np.float64)
+        data_different_dtype = data.geometry.allocate(0, dtype=np.float64)
         self.assertFalse(data == data_different_dtype)
 
 
@@ -1355,9 +1354,6 @@ def test_fill_dimension_ImageData(self, xp, device, raise_error, err_type):
         a = xp.ones((4,2))
         # default_labels = [ImageDimension["VERTICAL"], ImageDimension["HORIZONTAL_Y"], ImageDimension["HORIZONTAL_X"]]
 
-        data = u.as_array()
-        axis_number = u.get_dimension_axis('horizontal_y')
-
         u.fill(a, horizontal_y=0)
         np.testing.assert_array_equal(u.get_slice(horizontal_y=0).as_array(), a)
 

From 8283b878285ca9db4d2bb439b1cc6d1023900b60 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 14:33:01 +0100
Subject: [PATCH 14/25] added cil array_api_ compat

---
 .../Python/cil/framework/array_api_compat.py  | 64 +++++++++++++++++++
 .../Python/cil/framework/data_container.py    | 61 +++++-------------
 Wrappers/Python/cil/framework/image_data.py   |  7 +-
 3 files changed, 85 insertions(+), 47 deletions(-)
 create mode 100644 Wrappers/Python/cil/framework/array_api_compat.py

diff --git a/Wrappers/Python/cil/framework/array_api_compat.py b/Wrappers/Python/cil/framework/array_api_compat.py
new file mode 100644
index 0000000000..02d1a97de6
--- /dev/null
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -0,0 +1,64 @@
+# Updates to array API compatibility layer for CIL
+from array_api_compat import array_namespace
+
+def expand_dims(array, axis):
+    '''Expand dimensions of an array along specified axes.
+    
+    Parameters
+    ----------
+    array : array-like
+        The input array to expand.
+    axis : int or tuple of int
+        The axis or axes along which to expand the dimensions.
+        
+    Returns
+    -------
+    array-like
+        The array with expanded dimensions. It may be a new array or the same array with expanded dimensions.
+
+    Raises:
+    --------
+    IndexError If provided an invalid axis position, an IndexError should be raised.
+    
+    Notes:
+    This function recursively expands the dimensions of the input array along the specified axes if a list or tuple of ints is provided.
+    '''
+    xp = array_namespace(array)
+    
+    if isinstance(axis, int):
+        return xp.expand_dims(array, axis=axis)
+    axis = list(axis)
+    ax = axis.pop(0)
+    if len(axis) == 1:
+        axis = axis[0]
+    return expand_dims(xp.expand_dims(array, axis=ax), axis=axis)
+
+def squeeze(array, axis=None):
+    '''squeezes the array, removing all singleton dimensions recursively
+    
+    Parameters
+    ----------
+    array : array-like
+        The array to squeeze
+    axis : int or tuple of int, optional
+        The axis or axes to squeeze. If None, all singleton dimensions are removed.
+
+    Returns
+    -------
+    array-like
+        The squeezed array with all singleton dimensions removed. If the input array has no singleton dimensions, it is returned unchanged.
+    '''
+    xp = array_namespace(array)
+    # find and remove singleton dimensions
+    if axis is None:
+        s = xp.nonzero(xp.asarray(array.shape) == 1)[0]
+        axis = s.tolist()
+    if isinstance(axis, int):
+        return xp.squeeze(array, axis=axis)
+    # process from the largest axis to the smallest
+    axis = list(axis)
+    axis.sort(reverse=True)
+    ax = axis.pop(0)
+    if len(axis) == 1:
+        axis = axis[0]
+    return squeeze(xp.squeeze(array, axis=ax), axis=axis)
\ No newline at end of file
diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index a535cc711f..69e1956b55 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -26,23 +26,9 @@
 
 from .cilacc import cilacc
 from cil.utilities.multiprocessing import NUM_THREADS
+from .array_api_compat import squeeze as cil_squeeze
+from .array_api_compat import expand_dims as cil_expand_dims
 
-def squeeze_array(array):
-    '''squeezes the array, removing all singleton dimensions recursively
-    
-    Parameters
-    ----------
-    array : array-like
-        The array to squeeze
-    '''
-    xp = array_namespace(array)
-    # find and remove singleton dimensions
-    s = xp.nonzero(xp.asarray(array.shape) == 1)[0]
-    singletons = s.tolist()
-    if len(singletons) > 1:
-        # remove singleton dimension recursively
-        return squeeze_array(xp.squeeze(array, axis=singletons[0]))
-    return array
 
 
 class DataContainer(object):
@@ -264,52 +250,39 @@ def fill(self, array, **dimension):
         '''
         if id(array) == id(self.array):
             return
+        if issubclass(array.__class__ , DataContainer):
+            return self.fill(array.as_array(), **dimension)
         if dimension == {}:
-            if issubclass(array.__class__ , DataContainer):
-
-                try:
-                    if self.dimension_labels != array.dimension_labels:
-                        raise ValueError('Input array is not in the same order as destination array. Use "array.reorder()"')
-                except AttributeError:
-                    pass
-
-                if self.array.shape == array.shape:
-                    # numpy.copyto(self.array, array.array)
-                    xp = array_namespace(self.as_array())
-                    self.array.__setitem__(slice(None, None, None), array.array)
-                else:
-                    raise ValueError('Cannot fill with the provided array.' + \
-                                     'Expecting shape {0} got {1}'.format(
-                                     self.shape,array.shape))
-            else:
-                # raise TypeError('Can fill only with number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
-                import warnings
-                warnings.filterwarnings('error')
-                xp = array_namespace(self.as_array())
-                self.array.__setitem__(slice(None, None, None), array)
-                warnings.resetwarnings()
+            # raise TypeError('Can fill only with number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
+            import warnings
+            warnings.filterwarnings('error')
+            xp = array_namespace(self.as_array())
+            self.array.__setitem__(slice(None, None, None), array)
+            warnings.resetwarnings()
         else:
             slices = []
-            where = 0
+            where = [] * self.number_of_dimensions
             for i,el in enumerate(self.dimension_labels):
                 for k,v in dimension.items():
                     if el == k:
                         slices.append(slice(v,v+1,None))
-                        where = i
+                        where[i] = 1
                     else:
                         slices.append(slice(None, None, None))
             # give new shape to the array to insert so that it fits
             # the shape of the target data
             try:
                  
-                xp = array_namespace(array)
-                xp.expand_dims(array, axis=where)
+                array = cil_expand_dims(array, axis=where)
                 self.array.__setitem__(tuple(slices), array)
-                xp.squeeze(array, axis=where)
+                array = cil_squeeze(array, axis=where)
 
             except AttributeError as ae:
                 self.array.__setitem__(tuple(slices), array)
                 
+            except TypeError as ae:
+                self.array.__setitem__(tuple(slices), array)
+                
 
     def check_dimensions(self, other):
         return self.shape == other.shape
diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index 79a6a91109..ba59b5d1f5 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -17,7 +17,8 @@
 # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
 import numpy as np
 
-from .data_container import DataContainer, squeeze_array
+from .data_container import DataContainer
+from .array_api_compat import squeeze as cil_squeeze 
 from .labels import ImageDimension, Backend
 
 from array_api_compat import array_namespace # https://data-apis.org/array-api-compat/
@@ -61,7 +62,7 @@ def __init__(self,
         elif issubclass(type(array) , DataContainer):
             # this is a reference so we might make a mess modifying both ends
             array = array.as_array()
-            xp = array_namespace(array)
+            array = cil_squeeze(array)
         # elif issubclass(type(array) , np.ndarray):
             # remove singleton dimensions
             # array = np.squeeze(array)
@@ -69,7 +70,7 @@ def __init__(self,
             # Consider array as an object is compliant to the array API
             # https://docs.scipy.org/doc/scipy-1.15.2/dev/api-dev/array_api.html 
             # this might raise an exception but that's fine
-            array = squeeze_array(array)
+            array = cil_squeeze(array)
         # else:
         #     raise TypeError('array must be a CIL type DataContainer or np.ndarray got {}'.format(type(array)))
 

From b695c5bb55325a5d05d407fc85591b7c040a7206 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 15:27:34 +0100
Subject: [PATCH 15/25] fix test_fill_dimension_AcquisitionData

---
 Wrappers/Python/cil/framework/data_container.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 69e1956b55..aff8bccb4f 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -260,15 +260,14 @@ def fill(self, array, **dimension):
             self.array.__setitem__(slice(None, None, None), array)
             warnings.resetwarnings()
         else:
-            slices = []
-            where = [] * self.number_of_dimensions
+            slices = [slice(None, None, None)] * self.number_of_dimensions
+            where = []
             for i,el in enumerate(self.dimension_labels):
                 for k,v in dimension.items():
                     if el == k:
-                        slices.append(slice(v,v+1,None))
-                        where[i] = 1
-                    else:
-                        slices.append(slice(None, None, None))
+                        slices[i] = slice(v,v+1,None)
+                        where.append(i)
+
             # give new shape to the array to insert so that it fits
             # the shape of the target data
             try:

From a6c6e569c90c10f9530c01b1d6b41993514af94f Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 15:46:44 +0100
Subject: [PATCH 16/25] fixed DataContainer unittests

---
 Wrappers/Python/cil/framework/array_api_compat.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/Wrappers/Python/cil/framework/array_api_compat.py b/Wrappers/Python/cil/framework/array_api_compat.py
index 02d1a97de6..ab22d93fc5 100644
--- a/Wrappers/Python/cil/framework/array_api_compat.py
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -23,10 +23,13 @@ def expand_dims(array, axis):
     Notes:
     This function recursively expands the dimensions of the input array along the specified axes if a list or tuple of ints is provided.
     '''
+    print("axis", axis)
     xp = array_namespace(array)
     
     if isinstance(axis, int):
         return xp.expand_dims(array, axis=axis)
+    if len(axis) == 1:
+        return xp.expand_dims(array, axis=axis[0])
     axis = list(axis)
     ax = axis.pop(0)
     if len(axis) == 1:
@@ -53,8 +56,15 @@ def squeeze(array, axis=None):
     if axis is None:
         s = xp.nonzero(xp.asarray(array.shape) == 1)[0]
         axis = s.tolist()
+        if len(axis) == 1:
+            axis = axis[0]
+        elif len(axis) == 0:
+            return array
     if isinstance(axis, int):
         return xp.squeeze(array, axis=axis)
+    if len(axis) == 1:
+        return xp.squeeze(array, axis=axis[0])
+    
     # process from the largest axis to the smallest
     axis = list(axis)
     axis.sort(reverse=True)

From c9977e90ad460382a3a52b1c191c23b8753b20e8 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 16:12:32 +0100
Subject: [PATCH 17/25] removed prints

---
 Wrappers/Python/cil/framework/array_api_compat.py | 2 +-
 Wrappers/Python/cil/framework/data_container.py   | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/Wrappers/Python/cil/framework/array_api_compat.py b/Wrappers/Python/cil/framework/array_api_compat.py
index ab22d93fc5..d27261e141 100644
--- a/Wrappers/Python/cil/framework/array_api_compat.py
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -23,7 +23,6 @@ def expand_dims(array, axis):
     Notes:
     This function recursively expands the dimensions of the input array along the specified axes if a list or tuple of ints is provided.
     '''
-    print("axis", axis)
     xp = array_namespace(array)
     
     if isinstance(axis, int):
@@ -59,6 +58,7 @@ def squeeze(array, axis=None):
         if len(axis) == 1:
             axis = axis[0]
         elif len(axis) == 0:
+            # nothing to do
             return array
     if isinstance(axis, int):
         return xp.squeeze(array, axis=axis)
diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 0cb34853dd..f843331f33 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -185,7 +185,6 @@ def get_slice(self, **kw):
                     new_array = self.as_array()
                 new_array = new_array.take(indices=value, axis=axis)
 
-        print ("new_array.shape", new_array.shape)
         # xp = array_namespace(new_array)
         # new_array = xp.squeeze(new_array, axis=tuple(squeeze_axes))
         if new_array.ndim > 1:

From cbe1f3b8608a4cf0eb9cd120cf671e2ff0cdf2ce Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 16:28:30 +0100
Subject: [PATCH 18/25] temporarily disable some tests if torch is not
 installed

---
 Wrappers/Python/test/test_DataContainer.py | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 1f20a51b10..a167eb6c53 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -32,8 +32,11 @@
 import array_api_compat
 from unittest_parametrize import param, parametrize, ParametrizedTestCase
 import numpy
-import torch
-
+try:
+    import torch
+except ImportError:
+    torch = None
+import unittest
 
 log = logging.getLogger(__name__)
 initialise_tests()
@@ -50,7 +53,8 @@ def create_DataContainer(self, X,Y,Z, value=1, backend='numpy'):
         #print("a refcount " , sys.getrefcount(a))
         ds = DataContainer(a, False, ['X', 'Y', 'Z'])
         return ds
-
+    
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
@@ -79,6 +83,7 @@ def testGb_creation_nocopy(self):
         ds1 = ds.clone()
         self.assertNotEqual(aid(ds), aid(ds1))
 
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
@@ -89,6 +94,7 @@ def test_ndim(self, xp, raise_error, err_type):
         self.assertEqual(x_np.ndim, x_cil.ndim)
         self.assertEqual(3, x_cil.ndim)
 
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
@@ -107,6 +113,7 @@ def test_DataContainer_equal(self, xp, raise_error, err_type):
         data1 += 1
         self.assertFalse((data ==  data1).all())
 
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
@@ -146,6 +153,7 @@ def test_AcquisitionData_equal(self, xp, raise_error, err_type):
         data_different_labels.geometry.set_labels([AcquisitionDimension("ANGLE"), AcquisitionDimension("CHANNEL") ])
         self.assertFalse(data == data_different_labels)
 
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
@@ -182,6 +190,7 @@ def test_ImageData_equal(self, xp, device, raise_error, err_type):
         data_different_labels.geometry.set_labels([ImageDimension("VERTICAL"), ImageDimension("HORIZONTAL_X")])
         self.assertFalse(data == data_different_labels)
 
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
@@ -221,6 +230,7 @@ def testInlineAlgebra(self, xp, device, raise_error, err_type):
         np.testing.assert_array_almost_equal(ds.as_array(), b)
         # self.assertEqual(ds.as_array()[0][0][0], 1.)
 
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
@@ -1344,6 +1354,7 @@ def test_negation(self):
 
         np.testing.assert_array_equal(ds.as_array(), -a)
 
+    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),

From 1a8465920ee5629bac2734d834180f57a4305ce5 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 16:33:29 +0100
Subject: [PATCH 19/25] skip tests with torch but executes with numpy

---
 Wrappers/Python/test/test_DataContainer.py | 24 ++++++++++++++--------
 1 file changed, 16 insertions(+), 8 deletions(-)

diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index a167eb6c53..7fb0ecb5b4 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -54,12 +54,13 @@ def create_DataContainer(self, X,Y,Z, value=1, backend='numpy'):
         ds = DataContainer(a, False, ['X', 'Y', 'Z'])
         return ds
     
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
          ]) 
     def test_creation_nocopy(self, xp, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         shape = 2, 3, 4, 5
         size = xp.prod(xp.asarray(shape))
         a = xp.arange(size)
@@ -83,23 +84,25 @@ def testGb_creation_nocopy(self):
         ds1 = ds.clone()
         self.assertNotEqual(aid(ds), aid(ds1))
 
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
          ]) 
     def test_ndim(self, xp, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         x_np = xp.arange(0, 60).reshape(3,4,5)
         x_cil = DataContainer(x_np)
         self.assertEqual(x_np.ndim, x_cil.ndim)
         self.assertEqual(3, x_cil.ndim)
 
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
          ]) 
     def test_DataContainer_equal(self, xp, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         array = xp.linspace(-1, 1, 32, dtype=xp.float32).reshape(4, 8)
         data = DataContainer(array)
         data1 = data.copy()
@@ -113,12 +116,13 @@ def test_DataContainer_equal(self, xp, raise_error, err_type):
         data1 += 1
         self.assertFalse((data ==  data1).all())
 
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
          ]) 
     def test_AcquisitionData_equal(self, xp, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         array = xp.linspace(-1, 1, 32, dtype=xp.float32).reshape(4, 8)
         geom = AcquisitionGeometry.create_Parallel3D().set_panel((8, 4)).set_channels(1).set_angles([1])
         data = AcquisitionData(array, geometry=geom)
@@ -153,12 +157,13 @@ def test_AcquisitionData_equal(self, xp, raise_error, err_type):
         data_different_labels.geometry.set_labels([AcquisitionDimension("ANGLE"), AcquisitionDimension("CHANNEL") ])
         self.assertFalse(data == data_different_labels)
 
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
          ]) 
     def test_ImageData_equal(self, xp, device, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         array = xp.linspace(-1, 1, 32, dtype=xp.float32).reshape(4, 8)
         geom = ImageGeometry(voxel_num_x=8, voxel_num_y=4)
         data = ImageData(array, geometry=geom)
@@ -190,12 +195,13 @@ def test_ImageData_equal(self, xp, device, raise_error, err_type):
         data_different_labels.geometry.set_labels([ImageDimension("VERTICAL"), ImageDimension("HORIZONTAL_X")])
         self.assertFalse(data == data_different_labels)
 
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
          ]) 
     def testInlineAlgebra(self, xp, device, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         X, Y, Z = 8, 16, 32
         a = np.ones((X, Y, Z), dtype='float32')
         b = np.ones((X, Y, Z), dtype='float32')
@@ -230,12 +236,13 @@ def testInlineAlgebra(self, xp, device, raise_error, err_type):
         np.testing.assert_array_almost_equal(ds.as_array(), b)
         # self.assertEqual(ds.as_array()[0][0][0], 1.)
 
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
          ])
     def test_unary_operations(self, xp, device, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         X, Y, Z = 8, 16, 32
         a = -xp.ones((X, Y, Z), dtype=xp.float32)
         b = xp.ones((X, Y, Z), dtype=xp.float32)
@@ -1354,12 +1361,13 @@ def test_negation(self):
 
         np.testing.assert_array_equal(ds.as_array(), -a)
 
-    @unittest.skipIf(torch is None, "torch not available")
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
          ])
     def test_fill_dimension_ImageData(self, xp, device, raise_error, err_type):
+        if xp is None:
+            self.skipTest("torch not available")
         ig = ImageGeometry(2,3,4)
         u = ig.allocate(0)
         a = xp.ones((4,2))

From b4b3cda4053fd37ed1bf92fae1285a8da02d45a1 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 16:43:29 +0100
Subject: [PATCH 20/25] prevent AttributeError: module 'numpy' has no attribute
 'exceptions'

---
 Wrappers/Python/test/test_DataContainer.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 7fb0ecb5b4..5cdee426dd 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -746,7 +746,8 @@ def complex_allocate_geometry_test(self, geometry):
         try:
             data1.fill(r)
             self.assertTrue(False)
-        except numpy.exceptions.ComplexWarning as err:
+        # except numpy.exceptions.ComplexWarning as err:
+        except Exception as err:
             log.info(str(err))
             self.assertTrue(True)
 

From 97fbe3557c27867515acc9f68402d6a2a2daa445 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Wed, 4 Jun 2025 16:56:53 +0100
Subject: [PATCH 21/25] update the type of the data according to the parameter
 requested

---
 Wrappers/Python/cil/framework/acquisition_geometry.py | 5 +++--
 Wrappers/Python/cil/framework/image_geometry.py       | 5 +++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/Wrappers/Python/cil/framework/acquisition_geometry.py b/Wrappers/Python/cil/framework/acquisition_geometry.py
index c257e58517..2a1c04abd4 100644
--- a/Wrappers/Python/cil/framework/acquisition_geometry.py
+++ b/Wrappers/Python/cil/framework/acquisition_geometry.py
@@ -2173,8 +2173,9 @@ def allocate(self, value=0, **kwargs):
         :type dtype: numpy type, default numpy.float32
         '''
         dtype = kwargs.get('dtype', self.dtype)
-
-        out = AcquisitionData(geometry=self.copy(),
+        ng = self.copy()
+        ng.dtype = dtype
+        out = AcquisitionData(geometry=ng,
                               dtype=dtype,
                               suppress_warning=True)
 
diff --git a/Wrappers/Python/cil/framework/image_geometry.py b/Wrappers/Python/cil/framework/image_geometry.py
index 6ffc03ed66..5bf54d3e27 100644
--- a/Wrappers/Python/cil/framework/image_geometry.py
+++ b/Wrappers/Python/cil/framework/image_geometry.py
@@ -262,8 +262,9 @@ def allocate(self, value=0, **kwargs):
 
         if kwargs.get('dimension_labels', None) is not None:
             raise ValueError("Deprecated: 'dimension_labels' cannot be set with 'allocate()'. Use 'geometry.set_labels()' to modify the geometry before using allocate.")
-
-        out = ImageData(geometry=self.copy(),
+        ng = self.copy()
+        ng.dtype = dtype
+        out = ImageData(geometry=ng,
                             dtype=dtype,
                             suppress_warning=True)
 

From fdbdb3aee9b1d77cea01a6163d2c5d2d0aae10a6 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Mon, 9 Jun 2025 15:52:55 +0000
Subject: [PATCH 22/25] WIP added allclose

---
 .../Python/cil/framework/array_api_compat.py  | 23 +++++++++++++++-
 Wrappers/Python/cil/framework/image_data.py   | 20 ++++++++------
 Wrappers/Python/test/test_DataContainer.py    | 27 +++++++++++++++----
 3 files changed, 56 insertions(+), 14 deletions(-)

diff --git a/Wrappers/Python/cil/framework/array_api_compat.py b/Wrappers/Python/cil/framework/array_api_compat.py
index d27261e141..77c77e0887 100644
--- a/Wrappers/Python/cil/framework/array_api_compat.py
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -71,4 +71,25 @@ def squeeze(array, axis=None):
     ax = axis.pop(0)
     if len(axis) == 1:
         axis = axis[0]
-    return squeeze(xp.squeeze(array, axis=ax), axis=axis)
\ No newline at end of file
+    return squeeze(xp.squeeze(array, axis=ax), axis=axis)
+
+import pysnooper
+@pysnooper.snoop()
+def allclose(a, b, rtol=1e-5, atol=1e-6):
+    """
+    Check if two arrays are element-wise equal within a tolerance.allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False)[source]
+Returns True if two arrays are element-wise equal within a tolerance.
+
+The tolerance values are positive, typically very small numbers. 
+The relative difference (rtol * abs(b)) and the absolute difference atol are added together to compare against the absolute difference between a and b.
+    """
+    xp = array_namespace(a.as_array())
+    if array_namespace(b.as_array()) != xp:
+        raise TypeError('Can only compare arrays ' \
+        'with same namespace. Got {} and {}'.format(array_namespace(a), array_namespace(b)))
+    
+    diff = rtol * xp.abs(b) + atol
+    if xp.any(diff < xp.abs(a - b)):
+        print(f"Max difference: {diff.max()}")
+        return False
+    return True
diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index 738af7072f..bc687c9901 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -113,8 +113,9 @@ def __init__(self,
     def __eq__(self, other):
         '''
         Check if two ImageData objects are equal. This is done by checking if the geometry, data and dtype are equal.
-        Also, if the other object is a np.ndarray, it will check if the data and dtype are equal.
-        
+        Also, if the other object is a ndarray, it will check if the data and dtype are equal.
+        This will check equality as numpy allclose.
+
         Parameters
         ----------
         other: ImageData or np.ndarray
@@ -124,14 +125,17 @@ def __eq__(self, other):
         -------
         bool
             True if the two objects are equal, False otherwise.
-        '''
-
+        ''' 
+        from .array_api_compat import allclose as cil_allclose
         if isinstance(other, ImageData):
-            if np.array_equal(self.as_array(), other.as_array()) \
-                and self.geometry == other.geometry \
-                and self.dtype == other.dtype:
+            xp = array_namespace(self.array)
+            if self.geometry == other.geometry \
+                and self.dtype == other.dtype \
+                and self.shape == other.shape \
+                and cil_allclose(self.as_array(), other.as_array()):
                 return True 
-        elif np.array_equal(self.as_array(), other) and self.dtype==other.dtype:
+            return False
+        elif self.dtype==other.dtype and cil_allclose(self.as_array(), other):
             return True
         else:
             return False
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 5cdee426dd..4efbf1bc01 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -19,12 +19,12 @@
 import logging
 import sys
 
-import array_api_compat.numpy
 import numpy as np
 
 from cil.framework import (DataContainer, ImageGeometry, ImageData, VectorGeometry, AcquisitionData,
                            AcquisitionGeometry, BlockGeometry, VectorData)
 from cil.framework.labels import ImageDimension, AcquisitionDimension
+from cil.framework.array_api_compat import allclose as cil_allclose
 
 from testclass import CCPiTestClass
 from utils import initialise_tests
@@ -36,6 +36,11 @@
     import torch
 except ImportError:
     torch = None
+try:
+    import cupy
+except ImportError:
+    cupy = None
+
 import unittest
 
 log = logging.getLogger(__name__)
@@ -57,6 +62,7 @@ def create_DataContainer(self, X,Y,Z, value=1, backend='numpy'):
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
+         param(cupy, None, None, id="cupy"),
          ]) 
     def test_creation_nocopy(self, xp, raise_error, err_type):
         if xp is None:
@@ -87,6 +93,7 @@ def testGb_creation_nocopy(self):
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
+         param(cupy, None, None, id="cupy"),
          ]) 
     def test_ndim(self, xp, raise_error, err_type):
         if xp is None:
@@ -99,6 +106,7 @@ def test_ndim(self, xp, raise_error, err_type):
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
+         param(cupy, None, None, id="cupy"),
          ]) 
     def test_DataContainer_equal(self, xp, raise_error, err_type):
         if xp is None:
@@ -119,6 +127,7 @@ def test_DataContainer_equal(self, xp, raise_error, err_type):
     @parametrize("xp, raise_error, err_type", 
         [param(numpy, None, None, id="numpy"), 
          param(torch, None, None, id="torch"),
+         param(cupy, None, None, id="cupy"),
          ]) 
     def test_AcquisitionData_equal(self, xp, raise_error, err_type):
         if xp is None:
@@ -160,6 +169,7 @@ def test_AcquisitionData_equal(self, xp, raise_error, err_type):
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
+         param(cupy, None, None, None, id="cupy"),
          ]) 
     def test_ImageData_equal(self, xp, device, raise_error, err_type):
         if xp is None:
@@ -181,13 +191,17 @@ def test_ImageData_equal(self, xp, device, raise_error, err_type):
 
         # Check the equality of two ImageData with different shapes
         data_different_shape = data.copy()
-        data_different_shape.array = data_different_shape.array.reshape(8, 4)
+        # xp.reshape(data_different_shape.array, (8, 4), copy=False) does not seem to work
+        data_different_shape.array = xp.reshape(data_different_shape.array, (8, 4))
+
+
+        # data_different_shape.array = data_different_shape.array.reshape(8, 4)
 
         self.assertFalse(data == data_different_shape)
 
         # Check the equality of two ImageData with different dtypes
-        data_different_dtype = data.geometry.allocate(0, dtype=np.float64)
-        self.assertFalse(data == data_different_dtype)
+        data_different_dtype = data.geometry.allocate(0, dtype=xp.float64)
+        self.assertFalse(cil_allclose(data, data_different_dtype))
 
 
         # Check the equality of two ImageData with different labels
@@ -198,10 +212,11 @@ def test_ImageData_equal(self, xp, device, raise_error, err_type):
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
+         param(cupy, None, None, None, id="cupy"),
          ]) 
     def testInlineAlgebra(self, xp, device, raise_error, err_type):
         if xp is None:
-            self.skipTest("torch not available")
+            self.skipTest(f"xp not available")
         X, Y, Z = 8, 16, 32
         a = np.ones((X, Y, Z), dtype='float32')
         b = np.ones((X, Y, Z), dtype='float32')
@@ -239,6 +254,7 @@ def testInlineAlgebra(self, xp, device, raise_error, err_type):
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
+         param(cupy, None, None, None, id="cupy"),
          ])
     def test_unary_operations(self, xp, device, raise_error, err_type):
         if xp is None:
@@ -1365,6 +1381,7 @@ def test_negation(self):
     @parametrize("xp, device, raise_error, err_type", 
         [param(numpy, None, None, None, id="numpy"), 
          param(torch, None, None, None, id="torch_cpu"),
+         param(cupy, None, None, None, id="cupy"),
          ])
     def test_fill_dimension_ImageData(self, xp, device, raise_error, err_type):
         if xp is None:

From ba6194b303566249daa9cded4442eb4911efd888 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Mon, 16 Jun 2025 15:41:44 +0000
Subject: [PATCH 23/25] WIP for TotalVariation

---
 .../Python/cil/framework/array_api_compat.py  | 23 +++++
 .../Python/cil/framework/data_container.py    | 13 ++-
 Wrappers/Python/cil/framework/image_data.py   | 15 ++-
 .../Python/cil/framework/image_geometry.py    |  2 +-
 .../optimisation/functions/IndicatorBox.py    | 11 ++-
 .../optimisation/functions/TotalVariation.py  |  5 +-
 .../operators/FiniteDifferenceOperator.py     | 52 +++++-----
 experiments/cupy_TV.py                        | 97 +++++++++++++++++++
 8 files changed, 181 insertions(+), 37 deletions(-)
 create mode 100644 experiments/cupy_TV.py

diff --git a/Wrappers/Python/cil/framework/array_api_compat.py b/Wrappers/Python/cil/framework/array_api_compat.py
index 77c77e0887..02634f98b4 100644
--- a/Wrappers/Python/cil/framework/array_api_compat.py
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -1,5 +1,6 @@
 # Updates to array API compatibility layer for CIL
 from array_api_compat import array_namespace
+import sys
 
 def expand_dims(array, axis):
     '''Expand dimensions of an array along specified axes.
@@ -93,3 +94,25 @@ def allclose(a, b, rtol=1e-5, atol=1e-6):
         print(f"Max difference: {diff.max()}")
         return False
     return True
+
+def dtype_namespace(dtype):
+    """
+    Get the namespace of a given dtype.
+    
+    Parameters
+    ----------
+    dtype : data-type
+        The data type to check.
+        
+    Returns
+    -------
+    str
+        The namespace of the dtype.
+        
+    Raises
+    ------
+    TypeError
+        If the dtype is not recognized.
+    """
+    xp = sys.modules[dtype.__module__]
+    return xp
\ No newline at end of file
diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index f843331f33..34b78c1da2 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -480,7 +480,11 @@ def power(self, other, *args, **kwargs):
 
     def maximum(self, x2, *args, **kwargs):
         xp = array_namespace(self.array)
-        return self.pixel_wise_binary(xp.maximum, x2, *args, **kwargs)
+        try:
+            return self.pixel_wise_binary(xp.maximum, x2, *args, **kwargs)
+        except TypeError as te:
+            tmp = xp.ones_like(self.as_array()) * x2
+            return self.pixel_wise_binary(xp.maximum, tmp, *args, **kwargs)
 
     def minimum(self,x2, out=None, *args, **kwargs):
         xp = array_namespace(self.array)
@@ -693,7 +697,12 @@ def sqrt(self, *args,  **kwargs):
 
     def conjugate(self, *args,  **kwargs):
         xp = array_namespace(self.array)
-        return self.pixel_wise_unary(xp.conjugate, *args,  **kwargs)
+        try:
+            return self.pixel_wise_unary(xp.conjugate, *args,  **kwargs)
+        except AttributeError:
+            # torch has conj instead
+            import torch
+            return torch.conj(self.as_array())
 
     def exp(self, *args, **kwargs):
         '''Applies exp pixel-wise to the DataContainer'''
diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index bc687c9901..c1c6bf6635 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -19,6 +19,7 @@
 
 from .data_container import DataContainer
 from .array_api_compat import squeeze as cil_squeeze 
+from .array_api_compat import dtype_namespace
 from .labels import ImageDimension, Backend
 
 from array_api_compat import array_namespace # https://data-apis.org/array-api-compat/
@@ -61,7 +62,7 @@ def dimension_labels(self):
     def dimension_labels(self, val):
         if val is not None:
             raise ValueError("Unable to set the dimension_labels directly. Use geometry.set_labels() instead")
-
+    
     def __init__(self,
                  array = None,
                  deep_copy=False,
@@ -75,9 +76,15 @@ def __init__(self,
                     raise TypeError('dtype must match the array dtype got {} expected {}'.format(dtype, array.dtype))
         # TODO: change
         if array is None:
-            # defaults to a np array
-            xp = np
-            dtype = kwargs.get('dtype', xp.float32)
+            dtype = kwargs.get('dtype', None)
+            if dtype is not None:
+                # tries to infer the namespace from the dtype
+                xp = dtype_namespace(dtype)
+            else:
+                # defaults to a np array
+                xp = np
+                dtype = np.float32
+                
             array = xp.empty(geometry.shape, dtype=dtype)
         elif issubclass(type(array) , DataContainer):
             # this is a reference so we might make a mess modifying both ends
diff --git a/Wrappers/Python/cil/framework/image_geometry.py b/Wrappers/Python/cil/framework/image_geometry.py
index 5bf54d3e27..c31d683d93 100644
--- a/Wrappers/Python/cil/framework/image_geometry.py
+++ b/Wrappers/Python/cil/framework/image_geometry.py
@@ -270,7 +270,7 @@ def allocate(self, value=0, **kwargs):
 
         if isinstance(value, Number):
             # it's created empty, so we make it 0
-            out.array.fill(value)
+            out.fill(value)
         elif value in FillType:
             if value == FillType.RANDOM:
                 seed = kwargs.get('seed', None)
diff --git a/Wrappers/Python/cil/optimisation/functions/IndicatorBox.py b/Wrappers/Python/cil/optimisation/functions/IndicatorBox.py
index 3820b10c76..cf154b89d9 100644
--- a/Wrappers/Python/cil/optimisation/functions/IndicatorBox.py
+++ b/Wrappers/Python/cil/optimisation/functions/IndicatorBox.py
@@ -21,6 +21,7 @@
 import numba
 from cil.utilities import multiprocessing as cil_mp
 import logging
+from array_api_compat import array_namespace
 
 log = logging.getLogger(__name__)
 
@@ -278,11 +279,12 @@ class IndicatorBox_numpy(IndicatorBox):
     def evaluate(self, x):
         '''Evaluates IndicatorBox at x'''
 
-        if (np.all(x.as_array() >= self.lower)
-                and np.all(x.as_array() <= self.upper)):
+        xp = array_namespace(x.as_array())
+        if (xp.all(x.as_array() >= self.lower)
+                and xp.all(x.as_array() <= self.upper)):
             val = 0
         else:
-            val = np.inf
+            val = xp.inf
         return val
 
     def convex_conjugate(self, x):
@@ -290,7 +292,8 @@ def convex_conjugate(self, x):
         return x.maximum(0).sum()
 
     def _proximal(self, outarr):
-        np.clip(outarr,
+        xp = array_namespace(outarr)
+        xp.clip(outarr,
                 None if self.orig_lower is None else self.lower,
                 None if self.orig_upper is None else self.upper,
                 out=outarr)
diff --git a/Wrappers/Python/cil/optimisation/functions/TotalVariation.py b/Wrappers/Python/cil/optimisation/functions/TotalVariation.py
index 6c438c76a8..6caf0a6e9d 100644
--- a/Wrappers/Python/cil/optimisation/functions/TotalVariation.py
+++ b/Wrappers/Python/cil/optimisation/functions/TotalVariation.py
@@ -163,7 +163,8 @@ def __init__(self,
                  isotropic=True,
                  split=False,
                  strong_convexity_constant=0,
-                 warm_start=True):
+                 warm_start=True, 
+                 indicator_accelerated=False):
 
         super(TotalVariation, self).__init__(L=None)
 
@@ -188,7 +189,7 @@ def __init__(self,
             upper = np.inf
         self.lower = lower
         self.upper = upper
-        self.projection_C = IndicatorBox(lower, upper).proximal
+        self.projection_C = IndicatorBox(lower, upper, accelerated=indicator_accelerated).proximal
 
         # Setup GradientOperator as None. This is to avoid domain argument in the __init__
         self._gradient = None
diff --git a/Wrappers/Python/cil/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/cil/optimisation/operators/FiniteDifferenceOperator.py
index 5061ea024e..cec584987e 100644
--- a/Wrappers/Python/cil/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/cil/optimisation/operators/FiniteDifferenceOperator.py
@@ -21,6 +21,8 @@
 from cil.optimisation.operators import LinearOperator
 from cil.utilities.errors import InPlaceError
 
+from array_api_compat import array_namespace
+
 class FiniteDifferenceOperator(LinearOperator):
 
     r'''
@@ -90,6 +92,7 @@ def direct(self, x, out = None):
             raise InPlaceError(message="FiniteDifferenceOperator.direct cannot be used in place")
 
         x_asarr = x.as_array()
+        xp = array_namespace(x_asarr)
 
         outnone = False
         if out is None:
@@ -107,14 +110,14 @@ def direct(self, x, out = None):
         if self.method == 'forward':
 
             # interior nodes
-            np.subtract( x_asarr[tuple(self.get_slice(2, None))], \
+            xp.subtract( x_asarr[tuple(self.get_slice(2, None))], \
                              x_asarr[tuple(self.get_slice(1,-1))], \
                              out = outa[tuple(self.get_slice(1, -1))])
 
             if self.boundary_condition == 'Neumann':
 
                 # left boundary
-                np.subtract(x_asarr[tuple(self.get_slice(1,2))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(1,2))],\
                             x_asarr[tuple(self.get_slice(0,1))],
                             out = outa[tuple(self.get_slice(0,1))])
 
@@ -122,12 +125,12 @@ def direct(self, x, out = None):
             elif self.boundary_condition == 'Periodic':
 
                 # left boundary
-                np.subtract(x_asarr[tuple(self.get_slice(1,2))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(1,2))],\
                             x_asarr[tuple(self.get_slice(0,1))],
                             out = outa[tuple(self.get_slice(0,1))])
 
                 # right boundary
-                np.subtract(x_asarr[tuple(self.get_slice(0,1))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(0,1))],\
                             x_asarr[tuple(self.get_slice(-1,None))],
                             out = outa[tuple(self.get_slice(-1,None))])
 
@@ -141,26 +144,26 @@ def direct(self, x, out = None):
         elif self.method == 'backward':
 
             # interior nodes
-            np.subtract( x_asarr[tuple(self.get_slice(1, -1))], \
+            xp.subtract( x_asarr[tuple(self.get_slice(1, -1))], \
                              x_asarr[tuple(self.get_slice(0,-2))], \
                              out = outa[tuple(self.get_slice(1, -1))])
 
             if self.boundary_condition == 'Neumann':
 
                     # right boundary
-                    np.subtract( x_asarr[tuple(self.get_slice(-1, None))], \
+                    xp.subtract( x_asarr[tuple(self.get_slice(-1, None))], \
                                  x_asarr[tuple(self.get_slice(-2,-1))], \
                                  out = outa[tuple(self.get_slice(-1, None))])
 
             elif self.boundary_condition == 'Periodic':
 
                 # left boundary
-                np.subtract(x_asarr[tuple(self.get_slice(0,1))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(0,1))],\
                             x_asarr[tuple(self.get_slice(-1,None))],
                             out = outa[tuple(self.get_slice(0,1))])
 
                 # right boundary
-                np.subtract(x_asarr[tuple(self.get_slice(-1,None))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(-1,None))],\
                             x_asarr[tuple(self.get_slice(-2,-1))],
                             out = outa[tuple(self.get_slice(-1,None))])
 
@@ -175,7 +178,7 @@ def direct(self, x, out = None):
         elif self.method == 'centered':
 
             # interior nodes
-            np.subtract( x_asarr[tuple(self.get_slice(2, None))], \
+            xp.subtract( x_asarr[tuple(self.get_slice(2, None))], \
                              x_asarr[tuple(self.get_slice(0,-2))], \
                              out = outa[tuple(self.get_slice(1, -1))])
 
@@ -184,13 +187,13 @@ def direct(self, x, out = None):
             if self.boundary_condition == 'Neumann':
 
                 # left boundary
-                np.subtract( x_asarr[tuple(self.get_slice(1, 2))], \
+                xp.subtract( x_asarr[tuple(self.get_slice(1, 2))], \
                                  x_asarr[tuple(self.get_slice(0,1))], \
                                  out = outa[tuple(self.get_slice(0, 1))])
                 outa[tuple(self.get_slice(0, 1))] /=2.
 
                 # left boundary
-                np.subtract( x_asarr[tuple(self.get_slice(-1, None))], \
+                xp.subtract( x_asarr[tuple(self.get_slice(-1, None))], \
                                  x_asarr[tuple(self.get_slice(-2,-1))], \
                                  out = outa[tuple(self.get_slice(-1, None))])
                 outa[tuple(self.get_slice(-1, None))] /=2.
@@ -199,14 +202,14 @@ def direct(self, x, out = None):
                 pass
 
                # left boundary
-                np.subtract( x_asarr[tuple(self.get_slice(1, 2))], \
+                xp.subtract( x_asarr[tuple(self.get_slice(1, 2))], \
                                  x_asarr[tuple(self.get_slice(-1,None))], \
                                  out = outa[tuple(self.get_slice(0, 1))])
                 outa[tuple(self.get_slice(0, 1))] /= 2.
 
 
                 # left boundary
-                np.subtract( x_asarr[tuple(self.get_slice(0, 1))], \
+                xp.subtract( x_asarr[tuple(self.get_slice(0, 1))], \
                                  x_asarr[tuple(self.get_slice(-2,-1))], \
                                  out = outa[tuple(self.get_slice(-1, None))])
                 outa[tuple(self.get_slice(-1, None))] /= 2.
@@ -235,6 +238,7 @@ def adjoint(self, x, out=None):
         # Adjoint operation defined as
 
         x_asarr = x.as_array()
+        xp = array_namespace(x_asarr)
 
         outnone = False
         if out is None:
@@ -254,7 +258,7 @@ def adjoint(self, x, out=None):
         if self.method == 'forward':
 
             # interior nodes
-            np.subtract( x_asarr[tuple(self.get_slice(1, -1))], \
+            xp.subtract( x_asarr[tuple(self.get_slice(1, -1))], \
                              x_asarr[tuple(self.get_slice(0,-2))], \
                              out = outa[tuple(self.get_slice(1, -1))])
 
@@ -269,11 +273,11 @@ def adjoint(self, x, out=None):
             elif self.boundary_condition == 'Periodic':
 
                 # left boundary
-                np.subtract(x_asarr[tuple(self.get_slice(0,1))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(0,1))],\
                             x_asarr[tuple(self.get_slice(-1,None))],
                             out = outa[tuple(self.get_slice(0,1))])
                 # right boundary
-                np.subtract(x_asarr[tuple(self.get_slice(-1,None))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(-1,None))],\
                             x_asarr[tuple(self.get_slice(-2,-1))],
                             out = outa[tuple(self.get_slice(-1,None))])
 
@@ -287,7 +291,7 @@ def adjoint(self, x, out=None):
         elif self.method == 'backward':
 
             # interior nodes
-            np.subtract( x_asarr[tuple(self.get_slice(2, None))], \
+            xp.subtract( x_asarr[tuple(self.get_slice(2, None))], \
                              x_asarr[tuple(self.get_slice(1,-1))], \
                              out = outa[tuple(self.get_slice(1, -1))])
 
@@ -303,12 +307,12 @@ def adjoint(self, x, out=None):
             elif self.boundary_condition == 'Periodic':
 
                 # left boundary
-                np.subtract(x_asarr[tuple(self.get_slice(1,2))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(1,2))],\
                             x_asarr[tuple(self.get_slice(0,1))],
                             out = outa[tuple(self.get_slice(0,1))])
 
                 # right boundary
-                np.subtract(x_asarr[tuple(self.get_slice(0,1))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(0,1))],\
                             x_asarr[tuple(self.get_slice(-1,None))],
                             out = outa[tuple(self.get_slice(-1,None))])
 
@@ -323,7 +327,7 @@ def adjoint(self, x, out=None):
         elif self.method == 'centered':
 
             # interior nodes
-            np.subtract( x_asarr[tuple(self.get_slice(2, None))], \
+            xp.subtract( x_asarr[tuple(self.get_slice(2, None))], \
                              x_asarr[tuple(self.get_slice(0,-2))], \
                              out = outa[tuple(self.get_slice(1, -1))])
             outa[tuple(self.get_slice(1, -1))] /= 2.0
@@ -332,13 +336,13 @@ def adjoint(self, x, out=None):
             if self.boundary_condition == 'Neumann':
 
                 # left boundary
-                np.add(x_asarr[tuple(self.get_slice(0,1))],\
+                xp.add(x_asarr[tuple(self.get_slice(0,1))],\
                             x_asarr[tuple(self.get_slice(1,2))],
                             out = outa[tuple(self.get_slice(0,1))])
                 outa[tuple(self.get_slice(0,1))] /= 2.0
 
                 # right boundary
-                np.add(x_asarr[tuple(self.get_slice(-1,None))],\
+                xp.add(x_asarr[tuple(self.get_slice(-1,None))],\
                             x_asarr[tuple(self.get_slice(-2,-1))],
                             out = outa[tuple(self.get_slice(-1,None))])
 
@@ -348,13 +352,13 @@ def adjoint(self, x, out=None):
             elif self.boundary_condition == 'Periodic':
 
                 # left boundary
-                np.subtract(x_asarr[tuple(self.get_slice(1,2))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(1,2))],\
                             x_asarr[tuple(self.get_slice(-1,None))],
                             out = outa[tuple(self.get_slice(0,1))])
                 outa[tuple(self.get_slice(0,1))] /= 2.0
 
                 # right boundary
-                np.subtract(x_asarr[tuple(self.get_slice(0,1))],\
+                xp.subtract(x_asarr[tuple(self.get_slice(0,1))],\
                             x_asarr[tuple(self.get_slice(-2,-1))],
                             out = outa[tuple(self.get_slice(-1,None))])
                 outa[tuple(self.get_slice(-1,None))] /= 2.0
diff --git a/experiments/cupy_TV.py b/experiments/cupy_TV.py
new file mode 100644
index 0000000000..994895491e
--- /dev/null
+++ b/experiments/cupy_TV.py
@@ -0,0 +1,97 @@
+#%%
+from cil.utilities import dataexample
+from cil.optimisation.functions import TotalVariation
+from cil.plugins.ccpi_regularisation.functions import FGP_TV
+from cil.utilities.display import show2D
+
+import time
+
+N = 10
+num_iter = 1000
+isotropic = True
+
+#%%
+# get an example image
+
+data = {'numpy': dataexample.CAMERA.get(size=(128, 128))}
+geom = data['numpy'].geometry.copy()
+
+import numpy as np
+import torch as to
+import cupy as cp
+from cil.framework import ImageData
+
+
+data['torch'] = ImageData(to.from_numpy(data['numpy'].as_array()), 
+                          geometry=geom.copy())
+
+
+data['cupy'] = ImageData(cp.asarray(data['numpy'].as_array()), 
+                          geometry=geom.copy())
+
+data['c'] = data['numpy']
+
+#%% create a TotalVariation object
+TV = {}
+TV['cupy']  = TotalVariation(max_iteration=num_iter, isotropic=isotropic, backend='numpy')
+TV['torch'] = TotalVariation(max_iteration=num_iter, isotropic=isotropic, backend='numpy')
+TV['numpy'] = TotalVariation(max_iteration=num_iter, isotropic=isotropic, backend='numpy')
+TV['c'] = TotalVariation(max_iteration=num_iter, isotropic=isotropic, backend='c')
+
+
+
+#%%
+fgp = FGP_TV(max_iteration=num_iter, isotropic=isotropic, device='gpu')
+
+#%%
+d1 = fgp.proximal(data['numpy'], tau=1)
+t0 = time.time()
+for _ in range(N):
+    d1 = fgp.proximal(data['numpy'], tau=1, out=d1)
+t1 = time.time()
+dt_fgp = t1-t0
+print (f"Elapsed time: {dt_fgp:.6f} seconds")
+
+#%%
+# d3 = cTV.proximal(x, tau=1)
+
+#%%
+# cupy
+
+def run_TV(data, TV, backend, tau=1, N=10):
+    d = TV[backend].proximal(data[backend], tau=tau)
+    t0 = time.time()
+    for _ in range(N):
+        d = TV[backend].proximal(data[backend], tau=tau, out=d)
+    t1 = time.time()
+    dt = t1-t0
+    print (f"Elapsed time: {dt:.6f} seconds")
+    return ( d.as_array(), dt )
+
+
+#%%
+
+d2 , dt_2 = run_TV(data, TV, 'cupy', tau=1, N=N)
+#%%
+# d3 , dt_3 = run_TV(data, TV, 'torch', tau=1, N=N)
+# #%%
+# d4 , dt_4 = run_TV(data, TV, 'numpy', tau=1, N=N)
+# #%%
+# d5 , dt_5 = run_TV(data, TV, 'c', tau=1, N=N)
+
+# #%%
+# show2D(d4, title=f'np TotalVariation {dt_4}', origin='upper', cmap='inferno')
+
+# #%%
+# show2D([d1, d5, d4, d3.numpy(), cp.asnumpy(d2)], title=[f'FGP_TV {dt_fgp/N:.3f}s', 
+#                                  f'C TotalVariation {dt_5/N:.3f}s',
+#                                  f'np TotalVariation {dt_4/N:.3f}s',
+#                                  f'torch TotalVariation {dt_3/N:.3f}s',
+#                                  f'cp TotalVariation {dt_2/N:.3f}s',
+#                                  ], origin='upper', cmap='inferno',
+#                                  num_cols=2)
+# # %%
+# import array_api_compat
+# np_array = array_api_compat.to_numpy(data['cupy'].as_array())
+
+# # %%

From 012929343f5f5170058905851b60bb2f7fef3fc7 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Mon, 16 Jun 2025 15:57:27 +0000
Subject: [PATCH 24/25] WIP

---
 Wrappers/Python/cil/framework/data_container.py            | 1 -
 .../Python/cil/optimisation/functions/TotalVariation.py    | 7 +++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 34b78c1da2..3c7cf8513b 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -490,7 +490,6 @@ def minimum(self,x2, out=None, *args, **kwargs):
         xp = array_namespace(self.array)
         return self.pixel_wise_binary(xp.minimum, x2=x2, out=out, *args, **kwargs)
 
-
     def sapyb(self, a, y, b, out=None, num_threads=NUM_THREADS):
         '''performs a*self + b * y. Can be done in-place
 
diff --git a/Wrappers/Python/cil/optimisation/functions/TotalVariation.py b/Wrappers/Python/cil/optimisation/functions/TotalVariation.py
index 6caf0a6e9d..9079260ae6 100644
--- a/Wrappers/Python/cil/optimisation/functions/TotalVariation.py
+++ b/Wrappers/Python/cil/optimisation/functions/TotalVariation.py
@@ -271,7 +271,7 @@ def proximal(self, x, tau, out=None):
             tau *= strongly_convex_factor
 
         return solution
-
+    
     def _fista_on_dual_rof(self, x, tau, out=None):
         r""" Runs the Fast Gradient Projection (FGP) algorithm to solve the dual problem
         of the Total Variation Denoising problem (ROF).
@@ -306,7 +306,10 @@ def _fista_on_dual_rof(self, x, tau, out=None):
             tau.multiply(-self.regularisation_parameter, out=tau_reg_neg)
 
         if out is None:
-            out = self.gradient_operator.domain_geometry().allocate(0)
+            # out = self.gradient_operator.domain_geometry().allocate(0)
+            from cil.framework import ImageData
+            import cupy as cp
+            out = ImageData(cp.empty(x.shape, dtype=cp.float32), geometry=self.gradient_operator.domain_geometry().copy(), dtype=cp.float32)
 
         for k in range(self.iterations):
 

From 50c2ffa90ab2fcde5374d990fa3dad4245509713 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Mon, 16 Jun 2025 16:21:11 +0000
Subject: [PATCH 25/25] removed pysnooper

---
 Wrappers/Python/cil/framework/array_api_compat.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Wrappers/Python/cil/framework/array_api_compat.py b/Wrappers/Python/cil/framework/array_api_compat.py
index 02634f98b4..24fd594313 100644
--- a/Wrappers/Python/cil/framework/array_api_compat.py
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -74,8 +74,6 @@ def squeeze(array, axis=None):
         axis = axis[0]
     return squeeze(xp.squeeze(array, axis=ax), axis=axis)
 
-import pysnooper
-@pysnooper.snoop()
 def allclose(a, b, rtol=1e-5, atol=1e-6):
     """
     Check if two arrays are element-wise equal within a tolerance.allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False)[source]