[FEAT] Update astromodels via gammapy

.readthedocs.yaml

@@ -12,7 +12,6 @@ build:
       - sphinx-apidoc . -o docs/api
       - bash docs/download_notebooks.sh
 
-
 # Build documentation in the "docs/" directory with Sphinx
 sphinx:
   configuration: docs/conf.py

docs/index.rst

@@ -44,11 +44,13 @@ To get started check out the `get-started`_ section.
 
    intro
    notebooks/converter.ipynb
+   notebooks/am_in_gp.ipynb
 
    api/API.rst
 
 .. nblinkgallery::
 
    notebooks/crab_spectrum.ipynb
    notebooks/converter.ipynb
+   notebooks/am_in_gp.ipynb
 

docs/md_docs/am_in_gp.md

@@ -0,0 +1,149 @@
+---
+jupyter:
+  jupytext:
+    text_representation:
+      extension: .md
+      format_name: markdown
+      format_version: '1.3'
+      jupytext_version: 1.19.3
+  kernelspec:
+    display_name: Python (threeml)
+    language: python
+    name: python3
+---
+
+# Using astromodels in gammapy
+
+In case you do not want to use `threeML` for fitting but want to use `astromodels` for
+modelling you can use the [converter](./converter.ipynb) to do so.
+
+```python
+from gammapy_plugin.converter import AstromodelConverter
+from gammapy_plugin.models import SpectralModelConverted
+from astromodels.functions import Powerlaw
+from astromodels.core.model import Model
+from astromodels.sources import PointSource
+import astropy.units as u
+from astromodels.core.units import get_units
+
+get_units().energy = u.TeV # this is HIGHLY EXPERIMENTAL!!!
+pl = Powerlaw()
+ps = PointSource("crab", ra=83.63, dec=22.01, spectral_shape=pl)
+pl.piv = 1 * u.TeV
+pl.K = 1e-12 * u.Unit("TeV-1 cm-2 s-1")
+pl.index = -2
+model = Model(ps)
+conv = AstromodelConverter(model)
+```
+
+
+pl = Powerlaw()
+ps = PointSource("crab",ra=83.63, dec=22.01,spectral_shape = pl)
+
+model = Model(ps)
+conv = AstromodelConverter(model)
+
+
+Okay now go to all the gammapy setup
+```python
+from pathlib import Path
+import astropy.units as u
+from astropy.coordinates import Angle, SkyCoord
+from regions import CircleSkyRegion
+
+# %matplotlib inline
+import matplotlib.pyplot as plt
+from IPython.display import display
+from gammapy.data import DataStore
+from gammapy.datasets import (
+    Datasets,
+    FluxPointsDataset,
+    SpectrumDataset,
+)
+from gammapy.estimators import FluxPointsEstimator
+from gammapy.estimators.utils import resample_energy_edges
+from gammapy.makers import (
+    ReflectedRegionsBackgroundMaker,
+    SafeMaskMaker,
+    SpectrumDatasetMaker,
+)
+from gammapy.maps import MapAxis, RegionGeom, WcsGeom
+from gammapy.modeling import Fit
+from gammapy.modeling.models import (
+    ExpCutoffPowerLawSpectralModel,
+    SkyModel,
+    create_crab_spectral_model,
+)
+from gammapy.visualization import plot_spectrum_datasets_off_regions
+
+datastore = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")
+obs_ids = [23523, 23526, 23559, 23592]
+observations = datastore.get_observations(obs_ids)
+target_position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs")
+on_region_radius = Angle("0.11 deg")
+on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)
+exclusion_region = CircleSkyRegion(
+    center=SkyCoord(183.604, -8.708, unit="deg", frame="galactic"),
+    radius=0.5 * u.deg,
+)
+
+skydir = target_position.galactic
+geom = WcsGeom.create(
+    npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs"
+)
+
+exclusion_mask = ~geom.region_mask([exclusion_region])
+exclusion_mask.plot()
+plt.show()
+energy_axis = MapAxis.from_energy_bounds(
+    0.1, 40, nbin=10, per_decade=True, unit="TeV", name="energy"
+)
+energy_axis_true = MapAxis.from_energy_bounds(
+    0.05, 100, nbin=20, per_decade=True, unit="TeV", name="energy_true"
+)
+
+geom = RegionGeom.create(region=on_region, axes=[energy_axis])
+dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true)
+
+dataset_maker = SpectrumDatasetMaker(
+    containment_correction=True, selection=["counts", "exposure", "edisp"]
+)
+bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
+safe_mask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
+
+datasets = Datasets()
+
+for obs_id, observation in zip(obs_ids, observations):
+    dataset = dataset_maker.run(dataset_empty.copy(name=str(obs_id)), observation)
+    dataset_on_off = bkg_maker.run(dataset, observation)
+    dataset_on_off = safe_mask_maker.run(dataset_on_off, observation)
+    datasets.append(dataset_on_off)
+
+print(datasets)
+plt.figure()
+ax = exclusion_mask.plot()
+on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="k")
+plot_spectrum_datasets_off_regions(ax=ax, datasets=datasets)
+plt.show()
+
+
+datasets.models = [conv.gammapy_models[0]]
+
+
+fit_joint = Fit()
+result_joint = fit_joint.run(datasets=datasets)
+```
+```python
+print(result_joint)
+```
+
+```python
+display(result_joint.models.to_parameters_table())
+```
+
+```python
+ax_spectrum, ax_residuals = datasets[0].plot_fit()
+ax_spectrum.set_ylim(0.1, 40)
+plt.show()
+```
+

