Moving norm_calculation functions to utilities

Author: niharika2999

pysensors/optimizers/_gqr.py

@@ -8,6 +8,7 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 import pysensors as ps
+from matplotlib.patches import Circle
 
 
 class GQR(QR):
@@ -27,7 +28,7 @@ class GQR(QR):
 
     @ authors: Niharika Karnik (@nkarnik2999), Mohammad Abdo (@Jimmy-INL), and Krithika Manohar (@kmanohar)
     """
-    def __init__(self,idx_constrained,n_sensors,n_const_sensors,all_sensors,constraint_option):
+    def __init__(self,idx_constrained,n_sensors,n_const_sensors,all_sensors,constraint_option,nx,ny,r):
         """
         Attributes
         ----------
@@ -47,11 +48,15 @@ def __init__(self,idx_constrained,n_sensors,n_const_sensors,all_sensors,constrai
         """
         self.pivots_ = None
         self.optimality = None
+        
         self.constrainedIndices = idx_constrained
         self.nSensors = n_sensors
         self.nConstrainedSensors = n_const_sensors
         self.all_sensorloc = all_sensors
         self.constraint_option = constraint_option
+        self._nx = nx
+        self._ny = ny
+        self._r = r
 
     def fit(
         self,
@@ -92,16 +97,33 @@ def fit(
             # Norm of each column
             dlens = np.sqrt(np.sum(np.abs(r) ** 2, axis=0))
             if self.constraint_option == "max_n_const_sensors" :
-                dlens_updated = norm_calc_max_n_const_sensors(self.constrainedIndices,dlens,p,j, self.nConstrainedSensors,self.all_sensorloc,self.nSensors) 
+                dlens_updated = ps.utils._norm_calc.norm_calc_max_n_const_sensors(self.constrainedIndices,dlens,p,j, self.nConstrainedSensors,self.all_sensorloc,self.nSensors) 
+                i_piv = np.argmax(dlens_updated)
+                dlen = dlens_updated[i_piv]
             elif self.constraint_option == "exact_n_const_sensors" : 
-                dlens_updated = norm_calc_exact_n_const_sensors(self.constrainedIndices,dlens,p,j,self.nConstrainedSensors)
+                dlens_updated = ps.utils._norm_calc.norm_calc_exact_n_const_sensors(self.constrainedIndices,dlens,p,j,self.nConstrainedSensors)
+                i_piv = np.argmax(dlens_updated)
+                dlen = dlens_updated[i_piv]
             elif self.constraint_option == "predetermined_end":
-                dlens_updated = predetermined_norm_calc(self.constrainedIndices, dlens, p, j, self.nConstrainedSensors, self.nSensors)
+                dlens_updated = ps.utils._norm_calc.predetermined_norm_calc(self.constrainedIndices, dlens, p, j, self.nConstrainedSensors, self.nSensors)
+                i_piv = np.argmax(dlens_updated)
+                dlen = dlens_updated[i_piv]
+            elif self.constraint_option == "radii_constraints":
+                if j == 0:
+                    dlens_old = dlens 
+                    i_piv = np.argmax(dlens)
+                    dlen = dlens[i_piv]
+                else:
+
+                    dlens_updated = ps.utils._norm_calc.f_radii_constraint(j,dlens,dlens_old,p,self._nx,self._ny,self._r) #( self.radius,self._nx,self._ny,self.all_sensorloc,dlens,p,j) 
+                    i_piv = np.argmax(dlens_updated)
+                    dlen = dlens_updated[i_piv]
+                    dlens_old = dlens_updated
 
             # Choose pivot
-            i_piv = np.argmax(dlens_updated)
+            # i_piv = np.argmax(dlens_updated)
 
-            dlen = dlens_updated[i_piv]
+            # dlen = dlens_updated[i_piv]
 
             if dlen > 0:
                 u = r[:, i_piv] / dlen
@@ -128,116 +150,6 @@ def fit(
 
         return self
 
-## TODO: why not a part of the class?
-
-#function for mapping sensor locations with constraints
-def norm_calc_exact_n_const_sensors(lin_idx, dlens, piv, j, n_const_sensors): ##Will first force sensors into constrained region
-    #num_sensors should be fixed for each custom constraint (for now)
-    #num_sensors must be <= size of constraint region
-    """
-    Function for mapping constrained sensor locations with the QR procedure.
-
-    Parameters
-        ----------
-        lin_idx: np.ndarray, shape [No. of constrained locations]
-            Array which contains the constrained locationsof the grid in terms of column indices of basis_matrix.
-        dlens: np.ndarray, shape [Variable based on j]
-            Array which contains the norm of columns of basis matrix.
-        piv: np.ndarray, shape [n_features]
-            Ranked list of sensor locations.
-        n_const_sensors: int,
-            Number of sensors to be placed in the constrained area.
-        j: int,
-            Iterative variable in the QR algorithm.
-
-        Returns
-        -------
-        dlens : np.darray, shape [Variable based on j] with constraints mapped into it.
-    """
-    if j < n_const_sensors: # force sensors into constraint region
-        #idx = np.arange(dlens.shape[0])
-        #dlens[np.delete(idx, lin_idx)] = 0
-
-        didx = np.isin(piv[j:],lin_idx,invert=True)
-        dlens[didx] = 0
-    else:
-        didx = np.isin(piv[j:],lin_idx,invert=False)
-        dlens[didx] = 0
-    return dlens
-
-def norm_calc_max_n_const_sensors(lin_idx, dlens, piv, j, const_sensors,all_sensors,n_sensors): ##Optimal sensor placement with constraints (will place sensors in  the order of QR)
-    """
-    Function for mapping constrained sensor locations with the QR procedure (Optimally).
-
-    Parameters
-        ----------
-        lin_idx: np.ndarray, shape [No. of constrained locations]
-            Array which contains the constrained locationsof the grid in terms of column indices of basis_matrix.
-        dlens: np.ndarray, shape [Variable based on j]
-            Array which contains the norm of columns of basis matrix.
-        piv: np.ndarray, shape [n_features]
-            Ranked list of sensor locations.
-        j: int,
-            Iterative variable in the QR algorithm.
-        const_sensors: int,
-            Number of sensors to be placed in the constrained area.
-        all_sensors: np.ndarray, shape [n_features]
-            Ranked list of sensor locations.
-        n_sensors: integer,
-            Total number of sensors
-
-        Returns
-        -------
-        dlens : np.darray, shape [Variable based on j] with constraints mapped into it.
-    """
-    counter = 0
-    mask = np.isin(all_sensors,lin_idx,invert=False)
-    const_idx = all_sensors[mask]
-    updated_lin_idx = const_idx[const_sensors:]
-    for i in range(n_sensors):
-        if np.isin(all_sensors[i],lin_idx,invert=False):
-            counter += 1
-            if counter < const_sensors:
-                dlens = dlens
-            else:
-                didx = np.isin(piv[j:],updated_lin_idx,invert=False)
-                dlens[didx] = 0
-    return dlens
-
-def predetermined_norm_calc(lin_idx, dlens, piv, j, n_const_sensors, n_sensors):
-    """
-    Function for mapping constrained sensor locations with the QR procedure.
-
-    Parameters
-        ----------
-        lin_idx: np.ndarray, shape [No. of constrained locations]
-            Array which contains the constrained locationsof the grid in terms of column indices of basis_matrix.
-        dlens: np.ndarray, shape [Variable based on j]
-            Array which contains the norm of columns of basis matrix.
-        piv: np.ndarray, shape [n_features]
-            Ranked list of sensor locations.
-        n_const_sensors: int,
-            Number of sensors to be placed in the constrained area.
-        j: int,
-            Iterative variable in the QR algorithm.
-
-        Returns
-        -------
-        dlens : np.darray, shape [Variable based on j] with constraints mapped into it.
-    """
-    if (n_sensors - n_const_sensors) <= j <= n_sensors: # force sensors into constraint region
-        #idx = np.arange(dlens.shape[0])
-        #dlens[np.delete(idx, lin_idx)] = 0
-
-        didx = np.isin(piv[j:],lin_idx,invert=True)
-        dlens[didx] = 0
-    else:
-        didx = np.isin(piv[j:],lin_idx,invert=False)
-        dlens[didx] = 0
-    return dlens
-
-
-
 if __name__ == '__main__':
     faces = datasets.fetch_olivetti_faces(shuffle=True)
     X = faces.data
@@ -277,8 +189,9 @@ def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):
     #Find all sensor locations using built in QR optimizer
     #max_const_sensors = 230
     n_const_sensors = 3
-    n_sensors = 10
-    n_modes = 11
+    n_sensors = 30
+    n_modes = 50
+    r = 5
     basis = ps.basis.SVD(n_basis_modes=n_modes)
     optimizer  = ps.optimizers.QR()
     model = ps.SSPOR(basis = basis, optimizer=optimizer, n_sensors=n_sensors)
@@ -311,7 +224,7 @@ def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):
     # plt.title('Constrained region');
 
     ## Fit the dataset with the optimizer GQR
-    optimizer1 = GQR(sensors_constrained,n_sensors,n_const_sensors,all_sensors, constraint_option = "exact_n_const_sensors")
+    optimizer1 = GQR(sensors_constrained,n_sensors,n_const_sensors,all_sensors, constraint_option = "radii_constraints",nx = nx, ny = ny, r = r)
     model1 = ps.SSPOR(basis = basis, optimizer = optimizer1, n_sensors = n_sensors)
     model1.fit(X)
     all_sensors1 = model1.get_all_sensors()
@@ -322,13 +235,26 @@ def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):
     #yConstrained = np.floor(top_sensors[:n_const_sensors]/np.sqrt(n_features))
     #xConstrained = np.mod(top_sensors[:n_const_sensors],np.sqrt(n_features))
 
+    # img = np.zeros(n_features)
+    # img[top_sensors] = 16
+    # #plt.plot(xConstrained,yConstrained,'*r')
+    # plt.plot([xmin,xmin],[ymin,ymax],'r')
+    # plt.plot([xmin,xmax],[ymax,ymax],'r')
+    # plt.plot([xmax,xmax],[ymin,ymax],'r')
+    # plt.plot([xmin,xmax],[ymin,ymin],'r')
+    # plt.imshow(img.reshape(image_shape),cmap=plt.cm.binary)
+    # plt.title('n_sensors = {}, n_constr_sensors = {}'.format(n_sensors,n_const_sensors))
+    # plt.show()
+
     img = np.zeros(n_features)
     img[top_sensors] = 16
-    #plt.plot(xConstrained,yConstrained,'*r')
-    plt.plot([xmin,xmin],[ymin,ymax],'r')
-    plt.plot([xmin,xmax],[ymax,ymax],'r')
-    plt.plot([xmax,xmax],[ymin,ymax],'r')
-    plt.plot([xmin,xmax],[ymin,ymin],'r')
-    plt.imshow(img.reshape(image_shape),cmap=plt.cm.binary)
-    plt.title('n_sensors = {}, n_constr_sensors = {}'.format(n_sensors,n_const_sensors))
+    fig,ax = plt.subplots(1)
+    ax.set_aspect('equal')
+    ax.imshow(img.reshape(image_shape),cmap=plt.cm.binary)
+    print(top_sensors)
+    top_sensors_grid = np.unravel_index(top_sensors, (nx,ny))
+    # figure, axes = plt.subplots()
+    for i in range(len(top_sensors_grid[0])):
+        circ = Circle( (top_sensors_grid[1][i], top_sensors_grid[0][i]), r ,color='r',fill = False )
+        ax.add_patch(circ)
     plt.show()

pysensors/utils/__init__.py

@@ -7,6 +7,10 @@
 from ._constraints import functional_constraints
 from ._validation import determinant
 from ._validation import relative_reconstruction_error
+from ._norm_calc import norm_calc_exact_n_const_sensors
+from ._norm_calc import norm_calc_max_n_const_sensors
+from ._norm_calc import predetermined_norm_calc
+from ._norm_calc import f_radii_constraint
 
 __all__ = [
     "constrained_binary_solve",
@@ -17,5 +21,9 @@
     "box_constraints",
     "functional_constraints",
     "determinant",
-    "relative_reconstruction_error"
+    "relative_reconstruction_error",
+    "norm_calc_exact_n_const_sensors",
+    "norm_calc_max_n_const_sensors",
+    "predetermined_norm_calc",
+    "f_radii_constraint"
 ]