Merge pull request #17 from niharika2999/niharika-krithika-mohammad-regionConstraints

Adding region constraints to pysensors PR

Author: kmanohar

examples/spatially_constrained_qr.ipynb

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

@@ -0,0 +1,91 @@
+"""
+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):
+        '''
+            kwargs : Not defined but added to remain consistent with prior basis functions. 
+        '''
+        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/basis/_custom.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/classification/_sspoc.py

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

pysensors/optimizers/__init__.py

@@ -0,0 +1,117 @@
+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.
+        
+        Niharika Karnik, Mohammad G. Abdo, Carlos E. Estrada Perez, Jun Soo Yoo, Joshua J. Cogliati, Richard S. Skifton, 
+        Pattrick Calderoni, Steven L. Brunton, and Krithika Manohar. 
+        Optimal Sensor Placement with Adaptive Constraints for Nuclear Digital Twins. 2023. arXiv: 2306 . 13637 [math.OC].
+
+    @ 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
+
+        # 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
+            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)
+            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/optimizers/_gqr.py

@@ -1,10 +1,25 @@
 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 ._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",
+    "box_constraints",
+    "functional_constraints",
+    "exact_n",
+    "max_n",
+    "predetermined",
+    "determinant",
+    "relative_reconstruction_error"
 ]

pysensors/utils/__init__.py

@@ -0,0 +1,72 @@
+
+"""
+Various utility functions for mapping constrained sensors locations with the column indices for class GQR.
+"""
+
+import numpy as np
+
+
+def get_constraind_sensors_indices(x_min, x_max, y_min, y_max, nx, ny, all_sensors):
+    """
+    Function for mapping constrained sensor locations on the grid with the column indices of the basis_matrix.
+
+    Parameters
+    ----------
+    x_min: int, lower bound for the x-axis constraint
+    x_max : int, upper bound for the x-axis constraint
+    y_min : int, lower bound for the y-axis constraint
+    y_max : int, upper bound for the y-axis constraint
+    nx : int, image pixel (x dimensions of the grid)
+    ny : int, image pixel (y dimensions of the grid)
+    all_sensors : np.ndarray, shape [n_features], ranked list of sensor locations.
+
+    Returns
+    -------
+    idx_constrained : np.darray, shape [No. of constrained locations], array which contains the constrained
+        locations of the grid in terms of column indices of basis_matrix.
+    """
+    n_features = len(all_sensors)
+    image_size = int(np.sqrt(n_features))
+    a = np.unravel_index(all_sensors, (nx,ny))
+    constrained_sensorsx = []
+    constrained_sensorsy = []
+    for i in range(n_features):
+        if (a[0][i] >= x_min and a[0][i] <= x_max) and (a[1][i] >= y_min and a[1][i] <= y_max):
+            constrained_sensorsx.append(a[0][i])
+            constrained_sensorsy.append(a[1][i])
+
+    constrained_sensorsx = np.array(constrained_sensorsx)
+    constrained_sensorsy = np.array(constrained_sensorsy)
+    constrained_sensors_array = np.stack((constrained_sensorsy, constrained_sensorsx), axis=1)
+    constrained_sensors_tuple = np.transpose(constrained_sensors_array)
+    if len(constrained_sensorsx) == 0: ##Check to handle condition when number of sensors in the constrained region = 0
+        idx_constrained = []
+    else:
+        idx_constrained = np.ravel_multi_index(constrained_sensors_tuple, (nx,ny))
+    return idx_constrained
+
+def get_constrained_sensors_indices_linear(x_min, x_max, y_min, y_max,df):
+    """
+    Function for obtaining constrained column indices from already existing linear sensor locations on the grid.
+
+    Parameters
+    ----------
+    x_min: int, lower bound for the x-axis constraint
+    x_max : int, upper bound for the x-axis constraint
+    y_min : int, lower bound for the y-axis constraint
+    y_max : int, upper bound for the y-axis constraint
+    df : pandas.DataFrame, a dataframe containing the features  and samples
+
+    Returns
+    -------
+    idx_constrained : np.darray, shape [No. of constrained locations], array which contains the constrained
+        locations of the grid in terms of column indices of basis_matrix.
+    """
+    x = df['X (m)'].to_numpy()
+    n_features = x.shape[0]
+    y = df['Y (m)'].to_numpy()
+    idx_constrained = []
+    for i in range(n_features):
+        if (x[i] >= x_min and x[i] <= x_max) and (y[i] >= y_min and y[i] <= y_max):
+            idx_constrained.append(i)
+    return idx_constrained