From f800f7d4a8cd0fda415b769e4b27d94855f81f96 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Thu, 10 Apr 2025 00:01:53 +0200 Subject: [PATCH 1/4] Alternative solution to #5738 Signed-off-by: Quentin Remy --- gammapy/modeling/models/spectral.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py index 4e786b1cfd..033df9290e 100644 --- a/gammapy/modeling/models/spectral.py +++ b/gammapy/modeling/models/spectral.py @@ -491,6 +491,7 @@ def plot_error( sed_type="dnde", energy_power=0, n_points=100, + method="lin", **kwargs, ): """Plot spectral model error band. @@ -560,8 +561,16 @@ def plot_error( ax.yaxis.set_units(DEFAULT_UNIT[sed_type] * energy.unit**energy_power) flux, flux_err = self._get_plot_flux(sed_type=sed_type, energy=energy) - y_lo = scale_plot_flux(flux - flux_err, energy_power).quantity[:, 0, 0] - y_hi = scale_plot_flux(flux + flux_err, energy_power).quantity[:, 0, 0] + + if method == "log": + flux_err_log = flux_err.data / flux.data + y_lo = flux / np.exp(flux_err_log) + y_hi = flux * np.exp(flux_err_log) + y_lo = scale_plot_flux(y_lo, energy_power).quantity[:, 0, 0] + y_hi = scale_plot_flux(y_hi, energy_power).quantity[:, 0, 0] + else: + y_lo = scale_plot_flux(flux - flux_err, energy_power).quantity[:, 0, 0] + y_hi = scale_plot_flux(flux + flux_err, energy_power).quantity[:, 0, 0] with quantity_support(): ax.fill_between(energy.center, y_lo, y_hi, **kwargs) From c886800c55027ced23afa88d5bc97a2781c3918a Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Thu, 10 Apr 2025 02:12:42 +0200 Subject: [PATCH 2/4] implement the samples method Signed-off-by: Quentin Remy --- gammapy/modeling/models/spectral.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py index 033df9290e..6c853dd84d 100644 --- a/gammapy/modeling/models/spectral.py +++ b/gammapy/modeling/models/spectral.py @@ -491,6 +491,7 @@ def plot_error( sed_type="dnde", energy_power=0, n_points=100, + n_samples=10000, method="lin", **kwargs, ): @@ -568,10 +569,24 @@ def plot_error( y_hi = flux * np.exp(flux_err_log) y_lo = scale_plot_flux(y_lo, energy_power).quantity[:, 0, 0] y_hi = scale_plot_flux(y_hi, energy_power).quantity[:, 0, 0] - else: + elif method == "lin": y_lo = scale_plot_flux(flux - flux_err, energy_power).quantity[:, 0, 0] y_hi = scale_plot_flux(flux + flux_err, energy_power).quantity[:, 0, 0] - + elif method == "samples": + samples = np.random.multivariate_normal( + self.parameters.value, self.covariance, n_samples + ) + f_samples = np.array( + [ + self.evaluate(energy.center.value, *samples[k, :]) + for k in range(n_samples) + ] + ) + f_unit = self(1.0 * energy.center.unit).unit + y_lo = flux.copy(data=np.percentile(f_samples, 16, axis=0), unit=f_unit) + y_hi = flux.copy(data=np.percentile(f_samples, 84, axis=0), unit=f_unit) + y_lo = scale_plot_flux(y_lo, energy_power).quantity[:, 0, 0] + y_hi = scale_plot_flux(y_hi, energy_power).quantity[:, 0, 0] with quantity_support(): ax.fill_between(energy.center, y_lo, y_hi, **kwargs) From fc1a7334554645fd92b6aad6de24c9f3b1d1f3a9 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Thu, 10 Apr 2025 11:01:20 +0200 Subject: [PATCH 3/4] support for sed_type Signed-off-by: Quentin Remy --- gammapy/modeling/models/spectral.py | 39 +++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 8 deletions(-) diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py index 6c853dd84d..7673d54003 100644 --- a/gammapy/modeling/models/spectral.py +++ b/gammapy/modeling/models/spectral.py @@ -538,6 +538,7 @@ def plot_error( ``n_points`` bins between the given bounds. """ from gammapy.estimators.map.core import DEFAULT_UNIT + from gammapy.estimators import FluxMaps if self.is_norm_spectral_model: sed_type = "norm" @@ -576,17 +577,39 @@ def plot_error( samples = np.random.multivariate_normal( self.parameters.value, self.covariance, n_samples ) - f_samples = np.array( - [ - self.evaluate(energy.center.value, *samples[k, :]) - for k in range(n_samples) - ] + f = self(energy.center) + spl_ax = MapAxis(range(n_samples), node_type="center", name="sample") + f_samples = RegionNDMap.create(region=None, axes=[energy, spl_ax]) + f_samples.quantity = ( + np.array( + [ + self.evaluate(energy.center.value, *samples[k, :]) + for k in range(n_samples) + ] + ) + * f.unit + ) + + f_samples_map = FluxMaps.from_maps( + dict(dnde=f_samples), sed_type="dnde", reference_model=self + ) + f_samples = f_samples_map[sed_type] + + y_lo = RegionNDMap.create( + region=None, + axes=[energy], + unit=f_samples.unit, + data=np.percentile(f_samples, 16, axis=0), + ) + y_hi = RegionNDMap.create( + region=None, + axes=[energy], + unit=f_samples.unit, + data=np.percentile(f_samples, 84, axis=0), ) - f_unit = self(1.0 * energy.center.unit).unit - y_lo = flux.copy(data=np.percentile(f_samples, 16, axis=0), unit=f_unit) - y_hi = flux.copy(data=np.percentile(f_samples, 84, axis=0), unit=f_unit) y_lo = scale_plot_flux(y_lo, energy_power).quantity[:, 0, 0] y_hi = scale_plot_flux(y_hi, energy_power).quantity[:, 0, 0] + with quantity_support(): ax.fill_between(energy.center, y_lo, y_hi, **kwargs) From 7d091cda432f9c1da7c7f1c7241649fc3766813d Mon Sep 17 00:00:00 2001 From: Mireia Date: Thu, 10 Apr 2025 16:51:25 +0100 Subject: [PATCH 4/4] Fixes to PR --- gammapy/modeling/models/spectral.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py index 7673d54003..7bc7512358 100644 --- a/gammapy/modeling/models/spectral.py +++ b/gammapy/modeling/models/spectral.py @@ -577,13 +577,17 @@ def plot_error( samples = np.random.multivariate_normal( self.parameters.value, self.covariance, n_samples ) + # Reapply units to the sampled parameters + parameter_units = [par.unit for par in self.parameters] + samples = np.array([sample * parameter_units for sample in samples]) + f = self(energy.center) spl_ax = MapAxis(range(n_samples), node_type="center", name="sample") f_samples = RegionNDMap.create(region=None, axes=[energy, spl_ax]) f_samples.quantity = ( np.array( [ - self.evaluate(energy.center.value, *samples[k, :]) + self.evaluate(energy.center, *samples[k, :]) for k in range(n_samples) ] )