WIP: Functional constraints

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/basis/__init__.py

@@ -1,6 +1,6 @@
 from ._identity import Identity
 from ._random_projection import RandomProjection
 from ._svd import SVD
+from ._custom import Custom
 
-
-__all__ = ["Identity", "SVD", "RandomProjection"]
+__all__ = ["Identity", "SVD", "RandomProjection","Custom"]

pysensors/basis/_custom.py

@@ -0,0 +1,88 @@
+"""
+custom mode basis class.
+"""
+from ._base import InvertibleBasis
+from ._base import MatrixMixin
+
+class Custom(InvertibleBasis, MatrixMixin):
+    """
+    Use a custom transformation to map input features to
+    custom modes.
+
+    Assumes the data has already been centered (to have mean 0).
+
+    Parameters
+    ----------
+    n_basis_modes : int, optional (default 10)
+        Number of basis modes to retain. Cannot be larger than
+        the number of features ``n_features``, or the number of examples
+        ``n_examples``.
+    U: The custom basis matrix
+
+    Attributes
+    ----------
+    basis_matrix_ : numpy ndarray, shape (n_features, n_basis_modes)
+        The top n_basis_modes left singular vectors of the training data.
+
+    """
+
+    def __init__(self, U, n_basis_modes=10, **kwargs):
+        if isinstance(n_basis_modes, int) and n_basis_modes > 0:
+            super(Custom, self).__init__()#n_components=n_basis_modes, **kwargs
+            self._n_basis_modes = n_basis_modes
+            self.custom_basis_ = U
+        else:
+            raise ValueError("n_basis_modes must be a positive integer.")
+
+    def fit(self, X):
+        """
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, n_features)
+            The training data.
+
+        Returns
+        -------
+        self : instance
+        """
+        # self.basis_matrix_ = self.custom_basis_[:,: self.n_basis_modes] @ self.custom_basis_[:,: self.n_basis_modes].T @ X[: self.n_basis_modes, :].T.copy()
+        # self.basis_matrix_ = self.custom_basis_ @ self.custom_basis_.T @ X[: self.n_basis_modes, :].T.copy()
+        self.basis_matrix_ = self.custom_basis_[:,:self.n_basis_modes]
+        # self.basis_matrix_ = (X @ self.custom_basis_[:,:self.n_basis_modes] @ self.custom_basis_[:,:self.n_basis_modes].T)[:self.n_basis_modes,:].T
+
+        # self.basis_matrix_ = ((X @ self.custom_basis_).T)[:,:self.n_basis_modes]
+        # self.basis_matrix_ = ((X @ self.custom_basis_ @ self.custom_basis_.T).T)[:,:self.n_basis_modes]
+        return self
+
+    def matrix_inverse(self, n_basis_modes=None):
+        """
+        Get the inverse matrix mapping from measurement space to
+        coordinates with respect to the basis.
+
+        Note that this is not the inverse of the matrix returned by
+        ``self.matrix_representation``. It is the (pseudo) inverse of
+        the matrix whose columns are the basis modes.
+
+        Parameters
+        ----------
+        n_basis_modes : positive int, optional (default None)
+            Number of basis modes to be used to compute inverse.
+
+        Returns
+        -------
+        B : numpy ndarray, shape (n_basis_modes, n_features)
+            The inverse matrix.
+        """
+        n_basis_modes = self._validate_input(n_basis_modes)
+
+        return self.basis_matrix_[:, :n_basis_modes].T
+
+    @property
+    def n_basis_modes(self):
+        """Number of basis modes."""
+        return self._n_basis_modes
+
+    @n_basis_modes.setter
+    def n_basis_modes(self, n_basis_modes):
+        self._n_basis_modes = n_basis_modes
+        self.n_components = n_basis_modes

pysensors/classification/_sspoc.py