gammapy_plugin/test/test_gammapy_fit.py

@@ -0,0 +1,80 @@
+import astropy.units as u
+import numpy as np
+import pytest
+from astromodels.core.model import Model
+from astromodels.functions import Powerlaw
+from astromodels.sources import PointSource
+from gammapy.modeling import Fit
+from gammapy.modeling.models import PowerLawSpectralModel, SkyModel
+
+from gammapy_plugin.converter import AstromodelConverter
+
+
+def test_gammapy_fit(crab_test_data):
+    datasets = crab_test_data.copy()
+    datasets_gp = crab_test_data.copy()
+    np.random.seed(1234)
+
+    pl = Powerlaw()
+    ps = PointSource("crab", spectral_shape=pl, ra=83.63, dec=22.01)
+    model = Model(ps)
+    conv = AstromodelConverter(model)
+
+    datasets.models = [conv.gammapy_models[0]]  # using ReflectedRegionsBackground
+
+    fit = Fit()
+    result = fit.run(datasets=datasets)
+
+    spectral_model = PowerLawSpectralModel(
+        amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"),
+        index=2,
+        reference=1 * u.TeV,
+    )
+    model_gp = SkyModel(spectral_model=spectral_model, name="crab")
+
+    datasets_gp.models = [model_gp]
+
+    fit_joint = Fit()
+    result_joint = fit_joint.run(datasets=datasets_gp)
+
+    assert np.isclose(
+        result.models[0].parameters["crab.spectrum.main.Powerlaw.index"].value,
+        -result_joint.models[0].parameters["index"].value,
+        rtol=0.1,
+    )
+
+
+def test_xspec_wrapping(crab_test_data):
+
+    pytest.importorskip("xspec")
+
+    from astromodels.xspec import XS_powerlaw
+
+    datasets = crab_test_data.copy()
+    datasets = datasets.stack_reduce()
+    datasets_gp = crab_test_data.copy()
+    datasets_gp = datasets_gp.stack_reduce()
+    np.random.seed(1234)
+    pl = XS_powerlaw()
+    ps = PointSource("crab", spectral_shape=pl, ra=83.633, dec=22.014)
+    pl.phoindex = 2
+    pl.phoindex.min_value = -np.nan
+    pl.phoindex.max_value = np.nan
+
+    model = Model(ps)
+    conv = AstromodelConverter(model)
+    datasets.models = [conv.gammapy_models[0]]  # using ReflectedRegionsBackground
+    fit = Fit()
+    result = fit.run(datasets=datasets)
+
+    spectral_model = PowerLawSpectralModel(
+        amplitude=100 * u.Unit("cm-2 s-1 keV-1"),
+        index=2,
+        reference=1 * u.keV,
+    )
+    model_gp = SkyModel(spectral_model=spectral_model, name="crab")
+
+    datasets_gp.models = [model_gp]
+
+    fit_joint = Fit()
+    result_joint = fit_joint.run(datasets=datasets_gp)