Modifying radii_constraints code and gqr now works as a script for radii constraints

Author: Niharika Karnik

pysensors/optimizers/_gqr.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pysensors
 
 from pysensors.optimizers._qr import QR
 
@@ -110,17 +111,18 @@ def fit(
                 i_piv = np.argmax(dlens_updated)
                 dlen = dlens_updated[i_piv]
             elif self.constraint_option == "radii_constraints":
-                if j == 0:
-                    dlens_old = dlens 
+               
+                if j == 0: 
                     i_piv = np.argmax(dlens)
                     dlen = dlens[i_piv]
+                    dlens_old = dlens
                 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) 
+                    dlens_updated = ps.utils._norm_calc.f_radii_constraint(j,dlens,dlens_old,p,self._nx,self._ny,self._r,self.all_sensorloc, self.nSensors) #( 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)
 
@@ -146,30 +148,30 @@ def fit(
             R[j + 1 :, j] = 0
 
         self.pivots_ = p
-        self.optimality = np.trace(np.real(R))
-        print("The trace(R) = {}".format(self.optimality))
-
+        
         return self
 
 if __name__ == '__main__':
     faces = datasets.fetch_olivetti_faces(shuffle=True)
     X = faces.data
 
-    n_samples, n_features = X.shape
+    # n_samples, n_features = X.shape
+    X_small = X[:,:256]
+    n_samples, n_features = X_small.shape
     print('Number of samples:', n_samples)
     print('Number of features (sensors):', n_features)
 
     # Global centering
-    X = X - X.mean(axis=0)
+    X_small = X_small - X_small.mean(axis=0)
 
     # Local centering
-    X -= X.mean(axis=1).reshape(n_samples, -1)
+    X_small -= X_small.mean(axis=1).reshape(n_samples, -1)
 
     n_row, n_col = 2, 3
     n_components = n_row * n_col
-    image_shape = (64, 64)
-    nx = 64
-    ny = 64
+    image_shape = (16, 16)
+    nx = 16
+    ny = 16
 
     def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):
         '''Function for plotting faces'''
@@ -190,29 +192,29 @@ 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_sensors = 4
     n_modes = 40
-    r = 5
-    dmd = DMD(svd_rank=0,exact=True,opt=False)
-    dmd.fit(X.T)
-    U = dmd.modes.real
-    np.shape(U)
-    max_basis_modes = 200
-
-    model_dmd_unconstrained = ps.SSPOR(n_sensors=n_sensors, basis=ps.basis.Custom(n_basis_modes=n_modes, U=U))
-    model_dmd_unconstrained.fit(X)
-    basis_matrix_dmd = model_dmd_unconstrained.basis_matrix_
-
-    all_sensors_dmd_unconstrained = model_dmd_unconstrained.get_all_sensors()
-    top_sensors_dmd_unconstrained = model_dmd_unconstrained.get_selected_sensors()
-    optimality_dmd = ps.utils._validation.determinant(top_sensors_dmd_unconstrained, n_features, basis_matrix_dmd)
+    r = 2
+    # dmd = DMD(svd_rank=0,exact=True,opt=False)
+    # dmd.fit(X.T)
+    # U = dmd.modes.real
+    # np.shape(U)
+    # max_basis_modes = 200
+
+    # model_dmd_unconstrained = ps.SSPOR(n_sensors=n_sensors, basis=ps.basis.Custom(n_basis_modes=n_modes, U=U))
+    # model_dmd_unconstrained.fit(X)
+    # basis_matrix_dmd = model_dmd_unconstrained.basis_matrix_
+
+    # all_sensors_dmd_unconstrained = model_dmd_unconstrained.get_all_sensors()
+    # top_sensors_dmd_unconstrained = model_dmd_unconstrained.get_selected_sensors()
+    # optimality_dmd = ps.utils._validation.determinant(top_sensors_dmd_unconstrained, n_features, basis_matrix_dmd)
     # print(optimality0)
-    # basis = ps.basis.SVD(n_basis_modes=n_modes)
-    # optimizer  = ps.optimizers.QR()
-    # model = ps.SSPOR(basis = basis, optimizer=optimizer, n_sensors=n_sensors)
-    # model.fit(X)
-    # top_sensors0 = model.get_selected_sensors()
-    # all_sensors = model.get_all_sensors()
+    #basis = ps.basis.SVD(n_basis_modes=n_modes)
+    optimizer  = ps.optimizers.QR()
+    model = ps.SSPOR(optimizer=optimizer, n_sensors=n_sensors)
+    model.fit(X_small)
+    top_sensors0 = model.get_selected_sensors()
+    all_sensors = model.get_all_sensors()
     # basis_matrix0 = model.basis_matrix_
     # optimality0 = ps.utils._validation.determinant(top_sensors0, n_features, basis_matrix0)
     # print(optimality0)
@@ -221,10 +223,11 @@ def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):
     xmax = 64
     ymin = 10
     ymax = 30
-    sensors_constrained = ps.utils._constraints.get_constraind_sensors_indices(xmin,xmax,ymin,ymax,nx,ny,all_sensors_dmd_unconstrained) #Constrained column indices
+    sensors_constrained = ps.utils._constraints.get_constraind_sensors_indices(xmin,xmax,ymin,ymax,nx,ny,all_sensors) #Constrained column indices
     # didx = np.isin(all_sensors,sensors_constrained,invert=False)
     # const_index = np.nonzero(didx)
     # j =
+    
 
 
     ##Plotting the constrained region
