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])
              error = np.nan_to_num(par.error ** 2, nan=1)
              covar = Covariance(pars, data=[[error]])
              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, full_output=True):
          """Create dict for YAML serialisation"""
          tag = self.tag[0] if isinstance(self.tag, list) else self.tag
          params = self.parameters.to_dict()
  
          if full_output is False:
              base = self.__class__
              names = self.parameters.names
              for k, name in enumerate(names):
                  init = base.__dict__[name].to_dict()
                  for item in ["min", "max", "frozen
  
  [... 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])
              error = np.nan_to_num(par.error ** 2, nan=1)
              covar = Covariance(pars, data=[[error]])
              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
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/core.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Cube models (axes: lon, lat, energy)."""
  import copy
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.modeling.parameter import _get_parameters_str
  from gammapy.utils.scripts import make_name, make_path
  from gammapy.utils.fits import LazyFitsData, HDULocation
  from .core import Model, Models
  from .spatial import SpatialModel
  from .spectral import SpectralModel, PowerLawNormSpectralModel
  from .temporal import TemporalModel
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([self, *other])
          elif isinstance(other, (SkyModel, BackgroundModel)):
              return Models([self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __call__(self, lon, lat, energy, time=None):
          return self.evaluate(lon, lat, energy, time)
  
      def evaluate_geom(self, geom, gti=None):
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat, coords["energy_true"])
  
  
  class SkyModel(SkyModelBase):
      """Sky model component.
  
      This model represents a factorised sky model.
      It has `~gammapy.modeling.Parameters`
      combining the spatial and spectral parameters.
  
      Parameters
      ----------
      spectral_model : `~gammapy.modeling.models.SpectralModel`
          Spectral model
      spatial_model : `~gammapy.modeling.models.SpatialModel`
          Spatial model (must be normalised to integrate to 1)
      temporal_model : `~gammapy.modeling.models.temporalModel`
          Temporal model
      name : str
          Model identifier
      apply_irf : dict
          Dictionary declaring which IRFs should be applied to this model. Options
          are {"exposure": True, "psf": True, "edisp": True}
      datasets_names : list of str
          Which datasets this model is applied to.
      """
  
      tag = "SkyModel"
      _apply_irf_default = {"exposure": True, "psf": True, "edisp": True}
  
      def __init__(
          self,
          spectral_model,
          spatial_model=None,
          temporal_model=None,
          name=None,
          apply_irf=None,
          datasets_names=None,
      ):
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          self.temporal_model = temporal_model
          self._name = make_name(name)
  
          if apply_irf is None:
              apply_irf = self._apply_irf_default.copy()
  
          self.apply_irf = apply_irf
          self.datasets_names = datasets_names
          self._check_unit()
          super().__init__()
  
      @property
      def _models(self):
          models = self.spectral_model, self.spatial_model, self.temporal_model
          return [model for model in models if model is not None]
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/cube.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Cube models (axes: lon, lat, energy)."""
  import copy
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.modeling.parameter import _get_parameters_str
  from gammapy.utils.scripts import make_name, make_path
  from gammapy.utils.fits import LazyFitsData, HDULocation
  from .core import Model, Models
  from .spatial import SpatialModel
  from .spectral import SpectralModel, PowerLawNormSpectralModel
  from .temporal import TemporalModel
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([self, *other])
          elif isinstance(other, (SkyModel, BackgroundModel)):
              return Models([self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __call__(self, lon, lat, energy, time=None):
          return self.evaluate(lon, lat, energy, time)
  
      def evaluate_geom(self, geom, gti=None):
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat, coords["energy_true"])
  
  
  class SkyModel(SkyModelBase):
      """Sky model component.
  
      This model represents a factorised sky model.
      It has `~gammapy.modeling.Parameters`
      combining the spatial and spectral parameters.
  
      Parameters
      ----------
      spectral_model : `~gammapy.modeling.models.SpectralModel`
          Spectral model
      spatial_model : `~gammapy.modeling.models.SpatialModel`
          Spatial model (must be normalised to integrate to 1)
      temporal_model : `~gammapy.modeling.models.temporalModel`
          Temporal model
      name : str
          Model identifier
      apply_irf : dict
          Dictionary declaring which IRFs should be applied to this model. Options
          are {"exposure": True, "psf": True, "edisp": True}
      datasets_names : list of str
          Which datasets this model is applied to.
      """
  
      tag = "SkyModel"
      _apply_irf_default = {"exposure": True, "psf": True, "edisp": True}
  
      def __init__(
          self,
          spectral_model,
          spatial_model=None,
          temporal_model=None,
          name=None,
          apply_irf=None,
          datasets_names=None,
      ):
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          self.temporal_model = temporal_model
          self._name = make_name(name)
  
          if apply_irf is None:
              apply_irf = self._apply_irf_default.copy()
  
          self.apply_irf = apply_irf
          self.datasets_names = datasets_names
          self._check_unit()
          super().__init__()
  
      @property
      def _models(self):
          models = self.spectral_model, self.spatial_model, self.temporal_model
          return [model for model in models if model is not None]
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/cube.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.utils import lazyproperty
  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
  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."""
  
      _type = "spatial"
  
      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, energy=None):
          """Call evaluate method"""
          kwargs = {par.name: par.quantity for par in self.parameters}
  
          if energy is None and self.is_energy_dependent:
              raise ValueError("Missing energy value for evaluation")
  
          if energy is not None:
              kwargs["energy"] = energy
  
          return self.evaluate(lon, lat, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      # TODO: make this a hard-coded class attribute?
      @lazyproperty
      def is_energy_dependent(self):
          varnames = self.evaluate.__code__.co_varnames
          return "energy" in varnames
  
      @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.cov
  
  [... 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.utils import lazyproperty
  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
  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."""
      _type = "spatial"
  
      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, energy=None):
          """Call evaluate method"""
          kwargs = {par.name: par.quantity for par in self.parameters}
  
          if energy is None and self.is_energy_dependent:
              raise ValueError("Missing energy value for evaluation")
  
          if energy is not None:
              kwargs["energy"] = energy
  
          return self.evaluate(lon, lat, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      # TODO: make this a hard-coded class attribute?
      @lazyproperty
      def is_energy_dependent(self):
          varnames = self.evaluate.__code__.co_varnames
          return "energy" in varnames
  
      @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
  
  [... 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 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."""
  
      _type = "spectral"
  
      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)
  
      @property
      def type(self):
          return self._type
  
      @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 __mul__(self, other):
          if isinstance(other, SpectralModel):
              return CompoundSpectralModel(self, other, operator.mul)
          else:
              raise TypeError(f"Multiplication invalid for type {other!r}")
  
      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
  
  [... 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 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."""
      _type = "spectral"
  
      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)
  
      @property
      def type(self):
          return self._type
  
      @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 __mul__(self, other):
          if isinstance(other, SpectralModel):
              return CompoundSpectralModel(self, other, operator.mul)
          else:
              raise TypeError(f"Multiplication invalid for type {other!r}")
  
      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
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/spectral.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Time-dependent models."""
  import numpy as np
  import scipy.interpolate
  from astropy import units as u
  from astropy.table import Table
  from astropy.time import Time
  from astropy.utils import lazyproperty
  from gammapy.modeling import Parameter
  from gammapy.utils.random import InverseCDFSampler, get_random_state
  from gammapy.utils.scripts import make_path
  from gammapy.utils.time import time_ref_from_dict
  from .core import Model
  
  
  # TODO: make this a small ABC to define a uniform interface.
  class TemporalModel(Model):
      """Temporal model base class.
      evaluates on  astropy.time.Time objects"""
  
      _type = "temporal"
  
      def __call__(self, time):
          """Evaluate model
  
          Parameters
          ----------
          time : `~astropy.time.Time`
              Time object
          """
          kwargs = {par.name: par.quantity for par in self.parameters}
          time = u.Quantity(time.mjd, "day")
          return self.evaluate(time, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def time_sum(t_min, t_max):
          """
          Total time between t_min and t_max
  
          Parameters
          ----------
          t_min, t_max: `~astropy.time.Time`
              Lower and upper bound of integration range
  
          Returns
          -------
          time_sum : `~astropy.time.TimeDelta`
              Summed time in the intervals.
  
          """
          diff = t_max - t_min
          # TODO: this is a work-around for https://github.com/astropy/astropy/issues/10501
          return u.Quantity(np.sum(diff.to_value("day")), "day")
  
      def plot(self, time_range, ax=None):
          """
          Plot Temporal Model.
  
          Parameters
          ----------
          time_range : `~astropy.time.Time`
              times to plot the model
          ax : `~matplotlib.axes.Axes`, optional
              axis
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              axis
          """
  
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          t_min, t_max = time_range
          n_value = 100
          delta = t_max - t_min
          times = t_min + delta * np.linspace(0, 1, n_value)
          val = self(times)
          ax.plot(times.mjd, val)
          return ax
  
  
  class ConstantTemporalModel(TemporalModel):
      """Constant temporal model."""
  
      tag = ["ConstantTemporalModel", "const"]
  
      @staticmethod
      def evaluate(time):
          """Evaluate at given times."""
          return np.ones(time.shape)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : `~astropy.units.Quantity`
              Integrated flux norm on the given time intervals
          """
          return (t_max - t_min) / self.time_sum(t_min, t_max)
  
      @staticmethod
      def sample_time(n_events, t_min, t_max, random_state=0):
          """Sample arrival times of events.
  
          Parameters
          ----------
          n_events : int
              Number of events to sample.
          t_min : `~astropy.time.Time`
              Start time of the sampling
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/temporal.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Time-dependent models."""
  import numpy as np
  import scipy.interpolate
  from astropy import units as u
  from astropy.table import Table
  from astropy.time import Time
  from astropy.utils import lazyproperty
  from gammapy.modeling import Parameter
  from gammapy.utils.random import InverseCDFSampler, get_random_state
  from gammapy.utils.scripts import make_path
  from gammapy.utils.time import time_ref_from_dict
  from .core import Model
  
  
  # TODO: make this a small ABC to define a uniform interface.
  class TemporalModel(Model):
      """Temporal model base class.
      evaluates on  astropy.time.Time objects"""
      _type = "temporal"
  
      def __call__(self, time):
          """Evaluate model
  
          Parameters
          ----------
          time : `~astropy.time.Time`
              Time object
          """
          kwargs = {par.name: par.quantity for par in self.parameters}
          time = u.Quantity(time.mjd, "day")
          return self.evaluate(time, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def time_sum(t_min, t_max):
          """
          Total time between t_min and t_max
  
          Parameters
          ----------
          t_min, t_max: `~astropy.time.Time`
              Lower and upper bound of integration range
  
          Returns
          -------
          time_sum : `~astropy.time.TimeDelta`
              Summed time in the intervals.
  
          """
          diff = t_max - t_min
          # TODO: this is a work-around for https://github.com/astropy/astropy/issues/10501
          return u.Quantity(np.sum(diff.to_value("day")), "day")
  
      def plot(self, time_range, ax=None):
          """
          Plot Temporal Model.
  
          Parameters
          ----------
          time_range : `~astropy.time.Time`
              times to plot the model
          ax : `~matplotlib.axes.Axes`, optional
              axis
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              axis
          """
  
  
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          t_min, t_max = time_range
          n_value = 100
          delta = (t_max - t_min)
          times = t_min + delta * np.linspace(0, 1, n_value)
          val = self(times)
          ax.plot(times.mjd, val)
          return ax
  
  class ConstantTemporalModel(TemporalModel):
      """Constant temporal model."""
  
      tag = ["ConstantTemporalModel", "const"]
  
      @staticmethod
      def evaluate(time):
          """Evaluate at given times."""
          return np.ones(time.shape)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : `~astropy.units.Quantity`
              Integrated flux norm on the given time intervals
          """
          return (t_max - t_min) / self.time_sum(t_min, t_max)
  
      @staticmethod
      def sample_time(n_events, t_min, t_max, random_state=0):
          """Sample arrival times of events.
  
          Parameters
          ----------
          n_events : int
              Number of events to sample.
          t_min : `~astropy.time.Time`
              Start time of
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/temporal.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
  
  [... 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
  
  [... 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.time import Time
  from gammapy.data.gti import GTI
  from gammapy.datasets.map import MapEvaluator
  from gammapy.irf import EDispKernel, PSFKernel
  from gammapy.maps import Map, MapAxis, WcsGeom
  from gammapy.modeling.models import (
      BackgroundModel,
      ConstantSpectralModel,
      ConstantTemporalModel,
      GaussianSpatialModel,
      Models,
      PointSpatialModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      TemplateSpatialModel,
      CompoundSpectralModel,
      SkyModel,
      create_fermi_isotropic_diffuse_model,
  )
  from gammapy.utils.testing import requires_data
  
  
  @pytest.fixture(scope="session")
  def sky_model():
      spatial_model = GaussianSpatialModel(
          lon_0="3 deg", lat_0="4 deg", sigma="3 deg", frame="galactic"
      )
      spectral_model = PowerLawSpectralModel(
          index=2, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      temporal_model = ConstantTemporalModel()
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
          temporal_model=temporal_model,
          name="source-1",
      )
  
  
  @pytest.fixture(scope="session")
  def gti():
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      t_ref = Time(55555, format="mjd")
      gti = GTI.create(start, stop, reference_time=t_ref)
      return gti
  
  
  @pytest.fixture(scope="session")
  def diffuse_model():
      axis = MapAxis.from_nodes([0.1, 100], name="energy_true", unit="TeV", interp="log")
      m = Map.create(
          npix=(4, 3), binsz=2, axes=[axis], unit="cm-2 s-1 MeV-1 sr-1", frame="galactic"
      )
      m.data += 42
      spatial_model = TemplateSpatialModel(m, normalize=False)
      return SkyModel(PowerLawNormSpectralModel(), spatial_model)
  
  
  @pytest.fixture(scope="session")
  def geom():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 3), unit=u.TeV, name="energy")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def geom_true():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 4), unit=u.TeV, name="energy_true")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def exposure(geom_true):
      m = Map.from_geom(geom_true)
      m.quantity = np.ones(geom_true.data_shape) * u.Quantity("100 m2 s")
      m.data[1] *= 10
      return m
  
  
  @pytest.fixture(scope="session")
  def background(geom):
      m = Map.from_geom(geom)
     
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_cube.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.time import Time
  from gammapy.data.gti import GTI
  from gammapy.datasets.map import MapEvaluator
  from gammapy.irf import EDispKernel, PSFKernel
  from gammapy.maps import Map, MapAxis, WcsGeom
  from gammapy.modeling.models import (
      BackgroundModel,
      ConstantSpectralModel,
      ConstantTemporalModel,
      GaussianSpatialModel,
      Models,
      PointSpatialModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      TemplateSpatialModel,
      CompoundSpectralModel,
      SkyModel,
      create_fermi_isotropic_diffuse_model,
  )
  from gammapy.utils.testing import requires_data
  
  
  @pytest.fixture(scope="session")
  def sky_model():
      spatial_model = GaussianSpatialModel(
          lon_0="3 deg", lat_0="4 deg", sigma="3 deg", frame="galactic"
      )
      spectral_model = PowerLawSpectralModel(
          index=2, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      temporal_model = ConstantTemporalModel()
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
          temporal_model=temporal_model,
          name="source-1",
      )
  
  
  @pytest.fixture(scope="session")
  def gti():
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      t_ref = Time(55555, format="mjd")
      gti = GTI.create(start, stop, reference_time=t_ref)
      return gti
  
  
  @pytest.fixture(scope="session")
  def diffuse_model():
      axis = MapAxis.from_nodes([0.1, 100], name="energy_true", unit="TeV", interp="log")
      m = Map.create(
          npix=(4, 3), binsz=2, axes=[axis], unit="cm-2 s-1 MeV-1 sr-1", frame="galactic"
      )
      m.data += 42
      spatial_model = TemplateSpatialModel(m, normalize=False)
      return SkyModel(
          PowerLawNormSpectralModel(),
          spatial_model
      )
  
  
  @pytest.fixture(scope="session")
  def geom():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 3), unit=u.TeV, name="energy")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def geom_true():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 4), unit=u.TeV, name="energy_true")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def exposure(geom_true):
      m = Map.from_geom(geom_true)
      m.quantity = np.ones(geom_true.data_shape) * u.Quantity("100 m2 s")
      m.data[1] *= 10
      return m
  
  
  @pytest.fixture(scope="session")
  def background(geom):
      m = Map.from_
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_cube.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,
      PowerLawSpectralModel,
      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 model
  
  [... 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 "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
  
  [... 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
  import operator
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.maps import MapAxis
  from gammapy.modeling.models import (
      SPECTRAL_MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      ConstantSpectralModel,
      CompoundSpectralModel,
      ExpCutoffPowerLaw3FGLSpectralModel,
      ExpCutoffPowerLawSpectralModel,
      ExpCutoffPowerLawNormSpectralModel,
      GaussianSpectralModel,
      LogParabolaSpectralModel,
      LogParabolaNormSpectralModel,
      NaimaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      BrokenPowerLawSpectralModel,
      SmoothBrokenPowerLawSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
      TemplateSpectralModel,
  )
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      mpl_plot_check,
      requires_data,
      requires_dependency,
  )
  
  
  def table_model():
      energy = MapAxis.from_energy_bounds(0.1 * u.TeV, 100 * u.TeV, 1000).center
  
      model = PowerLawSpectralModel(
          index=2.3, amplitude="4 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      dnde = model(energy)
  
      return TemplateSpectralModel(energy, dnde)
  
  
  TEST_MODELS = [
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.6, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV cm-2 s-1"),
      ),
      dict(
          name="norm-powerlaw",
          model=PowerLawNormSpectralModel(
              tilt=2 * u.Unit(""), norm=4.0 * u.Unit(""), reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, ""),
          integral_1_10TeV=u.Quantity(3.6, "TeV"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV2"),
      ),
      dict(
          name="powerlaw2",
          model=PowerLaw2S
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_spectral.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  import operator
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.maps import MapAxis
  from gammapy.modeling.models import (
      SPECTRAL_MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      ConstantSpectralModel,
      CompoundSpectralModel,
      ExpCutoffPowerLaw3FGLSpectralModel,
      ExpCutoffPowerLawSpectralModel,
      ExpCutoffPowerLawNormSpectralModel,
      GaussianSpectralModel,
      LogParabolaSpectralModel,
      LogParabolaNormSpectralModel,
      NaimaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      BrokenPowerLawSpectralModel,
      SmoothBrokenPowerLawSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
      TemplateSpectralModel,
  )
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      mpl_plot_check,
      requires_data,
      requires_dependency,
  )
  
  
  def table_model():
      energy = MapAxis.from_energy_bounds(0.1 * u.TeV, 100 * u.TeV, 1000).center
  
      model = PowerLawSpectralModel(
          index=2.3, amplitude="4 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      dnde = model(energy)
  
      return TemplateSpectralModel(energy, dnde)
  
  
  TEST_MODELS = [
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.6, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV cm-2 s-1"),
      ),
      dict(
          name="norm-powerlaw",
          model=PowerLawNormSpectralModel(
              tilt=2 * u.Unit(""), norm=4.0 * u.Unit(""), reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, ""),
          integral_1_10TeV=u.Quantity(3.6, "TeV"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV2"),
      ),
      dict(
          name="powerlaw2",
          model=PowerLaw2S
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_spectral.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.table import Table
  from astropy.time import Time
  from gammapy.data.gti import GTI
  from gammapy.modeling.models import (
      ConstantTemporalModel,
      ExpDecayTemporalModel,
      GaussianTemporalModel,
      LightCurveTemplateTemporalModel,
      PowerLawSpectralModel,
      SkyModel,
  )
  from gammapy.utils.scripts import make_path
  from gammapy.utils.testing import requires_data, requires_dependency, mpl_plot_check
  
  
  # TODO: add light-curve test case from scratch
  # only use the FITS one for I/O (or not at all)
  @pytest.fixture(scope="session")
  def light_curve():
      path = "$GAMMAPY_DATA/tests/models/light_curve/lightcrv_PKSB1222+216.fits"
      return LightCurveTemplateTemporalModel.read(path)
  
  
  @requires_data()
  def test_light_curve_str(light_curve):
      ss = str(light_curve)
      assert "LightCurveTemplateTemporalModel" in ss
  
  
  @requires_data()
  def test_light_curve_evaluate(light_curve):
      t = Time(59500, format="mjd")
      val = light_curve(t)
      assert_allclose(val, 0.015512, rtol=1e-5)
  
      t = Time(46300, format="mjd")
      val = light_curve.evaluate(t.mjd, ext=3)
      assert_allclose(val, 0.01551196, rtol=1e-5)
  
  
  def rate(x, c="1e4 s"):
      c = u.Quantity(c)
      return np.exp(-x / c)
  
  
  def ph_curve(x, amplitude=0.5, x0=0.01):
      return 100.0 + amplitude * np.sin(2 * np.pi * (x - x0) / 1.0)
  
  
  def test_time_sampling(tmp_path):
      time = np.arange(0, 10, 0.06) * u.hour
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = rate(time)
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model = LightCurveTemplateTemporalModel(table)
  
      filename = str(make_path(tmp_path / "tmp.fits"))
      temporal_model.write(path=filename)
      model_read = temporal_model.read(filename)
      assert temporal_model.filename == filename
      assert model_read.filename == filename
      assert_allclose(model_read.table["TIME"].quantity.value, time.value)
  
      t_ref = "2010-01-01T00:00:00"
      t_min = "2010-01-01T00:00:00"
      t_max = "2010-01-01T08:00:00"
  
      sampler = temporal_model.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0, t_delta="10 min"
      )
  
      sampler = u.Quantity((sampler - Time(t_ref)).sec, "s")
      assert len(sampler) == 2
      assert_allclose(sampler.value, [12661.65802564, 7826.92991], rtol=1e-5)
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = np.ones(len(time
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_temporal.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.table import Table
  from astropy.time import Time
  from gammapy.data.gti import GTI
  from gammapy.modeling.models import (
      ConstantTemporalModel,
      ExpDecayTemporalModel,
      GaussianTemporalModel,
      LightCurveTemplateTemporalModel,
      PowerLawSpectralModel,
      SkyModel,
  )
  from gammapy.utils.scripts import make_path
  from gammapy.utils.testing import requires_data, requires_dependency, mpl_plot_check
  
  
  # TODO: add light-curve test case from scratch
  # only use the FITS one for I/O (or not at all)
  @pytest.fixture(scope="session")
  def light_curve():
      path = "$GAMMAPY_DATA/tests/models/light_curve/lightcrv_PKSB1222+216.fits"
      return LightCurveTemplateTemporalModel.read(path)
  
  
  @requires_data()
  def test_light_curve_str(light_curve):
      ss = str(light_curve)
      assert "LightCurveTemplateTemporalModel" in ss
  
  
  @requires_data()
  def test_light_curve_evaluate(light_curve):
      t = Time(59500, format="mjd")
      val = light_curve(t)
      assert_allclose(val, 0.015512, rtol=1e-5)
  
      t = Time(46300, format="mjd")
      val = light_curve.evaluate(t.mjd, ext=3)
      assert_allclose(val, 0.01551196, rtol=1e-5)
  
  
  def rate(x, c="1e4 s"):
      c = u.Quantity(c)
      return np.exp(-x / c)
  
  
  def ph_curve(x, amplitude=0.5, x0=0.01):
      return 100.0 + amplitude * np.sin(2 * np.pi * (x - x0) / 1.0)
  
  
  def test_time_sampling(tmp_path):
      time = np.arange(0, 10, 0.06) * u.hour
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = rate(time)
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model = LightCurveTemplateTemporalModel(table)
  
      filename = str(make_path(tmp_path / "tmp.fits"))
      temporal_model.write(path=filename)
      model_read = temporal_model.read(filename)
      assert temporal_model.filename == filename
      assert model_read.filename == filename
      assert_allclose(model_read.table["TIME"].quantity.value, time.value)
  
      t_ref = "2010-01-01T00:00:00"
      t_min = "2010-01-01T00:00:00"
      t_max = "2010-01-01T08:00:00"
  
      sampler = temporal_model.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0, t_delta="10 min"
      )
  
      sampler = u.Quantity((sampler - Time(t_ref)).sec, "s")
      assert len(sampler) == 2
      assert_allclose(sampler.value, [12661.65802564, 7826.92991], rtol=1e-5)
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = np.ones(len(time
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/tests/test_temporal.py