Mods to the constraint class

README.rst

@@ -55,6 +55,59 @@ Use the ``predict`` method to reconstruct a new function sampled at the chosen s
   :alt: A plot showing the function to be reconstructed, the learned sensor locations, and the reconstruction.
   :figclass: align-center
 
+Reconstruction with constraints
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+In most engineering applications, certain areas within the region of interest might allow a limited number of sensors or none at all.
+We develop a data-driven technique that incorporates constraints into an optimization framework for sensor placement, with the primary objective 
+of minimizing reconstruction errors under noisy sensor measurements. 
+
+This work has been implemented in the general QR optimizer for sensor selection. 
+This is an extension that requires a more intrusive access to the QR optimizer to facilitate a more adaptive optimization. It is a generalized version of cost constraints
+in the sense that users can allow `n_const_sensors` in the constrained area. If n = 0 this converges to the CCQR results. If there is 
+no constrained region it should converge to the results from QR optimizer.
+
+To implement constrained sensing we initialize the optimizer GQR and provide it additional kwargs such as the constrained region, number of allowable 
+sensors in the constrained region and the type of constraint. 
+
+Three strategies to deal with constraints are currently developed: 
+
+* ``max_n`` - Number of sensors in the constrained region should be less than or equal to the allowable constrained sensors.
+
+* ``exact_n`` - Number of sensors in the constrained region should be exactly equal to the allowable constrained sensors.
+
+* ``predetermined`` - A number of sensor locations are predetermined and the aim is to optimize the rest.
+
+.. code-block:: python
+
+  optimizer_exact = ps.optimizers.GQR()
+  opt_exact_kws={'idx_constrained':sensors_constrained,
+          'n_sensors':n_sensors,
+          'n_const_sensors':n_const_sensors,
+          'all_sensors':all_sensors,
+          'constraint_option':"exact_n"}
+
+We have further provided functions to compute the sensors in the constrained regions. For example if the user provides the center and radius of a circular
+constrained region, the constraints in utils compute the constrained sensor indices. Direct constraint plotting capabilities have also been developed. 
+
+The constrained shapes currently implemented are: 
+
+* ``Circle``
+
+* ``Cylinder``
+
+* ``Line``
+
+* ``Parabola``
+
+* ``Ellipse``
+
+* ``Polygon``
+
+* ``UserDefinedConstraints``
+
+  - This type of constraint has the ability to take in either a function from the user or a 
+  .py file which contains a functional definition of the constrained region.
+
 Classification
 ^^^^^^^^^^^^^^
 Classification is the problem of predicting which category an example belongs to, given a set of training data (e.g. determining whether digital photos are of dogs or cats).

examples/README.rst

@@ -37,6 +37,17 @@ ocean at any given point.
 Reproduces an example from `Manohar et al. (2018) <https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=8361090>`_
 where sensor locations are learned for a monomial basis for the task of reconstruction.
 
+`Spatial constraints example <https://python-sensors.readthedocs.io/en/latest/examples/spatially_constrained_qr.html>`_
+----------------------------------------------------------------------------------------------------
+Sensor locations are learned for a constrained sensing problem for the task of reconstruction. For further details of constrained sensing refer 
+`Karnik et al. (2024) <https://ieeexplore.ieee.org/abstract/document/10453459>`_
+
+`Functional constraints example <https://python-sensors.readthedocs.io/en/latest/examples/functional_constraints.html>`_
+----------------------------------------------------------------------------------------------------
+Sensor locations are learned for various shapes of constrained regions. For further details refer to 
+`Karnik et al. (2024) <https://www.mdpi.com/1996-1073/17/13/3355>`_
+
+
 
 Full table of contents
 ----------------------

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",