From 7bbb84a69c9793fa0bda2987212781a6b904cabc Mon Sep 17 00:00:00 2001 From: Tobias Preis Date: Thu, 18 Jun 2026 18:18:37 +0200 Subject: [PATCH 1/7] add: framework for display_model --- gammapy_plugin/gammapy_like.py | 71 +++++++++++++++++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/gammapy_plugin/gammapy_like.py b/gammapy_plugin/gammapy_like.py index c0aa438..cf99892 100644 --- a/gammapy_plugin/gammapy_like.py +++ b/gammapy_plugin/gammapy_like.py @@ -1,5 +1,5 @@ import logging -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Union, Optional, Dict, Any import numpy as np from astromodels.core.model import Model @@ -16,6 +16,7 @@ if TYPE_CHECKING: from threeML.analysis_results import _AnalysisResults + from threeML.io.plotting.data_residual_plot import ResidualPlot __all__ = ["GammapyLike"] @@ -271,3 +272,71 @@ def gammapy_model(self): def frame(self) -> str: """Coordinate Frame of the plugin.""" return self._frame + + def display_model( + self, + data_color: str = "k", + model_color: str = "r", + background_color: str = "b", + step: bool = True, + show_data: bool = True, + show_residuals: bool = True, + ratio_residuals: bool = False, + show_legend: bool = True, + min_rate: Union[int, float] = 1e-99, + model_label: Optional[str] = None, + model_kwargs: Optional[Dict[str, Any]] = None, + data_kwargs: Optional[Dict[str, Any]] = None, + background_label: Optional[str] = None, + background_kwargs: Optional[Dict[str, Any]] = None, + source_only: bool = True, + show_background: bool = False, + **kwargs, + ) -> "ResidualPlot": + + from threeML.io.plotting.data_residual_plot import ResidualPlot + + # TODO: have to ensure this is only for specturm data set valid + + residual_plot = ResidualPlot( + show_residuals=True, + ratio_residuals=True, + **kwargs, + ) + + # compute the values for the plotting + # TODO: this is only the case for 1D datasets + + y_unweighted = self._datasets[0].counts.data.reshape(-1) + x = self._datasets[0].counts.geom.axes["energy"].as_plot_center.to("keV").value + xerr = [ + self._datasets[0].counts.geom.axes["energy"].as_plot_xerr[i].to("keV").value + for i in [0, 1] + ] + + bins = ( + self._datasets[0].counts.geom.axes["energy"].as_plot_edges.to("keV").value + ) + widths = np.diff(bins) + y = y_unweighted / widths + + residuals = ( + self._datasets[0].counts.get_spectrum() - self._datasets[0].npred() + ) / self._datasets[0].npred() + + residual_plot.add_data( + x, + y, + residuals, + xerr=xerr, + label=self._name, + show_data=show_data, + ) + + return residual_plot.finalize( + xlabel="Energy\n(keV)", + ylabel="Counts/keV", + xscale="log", + yscale="log", + show_legend=show_legend, + ) From 6e117fc96a7cb0a1c8ad9ecbac79d99594b9caf0 Mon Sep 17 00:00:00 2001 From: Tobias Preis Date: Thu, 18 Jun 2026 18:18:49 +0200 Subject: [PATCH 2/7] refactor: imports --- gammapy_plugin/gammapy_like.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy_plugin/gammapy_like.py b/gammapy_plugin/gammapy_like.py index cf99892..41fe6f2 100644 --- a/gammapy_plugin/gammapy_like.py +++ b/gammapy_plugin/gammapy_like.py @@ -1,5 +1,5 @@ import logging -from typing import TYPE_CHECKING, Union, Optional, Dict, Any +from typing import TYPE_CHECKING, Any, Dict, Optional, Union import numpy as np from astromodels.core.model import Model From 2a78d1dee70940f3b868af746d48c675428be018 Mon Sep 17 00:00:00 2001 From: Tobias Preis Date: Thu, 18 Jun 2026 19:06:53 +0200 Subject: [PATCH 3/7] add: first working display_model --- gammapy_plugin/gammapy_like.py | 74 +++++++++++++--------------------- 1 file changed, 28 insertions(+), 46 deletions(-) diff --git a/gammapy_plugin/gammapy_like.py b/gammapy_plugin/gammapy_like.py index 41fe6f2..4bc7de4 100644 --- a/gammapy_plugin/gammapy_like.py +++ b/gammapy_plugin/gammapy_like.py @@ -275,22 +275,6 @@ def frame(self) -> str: def display_model( self, - data_color: str = "k", - model_color: str = "r", - background_color: str = "b", - step: bool = True, - show_data: bool = True, - show_residuals: bool = True, - ratio_residuals: bool = False, - show_legend: bool = True, - min_rate: Union[int, float] = 1e-99, - model_label: Optional[str] = None, - model_kwargs: Optional[Dict[str, Any]] = None, - data_kwargs: Optional[Dict[str, Any]] = None, - background_label: Optional[str] = None, - background_kwargs: Optional[Dict[str, Any]] = None, - source_only: bool = True, - show_background: bool = False, **kwargs, ) -> "ResidualPlot": @@ -303,40 +287,38 @@ def display_model( ratio_residuals=True, **kwargs, ) - - # compute the values for the plotting - # TODO: this is only the case for 1D datasets - - y_unweighted = self._datasets[0].counts.data.reshape(-1) - x = self._datasets[0].counts.geom.axes["energy"].as_plot_center.to("keV").value - xerr = [ - self._datasets[0].counts.geom.axes["energy"].as_plot_xerr[i].to("keV").value - for i in [0, 1] - ] - - bins = ( - self._datasets[0].counts.geom.axes["energy"].as_plot_edges.to("keV").value - ) - widths = np.diff(bins) - y = y_unweighted / widths - - residuals = ( - self._datasets[0].counts.get_spectrum() - self._datasets[0].npred() - ) / self._datasets[0].npred() - - residual_plot.add_data( - x, - y, - residuals, - xerr=xerr, - label=self._name, - show_data=show_data, - ) + for i in range(len(self._datasets)): + y_unweighted = self._datasets[i].counts.data.reshape(-1) + x = self._datasets[i].counts.geom.axes["energy"].as_plot_center.to("keV").value + xerr = [ + self._datasets[i].counts.geom.axes["energy"].as_plot_xerr[j].to("keV").value + for j in [0, 1] + ] + + bins = ( + self._datasets[i].counts.geom.axes["energy"].as_plot_edges.to("keV").value + ) + widths = np.diff(bins) + y = y_unweighted / widths + + residuals = ( + self._datasets[i].counts.get_spectrum() - self._datasets[i].npred() + ) / self._datasets[i].npred() + residuals = residuals.data.reshape(-1) + + residual_plot.add_data( + x, + y, + residuals, + xerr=xerr, + label=self._datasets[i].name, + show_data=kwargs.get("show_data",True), + ) return residual_plot.finalize( xlabel="Energy\n(keV)", ylabel="Counts/keV", xscale="log", yscale="log", - show_legend=show_legend, + show_legend=kwargs.get("show_legend",True), ) From b153bb67544077d473761b2174b9740cd90c9b90 Mon Sep 17 00:00:00 2001 From: Tobias Preis Date: Fri, 19 Jun 2026 14:27:22 +0200 Subject: [PATCH 4/7] add: expected count rate --- gammapy_plugin/gammapy_like.py | 36 +++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/gammapy_plugin/gammapy_like.py b/gammapy_plugin/gammapy_like.py index 4bc7de4..4b0062f 100644 --- a/gammapy_plugin/gammapy_like.py +++ b/gammapy_plugin/gammapy_like.py @@ -283,27 +283,39 @@ def display_model( # TODO: have to ensure this is only for specturm data set valid residual_plot = ResidualPlot( - show_residuals=True, - ratio_residuals=True, **kwargs, ) for i in range(len(self._datasets)): + y_unweighted = self._datasets[i].counts.data.reshape(-1) - x = self._datasets[i].counts.geom.axes["energy"].as_plot_center.to("keV").value + x = ( + self._datasets[i] + .counts.geom.axes["energy"] + .as_plot_center.to("keV") + .value + ) xerr = [ - self._datasets[i].counts.geom.axes["energy"].as_plot_xerr[j].to("keV").value + self._datasets[i] + .counts.geom.axes["energy"] + .as_plot_xerr[j] + .to("keV") + .value for j in [0, 1] ] bins = ( - self._datasets[i].counts.geom.axes["energy"].as_plot_edges.to("keV").value + self._datasets[i] + .counts.geom.axes["energy"] + .as_plot_edges.to("keV") + .value ) widths = np.diff(bins) y = y_unweighted / widths residuals = ( - self._datasets[i].counts.get_spectrum() - self._datasets[i].npred() - ) / self._datasets[i].npred() + self._datasets[i].counts.get_spectrum() + - self._datasets[i].npred().get_spectrum() + ) / self._datasets[i].npred().get_spectrum() residuals = residuals.data.reshape(-1) residual_plot.add_data( @@ -312,7 +324,13 @@ def display_model( residuals, xerr=xerr, label=self._datasets[i].name, - show_data=kwargs.get("show_data",True), + show_data=kwargs.get("show_data", True), + ) + + residual_plot.add_model( + x, + self._datasets[i].npred().get_spectrum().data.reshape(-1) / widths, + label=kwargs.get("model_label", "Expected"), ) return residual_plot.finalize( @@ -320,5 +338,5 @@ def display_model( ylabel="Counts/keV", xscale="log", yscale="log", - show_legend=kwargs.get("show_legend",True), + show_legend=kwargs.get("show_legend", True), ) From 46101b900d08df59792137bd6a3882e451c4914d Mon Sep 17 00:00:00 2001 From: Tobias Preis Date: Fri, 19 Jun 2026 14:30:37 +0200 Subject: [PATCH 5/7] add: display_spectrum_model_counts --- examples/example_crab-spectrum_fit.md | 117 +------------------------- 1 file changed, 3 insertions(+), 114 deletions(-) diff --git a/examples/example_crab-spectrum_fit.md b/examples/example_crab-spectrum_fit.md index 4477765..ebad819 100644 --- a/examples/example_crab-spectrum_fit.md +++ b/examples/example_crab-spectrum_fit.md @@ -7,9 +7,8 @@ jupyter: format_version: '1.3' jupytext_version: 1.19.1 kernelspec: - display_name: docs_env + display_name: Python 3 language: python - name: docs_env --- ```python @@ -38,6 +37,7 @@ from gammapy_plugin.converter import AstromodelConverter from gammapy_plugin.gammapy_like import GammapyLike from gammapy_plugin.test.utils import get_close +from threeML.io.plotting.post_process_data_plots import display_spectrum_model_counts # get_units().energy = u.TeV ``` @@ -120,50 +120,6 @@ jl.fit() res = jl.results ``` -```python -import numpy as np -import matplotlib.pyplot as plt - -e = np.geomspace(1e-1, 1e2, 1000) * u.TeV -``` - -```python -def mohrmann_weird_spec(E, N0, E0, gamma1, gamma2, Ebreak, beta): - return ( - N0 - * np.power(E / E0, -gamma1) - * np.power(1 + np.power(E / Ebreak, (gamma2 - gamma1) / beta), -beta) - ) -``` - -```python -plt.plot(e, e**2 * ps(e).to("TeV-1 cm-2 s-1"), label="this shit") -plt.plot( - e, - e**2 - * mohrmann_weird_spec( - e, - 3.35 * 1e-10 * u.Unit("TeV-1 s-1 cm-2"), - 1.0 * u.TeV, - 1.61, - 2.95, - 0.33 * u.TeV, - 1.73, - ), - label="mohrmann", -) -plt.fill_between( - np.geomspace(0.4 * u.TeV, 50 * u.TeV, 100), - np.geomspace(0.4 * u.TeV, 50 * u.TeV, 100) ** 2 - * ps(np.geomspace(0.4, 50, 100) * u.TeV).to("TeV -1 cm-2 s-1"), - alpha=0.2, -) -# plt.ylim(3*10e-15,3*10e-11) -plt.xscale("log") -plt.yscale("log") -plt.legend() -``` - ```python plot_spectra( res, energy_unit=u.TeV, flux_unit=u.Unit("TeV-1 cm-2 s-1"), ene_min=0.4, ene_max=50 @@ -171,74 +127,7 @@ plot_spectra( ``` ```python -?plot_spectra -``` - -```python -dir(res) -``` - -```python -res.plot_chains() -``` - -```python -res.write_to("crappycrabby.fits", overwrite=True) -``` - -```python -from astropy.io import fits as fits -``` - -```python -with fits.open("crappycrabby.fits") as f: - print(*[f"{x}\n" for x in f[1].header["MODEL"].split("_NEWLINE_")]) -``` - -```python -from threeML import load_analysis_results -from astromodels import * - -# get_units().energy = u.keV -# get_units().to_dict() -``` - -```python -ar = load_analysis_results("crappycrabby.fits") -``` - -```python - -``` - -```python -ar -``` - -```python -plot_spectra(ar, ene_min=1e8, ene_max=1e10) -``` - -```python -ar.optimized_model.crab -``` - -```python -threeML_config.model_plot.point_source_plot["flux_unit"] = "1/(TeV s cm2)" -``` - -```python - -``` - -```python -threeML_config.model_plot.point_source_plot["ene_unit"] = "TeV" -``` - -```python -import numpy as np - -np.linspace(1 * u.TeV, 10 * u.TeV, 10) +display_spectrum_model_counts(jl) ``` ```python From b32c1f6acea5fab28a4a31d5597707993630a682 Mon Sep 17 00:00:00 2001 From: Tobias Preis Date: Fri, 19 Jun 2026 14:30:53 +0200 Subject: [PATCH 6/7] chore: update jupyter kernel descriptio --- examples/example_extended_source_fov_bkg.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/example_extended_source_fov_bkg.md b/examples/example_extended_source_fov_bkg.md index 5eaa552..52593cf 100644 --- a/examples/example_extended_source_fov_bkg.md +++ b/examples/example_extended_source_fov_bkg.md @@ -7,9 +7,8 @@ jupyter: format_version: '1.3' jupytext_version: 1.19.1 kernelspec: - display_name: Python (threeML) + display_name: Python 3 language: python - name: threeml --- # 3D Analysis of an Extended Source From 008b672a892014e68d42a06e4b169c9706d01476 Mon Sep 17 00:00:00 2001 From: Tobias Preis Date: Tue, 23 Jun 2026 15:43:48 +0200 Subject: [PATCH 7/7] add: support for multi spec components --- gammapy_plugin/models/spectral.py | 46 +++++++++++++++---------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/gammapy_plugin/models/spectral.py b/gammapy_plugin/models/spectral.py index 39931ed..64d1fb6 100644 --- a/gammapy_plugin/models/spectral.py +++ b/gammapy_plugin/models/spectral.py @@ -30,6 +30,7 @@ def __init__(self, function: Union[Function, list], **kwargs) -> None: else: self._x_unit = self._astromodel_function.x_unit self._y_unit = self._astromodel_function.y_unit + self._components_parameters = self._astromodel_function.parameters.values() elif isinstance(function, list): for f in function: assert isinstance( @@ -42,21 +43,27 @@ def __init__(self, function: Union[Function, list], **kwargs) -> None: self._components_parameters = [] for f in self._astromodel_function: - if x_unit is None and f.x_unit is not None: - x_unit = f.x_unit - elif x_unit is not None and f.x_unit is not None: - assert x_unit == f.x_unit, "Component x_unit not matching" - # TODO maybe transform also possible need to check + if not hasattr(f, "_requested_x_unit"): + if x_unit is None and f.x_unit is not None: + x_unit = f.x_unit + elif x_unit is not None and f.x_unit is not None: + assert x_unit == f.x_unit, "Component x_unit not matching" + # TODO: maybe transform also possible need to check + else: + raise ValueError(f"Your Component {f.name} has no x_unit") else: - raise ValueError(f"Your Component {f.name} has no x_unit") - - if y_unit is None and f.y_unit is not None: - y_unit = f.y_unit - elif y_unit is not None and f.y_unit is not None: - assert y_unit == f.y_unit, "Component y_unit not matching" - # TODO maybe transform also possible need to check - else: # pragma: no cover - raise ValueError(f"Your Component {f.name} has no y_unit") + x_unit = f._requested_x_unit + if not hasattr(f, "_requested_y_unit"): + if y_unit is None and f.y_unit is not None: + y_unit = f.y_unit + elif y_unit is not None and f.y_unit is not None: + assert y_unit == f.y_unit, "Component y_unit not matching" + # TODO maybe transform also possible need to check + else: # pragma: no cover + raise ValueError(f"Your Component {f.name} has no y_unit") + else: + y_unit = f._requested_y_unit + for p in f.parameters.values(): self._components_parameters.append(p) self._x_unit = x_unit @@ -83,11 +90,7 @@ def _setup_parameters(self): # needed later for correctly evaluating the function self._mapping = {} self._mapping_free = {} - parameter_dict = ( - self._astromodel_function.parameters.values() - if self._components_parameters is None - else self._components_parameters - ) + parameter_dict = self._components_parameters for v in parameter_dict: vmin = np.nan vmax = np.nan @@ -118,7 +121,6 @@ def evaluate(self, energy, **kwargs): shape = energy.shape energy = energy.flatten() if self._components is not None: - vals = [] if shape is None: for i in range(self._components): @@ -131,10 +133,6 @@ def evaluate(self, energy, **kwargs): return sum(vals) else: - kwargs_mapped = {} - for k, v in kwargs.items(): - if self._astromodel_function.path in k: - kwargs_mapped[k.split(f"{self._astromodel_function.path}.")[1]] = v if shape is None: return self._astromodel_function(energy) else: