fix: adding inequality pairs

any data with 'pair_id' is interpreted as an inequality pair data point.
Inequality pairs are added with a step threshold

Author: lachlangrose

LoopStructural/interpolators/_discrete_interpolator.py

@@ -285,39 +285,66 @@ def add_value_inequality_constraints(self, w: float = 1.0):
             cols = self.support.elements[element[inside]]
             self.add_inequality_constraints_to_matrix(a, points[:, 3:5], cols, 'inequality_value')
 
-    def add_inequality_pairs_constraints(self, w: float = 1.0):
-        pairs = self.get_inequality_pairs_constraints()
-        if pairs['upper'].shape[0] == 0 or pairs['lower'].shape[0] == 0:
-            return
-        upper_interpolation = self.support.get_element_for_location(pairs['upper'])
-        lower_interpolation = self.support.get_element_for_location(pairs['lower'])
-        ij = np.array(
-            [
-                *np.meshgrid(
-                    np.arange(0, int(upper_interpolation[3].sum()), dtype=int),
-                    np.arange(0, int(lower_interpolation[3].sum()), dtype=int),
+    def add_inequality_pairs_constraints(
+        self, w: float = 1.0, upper_bound=np.finfo(float).eps, lower_bound=-np.inf
+    ):
+
+        points = self.get_inequality_pairs_constraints()
+        if points.shape[0] > 0:
+
+            # assemble a list of pairs in the model
+            # this will make pairs even across stratigraphic boundaries
+            # TODO add option to only add stratigraphic pairs
+            pairs = {}
+            k = 0
+            for i in np.unique(points[:, self.support.dimension]):
+                for j in np.unique(points[:, self.support.dimension]):
+                    if i == j:
+                        continue
+                    if tuple(sorted([i, j])) not in pairs:
+                        pairs[tuple(sorted([i, j]))] = k
+                        k += 1
+            pairs = list(pairs.keys())
+            for pair in pairs:
+                upper_points = points[points[:, self.support.dimension] == pair[0]]
+                lower_points = points[points[:, self.support.dimension] == pair[1]]
+
+                upper_interpolation = self.support.get_element_for_location(upper_points)
+                lower_interpolation = self.support.get_element_for_location(lower_points)
+                ij = np.array(
+                    [
+                        *np.meshgrid(
+                            np.arange(0, int(upper_interpolation[3].sum()), dtype=int),
+                            np.arange(0, int(lower_interpolation[3].sum()), dtype=int),
+                        )
+                    ],
+                    dtype=int,
                 )
-            ],
-            dtype=int,
-        )
 
-        ij = ij.reshape(2, -1).T
-        rows = np.arange(0, ij.shape[0], dtype=int)
-        rows = np.tile(rows, (upper_interpolation[1].shape[-1], 1)).T
-        rows = np.hstack([rows, rows])
-        a = upper_interpolation[1][upper_interpolation[3]][ij[:, 0]]  # np.ones(ij.shape[0])
-        a = np.hstack([a, -lower_interpolation[1][lower_interpolation[3]][ij[:, 1]]])
-        cols = np.hstack(
-            [
-                self.support.elements[upper_interpolation[2][upper_interpolation[3]][ij[:, 0]]],
-                self.support.elements[lower_interpolation[2][lower_interpolation[3]][ij[:, 1]]],
-            ]
-        )
+                ij = ij.reshape(2, -1).T
+                rows = np.arange(0, ij.shape[0], dtype=int)
+                rows = np.tile(rows, (upper_interpolation[1].shape[-1], 1)).T
+                rows = np.hstack([rows, rows])
+                a = upper_interpolation[1][upper_interpolation[3]][ij[:, 0]]  # np.ones(ij.shape[0])
+                a = np.hstack([a, -lower_interpolation[1][lower_interpolation[3]][ij[:, 1]]])
+                cols = np.hstack(
+                    [
+                        self.support.elements[
+                            upper_interpolation[2][upper_interpolation[3]][ij[:, 0]]
+                        ],
+                        self.support.elements[
+                            lower_interpolation[2][lower_interpolation[3]][ij[:, 1]]
+                        ],
+                    ]
+                )
 
-        bounds = np.zeros((ij.shape[0], 2))
-        bounds[:, 0] = np.finfo('float').eps
-        bounds[:, 1] = 1e10
-        self.add_inequality_constraints_to_matrix(a, bounds, cols, 'inequality_pairs')
+                bounds = np.zeros((ij.shape[0], 2))
+                bounds[:, 0] = lower_bound
+                bounds[:, 1] = upper_bound
+
+                self.add_inequality_constraints_to_matrix(
+                    a, bounds, cols, f'inequality_pairs_{pair[0]}_{pair[1]}'
+                )
 
     def add_inequality_feature(
         self,
@@ -486,7 +513,7 @@ def build_inequality_matrix(self):
         if len(mats) == 0:
             return None, None
         Q = sparse.vstack(mats)
-        bounds = np.hstack(bounds)
+        bounds = np.vstack(bounds)
         return Q, bounds
 
     def solve_system(
@@ -562,6 +589,12 @@ def solve_system(
             return True
         elif solver == 'admm':
             logger.info("Solving using admm")
+
+            if 'x0' in solver_kwargs:
+                x0 = solver_kwargs['x0'](self.support)
+            else:
+                x0 = np.zeros(A.shape[1])
+            solver_kwargs.pop('x0', None)
             if Q is None:
                 logger.warning("No inequality constraints, using lsmr")
                 return self.solve_system('lsmr', solver_kwargs)
@@ -579,7 +612,7 @@ def solve_system(
                     b,
                     Q,
                     bounds,
-                    x0=solver_kwargs.pop('x0', np.zeros(A.shape[1])),
+                    x0=x0,
                     admm_weight=solver_kwargs.pop('admm_weight', 0.01),
                     nmajor=solver_kwargs.pop('nmajor', 200),
                     linsys_solver_kwargs=solver_kwargs,

LoopStructural/interpolators/_geological_interpolator.py

@@ -186,13 +186,11 @@ def set_value_inequality_constraints(self, points: np.ndarray):
         self.data["inequality"] = points
         self.up_to_date = False
 
-    def set_inequality_pairs_constraints(self, lower_points: np.ndarray, upper_points: np.ndarray):
-        if lower_points.shape[1] < 3:
-            raise ValueError("Inequality pairs constraints must at least have X,Y,Z")
-        if upper_points.shape[1] < 3:
-            raise ValueError("Inequality pairs constraints must at least have X,Y,Z")
-        self.data["inequality_pairs_upper"] = upper_points
-        self.data['inequality_pairs_lower'] = lower_points
+    def set_inequality_pairs_constraints(self, points: np.ndarray):
+        if points.shape[1] < 4:
+            raise ValueError("Inequality pairs constraints must at least have X,Y,Z,rock_id")
+
+        self.data["inequality_pairs"] = points
         self.up_to_date = False
 
     def get_value_constraints(self):
@@ -256,10 +254,7 @@ def get_inequality_value_constraints(self):
         return self.data["inequality"]
 
     def get_inequality_pairs_constraints(self):
-        return {
-            'lower': self.data["inequality_pairs_lower"],
-            'upper': self.data["inequality_pairs_upper"],
-        }
+        return self.data["inequality_pairs"]
 
     # @abstractmethod
     def setup(self, **kwargs):
@@ -341,8 +336,7 @@ def clean(self):
             "tangent": np.zeros((0, 7)),
             "interface": np.zeros((0, 5)),
             "inequality": np.zeros((0, 6)),
-            "inequality_pairs_lower": np.zeros((0, 3)),
-            "inequality_pairs_upper": np.zeros((0, 3)),
+            "inequality_pairs": np.zeros((0, 4)),
         }
         self.up_to_date = False
         self.n_g = 0

LoopStructural/modelling/features/builders/_geological_feature_builder.py

@@ -18,6 +18,7 @@
     tangent_vec_names,
     interface_name,
     inequality_name,
+    pairs_name,
 )
 from ....modelling.features import GeologicalFeature
 from ....modelling.features.builders import BaseBuilder
@@ -283,8 +284,11 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, *
         mask = np.all(~np.isnan(data.loc[:, inequality_name()].to_numpy(float)), axis=1)
         if mask.sum() > 0:
             inequality_data = data.loc[mask, xyz_names() + inequality_name()].to_numpy(float)
-            self.interpolator.set_inequality_constraints(inequality_data)
-
+            self.interpolator.set_value_inequality_constraints(inequality_data)
+        mask = np.all(~np.isnan(data.loc[:, pairs_name()].to_numpy(float)), axis=1)
+        if mask.sum() > 0:
+            pairs_data = data.loc[mask, xyz_names() + pairs_name()].to_numpy(float)
+            self.interpolator.set_inequality_pairs_constraints(pairs_data)
         self.data_added = True
         self._up_to_date = False
 

LoopStructural/utils/helper.py

@@ -267,6 +267,10 @@ def polarity_name():
     return ["polarity"]
 
 
+def pairs_name():
+    return ["pair_id"]
+
+
 def all_heading():
     return (
         xyz_names()
@@ -280,6 +284,7 @@ def all_heading():
         + interface_name()
         + polarity_name()
         + inequality_name()
+        + pairs_name()
     )