@@ -259,7 +259,7 @@ def fit(
 
         if self.threshold is None:
             # Chosen as in Brunton et al. (2016)
-            threshold = np.sqrt(np.sum(s**2)) / (
+            threshold = np.sqrt(np.sum(s ** 2)) / (
                 2 * self.basis_matrix_inverse_.shape[0] * n_classes
             )
         else:

pysensors/optimizers/__init__.py

@@ -1,8 +1,9 @@
 from ._ccqr import CCQR
 from ._qr import QR
-
+from ._gqr import GQR
 
 __all__ = [
     "CCQR",
     "QR",
+    "GQR"
 ]

pysensors/optimizers/_gqr.py

@@ -0,0 +1,125 @@
+import numpy as np
+import pysensors
+
+from pysensors.optimizers._qr import QR
+
+import matplotlib.pyplot as plt
+from sklearn import datasets
+from sklearn import metrics
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+import pysensors as ps
+from matplotlib.patches import Circle
+from pysensors.utils._norm_calc import returnInstance as normCalcReturnInstance
+
+class GQR(QR):
+    """
+    General QR optimizer for sensor selection.
+    Ranks sensors in descending order of "importance" based on
+    reconstruction accuracy. This is an extension that requires a more intrusive
+    access to the QR optimizer to facilitate a more adaptive optimization. This 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. and if no constrained region it should converge to the results from QR optimizer.
+
+    See the following reference for more information
+        Manohar, Krithika, et al.
+        "Data-driven sparse sensor placement for reconstruction:
+        Demonstrating the benefits of exploiting known patterns."
+        IEEE Control Systems Magazine 38.3 (2018): 63-86.
+
+    @ authors: Niharika Karnik (@nkarnik2999), Mohammad Abdo (@Jimmy-INL), and Krithika Manohar (@kmanohar)
+    """
+    def __init__(self):
+        """
+        Attributes
+        ----------
+        pivots_ : np.ndarray, shape [n_features]
+            Ranked list of sensor locations.
+        idx_constrained : np.ndarray, shape [No. of constrained locations]
+            Column Indices of the sensors in the constrained locations.
+        n_sensors : integer,
+            Total number of sensors
+        n_const_sensors : integer,
+            Total number of sensors required by the user in the constrained region.
+        all_sensors : np.ndarray, shape [n_features]
+            Optimally placed list of sensors obtained from QR pivoting algorithm.
+        constraint_option : string,
+            max_n_const_sensors : The number of sensors in the constrained region should be less than or equal to n_const_sensors.
+            exact_n_const_sensors : The number of sensors in the constrained region should be exactly equal to n_const_sensors.
+        """
+        self.pivots_ = None
+        self.idx_constrained = []
+        self.n_sensors = None
+        self.n_const_sensors = 0
+        self.all_sensors = []
+        self.constraint_option = ''
+        self.nx = None
+        self.ny = None
+        self.r = 1
+
+    def fit(self,basis_matrix,**optimizer_kws):
+        """
+        Parameters
+        ----------
+        basis_matrix: np.ndarray, shape [n_features, n_samples]
+            Matrix whose columns are the basis vectors in which to
+            represent the measurement data.
+        optimizer_kws: dictionary, optional
+            Keyword arguments to be passed to the qr method.
+
+        Returns
+        -------
+        self: a fitted :class:`pysensors.optimizers.GQR` instance
+        """
+        [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
+
+        ## Assertions and checks:
+        # if self.n_sensors > n_features - max_const_sensors + self.nConstrainedSensors:
+        #     raise IOError ("n_sensors cannot be larger than n_features - all possible locations in the constrained area + allowed constrained sensors")
+        # if self.n_sensors > n_samples + self.nConstrainedSensors: ## Handling zero constraint?
+        #     raise IOError ("Currently n_sensors should be less than min(number of samples, number of modes) + number of constrained sensors,\
+        #                    got: n_sensors = {}, n_samples + const_sensors = {} + {} = {}".format(self.n_sensors,n_samples,self.nConstrainedSensors,n_samples+self.nConstrainedSensors))
+
+        # Initialize helper variables
+        R = basis_matrix.conj().T.copy()
+        p = np.arange(n_features)
+        k = min(n_samples, n_features)
+
+
+        for j in range(k):
+            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_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:
+                u = r[:, i_piv] / dlen
+                u[0] += np.sign(u[0]) + (u[0] == 0)
+                u /= np.sqrt(abs(u[0]))
+            else:
+                u = r[:, i_piv]
+                u[0] = np.sqrt(2)
+
+            # Track column pivots
+            i_piv += j # true permutation index is i_piv shifted by the iteration counter j
+            p[[j, i_piv]] = p[[i_piv, j]]
+
+            # Switch columns
+            R[:, [j, i_piv]] = R[:, [i_piv, j]]
+
+            # Apply reflector
+            R[j:, j:] -= np.outer(u, np.dot(u, R[j:, j:]))
+            R[j + 1 :, j] = 0
+        self.pivots_ = p
+        return self

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,10 +1,50 @@
 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_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
+from ._validation import determinant
+from ._validation import relative_reconstruction_error
 
 __all__ = [
     "constrained_binary_solve",
     "constrained_multiclass_solve",
     "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",
+    "determinant",
+    "relative_reconstruction_error"
 ]