From bc9531d35eebed137fa944d7cd6e00f5307dd5d3 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 17 Feb 2026 16:52:16 +0000
Subject: [PATCH 01/14] attempt to reconcile current master to branch with
 pytorch

---
 .../Python/cil/framework/data_container.py    | 272 ++++++++++++------
 1 file changed, 177 insertions(+), 95 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index f7321550e4..741890b844 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -27,7 +27,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
+from array_api_compat import array_namespace
 
 class DataContainer(object):
     '''Generic class to hold data
@@ -91,15 +93,12 @@ def size(self):
     __container_priority__ = 1
     def __init__ (self, array, deep_copy=True, dimension_labels=None,
                   **kwargs):
-        if type(array) == numpy.ndarray:
-            if deep_copy:
-                self.array = array.copy()
-            else:
-                self.array = array
+        if deep_copy:
+            # self.array = array.copy()
+            xp = array_namespace(array)
+            self.array = xp.asarray(array, copy=True)
         else:
-            raise TypeError('Array must be NumpyArray, passed {0}'\
-                            .format(type(array)))
-
+            self.array = array
         #Don't set for derived classes
         if type(self) is DataContainer:
             self.dimension_labels = dimension_labels
@@ -155,6 +154,11 @@ def as_array(self):
         '''
         return self.array
 
+    def asarray(self):
+        '''
+        Alias to as_array to be compatible with array API. Returns the pointer to the array.
+        '''
+        return self.as_array()
 
     def get_slice(self, **kw):
         '''
@@ -210,7 +214,8 @@ def reorder(self, order):
         for i, axis in enumerate(order):
             new_order[i] = self.dimension_labels.index(axis)
             dimension_labels_new[i] = axis
-
+        
+        # FIXME: works only with NumPy arrays
         self.array = numpy.ascontiguousarray(numpy.transpose(self.array, new_order))
 
         if self.geometry is None:
@@ -293,29 +298,21 @@ def fill(self, array, **kwargs):
         else:
             index = (slice(None),) * self.array.ndim
 
+        if dimension_labels is not None:
+            indexed_dimension_labels = [label for label in self.dimension_labels if label not in dimension]
+            if tuple(indexed_dimension_labels) != array.dimension_labels:
+                raise ValueError('Input array is not in the same order as destination array. Use "array.reorder()"')
+
         if id(array) == id(self.array):
             return
         
-        if isinstance(array, numpy.ndarray):
-            numpy.copyto(self.array[index], array)
+        # if isinstance(array, numpy.ndarray):
+        #     numpy.copyto(self.array[index], array)
         
-        elif isinstance(array, Number):
-            self.array[index] = array
+        # elif isinstance(array, Number):
+        #     self.array[index] = array
         
-        elif issubclass(array.__class__ , DataContainer):
-            if dimension_labels is not None:
-                indexed_dimension_labels = [label for label in self.dimension_labels if label not in dimension]
-                if tuple(indexed_dimension_labels) != array.dimension_labels:
-                    raise ValueError('Input array is not in the same order as destination array. Use "array.reorder()"')
-            
-            if self.array[index].shape == array.shape:
-                numpy.copyto(self.array[index], array.array)
-            else:
-                raise ValueError('Cannot fill with the provided array.' + \
-                                    'Expecting shape {0} got {1}'.format(
-                                    self.array[index].shape, array.shape))
-            
-        elif array in FillType:
+        if array in FillType:
 
             seed = kwargs.pop("seed", None)
             
@@ -363,9 +360,48 @@ def fill(self, array, **kwargs):
                     r = rng.integers(min_value, max_value, size=self.array[index].shape, dtype=numpy.int32).astype(self.dtype)
                 self.array[index] = r
 
+        # else:
+        #     raise TypeError('Can fill only with random method, number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
+        if issubclass(array.__class__ , DataContainer):
+            return self.fill(array.as_array(), **kwargs)
+
+
+        if dimension == {}:
+            import warnings
+            warnings.filterwarnings('error')
+            xp = array_namespace(self.as_array())
+            self.array.__setitem__(slice(None, None, None), array)
+            warnings.resetwarnings()
         else:
-            raise TypeError('Can fill only with random method, number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
 
+            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[i] = slice(v,v+1,None)
+                        where.append(i)
+
+            if self.array[index].shape == array.shape:
+                # numpy.copyto(self.array[index], array.array)
+                try:
+                 
+                    array = cil_expand_dims(array, axis=where)
+                    self.array.__setitem__(tuple(slices), array)
+                    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)
+                    
+            else:
+                raise ValueError('Cannot fill with the provided array.' + \
+                                    'Expecting shape {0} got {1}'.format(
+                                    self.array[index].shape, array.shape))
+            
+        
         if kwargs:
             warnings.warn(f"Unused keyword arguments: {kwargs}", stacklevel=2)
 
@@ -495,69 +531,71 @@ 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)))
-            
-            if self.geometry is None:
-                return type(self)(out,
-                    deep_copy=False,
-                    dimension_labels=self.dimension_labels)
+                out = pwop(self.as_array() , x2 , *args, **kwargs )
             
+            geom = self.geometry
+            if geom is not None:
+                geom = self.geometry.copy()
             return type(self)(out,
-                deep_copy=False,
-                geometry=self.geometry)
-
-
-        elif issubclass(type(out), DataContainer) and issubclass(type(x2), DataContainer):
-            if self.check_dimensions(out) and self.check_dimensions(x2):
-                kwargs['out'] = out.as_array()
-                pwop(self.as_array(), x2.as_array(), *args, **kwargs )
-                #return type(self)(out.as_array(),
-                #       deep_copy=False,
-                #       dimension_labels=self.dimension_labels,
-                #       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}")
-                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):
-            if self.array.shape == out.shape and self.array.dtype == out.dtype:
-                kwargs['out'] = out
-                pwop(self.as_array(), x2, *args, **kwargs)
-                #return type(self)(out,
-                #       deep_copy=False,
-                #       dimension_labels=self.dimension_labels,
-                #       geometry=self.geometry)
+                   deep_copy=False,
+                   dimension_labels=self.dimension_labels,
+                   geometry= None if self.geometry is None else self.geometry.copy(),
+                   suppress_warning=True)
         else:
-            raise ValueError(f"incompatible class: {pwop.__name__} {type(out)}")
+            # check the size and dimension of out
+            whatswrong = {}
+            if not self.check_dimensions(out):
+                whatswrong['out shape'] = out.shape
+            if not out.dtype == self.dtype:
+                whatswrong['out dtype'] = out.dtype
+            if not whatswrong == {}:
+                # report to the user what's wrong
+                msg = "Wrong size for data memory:\n"
+                for k,v in whatswrong.items():
+                    msg += f"{k} {v}\n"
+                raise ValueError(msg)
+            
+            if isinstance(x2, Number):
+                out = pwop(self.as_array() , x2 , *args, **kwargs )
+            else:
+                whatswrong = {}
+                if not self.check_dimensions(x2):
+                    whatswrong['x2 shape'] = x2.shape
+                if not out.dtype == self.dtype:
+                    whatswrong['x2 dtype'] = x2.dtype
+                if not whatswrong == {}:
+                    # report to the user what's wrong
+                    msg = "Wrong size for data memory:\n"
+                    for k,v in whatswrong.items():
+                        msg += f"{k} {v}\n"
+                    raise ValueError(msg)
+                if issubclass(x2.__class__ , DataContainer):
+                    out = pwop(self.as_array() , x2.as_array() , *args, **kwargs )
+                else:
+                    out = pwop(self.as_array() , x2 , *args, **kwargs )
+                return 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.sapyb(-1,self,1, out=kwargs.get('out', None))
-        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 \
@@ -565,16 +603,24 @@ def divide(self, other, *args, **kwargs):
             _out = other.divide(self, *args, **kwargs)
             _out.power(-1, out=_out)
             return _out
-        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)
+        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):
-        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):
@@ -744,44 +790,78 @@ 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):
-            if self.array.shape == out.shape and self.array.dtype == out.dtype:
-                kwargs['out'] = out
-                pwop(self.as_array(), *args, **kwargs)
+        
+        # elif issubclass(type(out), numpy.ndarray):
         else:
-            raise ValueError("incompatible class: {pwop.__name__} {type(out)}")
-
+            # check the size and dimension of out
+            whatswrong = {}
+            if not self.check_dimensions(out):
+                whatswrong['out shape'] = out.shape
+            if not out.dtype == self.dtype:
+                whatswrong['out dtype'] = out.dtype
+            if not whatswrong == {}:
+                # report to the user what's wrong
+                msg = "Wrong size for data memory:\n"
+                for k,v in whatswrong.items():
+                    msg += f"{k} {v}\n"
+                raise ValueError(msg)
+            
+            # FIXME is this necessary? Can one pass a numpy array in out?
+            kwargs['out'] = out
+            pwop(self.as_array(), *args, **kwargs)
+            
     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)
+        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'''
-        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):
         '''return the squared euclidean norm of the DataContainer viewed as a vector'''
-        #shape = self.shape
-        #size = reduce(lambda x,y:x*y, shape, 1)
-        #y = numpy.reshape(self.as_array(), (size, ))
         return self.dot(self)
-        #return self.dot(self)
+        
     def norm(self, **kwargs):
         '''return the euclidean norm of the DataContainer viewed as a vector'''
         return numpy.sqrt(self.squared_norm(**kwargs))
@@ -796,6 +876,7 @@ def dot(self, other, *args, **kwargs):
                     method))
 
         if self.shape == other.shape:
+            # FIXME this seems not to have been addressed
             if method == 'numpy':
                 return numpy.dot(self.as_array().ravel(), other.as_array().ravel().conjugate())
             elif method == 'reduce':
@@ -833,6 +914,7 @@ def _directional_reduction_unary(self, reduction_function, axis=None, out=None,
         scalar or ndarray
             The result of the unary function
         """
+        # FIXME Address this
         if axis is not None:
             axis = self.get_dimension_axis(axis)
 

From 0303cbf1cc04e4d82ea3820a762f6da18fd9c456 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Tue, 17 Feb 2026 17:04:33 +0000
Subject: [PATCH 02/14] a few updates

---
 .../Python/cil/framework/acquisition_data.py  |  14 +--
 .../Python/cil/framework/array_api_compat.py  | 116 ++++++++++++++++++
 Wrappers/Python/cil/framework/image_data.py   |  36 ++++--
 recipe/meta.yaml                              |   1 +
 4 files changed, 150 insertions(+), 17 deletions(-)
 create mode 100644 Wrappers/Python/cil/framework/array_api_compat.py

diff --git a/Wrappers/Python/cil/framework/acquisition_data.py b/Wrappers/Python/cil/framework/acquisition_data.py
index 308f0d5462..404e6362c4 100644
--- a/Wrappers/Python/cil/framework/acquisition_data.py
+++ b/Wrappers/Python/cil/framework/acquisition_data.py
@@ -22,6 +22,8 @@
 from .labels import AcquisitionDimension, Backend, AcquisitionType
 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):
@@ -83,17 +85,15 @@ def __init__(self,
         if array is None:
             if dtype is None:
                 dtype = geometry.dtype
-            array = numpy.empty(geometry.shape, dtype)
+            xp = array_api_compat.numpy
+            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)))
+            # remove singleton dimensions
+            array = array.squeeze()
 
         if array.shape != geometry.shape:
             raise ValueError('Shape mismatch got {} expected {}'.format(array.shape, geometry.shape))
@@ -118,7 +118,7 @@ def __eq__(self, other):
         bool
             True if the two objects are equal, False otherwise.
         '''
-
+        # FIXME: address this
         if isinstance(other, AcquisitionData):
             if numpy.array_equal(self.as_array(), other.as_array()) \
                 and self.geometry == other.geometry \
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..24fd594313
--- /dev/null
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -0,0 +1,116 @@
+# 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.
+    
+    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)
+    if len(axis) == 1:
+        return xp.expand_dims(array, axis=axis[0])
+    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 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)
+    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)
+    ax = axis.pop(0)
+    if len(axis) == 1:
+        axis = axis[0]
+    return squeeze(xp.squeeze(array, axis=ax), axis=axis)
+
+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
+
+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/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index 67cfbff5b7..2e8df2d486 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -21,6 +21,11 @@
 from .data_container import DataContainer
 from .labels import ImageDimension, Backend
 
+from .array_api_compat import squeeze as cil_squeeze 
+from .array_api_compat import dtype_namespace
+import array_api_compat
+from array_api_compat import array_namespace # https://data-apis.org/array-api-compat/
+
 class ImageData(DataContainer):
     """
     DataContainer for holding 2D or 3D image data
@@ -76,17 +81,23 @@ def __init__(self,
         if array is None:
             if dtype is None:
                 dtype = geometry.dtype
-            array = numpy.empty(geometry.shape, dtype=dtype)
-
+            xp = array_api_compat.numpy
+            array = xp.empty(geometry.shape, dtype=dtype)
+    
         elif issubclass(type(array) , DataContainer):
             array = array.as_array()
+            array = cil_squeeze(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:
-            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
+            array = cil_squeeze(array)
+        # else:
+        #     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))
@@ -115,12 +126,16 @@ def __eq__(self, other):
             True if the two objects are equal, False otherwise.
         '''
 
+        from .array_api_compat import allclose as cil_allclose
         if isinstance(other, ImageData):
-            if numpy.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 numpy.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
@@ -180,6 +195,7 @@ def apply_circular_mask(self, radius=0.99, in_place=True):
             If `in_place = False` returns a new ImageData object with the masked data
 
         """
+        # FIXME Address this
         ig = self.geometry
 
         # grid
diff --git a/recipe/meta.yaml b/recipe/meta.yaml
index 24d414883f..f17fdd4fb9 100644
--- a/recipe/meta.yaml
+++ b/recipe/meta.yaml
@@ -62,6 +62,7 @@ requirements:
     - cil-data >=22
     - tqdm
     - numba
+    - array-api-compat
 
   # version constraints of optional dependencies
   run_constrained:

From 2894fc790e240de094fba63b6e09c2b44c6706e2 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Thu, 19 Feb 2026 11:02:29 +0000
Subject: [PATCH 03/14] fix unittests

---
 .../Python/cil/framework/data_container.py    | 23 +++++++++++--------
 1 file changed, 14 insertions(+), 9 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 741890b844..8cce2ddbd6 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -298,10 +298,6 @@ def fill(self, array, **kwargs):
         else:
             index = (slice(None),) * self.array.ndim
 
-        if dimension_labels is not None:
-            indexed_dimension_labels = [label for label in self.dimension_labels if label not in dimension]
-            if tuple(indexed_dimension_labels) != array.dimension_labels:
-                raise ValueError('Input array is not in the same order as destination array. Use "array.reorder()"')
 
         if id(array) == id(self.array):
             return
@@ -363,6 +359,10 @@ def fill(self, array, **kwargs):
         # else:
         #     raise TypeError('Can fill only with random method, number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
         if issubclass(array.__class__ , DataContainer):
+            if dimension_labels is not None:
+                indexed_dimension_labels = [label for label in self.dimension_labels if label not in dimension]
+                if tuple(indexed_dimension_labels) != array.dimension_labels:
+                    raise ValueError('Input array is not in the same order as destination array. Use "array.reorder()"')
             return self.fill(array.as_array(), **kwargs)
 
 
@@ -524,7 +524,7 @@ def copy(self):
     ## binary operations
 
     def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
-        out = kwargs.get('out', None)
+        out = kwargs.pop('out', None)
 
         if out is None:
             if isinstance(x2, Number):
@@ -543,6 +543,9 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
                    geometry= None if self.geometry is None else self.geometry.copy(),
                    suppress_warning=True)
         else:
+            # assuming it is a CIL DataContainer we extract the actual data structure
+            outarr = out.as_array() 
+            kwargs['out'] = outarr
             # check the size and dimension of out
             whatswrong = {}
             if not self.check_dimensions(out):
@@ -557,7 +560,7 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
                 raise ValueError(msg)
             
             if isinstance(x2, Number):
-                out = pwop(self.as_array() , x2 , *args, **kwargs )
+                pwop(self.as_array() , x2 , *args, **kwargs )
             else:
                 whatswrong = {}
                 if not self.check_dimensions(x2):
@@ -571,10 +574,12 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
                         msg += f"{k} {v}\n"
                     raise ValueError(msg)
                 if issubclass(x2.__class__ , DataContainer):
-                    out = pwop(self.as_array() , x2.as_array() , *args, **kwargs )
+                    pwop(self.as_array() , x2.as_array() , *args, **kwargs )
                 else:
-                    out = pwop(self.as_array() , x2 , *args, **kwargs )
-                return out
+                    pwop(self.as_array() , x2 , *args, **kwargs )
+
+            out.fill(outarr)
+            return out
 
     def add(self, other, *args, **kwargs):
         if hasattr(other, '__container_priority__') and \

From fe33d21100dfcf06430a12f789c30557d1f9dc8d Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Thu, 19 Feb 2026 12:34:23 +0000
Subject: [PATCH 04/14] wip passing unittests

---
 .../Python/cil/framework/data_container.py    |  9 +++++++--
 Wrappers/Python/cil/framework/image_data.py   | 20 +++++++++----------
 Wrappers/Python/test/test_DataContainer.py    | 14 ++++++-------
 3 files changed, 24 insertions(+), 19 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 8cce2ddbd6..3899218646 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -355,6 +355,7 @@ def fill(self, array, **kwargs):
                 else:
                     r = rng.integers(min_value, max_value, size=self.array[index].shape, dtype=numpy.int32).astype(self.dtype)
                 self.array[index] = r
+            return self
 
         # else:
         #     raise TypeError('Can fill only with random method, number, numpy array or DataContainer and subclasses. Got {}'.format(type(array)))
@@ -550,11 +551,13 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
             whatswrong = {}
             if not self.check_dimensions(out):
                 whatswrong['out shape'] = out.shape
+                whatswrong['expected shape'] = self.shape
             if not out.dtype == self.dtype:
                 whatswrong['out dtype'] = out.dtype
+                whatswrong['expected dtype'] = self.dtype
             if not whatswrong == {}:
                 # report to the user what's wrong
-                msg = "Wrong size for data memory:\n"
+                msg = "Wrong size or type for data memory:\n"
                 for k,v in whatswrong.items():
                     msg += f"{k} {v}\n"
                 raise ValueError(msg)
@@ -565,11 +568,13 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
                 whatswrong = {}
                 if not self.check_dimensions(x2):
                     whatswrong['x2 shape'] = x2.shape
+                    whatswrong['expected shape'] = self.shape
                 if not out.dtype == self.dtype:
                     whatswrong['x2 dtype'] = x2.dtype
+                    whatswrong['expected dtype'] = self.dtype
                 if not whatswrong == {}:
                     # report to the user what's wrong
-                    msg = "Wrong size for data memory:\n"
+                    msg = "Wrong size or type for data memory:\n"
                     for k,v in whatswrong.items():
                         msg += f"{k} {v}\n"
                     raise ValueError(msg)
diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py
index 2e8df2d486..021a0eb1a0 100644
--- a/Wrappers/Python/cil/framework/image_data.py
+++ b/Wrappers/Python/cil/framework/image_data.py
@@ -15,7 +15,7 @@
 #
 # Authors:
 # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
-import numpy
+import numpy as np
 import warnings
 
 from .data_container import DataContainer
@@ -81,7 +81,7 @@ def __init__(self,
         if array is None:
             if dtype is None:
                 dtype = geometry.dtype
-            xp = array_api_compat.numpy
+            xp = np
             array = xp.empty(geometry.shape, dtype=dtype)
     
         elif issubclass(type(array) , DataContainer):
@@ -156,7 +156,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)
@@ -202,10 +202,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
@@ -217,17 +217,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
@@ -244,13 +244,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
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index cb8bcf3d32..a9749d5c4d 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -1060,15 +1060,15 @@ def test_multiply_out(self):
 
     def test_sapyb_datacontainer_f(self):
         #a vec, b vec
-
-        ig = ImageGeometry(10,10)
+        N,M=2,2
+        ig = ImageGeometry(N,M)
         d1 = ig.allocate(dtype=np.float32)
         d2 = ig.allocate(dtype=np.float32)
         a = ig.allocate(dtype=np.float32)
         b = ig.allocate(dtype=np.float32)
 
-        d1.fill(np.asarray(np.arange(1,101).reshape(10,10), dtype=np.float32))
-        d2.fill(np.asarray(np.arange(1,101).reshape(10,10), dtype=np.float32))
+        d1.fill(np.asarray(np.arange(1,N*M+1).reshape(N,M), dtype=np.float32))
+        d2.fill(np.asarray(np.arange(1,N*M+1).reshape(N,M), dtype=np.float32))
         a.fill(1.0/d1.as_array())
         b.fill(-1.0/d2.as_array())
 
@@ -1076,13 +1076,13 @@ def test_sapyb_datacontainer_f(self):
         # equals to 1 + -1 = 0
         out = d1.sapyb(a,d2,b)
         res = np.zeros_like(d1.as_array())
-        np.testing.assert_array_equal(res, out.as_array())
+        atol = np.finfo(np.float32).eps
+        np.testing.assert_allclose(res, out.as_array(), atol=atol)
 
         out.fill(0)
         d1.sapyb(a,d2,b, out)
         res = np.zeros_like(d1.as_array())
-        np.testing.assert_array_equal(res, out.as_array())
-
+        np.testing.assert_allclose(res, out.as_array(), atol=atol)
 
     def test_sapyb_scalar_f(self):
         # a,b scalar

From d253e0b9cb3d17ec2b83a6c0c977149784cde785 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Thu, 19 Feb 2026 14:36:05 +0000
Subject: [PATCH 05/14] wip

---
 .../Python/cil/framework/acquisition_data.py  | 10 +++---
 .../Python/cil/framework/data_container.py    | 32 ++++++++++---------
 2 files changed, 22 insertions(+), 20 deletions(-)

diff --git a/Wrappers/Python/cil/framework/acquisition_data.py b/Wrappers/Python/cil/framework/acquisition_data.py
index 404e6362c4..c03362ffc2 100644
--- a/Wrappers/Python/cil/framework/acquisition_data.py
+++ b/Wrappers/Python/cil/framework/acquisition_data.py
@@ -16,7 +16,7 @@
 # Authors:
 # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
 # Joshua DM Hellier (University of Manchester) [refactorer]
-import numpy
+import numpy as np
 import warnings
 
 from .labels import AcquisitionDimension, Backend, AcquisitionType
@@ -85,7 +85,7 @@ def __init__(self,
         if array is None:
             if dtype is None:
                 dtype = geometry.dtype
-            xp = array_api_compat.numpy
+            xp = np
             array = xp.empty(geometry.shape, dtype=dtype)
     
         elif issubclass(type(array) , DataContainer):
@@ -120,11 +120,11 @@ def __eq__(self, other):
         '''
         # FIXME: address this
         if isinstance(other, AcquisitionData):
-            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
@@ -149,7 +149,7 @@ def _get_slice(self, **kwargs):
             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
             kwargs['vertical'] = ind0
             out = DataContainer.get_slice(self, **kwargs)
diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 3899218646..25ef7f16be 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -364,6 +364,10 @@ def fill(self, array, **kwargs):
                 indexed_dimension_labels = [label for label in self.dimension_labels if label not in dimension]
                 if tuple(indexed_dimension_labels) != array.dimension_labels:
                     raise ValueError('Input array is not in the same order as destination array. Use "array.reorder()"')
+            # slice dimensions parameters have been stripped out of kwargs,
+            # Put them back for fill to use.
+            for k,v in dimension.items():
+                kwargs[k] = v
             return self.fill(array.as_array(), **kwargs)
 
 
@@ -383,26 +387,24 @@ def fill(self, array, **kwargs):
                         slices[i] = slice(v,v+1,None)
                         where.append(i)
 
-            if self.array[index].shape == array.shape:
-                # numpy.copyto(self.array[index], array.array)
-                try:
-                 
+            try:
+                # if it's an array try to fill, otherwise catch the error and try to set the item directly (e.g. for numbers)
+                if self.array[index].shape == array.shape:
+                    # numpy.copyto(self.array[index], array.array)     
                     array = cil_expand_dims(array, axis=where)
                     self.array.__setitem__(tuple(slices), array)
                     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)
-                    
-            else:
-                raise ValueError('Cannot fill with the provided array.' + \
+                else:
+                    raise ValueError('Cannot fill with the provided array.' + \
                                     'Expecting shape {0} got {1}'.format(
-                                    self.array[index].shape, array.shape))
-            
+                                    self.array[index].shape, array.shape))   
         
+            except AttributeError as ae:
+                self.array.__setitem__(tuple(slices), array)
+                    
+            except TypeError as ae:
+                self.array.__setitem__(tuple(slices), array)
+            
         if kwargs:
             warnings.warn(f"Unused keyword arguments: {kwargs}", stacklevel=2)
 

From fe4ae0e57ab034a33ff277016f6b77a0a5aa1c89 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Thu, 19 Feb 2026 16:33:25 +0000
Subject: [PATCH 06/14] wip tests

---
 .../Python/cil/framework/array_api_compat.py  | 21 ++++++++++++++-----
 .../Python/cil/framework/data_container.py    |  1 +
 2 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/Wrappers/Python/cil/framework/array_api_compat.py b/Wrappers/Python/cil/framework/array_api_compat.py
index 24fd594313..d830b0ec54 100644
--- a/Wrappers/Python/cil/framework/array_api_compat.py
+++ b/Wrappers/Python/cil/framework/array_api_compat.py
@@ -77,18 +77,29 @@ def squeeze(array, axis=None):
 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]
+    parameters:
+    a, b: DataContainer or array_like
 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:
+    try:
+        aarr = a.as_array()
+    except AttributeError:
+        aarr = a
+    try:
+        barr = b.as_array()
+    except AttributeError:
+        barr = b
+
+    xp = array_namespace(aarr)
+    if array_namespace(barr) != xp:
         raise TypeError('Can only compare arrays ' \
-        'with same namespace. Got {} and {}'.format(array_namespace(a), array_namespace(b)))
+        'with same namespace. Got {} and {}'.format(array_namespace(aarr), array_namespace(barr)))
     
-    diff = rtol * xp.abs(b) + atol
-    if xp.any(diff < xp.abs(a - b)):
+    diff = rtol * xp.abs(barr) + atol
+    if xp.any(diff < xp.abs(aarr - barr)):
         print(f"Max difference: {diff.max()}")
         return False
     return True
diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 25ef7f16be..610345f946 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -276,6 +276,7 @@ def fill(self, array, **kwargs):
         will copy the data in new_data into the data container.
         '''
         from cil.framework.labels import FillType
+        import warnings
         dimension = {}
 
         try:

From 155f42847966dba180f5e20e907ffdd20b8b9778 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Thu, 19 Feb 2026 18:06:02 +0000
Subject: [PATCH 07/14] add skip tests if no IPP

---
 Wrappers/Python/test/test_DataProcessor.py | 12 +++--
 Wrappers/Python/test/test_functions.py     | 60 ++++++++++++++++------
 Wrappers/Python/test/test_out_in_place.py  |  3 +-
 3 files changed, 53 insertions(+), 22 deletions(-)

diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py
index 32914873c2..b2758b7c83 100644
--- a/Wrappers/Python/test/test_DataProcessor.py
+++ b/Wrappers/Python/test/test_DataProcessor.py
@@ -459,7 +459,7 @@ def test_process_image_geometry(self):
 
         self.assertEqual(ig_gold, ig_out, msg="Binning image geometry with offset roi failed")
 
-
+    @unittest.skipUnless(has_ipp, "Requires IPP libraries")
     def test_bin_array_consistency(self):
 
         ig = ImageGeometry(64,32,16,channels=8)
@@ -496,7 +496,8 @@ def test_bin_array_consistency(self):
 
         numpy.testing.assert_allclose(binned_by_hand,binned_arr_numpy,atol=1e-6)
         numpy.testing.assert_allclose(binned_by_hand,binned_arr_acc,atol=1e-6)
-
+    
+    @unittest.skipUnless(has_ipp, "Requires IPP libraries")
     def test_bin_image_data(self):
         """
         Binning results tested with test_binning_cpp_ so this is checking wrappers with axis labels and geometry
@@ -545,6 +546,7 @@ def test_bin_image_data(self):
 
 
 
+    @unittest.skipUnless(has_ipp, "Requires IPP libraries")
     def test_bin_acquisition_data(self):
         """
         Binning results tested with test_binning_cpp_ so this is checking wrappers with axis labels and geometry
@@ -590,7 +592,7 @@ def test_bin_acquisition_data(self):
         self.assertEqual(binned_data.geometry, binned_by_hand.geometry)
 
 
-
+    @unittest.skipUnless(has_ipp, "Requires IPP libraries")
     def test_process_acquisition(self):
 
         arr=numpy.arange(24,dtype=numpy.float32).reshape(2,3,4)
@@ -624,7 +626,7 @@ def test_process_acquisition(self):
         self.assertEqual(data_out.geometry, geometry_gold,
         msg="Binner failed with geometry mismatch. Got:\n{0}\nExpected:\n{1}".format(data_out.geometry, geometry_gold))
 
-
+    @unittest.skipUnless(has_ipp, "Requires IPP libraries")
     def test_process_image(self):
 
         arr=numpy.arange(24,dtype=numpy.float32).reshape(2,3,4)
@@ -658,7 +660,7 @@ def test_process_image(self):
         self.assertEqual(data_out.geometry, geometry_gold,
         msg="Binner failed with geometry mismatch. Got:\n{0}\nExpected:\n{1}".format(data_out.geometry, geometry_gold))
 
-
+    @unittest.skipUnless(has_ipp, "Requires IPP libraries")
     def test_process_data_container(self):
 
         arr=numpy.arange(24,dtype=numpy.float32).reshape(2,3,4)
diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py
index 61e65c03b9..8545cd4a83 100644
--- a/Wrappers/Python/test/test_functions.py
+++ b/Wrappers/Python/test/test_functions.py
@@ -1123,23 +1123,29 @@ def test_soft_shrinkage(self):
     def soft_shrinkage_test(self, x):
         tau = 1.
         ret = soft_shrinkage(x, tau)
-        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()), atol=atol)
         tau = 2.
         ret = soft_shrinkage(x, tau)
-        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()), atol=atol)
         tau = -1.
         ret = soft_shrinkage(x, tau)
-        np.testing.assert_allclose(ret.as_array(), 2 * np.ones_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), 2 * np.ones_like(x.as_array()), atol=atol)
         tau = 1.
         ret = soft_shrinkage(-0.5 * x, tau)
-        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()), atol=atol)
         tau = 2.
         ret = soft_shrinkage(-0.5 *x, tau)
-        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()), atol=atol)
         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()))
+            atol = np.finfo(ret.as_array().dtype).eps
+            np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()), atol=atol)
         tau = -3j
         with self.assertRaises(ValueError):
             ret = soft_shrinkage(-0.5 *x, tau)
@@ -1147,42 +1153,64 @@ def soft_shrinkage_test(self, x):
         # tau np.ndarray
         tau = 1. * np.ones_like(x.as_array())
         ret = soft_shrinkage(x, tau)
+        atol = np.finfo(ret.as_array().dtype).eps
         np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()))
+        
         tau = 2.* np.ones_like(x.as_array())
         ret = soft_shrinkage(x, tau)
-        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()), atol=atol)
+        
         tau = -1.* np.ones_like(x.as_array())
         ret = soft_shrinkage(x, tau)
+        atol = np.finfo(ret.as_array().dtype).eps
         np.testing.assert_allclose(ret.as_array(), 2 * np.ones_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 * np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()), atol=atol)
+        
         tau = 2.* np.ones_like(x.as_array())
         ret = soft_shrinkage(-0.5 *x, tau)
-        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()), atol=atol)
+        
         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()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()), atol=atol)
 
         # tau DataContainer
         tau = 1. * x
         ret = soft_shrinkage(x, tau)
-        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()), atol=atol)
+        
         tau = 2. * x
         ret = soft_shrinkage(x, tau)
-        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), np.zeros_like(x.as_array()), atol=atol)
+        
         tau = -1. * x
         ret = soft_shrinkage(x, tau)
-        np.testing.assert_allclose(ret.as_array(), 2 * np.ones_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), 2 * np.ones_like(x.as_array()), atol=atol)
+        
         tau = 1. * x
         ret = soft_shrinkage(-0.5 * x, tau)
-        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()), atol=atol)
+        
         tau = 2. * x
         ret = soft_shrinkage(-0.5 *x, tau)
-        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1 * np.zeros_like(x.as_array()), atol=atol)
+        
         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()))
+        atol = np.finfo(ret.as_array().dtype).eps
+        np.testing.assert_allclose(ret.as_array(), -1.5 * np.ones_like(x.as_array()), atol=atol)
 
         np.testing.assert_allclose(ret.as_array().imag, np.zeros_like(ret.as_array().imag), atol=1e-6, rtol=1e-6)
 
diff --git a/Wrappers/Python/test/test_out_in_place.py b/Wrappers/Python/test/test_out_in_place.py
index 6961f081b1..1cc46c99a8 100644
--- a/Wrappers/Python/test/test_out_in_place.py
+++ b/Wrappers/Python/test/test_out_in_place.py
@@ -60,7 +60,7 @@
 
 from utils import  initialise_tests
 
-from utils import has_astra, has_tigre, has_nvidia
+from utils import has_astra, has_tigre, has_nvidia, has_ipp
 
 
 
@@ -521,6 +521,7 @@ def in_place_check(self, processor, data, data_array_index, *args):
             '\nFor geometry type: \n' + str(data.geometry)
             raise type(e)(error_message + '\n\n' + str(e))
 
+    @unittest.skipUnless(has_ipp, "IPP not installed")
     def test_out(self):
         """
         Tests to check output from Processors, including:

From f1bdaa5604412592ab497654053609fcff7f5a58 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Fri, 20 Feb 2026 09:21:40 +0000
Subject: [PATCH 08/14] raise exception only if backend raises

---
 .../Python/cil/framework/data_container.py    | 80 ++++++++++---------
 1 file changed, 44 insertions(+), 36 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index 610345f946..e8bf2aece5 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -550,42 +550,50 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
             # assuming it is a CIL DataContainer we extract the actual data structure
             outarr = out.as_array() 
             kwargs['out'] = outarr
-            # check the size and dimension of out
-            whatswrong = {}
-            if not self.check_dimensions(out):
-                whatswrong['out shape'] = out.shape
-                whatswrong['expected shape'] = self.shape
-            if not out.dtype == self.dtype:
-                whatswrong['out dtype'] = out.dtype
-                whatswrong['expected dtype'] = self.dtype
-            if not whatswrong == {}:
-                # report to the user what's wrong
-                msg = "Wrong size or type for data memory:\n"
-                for k,v in whatswrong.items():
-                    msg += f"{k} {v}\n"
-                raise ValueError(msg)
-            
+                        
             if isinstance(x2, Number):
-                pwop(self.as_array() , x2 , *args, **kwargs )
-            else:
-                whatswrong = {}
-                if not self.check_dimensions(x2):
-                    whatswrong['x2 shape'] = x2.shape
-                    whatswrong['expected shape'] = self.shape
-                if not out.dtype == self.dtype:
-                    whatswrong['x2 dtype'] = x2.dtype
-                    whatswrong['expected dtype'] = self.dtype
-                if not whatswrong == {}:
-                    # report to the user what's wrong
-                    msg = "Wrong size or type for data memory:\n"
-                    for k,v in whatswrong.items():
-                        msg += f"{k} {v}\n"
-                    raise ValueError(msg)
-                if issubclass(x2.__class__ , DataContainer):
-                    pwop(self.as_array() , x2.as_array() , *args, **kwargs )
-                else:
+                try:
                     pwop(self.as_array() , x2 , *args, **kwargs )
-
+                except Exception as bknerr:
+                    # check the size and dimension of out
+                    whatswrong = {}
+                    msg = ""
+                    if not self.check_dimensions(out):
+                        whatswrong['out shape'] = out.shape
+                        whatswrong['expected shape'] = self.shape
+                    if not out.dtype == self.dtype:
+                        whatswrong['out dtype'] = out.dtype
+                        whatswrong['expected dtype'] = self.dtype
+                    if not whatswrong == {}:
+                        # report to the user what's wrong
+                        msg = "Wrong size or type for data memory:\n"
+                        for k,v in whatswrong.items():
+                            msg += f"{k} {v}\n"
+                        print (msg)
+                    # report the error from the backend plus the error from the checks
+                    raise
+            else:
+                try:
+                    if issubclass(x2.__class__ , DataContainer):
+                        pwop(self.as_array() , x2.as_array() , *args, **kwargs )
+                    else:
+                        pwop(self.as_array() , x2 , *args, **kwargs )
+                except Exception as bknerr:
+                    whatswrong = {}
+                    if not self.check_dimensions(x2):
+                        whatswrong['x2 shape'] = x2.shape
+                        whatswrong['expected shape'] = self.shape
+                    if not out.dtype == self.dtype:
+                        whatswrong['x2 dtype'] = x2.dtype
+                        whatswrong['expected dtype'] = self.dtype
+                    if not whatswrong == {}:
+                        # report to the user what's wrong
+                        msg = "Wrong size or type for data memory:\n"
+                        for k,v in whatswrong.items():
+                            msg += f"{k} {v}\n"
+                        print(msg)
+                    # report the error from the backend plus the error from the checks
+                    raise
             out.fill(outarr)
             return out
 
@@ -670,9 +678,9 @@ 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("RuntimeError in sapyb, defaulting to Python due to: {}".format(rte))
             except TypeError as te:
-                warnings.warn("sapyb defaulting to Python due to: {}".format(te))
+                warnings.warn("TypeError in sapyb, defaulting to Python due to: {}".format(te))
             finally:
                 pass
 

From bb3a66688d41ce9d1ab3442a023693e6f539d9f2 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Fri, 20 Feb 2026 09:25:24 +0000
Subject: [PATCH 09/14] removed self handling of backend exception

---
 .../Python/cil/framework/data_container.py    | 45 +++----------------
 1 file changed, 5 insertions(+), 40 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index e8bf2aece5..d29cccd704 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -552,48 +552,13 @@ def pixel_wise_binary(self, pwop, x2, *args,  **kwargs):
             kwargs['out'] = outarr
                         
             if isinstance(x2, Number):
-                try:
                     pwop(self.as_array() , x2 , *args, **kwargs )
-                except Exception as bknerr:
-                    # check the size and dimension of out
-                    whatswrong = {}
-                    msg = ""
-                    if not self.check_dimensions(out):
-                        whatswrong['out shape'] = out.shape
-                        whatswrong['expected shape'] = self.shape
-                    if not out.dtype == self.dtype:
-                        whatswrong['out dtype'] = out.dtype
-                        whatswrong['expected dtype'] = self.dtype
-                    if not whatswrong == {}:
-                        # report to the user what's wrong
-                        msg = "Wrong size or type for data memory:\n"
-                        for k,v in whatswrong.items():
-                            msg += f"{k} {v}\n"
-                        print (msg)
-                    # report the error from the backend plus the error from the checks
-                    raise
             else:
-                try:
-                    if issubclass(x2.__class__ , DataContainer):
-                        pwop(self.as_array() , x2.as_array() , *args, **kwargs )
-                    else:
-                        pwop(self.as_array() , x2 , *args, **kwargs )
-                except Exception as bknerr:
-                    whatswrong = {}
-                    if not self.check_dimensions(x2):
-                        whatswrong['x2 shape'] = x2.shape
-                        whatswrong['expected shape'] = self.shape
-                    if not out.dtype == self.dtype:
-                        whatswrong['x2 dtype'] = x2.dtype
-                        whatswrong['expected dtype'] = self.dtype
-                    if not whatswrong == {}:
-                        # report to the user what's wrong
-                        msg = "Wrong size or type for data memory:\n"
-                        for k,v in whatswrong.items():
-                            msg += f"{k} {v}\n"
-                        print(msg)
-                    # report the error from the backend plus the error from the checks
-                    raise
+                if issubclass(x2.__class__ , DataContainer):
+                    pwop(self.as_array() , x2.as_array() , *args, **kwargs )
+                else:
+                    pwop(self.as_array() , x2 , *args, **kwargs )
+                
             out.fill(outarr)
             return out
 

From 2fa265ec6e97a3c979a3d3cb81b2cef246af4be3 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Fri, 20 Feb 2026 10:42:24 +0000
Subject: [PATCH 10/14] add test with parametrize

---
 Wrappers/Python/test/test_DataContainer.py | 58 +++++++++++++++++-----
 1 file changed, 45 insertions(+), 13 deletions(-)

diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index a9749d5c4d..99de63953e 100644
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -20,14 +20,26 @@
 import sys
 
 import numpy as np
+import numpy
+try:
+    import torch
+except ImportError:
+    torch = None
+try:
+    import cupy
+except ImportError:
+    cupy = None
 
 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
 
+from unittest_parametrize import ParametrizedTestCase, parametrize, param
+
 log = logging.getLogger(__name__)
 initialise_tests()
 
@@ -37,7 +49,7 @@ def aid(x):
     return x.as_array().__array_interface__['data'][0]
 
 
-class TestDataContainer(CCPiTestClass):
+class TestDataContainer(ParametrizedTestCase, CCPiTestClass):
     def create_DataContainer(self, X,Y,Z, value=1):
         a = value * np.ones((X, Y, Z), dtype='float32')
         #print("a refcount " , sys.getrefcount(a))
@@ -156,41 +168,61 @@ def test_ImageData_equal(self):
         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"),
+         param(cupy, None, None, None, id="cupy"),
+         ]) 
+    def testInlineAlgebra(self, xp, device, raise_error, err_type):
+        if xp is None:
+            self.skipTest(f"xp not available")
+
+        if xp == torch:
+            import array_api_compat.torch as xp
+        elif xp == cupy:
+            import array_api_compat.cupy as xp
+        elif xp == numpy:
+            import array_api_compat.numpy as xp
+        else:
+            raise ValueError(f"Unsupported xp: {xp}")
         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 += 2
         # self.assertEqual(ds.as_array()[0][0][0], 3.)
-        np.testing.assert_array_almost_equal(ds.as_array(), 3 * b)
+        
+        # np.testing.assert_array_almost_equal(ds.as_array(), 3 * b)
+        cil_allclose(ds.as_array(), 3 * b)
         ds -= 2
         # self.assertEqual(ds.as_array()[0][0][0], 1.)
-        np.testing.assert_array_almost_equal(ds.as_array(), b)
+        # np.testing.assert_array_almost_equal(ds.as_array(), b)
+        cil_allclose(ds.as_array(), b)
 
         ds *= 2
         # self.assertEqual(ds.as_array()[0][0][0], 2.)
-        np.testing.assert_array_almost_equal(ds.as_array(), b * 2)
+        cil_allclose(ds.as_array(), b * 2)
 
         ds /= 2
         # self.assertEqual(ds.as_array()[0][0][0], 1.)
-        np.testing.assert_array_almost_equal(ds.as_array(), b)
+        cil_allclose(ds.as_array(), b)
 
         ds1 = ds.copy()
         ds1 += 1
         ds += ds1
-        np.testing.assert_array_almost_equal(ds.as_array(), 3 * b)
+        cil_allclose(ds.as_array(), 3 * b)
         # self.assertEqual(ds.as_array()[0][0][0], 3.)
         ds -= ds1
-        np.testing.assert_array_almost_equal(ds.as_array(), b)
+        cil_allclose(ds.as_array(), b)
         # self.assertEqual(ds.as_array()[0][0][0], 1.)
         ds *= ds1
-        np.testing.assert_array_almost_equal(ds.as_array(), 2 * b)
+        cil_allclose(ds.as_array(), 2 * b)
         # self.assertEqual(ds.as_array()[0][0][0], 2.)
         ds /= ds1
-        np.testing.assert_array_almost_equal(ds.as_array(), b)
+        cil_allclose(ds.as_array(), b)
         # self.assertEqual(ds.as_array()[0][0][0], 1.)
 
+
     def test_unary_operations(self):
         X, Y, Z = 8, 16, 32
         a = -np.ones((X, Y, Z), dtype='float32')
@@ -1033,8 +1065,8 @@ def test_mean_direction(self):
         vd.dimension_labels = 'x'
         np.testing.assert_almost_equal(vd.mean(axis='x'), np.mean(vd))
 
-
     def test_multiply_out(self):
+        from cil.framework.array_api_compat import allclose as cil_allclose
         ig = ImageGeometry(10,11,12)
         u = ig.allocate()
         a = np.ones(u.shape)

From ace367e9abbb18ef3cf062043e6c573d85b0c617 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Sat, 21 Feb 2026 18:27:26 +0000
Subject: [PATCH 11/14] Update IndicatorBox

---
 .../Python/cil/optimisation/functions/IndicatorBox.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

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)

From 43b92497de4a059d34a1a13b6ac20744ebb0b780 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Sat, 21 Feb 2026 18:28:46 +0000
Subject: [PATCH 12/14] update get_slice

---
 Wrappers/Python/cil/framework/data_container.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index d29cccd704..f7ade972e6 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -180,7 +180,9 @@ def get_slice(self, **kw):
                 dimension_labels_list.remove(key)
                 if new_array is None:
                     new_array = self.as_array()
-                new_array = new_array.take(indices=value, axis=axis)
+                xp = array_api_compat.get_array_module(new_array)
+                # new_array = new_array.take(indices=value, axis=axis)
+                new_array = xp.take(new_array, indices=value, axis=axis)
 
         if new_array.ndim > 1:
             return DataContainer(new_array, False, dimension_labels_list)
@@ -266,6 +268,11 @@ def fill(self, array, **kwargs):
             required it should be passed directly as a kwarg.
             To fill random numbers using the earlier behaviour use `array='random_deprecated'` 
             or `array='random_int_deprecated'` 
+
+        dc.fill(some_data, vertical=1, horizontal_x=32)
+        will copy the data in `some_data` into the `dc` 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__
         
         Example
         -------
@@ -375,7 +382,6 @@ def fill(self, array, **kwargs):
         if dimension == {}:
             import warnings
             warnings.filterwarnings('error')
-            xp = array_namespace(self.as_array())
             self.array.__setitem__(slice(None, None, None), array)
             warnings.resetwarnings()
         else:

From a1eb356b8b7a16b6b3b4e6bbaf24a12214056a30 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Sat, 21 Feb 2026 19:12:39 +0000
Subject: [PATCH 13/14] add array0api-compat to requirements for GHA

---
 scripts/requirements-test-windows.yml | 1 +
 scripts/requirements-test.yml         | 1 +
 2 files changed, 2 insertions(+)

diff --git a/scripts/requirements-test-windows.yml b/scripts/requirements-test-windows.yml
index 2bfa84824f..64a32038dc 100644
--- a/scripts/requirements-test-windows.yml
+++ b/scripts/requirements-test-windows.yml
@@ -45,5 +45,6 @@ dependencies:
   - tqdm
   - zenodo_get >=1.6
   - pip
+  - array-api-compat
   - pip:
     - unittest-parametrize
diff --git a/scripts/requirements-test.yml b/scripts/requirements-test.yml
index a68116b2c2..d37428ecc7 100644
--- a/scripts/requirements-test.yml
+++ b/scripts/requirements-test.yml
@@ -46,5 +46,6 @@ dependencies:
   - tqdm
   - zenodo_get >=1.6
   - pip
+  - array-api-compat
   - pip:
     - unittest-parametrize

From f27d4ed42dd632eb65c7fa00b10348fd7f405006 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <14138589+paskino@users.noreply.github.com>
Date: Sat, 21 Feb 2026 23:09:53 +0000
Subject: [PATCH 14/14] restore code

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

diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py
index f7ade972e6..5090f48827 100644
--- a/Wrappers/Python/cil/framework/data_container.py
+++ b/Wrappers/Python/cil/framework/data_container.py
@@ -180,10 +180,8 @@ def get_slice(self, **kw):
                 dimension_labels_list.remove(key)
                 if new_array is None:
                     new_array = self.as_array()
-                xp = array_api_compat.get_array_module(new_array)
-                # new_array = new_array.take(indices=value, axis=axis)
-                new_array = xp.take(new_array, indices=value, axis=axis)
-
+                new_array = new_array.take(indices=value, axis=axis)
+                
         if new_array.ndim > 1:
             return DataContainer(new_array, False, dimension_labels_list)
         from .vector_data import VectorData