Merge pull request #188 from CUQI-DTU/sprint17_array_module

Move Data and CUQIarray classes to a separate module (array module)

cuqi/__init__.py

@@ -9,6 +9,7 @@
 from . import pde
 from . import problem
 from . import sampler
+from . import array
 from . import samples
 from . import solver
 from . import testproblem

cuqi/array/__init__.py

@@ -0,0 +1 @@
+from ._array import CUQIarray

cuqi/array/_array.py

@@ -0,0 +1,108 @@
+import numpy as np
+from cuqi.geometry import _DefaultGeometry
+
+
+class CUQIarray(np.ndarray):
+    """
+    A class to represent data arrays, subclassed from numpy array, along with geometry and plotting
+
+    Parameters
+    ----------
+    input_array : ndarray
+        A numpy array holding the parameter or function values. 
+    
+    is_par : bool, default True
+        Boolean flag whether input_array is to be interpreted as parameter (True) or function values (False).
+
+    geometry : cuqi.geometry.Geometry, default None
+        Contains the geometry related of the data
+
+    Attributes
+    ----------
+    funvals : CUQIarray
+        Returns itself as function values.
+
+    parameters : CUQIarray
+        Returns itself as parameters.
+
+    Methods
+    ----------
+    :meth:`plot`: Plots the data as function or parameters.
+    """
+
+    def __repr__(self) -> str: 
+        return "CUQIarray: NumPy array wrapped with geometry.\n" + \
+               "---------------------------------------------\n\n" + \
+            "Geometry:\n {}\n\n".format(self.geometry) + \
+            "Parameters:\n {}\n\n".format(self.is_par) + \
+            "Array:\n" + \
+            super().__repr__()
+
+    def __new__(cls, input_array, is_par=True, geometry=None):
+        # Input array is an already formed ndarray instance
+        # We first cast to be our class type
+        obj = np.asarray(input_array).view(cls)
+        # add the new attribute to the created instance
+        obj.is_par = is_par
+        if (not is_par) and (geometry is None):
+            raise ValueError("geometry cannot be none when initializing a CUQIarray as function values (with is_par False).")
+        if is_par and (obj.ndim>1):
+            raise ValueError("input_array cannot be multidimensional when initializing CUQIarray as parameter (with is_par True).")
+        if geometry is None:
+            geometry = _DefaultGeometry(grid=obj.__len__())
+        obj.geometry = geometry
+        # Finally, we must return the newly created object:
+        return obj
+
+    def __array_finalize__(self, obj):
+        # see InfoArray.__array_finalize__ for comments
+        if obj is None: return
+        self.is_par = getattr(obj, 'is_par', True)
+        self.geometry = getattr(obj, 'geometry', None)
+
+    @property
+    def funvals(self):
+        if self.is_par is True:
+            vals = self.geometry.par2fun(self)
+        else:
+            vals = self
+
+        if isinstance(vals, np.ndarray):
+            if vals.dtype == np.dtype('O'):
+                # if vals is of type np.ndarray, but the data type of the array
+                # is object (e.g. FEniCS function), then extract the object and
+                # return it. reshape(1) is needed to convert the shape from
+                # () to (1,).
+                return self.reshape(1)[0]
+            else:
+                # else, cast the np.ndarray to a CUQIarray
+                return type(self)(vals,is_par=False,geometry=self.geometry) #vals.view(np.ndarray)
+        else:
+            return vals 
+
+    @property
+    def parameters(self):
+        if self.is_par is False:
+            vals = self.geometry.fun2par(self)
+        else:
+            vals = self
+        return type(self)(vals,is_par=True,geometry=self.geometry)
+
+    def to_numpy(self):
+        """Return a numpy array of the CUQIarray data. If is_par is True, then 
+        the parameters are returned as numpy.ndarray. If is_par is False, then 
+        the function values are returned instead.
+        """
+        try:
+            return self.view(np.ndarray)
+        except:
+            raise ValueError(
+                f"Cannot convert {self.__class__.__name__} to numpy array")
+
+    def plot(self, plot_par=False, **kwargs):
+        if plot_par:
+            kwargs["is_par"]=True
+            return self.geometry.plot(self.parameters, plot_par=plot_par, **kwargs)
+        else:
+            kwargs["is_par"]=False
+            return self.geometry.plot(self.funvals, **kwargs)

cuqi/distribution/_distribution.py

@@ -5,7 +5,8 @@
 from functools import partial
 from cuqi.density import Density, EvaluatedDensity
 from cuqi.likelihood import Likelihood
-from cuqi.samples import Samples, CUQIarray
+from cuqi.samples import Samples
+from cuqi.array import CUQIarray
 from cuqi.geometry import _DefaultGeometry, Geometry
 from cuqi.utilities import infer_len, get_writeable_attributes, get_writeable_properties, get_non_default_args, get_indirect_variables
 import numpy as np # To be replaced by cuqi.array_api

cuqi/likelihood/_likelihood.py

@@ -21,7 +21,7 @@ class Likelihood(Density):
     distribution : ~cuqi.distribution.Distribution
         | Distribution to create likelihood from.
     
-    data : ~cuqi.samples.CUQIarray or array_like
+    data : ~cuqi.array.CUQIarray or array_like
         | Observation to create likelihood from.
 
     """

cuqi/model/_model.py

@@ -2,7 +2,8 @@
 from scipy.sparse import csc_matrix
 from scipy.sparse import hstack
 from scipy.linalg import solve
-from cuqi.samples import Samples, CUQIarray
+from cuqi.samples import Samples
+from cuqi.array import CUQIarray
 from cuqi.geometry import Geometry, _DefaultGeometry, _get_identity_geometries
 import cuqi
 import matplotlib.pyplot as plt
@@ -93,7 +94,7 @@ def _input2fun(self, x, geometry, is_par):
 
         Parameters
         ----------
-        x : ndarray or cuqi.samples.CUQIarray
+        x : ndarray or cuqi.array.CUQIarray
             The input value to be converted.
 
         geometry : cuqi.geometry.Geometry
@@ -105,7 +106,7 @@ def _input2fun(self, x, geometry, is_par):
 
         Returns
         -------
-        ndarray or cuqi.samples.CUQIarray
+        ndarray or cuqi.array.CUQIarray
             The input value represented as a function.
         """
         if type(x) is CUQIarray and not isinstance(x.geometry, _DefaultGeometry):
@@ -120,18 +121,18 @@ def _output2par(self, out, geometry, to_CUQIarray=False):
 
         Parameters
         ----------
-        out : ndarray or cuqi.samples.CUQIarray
+        out : ndarray or cuqi.array.CUQIarray
             The output value to be converted.
 
         geometry : cuqi.geometry.Geometry
             The geometry representing the argument `out`.
 
         to_CUQIarray : bool
-            If True, the output is wrapped as a cuqi.samples.CUQIarray.
+            If True, the output is wrapped as a cuqi.array.CUQIarray.
 
         Returns
         -------
-        ndarray or cuqi.samples.CUQIarray
+        ndarray or cuqi.array.CUQIarray
             The output value represented as parameters.
         """ 
         out = geometry.fun2par(out)
@@ -156,7 +157,7 @@ def _apply_func(self, func, func_range_geometry, func_domain_geometry, x, is_par
         func_domain_geometry : cuqi.geometry.Geometry
             The geometry representing the function `func` domain.
 
-        x : ndarray or cuqi.samples.CUQIarray
+        x : ndarray or cuqi.array.CUQIarray
             The input value to the operator.
 
         is_par : bool
@@ -165,7 +166,7 @@ def _apply_func(self, func, func_range_geometry, func_domain_geometry, x, is_par
 
         Returns
         -------
-        ndarray or cuqi.samples.CUQIarray
+        ndarray or cuqi.array.CUQIarray
             The output of the function `func` converted to parameters.
         """ 
         # If input x is Samples we apply func for each sample
@@ -211,7 +212,7 @@ def forward(self, *args, is_par=True, **kwargs):
 
         Parameters
         ----------
-        *args : ndarray or cuqi.samples.CUQIarray
+        *args : ndarray or cuqi.array.CUQIarray
             The model input.
 
         is_par : bool
@@ -223,7 +224,7 @@ def forward(self, *args, is_par=True, **kwargs):
 
         Returns
         -------
-        ndarray or cuqi.samples.CUQIarray
+        ndarray or cuqi.array.CUQIarray
             The model output. Always returned as parameters.
         """
 
@@ -426,12 +427,12 @@ def adjoint(self, y, is_par=True):
 
         Parameters
         ----------
-        y : ndarray or cuqi.samples.CUQIarray
+        y : ndarray or cuqi.array.CUQIarray
             The adjoint model input.
 
         Returns
         -------
-        ndarray or cuqi.samples.CUQIarray
+        ndarray or cuqi.array.CUQIarray
             The adjoint model output. Always returned as parameters.
         """
         return self._apply_func(self._adjoint_func,

cuqi/problem/_problem.py

@@ -11,7 +11,7 @@
 from cuqi.likelihood import Likelihood
 from cuqi.geometry import _DefaultGeometry
 from cuqi.utilities import ProblemInfo
-from cuqi.samples import CUQIarray
+from cuqi.array import CUQIarray
 
 from copy import copy
 
@@ -228,7 +228,7 @@ def ML(self, disp=True, x0=None) -> CUQIarray:
         x_ML, solver_info = self._solve_max_point(self.likelihood, disp=disp, x0=x0)
 
         # Wrap the result in a CUQIarray and add solver info
-        x_ML = cuqi.samples.CUQIarray(x_ML, geometry=self.likelihood.geometry)
+        x_ML = cuqi.array.CUQIarray(x_ML, geometry=self.likelihood.geometry)
         x_ML.info = solver_info
 
         return x_ML
@@ -285,7 +285,7 @@ def MAP(self, disp=True, x0=None) -> CUQIarray:
             x_MAP, solver_info = self._solve_max_point(self.posterior, disp=disp, x0=x0)
 
         # Wrap the result in a CUQIarray and add solver info
-        x_MAP = cuqi.samples.CUQIarray(x_MAP, geometry=self.posterior.geometry)
+        x_MAP = cuqi.array.CUQIarray(x_MAP, geometry=self.posterior.geometry)
         x_MAP.info = solver_info
         return x_MAP
 

cuqi/samples/__init__.py

@@ -1 +1 @@
-from ._samples import CUQIarray, Samples, Data
+from ._samples import Samples