diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8273c83 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +tests/__pycache__/* +*.pyc diff --git a/src/graph_laplacian.py b/src/graph_laplacian.py index d83917f..6e33108 100644 --- a/src/graph_laplacian.py +++ b/src/graph_laplacian.py @@ -1,14 +1,15 @@ -# tools +# tools import numpy as np import scipy import torch -from scipy.sparse import lil_matrix, block_diag,csr_array,diags,csr_matrix +from scipy.sparse import lil_matrix, block_diag, csr_array, diags, csr_matrix def create_laplacian_matrix(nx, ny=None): """ - Helper method to create the laplacian matrix for the laplacian regularization - + Helper method to create the laplacian matrix for the laplacian + regularization + Parameters ---------- :param nx: height of the original image @@ -24,23 +25,23 @@ def create_laplacian_matrix(nx, ny=None): """ if ny is None: ny = nx - assert(nx>1) - assert(ny>1) - #Blocks corresponding to the corner of the image (linking row elements) - top_block=lil_matrix((ny,ny),dtype=np.float32) - top_block.setdiag([2]+[3]*(ny-2)+[2]) - top_block.setdiag(-1,k=1) - top_block.setdiag(-1,k=-1) - #Blocks corresponding to the middle of the image (linking row elements) - mid_block=lil_matrix((ny,ny),dtype=np.float32) - mid_block.setdiag([3]+[4]*(ny-2)+[3]) - mid_block.setdiag(-1,k=1) - mid_block.setdiag(-1,k=-1) - #Construction of the diagonal of blocks - list_blocks=[top_block]+[mid_block]*(nx-2)+[top_block] - blocks=block_diag(list_blocks) - #Diagonals linking different rows - blocks.setdiag(-1,k=ny) + assert (nx > 1) + assert (ny > 1) + # Blocks corresponding to the corner of the image (linking row elements) + top_block = lil_matrix((ny, ny), dtype=np.float32) + top_block.setdiag([2] + [3] * (ny - 2) + [2]) + top_block.setdiag(-1, k=1) + top_block.setdiag(-1, k=-1) + # Blocks corresponding to the middle of the image (linking row elements) + mid_block = lil_matrix((ny, ny), dtype=np.float32) + mid_block.setdiag([3] + [4]*(ny - 2) + [3]) + mid_block.setdiag(-1, k=1) + mid_block.setdiag(-1, k=-1) + # Construction of the diagonal of blocks + list_blocks = [top_block] + [mid_block]*(nx-2) + [top_block] + blocks = block_diag(list_blocks) + # Diagonals linking different rows + blocks.setdiag(-1, k=ny) return blocks @@ -64,7 +65,7 @@ def delete_from_csr(mat, row_indices=[], col_indices=[]): row_mask[rows] = False col_mask = np.ones(mat.shape[1], dtype=bool) col_mask[cols] = False - return mat[row_mask][:,col_mask] + return mat[row_mask][:, col_mask] elif len(rows) > 0: mask = np.ones(mat.shape[0], dtype=bool) mask[rows] = False @@ -72,11 +73,12 @@ def delete_from_csr(mat, row_indices=[], col_indices=[]): elif len(cols) > 0: mask = np.ones(mat.shape[1], dtype=bool) mask[cols] = False - return mat[:,mask] + return mat[:, mask] else: return mat - -def remove_vertex(L,lst_rows=[],lst_cols=[]): + + +def remove_vertex(L, lst_rows=[], lst_cols=[]): """ Function that removes a vertex and adjust the graph laplacian matrix. """ @@ -84,76 +86,78 @@ def remove_vertex(L,lst_rows=[],lst_cols=[]): L_cut = L_cut - diags(L_cut.diagonal()) L_cut = L_cut - diags(L_cut.sum(axis=1).A1) assert (L_cut.sum(axis=1).A1 == np.zeros(L_cut.shape[0])).all() - return L_cut -def laplacian(Y,nx,ny=None): + +def laplacian(Y, nx, ny=None): """ Function that removes the vertices corresponding to Nan locations of tensor Y. Args: - Y: torch.tensor of the observations (Float64) - nx,ny: dimensions of the image (Int64) """ - # find Nan indices - Nan_indices = torch.where(torch.isnan(Y.ravel())==True)[0] + # find Nan indices + Nan_indices = torch.where(torch.isnan(Y.ravel()) is True)[0] # get list of indices list_idx = list(Nan_indices.detach().numpy()) # create Laplacian L = create_laplacian_matrix(nx, ny=ny) - L = remove_vertex(L,list_idx,list_idx) + L = remove_vertex(L, list_idx, list_idx) return L + def unnormalized_laplacian(y, nx, ny=None, method='inverse'): """Construct numpy array with non zeros weights and non zeros indices. - - Args: + + Args: - y: np.array for observations of size (size_x*size_y) - method: str indicating how to compute the weight of the unnormalized laplacian """ - # create laplacian matrix - lapl_tmp = laplacian(y,nx,ny) + # create laplacian matrix + lapl_tmp = laplacian(y, nx, ny) lapl_tmp.setdiag(np.zeros(nx*ny)) - + # select the non nan indices - y_tmp = y[torch.isnan(y)==False] + y_tmp = y[torch.isnan(y) is False] # store non zero indices idx_rows = np.array(lapl_tmp.nonzero()[0]) idx_cols = np.array(lapl_tmp.nonzero()[1]) - + # construct the set of weights nnz_w = np.zeros_like(idx_rows, dtype=np.float32) # construction of the non zeros weights if method == 'inverse': - nnz_w = 1/(np.abs(y_tmp[idx_rows] - y_tmp[idx_cols])+1e-4) + nnz_w = 1/(np.abs(y_tmp[idx_rows] - y_tmp[idx_cols]) + 1e-4) else: nnz_w = np.exp(-np.abs(y_tmp[idx_rows] - y_tmp[idx_cols])) - - # construct the non diagonal terms of the Laplacian + + # construct the non diagonal terms of the Laplacian lapl_nondiag = csr_array((nnz_w, (idx_rows, idx_cols)), shape=(lapl_tmp.shape[0], lapl_tmp.shape[0]), dtype=np.float32) # construct the diagonal terms of the Laplacian lapl_diag = diags(lapl_nondiag.sum(axis=0)) - + # construct the Laplacian L = lapl_diag - lapl_nondiag - + return L + def laplacian_chain(nb_vertices): """ Construct the Laplacian matrix of a chain. """ - L = np.zeros((nb_vertices,nb_vertices)) - - # First vertex - L[0,0] = 1 - L[0,1] = -1 - + L = np.zeros((nb_vertices, nb_vertices)) + + # First vertex + L[0, 0] = 1 + L[0, 1] = -1 + # Toeplitz matrix if nb_vertices > 2: first_row = torch.zeros(nb_vertices) @@ -165,34 +169,34 @@ def laplacian_chain(nb_vertices): first_col[0] = -1 D = scipy.linalg.toeplitz(first_col, r=first_row) - - L[1:nb_vertices-1,:] = D - + + L[1:nb_vertices-1, :] = D + # Last vertex - L[-1,-2] = -1 - L[-1,-1] = 1 - + L[-1, -2] = -1 + L[-1, -1] = 1 + return L + def unnormalized_laplacian_chain(y, nx, method='inverse'): """Construct numpy array with non zeros weights and non zeros indices. - - Args: + + Args: - y: np.array for aggregated observations of size (nb_bins) - method: str indicating how to compute the weight of the unnormalized laplacian """ - # create laplacian matrix + # create laplacian matrix lapl_tmp = csr_matrix(laplacian_chain(nx)) lapl_tmp.setdiag(np.zeros(nx*nx)) - + # select the non nan indices y_tmp = np.nan_to_num(y) - + # store non zero indices idx_rows = np.array(lapl_tmp.nonzero()[0]) idx_cols = np.array(lapl_tmp.nonzero()[1]) - - + # construct the set of weights nnz_w = np.zeros_like(idx_rows, dtype=np.float32) @@ -201,15 +205,14 @@ def unnormalized_laplacian_chain(y, nx, method='inverse'): nnz_w = 1/(np.abs(y_tmp[idx_rows] - y_tmp[idx_cols])) else: nnz_w = np.exp(-np.abs(y_tmp[idx_rows] - y_tmp[idx_cols])) - - - # construct the non diagonal terms of the Laplacian + + # construct the non diagonal terms of the Laplacian lapl_nondiag = csr_array((nnz_w, (idx_rows, idx_cols)), shape=(lapl_tmp.shape[0], lapl_tmp.shape[0]), dtype=np.float32) # construct the diagonal terms of the Laplacian lapl_diag = diags(lapl_nondiag.sum(axis=0)) - + # construct the Laplacian L = lapl_diag - lapl_nondiag - - return L \ No newline at end of file + + return L diff --git a/tests/test_background.py b/tests/test_background.py index f323dea..1513871 100644 --- a/tests/test_background.py +++ b/tests/test_background.py @@ -1,8 +1,6 @@ import unittest import torch import numpy as np - - import os import sys maindir = os.getcwd() @@ -11,7 +9,6 @@ sys.path.append(main_path+"/ds4ms/code/src") from background import background - class TestBackground(unittest.TestCase): def setUp(self): @@ -19,7 +16,7 @@ class TestBackground(unittest.TestCase): self.bg.load_data(verbose=True) self.bg.set_grid_volume(dqx=0.03, dqy=0.03, dE=0.08) self.bg.set_binned_data() - self.bg.set_radial_bins(max_radius=6.0, n_bins=10 + self.bg.set_radial_bins(max_radius=6.0, n_bins=10) self.bg.Ygrid = torch.tensor(np.random.rand(10, 10), dtype=torch.float64) def test_load_data(self): @@ -166,12 +163,14 @@ class TestBackground(unittest.TestCase): alpha_range = torch.tensor([1.0]) beta_range = torch.tensor([1.0]) mu_range = torch.tensor([1.0]) - result = self.bg.cross_validation(lambda_range=lambda_range, alpha_range=alpha_range, beta_range=beta_range, mu_range=mu_range, n_epochs=1, verbose=False) + result = self.bg.cross_validation(lambda_range=lambda_range, alpha_range=alpha_range, + beta_range=beta_range, mu_range=mu_range, n_epochs=1, verbose=False) self.assertIsNotNone(result) def test_compute_mask(self): result = self.bg.compute_mask(q=0.75, e_cut=None) self.assertIsNotNone(result) + if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() diff --git a/tests/test_graph_laplacian.py b/tests/test_graph_laplacian.py index e75a164..994204d 100644 --- a/tests/test_graph_laplacian.py +++ b/tests/test_graph_laplacian.py @@ -2,14 +2,15 @@ import unittest import numpy as np import torch from scipy.sparse import lil_matrix, csr_matrix, coo_matrix - - import os import sys maindir = os.getcwd() main_path = maindir[:maindir.find('ds4ms/code')] sys.path.append(main_path+"/ds4ms/code/src") -from graph_laplacian import create_laplacian_matrix, delete_from_csr, remove_vertex, laplacian, unnormalized_laplacian, laplacian_chain, unnormalized_laplacian_chain +from graph_laplacian import create_laplacian_matrix, delete_from_csr, \ + remove_vertex, laplacian, unnormalized_laplacian, laplacian_chain, \ + unnormalized_laplacian_chain + class TestGraphLaplacian(unittest.TestCase): @@ -66,5 +67,6 @@ class TestGraphLaplacian(unittest.TestCase): self.assertEqual(L.shape, (nx, nx)) self.assertIsInstance(L, csr_matrix) + if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main()