12 feat add support for threeml residualplot

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,125 +120,14 @@ 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
 )
 ```
 
 ```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

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

gammapy_plugin/gammapy_like.py

@@ -1,5 +1,5 @@
 import logging
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any, Dict, Optional, Union
 
 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,
+        **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(
+            **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
+            )
+            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().get_spectrum()
+            ) / self._datasets[i].npred().get_spectrum()
+            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),
+            )
+
+            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(
+            xlabel="Energy\n(keV)",
+            ylabel="Counts/keV",
+            xscale="log",
+            yscale="log",
+            show_legend=kwargs.get("show_legend", True),
+        )

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: