fix: ✨ adding interpolator factory

Author: lachlangrose

LoopStructural/interpolators/__init__.py

@@ -78,3 +78,14 @@ class InterpolatorType(IntEnum):
 from ._p1interpolator import P1Interpolator
 from ._p2interpolator import P2Interpolator
 from ._builders import get_interpolator
+
+interpolator_map = {
+    InterpolatorType.BASE: GeologicalInterpolator,
+    InterpolatorType.BASE_DISCRETE: DiscreteInterpolator,
+    InterpolatorType.FINITE_DIFFERENCE: FiniteDifferenceInterpolator,
+    InterpolatorType.DISCRETE_FOLD: DiscreteFoldInterpolator,
+    InterpolatorType.PIECEWISE_LINEAR: PiecewiseLinearInterpolator,
+    InterpolatorType.PIECEWISE_QUADRATIC: PiecewiseLinearInterpolator,
+    InterpolatorType.BASE_DATA_SUPPORTED: GeologicalInterpolator,
+    # InterpolatorType.SURFE: SurfeRBFInterpolator,
+}

LoopStructural/interpolators/_discrete_fold_interpolator.py

@@ -144,9 +144,7 @@ def add_fold_constraints(
             """
             np.random.shuffle(element_idx)
 
-            logger.info(
-                f"Adding fold orientation constraint to {self.propertyname} w = {fold_orientation}"
-            )
+            logger.info(f"Adding fold orientation constraint to w = {fold_orientation}")
             A = np.einsum(
                 "ij,ijk->ik",
                 deformed_orientation[element_idx[::step], :],
@@ -165,9 +163,7 @@ def add_fold_constraints(
             """
             np.random.shuffle(element_idx)
 
-            logger.info(
-                f"Adding fold axis constraint to {self.propertyname} w = {fold_axis_w}"
-            )
+            logger.info(f"Adding fold axis constraint to  w = {fold_axis_w}")
             A = np.einsum(
                 "ij,ijk->ik",
                 fold_axis[element_idx[::step], :],
@@ -188,7 +184,7 @@ def add_fold_constraints(
             np.random.shuffle(element_idx)
 
             logger.info(
-                f"Adding fold normalisation constraint to {self.propertyname} w = {fold_normalisation}"
+                f"Adding fold normalisation constraint to  w = {fold_normalisation}"
             )
             A = np.einsum(
                 "ij,ijk->ik", dgz[element_idx[::step], :], eg[element_idx[::step], :, :]
@@ -212,7 +208,7 @@ def add_fold_constraints(
             fold constant gradient
             """
             logger.info(
-                f"Adding fold regularisation constraint to {self.propertyname} w = {fold_regularisation[0]} {fold_regularisation[1]} {fold_regularisation[2]}"
+                f"Adding fold regularisation constraint to  w = {fold_regularisation[0]} {fold_regularisation[1]} {fold_regularisation[2]}"
             )
 
             idc, c, ncons = fold_cg(

LoopStructural/interpolators/_discrete_interpolator.py

@@ -21,7 +21,7 @@
 class DiscreteInterpolator(GeologicalInterpolator):
     """ """
 
-    def __init__(self, support):
+    def __init__(self, support, data={}, c=None, up_to_date=False):
         """
         Base class for a discrete interpolator e.g. piecewise linear or finite difference which is
         any interpolator that solves the system using least squares approximation
@@ -31,9 +31,14 @@ def __init__(self, support):
         support
             A discrete mesh with, nodes, elements, etc
         """
-        GeologicalInterpolator.__init__(self)
+        GeologicalInterpolator.__init__(self, data=data, up_to_date=up_to_date)
         self.B = []
         self.support = support
+        self.c = (
+            np.array(c)
+            if c is not None and np.array(c).shape[0] == self.support.n_nodes
+            else np.zeros(self.support.n_nodes)
+        )
         self.region_function = lambda xyz: np.ones(xyz.shape[0], dtype=bool)
         # self.region_map[self.region] = np.array(range(0,
         # len(self.region_map[self.region])))
@@ -95,21 +100,6 @@ def region_map(self):
         region_map[self.region] = np.array(range(0, len(region_map[self.region])))
         return region_map
 
-    def set_property_name(self, propertyname):
-        """
-        Set the property name attribute, this is usually used to
-        save the property on the support
-
-        Parameters
-        ----------
-        propertyname
-
-        Returns
-        -------
-
-        """
-        self.propertyname = propertyname
-
     def set_region(self, region=None):
         """
         Set the region of the support the interpolator is working on
@@ -724,7 +714,7 @@ def _solve_pyamg(self, A, B, tol=1e-12, x0=None, verb=False, **kwargs):
         logger.info("Solving using pyamg: tol {}".format(tol))
         return pyamg.solve(A, B, tol=tol, x0=x0, verb=verb)[: self.nx]
 
-    def _solve(self, solver="cg", **kwargs):
+    def solve_system(self, solver="cg", **kwargs):
         """
         Main entry point to run the solver and update the node value
         attribute for the
@@ -784,24 +774,18 @@ def _solve(self, solver="cg", **kwargs):
                 P, A, q, l, u, mkl=kwargs.get("mkl", False)
             )  # , **kwargs)
         # check solution is not nan
-        # self.support.properties[self.propertyname] = self.c
         if np.all(self.c == np.nan):
             self.valid = False
             logger.warning("Solver not run, no scalar field")
             return
         # if solution is all 0, probably didn't work
         if np.all(self.c[self.region] == 0):
             self.valid = False
-            logger.warning(
-                "No solution, {} scalar field 0. Add more data.".format(
-                    self.propertyname
-                )
-            )
+            logger.warning("No solution, scalar field 0. Add more data.")
+
             return
         self.valid = True
-        logging.info(
-            f"Solving interpolation: {self.propertyname} took: {time()-starttime}"
-        )
+        logging.info(f"Solving interpolation took: {time()-starttime}")
         self.up_to_date = True
 
     def update(self):
@@ -822,7 +806,7 @@ def update(self):
             return False
         if not self.up_to_date:
             self.setup_interpolator()
-            return self._solve(self.solver)
+            return self.solve_system(self.solver)
 
     def evaluate_value(self, evaluation_points: np.ndarray) -> np.ndarray:
         """Evaluate the value of the interpolator at location
@@ -864,3 +848,12 @@ def evaluate_gradient(self, evaluation_points: np.ndarray) -> np.ndarray:
         if evaluation_points.shape[0] > 0:
             return self.support.evaluate_gradient(evaluation_points, self.c)
         return np.zeros((0, 3))
+
+    def to_dict(self):
+        return {
+            "type": self.type.name,
+            "support": self.support.to_dict(),
+            "c": self.c,
+            **super().to_dict(),
+            # 'region_function':self.region_function,
+        }

LoopStructural/interpolators/_finite_difference_interpolator.py

@@ -17,7 +17,7 @@
 
 
 class FiniteDifferenceInterpolator(DiscreteInterpolator):
-    def __init__(self, grid):
+    def __init__(self, grid, data={}):
         """
         Finite difference interpolation on a regular cartesian grid
 
@@ -26,7 +26,7 @@ def __init__(self, grid):
         grid : StructuredGrid
         """
         self.shape = "rectangular"
-        DiscreteInterpolator.__init__(self, grid)
+        DiscreteInterpolator.__init__(self, grid, data=data)
         # default weights for the interpolation matrix are 1 in x,y,z and
         # 1/
         self.set_interpolation_weights(
@@ -52,7 +52,7 @@ def __init__(self, grid):
         # grid.step_vector[2]
         self.type = InterpolatorType.FINITE_DIFFERENCE
 
-    def _setup_interpolator(self, **kwargs):
+    def setup_interpolator(self, **kwargs):
         """
 
         Parameters
@@ -95,7 +95,6 @@ def _setup_interpolator(self, **kwargs):
                 self.assemble_inner(o[0], o[1])
         # otherwise just use defaults
         if "operators" not in kwargs:
-
             operator = Operator.Dxy_mask
             weight = (
                 self.interpolation_weights["dxy"] / 4
@@ -179,7 +178,7 @@ def add_vaue_constraints(self, w=1.0):
             )
             if np.sum(inside) <= 0:
                 logger.warning(
-                    f"{self.propertyname}: {np.sum(~inside)} \
+                    f"{np.sum(~inside)} \
                         value constraints not added: outside of model bounding box"
                 )
 
@@ -218,7 +217,7 @@ def add_inequality_constraints(self, w=1.0):
             )
             if np.sum(inside) <= 0:
                 logger.warning(
-                    f"{self.propertyname}: {np.sum(~inside)} \
+                    f"{np.sum(~inside)} \
                         value constraints not added: outside of model bounding box"
                 )
 
@@ -335,7 +334,7 @@ def add_gradient_constraints(self, w=1.0):
             )
             if np.sum(inside) <= 0:
                 logger.warning(
-                    f"{self.propertyname}: {np.sum(~inside)} \
+                    f" {np.sum(~inside)} \
                         norm constraints not added: outside of model bounding box"
                 )
 
@@ -400,7 +399,7 @@ def add_norm_constraints(self, w=1.0):
             )
             if np.sum(inside) <= 0:
                 logger.warning(
-                    f"{self.propertyname}: {np.sum(~inside)} \
+                    f"{np.sum(~inside)} \
                         norm constraints not added: outside of model bounding box"
                 )
             self.up_to_date = False
@@ -459,7 +458,7 @@ def add_gradient_orthogonal_constraints(
             )
             if np.sum(inside) <= 0:
                 logger.warning(
-                    f"{self.propertyname}: {np.sum(~inside)} \
+                    f"{np.sum(~inside)} \
                         gradient constraints not added: outside of model bounding box"
                 )
             self.up_to_date = False

LoopStructural/interpolators/_geological_interpolator.py

@@ -1,6 +1,8 @@
 """
 Base geological interpolator
 """
+from abc import ABC, ABCMeta, abstractmethod
+from LoopStructural.utils.exceptions import LoopTypeError
 from ..interpolators import InterpolatorType
 import numpy as np
 
@@ -9,22 +11,24 @@
 logger = getLogger(__name__)
 
 
-class GeologicalInterpolator:
+class GeologicalInterpolator(metaclass=ABCMeta):
     """
     Attributes
     ----------
     data : dict
         a dictionary with np.arrays for gradient, value, normal, tangent data
     """
 
-    def __init__(self):
+    @abstractmethod
+    def __init__(self, data={}, up_to_date=False):
         """
         This class is the base class for a geological interpolator and contains
         all of the main interface functions. Any class that is inheriting from
         this should be callable by using any of these functions. This will
         enable interpolators to be interchanged.
         """
-        self.data = {}  # None
+        self._data = {}
+        self.data = data  # None
         self.clean()  # init data structure
 
         self.n_g = 0
@@ -33,12 +37,22 @@ def __init__(self):
         self.n_t = 0
 
         self.type = InterpolatorType.BASE
-        self.up_to_date = False
+        self.up_to_date = up_to_date
         self.constraints = []
-        self.propertyname = "defaultproperty"
         self.__str = "Base Geological Interpolator"
         self.valid = False
 
+    @property
+    def data(self):
+        return self._data
+
+    @data.setter
+    def data(self, data):
+        if data is None:
+            data = {}
+        for k, v in data.items():
+            self._data[k] = np.array(v)
+
     def __str__(self):
         name = f"{self.type} \n"
         name += f"{self.n_g} gradient points\n"
@@ -48,23 +62,16 @@ def __str__(self):
         name += f"{self.n_g + self.n_i + self.n_n + self.n_t} total points\n"
         return name
 
+    def check_array(self, array: np.ndarray):
+        try:
+            return np.array(array)
+        except Exception as e:
+            raise LoopTypeError(str(e))
+
+    @abstractmethod
     def set_region(self, **kwargs):
         pass
 
-    def set_property_name(self, name):
-        """
-        Set the name of the interpolated property
-        Parameters
-        ----------
-        name : string
-            name of the property to be saved on a mesh
-
-        Returns
-        -------
-
-        """
-        self.propertyname = name
-
     def set_value_constraints(self, points: np.ndarray):
         """
 
@@ -78,6 +85,7 @@ def set_value_constraints(self, points: np.ndarray):
         -------
 
         """
+        points = self.check_array(points)
         if points.shape[1] < 4:
             raise ValueError("Value points must at least have X,Y,Z,val")
         self.data["value"] = points
@@ -208,34 +216,48 @@ def get_interface_constraints(self):
     def get_inequality_constraints(self):
         return self.data["inequality"]
 
+    # @abstractmethod
     def setup(self, **kwargs):
         """
         Runs all of the required setting up stuff
         """
-        self._setup_interpolator(**kwargs)
+        self.setup_interpolator(**kwargs)
 
+    @abstractmethod
     def setup_interpolator(self, **kwargs):
         """
         Runs all of the required setting up stuff
         """
-        self._setup_interpolator(**kwargs)
+        self.setup_interpolator(**kwargs)
 
+    @abstractmethod
     def solve_system(self, **kwargs):
         """
         Solves the interpolation equations
         """
         self._solve(**kwargs)
         self.up_to_date = True
 
+    @abstractmethod
     def update(self):
         return False
 
+    @abstractmethod
     def evaluate_value(self, locations: np.ndarray):
         raise NotImplementedError("evaluate_value not implemented")
 
+    @abstractmethod
     def evaluate_gradient(self, locations: np.ndarray):
         raise NotImplementedError("evaluate_gradient not implemented")
 
+    def to_dict(self):
+        return {
+            "type": self.type,
+            "data": self.data,
+            "up_to_date": self.up_to_date,
+            "valid": self.valid,
+        }
+
     def clean(self):
         """
         Removes all of the data from an interpolator