test: ✅ new test for equality constraints

Added a new test for equality constraints using FDI. Need to add more options for other solver.

Author: Lachlan Grose

LoopStructural/interpolators/discrete_interpolator.py

@@ -277,6 +277,7 @@ def add_equality_constraints(self, node_idx, values, name="undefined"):
         gi[self.region] = np.arange(0, self.nx)
         idc = gi[node_idx]
         outside = ~(idc == -1)
+        
         self.equal_constraints[name] = {
             "A": np.ones(idc[outside].shape[0]),
             "B": values[outside],
@@ -434,15 +435,22 @@ def build_matrix(self, square=True, damp=0.0, ie=False):
             # c are the node values and y are the
             # lagrange multipliers#
             nc = 0
+            a = []
+            rows = []
+            cols = []
+            b = []
             for c in self.equal_constraints.values():
                 b.extend((c["B"]).tolist())
+                aa = c["A"].flatten()
                 mask = aa == 0
-                a.extend(c["A"].flatten()[~mask].tolist())
+                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),
+                shape=(self.eq_const_c, self.nx),
                 dtype=float,
             ).tocsr()
 

tests/unit_tests/interpolator/test_solvers.py

@@ -93,10 +93,7 @@ def test_inequality_FDI_nodes():
     data['l'] = -3
     data['u'] = 10
     randind = np.arange(0,len(data))
-    # np.random.shuffle(randind)
-    # # data.loc[randind[:int(50*50*.5)],'val'] = np.nan
-    # np.random.shuffle(randind)
-    # data.loc[randind[:int(50*50*.1)],'val'] = 0
+
     origin = np.array([-0.1,-0.1,-0.1])
     maximum = np.array([1.1,1.1,1.1])
     nsteps = np.array([20,20,20])
@@ -106,18 +103,30 @@ def test_inequality_FDI_nodes():
     interpolator.set_value_constraints(data[['X','Y','Z','val','w']].to_numpy())
     interpolator.set_inequality_constraints(data[['X','Y','Z','l','u']].to_numpy())
     interpolator._setup_interpolator()
-    # col = np.arange(0,interpolator.nx,dtype=int)
-    # col = np.tile(col, (interpolator.nx, 1)).T
-    # interpolator.add_inequality_constraints_to_matrix(np.eye(interpolator.nx),
-    #                                                 np.zeros(interpolator.nx)-4,
-    #                                                 np.zeros(interpolator.nx)+np.inf,
-    #                                                 col
-    #                                                 )
     interpolator.solve_system(solver='osqp')
 
-    # print(np.sum(interpolator.evaluate_value(data[['X','Y','Z']].to_numpy())-data[['val']].to_numpy())/len(data))
-    # assert np.sum(interpolator.evaluate_value(data[['X','Y','Z']].to_numpy())-data[['val']].to_numpy())/len(data) < 0.5
+def test_equality_FDI_nodes():
+    xy = np.array(np.meshgrid(np.linspace(0,1,50),np.linspace(0,1,50))).T.reshape(-1,2)
+    xyz = np.hstack([xy,np.zeros((xy.shape[0],1))])
+    data = pd.DataFrame(xyz,columns=['X','Y','Z'])
+    data['val'] = np.sin(data['X'])
+    data['w'] = 1
+    data['feature_name'] = 'strati'
+    origin = np.array([-0.1,-0.1,-0.1])
+    maximum = np.array([1.1,1.1,1.1])
+    nsteps = np.array([20,20,20])
+    step_vector = (maximum-origin)/nsteps
+    grid = StructuredGrid(origin=origin,nsteps=nsteps,step_vector=step_vector)
+    interpolator = FDI(grid)
+    interpolator.set_value_constraints(data[['X','Y','Z','val','w']].to_numpy())
+
+    node_idx = np.arange(0,interpolator.nx)[interpolator.support.nodes[:,2]>.9]
+    interpolator.add_equality_constraints(node_idx,np.ones(node_idx.shape[0]),name='top')
+    interpolator._setup_interpolator()
+    interpolator.solve_system(solver='cg')
+
 if __name__ == '__main__':
-    test_inequality_FDI()
-    test_inequality_FDI_nodes()
-    test_FDI()
+    # test_inequality_FDI()
+    # test_inequality_FDI_nodes()
+    # test_FDI()
+    test_equality_FDI_nodes()