From c3f42837d5377ad75b0db808bf8214c5ebef9528 Mon Sep 17 00:00:00 2001 From: Hiroki Yoneda Date: Tue, 29 Oct 2024 17:42:31 +0100 Subject: [PATCH] Added unittest for MAP RL --- tests/image_deconvolution/conftest.py | 2 +- tests/image_deconvolution/test_algorithm.py | 158 ++++++++++++++++++-- 2 files changed, 150 insertions(+), 10 deletions(-) diff --git a/tests/image_deconvolution/conftest.py b/tests/image_deconvolution/conftest.py index 54b09384..0ddecc52 100644 --- a/tests/image_deconvolution/conftest.py +++ b/tests/image_deconvolution/conftest.py @@ -10,7 +10,7 @@ def dataset(): event_binned_data = Histogram.open(test_data.path / "test_event_histogram_galacticCDS.hdf5").project(["Em", "Phi", "PsiChi"]) - dict_bkg_binned_data = {"bkg": Histogram.open(test_data.path / "test_event_histogram_galacticCDS.hdf5").project(["Em", "Phi", "PsiChi"])} + dict_bkg_binned_data = {"bkg": Histogram.open(test_data.path / "test_event_histogram_galacticCDS.hdf5").project(["Em", "Phi", "PsiChi"]) * 0.1} precomputed_response = Histogram.open(test_data.path / "test_precomputed_response.h5") data = DataIF_COSI_DC2.load(name = "testdata_galacticCDS", diff --git a/tests/image_deconvolution/test_algorithm.py b/tests/image_deconvolution/test_algorithm.py index ea77a23a..79fc44dc 100644 --- a/tests/image_deconvolution/test_algorithm.py +++ b/tests/image_deconvolution/test_algorithm.py @@ -1,11 +1,15 @@ import pytest +import numpy as np from yayc import Configurator -from cosipy.image_deconvolution import RichardsonLucySimple, RichardsonLucy +from cosipy.image_deconvolution import RichardsonLucySimple, RichardsonLucy, MAP_RichardsonLucy def test_RicharsonLucySimple(dataset, model, mask, tmp_path): - parameter = Configurator({"iteration_max": 2, + num_iteration = 3 + + parameter = Configurator({"iteration_max": num_iteration, + "response_weighting": {"activate": True, "index": 0.5}, "minimum_flux": {"value": 0.0, "unit": "cm-2 s-1 sr-1"}, "background_normalization_optimization": {"activate": True, "range": {"bkg": [0.9, 1.1]}}, @@ -19,14 +23,18 @@ def test_RicharsonLucySimple(dataset, model, mask, tmp_path): algorithm.initialization() - algorithm.iteration() - algorithm.iteration() + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break algorithm.finalization() def test_RicharsonLucy(dataset, model, mask, tmp_path): - parameter = Configurator({"iteration_max": 2, + num_iteration = 3 + + parameter = Configurator({"iteration_max": num_iteration, "minimum_flux": {"value": 0.0, "unit": "cm-2 s-1 sr-1"}, "acceleration": {"activate": True, "alpha_max": 10.0}, "response_weighting": {"activate": True, "index": 0.5}, @@ -44,10 +52,14 @@ def test_RicharsonLucy(dataset, model, mask, tmp_path): algorithm.initialization() - algorithm.iteration() - algorithm.iteration() + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break algorithm.finalization() + + assert np.isclose(algorithm.results[-1]['log-likelihood'][0], 5259.242894137228) # wo/ acceleration and overwrite the directory parameter["acceleration:activate"] = False @@ -59,7 +71,135 @@ def test_RicharsonLucy(dataset, model, mask, tmp_path): algorithm.initialization() - algorithm.iteration() - algorithm.iteration() + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break + + algorithm.finalization() + + assert np.isclose(algorithm.results[-1]['log-likelihood'][0], 5226.135408675148) + +def test_MAP_RichardsonLucy(dataset, model, mask, tmp_path): + + num_iteration = 10 + + parameter = Configurator({"iteration_max": num_iteration, + "minimum_flux": {"value": 0.0, "unit": "cm-2 s-1 sr-1"}, + "response_weighting": {"activate": True, "index": 0.5}, + "background_normalization_optimization": {"activate": True, + "range": {"bkg": [0.9, 1.1]}}, + "stopping_criteria": {"statistics": "log-posterior", + "threshold": 1e-2}, + "prior": {"TSV" :{"coefficient": 1.e-10}, + "gamma":{"model":{"theta": {"value": np.inf, "unit": "cm-2 s-1 sr-1"}, + "k": {"value": 0.999}}, + "background": {"theta": {"value": np.inf}, "k": {"value": 1.0}} + } + }, + "save_results": {"activate": True, "directory": f"{str(tmp_path)}", "only_final_result": True} + }) + + # first run + algorithm = MAP_RichardsonLucy(initial_model = model, + dataset = dataset, + mask = mask, + parameter = parameter) + + algorithm.initialization() + + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break + + algorithm.finalization() + + assert np.isclose(algorithm.results[-1]['log-posterior'], 5298.642098311226) + + # background fixed + parameter["background_normalization_optimization:activate"] = False + + algorithm = MAP_RichardsonLucy(initial_model = model, + dataset = dataset, + mask = mask, + parameter = parameter) + + algorithm.initialization() + + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break + + algorithm.finalization() + + assert np.isclose(algorithm.results[-1]['log-posterior'], 4951.657349512739) + + # with large threshold + parameter["stopping_criteria:threshold"] = 1e10 + + algorithm = MAP_RichardsonLucy(initial_model = model, + dataset = dataset, + mask = mask, + parameter = parameter) + + algorithm.initialization() + + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break + + algorithm.finalization() + + assert len(algorithm.results) == 2 + + # with log-likelihood + parameter["stopping_criteria:statistics"] = "log-likelihood" + + algorithm = MAP_RichardsonLucy(initial_model = model, + dataset = dataset, + mask = mask, + parameter = parameter) + + algorithm.initialization() + + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break + + algorithm.finalization() + + assert np.isclose(algorithm.results[-1]['log-likelihood'][0], 4850.641606406938) + + # without gamma prior + parameter["stopping_criteria:statistics"] = 'log-posterior' + parameter["prior"] = {"TSV":{"coefficient": 1.e-10}} + + algorithm = MAP_RichardsonLucy(initial_model = model, + dataset = dataset, + mask = mask, + parameter = parameter) + + algorithm.initialization() + + for i in range(num_iteration): + stop = algorithm.iteration() + if stop: + break algorithm.finalization() + + assert np.isclose(algorithm.results[-1]['log-posterior'], 4850.604374561122) + + # wrong statistics + parameter["stopping_criteria:statistics"] = "likelihooooooooood!!!" + + with pytest.raises(ValueError) as e_info: + algorithm = MAP_RichardsonLucy(initial_model = model, + dataset = dataset, + mask = mask, + parameter = parameter) +