files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Fermi catalog and source classes."""
  import abc
  import warnings
  import numpy as np
  import astropy.units as u
  from astropy.table import Column, Table
  from astropy.time import Time
  from astropy.wcs import FITSFixedWarning
  from gammapy.estimators import FluxPoints, LightCurve
  from gammapy.modeling.models import (
      DiskSpatialModel,
      GaussianSpatialModel,
      Model,
      PointSpatialModel,
      SkyModel,
      TemplateSpatialModel,
  )
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.scripts import make_path
  from gammapy.utils.table import table_standardise_units_inplace
  from .core import SourceCatalog, SourceCatalogObject
  
  __all__ = [
      "SourceCatalogObject4FGL",
      "SourceCatalogObject3FGL",
      "SourceCatalogObject2FHL",
      "SourceCatalogObject3FHL",
      "SourceCatalog4FGL",
      "SourceCatalog3FGL",
      "SourceCatalog2FHL",
      "SourceCatalog3FHL",
  ]
  
  
  def compute_flux_points_ul(quantity, quantity_errp):
      """Compute UL value for fermi flux points.
  
      See https://arxiv.org/pdf/1501.02003.pdf (page 30)
      """
      return 2 * quantity_errp + quantity
  
  
  class SourceCatalogObjectFermiBase(SourceCatalogObject, abc.ABC):
      """Base class for Fermi-LAT catalogs."""
  
      asso = ["ASSOC1", "ASSOC2", "ASSOC_TEV", "ASSOC_GAM1", "ASSOC_GAM2", "ASSOC_GAM3"]
  
      def __str__(self):
          return self.info()
  
      def info(self, info="all"):
          """Summary info string.
  
          Parameters
          ----------
          info : {'all', 'basic', 'more', 'position', 'spectral','lightcurve'}
              Comma separated list of options
          """
          if info == "all":
              info = "basic,more,position,spectral,lightcurve"
  
          ss = ""
          ops = info.split(",")
          if "basic" in ops:
              ss += self._info_basic()
          if "more" in ops:
              ss += self._info_more()
          if "position" in ops:
              ss += self._info_position()
              if not self.is_pointlike:
                  ss += self._info_morphology()
          if "spectral" in ops:
              ss += self._info_spectral_fit()
              ss += self._info_spectral_points()
          if "lightcurve" in ops:
              ss += self._info_lightcurve()
          return ss
  
      def _info_basic(self):
          d = self.data
          keys = self.asso
          ss = "\n*** Basic info ***\n\n"
          ss += "Catalog row index (zero-based) : {}\n".format(self.row_index)
          ss += "{:<20s} : {}\n".format("Source name", self.name)
          if "Extended_Source_Name" in d:
              ss += "{:<20s} : {}\n".format("Extended name", d["Extended_Source_Name"])
  
          def get_nonentry_keys(keys):
              vals = [d[_].strip() for _ in keys]
              return ", ".join([_ for _ in vals if _ != ""])
  
          associations = get_nonentry_keys(keys)
          ss += "{:<16s} : {}\n".format("Associations", associations)
          try:
              ss += "{:<16s} : {:.3f}\n".format("ASSOC_PROB_BAY", d["ASSOC_PROB_BAY"])
              ss += "{:<16s} : {:.3f}\n".format("ASSOC_PROB_LR", d["ASSOC_PROB_LR"])
          except (KeyError):
              pass
          try:
              ss += "{:<16s} : {}\n".format("Class1", d["CLASS1"])
          except (KeyError):
              ss += "{:<16s} : {}\n".format("Class", d["CLASS"])
          try:
              ss += "{:<16s} : {}\n".format("Class2", d["CLASS2"])
          except (KeyError):
              pass
          ss += "{:<16s} : {}\n".format("TeVCat flag", d.get("TEVCAT_FLAG", "N/A"))
          return ss
  
      @abc.abstractmethod
      def _info_more(self):
          pass
  
      def _info_position(self):
          d = self.data
          ss = "\n*** Position info ***\n\n"
          ss += "{:<20s} : {:.3f}\n".format("RA", d["RAJ2000"])
          ss += "{:<20s} : {:.3f}\n".format("DEC", d["DEJ2000"])
          ss += "{:<20s} : {:.3f}\n".format("GLON", d["GLON"])
          ss += "{:<20s} : {:.3f}\n".format("GLAT", d["GLAT"])
  
          ss += "\n"
          ss += "{:<20s} : {:.4f}\n".format("Semimajor (68%)", d["Conf_68_SemiMajor"])
          ss += "{:<20s} : {:.4f}\n".format("Semiminor (68%)", d["Conf_68_SemiMinor"])
          ss += "{:<20s} : {:.2f}\n".format("Position angle (68%)", d["Conf_68_PosAng"])
          ss += "{:<20s} : {:.4f}\n".format("Semimajor (95%)", d["Conf_95_SemiMajor"])
          ss += "{:<20s} : {:.4f}\n".format("Semiminor (95%)", d["Conf_95_SemiMinor"])
          ss += "{:<20s} : {:.2f}\n".format("Position angle (95%)", d["Conf_95_PosAng"])
          ss += "{:<20s} : {:.0f}\n".format("ROI number", d["ROI_num"])
          return ss
  
      def _info_morphology(self):
          e = self.data_extended
          ss = "\n*** Extended source information ***\n\n"
          ss += "{:<16s} : {}\n".format("Model form", e["Model_Form"])
          ss += "{:<16s} : {:.4f}\n".format("Model semimajor", e["Model_SemiMajor"])
          ss += "{:<16s} : {:.4f}\n".format("Model semiminor", e["Model_SemiMinor"])
          ss += "{:<16s} : {:.4f}\n".format("Position angle", e["Model_PosAng"])
          try:
              ss += "{:<16s} : {}\n".format("Spatial function", e["Spatial_Function"])
          except KeyError:
              pass
          ss += "{:<16s} : {}\n\n".format("Spatial filename", e["Spatial_Filename"])
          return ss
  
      def _info_spectral_fit(self):
          return "\n"
  
      def _info_spectral_points(self):
          ss = "\n*** Spectral points ***\n\n"
          lines = self.flux_points.table_formatted.pformat(max_width=-1, max_lines=-1)
          ss += "\n".join(lines)
          return ss
  
      def _info_lightcurve(self):
          return "\n"
  
      @property
      def is_pointlike(self):
          return self.data["Extended_Source_Name"].strip() == ""
  
      # FIXME: this should be renamed `set_position_error`,
      # and `phi_0` isn't filled correctly, other parameters missing
      # see https://github.com/gammapy/gammapy/pull/2533#issuecomment-553329049
      def _set_spatial_errors(self, model):
          d = self.data
  
          if "Pos_err_68" in d:
              percent = 0.68
              semi_minor = d["Pos_err_68"]
              semi_major = d["Pos_err_68"]
              phi_0 = 0.0
          else:
              percent = 0.95
              semi_minor = d["Conf_95_SemiMinor"]
              semi_major = d["Conf_95_SemiMajor"]
              phi_0 = d["Conf_95_PosAng"]
  
          if np.isnan(phi_0):
              phi_0 = 0.0 * u.deg
  
          scale_1sigma = Gauss2DPDF().containment_radius(percent)
          lat_err = semi_major / scale_1sigma
          lon_err = semi_minor / scale_1sigma / np.cos(d["DEJ2000"])
  
          if  "TemplateSpatialModel" not in model.tag:
              model
  
  [... the rest is omitted ...]
files: gammapy/catalog/fermi.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Fermi catalog and source classes."""
  import abc
  import warnings
  import numpy as np
  import astropy.units as u
  from astropy.table import Column, Table
  from astropy.time import Time
  from astropy.wcs import FITSFixedWarning
  from gammapy.estimators import FluxPoints, LightCurve
  from gammapy.modeling.models import (
      DiskSpatialModel,
      GaussianSpatialModel,
      Model,
      PointSpatialModel,
      SkyModel,
      TemplateSpatialModel,
  )
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.scripts import make_path
  from gammapy.utils.table import table_standardise_units_inplace
  from .core import SourceCatalog, SourceCatalogObject
  
  __all__ = [
      "SourceCatalogObject4FGL",
      "SourceCatalogObject3FGL",
      "SourceCatalogObject2FHL",
      "SourceCatalogObject3FHL",
      "SourceCatalog4FGL",
      "SourceCatalog3FGL",
      "SourceCatalog2FHL",
      "SourceCatalog3FHL",
  ]
  
  
  def compute_flux_points_ul(quantity, quantity_errp):
      """Compute UL value for fermi flux points.
  
      See https://arxiv.org/pdf/1501.02003.pdf (page 30)
      """
      return 2 * quantity_errp + quantity
  
  
  class SourceCatalogObjectFermiBase(SourceCatalogObject, abc.ABC):
      """Base class for Fermi-LAT catalogs."""
  
      asso = ["ASSOC1", "ASSOC2", "ASSOC_TEV", "ASSOC_GAM1", "ASSOC_GAM2", "ASSOC_GAM3"]
  
      def __str__(self):
          return self.info()
  
      def info(self, info="all"):
          """Summary info string.
  
          Parameters
          ----------
          info : {'all', 'basic', 'more', 'position', 'spectral','lightcurve'}
              Comma separated list of options
          """
          if info == "all":
              info = "basic,more,position,spectral,lightcurve"
  
          ss = ""
          ops = info.split(",")
          if "basic" in ops:
              ss += self._info_basic()
          if "more" in ops:
              ss += self._info_more()
          if "position" in ops:
              ss += self._info_position()
              if not self.is_pointlike:
                  ss += self._info_morphology()
          if "spectral" in ops:
              ss += self._info_spectral_fit()
              ss += self._info_spectral_points()
          if "lightcurve" in ops:
              ss += self._info_lightcurve()
          return ss
  
      def _info_basic(self):
          d = self.data
          keys = self.asso
          ss = "\n*** Basic info ***\n\n"
          ss += "Catalog row index (zero-based) : {}\n".format(self.row_index)
          ss += "{:<20s} : {}\n".format("Source name", self.name)
          if "Extended_Source_Name" in d:
              ss += "{:<20s} : {}\n".format("Extended name", d["Extended_Source_Name"])
  
          def get_nonentry_keys(keys):
              vals = [d[_].strip() for _ in keys]
              return ", ".join([_ for _ in vals if _ != ""])
  
          associations = get_nonentry_keys(keys)
          ss += "{:<16s} : {}\n".format("Associations", associations)
          try:
              ss += "{:<16s} : {:.3f}\n".format("ASSOC_PROB_BAY", d["ASSOC_PROB_BAY"])
              ss += "{:<16s} : {:.3f}\n".format("ASSOC_PROB_LR", d["ASSOC_PROB_LR"])
          except (KeyError):
              pass
          try:
              ss += "{:<16s} : {}\n".format("Class1", d["CLASS1"])
          except (KeyError):
              ss += "{:<16s} : {}\n".format("Class", d["CLASS"])
          try:
              ss += "{:<16s} : {}\n".format("Class2", d["CLASS2"])
          except (KeyError):
              pass
          ss += "{:<16s} : {}\n".format("TeVCat flag", d.get("TEVCAT_FLAG", "N/A"))
          return ss
  
      @abc.abstractmethod
      def _info_more(self):
          pass
  
      def _info_position(self):
          d = self.data
          ss = "\n*** Position info ***\n\n"
          ss += "{:<20s} : {:.3f}\n".format("RA", d["RAJ2000"])
          ss += "{:<20s} : {:.3f}\n".format("DEC", d["DEJ2000"])
          ss += "{:<20s} : {:.3f}\n".format("GLON", d["GLON"])
          ss += "{:<20s} : {:.3f}\n".format("GLAT", d["GLAT"])
  
          ss += "\n"
          ss += "{:<20s} : {:.4f}\n".format("Semimajor (68%)", d["Conf_68_SemiMajor"])
          ss += "{:<20s} : {:.4f}\n".format("Semiminor (68%)", d["Conf_68_SemiMinor"])
          ss += "{:<20s} : {:.2f}\n".format("Position angle (68%)", d["Conf_68_PosAng"])
          ss += "{:<20s} : {:.4f}\n".format("Semimajor (95%)", d["Conf_95_SemiMajor"])
          ss += "{:<20s} : {:.4f}\n".format("Semiminor (95%)", d["Conf_95_SemiMinor"])
          ss += "{:<20s} : {:.2f}\n".format("Position angle (95%)", d["Conf_95_PosAng"])
          ss += "{:<20s} : {:.0f}\n".format("ROI number", d["ROI_num"])
          return ss
  
      def _info_morphology(self):
          e = self.data_extended
          ss = "\n*** Extended source information ***\n\n"
          ss += "{:<16s} : {}\n".format("Model form", e["Model_Form"])
          ss += "{:<16s} : {:.4f}\n".format("Model semimajor", e["Model_SemiMajor"])
          ss += "{:<16s} : {:.4f}\n".format("Model semiminor", e["Model_SemiMinor"])
          ss += "{:<16s} : {:.4f}\n".format("Position angle", e["Model_PosAng"])
          try:
              ss += "{:<16s} : {}\n".format("Spatial function", e["Spatial_Function"])
          except KeyError:
              pass
          ss += "{:<16s} : {}\n\n".format("Spatial filename", e["Spatial_Filename"])
          return ss
  
      def _info_spectral_fit(self):
          return "\n"
  
      def _info_spectral_points(self):
          ss = "\n*** Spectral points ***\n\n"
          lines = self.flux_points.table_formatted.pformat(max_width=-1, max_lines=-1)
          ss += "\n".join(lines)
          return ss
  
      def _info_lightcurve(self):
          return "\n"
  
      @property
      def is_pointlike(self):
          return self.data["Extended_Source_Name"].strip() == ""
  
      # FIXME: this should be renamed `set_position_error`,
      # and `phi_0` isn't filled correctly, other parameters missing
      # see https://github.com/gammapy/gammapy/pull/2533#issuecomment-553329049
      def _set_spatial_errors(self, model):
          d = self.data
  
          if "Pos_err_68" in d:
              percent = 0.68
              semi_minor = d["Pos_err_68"]
              semi_major = d["Pos_err_68"]
              phi_0 = 0.0
          else:
              percent = 0.95
              semi_minor = d["Conf_95_SemiMinor"]
              semi_major = d["Conf_95_SemiMajor"]
              phi_0 = d["Conf_95_PosAng"]
  
          if np.isnan(phi_0):
              phi_0 = 0.0 * u.deg
  
          scale_1sigma = Gauss2DPDF().containment_radius(percent)
          lat_err = semi_major / scale_1sigma
          lon_err = semi_minor / scale_1sigma / np.cos(d["DEJ2000"])
  
          if model.tag != "TemplateSpatialModel":
              model.parameters["
  
  [... the rest is omitted ...]
files: gammapy/catalog/fermi.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  from astropy import units as u
  from astropy.time import Time
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.catalog import (
      SourceCatalog2FHL,
      SourceCatalog3FGL,
      SourceCatalog3FHL,
      SourceCatalog4FGL,
  )
  from gammapy.modeling.models import (
      ExpCutoffPowerLaw3FGLSpectralModel,
      LogParabolaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      SuperExpCutoffPowerLaw3FGLSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
  )
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      assert_time_allclose,
      requires_data,
  )
  
  SOURCES_4FGL = [
      dict(
          idx=0,
          name="4FGL J0000.3-7355",
          str_ref_file="data/4fgl_J0000.3-7355.txt",
          spec_type=PowerLawSpectralModel,
          dnde=u.Quantity(2.9476e-11, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(5.3318e-12, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=3,
          name="4FGL J0001.5+2113",
          str_ref_file="data/4fgl_J0001.5+2113.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(2.8545e-8, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(1.3324e-9, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=7,
          name="4FGL J0002.8+6217",
          str_ref_file="data/4fgl_J0002.8+6217.txt",
          spec_type=SuperExpCutoffPowerLaw4FGLSpectralModel,
          dnde=u.Quantity(2.084e-09, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(1.0885e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=2718,
          name="4FGL J1409.1-6121e",
          str_ref_file="data/4fgl_J1409.1-6121e.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(1.3237202133031811e-12, "cm-2 s-1 MeV-1"),
          dnde_err=u.Quantity(4.513233455580648e-14, "cm-2 s-1 MeV-1"),
      ),
  ]
  
  SOURCES_3FGL = [
      dict(
          idx=0,
          name="3FGL J0000.1+6545",
          str_ref_file="data/3fgl_J0000.1+6545.txt",
          spec_type=PowerLawSpectralModel,
          dnde=u.Quantity(1.4351261e-9, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.1356270e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=4,
          name="3FGL J0001.4+2120",
          str_ref_file="data/3fgl_J0001.4+2120.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(8.3828599e-10, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.6713238e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=55,
          name="3FGL J0023.4+0923",
          str_ref_file="data/3fgl_J0023.4+0923.txt",
          spec_type=ExpCutoffPowerLaw3FGLSpectralModel,
          dnde=u.Quantity(1.8666925e-09, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.2068837e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=960,
          name="3FGL J0835.3-4510",
          str_ref_file="data/3fgl_J0835.3-4510.txt",
          spec_type=SuperExpCutoffPowerLaw3FGLSpectralModel,
          dnde=u.Quantity(1.6547128794756733e-06, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(1.6621504e-11, "cm-2 s-1 MeV-1"),
      ),
  ]
  
  SOURCES_2FHL = [
      dict(
          idx=221,
          name="2FHL J1445.1-0329",
          str_ref_file="data/2fhl_j1445.1-0329.txt",
          spec_type=PowerLaw2SpectralModel,
          dnde=u.Quantity(1.065463448091757e-10, "cm-2 s-1 TeV-1"),
          dnde_err=u.Quantity(4.9691205387540815e-11, "cm-2 s-1 TeV-1"),
      ),
      dict(
          idx=134,
          name="2FHL J0822.6-4250e",
          str_ref_file="data/2fhl_j0822.6-4250e.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(2.46548351696472e-10, "cm-2 s-1 TeV-1"),
          dnde_err=u.Quantity(9.771755529198772e-11, "cm-2 s-1 TeV-1"),
      ),
  ]
  
  SOURCES_3FHL = [
      dict(
          idx=352,
          name="3FHL J0534.5+2201",
          spec_type=PowerLawSpectralModel,
          dnde=u.Quantity(6.3848912826152664e-12, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.679593524691324e-13, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=1442,
          name="3FHL J2158.8-3013",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(2.056998292908196e-12, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(4.219030630302381e-13, "cm-2 s-1 GeV-1"),
      ),
  ]
  
  
  @requires_data()
  class TestFermi4FGLObject:
      @classmethod
      def setup_class(cls):
          cls.cat = SourceCatalog4FGL()
          cls.source_name = "4FGL J0534.5+2200"
          cls.source = cls.cat[cls.source_name]
  
      def test_name(self):
          assert self.source.name == self.source_name
  
      def test_row_index(self):
          assert self.source.row_index == 995
  
      @pytest.mark.parametrize("ref", SOURCES_4FGL, ids=lambda _: _["name"])
      def test_str(self, ref):
          actual = str(self.cat[ref["idx"]])
          expected = open(get_pkg_data_filename(ref["str_ref_file"])).read()
          assert actual == expected
  
      @pytest.mark.parametrize("ref", SOURCES_4FGL, ids=lambda _: _["name"])
      def test_spectral_model(self, ref):
          model = self.cat[ref["idx"]].spectral_model()
  
          e_ref = model.reference.quantity
          dnde, dnde_err = model.evaluate_error(e_ref)
          assert isinstance(model, ref["spec_type"])
          assert_quantity_allclose(dnde, ref["dnde"], rtol=1e-4)
          assert_quantity_allclose(dnde_err, ref["dnde_err"], rtol=1e-4)
  
      def test_spatial_model(self):
          model = self.cat["4FGL J0000.3-7355"].spatial_model()
          assert "PointSpatialModel" in model.tag
          assert model.frame == "icrs"
          p = model.parameters
          assert_allclose(p["lon_0"].value, 0.098
  
  [... the rest is omitted ...]
files: gammapy/catalog/tests/test_fermi.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  from astropy import units as u
  from astropy.time import Time
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.catalog import (
      SourceCatalog2FHL,
      SourceCatalog3FGL,
      SourceCatalog3FHL,
      SourceCatalog4FGL,
  )
  from gammapy.modeling.models import (
      ExpCutoffPowerLaw3FGLSpectralModel,
      LogParabolaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      SuperExpCutoffPowerLaw3FGLSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
  )
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      assert_time_allclose,
      requires_data,
  )
  
  SOURCES_4FGL = [
      dict(
          idx=0,
          name="4FGL J0000.3-7355",
          str_ref_file="data/4fgl_J0000.3-7355.txt",
          spec_type=PowerLawSpectralModel,
          dnde=u.Quantity(2.9476e-11, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(5.3318e-12, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=3,
          name="4FGL J0001.5+2113",
          str_ref_file="data/4fgl_J0001.5+2113.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(2.8545e-8, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(1.3324e-9, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=7,
          name="4FGL J0002.8+6217",
          str_ref_file="data/4fgl_J0002.8+6217.txt",
          spec_type=SuperExpCutoffPowerLaw4FGLSpectralModel,
          dnde=u.Quantity(2.084e-09, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(1.0885e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=2718,
          name="4FGL J1409.1-6121e",
          str_ref_file="data/4fgl_J1409.1-6121e.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(1.3237202133031811e-12, "cm-2 s-1 MeV-1"),
          dnde_err=u.Quantity(4.513233455580648e-14, "cm-2 s-1 MeV-1"),
      ),
  ]
  
  SOURCES_3FGL = [
      dict(
          idx=0,
          name="3FGL J0000.1+6545",
          str_ref_file="data/3fgl_J0000.1+6545.txt",
          spec_type=PowerLawSpectralModel,
          dnde=u.Quantity(1.4351261e-9, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.1356270e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=4,
          name="3FGL J0001.4+2120",
          str_ref_file="data/3fgl_J0001.4+2120.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(8.3828599e-10, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.6713238e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=55,
          name="3FGL J0023.4+0923",
          str_ref_file="data/3fgl_J0023.4+0923.txt",
          spec_type=ExpCutoffPowerLaw3FGLSpectralModel,
          dnde=u.Quantity(1.8666925e-09, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.2068837e-10, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=960,
          name="3FGL J0835.3-4510",
          str_ref_file="data/3fgl_J0835.3-4510.txt",
          spec_type=SuperExpCutoffPowerLaw3FGLSpectralModel,
          dnde=u.Quantity(1.6547128794756733e-06, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(1.6621504e-11, "cm-2 s-1 MeV-1"),
      ),
  ]
  
  SOURCES_2FHL = [
      dict(
          idx=221,
          name="2FHL J1445.1-0329",
          str_ref_file="data/2fhl_j1445.1-0329.txt",
          spec_type=PowerLaw2SpectralModel,
          dnde=u.Quantity(1.065463448091757e-10, "cm-2 s-1 TeV-1"),
          dnde_err=u.Quantity(4.9691205387540815e-11, "cm-2 s-1 TeV-1"),
      ),
      dict(
          idx=134,
          name="2FHL J0822.6-4250e",
          str_ref_file="data/2fhl_j0822.6-4250e.txt",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(2.46548351696472e-10, "cm-2 s-1 TeV-1"),
          dnde_err=u.Quantity(9.771755529198772e-11, "cm-2 s-1 TeV-1"),
      ),
  ]
  
  SOURCES_3FHL = [
      dict(
          idx=352,
          name="3FHL J0534.5+2201",
          spec_type=PowerLawSpectralModel,
          dnde=u.Quantity(6.3848912826152664e-12, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(2.679593524691324e-13, "cm-2 s-1 GeV-1"),
      ),
      dict(
          idx=1442,
          name="3FHL J2158.8-3013",
          spec_type=LogParabolaSpectralModel,
          dnde=u.Quantity(2.056998292908196e-12, "cm-2 s-1 GeV-1"),
          dnde_err=u.Quantity(4.219030630302381e-13, "cm-2 s-1 GeV-1"),
      ),
  ]
  
  
  @requires_data()
  class TestFermi4FGLObject:
      @classmethod
      def setup_class(cls):
          cls.cat = SourceCatalog4FGL()
          cls.source_name = "4FGL J0534.5+2200"
          cls.source = cls.cat[cls.source_name]
  
      def test_name(self):
          assert self.source.name == self.source_name
  
      def test_row_index(self):
          assert self.source.row_index == 995
  
      @pytest.mark.parametrize("ref", SOURCES_4FGL, ids=lambda _: _["name"])
      def test_str(self, ref):
          actual = str(self.cat[ref["idx"]])
          expected = open(get_pkg_data_filename(ref["str_ref_file"])).read()
          assert actual == expected
  
      @pytest.mark.parametrize("ref", SOURCES_4FGL, ids=lambda _: _["name"])
      def test_spectral_model(self, ref):
          model = self.cat[ref["idx"]].spectral_model()
  
          e_ref = model.reference.quantity
          dnde, dnde_err = model.evaluate_error(e_ref)
          assert isinstance(model, ref["spec_type"])
          assert_quantity_allclose(dnde, ref["dnde"], rtol=1e-4)
          assert_quantity_allclose(dnde_err, ref["dnde_err"], rtol=1e-4)
  
      def test_spatial_model(self):
          model = self.cat["4FGL J0000.3-7355"].spatial_model()
          assert model.tag == "PointSpatialModel"
          assert model.frame == "icrs"
          p = model.parameters
          assert_allclose(p["lon_0"].value, 0.0983
  
  [... the rest is omitted ...]
files: gammapy/catalog/tests/test_fermi.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import collections.abc
  import copy
  from os.path import split
  import numpy as np
  import astropy.units as u
  from astropy.table import Table
  import yaml
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.scripts import make_name, make_path
  
  
  def _set_link(shared_register, model):
      for param in model.parameters:
          name = param.name
          link_label = param._link_label_io
          if link_label is not None:
              if link_label in shared_register:
                  new_param = shared_register[link_label]
                  setattr(model, name, new_param)
              else:
                  shared_register[link_label] = param
      return shared_register
  
  
  __all__ = ["Model", "Models", "ProperModels"]
  
  
  class Model:
      """Model base class."""
  
      def __init__(self, **kwargs):
          # Copy default parameters from the class to the instance
          default_parameters = self.default_parameters.copy()
  
          for par in default_parameters:
              value = kwargs.get(par.name, par)
  
              if not isinstance(value, Parameter):
                  par.quantity = u.Quantity(value)
              else:
                  par = value
  
              setattr(self, par.name, par)
  
          self._covariance = Covariance(self.parameters)
  
      def __init_subclass__(cls, **kwargs):
          # Add parameters list on the model sub-class (not instances)
          cls.default_parameters = Parameters(
              [_ for _ in cls.__dict__.values() if isinstance(_, Parameter)]
          )
  
      @classmethod
      def from_parameters(cls, parameters, **kwargs):
          """Create model from parameter list
  
          Parameters
          ----------
          parameters : `Parameters`
              Parameters for init
  
          Returns
          -------
          model : `Model`
              Model instance
          """
          for par in parameters:
              kwargs[par.name] = par
          return cls(**kwargs)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance(self.parameters)
  
      @property
      def covariance(self):
          self._check_covariance()
          for par in self.parameters:
              pars = Parameters([par])
              covar = Covariance(pars, data=[[par.error ** 2]])
              self._covariance.set_subcovariance(covar)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for par in self.parameters:
              pars = Parameters([par])
              variance = self._covariance.get_subcovariance(pars)
              par.error = np.sqrt(variance)
  
      @property
      def parameters(self):
          """Parameters (`~gammapy.modeling.Parameters`)"""
          return Parameters(
              [getattr(self, name) for name in self.default_parameters.names]
          )
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      def to_dict(self):
          """Create dict for YAML serialisation"""
          tag = self.tag[0] if isinstance(self.tag, list) else self.tag
          return {"type": tag, "parameters": self.parameters.to_dict()}
  
      @classmethod
      def from_dict(cls, data):
          kwargs = {}
          parameters = Parameters.from_dict(data["parameters"])
  
          # TODO: this is a special case for spatial models, maybe better move to `SpatialModel` base class
          if "frame" in data:
              kwargs["frame"] = data["frame"]
  
          return cls.from_parameters(parameters, **kwargs)
  
      @staticmethod
      def create(tag, *args, **kwargs):
          """Create a model instance.
  
          Examples
          --------
          >>> from gammapy.modeling.models import Model
          >>> spectral_model = Model.create("PowerLaw2SpectralModel", amplitude="1e-10 cm-2 s-1", index=3)
          >>> type(spectral_model)
          gammapy.modeling.models.spectral.PowerLaw2SpectralModel
          """
          from . import MODEL_REGISTRY
  
          cls = MODEL_REGISTRY.get_cls(tag)
          return cls(*args, **kwargs)
  
      def __str__(self):
          return f"{self.__class__.__name__}\n\n{self.parameters.to_table()}"
  
  
  class Models(collections.abc.MutableSequence):
      """Sky model collection.
  
      Parameters
      ----------
      models : `SkyModel`, list of `SkyModel` or `Models`
          Sky models
      """
  
      def __init__(self, models=None):
          if models is None:
              models = []
  
          if isinstance(models, Models):
              models = models._models
          elif isinstance(models, Model):
              models = [models]
          elif not isinstance(models, list):
              raise TypeError(f"Invalid type: {models!r}")
  
          unique_names = []
          for model in models:
              if model.name in unique_names:
                  raise (ValueError("Model names must be unique"))
              unique_names.append(model.name)
  
          self._models = models
          self._covar_file = None
          self._covariance = Covariance(self.parameters)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in self._models]
              )
  
      @property
      def covariance(self):
          self._check_covariance()
  
          for model in self._models:
              self._covariance.set_subcovariance(model.covariance)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for model in self._models:
              subcovar = self._covariance.get_subcovariance(model.covariance.parameters)
              model.covariance = subcovar
  
      @property
      def parameters(self):
          return Parameters.from_stack([_.parameters for _ in self._models])
  
      @property
      def parameters_unique_names(self):
          from gammapy.modeling.models import SkyModel
  
          param_names = []
          for m in self._models:
              if isinstance(m, SkyModel):
                  for p in m.parameters:
                      if (
                          m.spectral_model is not None
                          and p in m.spectral_model.parameters
                      ):
                          tag = ".spectral."
                      elif (
                          m.spatial_model is not None and p in m.spatial_model.parameters
                      ):
                          tag = ".spatial."
                      elif (
                          m.temporal_model is not None
                          and p in m.temporal_model.parameters
                      ):
                          tag = ".temporal."
                      param_names.append(m.name + tag + p.name)
              else:
                  for p in m.parameters:
                      param_names.append(m.name + "." + p.name)
          return param_names
  
      @property
      def names(self):
          return [m.name for m in self._models]
  
      @classmethod
      def read(cls, filename):
          """Read from YAML file."""
          yaml_str = make_path(filename).read_text()
          path, filename = split(filename)
          return cls.from_yaml(yaml_str, path=path)
  
      @classmethod
      def from_yaml(cls, yaml_str, path=""):
          """Create from YAML string."""
          data = yaml.safe_load(yaml_str)
          return cls.from_dict(data, path=path)
  
      @classmethod
      def from_dict(cls, data, path=""):
          """Create from dict."""
          from . import MODEL_REGISTRY, SkyModel
  
          models = []
  
          for component in data["components"]:
              model = MODEL_REGISTRY.get_cls(component["type"]).from_dict(component)
              models.append(model)
  
          models = cls(models)
  
          if "covariance" in data:
              filename = data["covariance"]
              path = make_path(path)
              if not (path / filename).exists():
                  path, filename = split(filename)
  
              models.read_covariance(path, filename, format="ascii.fixed_width")
  
          shared_register = {}
          for model in models:
              if isinstance(model, SkyModel):
                  submodels = [
                      model.spectral_model,
                      model.spatial_model,
                      model.temporal_model,
                  ]
                  for submodel in submodels:
                      if submodel is not None:
                          shared_register = _set_link(shared_register, submodel)
              else:
                  shared_register = _set_link(shared_register, model)
          return models
  
      def write(self, path, overwrite=False, write_covariance=True):
          """Write to YAML file.
  
          Parameters
          ----------
          path : `pathlib.Path`
              path to write files
          overwrite : bool
              overwrite files
          write_covariance : bool
              save covariance or not
          """
  
          base_path, _ = split(path)
          path = make_path(path)
          base_path = make_path(base_path)
  
          if path.exists
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/core.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import collections.abc
  import copy
  from os.path import split
  import numpy as np
  import astropy.units as u
  from astropy.table import Table
  import yaml
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.scripts import make_name, make_path
  
  
  def _set_link(shared_register, model):
      for param in model.parameters:
          name = param.name
          link_label = param._link_label_io
          if link_label is not None:
              if link_label in shared_register:
                  new_param = shared_register[link_label]
                  setattr(model, name, new_param)
              else:
                  shared_register[link_label] = param
      return shared_register
  
  
  __all__ = ["Model", "Models", "ProperModels"]
  
  
  class Model:
      """Model base class."""
  
      def __init__(self, **kwargs):
          # Copy default parameters from the class to the instance
          default_parameters = self.default_parameters.copy()
  
          for par in default_parameters:
              value = kwargs.get(par.name, par)
  
              if not isinstance(value, Parameter):
                  par.quantity = u.Quantity(value)
              else:
                  par = value
  
              setattr(self, par.name, par)
  
          self._covariance = Covariance(self.parameters)
  
      def __init_subclass__(cls, **kwargs):
          # Add parameters list on the model sub-class (not instances)
          cls.default_parameters = Parameters(
              [_ for _ in cls.__dict__.values() if isinstance(_, Parameter)]
          )
  
      @classmethod
      def from_parameters(cls, parameters, **kwargs):
          """Create model from parameter list
  
          Parameters
          ----------
          parameters : `Parameters`
              Parameters for init
  
          Returns
          -------
          model : `Model`
              Model instance
          """
          for par in parameters:
              kwargs[par.name] = par
          return cls(**kwargs)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance(self.parameters)
  
      @property
      def covariance(self):
          self._check_covariance()
          for par in self.parameters:
              pars = Parameters([par])
              covar = Covariance(pars, data=[[par.error ** 2]])
              self._covariance.set_subcovariance(covar)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for par in self.parameters:
              pars = Parameters([par])
              variance = self._covariance.get_subcovariance(pars)
              par.error = np.sqrt(variance)
  
      @property
      def parameters(self):
          """Parameters (`~gammapy.modeling.Parameters`)"""
          return Parameters(
              [getattr(self, name) for name in self.default_parameters.names]
          )
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      def to_dict(self):
          """Create dict for YAML serialisation"""
          return {"type": self.tag, "parameters": self.parameters.to_dict()}
  
      @classmethod
      def from_dict(cls, data):
          kwargs = {}
          parameters = Parameters.from_dict(data["parameters"])
  
          # TODO: this is a special case for spatial models, maybe better move to `SpatialModel` base class
          if "frame" in data:
              kwargs["frame"] = data["frame"]
  
          return cls.from_parameters(parameters, **kwargs)
  
      @staticmethod
      def create(tag, *args, **kwargs):
          """Create a model instance.
  
          Examples
          --------
          >>> from gammapy.modeling.models import Model
          >>> spectral_model = Model.create("PowerLaw2SpectralModel", amplitude="1e-10 cm-2 s-1", index=3)
          >>> type(spectral_model)
          gammapy.modeling.models.spectral.PowerLaw2SpectralModel
          """
          from . import MODEL_REGISTRY
  
          cls = MODEL_REGISTRY.get_cls(tag)
          return cls(*args, **kwargs)
  
      def __str__(self):
          return f"{self.__class__.__name__}\n\n{self.parameters.to_table()}"
  
  
  class Models(collections.abc.MutableSequence):
      """Sky model collection.
  
      Parameters
      ----------
      models : `SkyModel`, list of `SkyModel` or `Models`
          Sky models
      """
  
      def __init__(self, models=None):
          if models is None:
              models = []
  
          if isinstance(models, Models):
              models = models._models
          elif isinstance(models, Model):
              models = [models]
          elif not isinstance(models, list):
              raise TypeError(f"Invalid type: {models!r}")
  
          unique_names = []
          for model in models:
              if model.name in unique_names:
                  raise (ValueError("Model names must be unique"))
              unique_names.append(model.name)
  
          self._models = models
          self._covar_file = None
          self._covariance = Covariance(self.parameters)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in self._models]
              )
  
      @property
      def covariance(self):
          self._check_covariance()
  
          for model in self._models:
              self._covariance.set_subcovariance(model.covariance)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for model in self._models:
              subcovar = self._covariance.get_subcovariance(model.covariance.parameters)
              model.covariance = subcovar
  
      @property
      def parameters(self):
          return Parameters.from_stack([_.parameters for _ in self._models])
  
      @property
      def parameters_unique_names(self):
          from gammapy.modeling.models import SkyModel
  
          param_names = []
          for m in self._models:
              if isinstance(m, SkyModel):
                  for p in m.parameters:
                      if (
                          m.spectral_model is not None
                          and p in m.spectral_model.parameters
                      ):
                          tag = ".spectral."
                      elif (
                          m.spatial_model is not None and p in m.spatial_model.parameters
                      ):
                          tag = ".spatial."
                      elif (
                          m.temporal_model is not None
                          and p in m.temporal_model.parameters
                      ):
                          tag = ".temporal."
                      param_names.append(m.name + tag + p.name)
              else:
                  for p in m.parameters:
                      param_names.append(m.name + "." + p.name)
          return param_names
  
      @property
      def names(self):
          return [m.name for m in self._models]
  
      @classmethod
      def read(cls, filename):
          """Read from YAML file."""
          yaml_str = make_path(filename).read_text()
          path, filename = split(filename)
          return cls.from_yaml(yaml_str, path=path)
  
      @classmethod
      def from_yaml(cls, yaml_str, path=""):
          """Create from YAML string."""
          data = yaml.safe_load(yaml_str)
          return cls.from_dict(data, path=path)
  
      @classmethod
      def from_dict(cls, data, path=""):
          """Create from dict."""
          from . import MODEL_REGISTRY, SkyModel
  
          models = []
  
          for component in data["components"]:
              model = MODEL_REGISTRY.get_cls(component["type"]).from_dict(component)
              models.append(model)
  
          models = cls(models)
  
          if "covariance" in data:
              filename = data["covariance"]
              path = make_path(path)
              if not (path / filename).exists():
                  path, filename = split(filename)
  
              models.read_covariance(path, filename, format="ascii.fixed_width")
  
          shared_register = {}
          for model in models:
              if isinstance(model, SkyModel):
                  submodels = [
                      model.spectral_model,
                      model.spatial_model,
                      model.temporal_model,
                  ]
                  for submodel in submodels:
                      if submodel is not None:
                          shared_register = _set_link(shared_register, submodel)
              else:
                  shared_register = _set_link(shared_register, model)
          return models
  
      def write(self, path, overwrite=False, write_covariance=True):
          """Write to YAML file.
  
          Parameters
          ----------
          path : `pathlib.Path`
              path to write files
          overwrite : bool
              overwrite files
          write_covariance : bool
              save covariance or not
          """
  
          base_path, _ = split(path)
          path = make_path(path)
          base_path = make_path(base_path)
  
          if path.exists() and not overwrite:
              raise IOError(f"File exists already: {path
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/core.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spatial models."""
  import logging
  import numpy as np
  import scipy.integrate
  import scipy.special
  import astropy.units as u
  from astropy.coordinates import Angle, SkyCoord
  from astropy.coordinates.angle_utilities import angular_separation, position_angle
  from regions import (
      CircleAnnulusSkyRegion,
      EllipseSkyRegion,
      PointSkyRegion,
      PolygonSkyRegion,
  )
  from gammapy.maps import Map, WcsGeom
  from gammapy.modeling import Parameter, Parameters
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  log = logging.getLogger(__name__)
  
  
  def compute_sigma_eff(lon_0, lat_0, lon, lat, phi, major_axis, e):
      """Effective radius, used for the evaluation of elongated models"""
      phi_0 = position_angle(lon_0, lat_0, lon, lat)
      d_phi = phi - phi_0
      minor_axis = Angle(major_axis * np.sqrt(1 - e ** 2))
  
      a2 = (major_axis * np.sin(d_phi)) ** 2
      b2 = (minor_axis * np.cos(d_phi)) ** 2
      denominator = np.sqrt(a2 + b2)
      sigma_eff = major_axis * minor_axis / denominator
      return minor_axis, sigma_eff
  
  
  class SpatialModel(Model):
      """Spatial model base class."""
  
      def __init__(self, **kwargs):
          frame = kwargs.pop("frame", "icrs")
          super().__init__(**kwargs)
          if not hasattr(self, "frame"):
              self.frame = frame
  
      def __call__(self, lon, lat):
          """Call evaluate method"""
          kwargs = {par.name: par.quantity for par in self.parameters}
          return self.evaluate(lon, lat, **kwargs)
  
      @property
      def position(self):
          """Spatial model center position"""
          lon = self.lon_0.quantity
          lat = self.lat_0.quantity
          return SkyCoord(lon, lat, frame=self.frame)
  
      @position.setter
      def position(self, skycoord):
          """Spatial model center position"""
          coord = skycoord.transform_to(self.frame)
          self.lon_0.quantity = coord.data.lon
          self.lat_0.quantity = coord.data.lat
  
      # TODO: get rid of this!
      _phi_0 = 0.0
  
      @property
      def phi_0(self):
          return self._phi_0
  
      @phi_0.setter
      def phi_0(self, phi_0=0.0):
          self._phi_0 = phi_0
  
      @property
      def position_error(self):
          """Get 95% containment position error as (`~regions.EllipseSkyRegion`)"""
          if self.covariance is None:
              return EllipseSkyRegion(
                  center=self.position,
                  height=np.nan * u.deg,
                  width=np.nan * u.deg,
                  angle=np.nan * u.deg,
              )
          pars = self.parameters
          sub_covar = self.covariance.get_subcovariance(["lon_0", "lat_0"]).data.copy()
          cos_lat = np.cos(self.lat_0.quantity.to_value("rad"))
          sub_covar[0, 0] *= cos_lat ** 2.0
          sub_covar[0, 1] *= cos_lat
          sub_covar[1, 0] *= cos_lat
          eig_vals, eig_vecs = np.linalg.eig(sub_covar)
          lon_err, lat_err = np.sqrt(eig_vals)
          y_vec = eig_vecs[:, 0]
          phi = (np.arctan2(y_vec[1], y_vec[0]) * u.rad).to("deg") + self.phi_0
          err = np.sort([lon_err, lat_err])
          scale_r95 = Gauss2DPDF().containment_radius(0.95)
          err *= scale_r95
          if err[1] == lon_err * scale_r95:
              phi += 90 * u.deg
              height = 2 * err[1] * pars["lon_0"].unit
              width = 2 * err[0] * pars["lat_0"].unit
          else:
              height = 2 * err[1] * pars["lat_0"].unit
              width = 2 * err[0] * pars["lon_0"].unit
          return EllipseSkyRegion(
              center=self.position, height=height, width=width, angle=phi
          )
  
      def evaluate_geom(self, geom):
          """Evaluate model on `~gammapy.maps.Geom`."""
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat)
  
      def integrate_geom(self, geom):
          """Integrate model on `~gammapy.maps.Geom`."""
          values = self.evaluate_geom(geom)
          data = values * geom.solid_angle()
          return Map.from_geom(geom=geom, data=data.value, unit=data.unit)
  
      def to_dict(self):
          """Create dict for YAML serilisation"""
          data = super().to_dict()
          data["frame"] = self.frame
          data["parameters"] = data.pop("parameters")
          return data
  
      def plot(self, ax=None, geom=None, **kwargs):
          """Plot spatial model.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          geom : `~gammapy.maps.WcsGeom`, optional
              Geom to use for plotting.
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          if self.evaluation_radius is None and geom is None:
              raise ValueError(
                  f"{self.__class__.__name__} requires geom to be defined for plotting."
              )
  
          if geom is None:
              width = 2 * max(self.evaluation_radius, 0.1 * u.deg)
              geom = WcsGeom.create(
                  skydir=self.position, frame=self.frame, width=width, binsz=0.02
              )
  
          data = self.evaluate_geom(geom)
          m = Map.from_geom(geom, data=data.value, unit=data.unit)
          _, ax, _ = m.plot(ax=ax, **kwargs)
  
          return ax
  
      def plot_error(self, ax=None, **kwargs):
          """Plot position error
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          # plot center position
          lon, lat = self.lon_0.value, self.lat_0.value
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("marker", "x")
          kwargs.setdefault("color", "red")
          kwargs.setdefault("label", "position")
  
          ax.scatter(lon, lat, transform=ax.get_transform(self.frame), **kwargs)
  
          # plot position error
          if not np.all(self.covariance.data == 0):
              region = self.position_error.to_pixel(ax.wcs)
              artist = region.as_artist(facecolor="none", edgecolor=kwargs["color"])
              ax.add_artist(artist)
  
          return ax
  
  
  class PointSpatialModel(SpatialModel):
      r"""Point Source.
  
      For more information see :ref:`point-spatial-model`.
  
      Parameters
      ----------
      lon_0, lat_0 : `~astropy.coordinates.Angle`
          Center position
      frame : {"icrs", "galactic"}
          Center position coordinate frame
      """
  
      tag = ["PointSpatialModel", "PS"]
      lon_0 = Parameter("lon_0", "0 deg")
      lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
  
      @property
      def evaluation_radius(self):
          """Evaluation radius (`~astropy.coordinates.Angle`).
  
          Set as zero degrees.
          """
          return 0 * u.deg
  
      @staticmethod
      def _grid_weights(x, y, x0, y0):
          """Compute 4-pixel weights such that centroid is preserved."""
          dx = np.abs(x - x0)
          dx = np.where(dx < 1, 1 - dx, 
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/spatial.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spatial models."""
  import logging
  import numpy as np
  import scipy.integrate
  import scipy.special
  import astropy.units as u
  from astropy.coordinates import Angle, SkyCoord
  from astropy.coordinates.angle_utilities import angular_separation, position_angle
  from regions import (
      CircleAnnulusSkyRegion,
      EllipseSkyRegion,
      PointSkyRegion,
      PolygonSkyRegion,
  )
  from gammapy.maps import Map, WcsGeom
  from gammapy.modeling import Parameter, Parameters
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  log = logging.getLogger(__name__)
  
  
  def compute_sigma_eff(lon_0, lat_0, lon, lat, phi, major_axis, e):
      """Effective radius, used for the evaluation of elongated models"""
      phi_0 = position_angle(lon_0, lat_0, lon, lat)
      d_phi = phi - phi_0
      minor_axis = Angle(major_axis * np.sqrt(1 - e ** 2))
  
      a2 = (major_axis * np.sin(d_phi)) ** 2
      b2 = (minor_axis * np.cos(d_phi)) ** 2
      denominator = np.sqrt(a2 + b2)
      sigma_eff = major_axis * minor_axis / denominator
      return minor_axis, sigma_eff
  
  
  class SpatialModel(Model):
      """Spatial model base class."""
  
      def __init__(self, **kwargs):
          frame = kwargs.pop("frame", "icrs")
          super().__init__(**kwargs)
          if not hasattr(self, "frame"):
              self.frame = frame
  
      def __call__(self, lon, lat):
          """Call evaluate method"""
          kwargs = {par.name: par.quantity for par in self.parameters}
          return self.evaluate(lon, lat, **kwargs)
  
      @property
      def position(self):
          """Spatial model center position"""
          lon = self.lon_0.quantity
          lat = self.lat_0.quantity
          return SkyCoord(lon, lat, frame=self.frame)
  
      @position.setter
      def position(self, skycoord):
          """Spatial model center position"""
          coord = skycoord.transform_to(self.frame)
          self.lon_0.quantity = coord.data.lon
          self.lat_0.quantity = coord.data.lat
  
      # TODO: get rid of this!
      _phi_0 = 0.0
  
      @property
      def phi_0(self):
          return self._phi_0
  
      @phi_0.setter
      def phi_0(self, phi_0=0.0):
          self._phi_0 = phi_0
  
      @property
      def position_error(self):
          """Get 95% containment position error as (`~regions.EllipseSkyRegion`)"""
          if self.covariance is None:
              return EllipseSkyRegion(
                  center=self.position,
                  height=np.nan * u.deg,
                  width=np.nan * u.deg,
                  angle=np.nan * u.deg,
              )
          pars = self.parameters
          sub_covar = self.covariance.get_subcovariance(["lon_0", "lat_0"]).data.copy()
          cos_lat = np.cos(self.lat_0.quantity.to_value("rad"))
          sub_covar[0, 0] *= cos_lat ** 2.0
          sub_covar[0, 1] *= cos_lat
          sub_covar[1, 0] *= cos_lat
          eig_vals, eig_vecs = np.linalg.eig(sub_covar)
          lon_err, lat_err = np.sqrt(eig_vals)
          y_vec = eig_vecs[:, 0]
          phi = (np.arctan2(y_vec[1], y_vec[0]) * u.rad).to("deg") + self.phi_0
          err = np.sort([lon_err, lat_err])
          scale_r95 = Gauss2DPDF().containment_radius(0.95)
          err *= scale_r95
          if err[1] == lon_err * scale_r95:
              phi += 90 * u.deg
              height = 2 * err[1] * pars["lon_0"].unit
              width = 2 * err[0] * pars["lat_0"].unit
          else:
              height = 2 * err[1] * pars["lat_0"].unit
              width = 2 * err[0] * pars["lon_0"].unit
          return EllipseSkyRegion(
              center=self.position, height=height, width=width, angle=phi
          )
  
      def evaluate_geom(self, geom):
          """Evaluate model on `~gammapy.maps.Geom`."""
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat)
  
      def integrate_geom(self, geom):
          """Integrate model on `~gammapy.maps.Geom`."""
          values = self.evaluate_geom(geom)
          data = values * geom.solid_angle()
          return Map.from_geom(geom=geom, data=data.value, unit=data.unit)
  
      def to_dict(self):
          """Create dict for YAML serilisation"""
          data = super().to_dict()
          data["frame"] = self.frame
          data["parameters"] = data.pop("parameters")
          return data
  
      def plot(self, ax=None, geom=None, **kwargs):
          """Plot spatial model.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          geom : `~gammapy.maps.WcsGeom`, optional
              Geom to use for plotting.
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          if self.evaluation_radius is None and geom is None:
              raise ValueError(
                  f"{self.__class__.__name__} requires geom to be defined for plotting."
              )
  
          if geom is None:
              width = 2 * max(self.evaluation_radius, 0.1 * u.deg)
              geom = WcsGeom.create(
                  skydir=self.position, frame=self.frame, width=width, binsz=0.02
              )
  
          data = self.evaluate_geom(geom)
          m = Map.from_geom(geom, data=data.value, unit=data.unit)
          _, ax, _ = m.plot(ax=ax, **kwargs)
  
          return ax
  
      def plot_error(self, ax=None, **kwargs):
          """Plot position error
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          # plot center position
          lon, lat = self.lon_0.value, self.lat_0.value
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("marker", "x")
          kwargs.setdefault("color", "red")
          kwargs.setdefault("label", "position")
  
          ax.scatter(lon, lat, transform=ax.get_transform(self.frame), **kwargs)
  
          # plot position error
          if not np.all(self.covariance.data == 0):
              region = self.position_error.to_pixel(ax.wcs)
              artist = region.as_artist(facecolor="none", edgecolor=kwargs["color"])
              ax.add_artist(artist)
  
          return ax
  
  
  class PointSpatialModel(SpatialModel):
      r"""Point Source.
  
      For more information see :ref:`point-spatial-model`.
  
      Parameters
      ----------
      lon_0, lat_0 : `~astropy.coordinates.Angle`
          Center position
      frame : {"icrs", "galactic"}
          Center position coordinate frame
      """
  
      tag = "PointSpatialModel"
      lon_0 = Parameter("lon_0", "0 deg")
      lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
  
      @property
      def evaluation_radius(self):
          """Evaluation radius (`~astropy.coordinates.Angle`).
  
          Set as zero degrees.
          """
          return 0 * u.deg
  
      @staticmethod
      def _grid_weights(x, y, x0, y0):
          """Compute 4-pixel weights such that centroid is preserved."""
          dx = np.abs(x - x0)
          dx = np.where(dx < 1, 1 - dx, 0)
  
         
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/spatial.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spectral models for Gammapy."""
  import operator
  import numpy as np
  import scipy.optimize
  import scipy.special
  import astropy.units as u
  from astropy import constants as const
  from astropy.table import Table
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.integrate import evaluate_integral_pwl, trapz_loglog
  from gammapy.utils.interpolation import ScaledRegularGridInterpolator
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  
  def integrate_spectrum(func, emin, emax, ndecade=100, intervals=False):
      """Integrate 1d function using the log-log trapezoidal rule.
  
      If scalar values for xmin and xmax are passed an oversampled grid is generated using the
      ``ndecade`` keyword argument. If xmin and xmax arrays are passed, no
      oversampling is performed and the integral is computed in the provided
      grid.
  
      Parameters
      ----------
      func : callable
          Function to integrate.
      emin : `~astropy.units.Quantity`
          Integration range minimum
      emax : `~astropy.units.Quantity`
          Integration range minimum
      ndecade : int, optional
          Number of grid points per decade used for the integration.
          Default : 100.
      intervals : bool, optional
          Return integrals in the grid not the sum, default: False
      """
      if emin.isscalar and emax.isscalar:
          energies = MapAxis.from_energy_bounds(
              emin=emin, emax=emax, nbin=ndecade, per_decade=True
          ).edges
      else:
          energies = edges_from_lo_hi(emin, emax)
  
      values = func(energies)
  
      integral = trapz_loglog(values, energies)
  
      if intervals:
          return integral
  
      return integral.sum()
  
  
  class SpectralModel(Model):
      """Spectral model base class."""
  
      def __call__(self, energy):
          kwargs = {par.name: par.quantity for par in self.parameters}
          kwargs = self._convert_evaluate_unit(kwargs, energy)
          return self.evaluate(energy, **kwargs)
  
      @staticmethod
      def _convert_evaluate_unit(kwargs_ref, energy):
          kwargs = {}
          for name, quantity in kwargs_ref.items():
              if quantity.unit.physical_type == "energy":
                  quantity = quantity.to(energy.unit)
              kwargs[name] = quantity
          return kwargs
  
      def __add__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.add)
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __sub__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.sub)
  
      def __rsub__(self, model):
          return self.__sub__(model)
  
      def _evaluate_gradient(self, energy, eps):
          n = len(self.parameters)
          f = self(energy)
          shape = (n, len(np.atleast_1d(energy)))
          df_dp = np.zeros(shape)
  
          for idx, parameter in enumerate(self.parameters):
              if parameter.frozen or eps[idx] == 0:
                  continue
  
              parameter.value += eps[idx]
              df = self(energy) - f
              df_dp[idx] = df.value / eps[idx]
  
              # Reset model to original parameter
              parameter.value -= eps[idx]
  
          return df_dp
  
      def evaluate_error(self, energy, epsilon=1e-4):
          """Evaluate spectral model with error propagation.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to evaluate
          epsilon : float
              Step size of the gradient evaluation. Given as a
              fraction of the parameter error.
  
          Returns
          -------
          dnde, dnde_error : tuple of `~astropy.units.Quantity`
              Tuple of flux and flux error.
          """
          p_cov = self.covariance
          eps = np.sqrt(np.diag(p_cov)) * epsilon
  
          df_dp = self._evaluate_gradient(energy, eps)
          f_cov = df_dp.T @ p_cov @ df_dp
          f_err = np.sqrt(np.diagonal(f_cov))
  
          q = self(energy)
          return u.Quantity([q.value, f_err], unit=q.unit)
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate spectral model numerically.
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE
  
          If array input for ``emin`` and ``emax`` is given you have to set
          ``intervals=True`` if you want the integral in each energy bin.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to :func:`~gammapy.utils.integrate.integrate_spectrum`
          """
          return integrate_spectrum(self, emin, emax, **kwargs)
  
      def energy_flux(self, emin, emax, **kwargs):
          r"""Compute energy flux in given energy range.
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          def f(x):
              return x * self(x)
  
          return integrate_spectrum(f, emin, emax, **kwargs)
  
      def plot(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model curve.
  
          kwargs are forwarded to `matplotlib.pyplot.plot`
  
          By default a log-log scaling of the axes is used, if you want to change
          the y axis scaling to linear you can use::
  
              from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel
              from astropy import units as u
  
              pwl = ExpCutoffPowerLawSpectralModel()
              ax = pwl.plot(energy_range=(0.1, 100) * u.TeV)
              ax.set_yscale('linear')
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          # evaluate model
          flux = self(energy).to(flux_unit)
  
          y = self._plot_scale_flux(energy, flux, energy_power)
  
          ax.plot(energy.value, y.value, **kwargs)
  
          self._plot_format_ax(ax, energy, y, energy_power)
          return ax
  
      def plot_error(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model error band.
  
          .. note::
  
              This method calls ``ax.set_yscale("log", nonposy='clip')`` and
              ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation.
              The additional argument ``nonposx='clip'`` avoids artefacts in the plot,
              when the error band extends to negative values (see also
              https://github.com/matplotlib/matplotlib/issues/8623
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/spectral.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spectral models for Gammapy."""
  import operator
  import numpy as np
  import scipy.optimize
  import scipy.special
  import astropy.units as u
  from astropy import constants as const
  from astropy.table import Table
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.integrate import evaluate_integral_pwl, trapz_loglog
  from gammapy.utils.interpolation import ScaledRegularGridInterpolator
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  
  def integrate_spectrum(func, emin, emax, ndecade=100, intervals=False):
      """Integrate 1d function using the log-log trapezoidal rule.
  
      If scalar values for xmin and xmax are passed an oversampled grid is generated using the
      ``ndecade`` keyword argument. If xmin and xmax arrays are passed, no
      oversampling is performed and the integral is computed in the provided
      grid.
  
      Parameters
      ----------
      func : callable
          Function to integrate.
      emin : `~astropy.units.Quantity`
          Integration range minimum
      emax : `~astropy.units.Quantity`
          Integration range minimum
      ndecade : int, optional
          Number of grid points per decade used for the integration.
          Default : 100.
      intervals : bool, optional
          Return integrals in the grid not the sum, default: False
      """
      if emin.isscalar and emax.isscalar:
          energies = MapAxis.from_energy_bounds(
              emin=emin, emax=emax, nbin=ndecade, per_decade=True
          ).edges
      else:
          energies = edges_from_lo_hi(emin, emax)
  
      values = func(energies)
  
      integral = trapz_loglog(values, energies)
  
      if intervals:
          return integral
  
      return integral.sum()
  
  
  class SpectralModel(Model):
      """Spectral model base class."""
  
      def __call__(self, energy):
          kwargs = {par.name: par.quantity for par in self.parameters}
          kwargs = self._convert_evaluate_unit(kwargs, energy)
          return self.evaluate(energy, **kwargs)
  
      @staticmethod
      def _convert_evaluate_unit(kwargs_ref, energy):
          kwargs = {}
          for name, quantity in kwargs_ref.items():
              if quantity.unit.physical_type == "energy":
                  quantity = quantity.to(energy.unit)
              kwargs[name] = quantity
          return kwargs
  
      def __add__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.add)
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __sub__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.sub)
  
      def __rsub__(self, model):
          return self.__sub__(model)
  
      def _evaluate_gradient(self, energy, eps):
          n = len(self.parameters)
          f = self(energy)
          shape = (n, len(np.atleast_1d(energy)))
          df_dp = np.zeros(shape)
  
          for idx, parameter in enumerate(self.parameters):
              if parameter.frozen or eps[idx] == 0:
                  continue
  
              parameter.value += eps[idx]
              df = self(energy) - f
              df_dp[idx] = df.value / eps[idx]
  
              # Reset model to original parameter
              parameter.value -= eps[idx]
  
          return df_dp
  
      def evaluate_error(self, energy, epsilon=1e-4):
          """Evaluate spectral model with error propagation.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to evaluate
          epsilon : float
              Step size of the gradient evaluation. Given as a
              fraction of the parameter error.
  
          Returns
          -------
          dnde, dnde_error : tuple of `~astropy.units.Quantity`
              Tuple of flux and flux error.
          """
          p_cov = self.covariance
          eps = np.sqrt(np.diag(p_cov)) * epsilon
  
          df_dp = self._evaluate_gradient(energy, eps)
          f_cov = df_dp.T @ p_cov @ df_dp
          f_err = np.sqrt(np.diagonal(f_cov))
  
          q = self(energy)
          return u.Quantity([q.value, f_err], unit=q.unit)
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate spectral model numerically.
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE
  
          If array input for ``emin`` and ``emax`` is given you have to set
          ``intervals=True`` if you want the integral in each energy bin.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to :func:`~gammapy.utils.integrate.integrate_spectrum`
          """
          return integrate_spectrum(self, emin, emax, **kwargs)
  
      def energy_flux(self, emin, emax, **kwargs):
          r"""Compute energy flux in given energy range.
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          def f(x):
              return x * self(x)
  
          return integrate_spectrum(f, emin, emax, **kwargs)
  
      def plot(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model curve.
  
          kwargs are forwarded to `matplotlib.pyplot.plot`
  
          By default a log-log scaling of the axes is used, if you want to change
          the y axis scaling to linear you can use::
  
              from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel
              from astropy import units as u
  
              pwl = ExpCutoffPowerLawSpectralModel()
              ax = pwl.plot(energy_range=(0.1, 100) * u.TeV)
              ax.set_yscale('linear')
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          # evaluate model
          flux = self(energy).to(flux_unit)
  
          y = self._plot_scale_flux(energy, flux, energy_power)
  
          ax.plot(energy.value, y.value, **kwargs)
  
          self._plot_format_ax(ax, energy, y, energy_power)
          return ax
  
      def plot_error(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model error band.
  
          .. note::
  
              This method calls ``ax.set_yscale("log", nonposy='clip')`` and
              ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation.
              The additional argument ``nonposx='clip'`` avoids artefacts in the plot,
              when the error band extends to negative values (see also
              https://github.com/matplotlib/matplotlib/issues/8623
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/spectral.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.modeling.models import Model, Parameter, Parameters
  from gammapy.datasets import Datasets
  from gammapy.utils.testing import requires_data
  
  
  class MyModel(Model):
      """Simple model example"""
  
      x = Parameter("x", 1, "cm")
      y = Parameter("y", 2)
  
  
  class CoModel(Model):
      """Compound model example"""
  
      norm = Parameter("norm", 42, "cm")
  
      def __init__(self, m1, m2, norm=norm.quantity):
          self.m1 = m1
          self.m2 = m2
          super().__init__(norm=norm)
  
      @property
      def parameters(self):
          return Parameters([self.norm]) + self.m1.parameters + self.m2.parameters
  
  
  class WrapperModel(Model):
      """Wrapper compound model.
  
      Dynamically generated parameters in `__init__`,
      and a parameter name conflict with the wrapped
      model, both have a parameter called "y".
      """
  
      def __init__(self, m1, a=1, y=99):
          self.m1 = m1
          a = Parameter("a", a)
          y = Parameter("y", y)
          self.default_parameters = Parameters([a, y])
          super().__init__(a=a, y=y)
  
      @property
      def parameters(self):
          return Parameters([self.a, self.y]) + self.m1.parameters
  
  
  def test_model_class():
      assert isinstance(MyModel.parameters, property)
      assert MyModel.x.name == "x"
      assert MyModel.default_parameters["x"] is MyModel.x
  
  
  def test_model_class_par_init():
      x = Parameter("x", 4, "cm")
      y = Parameter("y", 10)
  
      model = MyModel(x=x, y=y)
  
      assert x is model.x
      assert y is model.y
  
  
  def test_model_init():
      m = MyModel()
      assert m.x.name == "x"
      assert m.x.value == 1
      assert m.x is m.parameters[0]
      assert m.y is m.parameters[1]
      assert m.parameters is not MyModel.default_parameters
  
      m = MyModel(x="99 cm")
      assert m.x.value == 99
      assert m.y.value == 2
  
      # Currently we always convert to the default unit of a parameter
      # TODO: discuss if this is the behaviour we want, or if we instead
      # should change to the user-set unit, as long as it's compatible
      m = MyModel(x=99 * u.m)
      assert_allclose(m.x.value, 9900)
      assert m.x.unit == "cm"
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99)
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99 * u.s)
  
  
  def test_wrapper_model():
      outer = MyModel()
      m = WrapperModel(outer)
  
      assert isinstance(m.a, Parameter)
      assert m.y.value == 99
  
      assert m.parameters.names == ["a", "y", "x", "y"]
  
  
  def test_model_parameter():
      m = MyModel(x="99 cm")
      assert isinstance(m.x, Parameter)
      assert m.x.value == 99
      assert m.x.unit == "cm"
  
      with pytest.raises(TypeError):
          m.x = 99
  
      with pytest.raises(TypeError):
          m.x = 99 * u.cm
  
  
  # TODO: implement parameter linking. Not working ATM!
  @pytest.mark.xfail()
  def test_model_parameter_link():
      # Assigning a parameter should create a link
      m = MyModel()
      par = MyModel.x.copy()
      m.x = par
      assert isinstance(m.x, Parameter)
      assert m.x is par
      # model.parameters should be in sync with attributes
      assert m.x is m.parameters["x"]
  
  
  def test_model_copy():
      m = MyModel()
  
      m2 = m.copy()
  
      # Models should be independent
      assert m.parameters is not m2.parameters
      assert m.parameters[0] is not m2.parameters[0]
  
  
  def test_model_create():
      spectral_model = Model.create(
          "PowerLaw2SpectralModel", amplitude="1e-10 cm-2 s-1", index=3
      )
      assert "PowerLaw2SpectralModel" in spectral_model.tag
      assert_allclose(spectral_model.index.value, 3)
  
  
  def test_compound_model():
      m1 = MyModel()
      m2 = MyModel(x=10 * u.cm, y=20)
      m = CoModel(m1, m2)
      assert len(m.parameters) == 5
      assert m.parameters.names == ["norm", "x", "y", "x", "y"]
      assert_allclose(m.parameters.values, [42, 1, 2, 10, 20])
  
  
  def test_parameter_link_init():
      m1 = MyModel()
      m2 = MyModel(y=m1.y)
  
      assert m1.y is m2.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  def test_parameter_link():
      m1 = MyModel()
      m2 = MyModel()
  
      m2.y = m1.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  @requires_data()
  def test_models_management(tmp_path):
      path = "$GAMMAPY_DATA/tests/models"
      filedata = "gc_example_datasets.yaml"
      filemodel = "gc_example_models.yaml"
  
      datasets = Datasets.read(path, filedata, filemodel)
  
      model1 = datasets.models[0].copy(name="model1", datasets_names=None)
      model2 = datasets.models[0].copy(name="model2", datasets_names=[datasets[1].name])
      model3 = datasets.models[0].copy(name="model3", datasets_names=[datasets[0].name])
  
      model1b = datasets.models[0].copy(name="model1", datasets_names=None)
      model1b.spectral_model.amplitude.value *= 2
  
      names0 = datasets[0].models.names
      names1 = datasets[1].models.names
  
      datasets[0].models.append(model1)
      _ = datasets[0].models + model2
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[0].models["model1"].datasets_names == None
      assert datasets[0].models["model2"].datasets_names == [
          datasets[1].name,
          datasets[0].name,
      ]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      # TODO consistency check at datasets level ?
      # or force same Models for each dataset._models on datasets init ?
      # here we have the right behavior: model1 and model2 are also added to dataset1
      # because serialization create a global model object shared by all datasets
      # if that was not the case we could have inconsistancies
      # such as model1.datasets_names == None added only to dataset1
      # user can still create such inconsistancies if they define datasets
      # with diferent Models objects for each dataset.
  
      del datasets[0].models["model1"]
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.remove(model2)
      assert datasets[0].models.names == names0
  
      datasets.models.append(model2)
      assert model2 in datasets.models
      assert model2 in datasets[1].models
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.extend([model1, model3])
      assert datasets[0].models.names == names0 + ["model2", "model1", "model3"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      assert datasets[0].models.names == names0
      assert datasets[1].models.names == names1
      datasets.models.extend([model1, model2, model3])
      assert datasets[0].models.names == names0 + ["model1", "model2", "model3"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      _ = datasets.models + [model1, model2]
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      datasets[0].models["model2"] = model3
      assert datasets[0].models.names == names0 + ["model1", "model3"]
      assert datasets[1].models
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_core.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.modeling.models import Model, Parameter, Parameters
  from gammapy.datasets import Datasets
  from gammapy.utils.testing import requires_data
  
  
  class MyModel(Model):
      """Simple model example"""
  
      x = Parameter("x", 1, "cm")
      y = Parameter("y", 2)
  
  
  class CoModel(Model):
      """Compound model example"""
  
      norm = Parameter("norm", 42, "cm")
  
      def __init__(self, m1, m2, norm=norm.quantity):
          self.m1 = m1
          self.m2 = m2
          super().__init__(norm=norm)
  
      @property
      def parameters(self):
          return Parameters([self.norm]) + self.m1.parameters + self.m2.parameters
  
  
  class WrapperModel(Model):
      """Wrapper compound model.
  
      Dynamically generated parameters in `__init__`,
      and a parameter name conflict with the wrapped
      model, both have a parameter called "y".
      """
  
      def __init__(self, m1, a=1, y=99):
          self.m1 = m1
          a = Parameter("a", a)
          y = Parameter("y", y)
          self.default_parameters = Parameters([a, y])
          super().__init__(a=a, y=y)
  
      @property
      def parameters(self):
          return Parameters([self.a, self.y]) + self.m1.parameters
  
  
  def test_model_class():
      assert isinstance(MyModel.parameters, property)
      assert MyModel.x.name == "x"
      assert MyModel.default_parameters["x"] is MyModel.x
  
  
  def test_model_class_par_init():
      x = Parameter("x", 4, "cm")
      y = Parameter("y", 10)
  
      model = MyModel(x=x, y=y)
  
      assert x is model.x
      assert y is model.y
  
  
  def test_model_init():
      m = MyModel()
      assert m.x.name == "x"
      assert m.x.value == 1
      assert m.x is m.parameters[0]
      assert m.y is m.parameters[1]
      assert m.parameters is not MyModel.default_parameters
  
      m = MyModel(x="99 cm")
      assert m.x.value == 99
      assert m.y.value == 2
  
      # Currently we always convert to the default unit of a parameter
      # TODO: discuss if this is the behaviour we want, or if we instead
      # should change to the user-set unit, as long as it's compatible
      m = MyModel(x=99 * u.m)
      assert_allclose(m.x.value, 9900)
      assert m.x.unit == "cm"
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99)
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99 * u.s)
  
  
  def test_wrapper_model():
      outer = MyModel()
      m = WrapperModel(outer)
  
      assert isinstance(m.a, Parameter)
      assert m.y.value == 99
  
      assert m.parameters.names == ["a", "y", "x", "y"]
  
  
  def test_model_parameter():
      m = MyModel(x="99 cm")
      assert isinstance(m.x, Parameter)
      assert m.x.value == 99
      assert m.x.unit == "cm"
  
      with pytest.raises(TypeError):
          m.x = 99
  
      with pytest.raises(TypeError):
          m.x = 99 * u.cm
  
  
  # TODO: implement parameter linking. Not working ATM!
  @pytest.mark.xfail()
  def test_model_parameter_link():
      # Assigning a parameter should create a link
      m = MyModel()
      par = MyModel.x.copy()
      m.x = par
      assert isinstance(m.x, Parameter)
      assert m.x is par
      # model.parameters should be in sync with attributes
      assert m.x is m.parameters["x"]
  
  
  def test_model_copy():
      m = MyModel()
  
      m2 = m.copy()
  
      # Models should be independent
      assert m.parameters is not m2.parameters
      assert m.parameters[0] is not m2.parameters[0]
  
  
  def test_model_create():
      spectral_model = Model.create(
          "PowerLaw2SpectralModel", amplitude="1e-10 cm-2 s-1", index=3
      )
      assert spectral_model.tag == "PowerLaw2SpectralModel"
      assert_allclose(spectral_model.index.value, 3)
  
  
  def test_compound_model():
      m1 = MyModel()
      m2 = MyModel(x=10 * u.cm, y=20)
      m = CoModel(m1, m2)
      assert len(m.parameters) == 5
      assert m.parameters.names == ["norm", "x", "y", "x", "y"]
      assert_allclose(m.parameters.values, [42, 1, 2, 10, 20])
  
  
  def test_parameter_link_init():
      m1 = MyModel()
      m2 = MyModel(y=m1.y)
  
      assert m1.y is m2.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  def test_parameter_link():
      m1 = MyModel()
      m2 = MyModel()
  
      m2.y = m1.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  @requires_data()
  def test_models_management(tmp_path):
      path = "$GAMMAPY_DATA/tests/models"
      filedata = "gc_example_datasets.yaml"
      filemodel = "gc_example_models.yaml"
  
      datasets = Datasets.read(path, filedata, filemodel)
  
      model1 = datasets.models[0].copy(name="model1", datasets_names=None)
      model2 = datasets.models[0].copy(name="model2", datasets_names=[datasets[1].name])
      model3 = datasets.models[0].copy(name="model3", datasets_names=[datasets[0].name])
  
      model1b = datasets.models[0].copy(name="model1", datasets_names=None)
      model1b.spectral_model.amplitude.value *= 2
  
      names0 = datasets[0].models.names
      names1 = datasets[1].models.names
  
      datasets[0].models.append(model1)
      _ = datasets[0].models + model2
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[0].models["model1"].datasets_names == None
      assert datasets[0].models["model2"].datasets_names == [
          datasets[1].name,
          datasets[0].name,
      ]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      # TODO consistency check at datasets level ?
      # or force same Models for each dataset._models on datasets init ?
      # here we have the right behavior: model1 and model2 are also added to dataset1
      # because serialization create a global model object shared by all datasets
      # if that was not the case we could have inconsistancies
      # such as model1.datasets_names == None added only to dataset1
      # user can still create such inconsistancies if they define datasets
      # with diferent Models objects for each dataset.
  
      del datasets[0].models["model1"]
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.remove(model2)
      assert datasets[0].models.names == names0
  
      datasets.models.append(model2)
      assert model2 in datasets.models
      assert model2 in datasets[1].models
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.extend([model1, model3])
      assert datasets[0].models.names == names0 + ["model2", "model1", "model3"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      assert datasets[0].models.names == names0
      assert datasets[1].models.names == names1
      datasets.models.extend([model1, model2, model3])
      assert datasets[0].models.names == names0 + ["model1", "model2", "model3"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      _ = datasets.models + [model1, model2]
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      datasets[0].models["model2"] = model3
      assert datasets[0].models.names == names0 + ["model1", "model3"]
      assert datasets[1].models.n
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_core.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.table import Table
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling.models import (
      MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      BackgroundModel,
      Model,
      Models,
  )
  from gammapy.utils.scripts import read_yaml, write_yaml
  from gammapy.utils.testing import requires_data
  
  
  @requires_data()
  def test_dict_to_skymodels():
      filename = get_pkg_data_filename("data/examples.yaml")
      models_data = read_yaml(filename)
      models = Models.from_dict(models_data)
  
      assert len(models) == 5
  
      model0 = models[0]
      assert isinstance(model0, BackgroundModel)
      assert model0.name == "background_irf"
  
      model0 = models[1]
      assert "ExpCutoffPowerLawSpectralModel" in model0.spectral_model.tag
      assert "PointSpatialModel" in model0.spatial_model.tag
  
      pars0 = model0.parameters
      assert pars0["index"].value == 2.1
      assert pars0["index"].unit == ""
      assert np.isnan(pars0["index"].max)
      assert np.isnan(pars0["index"].min)
      assert not pars0["index"].frozen
  
      assert pars0["lon_0"].value == -0.5
      assert pars0["lon_0"].unit == "deg"
      assert pars0["lon_0"].max == 180.0
      assert pars0["lon_0"].min == -180.0
      assert pars0["lon_0"].frozen
  
      assert pars0["lat_0"].value == -0.0005
      assert pars0["lat_0"].unit == "deg"
      assert pars0["lat_0"].max == 90.0
      assert pars0["lat_0"].min == -90.0
      assert pars0["lat_0"].frozen
  
      assert pars0["lambda_"].value == 0.006
      assert pars0["lambda_"].unit == "TeV-1"
      assert np.isnan(pars0["lambda_"].min)
      assert np.isnan(pars0["lambda_"].max)
  
      model1 = models[2]
      assert "PL" in model1.spectral_model.tag
      assert "PowerLawSpectralModel" in model1.spectral_model.tag
      assert "DiskSpatialModel" in model1.spatial_model.tag
      assert "disk" in model1.spatial_model.tag
      assert "LightCurveTemplateTemporalModel" in model1.temporal_model.tag
  
      pars1 = model1.parameters
      assert pars1["index"].value == 2.2
      assert pars1["index"].unit == ""
      assert pars1["lat_0"].scale == 1.0
      assert pars1["lat_0"].factor == pars1["lat_0"].value
  
      assert np.isnan(pars1["index"].max)
      assert np.isnan(pars1["index"].min)
  
      assert pars1["r_0"].unit == "deg"
  
      model2 = models[3]
      assert_allclose(model2.spectral_model.energy.data, [34.171, 44.333, 57.517])
      assert model2.spectral_model.energy.unit == "MeV"
      assert_allclose(
          model2.spectral_model.values.data, [2.52894e-06, 1.2486e-06, 6.14648e-06]
      )
      assert model2.spectral_model.values.unit == "1 / (cm2 MeV s sr)"
  
      assert "TemplateSpectralModel" in model2.spectral_model.tag
      assert "TemplateSpatialModel" in model2.spatial_model.tag
  
      assert model2.spatial_model.parameters["norm"].value == 1.0
      assert not model2.spatial_model.normalize
      assert model2.spectral_model.parameters["norm"].value == 2.1
  
      # TODO problem of duplicate parameter name between TemplateSpatialModel and TemplateSpectralModel
      # assert model2.parameters["norm"].value == 2.1 # fail
  
  
  @requires_data()
  def test_sky_models_io(tmp_path):
      # TODO: maybe change to a test case where we create a model programatically?
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      models.covariance = np.eye(len(models.parameters))
      models.write(tmp_path / "tmp.yaml")
      models = Models.read(tmp_path / "tmp.yaml")
      assert models._covar_file == "tmp_covariance.dat"
      assert_allclose(models.covariance.data, np.eye(len(models.parameters)))
      assert_allclose(models.parameters["lat_0"].min, -90.0)
  
      # TODO: not sure if we should just round-trip, or if we should
      # check YAML file content (e.g. against a ref file in the repo)
      # or check serialised dict content
  
  
  @requires_data()
  def test_absorption_io(tmp_path):
      dominguez = Absorption.read_builtin("dominguez")
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel"),
          absorption=dominguez,
          redshift=0.5,
      )
      assert len(model.parameters) == 5
  
      model_dict = model.to_dict()
      parnames = [_["name"] for _ in model_dict["parameters"]]
      assert parnames == ["redshift", "alpha_norm"]
  
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert new_model.redshift.value == 0.5
      assert new_model.alpha_norm.name == "alpha_norm"
      assert new_model.alpha_norm.value == 1
      assert "PowerLawSpectralModel" in new_model.spectral_model.tag
      assert_allclose(new_model.absorption.energy, dominguez.energy)
      assert_allclose(new_model.absorption.param, dominguez.param)
      assert len(new_model.parameters) == 5
  
      test_absorption = Absorption(
          u.Quantity(range(3), "keV"),
          u.Quantity(range(2), ""),
          u.Quantity(np.ones((2, 3)), ""),
      )
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel"),
          absorption=test_absorption,
          redshift=0.5,
      )
      model_dict = model.to_dict()
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert_allclose(new_model.absorption.energy, test_absorption.energy)
      assert_allclose(new_model.absorption.param, test_absorption.param)
  
      write_yaml(model_dict, tmp_path / "tmp.yaml")
      read_yaml(tmp_path / "tmp.yaml")
  
  
  def make_all_models():
      """Make an instance of each model, for testing."""
      yield Model.create("ConstantSpatialModel")
      map_constantmodel = Map.create(npix=(10, 20), unit="sr-1")
      yield Model.create("TemplateSpatialModel", map=map_constantmodel)
      yield Model.create("DiskSpatialModel", lon_0="1 deg", lat_0="2 deg", r_0="3 deg")
      yield Model.create(
          "GaussianSpatialModel", lon_0="1 deg", lat_0="2 deg", sigma="3 deg"
      )
      yield Model.create("PointSpatialModel", lon_0="1 deg", lat_0="2 deg")
      yield Model.create(
          "ShellSpatialModel", lon_0="1 deg", lat_0="2 deg", radius="3 deg", width="4 deg"
      )
      yield Model.create("ConstantSpectralModel", const="99 cm-2 s-1 TeV-1")
      # TODO: yield Model.create("CompoundSpectralModel")
      yield Model.create("PowerLawSpectralModel")
      yield Model.create("PowerLaw2SpectralModel")
      yield Model.create("ExpCutoffPowerLawSpectralModel")
      yield Model.create("ExpCutoffPowerLaw3FGLSpectralModel")
      yield Model.create("SuperExpCutoffPowerLaw3FGLSpectralModel")
      yield Model.create("SuperExpCutoffPowerLaw4FGLSpectralModel")
      yield Model.create("LogParabolaSpectralModel")
      yield Model.create(
          "TemplateSpectralModel", energy=[1, 2] * u.cm, values=[3, 4] * u.cm
      )  # TODO: add unit validation?
      yield Model.create("GaussianSpectralModel")
      # TODO: yield Model.create("AbsorbedSpectralModel")
     
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_io.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.table import Table
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling.models import (
      MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      BackgroundModel,
      Model,
      Models,
  )
  from gammapy.utils.scripts import read_yaml, write_yaml
  from gammapy.utils.testing import requires_data
  
  
  @requires_data()
  def test_dict_to_skymodels():
      filename = get_pkg_data_filename("data/examples.yaml")
      models_data = read_yaml(filename)
      models = Models.from_dict(models_data)
  
      assert len(models) == 5
  
      model0 = models[0]
      assert isinstance(model0, BackgroundModel)
      assert model0.name == "background_irf"
  
      model0 = models[1]
      assert model0.spectral_model.tag == "ExpCutoffPowerLawSpectralModel"
      assert model0.spatial_model.tag == "PointSpatialModel"
  
      pars0 = model0.parameters
      assert pars0["index"].value == 2.1
      assert pars0["index"].unit == ""
      assert np.isnan(pars0["index"].max)
      assert np.isnan(pars0["index"].min)
      assert not pars0["index"].frozen
  
      assert pars0["lon_0"].value == -0.5
      assert pars0["lon_0"].unit == "deg"
      assert pars0["lon_0"].max == 180.0
      assert pars0["lon_0"].min == -180.0
      assert pars0["lon_0"].frozen
  
      assert pars0["lat_0"].value == -0.0005
      assert pars0["lat_0"].unit == "deg"
      assert pars0["lat_0"].max == 90.0
      assert pars0["lat_0"].min == -90.0
      assert pars0["lat_0"].frozen
  
      assert pars0["lambda_"].value == 0.006
      assert pars0["lambda_"].unit == "TeV-1"
      assert np.isnan(pars0["lambda_"].min)
      assert np.isnan(pars0["lambda_"].max)
  
      model1 = models[2]
      assert model1.spectral_model.tag == "PowerLawSpectralModel"
      assert model1.spatial_model.tag == "DiskSpatialModel"
      assert model1.temporal_model.tag == "LightCurveTemplateTemporalModel"
  
      pars1 = model1.parameters
      assert pars1["index"].value == 2.2
      assert pars1["index"].unit == ""
      assert pars1["lat_0"].scale == 1.0
      assert pars1["lat_0"].factor == pars1["lat_0"].value
  
      assert np.isnan(pars1["index"].max)
      assert np.isnan(pars1["index"].min)
  
      assert pars1["r_0"].unit == "deg"
  
      model2 = models[3]
      assert_allclose(model2.spectral_model.energy.data, [34.171, 44.333, 57.517])
      assert model2.spectral_model.energy.unit == "MeV"
      assert_allclose(
          model2.spectral_model.values.data, [2.52894e-06, 1.2486e-06, 6.14648e-06]
      )
      assert model2.spectral_model.values.unit == "1 / (cm2 MeV s sr)"
  
      assert model2.spectral_model.tag == "TemplateSpectralModel"
      assert model2.spatial_model.tag == "TemplateSpatialModel"
  
      assert model2.spatial_model.parameters["norm"].value == 1.0
      assert not model2.spatial_model.normalize
      assert model2.spectral_model.parameters["norm"].value == 2.1
  
      # TODO problem of duplicate parameter name between TemplateSpatialModel and TemplateSpectralModel
      # assert model2.parameters["norm"].value == 2.1 # fail
  
  
  @requires_data()
  def test_sky_models_io(tmp_path):
      # TODO: maybe change to a test case where we create a model programatically?
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      models.covariance = np.eye(len(models.parameters))
      models.write(tmp_path / "tmp.yaml")
      models = Models.read(tmp_path / "tmp.yaml")
      assert models._covar_file == "tmp_covariance.dat"
      assert_allclose(models.covariance.data, np.eye(len(models.parameters)))
      assert_allclose(models.parameters["lat_0"].min, -90.0)
  
      # TODO: not sure if we should just round-trip, or if we should
      # check YAML file content (e.g. against a ref file in the repo)
      # or check serialised dict content
  
  
  @requires_data()
  def test_absorption_io(tmp_path):
      dominguez = Absorption.read_builtin("dominguez")
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel"),
          absorption=dominguez,
          redshift=0.5,
      )
      assert len(model.parameters) == 5
  
      model_dict = model.to_dict()
      parnames = [_["name"] for _ in model_dict["parameters"]]
      assert parnames == ["redshift", "alpha_norm"]
  
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert new_model.redshift.value == 0.5
      assert new_model.alpha_norm.name == "alpha_norm"
      assert new_model.alpha_norm.value == 1
      assert new_model.spectral_model.tag == "PowerLawSpectralModel"
      assert_allclose(new_model.absorption.energy, dominguez.energy)
      assert_allclose(new_model.absorption.param, dominguez.param)
      assert len(new_model.parameters) == 5
  
      test_absorption = Absorption(
          u.Quantity(range(3), "keV"),
          u.Quantity(range(2), ""),
          u.Quantity(np.ones((2, 3)), ""),
      )
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel"),
          absorption=test_absorption,
          redshift=0.5,
      )
      model_dict = model.to_dict()
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert_allclose(new_model.absorption.energy, test_absorption.energy)
      assert_allclose(new_model.absorption.param, test_absorption.param)
  
      write_yaml(model_dict, tmp_path / "tmp.yaml")
      read_yaml(tmp_path / "tmp.yaml")
  
  
  def make_all_models():
      """Make an instance of each model, for testing."""
      yield Model.create("ConstantSpatialModel")
      map_constantmodel = Map.create(npix=(10, 20), unit="sr-1")
      yield Model.create("TemplateSpatialModel", map=map_constantmodel)
      yield Model.create("DiskSpatialModel", lon_0="1 deg", lat_0="2 deg", r_0="3 deg")
      yield Model.create(
          "GaussianSpatialModel", lon_0="1 deg", lat_0="2 deg", sigma="3 deg"
      )
      yield Model.create("PointSpatialModel", lon_0="1 deg", lat_0="2 deg")
      yield Model.create(
          "ShellSpatialModel", lon_0="1 deg", lat_0="2 deg", radius="3 deg", width="4 deg"
      )
      yield Model.create("ConstantSpectralModel", const="99 cm-2 s-1 TeV-1")
      # TODO: yield Model.create("CompoundSpectralModel")
      yield Model.create("PowerLawSpectralModel")
      yield Model.create("PowerLaw2SpectralModel")
      yield Model.create("ExpCutoffPowerLawSpectralModel")
      yield Model.create("ExpCutoffPowerLaw3FGLSpectralModel")
      yield Model.create("SuperExpCutoffPowerLaw3FGLSpectralModel")
      yield Model.create("SuperExpCutoffPowerLaw4FGLSpectralModel")
      yield Model.create("LogParabolaSpectralModel")
      yield Model.create(
          "TemplateSpectralModel", energy=[1, 2] * u.cm, values=[3, 4] * u.cm
      )  # TODO: add unit validation?
      yield Model.create("GaussianSpectralModel")
      # TODO: yield Model.create("AbsorbedSpectralModel")
      # TODO: yield Model.create("NaimaSpectralModel")
      # TODO: yield Model.create("ScaleSpectralModel")
      yield Model.create("ConstantTemporalModel
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_io.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  __all__ = ["Registry"]
  
  
  class Registry(list):
      """Registry class."""
  
      def get_cls(self, tag):
          for cls in self:
              if hasattr(cls, "tag") and tag in cls.tag:
                  return cls
          raise KeyError(f"No model found with tag: {tag!r}")
  
      def __str__(self):
          info = "Registry\n"
          info += "--------\n\n"
          tags = [_.tag[0] if isinstance(_.tag, list) else _.tag for _ in self]
          len_max = max([len(tag) for tag in tags])
  
          for item in self:
              info += f"\t{item.tag:{len_max}s}: {item.__name__}\n"
  
          return info.expandtabs(tabsize=2)
  
  
  [... the rest is omitted ...]
files: gammapy/utils/registry.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  
  __all__ = ["Registry"]
  
  
  class Registry(list):
      """Registry class."""
  
      def get_cls(self, tag):
          for cls in self:
              if hasattr(cls, "tag") and cls.tag == tag:
                  return cls
          raise KeyError(f"No model found with tag: {tag!r}")
  
      def __str__(self):
          info = "Registry\n"
          info += "--------\n\n"
  
          len_max = max([len(_.tag) for _ in self])
  
          for item in self:
              info += f"\t{item.tag:{len_max}s}: {item.__name__}\n"
  
          return info.expandtabs(tabsize=2)
  
  
  [... the rest is omitted ...]
files: gammapy/utils/registry.py