Interpolator API (2D fix, inequalities, surfe)

* docs: updating readme

* fix: adding plot method to loop objects

* fix: default solver to cg

* fix: adding inequality pairs to interpolator api

* fix: linting

* adding fit and eval gradient

* move default operators to the support for 2d/3d

* remove property wraper from structured tertra vtk

* add warning if ineq points are outside support

* pass inequality to loop interpolator

* add plot option for interpolator api

* make weights all 1. and scale in operator creation

* incorrect slicing fixed

* check constraints shape using interpolator dimensions

* remove non linear interpolator for now

* set bb step vector/nelements for 2d

* change neighbour indexing to order by global index.

This makes the matrix diagonal.

* updating fdi for 2d interpolation using support dimension to infer the dimension

* is row scaling the interpolation matrix sensible?

* updating interpolation api to use inequalities pairs/values

* fix: adding surfe to interpolator api

* fix: change argument for p1 interpolator

* fix: pass constraint name to add to matrix

* Inequalities are causing seg fault

* pass entire vector for strike dip vector

LoopStructural/datatypes/_bounding_box.py

@@ -148,7 +148,13 @@ def nelements(self, nelements: Union[int, float]):
         ele_vol = box_vol / nelements
         # calculate the step vector of a regular cube
         step_vector = np.zeros(self.dimensions)
-        step_vector[:] = ele_vol ** (1.0 / 3.0)
+        if self.dimensions == 2:
+            step_vector[:] = ele_vol ** (1.0 / 2.0)
+        elif self.dimensions == 3:
+            step_vector[:] = ele_vol ** (1.0 / 3.0)
+        else:
+            logger.warning("Can only set nelements for 2d or 3D bounding box")
+            return
         # number of steps is the length of the box / step vector
         nsteps = np.ceil((self.maximum - self.origin) / step_vector).astype(int)
         self.nsteps = nsteps
@@ -270,7 +276,10 @@ def with_buffer(self, buffer: float = 0.2) -> BoundingBox:
         origin = self.origin - buffer * (self.maximum - self.origin)
         maximum = self.maximum + buffer * (self.maximum - self.origin)
         return BoundingBox(
-            origin=origin, maximum=maximum, global_origin=self.global_origin + origin
+            origin=origin,
+            maximum=maximum,
+            global_origin=self.global_origin + origin,
+            dimensions=self.dimensions,
         )
 
     def get_value(self, name):

LoopStructural/datatypes/_point.py

@@ -3,6 +3,9 @@
 
 from typing import Optional, Union
 import io
+from LoopStructural.utils import getLogger
+
+logger = getLogger(__name__)
 
 
 @dataclass
@@ -29,6 +32,20 @@ def vtk(self):
         points["values"] = self.values
         return points
 
+    def plot(self, pyvista_kwargs={}):
+        """Calls pyvista plot on the vtk object
+
+        Parameters
+        ----------
+        pyvista_kwargs : dict, optional
+            kwargs passed to pyvista.DataSet.plot(), by default {}
+        """
+        try:
+            self.vtk().plot(**pyvista_kwargs)
+            return
+        except ImportError:
+            logger.error("pyvista is required for vtk")
+
     def save(self, filename: Union[str, io.StringIO], ext=None):
         if isinstance(filename, io.StringIO):
             if ext is None:
@@ -123,6 +140,20 @@ def vtk(self, geom='arrow', scale=1.0, scale_function=None, tolerance=0.05):
         # Perform the glyph
         return points.glyph(orient="vectors", geom=geom, tolerance=tolerance)
 
+    def plot(self, pyvista_kwargs={}):
+        """Calls pyvista plot on the vtk object
+
+        Parameters
+        ----------
+        pyvista_kwargs : dict, optional
+            kwargs passed to pyvista.DataSet.plot(), by default {}
+        """
+        try:
+            self.vtk().plot(**pyvista_kwargs)
+            return
+        except ImportError:
+            logger.error("pyvista is required for vtk")
+
     def save(self, filename):
         filename = str(filename)
         ext = filename.split('.')[-1]

LoopStructural/datatypes/_structured_grid.py

@@ -1,6 +1,9 @@
 from typing import Dict
 import numpy as np
 from dataclasses import dataclass
+from LoopStructural.utils import getLogger
+
+logger = getLogger(__name__)
 
 
 @dataclass
@@ -46,6 +49,20 @@ def vtk(self):
             grid.cell_data[name] = data.flatten(order="F")
         return grid
 
