Merge branch 'inequality'

Author: Lachlan Grose

.github/workflows/windows_CI.yml

@@ -18,7 +18,7 @@ jobs:
       - name: Installing dependencies
         shell: bash -l {0}
         run: |
-          conda install -c conda-forge cython numpy scipy scikit-image scikit-learn pyamg flake8 pytest networkx
+          conda install -c conda-forge cython numpy scipy scikit-image scikit-learn pyamg flake8 pytest networkx osqp
       - name: Checking formatting of code
         shell: bash -l {0}
         run: |

LoopStructural/interpolators/__init__.py

@@ -4,20 +4,24 @@
 """
 from enum import IntEnum
 
+
 class InterpolatorType(IntEnum):
-    '''
+    """
     Enum for the different interpolator types
-    
+
     1-9 should cover interpolators with supports
     9+ are data supported
-    '''
+    """
+
     BASE = 0
     BASE_DISCRETE = 1
     FINITE_DIFFERENCE = 2
     DISCRETE_FOLD = 3
     PIECEWISE_LINEAR = 4
     BASE_DATA_SUPPORTED = 10
     SURFE = 11
+
+
 from LoopStructural.interpolators.geological_interpolator import GeologicalInterpolator
 from LoopStructural.interpolators.discrete_interpolator import DiscreteInterpolator
 from LoopStructural.interpolators.structured_tetra import TetMesh
@@ -32,6 +36,3 @@ class InterpolatorType(IntEnum):
 from LoopStructural.interpolators.discrete_fold_interpolator import (
     DiscreteFoldInterpolator,
 )
-
-
-

LoopStructural/interpolators/base_structured_3d_support.py

@@ -27,7 +27,7 @@ def __init__(
         # we use property decorators to update these when different parts of
         # the geometry need to change
         # inisialise the private attributes
-        if np.any(step_vector==0):
+        if np.any(step_vector == 0):
             logger.warning(f"Step vector {step_vector} has zero values")
         self._nsteps = np.array(nsteps, dtype=int)
         self._step_vector = np.array(step_vector)
@@ -74,7 +74,7 @@ def step_vector(self, step_vector):
         newsteps = self._nsteps / change_factor
         self._nsteps = np.ceil(newsteps).astype(int)
         self._step_vector = step_vector
-        
+
     @property
     def origin(self):
         return self._origin
@@ -143,7 +143,7 @@ def nodes(self):
         max = self.origin + self.nsteps_cells * self.step_vector
         if np.any(np.isnan(self.nsteps)):
             raise ValueError("Cannot resize mesh nsteps is NaN")
-        if  np.any(np.isnan(self.origin)):
+        if np.any(np.isnan(self.origin)):
             raise ValueError("Cannot resize mesh origin is NaN")
 
         x = np.linspace(self.origin[0], max[0], self.nsteps[0])

LoopStructural/interpolators/discrete_fold_interpolator.py

@@ -6,9 +6,7 @@
 import numpy as np
 from LoopStructural.interpolators.cython.dsi_helper import fold_cg
 
-from LoopStructural.interpolators import (
-    PiecewiseLinearInterpolator, InterpolatorType
-)
+from LoopStructural.interpolators import PiecewiseLinearInterpolator, InterpolatorType
 
 from LoopStructural.utils import getLogger
 

LoopStructural/interpolators/discrete_interpolator.py

@@ -138,7 +138,6 @@ def reset(self):
 
         """
         logger.debug("Resetting interpolation constraints")
-        
 
     def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"):
         """
@@ -255,7 +254,7 @@ def remove_constraints_from_least_squares(
 
         pass
 
-    def add_equality_constraints(self, node_idx, values,name="undefined"):
+    def add_equality_constraints(self, node_idx, values, name="undefined"):
         """
         Adds hard constraints to the least squares system. For now this just
         sets
@@ -283,8 +282,7 @@ def add_equality_constraints(self, node_idx, values,name="undefined"):
             "B": values[outside].tolist(),
             "col": idc[outside].tolist(),
             # "w": w,
-            "row": np.arange(self.eq_const_c, self.eq_const_c+idc[outside].shape[0])
-            
+            "row": np.arange(self.eq_const_c, self.eq_const_c + idc[outside].shape[0]),
         }
         self.eq_const_c += idc[outside].shape[0]
         # ,'C':np.ones(idc[outside].shape[0]).tolist(),}
@@ -317,20 +315,38 @@ def add_inequality_constraints_to_matrix(self, A, l, u, idc, name="undefined"):
 
         """
         # map from mesh node index to region node index
-        gi = np.zeros(self.support.n_nodes)
+        gi = np.zeros(self.support.n_nodes,dtype=int)
         gi[:] = -1
-        gi[self.region] = np.arange(0, self.nx)
+        gi[self.region] = np.arange(0, self.nx,dtype=int)
         idc = gi[idc]
-        rows = np.arange(self.ineq_const_c, self.ineq_const_c+idc.shape[0])
-        rows =  np.tile(rows, (A.shape[-1], 1)).T
-        self.ineq_constraints[name] = {
-            "A": A,
-            "l": l,
-            "col": idc,
-            "u": u,
-            "row": rows
-        }
-        self.ineq_const_c+=idc.shape[0]
+        rows = np.arange(self.ineq_const_c, self.ineq_const_c + idc.shape[0])
+        rows = np.tile(rows, (A.shape[-1], 1)).T
+        self.ineq_constraints[name] = {"A": A, "l": l, "col": idc, "u": u, "row": rows}
+        self.ineq_const_c += idc.shape[0]
+
+    def add_inequality_feature(self, feature, lower=True, mask=None):
+
+        # add inequality value for the nodes of the mesh
+        # flag lower determines whether the feature is a lower bound or upper bound
+        # mask is just a boolean array determining which nodes to apply it to
+
+        value = feature(self.support.nodes)
+        if mask is None:
+            mask = np.ones(value.shape[0], dtype=bool)
+        l = np.zeros(value.shape[0]) - np.inf
+        u = np.zeros(value.shape[0]) + np.inf
+        mask = np.logical_and(mask, ~np.isnan(value))
+        if lower:
+            l[mask] = value[mask]
+        if lower == False:
+            u[mask] = value[mask]
+
+        self.add_inequality_constraints_to_matrix(
+            np.ones((value.shape[0], 1)),
+            l,
+            u,
+            np.arange(0, self.nx, dtype=int),
+        )
 
     def add_tangent_constraints(self, w=1.0):
         """
@@ -426,20 +442,21 @@ def build_matrix(self, square=True, damp=0.0, ie=False):
             nc = 0
             for c in self.equal_constraints.values():
                 aa = (c["A"]).flatten()
-                b.extend((c["B"] ).tolist())
+                b.extend((c["B"]).tolist())
                 mask = aa == 0
                 a.extend(aa[~mask].tolist())
                 rows.extend(c["row"].flatten()[~mask].tolist())
                 cols.extend(c["col"].flatten()[~mask].tolist())
             C = coo_matrix(
-                
-                    (np.array(a), (np.array(rows), cols)), shape=(self.eq_const_c, self.nx), dtype=float
-                ).tocsr()
-                
+                (np.array(a), (np.array(rows), cols)),
+                shape=(self.eq_const_c_, self.nx),
+                dtype=float,
+            ).tocsr()
+
             d = np.array(b)
             ATA = bmat([[ATA, C.T], [C, None]])
             ATB = np.hstack([ATB, d])
-        
+
         if isinstance(damp, bool):
             if damp == True:
                 damp = np.finfo("float").eps
@@ -449,7 +466,7 @@ def build_matrix(self, square=True, damp=0.0, ie=False):
             logger.info("Adding eps to matrix diagonal")
             ATA += eye(ATA.shape[0]) * damp
         if len(self.ineq_constraints) > 0 and ie:
-            print('using inequality constraints')
+            print("using inequality constraints")
             a = []
             l = []
             u = []
@@ -465,16 +482,19 @@ def build_matrix(self, square=True, damp=0.0, ie=False):
                 rows.extend(c["row"].flatten()[~mask].tolist())
                 cols.extend(c["col"].flatten()[~mask].tolist())
             Aie = coo_matrix(
-            (np.array(a), (np.array(rows), cols)), shape=(self.ineq_const_c, self.nx), dtype=float
-            ).tocsc()    # .tocsr()
+                (np.array(a), (np.array(rows), cols)),
+                shape=(self.ineq_const_c, self.nx),
+                dtype=float,
+            ).tocsc()  # .tocsr()
 
             uie = np.array(u)
             lie = np.array(l)
 
             return ATA, ATB, Aie.T.dot(Aie), Aie.T.dot(uie), Aie.T.dot(lie)
         return ATA, ATB
-    def _solve_osqp(self, P, A, q, l, u):
-        
+
+    def _solve_osqp(self, P, A, q, l, u,mkl=False):
+
         try:
             import osqp
         except ImportError:
@@ -502,13 +522,21 @@ def _solve_osqp(self, P, A, q, l, u):
         # # Solve problem
         # res = prob.solve()
 
-
         # Create an OSQP object
         prob = osqp.OSQP()
 
         # Setup workspace
         # osqp likes csc matrices
-        prob.setup(P.tocsc(), np.array(q), A.tocsc(), np.array(u), np.array(l))
+        linsys_solver='qdldl'
+        if mkl:
+            linsys_solver='mkl pardiso'
+        
+        try:
+            prob.setup(P.tocsc(), np.array(q), A.tocsc(), np.array(u), np.array(l),linsys_solver=linsys_solver)
+        except ValueError:
+            if mkl:
+                logger.error('MKL solver library path not correct. Please add to LD_LIBRARY_PATH') 
+                raise LoopImportError("Cannot import MKL pardiso")
         res = prob.solve()
         return res.x
 
@@ -676,8 +704,8 @@ def _solve(self, solver="cg", **kwargs):
             damp = True
         if solver == "lsqr":
             A, B = self.build_matrix(False)
-        elif solver =='osqp':
-            P, q, A, l, u = self.build_matrix(True,ie=True)
+        elif solver == "osqp":
+            P, q, A, l, u = self.build_matrix(True, ie=True)
         else:
             A, B = self.build_matrix(damp=damp)
 
@@ -702,8 +730,8 @@ def _solve(self, solver="cg", **kwargs):
         if solver == "external":
             logger.warning("Using external solver")
             self.c[self.region] = kwargs["external"](A, B)[: self.nx]
-        if solver == 'osqp':
-            self.c[self.region] = self._solve_osqp(P, A, q, l, u)#, **kwargs)
+        if solver == "osqp":
+            self.c[self.region] = self._solve_osqp(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):

LoopStructural/interpolators/finite_difference_interpolator.py

@@ -132,6 +132,7 @@ def _setup_interpolator(self, **kwargs):
         self.add_vaue_constraints(self.interpolation_weights["cpw"])
         self.add_tangent_constraints(self.interpolation_weights["tpw"])
         self.add_interface_constraints(self.interpolation_weights["ipw"])
+        self.add_inequality_constraints()
 
     def copy(self):
         """
@@ -161,12 +162,11 @@ def add_vaue_constraints(self, w=1.0):
             node_idx, inside = self.support.position_to_cell_corners(points[:, :3])
             # print(points[inside,:].shape)
 
-            gi = np.zeros(self.support.n_nodes)
+            gi = np.zeros(self.support.n_nodes, dtype=int)
             gi[:] = -1
-            gi[self.region] = np.arange(0, self.nx)
+            gi[self.region] = np.arange(0, self.nx, dtype=int)
             idc = np.zeros(node_idx.shape)
             idc[:] = -1
-
             idc[inside, :] = gi[node_idx[inside, :]]
             inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
             a = self.support.position_to_dof_coefs(points[inside, :3])
@@ -179,8 +179,41 @@ def add_vaue_constraints(self, w=1.0):
                 w=w * points[inside, 4],
                 name="value",
             )
-            if np.sum(inside)<=0:
-                logger.warning(f"{self.propertyname}: {np.sum(~inside)} value constraints not added: outside of model bounding box")
+            if np.sum(inside) <= 0:
+                logger.warning(
+                    f"{self.propertyname}: {np.sum(~inside)} value constraints not added: outside of model bounding box"
+                )
+
+    def add_inequality_constraints(self, w=1.0):
+        points = self.get_inequality_constraints()
+        # check that we have added some points
+        if points.shape[0] > 0:
+            node_idx, inside = self.support.position_to_cell_corners(points[:, :3])
+            # print(points[inside,:].shape)
+
+            gi = np.zeros(self.support.n_nodes, dtype=int)
+            gi[:] = -1
+            gi[self.region] = np.arange(0, self.nx, dtype=int)
+            idc = np.zeros(node_idx.shape, dtype=int)
+            idc[:] = -1
+
+            idc[inside, :] = gi[node_idx[inside, :]]
+            inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
+            a = self.support.position_to_dof_coefs(points[inside, :3])
+            # a*=w
+            # a/=np.product(self.support.step_vector)
+            self.add_inequality_constraints_to_matrix(
+                a.T,
+                points[inside, 3],
+                points[inside, 4],
+                idc[inside, :],
+                name="value_inequality",
+            )
+            if np.sum(inside) <= 0:
+                logger.warning(
+                    f"{self.propertyname}: {np.sum(~inside)} value constraints not added: outside of model bounding box"
+                )
+
     def add_interface_constraints(
         self, w=1.0
     ):  # for now weight all value points the same
@@ -245,8 +278,6 @@ def add_interface_constraints(
                     name="interface_{}".format(unique_id),
                 )
 
-       
-
     def add_gradient_constraints(self, w=1.0):
         """
 
@@ -297,8 +328,11 @@ def add_gradient_constraints(self, w=1.0):
             self.add_constraints_to_least_squares(
                 A, B, idc[inside, :], w=w * self.vol, name="gradient"
             )
-            if np.sum(inside)<=0:
-                logger.warning(f"{self.propertyname}: {np.sum(~inside)} norm constraints not added: outside of model bounding box")
+            if np.sum(inside) <= 0:
+                logger.warning(
+                    f"{self.propertyname}: {np.sum(~inside)} norm constraints not added: outside of model bounding box"
+                )
+
     def add_norm_constraints(self, w=1.0):
         """
         Add constraints to control the norm of the gradient of the scalar field
@@ -358,8 +392,10 @@ def add_norm_constraints(self, w=1.0):
                 w=w * self.vol,
                 name="norm",
             )
-            if np.sum(inside)<=0:
-                logger.warning(f"{self.propertyname}: {np.sum(~inside)} norm constraints not added: outside of model bounding box")
+            if np.sum(inside) <= 0:
+                logger.warning(
+                    f"{self.propertyname}: {np.sum(~inside)} norm constraints not added: outside of model bounding box"
+                )
 
     def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
         """
@@ -410,8 +446,10 @@ def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
             self.add_constraints_to_least_squares(
                 A, B, idc[inside, :], w=w * self.vol, name="gradient orthogonal"
             )
-            if np.sum(inside)<=0:
-                logger.warning(f"{self.propertyname}: {np.sum(~inside)} gradient constraints not added: outside of model bounding box")
+            if np.sum(inside) <= 0:
+                logger.warning(
+                    f"{self.propertyname}: {np.sum(~inside)} gradient constraints not added: outside of model bounding box"
+                )
 
     def add_regularisation(self, operator, w=0.1):
         """