From fd12bc9e03d902f1e69e36a30ba6d524b618740a Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 11 Jul 2024 17:37:57 +0200 Subject: [PATCH 01/10] Run convolve using Cupy if GPU is available --- invisible_cities/reco/deconv_functions.py | 76 ++++++++++++++++++----- 1 file changed, 61 insertions(+), 15 deletions(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index 33504cad7f..136df9fa07 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -10,6 +10,24 @@ from scipy.signal import convolve from scipy.spatial.distance import cdist +# Check if there is a GPU available +USING_GPU = False + +try: + import cupy as cp + import cupyx as cpx + from cupyx.scipy.signal import convolve + from cupyx.scipy.signal import fftconvolve + from cupyx.scipy.ndimage import convolve + + if cp.cuda.runtime.getDeviceCount(): + USING_GPU = True + print("GPUs available:", cp.cuda.runtime.getDeviceCount()) + else: + raise ImportError +except ImportError: + print("GPUs not available") + from ..core .core_functions import shift_to_bin_centers from ..core .core_functions import in_range @@ -291,25 +309,53 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0.): else: convolve_method = convolve - image = image.astype(float) - psf = psf.astype(float) - im_deconv = 0.5 * np.ones(image.shape) - s = slice(None, None, -1) - psf_mirror = psf[(s,) * psf.ndim] ### Allow for n-dim mirroring. - eps = np.finfo(image.dtype).eps ### Protection against 0 value - ref_image = image/image.max() - - for i in range(iterations): - x = convolve_method(im_deconv, psf, 'same') - np.place(x, x==0, eps) ### Protection against 0 value + if USING_GPU: + xp = cp + else: + xp = np + + # Convert image and psf to the appropriate array type (float) + image = xp.asarray(image, dtype=xp.float32) + psf = xp.asarray(psf, dtype=xp.float32) + + # Initialize im_deconv + im_deconv = 0.5 * xp.ones(image.shape, dtype=xp.float32) + + # Allow for n-dimensional mirroring + s = slice(None, None, -1) + psf_mirror = psf[(s,) * psf.ndim] + + # Protection against 0 value + eps = xp.finfo(image.dtype).eps + + # Normalize the reference image + ref_image = image / image.max() + + for _ in range(iterations): + x = convolve_method(im_deconv, psf, mode='same') + xp.place(x, x == 0, eps) relative_blur = image / x im_deconv *= convolve_method(relative_blur, psf_mirror, 'same') - with np.errstate(divide='ignore', invalid='ignore'): - rel_diff = np.sum(np.divide(((im_deconv/im_deconv.max() - ref_image)**2), ref_image)) - if rel_diff < iter_thr: ### Break if a given threshold is reached. + im_deconv_max = im_deconv.max() + im_deconv_ratio = im_deconv / im_deconv_max + + # This is needed because Cupy does not have a full errstate implementation + if USING_GPU: + rel_diff_array = (im_deconv_ratio - ref_image) ** 2 + with cpx.errstate(): + rel_diff_array = xp.where(ref_image != 0, rel_diff_array / ref_image, 0) + rel_diff = xp.sum(rel_diff_array) + else: + with np.errstate(divide='ignore', invalid='ignore'): + rel_diff = xp.sum(((im_deconv_ratio - ref_image) ** 2) / ref_image) + + if rel_diff < iter_thr: break - ref_image = im_deconv/im_deconv.max() + ref_image = im_deconv_ratio + + if USING_GPU: + im_deconv = cp.asnumpy(im_deconv) return im_deconv From 7b6b5336a0a39d52e2711da5e7cb2e828c6e022d Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 11 Jul 2024 18:06:08 +0200 Subject: [PATCH 02/10] fix types --- invisible_cities/reco/deconv_functions.py | 34 ++++++++--------------- 1 file changed, 11 insertions(+), 23 deletions(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index 136df9fa07..f52904689d 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -18,7 +18,6 @@ import cupyx as cpx from cupyx.scipy.signal import convolve from cupyx.scipy.signal import fftconvolve - from cupyx.scipy.ndimage import convolve if cp.cuda.runtime.getDeviceCount(): USING_GPU = True @@ -315,45 +314,34 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0.): xp = np # Convert image and psf to the appropriate array type (float) - image = xp.asarray(image, dtype=xp.float32) - psf = xp.asarray(psf, dtype=xp.float32) - - # Initialize im_deconv - im_deconv = 0.5 * xp.ones(image.shape, dtype=xp.float32) - - # Allow for n-dimensional mirroring - s = slice(None, None, -1) + image = xp.asarray(image, dtype=float) + psf = xp.asarray(psf, dtype=float) + im_deconv = 0.5 * xp.ones(image.shape) + s = slice(None, None, -1) psf_mirror = psf[(s,) * psf.ndim] + eps = xp.finfo(image.dtype).eps + ref_image = image / image.max() - # Protection against 0 value - eps = xp.finfo(image.dtype).eps - - # Normalize the reference image - ref_image = image / image.max() - for _ in range(iterations): - x = convolve_method(im_deconv, psf, mode='same') + x = convolve_method(im_deconv, psf, 'same') xp.place(x, x == 0, eps) relative_blur = image / x im_deconv *= convolve_method(relative_blur, psf_mirror, 'same') - im_deconv_max = im_deconv.max() - im_deconv_ratio = im_deconv / im_deconv_max - # This is needed because Cupy does not have a full errstate implementation if USING_GPU: - rel_diff_array = (im_deconv_ratio - ref_image) ** 2 + rel_diff_array = (im_deconv/im_deconv.max() - ref_image) ** 2 with cpx.errstate(): - rel_diff_array = xp.where(ref_image != 0, rel_diff_array / ref_image, 0) + rel_diff_array = xp.where(ref_image != 0, xp.divide(rel_diff_array, ref_image), 0) rel_diff = xp.sum(rel_diff_array) else: with np.errstate(divide='ignore', invalid='ignore'): - rel_diff = xp.sum(((im_deconv_ratio - ref_image) ** 2) / ref_image) + rel_diff = xp.sum(xp.divide(((im_deconv/im_deconv.max() - ref_image)**2), ref_image)) if rel_diff < iter_thr: break - ref_image = im_deconv_ratio + ref_image = im_deconv/im_deconv.max() if USING_GPU: im_deconv = cp.asnumpy(im_deconv) From 37c27c72e18c9ecf7c5b3f050f6c81fe9d0f6483 Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 11 Jul 2024 18:07:57 +0200 Subject: [PATCH 03/10] fix comments --- invisible_cities/reco/deconv_functions.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index f52904689d..c909a0b714 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -318,13 +318,13 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0.): psf = xp.asarray(psf, dtype=float) im_deconv = 0.5 * xp.ones(image.shape) s = slice(None, None, -1) - psf_mirror = psf[(s,) * psf.ndim] - eps = xp.finfo(image.dtype).eps - ref_image = image / image.max() + psf_mirror = psf[(s,) * psf.ndim] ### Allow for n-dim mirroring. + eps = xp.finfo(image.dtype).eps ### Protection against 0 value + ref_image = image/image.max() for _ in range(iterations): x = convolve_method(im_deconv, psf, 'same') - xp.place(x, x == 0, eps) + xp.place(x, x==0, eps) ### Protection against 0 value relative_blur = image / x im_deconv *= convolve_method(relative_blur, psf_mirror, 'same') From 62041f72082aa94ac7cf8fd532c1cc6c67f2d8ad Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 11 Jul 2024 18:08:50 +0200 Subject: [PATCH 04/10] fix comments --- invisible_cities/reco/deconv_functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index c909a0b714..b24ee6ffb1 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -338,7 +338,7 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0.): with np.errstate(divide='ignore', invalid='ignore'): rel_diff = xp.sum(xp.divide(((im_deconv/im_deconv.max() - ref_image)**2), ref_image)) - if rel_diff < iter_thr: + if rel_diff < iter_thr: ### Break if a given threshold is reached. break ref_image = im_deconv/im_deconv.max() From adade004a2971dde2c637ff15799c4e9828d6afb Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 12 Jun 2025 12:49:53 +0200 Subject: [PATCH 05/10] cache gpu test and add use_gpu flag for beersheba --- invisible_cities/cities/beersheba.py | 9 ++-- invisible_cities/reco/deconv_functions.py | 60 +++++++++++++++-------- 2 files changed, 45 insertions(+), 24 deletions(-) diff --git a/invisible_cities/cities/beersheba.py b/invisible_cities/cities/beersheba.py index f5c15d629c..fe4c242ea0 100644 --- a/invisible_cities/cities/beersheba.py +++ b/invisible_cities/cities/beersheba.py @@ -133,7 +133,8 @@ def deconvolve_signal(det_db : pd.DataFrame, n_dim : Optional[int]=2, cut_type : Optional[CutType]=CutType.abs, inter_method : Optional[InterpolationMethod]=InterpolationMethod.cubic, - n_iterations_g : Optional[int]=0): + n_iterations_g : Optional[int]=0, + use_gpu : Optional[bool]=False): """ Applies Lucy Richardson deconvolution to SiPM response with a @@ -172,7 +173,7 @@ def deconvolve_signal(det_db : pd.DataFrame, det_grid = [np.arange(det_db[var].min() + bs/2, det_db[var].max() - bs/2 + np.finfo(np.float32).eps, bs) for var, bs in zip(dimensions, bin_size)] deconvolution = deconvolve(n_iterations, iteration_tol, - sample_width, det_grid , inter_method) + sample_width, det_grid, inter_method, use_gpu) if not isinstance(energy_type , HitEnergy ): raise ValueError(f'energy_type {energy_type} is not a valid energy type.') @@ -211,7 +212,7 @@ def deconvolve_hits(df, z): psf_cols = psf.loc[:, cols] gaus = dist.pdf(psf_cols.values) psf = gaus.reshape(psf_cols.nunique()) - deconv_image = nan_to_num(richardson_lucy(deconv_image, psf, n_iterations_g, iteration_tol)) + deconv_image = nan_to_num(richardson_lucy(deconv_image, psf, n_iterations_g, iteration_tol, use_gpu)) return create_deconvolution_df(df, deconv_image.flatten(), pos, cut_type, e_cut, n_dim) @@ -439,6 +440,8 @@ def beersheba( files_in : OneOrManyFiles 'cubic' not supported for 3D deconvolution. n_iterations_g : int Number of Lucy-Richardson iterations for gaussian in 'separate mode' + use_gpu : bool + Whether to use the GPU for the deconvolution. ---------- Input diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index b24ee6ffb1..4dac2ea26e 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -1,5 +1,6 @@ import numpy as np import pandas as pd +import warnings from typing import List from typing import Tuple @@ -13,25 +14,16 @@ # Check if there is a GPU available USING_GPU = False -try: - import cupy as cp - import cupyx as cpx - from cupyx.scipy.signal import convolve - from cupyx.scipy.signal import fftconvolve - if cp.cuda.runtime.getDeviceCount(): - USING_GPU = True - print("GPUs available:", cp.cuda.runtime.getDeviceCount()) - else: - raise ImportError -except ImportError: - print("GPUs not available") from ..core .core_functions import shift_to_bin_centers from ..core .core_functions import in_range from .. types.symbols import InterpolationMethod +from functools import lru_cache + + def cut_and_redistribute_df(cut_condition : str, variables : List[str]=[]) -> Callable: @@ -215,7 +207,8 @@ def deconvolve(n_iterations : int, iteration_tol : float, sample_width : List[float], det_grid : List[np.ndarray], - inter_method : InterpolationMethod = InterpolationMethod.cubic + inter_method : InterpolationMethod = InterpolationMethod.cubic, + use_gpu : bool = False ) -> Callable: """ Deconvolves a given set of data (sensor position and its response) @@ -232,7 +225,8 @@ def deconvolve(n_iterations : int, iteration_tol : Stopping threshold (difference between iterations). sample_width : Sampling size of the sensors. det_grid : xy-coordinates of the detector grid, to interpolate on them - inter_method : Interpolation method. + inter_method : Interpolation method. Default is cubic. + use_gpu : Use GPU for the deconvolution. Default is False. Returns ---------- @@ -251,13 +245,37 @@ def deconvolve(data : Tuple[np.ndarray, ...], columns = var_name[:len(data)] psf_deco = psf.factor.values.reshape(psf.loc[:, columns].nunique().values) deconv_image = np.nan_to_num(richardson_lucy(inter_signal, psf_deco, - n_iterations, iteration_tol)) + n_iterations, iteration_tol, use_gpu)) return deconv_image, inter_pos return deconvolve -def richardson_lucy(image, psf, iterations=50, iter_thr=0.): + +@lru_cache +def is_gpu_available() -> bool: + ''' + Check if GPUs are available, import the necessary libraries and + return True if they are available, False otherwise. + ''' + try: + import cupy as cp + import cupyx as cpx + from cupyx.scipy.signal import convolve + from cupyx.scipy.signal import fftconvolve + + if cp.cuda.runtime.getDeviceCount(): + print("GPUs available:", cp.cuda.runtime.getDeviceCount()) + return True + else: + warnings.warn("use_gpu was set to True but no GPUs are available") + return False + except ImportError: + warnings.warn("use_gpu was set to True but cupy is not installed") + return False + + +def richardson_lucy(image, psf, iterations=50, iter_thr=0., use_gpu=False): """Richardson-Lucy deconvolution (modification from scikit-image package). The modification adds a value=0 protection, the possibility to stop iterating @@ -308,10 +326,10 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0.): else: convolve_method = convolve - if USING_GPU: - xp = cp - else: - xp = np + xp = np + if use_gpu: + if is_gpu_available(): + xp = cp # Convert image and psf to the appropriate array type (float) image = xp.asarray(image, dtype=float) @@ -329,7 +347,7 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0.): im_deconv *= convolve_method(relative_blur, psf_mirror, 'same') # This is needed because Cupy does not have a full errstate implementation - if USING_GPU: + if use_gpu and is_gpu_available(): rel_diff_array = (im_deconv/im_deconv.max() - ref_image) ** 2 with cpx.errstate(): rel_diff_array = xp.where(ref_image != 0, xp.divide(rel_diff_array, ref_image), 0) From 99544d6241a1c02f4ea7c973d5c01b7b6c364616 Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 12 Jun 2025 13:28:01 +0200 Subject: [PATCH 06/10] improve gpu check --- invisible_cities/reco/deconv_functions.py | 39 ++++++++++++----------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index 4dac2ea26e..70adc62a47 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -12,9 +12,14 @@ from scipy.spatial.distance import cdist # Check if there is a GPU available -USING_GPU = False - - +try: + import cupy as cp + import cupyx as cpx + from cupyx.scipy.signal import convolve + from cupyx.scipy.signal import fftconvolve + import cupyx.scipy.ndimage as cpx_ndi +except ImportError: + warnings.warn("Impossible to import cupy. Computations will be done in CPU.") from ..core .core_functions import shift_to_bin_centers from ..core .core_functions import in_range @@ -254,24 +259,20 @@ def deconvolve(data : Tuple[np.ndarray, ...], @lru_cache def is_gpu_available() -> bool: - ''' - Check if GPUs are available, import the necessary libraries and - return True if they are available, False otherwise. - ''' + """Check if a GPU is available for computations. + Returns + ------- + bool + True if a GPU is available, False otherwise. + """ try: - import cupy as cp - import cupyx as cpx - from cupyx.scipy.signal import convolve - from cupyx.scipy.signal import fftconvolve - if cp.cuda.runtime.getDeviceCount(): print("GPUs available:", cp.cuda.runtime.getDeviceCount()) return True else: - warnings.warn("use_gpu was set to True but no GPUs are available") + warnings.warn("Cupy is installed bu no GPUs are available. Computations will be done in CPU.") return False - except ImportError: - warnings.warn("use_gpu was set to True but cupy is not installed") + except NameError: return False @@ -327,9 +328,9 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0., use_gpu=False): convolve_method = convolve xp = np - if use_gpu: - if is_gpu_available(): - xp = cp + using_gpu = use_gpu and is_gpu_available() + if using_gpu: + xp = cp # Convert image and psf to the appropriate array type (float) image = xp.asarray(image, dtype=float) @@ -361,7 +362,7 @@ def richardson_lucy(image, psf, iterations=50, iter_thr=0., use_gpu=False): ref_image = im_deconv/im_deconv.max() - if USING_GPU: + if using_gpu: im_deconv = cp.asnumpy(im_deconv) return im_deconv From 39dad459abd58e5ed7947f7d7789554ce502fe08 Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 12 Jun 2025 16:33:27 +0200 Subject: [PATCH 07/10] fix syntax error --- invisible_cities/cities/beersheba.py | 2 +- invisible_cities/reco/deconv_functions.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/invisible_cities/cities/beersheba.py b/invisible_cities/cities/beersheba.py index 952b95f0bd..47c6645f59 100644 --- a/invisible_cities/cities/beersheba.py +++ b/invisible_cities/cities/beersheba.py @@ -138,7 +138,7 @@ def deconvolve_signal(det_db : pd.DataFrame, cut_type : Optional[CutType]=CutType.abs, inter_method : Optional[InterpolationMethod]=InterpolationMethod.cubic, n_iterations_g : Optional[int]=0, - use_gpu : Optional[bool]=False): + use_gpu : Optional[bool]=False): """ Applies Lucy Richardson deconvolution to SiPM response with a given set of PSFs and parameters. diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index c7bdc6f40c..7938f41ecb 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -327,8 +327,8 @@ def deconvolve(n_iterations : int, satellite_max_size : int, e_cut : float, cut_type : Optional[CutType] = CutType.abs, - inter_method : InterpolationMethod = InterpolationMethod.cubic - use_gpu : bool = False + inter_method : InterpolationMethod = InterpolationMethod.cubic, + use_gpu : Optional[bool] = False ) -> Callable: """ Deconvolves a given set of data (sensor position and its response) From 0ef7b53362342991c76998a8c5e3298dba101c38 Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 12 Jun 2025 16:52:04 +0200 Subject: [PATCH 08/10] fix keyword arg --- invisible_cities/cities/beersheba.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/invisible_cities/cities/beersheba.py b/invisible_cities/cities/beersheba.py index 47c6645f59..493c2b534f 100644 --- a/invisible_cities/cities/beersheba.py +++ b/invisible_cities/cities/beersheba.py @@ -187,7 +187,7 @@ def deconvolve_signal(det_db : pd.DataFrame, deconvolution = deconvolve(n_iterations, iteration_tol, sample_width, det_grid, **satellite_params, - inter_method = inter_method, use_gpu) + inter_method=inter_method, use_gpu=use_gpu) if not isinstance(energy_type , HitEnergy ): raise ValueError(f'energy_type {energy_type} is not a valid energy type.') @@ -230,7 +230,7 @@ def deconvolve_hits(df, z): deconv_image = nan_to_num(richardson_lucy(deconv_image, psf, iterations = n_iterations_g, iter_thr = iteration_tol, - **satellite_params, use_gpu)) + **satellite_params, use_gpu=use_gpu)) return create_deconvolution_df(df, deconv_image.flatten(), pos, cut_type, e_cut, n_dim) From 93ba84e0e3967e40193f5bffea805ba93ec9309c Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Thu, 12 Jun 2025 16:58:05 +0200 Subject: [PATCH 09/10] fix undeclared variable... --- invisible_cities/reco/deconv_functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index 7938f41ecb..9bced484ca 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -473,7 +473,7 @@ def richardson_lucy(image, psf, satellite_start_iter, satellite_max_size, e_cut, eps = xp.finfo(image.dtype).eps ### Protection against 0 value ref_image = image/image.max() - for _ in range(iterations): + for i in range(iterations): x = convolve_method(im_deconv, psf, 'same') xp.place(x, x==0, eps) ### Protection against 0 value relative_blur = image / x From f432109de9e3e912237c914a01f20cfb828984ed Mon Sep 17 00:00:00 2001 From: Stefano Roberto Soleti Date: Mon, 16 Jun 2025 16:43:38 +0200 Subject: [PATCH 10/10] spelling mistake Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- invisible_cities/reco/deconv_functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index 9bced484ca..6038198ae4 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -391,7 +391,7 @@ def is_gpu_available() -> bool: print("GPUs available:", cp.cuda.runtime.getDeviceCount()) return True else: - warnings.warn("Cupy is installed bu no GPUs are available. Computations will be done in CPU.") + warnings.warn("Cupy is installed but no GPUs are available. Computations will be done in CPU.") return False except NameError: return False