@@ -241,49 +244,51 @@ 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 = "max_n_const_sensors",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()
-    # basis_matrix = model1.basis_matrix_
-    # top_sensors = model1.get_selected_sensors()
-    # print(top_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( optimizer = optimizer1, n_sensors = n_sensors)
+    model1.fit(X_small)
+    all_sensors1 = model1.get_all_sensors()
+    basis_matrix = model1.basis_matrix_
+    top_sensors = model1.get_selected_sensors()
+    print(top_sensors)
     # optimality = ps.utils._validation.determinant(top_sensors, n_features, basis_matrix)
     # print(optimality)
-    optimizer_dmd_constrained = ps.optimizers.GQR(sensors_constrained,n_sensors,n_const_sensors,all_sensors_dmd_unconstrained,constraint_option = "exact_n_const_sensors",nx = nx, ny = ny, r = r)
-    model_dmd_constrained = ps.SSPOR(n_sensors=n_sensors, basis=ps.basis.Custom(n_basis_modes=n_modes, U=U), optimizer = optimizer_dmd_constrained)
-    model_dmd_constrained.fit(X)
-    all_sensors_dmd_constrained = model_dmd_constrained.get_all_sensors()
-
-    top_sensors_dmd_constrained = model_dmd_constrained.get_selected_sensors()
-    basis_matrix_dmd_constrained = model_dmd_constrained.basis_matrix_
-    optimality = ps.utils._validation.determinant(top_sensors_dmd_constrained, n_features, basis_matrix_dmd_constrained)
+    # optimizer_dmd_constrained = ps.optimizers.GQR(sensors_constrained,n_sensors,n_const_sensors,all_sensors_dmd_unconstrained,constraint_option = "exact_n_const_sensors",nx = nx, ny = ny, r = r)
+    # model_dmd_constrained = ps.SSPOR(n_sensors=n_sensors, basis=ps.basis.Custom(n_basis_modes=n_modes, U=U), optimizer = optimizer_dmd_constrained)
+    # model_dmd_constrained.fit(X)
+    # all_sensors_dmd_constrained = model_dmd_constrained.get_all_sensors()
+
+    # top_sensors_dmd_constrained = model_dmd_constrained.get_selected_sensors()
+    # basis_matrix_dmd_constrained = model_dmd_constrained.basis_matrix_
+    # optimality = ps.utils._validation.determinant(top_sensors_dmd_constrained, n_features, basis_matrix_dmd_constrained)
     # print(optimality)
 
     ## TODO: this can be done using ravel and unravel more elegantly
     #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_dmd_constrained] = 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
-    # 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)
+    # img[top_sensors_dmd_constrained] = 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
+    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/_norm_calc.py

@@ -109,48 +109,61 @@ def predetermined_norm_calc(lin_idx, dlens, piv, j, n_const_sensors, n_sensors):
         dlens[didx] = 0
     return dlens
 
-def f_radii_constraint(j,dlens,dlens_old,piv,nx,ny,r):
-    a = np.unravel_index(piv, (nx,ny))
-    n_features = len(piv)
+def f_radii_constraint(j,dlens,dlens_old,piv,nx,ny,r, all_sensors, n_sensors):
     if j == 1:
-        x_cord = a[0][j-1]
-        y_cord = a[1][j-1]
-        #print(x_cord, y_cord)
-        constrained_sensorsx = []
-        constrained_sensorsy = []
-        for i in range(n_features):
-            if ((a[0][i]-x_cord)**2 + (a[1][i]-y_cord)**2) < r**2: 
-                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)
-        idx_constrained = np.ravel_multi_index(constrained_sensors_tuple, (nx,ny))
-#         print(idx_constrained)
+        idx_constrained = get_constraind_sensors_indices_radii(j,piv,r, nx,ny, all_sensors)
+        print(idx_constrained)
         didx = np.isin(piv[j:],idx_constrained,invert= False)
         dlens[didx] = 0
         return dlens
-    else: 
+    else:
         result = np.where(dlens_old == 0)[0]
         result_list = result.tolist()
-        result_list = [x - 1 for x in result_list]
+        result_list = [x + (j-1) for x in result_list]
         result_array = np.array(result_list)
-        x_cord = a[0][j-1]
-        y_cord = a[1][j-1]
-        #print(x_cord, y_cord)
-        constrained_sensorsx = []
-        constrained_sensorsy = []
-        for i in range(n_features):
-            if ((a[0][i]-x_cord)**2 + (a[1][i]-y_cord)**2) < r**2: 
-                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)
-        idx_constrained = np.ravel_multi_index(constrained_sensors_tuple, (nx,ny))
-        t = np.concatenate((idx_constrained,result_array), axis = 0)
+        print(result_array)
+       
+        idx_constrained1 = get_constraind_sensors_indices_radii(j,piv,r, nx,ny, all_sensors)
+        t = np.concatenate((idx_constrained1,result_array), axis = 0)
         didx = np.isin(piv[j:],t,invert= False)
         dlens[didx] = 0
-        return dlens
+        return dlens
+
+def get_constraind_sensors_indices_radii(j,piv,r, nx,ny, all_sensors):
+    """
+    Function for mapping constrained sensor locations on the grid with the column indices of the basis_matrix.
+
+    Parameters
+        ----------
+        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 locationsof 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(piv, (nx,ny))
+    t = np.unravel_index(all_sensors, (nx,ny))
+    x_cord = a[0][j-1]
+    y_cord = a[1][j-1]
+    #print(x_cord,y_cord)
+    constrained_sensorsx = []
+    constrained_sensorsy = []
+    for i in range(n_features):
+        if ((t[0][i]-x_cord)**2 + (t[1][i]-y_cord)**2) < r**2:
+            constrained_sensorsx.append(t[0][i])
+            constrained_sensorsy.append(t[1][i])
+
+    constrained_sensorsx = np.array(constrained_sensorsx)
+    constrained_sensorsy = np.array(constrained_sensorsy)
+    constrained_sensors_array = np.stack((constrained_sensorsx, constrained_sensorsy), 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
+