+    def plot(self, pyvista_kwargs={}):
+        """Calls pyvista plot on the vtk object
+
+        Parameters
+        ----------
+        pyvista_kwargs : dict, optional
+            kwargs passed to pyvista.DataSet.plot(), by default {}
+        """
+        try:
+            self.vtk().plot(**pyvista_kwargs)
+            return
+        except ImportError:
+            logger.error("pyvista is required for vtk")
+
     def merge(self, other):
         if not np.all(np.isclose(self.origin, other.origin)):
             raise ValueError("Origin of grids must be the same")

LoopStructural/datatypes/_surface.py

@@ -2,6 +2,9 @@
 from typing import Optional
 import numpy as np
 import io
+from LoopStructural.utils import getLogger
+
+logger = getLogger(__name__)
 
 
 @dataclass
@@ -85,6 +88,20 @@ def vtk(self):
                 surface.cell_data[k] = np.array(v)
         return surface
 
+    def plot(self, pyvista_kwargs={}):
+        """Calls pyvista plot on the vtk object
+
+        Parameters
+        ----------
+        pyvista_kwargs : dict, optional
+            kwargs passed to pyvista.DataSet.plot(), by default {}
+        """
+        try:
+            self.vtk().plot(**pyvista_kwargs)
+            return
+        except ImportError:
+            logger.error("pyvista is required for vtk")
+
     def to_dict(self, flatten=False):
         triangles = self.triangles
         vertices = self.vertices

LoopStructural/interpolators/__init__.py

@@ -52,6 +52,7 @@ class InterpolatorType(IntEnum):
     "P2": InterpolatorType.PIECEWISE_QUADRATIC,
     "P1": InterpolatorType.PIECEWISE_LINEAR,
     "DFI": InterpolatorType.DISCRETE_FOLD,
+    'surfe': InterpolatorType.SURFE,
 }
 from ..interpolators._geological_interpolator import GeologicalInterpolator
 from ..interpolators._discrete_interpolator import DiscreteInterpolator
@@ -79,6 +80,11 @@ class InterpolatorType(IntEnum):
 from ..interpolators._p2interpolator import P2Interpolator
 from ..interpolators._p1interpolator import P1Interpolator
 
+try:
+    from ..interpolators._surfe_wrapper import SurfeRBFInterpolator
+except ImportError:
+    logger.warning("Surfe is not installed, SurfeRBFInterpolator will not be available")
+    SurfeRBFInterpolator = None
 interpolator_map = {
     InterpolatorType.BASE: GeologicalInterpolator,
     InterpolatorType.BASE_DISCRETE: DiscreteInterpolator,
@@ -87,6 +93,7 @@ class InterpolatorType(IntEnum):
     InterpolatorType.PIECEWISE_LINEAR: P1Interpolator,
     InterpolatorType.PIECEWISE_QUADRATIC: P2Interpolator,
     InterpolatorType.BASE_DATA_SUPPORTED: GeologicalInterpolator,
+    InterpolatorType.SURFE: SurfeRBFInterpolator,
 }
 
 support_interpolator_map = {
@@ -100,6 +107,12 @@ class InterpolatorType(IntEnum):
         3: SupportType.P2UnstructuredTetMesh,
         2: SupportType.P2Unstructured2d,
     },
+    InterpolatorType.SURFE: {
+        3: SupportType.DataSupported,
+        2: SupportType.DataSupported,
+    },
 }
 
 from ._interpolator_factory import InterpolatorFactory
+
+# from ._api import LoopInterpolator

LoopStructural/interpolators/_api.py

@@ -19,6 +19,8 @@ def __init__(
         dimensions: int = 3,
         type=InterpolatorType.FINITE_DIFFERENCE,
         nelements: int = 1000,
+        interpolator_setup_kwargs={},
+        buffer: float = 0.2,
     ):
         """Scikitlearn like interface for LoopStructural interpolators
         useful for quickly building an interpolator to apply to a dataset
@@ -37,21 +39,22 @@ def __init__(
         nelements : int, optional
             degrees of freedom for interpolator, by default 1000
         """
+        logger.warning("LoopInterpolator is experimental and the API is subject to change")
         self.dimensions = dimensions
         self.type = "FDI"
         self.bounding_box = bounding_box
         self.interpolator: GeologicalInterpolator = InterpolatorFactory.create_interpolator(
-            type,
-            bounding_box,
-            nelements,
+            type, bounding_box, nelements, buffer=buffer
         )
+        self.interpolator_setup_kwargs = interpolator_setup_kwargs
 
     def fit(
         self,
         values: Optional[np.ndarray] = None,
         tangent_vectors: Optional[np.ndarray] = None,
         normal_vectors: Optional[np.ndarray] = None,
-        inequality_constraints: Optional[np.ndarray] = None,
+        inequality_value_constraints: Optional[np.ndarray] = None,
+        inequality_pairs_constraints: Optional[np.ndarray] = None,
     ):
         """_summary_
 
@@ -66,16 +69,18 @@ def fit(
         inequality_constraints : Optional[np.ndarray], optional
             _description_, by default None
         """
+
         if values is not None:
             self.interpolator.set_value_constraints(values)
         if tangent_vectors is not None:
             self.interpolator.set_tangent_constraints(tangent_vectors)
         if normal_vectors is not None:
             self.interpolator.set_normal_constraints(normal_vectors)
-        if inequality_constraints:
-            pass
-
-        self.interpolator.setup()
+        if inequality_value_constraints is not None:
+            self.interpolator.set_value_inequality_constraints(inequality_value_constraints)
+        if inequality_pairs_constraints is not None:
+            self.interpolator.set_inequality_pairs_constraints(inequality_pairs_constraints)
+        self.interpolator.setup(**self.interpolator_setup_kwargs)
 
     def evaluate_scalar_value(self, locations: np.ndarray) -> np.ndarray:
         """Evaluate the value of the interpolator at locations
@@ -114,30 +119,93 @@ def fit_and_evaluate_value(
         values: Optional[np.ndarray] = None,
         tangent_vectors: Optional[np.ndarray] = None,
         normal_vectors: Optional[np.ndarray] = None,
-        inequality_constraints: Optional[np.ndarray] = None,
+        inequality_value_constraints: Optional[np.ndarray] = None,
+        inequality_pairs_constraints: Optional[np.ndarray] = None,
     ):
         # get locations
         self.fit(
             values=values,
             tangent_vectors=tangent_vectors,
             normal_vectors=normal_vectors,
-            inequality_constraints=inequality_constraints,
+            inequality_value_constraints=inequality_value_constraints,
+            inequality_pairs_constraints=inequality_pairs_constraints,
         )
         locations = self.interpolator.get_data_locations()
-        return self.evalute_scalar_value(locations)
+        return self.evaluate_scalar_value(locations)
 
     def fit_and_evaluate_gradient(
         self,
         values: Optional[np.ndarray] = None,
         tangent_vectors: Optional[np.ndarray] = None,
         normal_vectors: Optional[np.ndarray] = None,
-        inequality_constraints: Optional[np.ndarray] = None,
+        inequality_value_constraints: Optional[np.ndarray] = None,
+        inequality_pairs_constraints: Optional[np.ndarray] = None,
     ):
         self.fit(
             values=values,
             tangent_vectors=tangent_vectors,
             normal_vectors=normal_vectors,
-            inequality_constraints=inequality_constraints,
+            inequality_value_constraints=inequality_value_constraints,
+            inequality_pairs_constraints=inequality_pairs_constraints,
         )
         locations = self.interpolator.get_data_locations()
         return self.evaluate_gradient(locations)
+
+    def fit_and_evaluate_value_and_gradient(
+        self,
+        values: Optional[np.ndarray] = None,
+        tangent_vectors: Optional[np.ndarray] = None,
+        normal_vectors: Optional[np.ndarray] = None,
+        inequality_value_constraints: Optional[np.ndarray] = None,
+        inequality_pairs_constraints: Optional[np.ndarray] = None,
+    ):
+        self.fit(
+            values=values,
+            tangent_vectors=tangent_vectors,
+            normal_vectors=normal_vectors,
+            inequality_value_constraints=inequality_value_constraints,
+            inequality_pairs_constraints=inequality_pairs_constraints,
+        )
+        locations = self.interpolator.get_data_locations()
+        return self.evaluate_scalar_value(locations), self.evaluate_gradient(locations)
+
+    def plot(self, ax=None, **kwargs):
+        """Plots a 2d map scalar field or 3d pyvista plot
+
+        Parameters
+        ----------
+        ax : matplotlib axes, optional
+            The axes you want to add the plot to, otherwise it will make one, by default None
+
+        Returns
+        -------
+        _type_
+            _description_
+        """
+        if self.dimensions == 3:
+            vtkgrid = self.interpolator.support.vtk()
+            vtkgrid['val'] = self.interpolator.c
+            vtkgrid.plot(**kwargs)
+            return vtkgrid
+        elif self.dimensions == 2:
+            if ax is None:
+                import matplotlib.pyplot as plt
+
+                fig, ax = plt.subplots()
+            val = self.interpolator.c
+            val = np.rot90(val.reshape(self.interpolator.support.nsteps, order='F'), 3)
+            ax.imshow(
+                val,
+                origin='lower',
+                extent=[
+                    self.bounding_box.origin[0],
+                    self.bounding_box.maximum[0],
+                    self.bounding_box.origin[1],
+                    self.bounding_box.maximum[1],
+                ],
+                **kwargs,
+            )
+            return val, ax
+
+    # def isovalue(self, value: float, **kwargs):
+    #     self.interpolator.isovalue(value, **kwargs)