Adding a constraints class for the GQR optimizer

examples/twistParabolicConstraint.py

@@ -0,0 +1,44 @@
+
+import numpy as np
+
+def twistParabolicConstraint(x,y,**kwargs):
+    """
+    Function for evaluating constrained sensor locations on the grid by returning a negative value (if index is constrained) and a positive value (if index is unconstrained).
+
+    Parameters
+    ----------
+    x: float, x coordinate of all grid-points considered for sensor placement
+    y : float, y coordinate of all grid-points considered for sensor placement
+
+    **kwargs : h : float, x-coordinate of the vertex of the parabola we want to be constrained;
+               k : float, y-coordinate of the vertex of the parabola we want to be constrained;
+               a : float, The x-coordinate of the focus of the parabola.
+
+    Returns
+    -------
+    g : np.darray, shape [No. of grid points], 
+        A boolean array for every single grid point based on whether the grid point lies in the constrained region or not
+    """
+    # make sure the length of x is the same as y
+    assert len(x) == len(y)
+    if ('h' not in kwargs.keys()) or (kwargs['h'] == None) :
+        kwargs['h'] = 0.025
+    if ('k' not in kwargs.keys()) or (kwargs['k'] == None) :
+        kwargs['k'] = 0
+    if ('a' not in kwargs.keys()) or (kwargs['a'] == None) :
+        kwargs['a'] = 100
+    # initialize the constraint evaluation function g
+    g1 = np.zeros(len(x),dtype=float) - 1
+    g2 = np.zeros(len(x),dtype=float) - 1
+    g = np.zeros(len(x),dtype=float)
+    # loop over all given location and check if they violate the constraints
+    # make sure the retuned value is negative if it is in the constrained area and positive otherwise
+    for i in range(len(x)):
+        # circle of center (h,k)=(0,0) and radius 2.5 cm
+        g1[i] = (kwargs['a']*(x[i]-kwargs['h'])**2) - (y[i]-kwargs['k'])
+        #
+        # Second constraint:
+        g2[i] = y[i] - 0.2
+        if bool(g1[i]>=0) == bool(g2[i]>=0):
+            g[i] = (bool(g1[i]>=0) and bool(g2[i]>=0))-1
+    return g

examples/userExplicitConstraint1.py

@@ -0,0 +1,10 @@
+import numpy as np
+
+def userExplicitConstraint1(x,y,**kwargs):
+  '''
+  '''
+  assert len(x) == len(y)
+  g = np.zeros(len(x),dtype=float) 
+  for i in range(len(x)):
+    g[i] = ((x[i]-30)**2 + (y[i]-40)**2) - 5**2
+  return g

examples/userExplicitConstraint2.py

@@ -0,0 +1,4 @@
+
+def userExplicitConstraint2(x,y,**kwargs):
+  g = (x-20)**2 + (y-10)**2 - 49
+  return g

pysensors/optimizers/_gqr.py

@@ -74,7 +74,7 @@ def fit(self,basis_matrix,**optimizer_kws):
         [setattr(self,name,optimizer_kws.get(name,getattr(self,name))) for name in optimizer_kws.keys()]
         self._norm_calc_Instance = normCalcReturnInstance(self, self.constraint_option)
         n_features, n_samples = basis_matrix.shape  # We transpose basis_matrix below
-        max_const_sensors = len(self.idx_constrained) # Maximum number of sensors allowed in the constrained region
+        # max_const_sensors = len(self.idx_constrained) # Maximum number of sensors allowed in the constrained region
 
         ## Assertions and checks:
         # if self.n_sensors > n_features - max_const_sensors + self.nConstrainedSensors:
@@ -93,9 +93,14 @@ def fit(self,basis_matrix,**optimizer_kws):
             r = R[j:, j:]
 
             # Norm of each column
+            if j == 0:
+                dlens_old = np.sqrt(np.sum(np.abs(r) ** 2, axis=0))
+            else:
+                dlens_old = dlens
             dlens = np.sqrt(np.sum(np.abs(r) ** 2, axis=0))
-            dlens_updated = self._norm_calc_Instance(self.idx_constrained, dlens, p, j, self.n_const_sensors, dlens_old=dlens, all_sensors=self.all_sensors, n_sensors=self.n_sensors, nx=self.nx, ny=self.ny, r=self.r)
-            i_piv = np.argmax(dlens_updated)
+            dlens_updated = self._norm_calc_Instance(self.idx_constrained, dlens, p, j, self.n_const_sensors, dlens_old=dlens_old, all_sensors=self.all_sensors, n_sensors=self.n_sensors, nx=self.nx, ny=self.ny, r=self.r)
+            # i_piv = np.argmax(dlens_updated)
+            i_piv = np.where(dlens_updated==dlens_updated.max())[0][0]
             dlen = dlens_updated[i_piv]
 
             if dlen > 0:

pysensors/reconstruction/_sspor.py

@@ -9,6 +9,7 @@
 from ..basis import Identity
 from ..optimizers import CCQR
 from ..optimizers import QR
+from ..optimizers import GQR
 from ..utils import validate_input
 
 
@@ -510,7 +511,7 @@ def _validate_n_sensors(self):
         # If n_sensors exceeds n_samples, the cost-constrained QR algorithm may
         # place sensors in constrained areas.
         if (
-            isinstance(self.optimizer, CCQR)
+            (isinstance(self.optimizer, CCQR) or isinstance(self.optimizer, QR) or isinstance(self.optimizer, GQR))
             and self.n_sensors > self.basis_matrix_.shape[1]
         ):
             warnings.warn(

pysensors/utils/__init__.py

@@ -1,8 +1,22 @@
 from ._base import validate_input
 from ._optimizers import constrained_binary_solve
 from ._optimizers import constrained_multiclass_solve
-from ._constraints import get_constraind_sensors_indices
-from ._constraints import get_constrained_sensors_indices_linear
+from ._constraints import get_constrained_sensors_indices
+from ._constraints import get_constrained_sensors_indices_dataframe
+from ._constraints import BaseConstraint
+from ._constraints import Circle
+from ._constraints import Cylinder
+from ._constraints import Line
+from ._constraints import Ellipse
+from ._constraints import Parabola
+from ._constraints import Polygon
+from ._constraints import UserDefinedConstraints
+from ._constraints import check_constraints
+from ._constraints import constraints_eval
+
+from ._constraints import load_functional_constraints
+from ._constraints import get_coordinates_from_indices
+from ._constraints import get_indices_from_coordinates
 from ._norm_calc import exact_n
 from ._norm_calc import max_n
 from ._norm_calc import predetermined
@@ -15,8 +29,19 @@
     "validate_input",
     "get_constraind_sensors_indices",
     "get_constrained_sensors_indices_linear",
+    "BaseConstraint",
+    "Circle",
+    "Cylinder"
+    "Line",
+    "Parabola",
+    "Polygon",
+    "Ellipse",
+    "UserDefinedConstraints"
     "box_constraints",
+    "constraints_eval",
     "functional_constraints",
+    "get_coordinates_from_indices",
+    "get_indices_from_coordinates",
     "exact_n",
     "max_n",
     "predetermined",