From 58d325ac3ce601e6d5406ecd924c923dc231b4f9 Mon Sep 17 00:00:00 2001 From: Gonzalo Martinez Lema Date: Wed, 23 Oct 2024 16:41:04 +0200 Subject: [PATCH 1/7] Test all interpolation methods in `deconv_functions` --- .../reco/deconv_functions_test.py | 37 ++++++++++++++++--- 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/invisible_cities/reco/deconv_functions_test.py b/invisible_cities/reco/deconv_functions_test.py index 8841b7fea..73996674d 100644 --- a/invisible_cities/reco/deconv_functions_test.py +++ b/invisible_cities/reco/deconv_functions_test.py @@ -6,6 +6,7 @@ from pytest import mark from pytest import raises +from pytest import fixture from hypothesis import given from hypothesis.strategies import floats @@ -140,7 +141,17 @@ def test_interpolate_signal(): assert np.allclose(grid , sorted(set(inter_position[1]))) -def test_deconvolution_input(data_hdst, data_hdst_deconvolved): +@fixture(scope="session") +def new_grid_1mm(): + bin_sizes = (1., 1.) + det_db = DataSiPM('new', 0) + 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("X Y".split(), bin_sizes)] + return det_grid + + +def test_deconvolution_input(data_hdst, data_hdst_deconvolved, new_grid_1mm): ref_interpolation = np.load(data_hdst_deconvolved) hdst = load_dst(data_hdst, 'RECO', 'Events') @@ -148,10 +159,7 @@ def test_deconvolution_input(data_hdst, data_hdst_deconvolved): h = h.groupby(['X', 'Y']).Q.sum().reset_index() h = h[h.Q > 40] - det_db = DataSiPM('new', 0) - 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(['X', 'Y'], [1., 1.])] - interpolator = deconvolution_input([10., 10.], det_grid, InterpolationMethod.cubic) + interpolator = deconvolution_input([10., 10.], new_grid_1mm, InterpolationMethod.cubic) inter = interpolator((h.X, h.Y), h.Q) assert np.allclose(ref_interpolation['e_inter'], inter[0]) @@ -159,8 +167,25 @@ def test_deconvolution_input(data_hdst, data_hdst_deconvolved): assert np.allclose(ref_interpolation['y_inter'], inter[1][1]) +@mark.parametrize("interp_method", InterpolationMethod) +def test_deconvolution_input_interpolation_method(data_hdst, new_grid_1mm, interp_method): + """ + Check that it runs with all interpolation methods + """ + hdst = load_dst(data_hdst, 'RECO', 'Events') + hdst = hdst[(hdst.event == 3021916) & (hdst.npeak == 0) & (hdst.Q>40)] + + deconvolve = deconvolution_input([10., 10.], new_grid_1mm, interp_method) + output = deconvolve((hdst.X, hdst.Y), hdst.Q) + + assert output[0].shape == (70, 100) + assert len(output[1]) == 2 + assert output[1][0].shape == (7000,) + assert output[1][1].shape == (7000,) + + @mark.parametrize("interp_method", InterpolationMethod.__members__) -def test_deconvolution_input_interpolation_method(data_hdst, data_hdst_deconvolved, interp_method): +def test_deconvolution_input_wrong_interpolation_method_raises(data_hdst, data_hdst_deconvolved, interp_method): with raises(ValueError): deconvolution_input([10., 10.], [1., 1.], interp_method) From d10a611a37d11c03c550d7965b52c23b3a37953c Mon Sep 17 00:00:00 2001 From: Gonzalo Martinez Lema Date: Wed, 23 Oct 2024 11:40:36 +0200 Subject: [PATCH 2/7] Refactor and fix `deconvolution_input` for the case of no interpolation The no-interpolation case did not generate images with the same size as the other interpolation methods. --- invisible_cities/reco/deconv_functions.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/invisible_cities/reco/deconv_functions.py b/invisible_cities/reco/deconv_functions.py index dd0bcd3cf..799d8e08d 100644 --- a/invisible_cities/reco/deconv_functions.py +++ b/invisible_cities/reco/deconv_functions.py @@ -14,6 +14,7 @@ from scipy import ndimage as ndi from ..core .core_functions import shift_to_bin_centers +from .. core.core_functions import binedges_from_bincenters from ..core .core_functions import in_range from ..core .configure import check_annotations @@ -210,15 +211,17 @@ def deconvolution_input(data : Tuple[np.ndarray, ...], weight : np.ndarray ) -> Tuple[np.ndarray, Tuple[np.ndarray, ...]]: - ranges = [[coord.min() - 1.5 * sw, coord.max() + 1.5 * sw] for coord, sw in zip(data, sample_width)] - if inter_method in (InterpolationMethod.linear, InterpolationMethod.cubic, InterpolationMethod.nearest): - allbins = [np.arange(rang[0], rang[1] + np.finfo(np.float32).eps, sw) for rang, sw in zip(ranges, sample_width)] - Hs, edges = np.histogramdd(data, bins=allbins, normed=False, weights=weight) - elif inter_method is InterpolationMethod.nointerpolation: + eps = np.finfo(np.float32).eps + ranges = [ [ coord.min() - 1.5 * sw + , coord.max() + 1.5 * sw + eps] + for coord, sw in zip(data, sample_width) ] + if inter_method is InterpolationMethod.nointerpolation: allbins = [grid[in_range(grid, *rang)] for rang, grid in zip(ranges, det_grid)] + allbins = [binedges_from_bincenters(bins) for bins in allbins] Hs, edges = np.histogramdd(data, bins=allbins, normed=False, weights=weight) else: - raise ValueError(f'inter_method {inter_method} is not a valid interpolatin mode.') + allbins = [np.arange(*rang, sw) for rang, sw in zip(ranges, sample_width)] + Hs, edges = np.histogramdd(data, bins=allbins, normed=False, weights=weight) inter_points = np.meshgrid(*(shift_to_bin_centers(edge) for edge in edges), indexing='ij') inter_points = tuple (inter_p.flatten() for inter_p in inter_points) From 422783f70753f71e01d824ac56d936006f81f815 Mon Sep 17 00:00:00 2001 From: Gonzalo Martinez Lema Date: Wed, 23 Oct 2024 18:32:34 +0200 Subject: [PATCH 3/7] Remove deprecated test --- invisible_cities/reco/deconv_functions_test.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/invisible_cities/reco/deconv_functions_test.py b/invisible_cities/reco/deconv_functions_test.py index 73996674d..63d29fb12 100644 --- a/invisible_cities/reco/deconv_functions_test.py +++ b/invisible_cities/reco/deconv_functions_test.py @@ -184,12 +184,6 @@ def test_deconvolution_input_interpolation_method(data_hdst, new_grid_1mm, inter assert output[1][1].shape == (7000,) -@mark.parametrize("interp_method", InterpolationMethod.__members__) -def test_deconvolution_input_wrong_interpolation_method_raises(data_hdst, data_hdst_deconvolved, interp_method): - with raises(ValueError): - deconvolution_input([10., 10.], [1., 1.], interp_method) - - def test_deconvolve(data_hdst, data_hdst_deconvolved): ref_interpolation = np.load (data_hdst_deconvolved) From d00d4589f3ba736f1ede670cce69877534c5a6a4 Mon Sep 17 00:00:00 2001 From: Gonzalo Martinez Lema Date: Sat, 2 Nov 2024 12:44:58 +0100 Subject: [PATCH 4/7] Explain tests in docstrings --- invisible_cities/reco/deconv_functions_test.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/invisible_cities/reco/deconv_functions_test.py b/invisible_cities/reco/deconv_functions_test.py index 63d29fb12..c4f67997b 100644 --- a/invisible_cities/reco/deconv_functions_test.py +++ b/invisible_cities/reco/deconv_functions_test.py @@ -151,7 +151,11 @@ def new_grid_1mm(): return det_grid -def test_deconvolution_input(data_hdst, data_hdst_deconvolved, new_grid_1mm): +def test_deconvolution_input_exact_result(data_hdst, data_hdst_deconvolved, new_grid_1mm): + """ + Compare the output of the deconvolution with a reference output + stored in a file. + """ ref_interpolation = np.load(data_hdst_deconvolved) hdst = load_dst(data_hdst, 'RECO', 'Events') @@ -170,7 +174,9 @@ def test_deconvolution_input(data_hdst, data_hdst_deconvolved, new_grid_1mm): @mark.parametrize("interp_method", InterpolationMethod) def test_deconvolution_input_interpolation_method(data_hdst, new_grid_1mm, interp_method): """ - Check that it runs with all interpolation methods + Check that it runs with all interpolation methods. + Select one event/peak from input and apply a cut to obtain an + input image that results in round numbers. """ hdst = load_dst(data_hdst, 'RECO', 'Events') hdst = hdst[(hdst.event == 3021916) & (hdst.npeak == 0) & (hdst.Q>40)] From b01007b62b6e35234249a657c6d116d8881196cb Mon Sep 17 00:00:00 2001 From: Gonzalo Martinez Lema Date: Sat, 2 Nov 2024 12:46:10 +0100 Subject: [PATCH 5/7] Refactor `test_deconvolution_input_interpolation_method` --- .../reco/deconv_functions_test.py | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/invisible_cities/reco/deconv_functions_test.py b/invisible_cities/reco/deconv_functions_test.py index c4f67997b..8b5d28c64 100644 --- a/invisible_cities/reco/deconv_functions_test.py +++ b/invisible_cities/reco/deconv_functions_test.py @@ -179,15 +179,18 @@ def test_deconvolution_input_interpolation_method(data_hdst, new_grid_1mm, inter input image that results in round numbers. """ hdst = load_dst(data_hdst, 'RECO', 'Events') - hdst = hdst[(hdst.event == 3021916) & (hdst.npeak == 0) & (hdst.Q>40)] - - deconvolve = deconvolution_input([10., 10.], new_grid_1mm, interp_method) - output = deconvolve((hdst.X, hdst.Y), hdst.Q) - - assert output[0].shape == (70, 100) - assert len(output[1]) == 2 - assert output[1][0].shape == (7000,) - assert output[1][1].shape == (7000,) + hdst = hdst[(hdst.event == 3021916) & (hdst.npeak == 0)] + hdst = hdst.groupby(list("XY")).Q.sum().reset_index().loc[lambda df: df.Q>40] + + interpolator = deconvolution_input([10., 10.], new_grid_1mm, interp_method) + output = interpolator((hdst.X, hdst.Y), hdst.Q) + + shape = 120, 120 + nelements = shape[0] * shape[1] + assert output[0].shape == shape + assert len(output[1]) == 2 + assert output[1][0].shape == (nelements,) + assert output[1][1].shape == (nelements,) def test_deconvolve(data_hdst, data_hdst_deconvolved): From f41676ce1edf5f59d4b882e7a43650a92b59ff3d Mon Sep 17 00:00:00 2001 From: Gonzalo Martinez Lema Date: Sat, 2 Nov 2024 12:58:57 +0100 Subject: [PATCH 6/7] Extract data preparation into a fixture --- .../reco/deconv_functions_test.py | 31 ++++++++++++------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/invisible_cities/reco/deconv_functions_test.py b/invisible_cities/reco/deconv_functions_test.py index 8b5d28c64..ec41fed81 100644 --- a/invisible_cities/reco/deconv_functions_test.py +++ b/invisible_cities/reco/deconv_functions_test.py @@ -151,20 +151,29 @@ def new_grid_1mm(): return det_grid -def test_deconvolution_input_exact_result(data_hdst, data_hdst_deconvolved, new_grid_1mm): +@fixture(scope="session") +def data_hdst_first_peak(data_hdst): + event = 3021916 + peak = 0 + hdst = ( load_dst(data_hdst, 'RECO', 'Events') + .groupby("event npeak".split()) + .get_group((event, peak)) + .groupby(['X', 'Y']).Q.sum().reset_index() + .loc[lambda df: df.Q > 40] + ) + return hdst + + +def test_deconvolution_input_exact_result(data_hdst_first_peak, data_hdst_deconvolved, new_grid_1mm): """ Compare the output of the deconvolution with a reference output stored in a file. """ + hdst = data_hdst_first_peak ref_interpolation = np.load(data_hdst_deconvolved) - hdst = load_dst(data_hdst, 'RECO', 'Events') - h = hdst[(hdst.event == 3021916) & (hdst.npeak == 0)] - h = h.groupby(['X', 'Y']).Q.sum().reset_index() - h = h[h.Q > 40] - interpolator = deconvolution_input([10., 10.], new_grid_1mm, InterpolationMethod.cubic) - inter = interpolator((h.X, h.Y), h.Q) + inter = interpolator((hdst.X, hdst.Y), hdst.Q) assert np.allclose(ref_interpolation['e_inter'], inter[0]) assert np.allclose(ref_interpolation['x_inter'], inter[1][0]) @@ -172,16 +181,13 @@ def test_deconvolution_input_exact_result(data_hdst, data_hdst_deconvolved, new_ @mark.parametrize("interp_method", InterpolationMethod) -def test_deconvolution_input_interpolation_method(data_hdst, new_grid_1mm, interp_method): +def test_deconvolution_input_interpolation_method(data_hdst_first_peak, new_grid_1mm, interp_method): """ Check that it runs with all interpolation methods. Select one event/peak from input and apply a cut to obtain an input image that results in round numbers. """ - hdst = load_dst(data_hdst, 'RECO', 'Events') - hdst = hdst[(hdst.event == 3021916) & (hdst.npeak == 0)] - hdst = hdst.groupby(list("XY")).Q.sum().reset_index().loc[lambda df: df.Q>40] - + hdst = data_hdst_first_peak interpolator = deconvolution_input([10., 10.], new_grid_1mm, interp_method) output = interpolator((hdst.X, hdst.Y), hdst.Q) @@ -193,6 +199,7 @@ def test_deconvolution_input_interpolation_method(data_hdst, new_grid_1mm, inter assert output[1][1].shape == (nelements,) +#TODO: use data_hdst_first_peak and new_grid_1mm here too, if possible def test_deconvolve(data_hdst, data_hdst_deconvolved): ref_interpolation = np.load (data_hdst_deconvolved) From 850f8b592e28866ee1168c92d0def176c21177ec Mon Sep 17 00:00:00 2001 From: Gonzalo Martinez Lema Date: Sat, 2 Nov 2024 13:18:24 +0100 Subject: [PATCH 7/7] Recover `test_deconvolution_input_wrong_interpolation_method_raises` --- invisible_cities/reco/deconv_functions_test.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/invisible_cities/reco/deconv_functions_test.py b/invisible_cities/reco/deconv_functions_test.py index ec41fed81..b9c9b584f 100644 --- a/invisible_cities/reco/deconv_functions_test.py +++ b/invisible_cities/reco/deconv_functions_test.py @@ -199,6 +199,12 @@ def test_deconvolution_input_interpolation_method(data_hdst_first_peak, new_grid assert output[1][1].shape == (nelements,) +@mark.parametrize("interp_method", InterpolationMethod.__members__) +def test_deconvolution_input_wrong_interpolation_method_raises(interp_method): + with raises(ValueError): + deconvolution_input([10., 10.], [1., 1.], interp_method) + + #TODO: use data_hdst_first_peak and new_grid_1mm here too, if possible def test_deconvolve(data_hdst, data_hdst_deconvolved): ref_interpolation = np.load (data_hdst_deconvolved)