diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 8e4ca36..99ca186 100644 --- a/.readthedocs.yaml +++ b/.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 diff --git a/docs/index.rst b/docs/index.rst index 737b3eb..d38b00f 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -44,6 +44,7 @@ To get started check out the `get-started`_ section. intro notebooks/converter.ipynb + notebooks/am_in_gp.ipynb api/API.rst @@ -51,4 +52,5 @@ To get started check out the `get-started`_ section. notebooks/crab_spectrum.ipynb notebooks/converter.ipynb + notebooks/am_in_gp.ipynb diff --git a/docs/md_docs/am_in_gp.md b/docs/md_docs/am_in_gp.md new file mode 100644 index 0000000..7e02932 --- /dev/null +++ b/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() +``` + diff --git a/gammapy_plugin/test/test_gammapy_fit.py b/gammapy_plugin/test/test_gammapy_fit.py new file mode 100644 index 0000000..6754960 --- /dev/null +++ b/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)