diff --git a/README.md b/README.md index 36193d5..6b2cdf6 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,54 @@ # Accuracy Bounds -Install project with +This repository contains the official python-based implementation accompanying the paper: +**“Computable Sharp Accuracy Bounds for Inverse Problems”** ([arXiv:2510.10229](https://arxiv.org/abs/2510.10229)) -"" -pip install -e . -"" +If you use this software in your work, please cite our [preprint](https://arxiv.org/abs/2510.10229): + +```bibtex +@article{gottschling2025average, + title={Average Kernel Sizes--Computable Sharp Accuracy Bounds for Inverse Problems}, + author={Gottschling, Nina M and Iagaru, David and Gawlikowski, Jakob and Sgouralis, Ioannis}, + journal={arXiv preprint arXiv:2510.10229}, + year={2025} +} +``` + +For the definition of the worst-case kernel size, see [preprint](https://arxiv.org/abs/2311.16898): + +```bibtex +@article{gottschling2023existence, + title={On the existence of optimal multi-valued decoders and their accuracy bounds for undersampled inverse problems}, + author={Gottschling, Nina Maria and Campodonico, Paolo and Antun, Vegard and Hansen, Anders C}, + journal={arXiv preprint arXiv:2311.16898}, + year={2023} +} +``` +## Installation +This project has been tested under **Python 3.7** on a Unix development environment. + +### 1. Clone the repository and create a virtual environment and change into the repository folder. +``` +git clone https://github.com/nm19000/AccuracyBounds.git +cd AccuracyBounds +``` + +### 2. Create and activate a virtual environment +``` +python -m venv venv +source venv/bin/activate +``` +### 3. Install the project +#### 3.1 Lightweight Version (without PyTorch) +``` +pip install -e . +``` +#### 3.2 Full version (with PyTorch) +``` +pip install -e .[torch] +``` # Average Kernel Sizes - Computable Sharp Accuracy Bounds for Inverse Problems @@ -14,7 +56,6 @@ The reconstruction of an unknown quantity from noisy measurements is a mathemati ## Accuracy Bounds for Inverse Problems - Computation of worst-case and average kernel size for an inverse problem with noise of the form: $$ @@ -27,38 +68,15 @@ The lower bound to the worst-case error of any approximate inverse map is the wo # Testing -The algorithms for computing the worst-case and average kernel size are tested against linear algebra examples, where the worst-case and average kernel size can be calculated analytically in the limit of infinite datapoints. Thus, we ensure that the implemented algorithms compute the correct quantities. The algorithm versions that run with cuda can be tested with test/test_toy_example_cuda.py. The algorithm versions that only run with numpy can be tested with test/test_toy_example_np.py. +The algorithms for computing the worst-case and average kernel size are tested against linear algebra examples, where the worst-case and average kernel size can be calculated analytically in the limit of infinite datapoints. Thus, we ensure that the implemented algorithms compute the correct quantities. The algorithm versions that run with cuda can be tested with ``test/test_toy_example_torch.py``. The algorithm versions that only run with numpy can be tested with ``test/test_toy_example_np.py``. All tests in the project can be automatically run by calling ``pytest``. For computing the feasible sets there are two versions to compute these available: in terms of list and in terms of feasible appartenance matrices that allocate data points to feasible sets. Both versions are suitable for forward models with additive noise. -For interactive testing the algorithms to compute the worst-case and average kernel size with numpy and cuda please see test/playground/testing.ipynb. - -# References -If you use this software in your work, please cite our [preprint](https://arxiv.org/abs/2510.10229): - -```bibtex -@article{gottschling2025average, - title={Average Kernel Sizes--Computable Sharp Accuracy Bounds for Inverse Problems}, - author={Gottschling, Nina M and Iagaru, David and Gawlikowski, Jakob and Sgouralis, Ioannis}, - journal={arXiv preprint arXiv:2510.10229}, - year={2025} -} -``` - -For the definition of the worst-case kernel size, see [preprint](https://arxiv.org/abs/2311.16898): - -```bibtex -@article{gottschling2023existence, - title={On the existence of optimal multi-valued decoders and their accuracy bounds for undersampled inverse problems}, - author={Gottschling, Nina Maria and Campodonico, Paolo and Antun, Vegard and Hansen, Anders C}, - journal={arXiv preprint arXiv:2311.16898}, - year={2023} -} -``` +For interactive testing the algorithms to compute the worst-case and average kernel size with numpy and cuda please see ``test/playground/testing.ipynb``. ## Experiments from [Paper](https://arxiv.org/abs/2510.10229) - We demonstrate the validity of the algorithms on two inverse problems from different domains: fluorescence localization microscopy and super-resolution of multi-spectral satellite data. The code for generating the data for reproducing the localization microscopy experiments can be found in examples/data/localization_microscopy/data_A/source_code. The data used for the localization microscopy experiments can be found in examples/data/localization_microscopy. To reproduce the tabular data for the figures and tables in this manuscript, the average kernel size and loss computations, can be found in examples/example_localizationmicroscopy.ipynb. The instructions for reproducing the satellite data super-resolution experiments are to be found below. + We demonstrate the validity of the algorithms on two inverse problems from different domains: fluorescence localization microscopy and super-resolution of multi-spectral satellite data. The code for generating the data for reproducing the localization microscopy experiments can be found in ``examples/data/localization_microscopy/data_A/source_code``. The data used for the localization microscopy experiments can be found in examples/data/localization_microscopy. To reproduce the tabular data for the figures and tables in this manuscript, the average kernel size and loss computations, can be found in ``examples/example_localizationmicroscopy.ipynb``. The instructions for reproducing the satellite data super-resolution experiments are to be found below. ### Satellite Data Super Resolution Experiments @@ -104,37 +122,43 @@ The data used in our experiments can be downloaded from the indications on [hugg ### Preliminary Kernel Size Computations Run the command - " python test/S2_SR/Kernelsize_computations.py " - with the desired options specified in the file Kernelsize_computations.py +``` +python test/S2_SR/Kernelsize_computations.py +``` +with the desired options specified in the file ``Kernelsize_computations.py`` ### Preliminary Operator Computations To calculate the downsampling operator and its kernel projection in matrix form, run -" python examples/S2_SR/op_testing.py " +``` +python examples/S2_SR/op_testing.py +``` Manually change the values of the following variables inside the python file to enable or discard the corresponding computation or visualization : +``` plot_sparsity = False # To plot the sparsity pattern of the operators \ check_DSOp = False # To check that the downsampling operator uner matrix is correctly computed \ computeDS = False # To compute the downsampling operator under matrix form \ compute_P_null = False # To compute the kenrel projection \ check_P_null = False # To check that the kernel projection operator uner matrix is correctly computed\ scale_plot = False # To plot some satellite images with the scale bars \ +``` ### Run the experiments After having run the preliminary kernel size and kernel projection operator computations, the experiments are ready to be reproduced with the command -"" +``` python test/S2_SR/experiments.py -"" +``` Adjust the following parameters inside the python file: - +``` DSHR = True # Whether the lower resolution image is the downsampled version of the high resolution image (we run the experiments with DSHR = true)\ -light_loading = False # Whether you use the light dataloader or you want to use stored patches. We run the experiments with the parameter set to false, but it is recommended for more memory and speed efficiency. Warning : it has to correspond with the value of the --light_load parameter in the kernel size computations . If light_loading is set to false, you will need to generate a dataset where each file corresponds to a patch. This can be done using the function build_S2_patched_dataset_DSHR or build_S2_patched_dataset in the utils.py file. Warning such a patched dataset may contain more than 100 000 patches files for 119 full sized images. It is therefore not recommended. If activated, use the dataset SRDataset_perimg_lightload instead of the dataset SRDataset_perimg \ +light_loading = False # Whether you use the light dataloader or you want to use stored patches. We run the experiments with the parameter set to false, but it is recommended for more memory and speed efficiency. Warning : it has to correspond with the value of the ``--light_load`` parameter in the kernel size computations . If light_loading is set to false, you will need to generate a dataset where each file corresponds to a patch. This can be done using the function ``build_S2_patched_dataset_DSHR`` or ``build_S2_patched_dataset`` in the ``utils_torch.py`` file. Warning such a patched dataset may contain more than 100 000 patches files for 119 full sized images. It is therefore not recommended. If activated, use the dataset SRDataset_perimg_lightload instead of the dataset SRDataset_perimg \ PS_X = 16 # Patch size in high resolution (has to correspond with the kernel size computations) \ PS_Y = PS_X//4 \ p_norm = 2 # Defines the used norm among the $L^p$ norms \ @@ -143,7 +167,7 @@ noise_level_KS = 4000 # has to Correspond with the preliminary computations of t preload_feas_info = True # Preload or not the feasible information from the feasible appartenance matrix. It has to be activated the first time so that the json file grouping the essential information from the feasible appartenance matrix can be saved (it may take more than 30 minutes for 100 000 patches in the dataset). \ pred_type = 'bicub' # The model used for the predictions ('diff' for the opensr-modelm 'bilin' for bilinear interpolation and 'bicub' for bicubic interpolation) \ p = 2 # The p parameter for kernel size and the loss computations - +``` Modify the paths root_folder, feas_app_lightl_path etc manually in the file, according the needs. diff --git a/pyproject.toml b/pyproject.toml index d03e1f7..4f795e6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,12 +17,16 @@ dependencies = [ "click", "pdoc", "numpy", - "torch", "tqdm", "joblib", "scipy", ] +[project.optional-dependencies] +torch = [ + "torch", +] + [build-system] requires = [ "setuptools", diff --git a/src/accuracy_bounds/inverseproblems/utils.py b/src/accuracy_bounds/inverseproblems/utils.py index 997e28c..a5073a3 100644 --- a/src/accuracy_bounds/inverseproblems/utils.py +++ b/src/accuracy_bounds/inverseproblems/utils.py @@ -1,87 +1,4 @@ import numpy as np -from multiprocessing import Pool, cpu_count -from tqdm import tqdm -from scipy.sparse import coo_matrix -from scipy.sparse.linalg import svds -import torch -import scipy.sparse as sp - - - -#To calculate operators in the form of matrices (suits better to big operators) -class MatrixOpCalculator: - def __init__(self, n_in, n_out, Operator, num_workers = None, singular_threshold_ratio = 0.001): - self.n_in = n_in - self.n_out = n_out - self.operator = Operator - self.num_workers=num_workers - if num_workers is None: - self.num_workers = cpu_count() - self.singular_threshold_ratio = singular_threshold_ratio - - # Compute column i of the sparse matrix - def compute_column(self,i): - A = self.operator - e_i = np.zeros(self.n_in) - e_i[i] = 1.0 - col = A(e_i) - - # Find non-zero entries - row_idx = np.nonzero(col)[0] - data = col[row_idx] - col_idx = np.full_like(row_idx, i) - - return row_idx, col_idx, data - - def build_sparse_matrix_parallel(self): - n_in = self.n_in - n_out = self.n_out - - with Pool(self.num_workers) as pool: - # Run in parallel - results = list(tqdm(pool.imap(self.compute_column, range(n_in)), total=n_in)) - - # Collect results - row_indices = [] - col_indices = [] - data = [] - - for r, c, d in results: - row_indices.extend(r) - col_indices.extend(c) - data.extend(d) - - # Create sparse matrix in COO format first, then convert - A_sparse = coo_matrix((data, (row_indices, col_indices)), shape=(n_out, n_in)).tocsc() - return A_sparse - - - def get_range_space_basis(self, A_sparse, sigma_threshold_ratio = 0.001): - # Compute a few smallest singular values - p,q = A_sparse.shape - kmax = int(min(p,q)-1) - #t0 = time.time() - umat, sing, vt = svds(A_sparse, k =kmax ) # smallest magnitude - #t1 = time.time() - #print(f'Took {t1-t0:2f} seconds to compute the SVD of A for kmax = {kmax}') - - # Null space basis vectors: columns of vt.T corresponding to near-zero singular values - threshold = sigma_threshold_ratio * np.max(sing) - range_space_basis = vt[sing >= threshold].T - #print(f'Sigma 1 = {np.max(sing)}') - return range_space_basis - - - - def make_null_projection_operator(self,range_basis): - n = range_basis.shape[0] - - #def matvec(x): - # return null_basis @ (null_basis.T @ x) - - #return LinearOperator((n, n), matvec=matvec, dtype=null_basis.dtype) - return np.eye(n)- range_basis.dot(range_basis.T) - # Function to apply matrix transformation A to points def apply_forwardmodel(A, points): @@ -107,147 +24,4 @@ def projection_nullspace(A, x): project_ns = projection_nullspace_operator(A) x_ns = np.dot(project_ns,x) - return x_ns - -def offset_csr_block(local_csr: torch.Tensor, i0: int, j0: int, global_shape): - """ - Place a CSR 'local_csr' block at (row=i0, col=j0) inside a larger CSR tensor of shape 'global_shape'. - - local_csr: torch.sparse_csr_tensor of shape (h, w) - i0, j0: top-left offsets (rows, cols) in the global matrix - global_shape: (H, W) - returns: torch.sparse_csr_tensor of shape (H, W) - """ - assert local_csr.layout == torch.sparse_csr, "local_csr must be CSR" - H, W = global_shape - h, w = local_csr.shape - if not (0 <= i0 <= H - h): - raise ValueError(f"Row offset i0={i0} with block height {h} exceeds H={H}") - if not (0 <= j0 <= W - w): - raise ValueError(f"Col offset j0={j0} with block width {w} exceeds W={W}") - - crow = local_csr.crow_indices() - col = local_csr.col_indices() - val = local_csr.values() - - # shift columns - col_off = col + j0 - - nnz = val.numel() - device = local_csr.device - itype = crow.dtype - - # Build global crow: zeros up to i0, then local crow, then constant nnz afterwards - crow_g = torch.empty(H + 1, dtype=itype, device=device) - if i0 > 0: - crow_g[:i0] = 0 - crow_g[i0:i0 + h + 1] = crow - if i0 + h < H: - crow_g[i0 + h + 1:] = nnz - - return torch.sparse_csr_tensor(crow_g, col_off, val, size=global_shape, - device=device, dtype=local_csr.dtype) - - - -def insert_no_overlap_keep_A(A_coo, B_csr): - A = A_coo.coalesce() - B = B_csr.to_sparse_coo() - I_A, V_A = A.indices(), A.values() - I_B, V_B = B.indices(), B.values() - # linearize indices to detect overlaps - lin_A = I_A[0] * A.shape[1] + I_A[1] - lin_B = I_B[0] * A.shape[1] + I_B[1] - # keep only B entries not in A - mask = ~torch.isin(lin_B, lin_A) - I = torch.cat([I_A, I_B[:, mask]], dim=1) - V = torch.cat([V_A, V_B[mask]], dim=0) - return torch.sparse_coo_tensor(I, V, A.shape, device=A.device, dtype=A.dtype) - -def sparse_block(A, i0, i1, j0, j1, out_layout="coo"): - """ - Return A[i0:i1, j0:j1] as a sparse tensor without densifying. - Works for A in COO or CSR. out_layout: 'coo' or 'csr' - """ - Acoo = A.to_sparse_coo().coalesce() - r, c = Acoo.indices() - v = Acoo.values() - - m = (r >= i0) & (r < i1) & (c >= j0) & (c < j1) - r2 = r[m] - i0 - c2 = c[m] - j0 - v2 = v[m] - - B = torch.sparse_coo_tensor( - torch.stack([r2, c2], dim=0), - v2, - (i1 - i0, j1 - j0), - device=A.device, dtype=A.dtype - ).coalesce() - - return B if out_layout == "coo" else B.to_sparse_csr() - - - - -def torch_csr_to_scipy(A: torch.Tensor) -> sp.csr_matrix: - """"Convert matrix/linear forward model, as torch tensor, to scipy sparse matrix. - Arguments: - - A: torch tensor as matrix/forward model. - Returns: - - Scipy sparse matrix. - """ - assert A.layout == torch.sparse_csr, f"A must be sparse_csr {type(A.layout)}" - m, n = A.shape - indptr = A.crow_indices().detach().cpu().numpy() # length m+1 - indices = A.col_indices().detach().cpu().numpy() # length nnz - data = A.values().detach().cpu().numpy() # length nnz - - # (optional) SciPy prefers int32 for indices - if indices.dtype != np.int32: indices = indices.astype(np.int32, copy=False) - if indptr.dtype != np.int32: indptr = indptr.astype(np.int32, copy=False) - - sp_matrix = sp.csr_matrix((data, indices, indptr), shape=(m, n)) - - return sp_matrix - -def torch_coo_to_scipy(A: torch.Tensor, to='csr'): - A = A.to_sparse_coo().coalesce() - m, n = A.shape - ij = A.indices().detach().cpu().numpy() # shape [2, nnz] - data = A.values().detach().cpu().numpy() - row, col = ij[0], ij[1] - M = sp.coo_matrix((data, (row, col)), shape=(m, n)) - return M.tocsr() if to == 'csr' else M - -def torch_sparse_to_scipy_csr(A: torch.Tensor) -> sp.csr_matrix: - """ - Convert a PyTorch tensor (CSR/COO/dense) to SciPy csr_matrix. - - No densifying for sparse inputs. - - Handles CUDA tensors by moving index/data arrays to CPU. - """ - m, n = A.shape - - # Dense (strided) tensor - if A.layout == torch.strided: - return sp.csr_matrix(A.detach().cpu().numpy()) - - # Try CSR fast-path - try: - indptr = A.crow_indices() # only valid for CSR - indices = A.col_indices() - data = A.values() - except (AttributeError, RuntimeError): - # Not CSR → go through COO safely - Acoo = A.to_sparse_coo().coalesce() - ij = Acoo.indices().detach().cpu().numpy() # [2, nnz] - row, col = ij[0].astype(np.int32, copy=False), ij[1].astype(np.int32, copy=False) - dat = Acoo.values().detach().cpu().numpy() - return sp.coo_matrix((dat, (row, col)), shape=(m, n)).tocsr() - else: - # CSR build (ensure int32 indices for SciPy) - indptr_np = indptr.detach().cpu().numpy().astype(np.int32, copy=False) - indices_np = indices.detach().cpu().numpy().astype(np.int32, copy=False) - data_np = data.detach().cpu().numpy() - return sp.csr_matrix((data_np, indices_np, indptr_np), shape=(m, n)) - \ No newline at end of file + return x_ns \ No newline at end of file diff --git a/src/accuracy_bounds/inverseproblems/utils_torch.py b/src/accuracy_bounds/inverseproblems/utils_torch.py new file mode 100644 index 0000000..14115df --- /dev/null +++ b/src/accuracy_bounds/inverseproblems/utils_torch.py @@ -0,0 +1,226 @@ +import numpy as np +from multiprocessing import Pool, cpu_count +from tqdm import tqdm +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import svds +import torch +import scipy.sparse as sp + + +#To calculate operators in the form of matrices (suits better to big operators) +class MatrixOpCalculator: + def __init__(self, n_in, n_out, Operator, num_workers = None, singular_threshold_ratio = 0.001): + self.n_in = n_in + self.n_out = n_out + self.operator = Operator + self.num_workers=num_workers + if num_workers is None: + self.num_workers = cpu_count() + self.singular_threshold_ratio = singular_threshold_ratio + + # Compute column i of the sparse matrix + def compute_column(self,i): + A = self.operator + e_i = np.zeros(self.n_in) + e_i[i] = 1.0 + col = A(e_i) + + # Find non-zero entries + row_idx = np.nonzero(col)[0] + data = col[row_idx] + col_idx = np.full_like(row_idx, i) + + return row_idx, col_idx, data + + def build_sparse_matrix_parallel(self): + n_in = self.n_in + n_out = self.n_out + + with Pool(self.num_workers) as pool: + # Run in parallel + results = list(tqdm(pool.imap(self.compute_column, range(n_in)), total=n_in)) + + # Collect results + row_indices = [] + col_indices = [] + data = [] + + for r, c, d in results: + row_indices.extend(r) + col_indices.extend(c) + data.extend(d) + + # Create sparse matrix in COO format first, then convert + A_sparse = coo_matrix((data, (row_indices, col_indices)), shape=(n_out, n_in)).tocsc() + return A_sparse + + + def get_range_space_basis(self, A_sparse, sigma_threshold_ratio = 0.001): + # Compute a few smallest singular values + p,q = A_sparse.shape + kmax = int(min(p,q)-1) + #t0 = time.time() + umat, sing, vt = svds(A_sparse, k =kmax ) # smallest magnitude + #t1 = time.time() + #print(f'Took {t1-t0:2f} seconds to compute the SVD of A for kmax = {kmax}') + + # Null space basis vectors: columns of vt.T corresponding to near-zero singular values + threshold = sigma_threshold_ratio * np.max(sing) + range_space_basis = vt[sing >= threshold].T + #print(f'Sigma 1 = {np.max(sing)}') + return range_space_basis + + + + def make_null_projection_operator(self,range_basis): + n = range_basis.shape[0] + + #def matvec(x): + # return null_basis @ (null_basis.T @ x) + + #return LinearOperator((n, n), matvec=matvec, dtype=null_basis.dtype) + return np.eye(n)- range_basis.dot(range_basis.T) + + +def offset_csr_block(local_csr: torch.Tensor, i0: int, j0: int, global_shape): + """ + Place a CSR 'local_csr' block at (row=i0, col=j0) inside a larger CSR tensor of shape 'global_shape'. + + local_csr: torch.sparse_csr_tensor of shape (h, w) + i0, j0: top-left offsets (rows, cols) in the global matrix + global_shape: (H, W) + returns: torch.sparse_csr_tensor of shape (H, W) + """ + assert local_csr.layout == torch.sparse_csr, "local_csr must be CSR" + H, W = global_shape + h, w = local_csr.shape + if not (0 <= i0 <= H - h): + raise ValueError(f"Row offset i0={i0} with block height {h} exceeds H={H}") + if not (0 <= j0 <= W - w): + raise ValueError(f"Col offset j0={j0} with block width {w} exceeds W={W}") + + crow = local_csr.crow_indices() + col = local_csr.col_indices() + val = local_csr.values() + + # shift columns + col_off = col + j0 + + nnz = val.numel() + device = local_csr.device + itype = crow.dtype + + # Build global crow: zeros up to i0, then local crow, then constant nnz afterwards + crow_g = torch.empty(H + 1, dtype=itype, device=device) + if i0 > 0: + crow_g[:i0] = 0 + crow_g[i0:i0 + h + 1] = crow + if i0 + h < H: + crow_g[i0 + h + 1:] = nnz + + return torch.sparse_csr_tensor(crow_g, col_off, val, size=global_shape, + device=device, dtype=local_csr.dtype) + + + +def insert_no_overlap_keep_A(A_coo, B_csr): + A = A_coo.coalesce() + B = B_csr.to_sparse_coo() + I_A, V_A = A.indices(), A.values() + I_B, V_B = B.indices(), B.values() + # linearize indices to detect overlaps + lin_A = I_A[0] * A.shape[1] + I_A[1] + lin_B = I_B[0] * A.shape[1] + I_B[1] + # keep only B entries not in A + mask = ~torch.isin(lin_B, lin_A) + I = torch.cat([I_A, I_B[:, mask]], dim=1) + V = torch.cat([V_A, V_B[mask]], dim=0) + return torch.sparse_coo_tensor(I, V, A.shape, device=A.device, dtype=A.dtype) + +def sparse_block(A, i0, i1, j0, j1, out_layout="coo"): + """ + Return A[i0:i1, j0:j1] as a sparse tensor without densifying. + Works for A in COO or CSR. out_layout: 'coo' or 'csr' + """ + Acoo = A.to_sparse_coo().coalesce() + r, c = Acoo.indices() + v = Acoo.values() + + m = (r >= i0) & (r < i1) & (c >= j0) & (c < j1) + r2 = r[m] - i0 + c2 = c[m] - j0 + v2 = v[m] + + B = torch.sparse_coo_tensor( + torch.stack([r2, c2], dim=0), + v2, + (i1 - i0, j1 - j0), + device=A.device, dtype=A.dtype + ).coalesce() + + return B if out_layout == "coo" else B.to_sparse_csr() + + + + +def torch_csr_to_scipy(A: torch.Tensor) -> sp.csr_matrix: + """"Convert matrix/linear forward model, as torch tensor, to scipy sparse matrix. + Arguments: + - A: torch tensor as matrix/forward model. + Returns: + - Scipy sparse matrix. + """ + assert A.layout == torch.sparse_csr, f"A must be sparse_csr {type(A.layout)}" + m, n = A.shape + indptr = A.crow_indices().detach().cpu().numpy() # length m+1 + indices = A.col_indices().detach().cpu().numpy() # length nnz + data = A.values().detach().cpu().numpy() # length nnz + + # (optional) SciPy prefers int32 for indices + if indices.dtype != np.int32: indices = indices.astype(np.int32, copy=False) + if indptr.dtype != np.int32: indptr = indptr.astype(np.int32, copy=False) + + sp_matrix = sp.csr_matrix((data, indices, indptr), shape=(m, n)) + + return sp_matrix + +def torch_coo_to_scipy(A: torch.Tensor, to='csr'): + A = A.to_sparse_coo().coalesce() + m, n = A.shape + ij = A.indices().detach().cpu().numpy() # shape [2, nnz] + data = A.values().detach().cpu().numpy() + row, col = ij[0], ij[1] + M = sp.coo_matrix((data, (row, col)), shape=(m, n)) + return M.tocsr() if to == 'csr' else M + +def torch_sparse_to_scipy_csr(A: torch.Tensor) -> sp.csr_matrix: + """ + Convert a PyTorch tensor (CSR/COO/dense) to SciPy csr_matrix. + - No densifying for sparse inputs. + - Handles CUDA tensors by moving index/data arrays to CPU. + """ + m, n = A.shape + + # Dense (strided) tensor + if A.layout == torch.strided: + return sp.csr_matrix(A.detach().cpu().numpy()) + + # Try CSR fast-path + try: + indptr = A.crow_indices() # only valid for CSR + indices = A.col_indices() + data = A.values() + except (AttributeError, RuntimeError): + # Not CSR → go through COO safely + Acoo = A.to_sparse_coo().coalesce() + ij = Acoo.indices().detach().cpu().numpy() # [2, nnz] + row, col = ij[0].astype(np.int32, copy=False), ij[1].astype(np.int32, copy=False) + dat = Acoo.values().detach().cpu().numpy() + return sp.coo_matrix((dat, (row, col)), shape=(m, n)).tocsr() + else: + # CSR build (ensure int32 indices for SciPy) + indptr_np = indptr.detach().cpu().numpy().astype(np.int32, copy=False) + indices_np = indices.detach().cpu().numpy().astype(np.int32, copy=False) + data_np = data.detach().cpu().numpy() + return sp.csr_matrix((data_np, indices_np, indptr_np), shape=(m, n)) + \ No newline at end of file diff --git a/test/test_toy_example_cuda.py b/test/test_toy_example_torch.py similarity index 79% rename from test/test_toy_example_cuda.py rename to test/test_toy_example_torch.py index 5f03df7..74f6f99 100644 --- a/test/test_toy_example_cuda.py +++ b/test/test_toy_example_torch.py @@ -1,15 +1,22 @@ -import numpy as np -import torch -from torch.utils.data import DataLoader -from playground.generator_functions import random_uni_points_in_ball -from accuracy_bounds.inverseproblems.feasible_sets import compute_feasible_set_linear_forwardmodel -from accuracy_bounds.inverseproblems.kersize_compute import worstcase_kernelsize, worstcase_kernelsize_sym, average_kernelsize, average_kernelsize_sym -from accuracy_bounds.inverseproblems.utils import apply_forwardmodel -from accuracy_bounds.inverseproblems.kersize_compute_dataloader import worstcase_kernelsize_appartenance, average_kernelsize_appartenance, target_distances_cuda_V2 -from accuracy_bounds.inverseproblems.feasible_sets_dataloader import feasible_appartenance_additive_noise_dataloader_cuda, feasible_appartenance_additive_noise_cuda -from accuracy_bounds.inverseproblems.utils import torch_sparse_to_scipy_csr, torch_csr_to_scipy +import importlib.util +import pytest +has_torch = importlib.util.find_spec("torch") is not None +if has_torch: + import numpy as np + import torch + from torch.utils.data import DataLoader + from playground.generator_functions import random_uni_points_in_ball + from accuracy_bounds.inverseproblems.feasible_sets import compute_feasible_set_linear_forwardmodel + from accuracy_bounds.inverseproblems.kersize_compute import worstcase_kernelsize, worstcase_kernelsize_sym, average_kernelsize, average_kernelsize_sym + from accuracy_bounds.inverseproblems.utils import apply_forwardmodel + from accuracy_bounds.inverseproblems.kersize_compute_dataloader import worstcase_kernelsize_appartenance, average_kernelsize_appartenance, target_distances_cuda_V2 + from accuracy_bounds.inverseproblems.feasible_sets_dataloader import feasible_appartenance_additive_noise_dataloader_cuda, feasible_appartenance_additive_noise_cuda + from accuracy_bounds.inverseproblems.utils import torch_sparse_to_scipy_csr, torch_csr_to_scipy + + +@pytest.mark.skipif(not has_torch, reason="requires torch") def test_worstcase_cuda(): num_points = 3000 radius = 2 @@ -54,6 +61,7 @@ def test_worstcase_cuda(): assert error < 0.3, f"Analytic Worstcase Kernel (Cuda) Error: {error}" +@pytest.mark.skipif(not has_torch, reason="requires torch") def test_average_kersize_cuda(): num_points = 3000 radius = 2 @@ -97,6 +105,7 @@ def test_average_kersize_cuda(): assert error < 0.3, f"Analytic Average Kernel (Cuda) Error: {error}" +@pytest.mark.skipif(not has_torch, reason="requires torch") def test_feas_appartance_w_and_wo_dataloader(): num_points = 3000 @@ -129,7 +138,4 @@ def _to_dense(x): error_tensor = (a - b).abs() error = float(error_tensor.mean().item()) # scalar Python float - assert error < 0.3, f"Feasible Appearance Error: {error}" - - -test_feas_appartance_w_and_wo_dataloader() \ No newline at end of file + assert error < 0.3, f"Feasible Appearance Error: {error}" \ No newline at end of file