files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import logging
  import numpy as np
  import astropy.units as u
  from astropy.io import fits
  from astropy.nddata.utils import NoOverlapError
  from astropy.table import Table
  from astropy.utils import lazyproperty
  from regions import CircleSkyRegion
  from gammapy.cube.edisp_map import EDispMap
  from gammapy.cube.psf_kernel import PSFKernel
  from gammapy.cube.psf_map import PSFMap
  from gammapy.data import GTI
  from gammapy.irf import EDispKernel, EffectiveAreaTable
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling import Dataset, Parameters
  from gammapy.modeling.models import BackgroundModel, Models
  from gammapy.spectrum import SpectrumDataset, SpectrumDatasetOnOff
  from gammapy.stats import cash, cash_sum_cython, wstat
  from gammapy.utils.random import get_random_state
  from gammapy.utils.scripts import make_name, make_path
  from .exposure import _map_spectrum_weight
  
  __all__ = ["MapDataset", "MapDatasetOnOff"]
  
  log = logging.getLogger(__name__)
  
  CUTOUT_MARGIN = 0.1 * u.deg
  RAD_MAX = 0.66
  RAD_AXIS_DEFAULT = MapAxis.from_bounds(
      0, RAD_MAX, nbin=66, node_type="edges", name="theta", unit="deg"
  )
  MIGRA_AXIS_DEFAULT = MapAxis.from_bounds(
      0.2, 5, nbin=48, node_type="edges", name="migra"
  )
  
  BINSZ_IRF_DEFAULT = 0.2
  
  
  class MapDataset(Dataset):
      """Perform sky model likelihood fit on maps.
  
      Parameters
      ----------
      models : `~gammapy.modeling.models.Models`
          Source sky models.
      counts : `~gammapy.maps.WcsNDMap`
          Counts cube
      exposure : `~gammapy.maps.WcsNDMap`
          Exposure cube
      mask_fit : `~gammapy.maps.WcsNDMap`
          Mask to apply to the likelihood for fitting.
      psf : `~gammapy.cube.PSFKernel` or `~gammapy.cube.PSFMap`
          PSF kernel
      edisp : `~gammapy.irf.EDispKernel` or `~gammapy.cube.EDispMap`
          Energy dispersion kernel
      evaluation_mode : {"local", "global"}
          Model evaluation mode.
          The "local" mode evaluates the model components on smaller grids to save computation time.
          This mode is recommended for local optimization algorithms.
          The "global" evaluation mode evaluates the model components on the full map.
          This mode is recommended for global optimization algorithms.
      mask_safe : `~gammapy.maps.WcsNDMap`
          Mask defining the safe data range.
      gti : `~gammapy.data.GTI`
          GTI of the observation or union of GTI if it is a stacked observation
      """
  
      stat_type = "cash"
      tag = "MapDataset"
  
      def __init__(
          self,
          models=None,
          counts=None,
          exposure=None,
          mask_fit=None,
          psf=None,
          edisp=None,
          name=None,
          evaluation_mode="local",
          mask_safe=None,
          gti=None,
      ):
          if mask_fit is not None and mask_fit.data.dtype != np.dtype("bool"):
              raise ValueError("mask data must have dtype bool")
  
          if mask_safe is not None and mask_safe.data.dtype != np.dtype("bool"):
              raise ValueError("mask data must have dtype bool")
  
          self._name = make_name(name)
          self._evaluators = {}
  
          self.background_model = None
          self.evaluation_mode = evaluation_mode
          self.counts = counts
          self.exposure = exposure
          self.mask_fit = mask_fit
          self.psf = psf
          self.edisp = edisp
          self.mask_safe = mask_safe
          self.models = models
          self.gti = gti
  
          # check whether a reference geom is defined
          _ = self._geom
  
      @property
      def name(self):
          return self._name
  
      def __str__(self):
          str_ = f"{self.__class__.__name__}\n"
          str_ += "-" * len(self.__class__.__name__) + "\n"
          str_ += "\n"
  
          str_ += "\t{:32}: {} \n\n".format("Name", self.name)
  
          counts = np.nan
          if self.counts is not None:
              counts = np.sum(self.counts.data)
          str_ += "\t{:32}: {:.0f} \n".format("Total counts", counts)
  
          npred = np.nan
          if self.models is not None:
              npred = np.sum(self.npred().data)
          str_ += "\t{:32}: {:.2f}\n".format("Total predicted counts", npred)
  
          background = np.nan
          if self.background_model is not None:
              background = np.sum(self.background_model.evaluate().data)
          str_ += "\t{:32}: {:.2f}\n\n".format("Total background counts", background)
  
          exposure_min, exposure_max, exposure_unit = np.nan, np.nan, ""
          if self.exposure is not None:
              if self.mask_safe is not None:
                  mask = self.mask_safe.reduce_over_axes(np.logical_or).data
                  if not mask.any():
                      mask = None
              else:
                  mask = None
              exposure_min = np.min(self.exposure.data[..., mask])
              exposure_max = np.max(self.exposure.data[..., mask])
              exposure_unit = self.exposure.unit
  
          str_ += "\t{:32}: {:.2e} {}\n".format(
              "Exposure min", exposure_min, exposure_unit
          )
          str_ += "\t{:32}: {:.2e} {}\n\n".format(
              "Exposure max", exposure_max, exposure_unit
          )
  
          # data section
          n_bins = 0
          if self.counts is not None:
              n_bins = self.counts.data.size
          str_ += "\t{:32}: {} \n".format("Number of total bins", n_bins)
  
          n_fit_bins = 0
          if self.mask is not None:
              n_fit_bins = np.sum(self.mask.data)
          str_ += "\t{:32}: {} \n\n".format("Number of fit bins", n_fit_bins)
  
          # likelihood section
          str_ += "\t{:32}: {}\n".format("Fit statistic type", self.stat_type)
  
          stat = np.nan
          if self.counts is not None and self.models is not None:
              stat = self.stat_sum()
          str_ += "\t{:32}: {:.2f}\n\n".format("Fit statistic value (-2 log(L))", stat)
  
          # model section
          n_models = 0
          if self.models is not None:
              n_models = len(self.models)
  
          str_ += "\t{:32}: {} \n".format("Number of models", n_models)
          str_ += "\t{:32}: {}\n".format("Number of parameters", len(self.parameters))
          str_ += "\t{:32}: {}\n\n".format(
              "Number of free parameters", len(self.parameters.free_parameters)
          )
  
          if self.models is not None:
              str_ += "\t" + "\n\t".join(str(self.models).split("\n")[2:])
  
          return str_.expandtabs(tabsize=2)
  
      @property
      def models(self):
          """Models (`~gammapy.modeling.models.Models`)."""
          return self._models
  
      @models.setter
      def models(self, value):
          if value is not None:
              self._models = Models(value)
          else:
              self._models = None
  
          if self.models is not None:
              for model in self.models:
                  if isinstance(model, BackgroundModel):
                      self.background_model = model
                      break
              else:
                  log.warning(f"No background model defined for dataset {self.name}")
  
      @property
      def evaluators(self):
          """Model evaluators"""
          # this call is needed to trigger the setup of the evaluators
          if not self._evaluators:
              self.npred()
  
          return self._evaluators
  
      @property
      def parameters(self):
          """List of parameters (`~gammapy.modeling.Parameters`)"""
          return self.models.parameters
  
      @property
      def _geom(self):
          """Main analysis geometry"""
          if self.counts is not None:
              return self.counts.geom
          elif self.background_model is not None:
              return self.background_model.map.geom
          elif self.mask_safe is not None:
              return self.mask_safe.geom
          elif self.mask_fit is not None:
              return self.mask_fit.geom
          else:
              raise ValueError(
                  "Either 'counts', 'background_model', 'mask_fit'"
                  " or 'mask_safe' must be defined."
              )
  
      @property
      def data_shape(self):
          """Shape of the counts or background data (tuple)"""
          return self._geom.data_shape
  
      def npred(self):
          """Predicted source and background counts (`~gammapy.maps.Map`)."""
          npred_total = Map.from_geom(self._geom, dtype=float)
  
          if self.models:
              for model in self.models:
                  evaluator = self._evaluators.get(model.name)
  
                  if evaluator is None:
                      evaluator = MapEvaluator(
                          model=model, evaluation_mode=self.evaluation_mode
                      )
                      self._evaluators[model.name] = evaluator
  
                  # if the model component drifts out of its support the evaluator has
                  # has to be updated
                  if evaluator.needs_update:
                      evaluator.update(self.exposure, self.psf, self.edisp, self._geom)
  
                  if evaluator.contributes:
                      npred = evaluator.compute_npred()
                      npred_total.stack(npred)
  
          return npred_total
  
      @classmethod
      def from_geoms(
          cls,
          geom,
          geom_exposure,
          geom_psf,
          geom_edisp,
          reference_time="2000-01-01",
          name=None,
          **kwargs,
      ):
          """
          Create a MapDataset object with zero filled maps according to the specified geometries
  
          Parameters
          ----------
          geom : `Geom`
              geometry for the counts and background maps
          geom_exposure : `Geom`
              geometry for the exposure map
          geom_psf : `Geom`
              geometry for the psf map
          geom_edisp : `Geom`
              geometry for the energy dispersion map
          reference_time : `~astropy.time.Time`
              the reference time to use in GTI definition
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          empty_maps : `MapDataset`
              A MapDataset containing zero filled maps
          """
          kwargs = kwargs.copy()
          kwargs["name"] = name
          kwargs["counts"] = Map.from_geom(geom, unit="")
  
          background = Map.from_geom(geom, unit="")
          kwargs["models"] = Models([BackgroundModel(background)])
          kwargs["exposure"] = Map.from_geom(geom_exposure, unit="m2 s")
          kwargs["edisp"] = EDispMap.from_geom(geom_edisp)
          kwargs["psf"] = PSFMap.from_geom(geom_psf)
          kwargs["gti"] = GTI.create([] * u.s, [] * u.s, reference_time=reference_time)
          kwargs["mask_safe"] = Map.from_geom(geom, unit="", dtype=bool)
  
          return cls(**kwargs)
  
      @classmethod
      def create(
          cls,
          geom,
          energy_axis_true=None,
          migra_axis=None,
          rad_axis=None,
          binsz_irf=None,
          reference_time="2000-01-01",
          name=None,
          **kwargs,
      ):
          """Create a MapDataset object with zero filled maps.
  
          Parameters
          ----------
          geom : `~gammapy.maps.WcsGeom`
              Reference target geometry in reco energy, used for counts and background maps
          energy_axis_true : `~gammapy.maps.MapAxis`
              True energy axis used for IRF maps
          migra_axis : `~gammapy.maps.MapAxis`
              Migration axis for the energy dispersion map
          rad_axis : `~gammapy.maps.MapAxis`
              Rad axis for the psf map
          binsz_irf : float
              IRF Map pixel size in degrees.
          reference_time : `~astropy.time.Time`
              the reference time to use in GTI definition
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          empty_maps : `MapDataset`
              A MapDataset containing zero filled maps
          """
          migra_axis = migra_axis or MIGRA_AXIS_DEFAULT
          rad_axis = rad_axis or RAD_AXIS_DEFAULT
          energy_axis_true = energy_axis_true or geom.get_axis_by_name("energy")
          binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT
  
          geom_image = geom.to_image()
          geom_exposure = geom_image.to_cube([energy_axis_true])
          geom_irf = geom_image.to_binsz(binsz=binsz_irf)
          geom_psf = geom_irf.to_cube([rad_axis, energy_axis_true])
          geom_edisp = geom_irf.to_cube([migra_axis, energy_axis_true])
  
          return cls.from_geoms(
              geom,
              geom_exposure,
              geom_psf,
              geom_edisp,
              reference_time=reference_time,
              name=name,
              **kwargs,
          )
  
      def stack(self, other):
          """Stack another dataset in place.
  
          Parameters
          ----------
          other: `~gammapy.cube.MapDataset`
              Map dataset to be stacked with this one.
          """
  
          if self.counts and other.counts:
              self.counts *= self.mask_safe
              self.counts.stack(other.counts, weights=other.mask_safe)
  
          if self.exposure and other.exposure:
              mask_image = self.mask_safe.reduce_over_axes(func=np.logical_or)
              self.exposure *= mask_image.data
              # TODO: apply energy dependent mask to exposure. Does this require
              #  a mask_safe in true energy?
              mask_image_other = other.mask_safe.reduce_over_axes(func=np.logical_or)
              self.exposure.stack(other.exposure, weights=mask_image_other)
  
          if self.background_model and other.background_model:
              self.background_model.map *= self.mask_safe
              self.background_model.stack(other.background_model, other.mask_safe)
  
          if self.mask_safe is not None and other.mask_safe is not None:
              self.mask_safe.stack(other.mask_safe)
  
          if self.psf and other.psf:
              if isinstance(self.psf, PSFMap) and isinstance(other.psf, PSFMap):
                  mask_irf = self._mask_safe_irf(self.psf.psf_map, mask_image)
                  self.psf.psf_map *= mask_irf.data
                  self.psf.exposure_map *= mask_irf.data
  
                  mask_image_other = other.mask_safe.reduce_over_axes(func=np.logical_or)
                  mask_irf_other = self._mask_safe_irf(
                      other.psf.psf_map, mask_image_other
                  )
                  self.psf.stack(other.psf, weights=mask_irf_other)
              else:
                  raise ValueError("Stacking of PSF kernels not supported")
  
          if self.edisp and other.edisp:
              if isinstance(self.edisp, EDispMap) and isinstance(other.edisp, EDispMap):
                  mask_irf = self._mask_safe_irf(self.edisp.edisp_map, mask_image)
                  self.edisp.edisp_map *= mask_irf.data
                  self.edisp.exposure_map *= mask_irf.data
  
                  mask_image_other = other.mask_safe.reduce_over_axes(func=np.logical_or)
                  mask_irf_other = self._mask_safe_irf(
                      other.edisp.edisp_map, mask_image_other
                  )
                  self.edisp.stack(other.edisp, weights=mask_irf_other)
              else:
                  raise ValueError("Stacking of edisp kernels not supported")
  
          if self.gti and other.gti:
              self.gti = self.gti.stack(other.gti).union()
  
      @staticmethod
      def _mask_safe_irf(irf_map, mask):
          geom = irf_map.geom.to_image()
          coords = geom.get_coord()
          data = mask.get_by_coord(coords).astype(bool)
          return Map.from_geom(geom=geom, data=data)
  
      def stat_array(self):
          """Likelihood per bin given the current model parameters"""
          return cash(n_on=self.counts.data, mu_on=self.npred().data)
  
      def residuals(self, method="diff"):
          """Compute residuals map.
  
          Parameters
          ----------
          method: {"diff", "diff/model", "diff/sqrt(model)"}
              Method used to compute the residuals. Available options are:
                  - "diff" (default): data - model
                  - "diff/model": (data - model) / model
                  - "diff/sqrt(model)": (data - model) / sqrt(model)
  
          Returns
          -------
          residuals : `gammapy.maps.WcsNDMap`
              Residual map.
          """
          return self._compute_residuals(self.counts, self.npred(), method=method)
  
      def plot_residuals(
          self,
          method="diff",
          smooth_kernel="gauss",
          smooth_radius="0.1 deg",
          region=None,
          figsize=(12, 4),
          **kwargs,
      ):
          """
          Plot spatial and spectral residuals.
  
          The spectral residuals are extracted from the provided region, and the
          normalization used for the residuals computation can be controlled using
          the method parameter. If no region is passed, only the spatial
          residuals are shown.
  
          Parameters
          ----------
          method : {"diff", "diff/model", "diff/sqrt(model)"}
              Method used to compute the residuals, see `MapDataset.residuals()`
          smooth_kernel : {'gauss', 'box'}
              Kernel shape.
          smooth_radius: `~astropy.units.Quantity`, str or float
              Smoothing width given as quantity or float. If a float is given it
              is interpreted as smoothing width in pixels.
          region: `~regions.Region`
              Region (pixel or sky regions accepted)
          figsize : tuple
              Figure size used for the plotting.
          **kwargs : dict
              Keyword arguments passed to `~matplotlib.pyplot.imshow`.
  
          Returns
          -------
          ax_image, ax_spec : `~matplotlib.pyplot.Axes`,
              Image and spectrum axes.
          """
          import matplotlib.pyplot as plt
  
          fig = plt.figure(figsize=figsize)
  
          counts, npred = self.counts, self.npred()
  
          if self.mask is not None:
              counts = counts * self.mask
              npred = npred * self.mask
  
          counts_spatial = counts.sum_over_axes().smooth(
              width=smooth_radius, kernel=smooth_kernel
          )
          npred_spatial = npred.sum_over_axes().smooth(
              width=smooth_radius, kernel=smooth_kernel
          )
          spatial_residuals = self._compute_residuals(
              counts_spatial, npred_spatial, method
          )
  
          if self.mask_safe is not None:
              mask = self.mask_safe.reduce_over_axes(func=np.logical_or)
              spatial_residuals.data[~mask.data] = np.nan
  
          # If no region is provided, skip spectral residuals
          ncols = 2 if region is not None else 1
          ax_image = fig.add_subplot(1, ncols, 1, projection=spatial_residuals.geom.wcs)
          ax_spec = None
  
          kwargs.setdefault("cmap", "coolwarm")
          kwargs.setdefault("stretch", "linear")
          kwargs.setdefault("vmin", -5)
          kwargs.setdefault("vmax", 5)
          spatial_residuals.plot(ax=ax_image, add_cbar=True, **kwargs)
  
          # Spectral residuals
          if region:
              ax_spec = fig.add_subplot(1, 2, 2)
              counts_spec = counts.get_spectrum(region=region)
              npred_spec = npred.get_spectrum(region=region)
              residuals = self._compute_residuals(counts_spec, npred_spec, method)
              ax = residuals.plot()
              ax.axhline(0, color="black", lw=0.5)
  
              y_max = 2 * np.nanmax(residuals.data)
              plt.ylim(-y_max, y_max)
              label = self._residuals_labels[method]
              plt.ylabel(f"Residuals ({label})")
  
              # Overlay spectral extraction region on the spatial residuals
              pix_region = region.to_pixel(wcs=spatial_residuals.geom.wcs)
              pix_region.plot(ax=ax_image)
  
          return ax_image, ax_spec
  
      @lazyproperty
      def _counts_data(self):
          return self.counts.data.astype(float)
  
      def stat_sum(self):
          """Total likelihood given the current model parameters."""
          counts, npred = self._counts_data, self.npred().data
  
          if self.mask is not None:
              return cash_sum_cython(counts[self.mask.data], npred[self.mask.data])
          else:
              return cash_sum_cython(counts.ravel(), npred.ravel())
  
      def fake(self, random_state="random-seed"):
          """Simulate fake counts for the current model and reduced IRFs.
  
          This method overwrites the counts defined on the dataset object.
  
          Parameters
          ----------
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
                  Defines random number generator initialisation.
                  Passed to `~gammapy.utils.random.get_random_state`.
          """
          random_state = get_random_state(random_state)
          npred = self.npred()
          npred.data = random_state.poisson(npred.data)
          self.counts = npred
  
      def to_hdulist(self):
          """Convert map dataset to list of HDUs.
  
          Returns
          -------
          hdulist : `~astropy.io.fits.HDUList`
              Map dataset list of HDUs.
          """
          # TODO: what todo about the model and background model parameters?
          exclude_primary = slice(1, None)
  
          hdu_primary = fits.PrimaryHDU()
          hdulist = fits.HDUList([hdu_primary])
          if self.counts is not None:
              hdulist += self.counts.to_hdulist(hdu="counts")[exclude_primary]
  
          if self.exposure is not None:
              hdulist += self.exposure.to_hdulist(hdu="exposure")[exclude_primary]
  
          if self.background_model is not None:
              hdulist += self.background_model.map.to_hdulist(hdu="background")[
                  exclude_primary
              ]
  
          if self.edisp is not None:
              if isinstance(self.edisp, EDispKernel):
                  hdus = self.edisp.to_hdulist()
                  hdus["MATRIX"].name = "edisp_matrix"
                  hdus["EBOUNDS"].name = "edisp_matrix_ebounds"
                  hdulist.append(hdus["EDISP_MATRIX"])
                  hdulist.append(hdus["EDISP_MATRIX_EBOUNDS"])
              else:
                  hdulist += self.edisp.edisp_map.to_hdulist(hdu="EDISP")[exclude_primary]
                  hdulist += self.edisp.exposure_map.to_hdulist(hdu="edisp_exposure")[
                      exclude_primary
                  ]
  
          if self.psf is not None:
              if isinstance(self.psf, PSFKernel):
                  hdulist += self.psf.psf_kernel_map.to_hdulist(hdu="psf_kernel")[
                      exclude_primary
                  ]
              else:
                  hdulist += self.psf.psf_map.to_hdulist(hdu="psf")[exclude_primary]
                  hdulist += self.psf.exposure_map.to_hdulist(hdu="psf_exposure")[
                      exclude_primary
                  ]
  
          if self.mask_safe is not None:
              mask_safe_int = self.mask_safe.copy()
              mask_safe_int.data = mask_safe_int.data.astype(int)
              hdulist += mask_safe_int.to_hdulist(hdu="mask_safe")[exclude_primary]
  
          if self.mask_fit is not None:
              mask_fit_int = self.mask_fit.copy()
              mask_fit_int.data = mask_fit_int.data.astype(int)
              hdulist += mask_fit_int.to_hdulist(hdu="mask_fit")[exclude_primary]
  
          if self.gti is not None:
              hdulist.append(fits.BinTableHDU(self.gti.table, name="GTI"))
  
          return hdulist
  
      @classmethod
      def from_hdulist(cls, hdulist, name=None):
          """Create map dataset from list of HDUs.
  
          Parameters
          ----------
          hdulist : `~astropy.io.fits.HDUList`
              List of HDUs.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset.
          """
          kwargs = {"name": name}
  
          if "COUNTS" in hdulist:
              kwargs["counts"] = Map.from_hdulist(hdulist, hdu="counts")
  
          if "EXPOSURE" in hdulist:
              kwargs["exposure"] = Map.from_hdulist(hdulist, hdu="exposure")
  
          if "BACKGROUND" in hdulist:
              background_map = Map.from_hdulist(hdulist, hdu="background")
              kwargs["models"] = Models([BackgroundModel(background_map)])
  
          if "EDISP_MATRIX" in hdulist:
              kwargs["edisp"] = EDispKernel.from_hdulist(
                  hdulist, hdu1="EDISP_MATRIX", hdu2="EDISP_MATRIX_EBOUNDS"
              )
          if "EDISP" in hdulist:
              edisp_map = Map.from_hdulist(hdulist, hdu="edisp")
              exposure_map = Map.from_hdulist(hdulist, hdu="edisp_exposure")
              kwargs["edisp"] = EDispMap(edisp_map, exposure_map)
  
          if "PSF_KERNEL" in hdulist:
              psf_map = Map.from_hdulist(hdulist, hdu="psf_kernel")
              kwargs["psf"] = PSFKernel(psf_map)
          if "PSF" in hdulist:
              psf_map = Map.from_hdulist(hdulist, hdu="psf")
              exposure_map = Map.from_hdulist(hdulist, hdu="psf_exposure")
              kwargs["psf"] = PSFMap(psf_map, exposure_map)
  
          if "MASK_SAFE" in hdulist:
              mask_safe = Map.from_hdulist(hdulist, hdu="mask_safe")
              mask_safe.data = mask_safe.data.astype(bool)
              kwargs["mask_safe"] = mask_safe
  
          if "MASK_FIT" in hdulist:
              mask_fit = Map.from_hdulist(hdulist, hdu="mask_fit")
              mask_fit.data = mask_fit.data.astype(bool)
              kwargs["mask_fit"] = mask_fit
  
          if "GTI" in hdulist:
              gti = GTI(Table.read(hdulist, hdu="GTI"))
              kwargs["gti"] = gti
  
          return cls(**kwargs)
  
      def write(self, filename, overwrite=False):
          """Write map dataset to file.
  
          Parameters
          ----------
          filename : str
              Filename to write to.
          overwrite : bool
              Overwrite file if it exists.
          """
          self.to_hdulist().writeto(make_path(filename), overwrite=overwrite)
  
      @classmethod
      def read(cls, filename, name=None):
          """Read map dataset from file.
  
          Parameters
          ----------
          filename : str
              Filename to read from.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset.
          """
          with fits.open(make_path(filename), memmap=False) as hdulist:
              return cls.from_hdulist(hdulist, name=name)
  
      @classmethod
      def from_dict(cls, data, components, models):
          """Create from dicts and models list generated from YAML serialization."""
          dataset = cls.read(data["filename"], name=data["name"])
          bkg_name = data["background"]
          model_names = data["models"]
          models_list = [model for model in models if model.name in model_names]
          models = Models(models_list)
  
          for component in components["components"]:
              if component["type"] == "BackgroundModel":
                  if component["name"] == bkg_name:
                      if "filename" not in component:
                          component["map"] = dataset.background_model.map
                      background_model = BackgroundModel.from_dict(component)
                      models.append(background_model)
  
          dataset.models = models
          return dataset
  
      def to_dict(self, filename=""):
          """Convert to dict for YAML serialization."""
          if self.models is None:
              models = []
          else:
              models = [_.name for _ in self.models]
  
          return {
              "name": self.name,
              "type": self.tag,
              "models": models,
              "background": self.background_model.name,
              "filename": str(filename),
          }
  
      def to_spectrum_dataset(self, on_region, containment_correction=False, name=None):
          """Return a ~gammapy.spectrum.SpectrumDataset from on_region.
  
          Counts and background are summed in the on_region.
  
          Effective area is taken from the average exposure divided by the livetime.
          Here we assume it is the sum of the GTIs.
  
          The energy dispersion kernel is obtained at the on_region center.
          Only regions with centers are supported.
  
          The model is not exported to the ~gammapy.spectrum.SpectrumDataset.
          It must be set after the dataset extraction.
  
          Parameters
          ----------
          on_region : `~regions.SkyRegion`
              the input ON region on which to extract the spectrum
          containment_correction : bool
              Apply containment correction for point sources and circular on regions
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `~gammapy.spectrum.SpectrumDataset`
              the resulting reduced dataset
          """
          kwargs = {"gti": self.gti, "name": name}
  
          if self.gti is not None:
              kwargs["livetime"] = self.gti.time_sum
          else:
              raise ValueError("No GTI in `MapDataset`, cannot compute livetime")
  
          if self.counts is not None:
              kwargs["counts"] = self.counts.get_spectrum(on_region, np.sum)
  
          if self.background_model is not None:
              kwargs["background"] = self.background_model.evaluate().get_spectrum(
                  on_region, np.sum
              )
  
          if self.exposure is not None:
              exposure = self.exposure.get_spectrum(on_region, np.mean)
              kwargs["aeff"] = EffectiveAreaTable(
                  energy_lo=exposure.energy.edges[:-1],
                  energy_hi=exposure.energy.edges[1:],
                  data=exposure.quantity / kwargs["livetime"],
              )
  
          if containment_correction:
              if not isinstance(on_region, CircleSkyRegion):
                  raise TypeError(
                      "Containement correction is only supported for"
                      " `CircleSkyRegion`."
                  )
              elif self.psf is None or isinstance(self.psf, PSFKernel):
                  raise ValueError("No PSFMap set. Containement correction impossible")
              else:
                  psf = self.psf.get_energy_dependent_table_psf(on_region.center)
                  containment = psf.containment(
                      kwargs["aeff"].energy.center, on_region.radius
                  )
                  kwargs["aeff"].data.data *= containment.squeeze()
  
          if self.edisp is not None:
              if isinstance(self.edisp, EDispKernel):
                  edisp = self.edisp
              else:
                  axis = self._geom.get_axis_by_name("energy")
                  edisp = self.edisp.get_edisp_kernel(on_region.center, e_reco=axis.edges)
              kwargs["edisp"] = edisp
  
          return SpectrumDataset(**kwargs)
  
      def to_image(self, spectrum=None, name=None):
          """Create images by summing over the energy axis.
  
          Exposure is weighted with an assumed spectrum,
          resulting in a weighted mean exposure image.
  
          Currently the PSFMap and EdispMap are dropped from the
          resulting image dataset.
  
          Parameters
          ----------
          spectrum : `~gammapy.modeling.models.SpectralModel`
              Spectral model to compute the weights.
              Default is power-law with spectral index of 2.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset containing images.
          """
          kwargs = {}
          kwargs["name"] = make_name(name)
          kwargs["gti"] = self.gti
  
          if self.mask_safe is not None:
              mask_safe = self.mask_safe
              kwargs["mask_safe"] = mask_safe.reduce_over_axes(
                  func=np.logical_or, keepdims=True
              )
          else:
              mask_safe = 1
  
          if self.counts is not None:
              counts = self.counts * mask_safe
              kwargs["counts"] = counts.sum_over_axes(keepdims=True)
  
          if self.exposure is not None:
              exposure = _map_spectrum_weight(self.exposure, spectrum)
              kwargs["exposure"] = exposure.sum_over_axes(keepdims=True)
  
          if self.background_model is not None:
              background = self.background_model.evaluate() * mask_safe
              background = background.sum_over_axes(keepdims=True)
              kwargs["models"] = Models([BackgroundModel(background)])
  
          if self.psf is not None:
              kwargs["psf"] = self.psf.to_image(spectrum=spectrum, keepdims=True)
  
          return self.__class__(**kwargs)
  
      def cutout(self, position, width, mode="trim", name=None):
          """Cutout map dataset.
  
          Parameters
          ----------
          position : `~astropy.coordinates.SkyCoord`
              Center position of the cutout region.
          width : tuple of `~astropy.coordinates.Angle`
              Angular sizes of the region in (lon, lat) in that specific order.
              If only one value is passed, a square region is extracted.
          mode : {'trim', 'partial', 'strict'}
              Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          cutout : `MapDataset`
              Cutout map dataset.
          """
          name = make_name(name)
          kwargs = {"gti": self.gti, "name": name}
          cutout_kwargs = {"position": position, "width": width, "mode": mode}
  
          if self.counts is not None:
              kwargs["counts"] = self.counts.cutout(**cutout_kwargs)
  
          if self.exposure is not None:
              kwargs["exposure"] = self.exposure.cutout(**cutout_kwargs)
  
          if self.background_model is not None:
              kwargs["models"] = self.background_model.cutout(**cutout_kwargs, name=name)
  
          if self.edisp is not None:
              kwargs["edisp"] = self.edisp.cutout(**cutout_kwargs)
  
          if self.psf is not None:
              kwargs["psf"] = self.psf.cutout(**cutout_kwargs)
  
          if self.mask_safe is not None:
              kwargs["mask_safe"] = self.mask_safe.cutout(**cutout_kwargs)
  
          if self.mask_fit is not None:
              kwargs["mask_fit"] = self.mask_fit.cutout(**cutout_kwargs)
  
          return self.__class__(**kwargs)
  
  
  class MapDatasetOnOff(MapDataset):
      """Map dataset for on-off likelihood fitting.
  
      Parameters
      ----------
      models : `~gammapy.modeling.models.Models`
          Source sky models.
      counts : `~gammapy.maps.WcsNDMap`
          Counts cube
      counts_off : `~gammapy.maps.WcsNDMap`
          Ring-convolved counts cube
      acceptance : `~gammapy.maps.WcsNDMap`
          Acceptance from the IRFs
      acceptance_off : `~gammapy.maps.WcsNDMap`
          Acceptance off
      exposure : `~gammapy.maps.WcsNDMap`
          Exposure cube
      mask_fit : `~numpy.ndarray`
          Mask to apply to the likelihood for fitting.
      psf : `~gammapy.cube.PSFKernel`
          PSF kernel
      edisp : `~gammapy.irf.EDispKernel`
          Energy dispersion
      background_model : `~gammapy.modeling.models.BackgroundModel`
          Background model to use for the fit.
      evaluation_mode : {"local", "global"}
          Model evaluation mode.
          The "local" mode evaluates the model components on smaller grids to save computation time.
          This mode is recommended for local optimization algorithms.
          The "global" evaluation mode evaluates the model components on the full map.
          This mode is recommended for global optimization algorithms.
      mask_safe : `~numpy.ndarray`
          Mask defining the safe data range.
      gti : `~gammapy.data.GTI`
          GTI of the observation or union of GTI if it is a stacked observation
      name : str
          Name of the dataset.
  
      """
  
      stat_type = "wstat"
      tag = "MapDatasetOnOff"
  
      def __init__(
          self,
          models=None,
          counts=None,
          counts_off=None,
          acceptance=None,
          acceptance_off=None,
          exposure=None,
          mask_fit=None,
          psf=None,
          edisp=None,
          background_model=None,
          name=None,
          evaluation_mode="local",
          mask_safe=None,
          gti=None,
      ):
          if mask_fit is not None and mask_fit.dtype != np.dtype("bool"):
              raise ValueError("mask data must have dtype bool")
  
          self.evaluation_mode = evaluation_mode
          self.counts = counts
          self.counts_off = counts_off
  
          if np.isscalar(acceptance):
              acceptance = np.ones(self.data_shape) * acceptance
  
          if np.isscalar(acceptance_off):
              acceptance_off = np.ones(self.data_shape) * acceptance_off
  
          self.acceptance = acceptance
          self.acceptance_off = acceptance_off
          self.exposure = exposure
          self.background_model = None
          self.mask_fit = mask_fit
          self.psf = psf
          self.edisp = edisp
          self.models = models
          self._name = make_name(name)
          self.mask_safe = mask_safe
          self.gti = gti
  
      def __str__(self):
          str_ = super().__str__()
  
          counts_off = np.nan
          if self.counts_off is not None:
              counts_off = np.sum(self.counts_off.data)
          str_ += "\t{:32}: {:.0f} \n".format("Total counts_off", counts_off)
  
          acceptance = np.nan
          if self.acceptance is not None:
              acceptance = np.sum(self.acceptance.data)
          str_ += "\t{:32}: {:.0f} \n".format("Acceptance", acceptance)
  
          acceptance_off = np.nan
          if self.acceptance_off is not None:
              acceptance_off = np.sum(self.acceptance_off.data)
          str_ += "\t{:32}: {:.0f} \n".format("Acceptance off", acceptance_off)
  
          return str_.expandtabs(tabsize=4)
  
      @property
      def parameters(self):
          """List of parameters (`~gammapy.modeling.Parameters`)"""
          parameters = []
  
          if self.models:
              parameters += self.models.parameters
  
          return Parameters(parameters)
  
      @property
      def alpha(self):
          """Exposure ratio between signal and background regions"""
          alpha = self.acceptance / self.acceptance_off
          alpha.data = np.nan_to_num(alpha.data)
          return alpha
  
      @property
      def background(self):
          """Predicted background in the on region.
  
          Notice that this definition is valid under the assumption of cash statistic.
          """
          return self.alpha * self.counts_off
  
      @property
      def excess(self):
          """Excess (counts - alpha * counts_off)"""
          return self.counts.data - self.background.data
  
      def stat_array(self):
          """Likelihood per bin given the current model parameters"""
          mu_sig = self.npred().data
          on_stat_ = wstat(
              n_on=self.counts.data,
              n_off=self.counts_off.data,
              alpha=list(self.alpha.data),
              mu_sig=mu_sig,
          )
          return np.nan_to_num(on_stat_)
  
      @classmethod
      def from_geoms(
          cls,
          geom,
          geom_exposure,
          geom_psf,
          geom_edisp,
          reference_time="2000-01-01",
          name=None,
          **kwargs,
      ):
          """
          Create a MapDatasetOnOff object with zero filled maps according to the specified geometries
  
          Parameters
          ----------
          geom : `gammapy.maps.WcsGeom`
              geometry for the counts, counts_off, acceptance and acceptance_off maps
          geom_exposure : `gammapy.maps.WcsGeom`
              geometry for the exposure map
          geom_psf : `gammapy.maps.WcsGeom`
              geometry for the psf map
          geom_edisp : `gammapy.maps.WcsGeom`
              geometry for the energy dispersion map
          reference_time : `~astropy.time.Time`
              the reference time to use in GTI definition
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          empty_maps : `MapDatasetOnOff`
              A MapDatasetOnOff containing zero filled maps
          """
          kwargs = kwargs.copy()
          kwargs["name"] = name
  
          for key in ["counts", "counts_off", "acceptance", "acceptance_off"]:
              kwargs[key] = Map.from_geom(geom, unit="")
  
          kwargs["exposure"] = Map.from_geom(geom_exposure, unit="m2 s")
          kwargs["edisp"] = EDispMap.from_geom(geom_edisp)
          kwargs["psf"] = PSFMap.from_geom(geom_psf)
          kwargs["gti"] = GTI.create([] * u.s, [] * u.s, reference_time=reference_time)
          kwargs["mask_safe"] = Map.from_geom(geom, dtype=bool)
  
          return cls(**kwargs)
  
      @classmethod
      def from_map_dataset(
          cls, dataset, acceptance, acceptance_off, counts_off=None, name=None
      ):
          """Create spectrum dataseton off from another dataset.
  
          Parameters
          ----------
          dataset : `MapDataset`
              Spectrum dataset defining counts, edisp, aeff, livetime etc.
          acceptance : `Map`
              Relative background efficiency in the on region.
          acceptance_off : `Map`
              Relative background efficiency in the off region.
          counts_off : `Map`
              Off counts map . If the dataset provides a background model,
              and no off counts are defined. The off counts are deferred from
              counts_off / alpha.
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          dataset : `MapDatasetOnOff`
              Map dataset on off.
  
          """
          kwargs = {"name": name}
  
          if counts_off is None and dataset.background_model is not None:
              alpha = acceptance / acceptance_off
              kwargs["counts_off"] = dataset.background_model.evaluate() / alpha
  
          return cls(
              counts=dataset.counts,
              exposure=dataset.exposure,
              counts_off=counts_off,
              edisp=dataset.edisp,
              gti=dataset.gti,
              mask_safe=dataset.mask_safe,
              mask_fit=dataset.mask_fit,
              acceptance=acceptance,
              acceptance_off=acceptance_off,
              name=dataset.name,
              evaluation_mode=dataset.evaluation_mode,
          )
  
      @property
      def _is_stackable(self):
          """Check if the Dataset contains enough information to be stacked"""
          if (
              self.acceptance_off is None
              or self.acceptance is None
              or self.counts_off is None
          ):
              return False
          else:
              return True
  
      def stack(self, other):
          r"""Stack another dataset in place.
  
          The ``acceptance`` of the stacked dataset is normalized to 1,
          and the stacked ``acceptance_off`` is scaled so that:
  
          .. math::
              \alpha_\text{stacked} =
              \frac{1}{a_\text{off}} =
              \frac{\alpha_1\text{OFF}_1 + \alpha_2\text{OFF}_2}{\text{OFF}_1 + OFF_2}
  
          Parameters
          ----------
          other : `MapDatasetOnOff`
              Other dataset
          """
          if not isinstance(other, MapDatasetOnOff):
              raise TypeError("Incompatible types for MapDatasetOnOff stacking")
  
          if not self._is_stackable or not other._is_stackable:
              raise ValueError("Cannot stack incomplete MapDatsetOnOff.")
  
          # Factor containing: self.alpha * self.counts_off + other.alpha * other.counts_off
          tmp_factor = (self.alpha * self.counts_off).copy()
          tmp_factor.data[~self.mask_safe.data] = 0
          tmp_factor.stack(other.alpha * other.counts_off, weights=other.mask_safe)
  
          # Stack the off counts (in place)
          self.counts_off.data[~self.mask_safe.data] = 0
          self.counts_off.stack(other.counts_off, weights=other.mask_safe)
  
          self.acceptance_off = self.counts_off / tmp_factor
          self.acceptance.data = np.ones(self.data_shape)
  
          super().stack(other)
  
      def stat_sum(self):
          """Total likelihood given the current model parameters."""
          return Dataset.stat_sum(self)
  
      def fake(self, background_model, random_state="random-seed"):
          """Simulate fake counts (on and off) for the current model and reduced IRFs.
  
          This method overwrites the counts defined on the dataset object.
  
          Parameters
          ----------
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
                  Defines random number generator initialisation.
                  Passed to `~gammapy.utils.random.get_random_state`.
          """
          random_state = get_random_state(random_state)
          npred = self.npred()
          npred.data = random_state.poisson(npred.data)
  
          npred_bkg = background_model.copy()
          npred_bkg.data = random_state.poisson(npred_bkg.data)
  
          self.counts = npred + npred_bkg
  
          npred_off = background_model / self.alpha
          npred_off.data = random_state.poisson(npred_off.data)
          self.counts_off = npred_off
  
      def to_hdulist(self):
          """Convert map dataset to list of HDUs.
  
          Returns
          -------
          hdulist : `~astropy.io.fits.HDUList`
              Map dataset list of HDUs.
          """
          hdulist = super().to_hdulist()
          exclude_primary = slice(1, None)
  
          if self.counts_off is not None:
              hdulist += self.counts_off.to_hdulist(hdu="counts_off")[exclude_primary]
  
          if self.acceptance is not None:
              hdulist += self.acceptance.to_hdulist(hdu="acceptance")[exclude_primary]
  
          if self.acceptance_off is not None:
              hdulist += self.acceptance_off.to_hdulist(hdu="acceptance_off")[
                  exclude_primary
              ]
  
          return hdulist
  
      @classmethod
      def from_hdulist(cls, hdulist, name=None):
          """Create map dataset from list of HDUs.
  
          Parameters
          ----------
          hdulist : `~astropy.io.fits.HDUList`
              List of HDUs.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset.
          """
          kwargs = {}
          kwargs["name"] = name
          if "COUNTS" in hdulist:
              kwargs["counts"] = Map.from_hdulist(hdulist, hdu="counts")
  
          if "COUNTS_OFF" in hdulist:
              kwargs["counts_off"] = Map.from_hdulist(hdulist, hdu="counts_off")
  
          if "ACCEPTANCE" in hdulist:
              kwargs["acceptance"] = Map.from_hdulist(hdulist, hdu="acceptance")
  
          if "ACCEPTANCE_OFF" in hdulist:
              kwargs["acceptance_off"] = Map.from_hdulist(hdulist, hdu="acceptance_off")
  
          if "EXPOSURE" in hdulist:
              kwargs["exposure"] = Map.from_hdulist(hdulist, hdu="exposure")
  
          if "EDISP_MATRIX" in hdulist:
              kwargs["edisp"] = EDispKernel.from_hdulist(
                  hdulist, hdu1="EDISP_MATRIX", hdu2="EDISP_MATRIX_EBOUNDS"
              )
  
          if "PSF_KERNEL" in hdulist:
              psf_map = Map.from_hdulist(hdulist, hdu="psf_kernel")
              kwargs["psf"] = PSFKernel(psf_map)
  
          if "MASK_SAFE" in hdulist:
              mask_safe_map = Map.from_hdulist(hdulist, hdu="mask_safe")
              kwargs["mask_safe"] = mask_safe_map.data.astype(bool)
  
          if "MASK_FIT" in hdulist:
              mask_fit_map = Map.from_hdulist(hdulist, hdu="mask_fit")
              kwargs["mask_fit"] = mask_fit_map.data.astype(bool)
  
          if "GTI" in hdulist:
              gti = GTI(Table.read(hdulist, hdu="GTI"))
              kwargs["gti"] = gti
          return cls(**kwargs)
  
      def to_spectrum_dataset(self, on_region, containment_correction=False, name=None):
          """Return a ~gammapy.spectrum.SpectrumDatasetOnOff from on_region.
  
          Counts and OFF counts are summed in the on_region.
  
          Acceptance is the average of all acceptances while acceptance OFF
          is taken such that number of excess is preserved in the on_region.
  
          Effective area is taken from the average exposure divided by the livetime.
          Here we assume it is the sum of the GTIs.
  
          The energy dispersion kernel is obtained at the on_region center.
          Only regions with centers are supported.
  
          The model is not exported to the ~gammapy.spectrum.SpectrumDataset.
          It must be set after the dataset extraction.
  
          Parameters
          ----------
          on_region : `~regions.SkyRegion`
              the input ON region on which to extract the spectrum
          containment_correction : bool
              Apply containment correction for point sources and circular on regions
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `~gammapy.spectrum.SpectrumDatasetOnOff`
              the resulting reduced dataset
          """
          dataset = super().to_spectrum_dataset(on_region, containment_correction, name)
  
          kwargs = {}
          if self.counts_off is not None:
              kwargs["counts_off"] = self.counts_off.get_spectrum(on_region, np.sum)
  
          if self.acceptance is not None:
              kwargs["acceptance"] = self.acceptance.get_spectrum(on_region, np.mean)
              background = self.background.get_spectrum(on_region, np.sum)
              kwargs["acceptance_off"] = (
                  kwargs["acceptance"] * kwargs["counts_off"] / background
              )
  
          return SpectrumDatasetOnOff.from_spectrum_dataset(dataset=dataset, **kwargs)
  
      def cutout(self, position, width, mode="trim", name=None):
          """Cutout map dataset.
  
          Parameters
          ----------
          position : `~astropy.coordinates.SkyCoord`
              Center position of the cutout region.
          width : tuple of `~astropy.coordinates.Angle`
              Angular sizes of the region in (lon, lat) in that specific order.
              If only one value is passed, a square region is extracted.
          mode : {'trim', 'partial', 'strict'}
              Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          cutout : `MapDatasetOnOff`
              Cutout map dataset.
          """
          cutout_kwargs = {
              "position": position,
              "width": width,
              "mode": mode,
              "name": name,
          }
  
          cutout_dataset = super().cutout(**cutout_kwargs)
  
          del cutout_kwargs["name"]
  
          if self.counts_off is not None:
              cutout_dataset.counts_off = self.counts_off.cutout(**cutout_kwargs)
  
          if self.acceptance is not None:
              cutout_dataset.acceptance = self.acceptance.cutout(**cutout_kwargs)
  
          if self.acceptance_off is not None:
              cutout_dataset.acceptance_off = self.acceptance_off.cutout(**cutout_kwargs)
  
          return cutout_dataset
  
      def to_image(self, spectrum=None, name=None):
          """Create images by summing over the energy axis.
  
          Exposure is weighted with an assumed spectrum,
          resulting in a weighted mean exposure image.
  
          Currently the PSFMap and EdispMap are dropped from the
          resulting image dataset.
  
          Parameters
          ----------
          spectrum : `~gammapy.modeling.models.SpectralModel`
              Spectral model to compute the weights.
              Default is power-law with spectral index of 2.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDatasetOnOff`
              Map dataset containing images.
          """
          kwargs = {"name": name}
          dataset = super().to_image(spectrum, name)
  
          if self.mask_safe is not None:
              mask_safe = self.mask_safe
          else:
              mask_safe = 1
  
          if self.counts_off is not None:
              counts_off = self.counts_off * mask_safe
              kwargs["counts_off"] = counts_off.sum_over_axes(keepdims=True)
  
          if self.acceptance is not None:
              acceptance = self.acceptance * mask_safe
              kwargs["acceptance"] = acceptance.sum_over_axes(keepdims=True)
  
              background = self.background * mask_safe
              background = background.sum_over_axes(keepdims=True)
              kwargs["acceptance_off"] = (
                  kwargs["acceptance"] * kwargs["counts_off"] / background
              )
  
          return self.from_map_dataset(dataset, **kwargs)
  
  
  class MapEvaluator:
      """Sky model evaluation on maps.
  
      This evaluates a sky model on a 3D map and convolves with the IRFs,
      and returns a map of the predicted counts.
      Note that background counts are not added.
  
      For now, we only make it work for 3D WCS maps with an energy axis.
      No HPX, no other axes, those can be added later here or via new
      separate model evaluator classes.
  
      Parameters
      ----------
      model : `~gammapy.modeling.models.SkyModel`
          Sky model
      exposure : `~gammapy.maps.Map`
          Exposure map
      psf : `~gammapy.cube.PSFKernel`
          PSF kernel
      edisp : `~gammapy.irf.EDispKernel`
          Energy dispersion
      evaluation_mode : {"local", "global"}
          Model evaluation mode.
      """
  
      def __init__(
          self, model=None, exposure=None, psf=None, edisp=None, evaluation_mode="local"
      ):
          self.model = model
          self.exposure = exposure
          self.psf = psf
          self.edisp = edisp
          self.contributes = True
  
          if evaluation_mode not in {"local", "global"}:
              raise ValueError(f"Invalid evaluation_mode: {evaluation_mode!r}")
  
          self.evaluation_mode = evaluation_mode
  
          # TODO: this is preliminary solution until we have further unified the model handling
          if isinstance(model, BackgroundModel):
              self.compute_npred = model.evaluate
              self.evaluation_mode = "global"
  
      @property
      def geom(self):
          """True energy map geometry (`~gammapy.maps.Geom`)"""
          return self.exposure.geom
  
      @property
      def needs_update(self):
          """Check whether the model component has drifted away from its support."""
          if self.exposure is None:
              return True
          elif self.evaluation_mode == "global" or self.model.evaluation_radius is None:
              return False
          else:
              position = self.model.position
              separation = self._init_position.separation(position)
              update = separation > (self.model.evaluation_radius + CUTOUT_MARGIN)
          return update
  
      def update(self, exposure, psf, edisp, geom):
          """Update MapEvaluator, based on the current position of the model component.
  
          Parameters
          ----------
          exposure : `~gammapy.maps.Map`
              Exposure map.
          psf : `gammapy.cube.PSFMap`
              PSF map.
          edisp : `gammapy.cube.EDispMap`
              Edisp map.
          geom : `WcsGeom`
              Counts geom
          """
          log.debug("Updating model evaluator")
          # cache current position of the model component
  
          if isinstance(edisp, EDispMap):
              e_reco = geom.get_axis_by_name("energy").edges
              self.edisp = edisp.get_edisp_kernel(self.model.position, e_reco=e_reco)
          else:
              self.edisp = edisp
  
          if isinstance(psf, PSFMap):
              self.psf = psf.get_psf_kernel(self.model.position, geom=exposure.geom)
          else:
              self.psf = psf
  
          if self.evaluation_mode == "local" and self.model.evaluation_radius is not None:
              self._init_position = self.model.position
              if self.psf is not None:
                  psf_width = np.max(self.psf.psf_kernel_map.geom.width)
              else:
                  psf_width = 0 * u.deg
  
              width = psf_width + 2 * (self.model.evaluation_radius + CUTOUT_MARGIN)
              try:
                  self.exposure = exposure.cutout(
                      position=self.model.position, width=width
                  )
                  self.contributes = True
              except (NoOverlapError, ValueError):
                  self.contributes = False
          else:
              self.exposure = exposure
  
      def compute_dnde(self):
          """Compute model differential flux at map pixel centers.
  
          Returns
          -------
          model_map : `~gammapy.maps.Map`
              Sky cube with data filled with evaluated model values.
              Units: ``cm-2 s-1 TeV-1 deg-2``
          """
          return self.model.evaluate_geom(self.geom)
  
      def compute_flux(self):
          """Compute model integral flux over map pixel volumes.
  
          For now, we simply multiply dnde with bin volume.
          """
          dnde = self.compute_dnde()
          volume = self.geom.bin_volume()
          return dnde * volume
  
      def apply_exposure(self, flux):
          """Compute npred cube
  
          For now just divide flux cube by exposure
          """
          npred = (flux * self.exposure.quantity).to_value("")
          return Map.from_geom(self.geom, data=npred, unit="")
  
      def apply_psf(self, npred):
          """Convolve npred cube with PSF"""
          tmp = npred.convolve(self.psf)
          tmp.data[tmp.data < 0.0] = 0
          return tmp
  
      def apply_edisp(self, npred):
          """Convolve map data with energy dispersion.
  
          Parameters
          ----------
          npred : `~gammapy.maps.Map`
              Predicted counts in true energy bins
  
          Returns
          -------
          npred_reco : `~gammapy.maps.Map`
              Predicted counts in reco energy bins
          """
          return npred.apply_edisp(self.edisp)
  
      def compute_npred(self):
          """
          Evaluate model predicted counts.
  
          Returns
          -------
          npred : `~gammapy.maps.Map`
              Predicted counts on the map (in reco energy bins)
          """
          flux = self.compute_flux()
          npred = self.apply_exposure(flux)
          if self.psf is not None and self.model.processing["psf"] == True:
              npred = self.apply_psf(npred)
          if self.edisp is not None and self.model.processing["edisp"] == True:
              npred = self.apply_edisp(npred)
  
          return npred
  
files: gammapy/cube/fit.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import logging
  import numpy as np
  import astropy.units as u
  from astropy.io import fits
  from astropy.nddata.utils import NoOverlapError
  from astropy.table import Table
  from astropy.utils import lazyproperty
  from regions import CircleSkyRegion
  from gammapy.cube.edisp_map import EDispMap
  from gammapy.cube.psf_kernel import PSFKernel
  from gammapy.cube.psf_map import PSFMap
  from gammapy.data import GTI
  from gammapy.irf import EDispKernel, EffectiveAreaTable
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling import Dataset, Parameters
  from gammapy.modeling.models import BackgroundModel, Models
  from gammapy.spectrum import SpectrumDataset, SpectrumDatasetOnOff
  from gammapy.stats import cash, cash_sum_cython, wstat
  from gammapy.utils.random import get_random_state
  from gammapy.utils.scripts import make_name, make_path
  from .exposure import _map_spectrum_weight
  
  __all__ = ["MapDataset", "MapDatasetOnOff"]
  
  log = logging.getLogger(__name__)
  
  CUTOUT_MARGIN = 0.1 * u.deg
  RAD_MAX = 0.66
  RAD_AXIS_DEFAULT = MapAxis.from_bounds(
      0, RAD_MAX, nbin=66, node_type="edges", name="theta", unit="deg"
  )
  MIGRA_AXIS_DEFAULT = MapAxis.from_bounds(
      0.2, 5, nbin=48, node_type="edges", name="migra"
  )
  
  BINSZ_IRF_DEFAULT = 0.2
  
  
  class MapDataset(Dataset):
      """Perform sky model likelihood fit on maps.
  
      Parameters
      ----------
      models : `~gammapy.modeling.models.Models`
          Source sky models.
      counts : `~gammapy.maps.WcsNDMap`
          Counts cube
      exposure : `~gammapy.maps.WcsNDMap`
          Exposure cube
      mask_fit : `~gammapy.maps.WcsNDMap`
          Mask to apply to the likelihood for fitting.
      psf : `~gammapy.cube.PSFKernel` or `~gammapy.cube.PSFMap`
          PSF kernel
      edisp : `~gammapy.irf.EDispKernel` or `~gammapy.cube.EDispMap`
          Energy dispersion kernel
      evaluation_mode : {"local", "global"}
          Model evaluation mode.
          The "local" mode evaluates the model components on smaller grids to save computation time.
          This mode is recommended for local optimization algorithms.
          The "global" evaluation mode evaluates the model components on the full map.
          This mode is recommended for global optimization algorithms.
      mask_safe : `~gammapy.maps.WcsNDMap`
          Mask defining the safe data range.
      gti : `~gammapy.data.GTI`
          GTI of the observation or union of GTI if it is a stacked observation
      """
  
      stat_type = "cash"
      tag = "MapDataset"
  
      def __init__(
          self,
          models=None,
          counts=None,
          exposure=None,
          mask_fit=None,
          psf=None,
          edisp=None,
          name=None,
          evaluation_mode="local",
          mask_safe=None,
          gti=None,
      ):
          if mask_fit is not None and mask_fit.data.dtype != np.dtype("bool"):
              raise ValueError("mask data must have dtype bool")
  
          if mask_safe is not None and mask_safe.data.dtype != np.dtype("bool"):
              raise ValueError("mask data must have dtype bool")
  
          self._name = make_name(name)
          self._evaluators = {}
  
          self.background_model = None
          self.evaluation_mode = evaluation_mode
          self.counts = counts
          self.exposure = exposure
          self.mask_fit = mask_fit
          self.psf = psf
          self.edisp = edisp
          self.mask_safe = mask_safe
          self.models = models
          self.gti = gti
  
          # check whether a reference geom is defined
          _ = self._geom
  
      @property
      def name(self):
          return self._name
  
      def __str__(self):
          str_ = f"{self.__class__.__name__}\n"
          str_ += "-" * len(self.__class__.__name__) + "\n"
          str_ += "\n"
  
          str_ += "\t{:32}: {} \n\n".format("Name", self.name)
  
          counts = np.nan
          if self.counts is not None:
              counts = np.sum(self.counts.data)
          str_ += "\t{:32}: {:.0f} \n".format("Total counts", counts)
  
          npred = np.nan
          if self.models is not None:
              npred = np.sum(self.npred().data)
          str_ += "\t{:32}: {:.2f}\n".format("Total predicted counts", npred)
  
          background = np.nan
          if self.background_model is not None:
              background = np.sum(self.background_model.evaluate().data)
          str_ += "\t{:32}: {:.2f}\n\n".format("Total background counts", background)
  
          exposure_min, exposure_max, exposure_unit = np.nan, np.nan, ""
          if self.exposure is not None:
              if self.mask_safe is not None:
                  mask = self.mask_safe.reduce_over_axes(np.logical_or).data
                  if not mask.any():
                      mask = None
              else:
                  mask = None
              exposure_min = np.min(self.exposure.data[..., mask])
              exposure_max = np.max(self.exposure.data[..., mask])
              exposure_unit = self.exposure.unit
  
          str_ += "\t{:32}: {:.2e} {}\n".format(
              "Exposure min", exposure_min, exposure_unit
          )
          str_ += "\t{:32}: {:.2e} {}\n\n".format(
              "Exposure max", exposure_max, exposure_unit
          )
  
          # data section
          n_bins = 0
          if self.counts is not None:
              n_bins = self.counts.data.size
          str_ += "\t{:32}: {} \n".format("Number of total bins", n_bins)
  
          n_fit_bins = 0
          if self.mask is not None:
              n_fit_bins = np.sum(self.mask.data)
          str_ += "\t{:32}: {} \n\n".format("Number of fit bins", n_fit_bins)
  
          # likelihood section
          str_ += "\t{:32}: {}\n".format("Fit statistic type", self.stat_type)
  
          stat = np.nan
          if self.counts is not None and self.models is not None:
              stat = self.stat_sum()
          str_ += "\t{:32}: {:.2f}\n\n".format("Fit statistic value (-2 log(L))", stat)
  
          # model section
          n_models = 0
          if self.models is not None:
              n_models = len(self.models)
  
          str_ += "\t{:32}: {} \n".format("Number of models", n_models)
          str_ += "\t{:32}: {}\n".format("Number of parameters", len(self.parameters))
          str_ += "\t{:32}: {}\n\n".format(
              "Number of free parameters", len(self.parameters.free_parameters)
          )
  
          if self.models is not None:
              str_ += "\t" + "\n\t".join(str(self.models).split("\n")[2:])
  
          return str_.expandtabs(tabsize=2)
  
      @property
      def models(self):
          """Models (`~gammapy.modeling.models.Models`)."""
          return self._models
  
      @models.setter
      def models(self, value):
          if value is not None:
              self._models = Models(value)
          else:
              self._models = None
  
          if self.models is not None:
              for model in self.models:
                  if isinstance(model, BackgroundModel):
                      self.background_model = model
                      break
              else:
                  log.warning(f"No background model defined for dataset {self.name}")
  
      @property
      def evaluators(self):
          """Model evaluators"""
          # this call is needed to trigger the setup of the evaluators
          if not self._evaluators:
              self.npred()
  
          return self._evaluators
  
      @property
      def parameters(self):
          """List of parameters (`~gammapy.modeling.Parameters`)"""
          return self.models.parameters
  
      @property
      def _geom(self):
          """Main analysis geometry"""
          if self.counts is not None:
              return self.counts.geom
          elif self.background_model is not None:
              return self.background_model.map.geom
          elif self.mask_safe is not None:
              return self.mask_safe.geom
          elif self.mask_fit is not None:
              return self.mask_fit.geom
          else:
              raise ValueError(
                  "Either 'counts', 'background_model', 'mask_fit'"
                  " or 'mask_safe' must be defined."
              )
  
      @property
      def data_shape(self):
          """Shape of the counts or background data (tuple)"""
          return self._geom.data_shape
  
      def npred(self):
          """Predicted source and background counts (`~gammapy.maps.Map`)."""
          npred_total = Map.from_geom(self._geom, dtype=float)
  
          if self.models:
              for model in self.models:
                  evaluator = self._evaluators.get(model.name)
  
                  if evaluator is None:
                      evaluator = MapEvaluator(model=model, evaluation_mode=self.evaluation_mode)
                      self._evaluators[model.name] = evaluator
  
                  # if the model component drifts out of its support the evaluator has
                  # has to be updated
                  if evaluator.needs_update:
                      evaluator.update(self.exposure, self.psf, self.edisp, self._geom)
  
                  if evaluator.contributes:
                      npred = evaluator.compute_npred()
                      npred_total.stack(npred)
  
          return npred_total
  
      @classmethod
      def from_geoms(
          cls,
          geom,
          geom_exposure,
          geom_psf,
          geom_edisp,
          reference_time="2000-01-01",
          name=None,
          **kwargs,
      ):
          """
          Create a MapDataset object with zero filled maps according to the specified geometries
  
          Parameters
          ----------
          geom : `Geom`
              geometry for the counts and background maps
          geom_exposure : `Geom`
              geometry for the exposure map
          geom_psf : `Geom`
              geometry for the psf map
          geom_edisp : `Geom`
              geometry for the energy dispersion map
          reference_time : `~astropy.time.Time`
              the reference time to use in GTI definition
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          empty_maps : `MapDataset`
              A MapDataset containing zero filled maps
          """
          kwargs = kwargs.copy()
          kwargs["name"] = name
          kwargs["counts"] = Map.from_geom(geom, unit="")
  
          background = Map.from_geom(geom, unit="")
          kwargs["models"] = Models([BackgroundModel(background)])
          kwargs["exposure"] = Map.from_geom(geom_exposure, unit="m2 s")
          kwargs["edisp"] = EDispMap.from_geom(geom_edisp)
          kwargs["psf"] = PSFMap.from_geom(geom_psf)
          kwargs["gti"] = GTI.create([] * u.s, [] * u.s, reference_time=reference_time)
          kwargs["mask_safe"] = Map.from_geom(geom, unit="", dtype=bool)
  
          return cls(**kwargs)
  
      @classmethod
      def create(
          cls,
          geom,
          energy_axis_true=None,
          migra_axis=None,
          rad_axis=None,
          binsz_irf=None,
          reference_time="2000-01-01",
          name=None,
          **kwargs,
      ):
          """Create a MapDataset object with zero filled maps.
  
          Parameters
          ----------
          geom : `~gammapy.maps.WcsGeom`
              Reference target geometry in reco energy, used for counts and background maps
          energy_axis_true : `~gammapy.maps.MapAxis`
              True energy axis used for IRF maps
          migra_axis : `~gammapy.maps.MapAxis`
              Migration axis for the energy dispersion map
          rad_axis : `~gammapy.maps.MapAxis`
              Rad axis for the psf map
          binsz_irf : float
              IRF Map pixel size in degrees.
          reference_time : `~astropy.time.Time`
              the reference time to use in GTI definition
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          empty_maps : `MapDataset`
              A MapDataset containing zero filled maps
          """
          migra_axis = migra_axis or MIGRA_AXIS_DEFAULT
          rad_axis = rad_axis or RAD_AXIS_DEFAULT
          energy_axis_true = energy_axis_true or geom.get_axis_by_name("energy")
          binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT
  
          geom_image = geom.to_image()
          geom_exposure = geom_image.to_cube([energy_axis_true])
          geom_irf = geom_image.to_binsz(binsz=binsz_irf)
          geom_psf = geom_irf.to_cube([rad_axis, energy_axis_true])
          geom_edisp = geom_irf.to_cube([migra_axis, energy_axis_true])
  
          return cls.from_geoms(
              geom,
              geom_exposure,
              geom_psf,
              geom_edisp,
              reference_time=reference_time,
              name=name,
              **kwargs,
          )
  
      def stack(self, other):
          """Stack another dataset in place.
  
          Parameters
          ----------
          other: `~gammapy.cube.MapDataset`
              Map dataset to be stacked with this one.
          """
  
          if self.counts and other.counts:
              self.counts *= self.mask_safe
              self.counts.stack(other.counts, weights=other.mask_safe)
  
          if self.exposure and other.exposure:
              mask_image = self.mask_safe.reduce_over_axes(func=np.logical_or)
              self.exposure *= mask_image.data
              # TODO: apply energy dependent mask to exposure. Does this require
              #  a mask_safe in true energy?
              mask_image_other = other.mask_safe.reduce_over_axes(func=np.logical_or)
              self.exposure.stack(other.exposure, weights=mask_image_other)
  
          if self.background_model and other.background_model:
              self.background_model.map *= self.mask_safe
              self.background_model.stack(other.background_model, other.mask_safe)
  
          if self.mask_safe is not None and other.mask_safe is not None:
              self.mask_safe.stack(other.mask_safe)
  
          if self.psf and other.psf:
              if isinstance(self.psf, PSFMap) and isinstance(other.psf, PSFMap):
                  mask_irf = self._mask_safe_irf(self.psf.psf_map, mask_image)
                  self.psf.psf_map *= mask_irf.data
                  self.psf.exposure_map *= mask_irf.data
  
                  mask_image_other = other.mask_safe.reduce_over_axes(func=np.logical_or)
                  mask_irf_other = self._mask_safe_irf(
                      other.psf.psf_map, mask_image_other
                  )
                  self.psf.stack(other.psf, weights=mask_irf_other)
              else:
                  raise ValueError("Stacking of PSF kernels not supported")
  
          if self.edisp and other.edisp:
              if isinstance(self.edisp, EDispMap) and isinstance(other.edisp, EDispMap):
                  mask_irf = self._mask_safe_irf(self.edisp.edisp_map, mask_image)
                  self.edisp.edisp_map *= mask_irf.data
                  self.edisp.exposure_map *= mask_irf.data
  
                  mask_image_other = other.mask_safe.reduce_over_axes(func=np.logical_or)
                  mask_irf_other = self._mask_safe_irf(
                      other.edisp.edisp_map, mask_image_other
                  )
                  self.edisp.stack(other.edisp, weights=mask_irf_other)
              else:
                  raise ValueError("Stacking of edisp kernels not supported")
  
          if self.gti and other.gti:
              self.gti = self.gti.stack(other.gti).union()
  
      @staticmethod
      def _mask_safe_irf(irf_map, mask):
          geom = irf_map.geom.to_image()
          coords = geom.get_coord()
          data = mask.get_by_coord(coords).astype(bool)
          return Map.from_geom(geom=geom, data=data)
  
      def stat_array(self):
          """Likelihood per bin given the current model parameters"""
          return cash(n_on=self.counts.data, mu_on=self.npred().data)
  
      def residuals(self, method="diff"):
          """Compute residuals map.
  
          Parameters
          ----------
          method: {"diff", "diff/model", "diff/sqrt(model)"}
              Method used to compute the residuals. Available options are:
                  - "diff" (default): data - model
                  - "diff/model": (data - model) / model
                  - "diff/sqrt(model)": (data - model) / sqrt(model)
  
          Returns
          -------
          residuals : `gammapy.maps.WcsNDMap`
              Residual map.
          """
          return self._compute_residuals(self.counts, self.npred(), method=method)
  
      def plot_residuals(
          self,
          method="diff",
          smooth_kernel="gauss",
          smooth_radius="0.1 deg",
          region=None,
          figsize=(12, 4),
          **kwargs,
      ):
          """
          Plot spatial and spectral residuals.
  
          The spectral residuals are extracted from the provided region, and the
          normalization used for the residuals computation can be controlled using
          the method parameter. If no region is passed, only the spatial
          residuals are shown.
  
          Parameters
          ----------
          method : {"diff", "diff/model", "diff/sqrt(model)"}
              Method used to compute the residuals, see `MapDataset.residuals()`
          smooth_kernel : {'gauss', 'box'}
              Kernel shape.
          smooth_radius: `~astropy.units.Quantity`, str or float
              Smoothing width given as quantity or float. If a float is given it
              is interpreted as smoothing width in pixels.
          region: `~regions.Region`
              Region (pixel or sky regions accepted)
          figsize : tuple
              Figure size used for the plotting.
          **kwargs : dict
              Keyword arguments passed to `~matplotlib.pyplot.imshow`.
  
          Returns
          -------
          ax_image, ax_spec : `~matplotlib.pyplot.Axes`,
              Image and spectrum axes.
          """
          import matplotlib.pyplot as plt
  
          fig = plt.figure(figsize=figsize)
  
          counts, npred = self.counts, self.npred()
  
          if self.mask is not None:
              counts = counts * self.mask
              npred = npred * self.mask
  
          counts_spatial = counts.sum_over_axes().smooth(
              width=smooth_radius, kernel=smooth_kernel
          )
          npred_spatial = npred.sum_over_axes().smooth(
              width=smooth_radius, kernel=smooth_kernel
          )
          spatial_residuals = self._compute_residuals(
              counts_spatial, npred_spatial, method
          )
  
          if self.mask_safe is not None:
              mask = self.mask_safe.reduce_over_axes(func=np.logical_or)
              spatial_residuals.data[~mask.data] = np.nan
  
          # If no region is provided, skip spectral residuals
          ncols = 2 if region is not None else 1
          ax_image = fig.add_subplot(1, ncols, 1, projection=spatial_residuals.geom.wcs)
          ax_spec = None
  
          kwargs.setdefault("cmap", "coolwarm")
          kwargs.setdefault("stretch", "linear")
          kwargs.setdefault("vmin", -5)
          kwargs.setdefault("vmax", 5)
          spatial_residuals.plot(ax=ax_image, add_cbar=True, **kwargs)
  
          # Spectral residuals
          if region:
              ax_spec = fig.add_subplot(1, 2, 2)
              counts_spec = counts.get_spectrum(region=region)
              npred_spec = npred.get_spectrum(region=region)
              residuals = self._compute_residuals(counts_spec, npred_spec, method)
              ax = residuals.plot()
              ax.axhline(0, color="black", lw=0.5)
  
              y_max = 2 * np.nanmax(residuals.data)
              plt.ylim(-y_max, y_max)
              label = self._residuals_labels[method]
              plt.ylabel(f"Residuals ({label})")
  
              # Overlay spectral extraction region on the spatial residuals
              pix_region = region.to_pixel(wcs=spatial_residuals.geom.wcs)
              pix_region.plot(ax=ax_image)
  
          return ax_image, ax_spec
  
      @lazyproperty
      def _counts_data(self):
          return self.counts.data.astype(float)
  
      def stat_sum(self):
          """Total likelihood given the current model parameters."""
          counts, npred = self._counts_data, self.npred().data
  
          if self.mask is not None:
              return cash_sum_cython(counts[self.mask.data], npred[self.mask.data])
          else:
              return cash_sum_cython(counts.ravel(), npred.ravel())
  
      def fake(self, random_state="random-seed"):
          """Simulate fake counts for the current model and reduced IRFs.
  
          This method overwrites the counts defined on the dataset object.
  
          Parameters
          ----------
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
                  Defines random number generator initialisation.
                  Passed to `~gammapy.utils.random.get_random_state`.
          """
          random_state = get_random_state(random_state)
          npred = self.npred()
          npred.data = random_state.poisson(npred.data)
          self.counts = npred
  
      def to_hdulist(self):
          """Convert map dataset to list of HDUs.
  
          Returns
          -------
          hdulist : `~astropy.io.fits.HDUList`
              Map dataset list of HDUs.
          """
          # TODO: what todo about the model and background model parameters?
          exclude_primary = slice(1, None)
  
          hdu_primary = fits.PrimaryHDU()
          hdulist = fits.HDUList([hdu_primary])
          if self.counts is not None:
              hdulist += self.counts.to_hdulist(hdu="counts")[exclude_primary]
  
          if self.exposure is not None:
              hdulist += self.exposure.to_hdulist(hdu="exposure")[exclude_primary]
  
          if self.background_model is not None:
              hdulist += self.background_model.map.to_hdulist(hdu="background")[
                  exclude_primary
              ]
  
          if self.edisp is not None:
              if isinstance(self.edisp, EDispKernel):
                  hdus = self.edisp.to_hdulist()
                  hdus["MATRIX"].name = "edisp_matrix"
                  hdus["EBOUNDS"].name = "edisp_matrix_ebounds"
                  hdulist.append(hdus["EDISP_MATRIX"])
                  hdulist.append(hdus["EDISP_MATRIX_EBOUNDS"])
              else:
                  hdulist += self.edisp.edisp_map.to_hdulist(hdu="EDISP")[exclude_primary]
                  hdulist += self.edisp.exposure_map.to_hdulist(hdu="edisp_exposure")[
                      exclude_primary
                  ]
  
          if self.psf is not None:
              if isinstance(self.psf, PSFKernel):
                  hdulist += self.psf.psf_kernel_map.to_hdulist(hdu="psf_kernel")[
                      exclude_primary
                  ]
              else:
                  hdulist += self.psf.psf_map.to_hdulist(hdu="psf")[exclude_primary]
                  hdulist += self.psf.exposure_map.to_hdulist(hdu="psf_exposure")[
                      exclude_primary
                  ]
  
          if self.mask_safe is not None:
              mask_safe_int = self.mask_safe.copy()
              mask_safe_int.data = mask_safe_int.data.astype(int)
              hdulist += mask_safe_int.to_hdulist(hdu="mask_safe")[exclude_primary]
  
          if self.mask_fit is not None:
              mask_fit_int = self.mask_fit.copy()
              mask_fit_int.data = mask_fit_int.data.astype(int)
              hdulist += mask_fit_int.to_hdulist(hdu="mask_fit")[exclude_primary]
  
          if self.gti is not None:
              hdulist.append(fits.BinTableHDU(self.gti.table, name="GTI"))
  
          return hdulist
  
      @classmethod
      def from_hdulist(cls, hdulist, name=None):
          """Create map dataset from list of HDUs.
  
          Parameters
          ----------
          hdulist : `~astropy.io.fits.HDUList`
              List of HDUs.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset.
          """
          kwargs = {"name": name}
  
          if "COUNTS" in hdulist:
              kwargs["counts"] = Map.from_hdulist(hdulist, hdu="counts")
  
          if "EXPOSURE" in hdulist:
              kwargs["exposure"] = Map.from_hdulist(hdulist, hdu="exposure")
  
          if "BACKGROUND" in hdulist:
              background_map = Map.from_hdulist(hdulist, hdu="background")
              kwargs["models"] = Models([BackgroundModel(background_map)])
  
          if "EDISP_MATRIX" in hdulist:
              kwargs["edisp"] = EDispKernel.from_hdulist(
                  hdulist, hdu1="EDISP_MATRIX", hdu2="EDISP_MATRIX_EBOUNDS"
              )
          if "EDISP" in hdulist:
              edisp_map = Map.from_hdulist(hdulist, hdu="edisp")
              exposure_map = Map.from_hdulist(hdulist, hdu="edisp_exposure")
              kwargs["edisp"] = EDispMap(edisp_map, exposure_map)
  
          if "PSF_KERNEL" in hdulist:
              psf_map = Map.from_hdulist(hdulist, hdu="psf_kernel")
              kwargs["psf"] = PSFKernel(psf_map)
          if "PSF" in hdulist:
              psf_map = Map.from_hdulist(hdulist, hdu="psf")
              exposure_map = Map.from_hdulist(hdulist, hdu="psf_exposure")
              kwargs["psf"] = PSFMap(psf_map, exposure_map)
  
          if "MASK_SAFE" in hdulist:
              mask_safe = Map.from_hdulist(hdulist, hdu="mask_safe")
              mask_safe.data = mask_safe.data.astype(bool)
              kwargs["mask_safe"] = mask_safe
  
          if "MASK_FIT" in hdulist:
              mask_fit = Map.from_hdulist(hdulist, hdu="mask_fit")
              mask_fit.data = mask_fit.data.astype(bool)
              kwargs["mask_fit"] = mask_fit
  
          if "GTI" in hdulist:
              gti = GTI(Table.read(hdulist, hdu="GTI"))
              kwargs["gti"] = gti
  
          return cls(**kwargs)
  
      def write(self, filename, overwrite=False):
          """Write map dataset to file.
  
          Parameters
          ----------
          filename : str
              Filename to write to.
          overwrite : bool
              Overwrite file if it exists.
          """
          self.to_hdulist().writeto(make_path(filename), overwrite=overwrite)
  
      @classmethod
      def read(cls, filename, name=None):
          """Read map dataset from file.
  
          Parameters
          ----------
          filename : str
              Filename to read from.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset.
          """
          with fits.open(make_path(filename), memmap=False) as hdulist:
              return cls.from_hdulist(hdulist, name=name)
  
      @classmethod
      def from_dict(cls, data, components, models):
          """Create from dicts and models list generated from YAML serialization."""
          dataset = cls.read(data["filename"], name=data["name"])
          bkg_name = data["background"]
          model_names = data["models"]
          models_list = [model for model in models if model.name in model_names]
          models = Models(models_list)
  
          for component in components["components"]:
              if component["type"] == "BackgroundModel":
                  if component["name"] == bkg_name:
                      if "filename" not in component:
                          component["map"] = dataset.background_model.map
                      background_model = BackgroundModel.from_dict(component)
                      models.append(background_model)
  
          dataset.models = models
          return dataset
  
      def to_dict(self, filename=""):
          """Convert to dict for YAML serialization."""
          if self.models is None:
              models = []
          else:
              models = [_.name for _ in self.models]
  
          return {
              "name": self.name,
              "type": self.tag,
              "models": models,
              "background": self.background_model.name,
              "filename": str(filename),
          }
  
      def to_spectrum_dataset(self, on_region, containment_correction=False, name=None):
          """Return a ~gammapy.spectrum.SpectrumDataset from on_region.
  
          Counts and background are summed in the on_region.
  
          Effective area is taken from the average exposure divided by the livetime.
          Here we assume it is the sum of the GTIs.
  
          The energy dispersion kernel is obtained at the on_region center.
          Only regions with centers are supported.
  
          The model is not exported to the ~gammapy.spectrum.SpectrumDataset.
          It must be set after the dataset extraction.
  
          Parameters
          ----------
          on_region : `~regions.SkyRegion`
              the input ON region on which to extract the spectrum
          containment_correction : bool
              Apply containment correction for point sources and circular on regions
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `~gammapy.spectrum.SpectrumDataset`
              the resulting reduced dataset
          """
          kwargs = {"gti": self.gti, "name": name}
  
          if self.gti is not None:
              kwargs["livetime"] = self.gti.time_sum
          else:
              raise ValueError("No GTI in `MapDataset`, cannot compute livetime")
  
          if self.counts is not None:
              kwargs["counts"] = self.counts.get_spectrum(on_region, np.sum)
  
          if self.background_model is not None:
              kwargs["background"] = self.background_model.evaluate().get_spectrum(
                  on_region, np.sum
              )
  
          if self.exposure is not None:
              exposure = self.exposure.get_spectrum(on_region, np.mean)
              kwargs["aeff"] = EffectiveAreaTable(
                  energy_lo=exposure.energy.edges[:-1],
                  energy_hi=exposure.energy.edges[1:],
                  data=exposure.quantity / kwargs["livetime"],
              )
  
          if containment_correction:
              if not isinstance(on_region, CircleSkyRegion):
                  raise TypeError(
                      "Containement correction is only supported for"
                      " `CircleSkyRegion`."
                  )
              elif self.psf is None or isinstance(self.psf, PSFKernel):
                  raise ValueError("No PSFMap set. Containement correction impossible")
              else:
                  psf = self.psf.get_energy_dependent_table_psf(on_region.center)
                  containment = psf.containment(
                      kwargs["aeff"].energy.center, on_region.radius
                  )
                  kwargs["aeff"].data.data *= containment.squeeze()
  
          if self.edisp is not None:
              if isinstance(self.edisp, EDispKernel):
                  edisp = self.edisp
              else:
                  axis = self._geom.get_axis_by_name("energy")
                  edisp = self.edisp.get_edisp_kernel(on_region.center, e_reco=axis.edges)
              kwargs["edisp"] = edisp
  
          return SpectrumDataset(**kwargs)
  
      def to_image(self, spectrum=None, name=None):
          """Create images by summing over the energy axis.
  
          Exposure is weighted with an assumed spectrum,
          resulting in a weighted mean exposure image.
  
          Currently the PSFMap and EdispMap are dropped from the
          resulting image dataset.
  
          Parameters
          ----------
          spectrum : `~gammapy.modeling.models.SpectralModel`
              Spectral model to compute the weights.
              Default is power-law with spectral index of 2.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset containing images.
          """
          kwargs = {}
          kwargs["name"] = make_name(name)
          kwargs["gti"] = self.gti
  
          if self.mask_safe is not None:
              mask_safe = self.mask_safe
              kwargs["mask_safe"] = mask_safe.reduce_over_axes(
                  func=np.logical_or, keepdims=True
              )
          else:
              mask_safe = 1
  
          if self.counts is not None:
              counts = self.counts * mask_safe
              kwargs["counts"] = counts.sum_over_axes(keepdims=True)
  
          if self.exposure is not None:
              exposure = _map_spectrum_weight(self.exposure, spectrum)
              kwargs["exposure"] = exposure.sum_over_axes(keepdims=True)
  
          if self.background_model is not None:
              background = self.background_model.evaluate() * mask_safe
              background = background.sum_over_axes(keepdims=True)
              kwargs["models"] = Models([BackgroundModel(background)])
  
          if self.psf is not None:
              kwargs["psf"] = self.psf.to_image(spectrum=spectrum, keepdims=True)
  
          return self.__class__(**kwargs)
  
      def cutout(self, position, width, mode="trim", name=None):
          """Cutout map dataset.
  
          Parameters
          ----------
          position : `~astropy.coordinates.SkyCoord`
              Center position of the cutout region.
          width : tuple of `~astropy.coordinates.Angle`
              Angular sizes of the region in (lon, lat) in that specific order.
              If only one value is passed, a square region is extracted.
          mode : {'trim', 'partial', 'strict'}
              Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          cutout : `MapDataset`
              Cutout map dataset.
          """
          name = make_name(name)
          kwargs = {"gti": self.gti, "name": name}
          cutout_kwargs = {"position": position, "width": width, "mode": mode}
  
          if self.counts is not None:
              kwargs["counts"] = self.counts.cutout(**cutout_kwargs)
  
          if self.exposure is not None:
              kwargs["exposure"] = self.exposure.cutout(**cutout_kwargs)
  
          if self.background_model is not None:
              kwargs["models"] = self.background_model.cutout(
                  **cutout_kwargs, name=name
              )
  
          if self.edisp is not None:
              kwargs["edisp"] = self.edisp.cutout(**cutout_kwargs)
  
          if self.psf is not None:
              kwargs["psf"] = self.psf.cutout(**cutout_kwargs)
  
          if self.mask_safe is not None:
              kwargs["mask_safe"] = self.mask_safe.cutout(**cutout_kwargs)
  
          if self.mask_fit is not None:
              kwargs["mask_fit"] = self.mask_fit.cutout(**cutout_kwargs)
  
          return self.__class__(**kwargs)
  
  
  class MapDatasetOnOff(MapDataset):
      """Map dataset for on-off likelihood fitting.
  
      Parameters
      ----------
      models : `~gammapy.modeling.models.Models`
          Source sky models.
      counts : `~gammapy.maps.WcsNDMap`
          Counts cube
      counts_off : `~gammapy.maps.WcsNDMap`
          Ring-convolved counts cube
      acceptance : `~gammapy.maps.WcsNDMap`
          Acceptance from the IRFs
      acceptance_off : `~gammapy.maps.WcsNDMap`
          Acceptance off
      exposure : `~gammapy.maps.WcsNDMap`
          Exposure cube
      mask_fit : `~numpy.ndarray`
          Mask to apply to the likelihood for fitting.
      psf : `~gammapy.cube.PSFKernel`
          PSF kernel
      edisp : `~gammapy.irf.EDispKernel`
          Energy dispersion
      background_model : `~gammapy.modeling.models.BackgroundModel`
          Background model to use for the fit.
      evaluation_mode : {"local", "global"}
          Model evaluation mode.
          The "local" mode evaluates the model components on smaller grids to save computation time.
          This mode is recommended for local optimization algorithms.
          The "global" evaluation mode evaluates the model components on the full map.
          This mode is recommended for global optimization algorithms.
      mask_safe : `~numpy.ndarray`
          Mask defining the safe data range.
      gti : `~gammapy.data.GTI`
          GTI of the observation or union of GTI if it is a stacked observation
      name : str
          Name of the dataset.
  
      """
  
      stat_type = "wstat"
      tag = "MapDatasetOnOff"
  
      def __init__(
          self,
          models=None,
          counts=None,
          counts_off=None,
          acceptance=None,
          acceptance_off=None,
          exposure=None,
          mask_fit=None,
          psf=None,
          edisp=None,
          background_model=None,
          name=None,
          evaluation_mode="local",
          mask_safe=None,
          gti=None,
      ):
          if mask_fit is not None and mask_fit.dtype != np.dtype("bool"):
              raise ValueError("mask data must have dtype bool")
  
          self.evaluation_mode = evaluation_mode
          self.counts = counts
          self.counts_off = counts_off
  
          if np.isscalar(acceptance):
              acceptance = np.ones(self.data_shape) * acceptance
  
          if np.isscalar(acceptance_off):
              acceptance_off = np.ones(self.data_shape) * acceptance_off
  
          self.acceptance = acceptance
          self.acceptance_off = acceptance_off
          self.exposure = exposure
          self.background_model = None
          self.mask_fit = mask_fit
          self.psf = psf
          self.edisp = edisp
          self.models = models
          self._name = make_name(name)
          self.mask_safe = mask_safe
          self.gti = gti
  
      def __str__(self):
          str_ = super().__str__()
  
          counts_off = np.nan
          if self.counts_off is not None:
              counts_off = np.sum(self.counts_off.data)
          str_ += "\t{:32}: {:.0f} \n".format("Total counts_off", counts_off)
  
          acceptance = np.nan
          if self.acceptance is not None:
              acceptance = np.sum(self.acceptance.data)
          str_ += "\t{:32}: {:.0f} \n".format("Acceptance", acceptance)
  
          acceptance_off = np.nan
          if self.acceptance_off is not None:
              acceptance_off = np.sum(self.acceptance_off.data)
          str_ += "\t{:32}: {:.0f} \n".format("Acceptance off", acceptance_off)
  
          return str_.expandtabs(tabsize=4)
  
      @property
      def parameters(self):
          """List of parameters (`~gammapy.modeling.Parameters`)"""
          parameters = []
  
          if self.models:
              parameters += self.models.parameters
  
          return Parameters(parameters)
  
      @property
      def alpha(self):
          """Exposure ratio between signal and background regions"""
          alpha = self.acceptance / self.acceptance_off
          alpha.data = np.nan_to_num(alpha.data)
          return alpha
  
      @property
      def background(self):
          """Predicted background in the on region.
  
          Notice that this definition is valid under the assumption of cash statistic.
          """
          return self.alpha * self.counts_off
  
      @property
      def excess(self):
          """Excess (counts - alpha * counts_off)"""
          return self.counts.data - self.background.data
  
      def stat_array(self):
          """Likelihood per bin given the current model parameters"""
          mu_sig = self.npred().data
          on_stat_ = wstat(
              n_on=self.counts.data,
              n_off=self.counts_off.data,
              alpha=list(self.alpha.data),
              mu_sig=mu_sig,
          )
          return np.nan_to_num(on_stat_)
  
      @classmethod
      def from_geoms(
          cls,
          geom,
          geom_exposure,
          geom_psf,
          geom_edisp,
          reference_time="2000-01-01",
          name=None,
          **kwargs,
      ):
          """
          Create a MapDatasetOnOff object with zero filled maps according to the specified geometries
  
          Parameters
          ----------
          geom : `gammapy.maps.WcsGeom`
              geometry for the counts, counts_off, acceptance and acceptance_off maps
          geom_exposure : `gammapy.maps.WcsGeom`
              geometry for the exposure map
          geom_psf : `gammapy.maps.WcsGeom`
              geometry for the psf map
          geom_edisp : `gammapy.maps.WcsGeom`
              geometry for the energy dispersion map
          reference_time : `~astropy.time.Time`
              the reference time to use in GTI definition
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          empty_maps : `MapDatasetOnOff`
              A MapDatasetOnOff containing zero filled maps
          """
          kwargs = kwargs.copy()
          kwargs["name"] = name
  
          for key in ["counts", "counts_off", "acceptance", "acceptance_off"]:
              kwargs[key] = Map.from_geom(geom, unit="")
  
          kwargs["exposure"] = Map.from_geom(geom_exposure, unit="m2 s")
          kwargs["edisp"] = EDispMap.from_geom(geom_edisp)
          kwargs["psf"] = PSFMap.from_geom(geom_psf)
          kwargs["gti"] = GTI.create([] * u.s, [] * u.s, reference_time=reference_time)
          kwargs["mask_safe"] = Map.from_geom(geom, dtype=bool)
  
          return cls(**kwargs)
  
      @classmethod
      def from_map_dataset(
          cls, dataset, acceptance, acceptance_off, counts_off=None, name=None
      ):
          """Create spectrum dataseton off from another dataset.
  
          Parameters
          ----------
          dataset : `MapDataset`
              Spectrum dataset defining counts, edisp, aeff, livetime etc.
          acceptance : `Map`
              Relative background efficiency in the on region.
          acceptance_off : `Map`
              Relative background efficiency in the off region.
          counts_off : `Map`
              Off counts map . If the dataset provides a background model,
              and no off counts are defined. The off counts are deferred from
              counts_off / alpha.
          name : str
              Name of the returned dataset.
  
          Returns
          -------
          dataset : `MapDatasetOnOff`
              Map dataset on off.
  
          """
          kwargs = {"name": name}
  
          if counts_off is None and dataset.background_model is not None:
              alpha = acceptance / acceptance_off
              kwargs["counts_off"] = dataset.background_model.evaluate() / alpha
  
          return cls(
              counts=dataset.counts,
              exposure=dataset.exposure,
              counts_off=counts_off,
              edisp=dataset.edisp,
              gti=dataset.gti,
              mask_safe=dataset.mask_safe,
              mask_fit=dataset.mask_fit,
              acceptance=acceptance,
              acceptance_off=acceptance_off,
              name=dataset.name,
              evaluation_mode=dataset.evaluation_mode,
          )
  
      @property
      def _is_stackable(self):
          """Check if the Dataset contains enough information to be stacked"""
          if (
              self.acceptance_off is None
              or self.acceptance is None
              or self.counts_off is None
          ):
              return False
          else:
              return True
  
      def stack(self, other):
          r"""Stack another dataset in place.
  
          The ``acceptance`` of the stacked dataset is normalized to 1,
          and the stacked ``acceptance_off`` is scaled so that:
  
          .. math::
              \alpha_\text{stacked} =
              \frac{1}{a_\text{off}} =
              \frac{\alpha_1\text{OFF}_1 + \alpha_2\text{OFF}_2}{\text{OFF}_1 + OFF_2}
  
          Parameters
          ----------
          other : `MapDatasetOnOff`
              Other dataset
          """
          if not isinstance(other, MapDatasetOnOff):
              raise TypeError("Incompatible types for MapDatasetOnOff stacking")
  
          if not self._is_stackable or not other._is_stackable:
              raise ValueError("Cannot stack incomplete MapDatsetOnOff.")
  
          # Factor containing: self.alpha * self.counts_off + other.alpha * other.counts_off
          tmp_factor = (self.alpha * self.counts_off).copy()
          tmp_factor.data[~self.mask_safe.data] = 0
          tmp_factor.stack(other.alpha * other.counts_off, weights=other.mask_safe)
  
          # Stack the off counts (in place)
          self.counts_off.data[~self.mask_safe.data] = 0
          self.counts_off.stack(other.counts_off, weights=other.mask_safe)
  
          self.acceptance_off = self.counts_off / tmp_factor
          self.acceptance.data = np.ones(self.data_shape)
  
          super().stack(other)
  
      def stat_sum(self):
          """Total likelihood given the current model parameters."""
          return Dataset.stat_sum(self)
  
      def fake(self, background_model, random_state="random-seed"):
          """Simulate fake counts (on and off) for the current model and reduced IRFs.
  
          This method overwrites the counts defined on the dataset object.
  
          Parameters
          ----------
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
                  Defines random number generator initialisation.
                  Passed to `~gammapy.utils.random.get_random_state`.
          """
          random_state = get_random_state(random_state)
          npred = self.npred()
          npred.data = random_state.poisson(npred.data)
  
          npred_bkg = background_model.copy()
          npred_bkg.data = random_state.poisson(npred_bkg.data)
  
          self.counts = npred + npred_bkg
  
          npred_off = background_model / self.alpha
          npred_off.data = random_state.poisson(npred_off.data)
          self.counts_off = npred_off
  
      def to_hdulist(self):
          """Convert map dataset to list of HDUs.
  
          Returns
          -------
          hdulist : `~astropy.io.fits.HDUList`
              Map dataset list of HDUs.
          """
          hdulist = super().to_hdulist()
          exclude_primary = slice(1, None)
  
          if self.counts_off is not None:
              hdulist += self.counts_off.to_hdulist(hdu="counts_off")[exclude_primary]
  
          if self.acceptance is not None:
              hdulist += self.acceptance.to_hdulist(hdu="acceptance")[exclude_primary]
  
          if self.acceptance_off is not None:
              hdulist += self.acceptance_off.to_hdulist(hdu="acceptance_off")[
                  exclude_primary
              ]
  
          return hdulist
  
      @classmethod
      def from_hdulist(cls, hdulist, name=None):
          """Create map dataset from list of HDUs.
  
          Parameters
          ----------
          hdulist : `~astropy.io.fits.HDUList`
              List of HDUs.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDataset`
              Map dataset.
          """
          kwargs = {}
          kwargs["name"] = name
          if "COUNTS" in hdulist:
              kwargs["counts"] = Map.from_hdulist(hdulist, hdu="counts")
  
          if "COUNTS_OFF" in hdulist:
              kwargs["counts_off"] = Map.from_hdulist(hdulist, hdu="counts_off")
  
          if "ACCEPTANCE" in hdulist:
              kwargs["acceptance"] = Map.from_hdulist(hdulist, hdu="acceptance")
  
          if "ACCEPTANCE_OFF" in hdulist:
              kwargs["acceptance_off"] = Map.from_hdulist(hdulist, hdu="acceptance_off")
  
          if "EXPOSURE" in hdulist:
              kwargs["exposure"] = Map.from_hdulist(hdulist, hdu="exposure")
  
          if "EDISP_MATRIX" in hdulist:
              kwargs["edisp"] = EDispKernel.from_hdulist(
                  hdulist, hdu1="EDISP_MATRIX", hdu2="EDISP_MATRIX_EBOUNDS"
              )
  
          if "PSF_KERNEL" in hdulist:
              psf_map = Map.from_hdulist(hdulist, hdu="psf_kernel")
              kwargs["psf"] = PSFKernel(psf_map)
  
          if "MASK_SAFE" in hdulist:
              mask_safe_map = Map.from_hdulist(hdulist, hdu="mask_safe")
              kwargs["mask_safe"] = mask_safe_map.data.astype(bool)
  
          if "MASK_FIT" in hdulist:
              mask_fit_map = Map.from_hdulist(hdulist, hdu="mask_fit")
              kwargs["mask_fit"] = mask_fit_map.data.astype(bool)
  
          if "GTI" in hdulist:
              gti = GTI(Table.read(hdulist, hdu="GTI"))
              kwargs["gti"] = gti
          return cls(**kwargs)
  
      def to_spectrum_dataset(self, on_region, containment_correction=False, name=None):
          """Return a ~gammapy.spectrum.SpectrumDatasetOnOff from on_region.
  
          Counts and OFF counts are summed in the on_region.
  
          Acceptance is the average of all acceptances while acceptance OFF
          is taken such that number of excess is preserved in the on_region.
  
          Effective area is taken from the average exposure divided by the livetime.
          Here we assume it is the sum of the GTIs.
  
          The energy dispersion kernel is obtained at the on_region center.
          Only regions with centers are supported.
  
          The model is not exported to the ~gammapy.spectrum.SpectrumDataset.
          It must be set after the dataset extraction.
  
          Parameters
          ----------
          on_region : `~regions.SkyRegion`
              the input ON region on which to extract the spectrum
          containment_correction : bool
              Apply containment correction for point sources and circular on regions
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `~gammapy.spectrum.SpectrumDatasetOnOff`
              the resulting reduced dataset
          """
          dataset = super().to_spectrum_dataset(on_region, containment_correction, name)
  
          kwargs = {}
          if self.counts_off is not None:
              kwargs["counts_off"] = self.counts_off.get_spectrum(on_region, np.sum)
  
          if self.acceptance is not None:
              kwargs["acceptance"] = self.acceptance.get_spectrum(on_region, np.mean)
              background = self.background.get_spectrum(on_region, np.sum)
              kwargs["acceptance_off"] = (
                  kwargs["acceptance"] * kwargs["counts_off"] / background
              )
  
          return SpectrumDatasetOnOff.from_spectrum_dataset(dataset=dataset, **kwargs)
  
      def cutout(self, position, width, mode="trim", name=None):
          """Cutout map dataset.
  
          Parameters
          ----------
          position : `~astropy.coordinates.SkyCoord`
              Center position of the cutout region.
          width : tuple of `~astropy.coordinates.Angle`
              Angular sizes of the region in (lon, lat) in that specific order.
              If only one value is passed, a square region is extracted.
          mode : {'trim', 'partial', 'strict'}
              Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          cutout : `MapDatasetOnOff`
              Cutout map dataset.
          """
          cutout_kwargs = {
              "position": position,
              "width": width,
              "mode": mode,
              "name": name,
          }
  
          cutout_dataset = super().cutout(**cutout_kwargs)
  
          del cutout_kwargs["name"]
  
          if self.counts_off is not None:
              cutout_dataset.counts_off = self.counts_off.cutout(**cutout_kwargs)
  
          if self.acceptance is not None:
              cutout_dataset.acceptance = self.acceptance.cutout(**cutout_kwargs)
  
          if self.acceptance_off is not None:
              cutout_dataset.acceptance_off = self.acceptance_off.cutout(**cutout_kwargs)
  
          return cutout_dataset
  
      def to_image(self, spectrum=None, name=None):
          """Create images by summing over the energy axis.
  
          Exposure is weighted with an assumed spectrum,
          resulting in a weighted mean exposure image.
  
          Currently the PSFMap and EdispMap are dropped from the
          resulting image dataset.
  
          Parameters
          ----------
          spectrum : `~gammapy.modeling.models.SpectralModel`
              Spectral model to compute the weights.
              Default is power-law with spectral index of 2.
          name : str
              Name of the new dataset.
  
          Returns
          -------
          dataset : `MapDatasetOnOff`
              Map dataset containing images.
          """
          kwargs = {"name": name}
          dataset = super().to_image(spectrum, name)
  
          if self.mask_safe is not None:
              mask_safe = self.mask_safe
          else:
              mask_safe = 1
  
          if self.counts_off is not None:
              counts_off = self.counts_off * mask_safe
              kwargs["counts_off"] = counts_off.sum_over_axes(keepdims=True)
  
          if self.acceptance is not None:
              acceptance = self.acceptance * mask_safe
              kwargs["acceptance"] = acceptance.sum_over_axes(keepdims=True)
  
              background = self.background * mask_safe
              background = background.sum_over_axes(keepdims=True)
              kwargs["acceptance_off"] = (
                  kwargs["acceptance"] * kwargs["counts_off"] / background
              )
  
          return self.from_map_dataset(dataset, **kwargs)
  
  
  class MapEvaluator:
      """Sky model evaluation on maps.
  
      This evaluates a sky model on a 3D map and convolves with the IRFs,
      and returns a map of the predicted counts.
      Note that background counts are not added.
  
      For now, we only make it work for 3D WCS maps with an energy axis.
      No HPX, no other axes, those can be added later here or via new
      separate model evaluator classes.
  
      Parameters
      ----------
      model : `~gammapy.modeling.models.SkyModel`
          Sky model
      exposure : `~gammapy.maps.Map`
          Exposure map
      psf : `~gammapy.cube.PSFKernel`
          PSF kernel
      edisp : `~gammapy.irf.EDispKernel`
          Energy dispersion
      evaluation_mode : {"local", "global"}
          Model evaluation mode.
      """
  
      def __init__(
          self, model=None, exposure=None, psf=None, edisp=None, evaluation_mode="local"
      ):
          self.model = model
          self.exposure = exposure
          self.psf = psf
          self.edisp = edisp
          self.contributes = True
  
          if evaluation_mode not in {"local", "global"}:
              raise ValueError(f"Invalid evaluation_mode: {evaluation_mode!r}")
  
          self.evaluation_mode = evaluation_mode
  
          # TODO: this is preliminary solution until we have further unified the model handling
          if isinstance(model, BackgroundModel):
              self.compute_npred = model.evaluate
              self.evaluation_mode = "global"
  
      @property
      def geom(self):
          """True energy map geometry (`~gammapy.maps.Geom`)"""
          return self.exposure.geom
  
      @property
      def needs_update(self):
          """Check whether the model component has drifted away from its support."""
          if self.exposure is None:
              return True
          elif self.evaluation_mode == "global" or self.model.evaluation_radius is None:
              return False
          else:
              position = self.model.position
              separation = self._init_position.separation(position)
              update = separation > (self.model.evaluation_radius + CUTOUT_MARGIN)
          return update
  
      def update(self, exposure, psf, edisp, geom):
          """Update MapEvaluator, based on the current position of the model component.
  
          Parameters
          ----------
          exposure : `~gammapy.maps.Map`
              Exposure map.
          psf : `gammapy.cube.PSFMap`
              PSF map.
          edisp : `gammapy.cube.EDispMap`
              Edisp map.
          geom : `WcsGeom`
              Counts geom
          """
          log.debug("Updating model evaluator")
          # cache current position of the model component
  
          if isinstance(edisp, EDispMap):
              e_reco = geom.get_axis_by_name("energy").edges
              self.edisp = edisp.get_edisp_kernel(self.model.position, e_reco=e_reco)
          else:
              self.edisp = edisp
  
          if isinstance(psf, PSFMap):
              self.psf = psf.get_psf_kernel(self.model.position, geom=exposure.geom)
          else:
              self.psf = psf
  
          if self.evaluation_mode == "local" and self.model.evaluation_radius is not None:
              self._init_position = self.model.position
              if self.psf is not None:
                  psf_width = np.max(self.psf.psf_kernel_map.geom.width)
              else:
                  psf_width = 0 * u.deg
  
              width = psf_width + 2 * (self.model.evaluation_radius + CUTOUT_MARGIN)
              try:
                  self.exposure = exposure.cutout(
                      position=self.model.position, width=width
                  )
                  self.contributes = True
              except (NoOverlapError, ValueError):
                  self.contributes = False
          else:
              self.exposure = exposure
  
      def compute_dnde(self):
          """Compute model differential flux at map pixel centers.
  
          Returns
          -------
          model_map : `~gammapy.maps.Map`
              Sky cube with data filled with evaluated model values.
              Units: ``cm-2 s-1 TeV-1 deg-2``
          """
          return self.model.evaluate_geom(self.geom)
  
      def compute_flux(self):
          """Compute model integral flux over map pixel volumes.
  
          For now, we simply multiply dnde with bin volume.
          """
          dnde = self.compute_dnde()
          volume = self.geom.bin_volume()
          return dnde * volume
  
      def apply_exposure(self, flux):
          """Compute npred cube
  
          For now just divide flux cube by exposure
          """
          npred = (flux * self.exposure.quantity).to_value("")
          return Map.from_geom(self.geom, data=npred, unit="")
  
      def apply_psf(self, npred):
          """Convolve npred cube with PSF"""
          tmp = npred.convolve(self.psf)
          tmp.data[tmp.data < 0.0] = 0
          return tmp
  
      def apply_edisp(self, npred):
          """Convolve map data with energy dispersion.
  
          Parameters
          ----------
          npred : `~gammapy.maps.Map`
              Predicted counts in true energy bins
  
          Returns
          -------
          npred_reco : `~gammapy.maps.Map`
              Predicted counts in reco energy bins
          """
          return npred.apply_edisp(self.edisp)
  
      def compute_npred(self):
          """
          Evaluate model predicted counts.
  
          Returns
          -------
          npred : `~gammapy.maps.Map`
              Predicted counts on the map (in reco energy bins)
          """
          flux = self.compute_flux()
          npred = self.apply_exposure(flux)
          if self.psf is not None:
              npred = self.apply_psf(npred)
          if self.edisp is not None:
              npred = self.apply_edisp(npred)
  
          return npred
  
files: gammapy/cube/fit.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
  from gammapy.modeling import Parameter, Parameters
  from gammapy.modeling.parameter import _get_parameters_str
  from gammapy.utils.scripts import make_name, make_path
  from .core import Model, Models
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      processing = {"psf": 1, "edisp": 1}
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([self, *other])
          elif isinstance(other, (SkyModel, SkyDiffuseCube)):
              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):
          return self.evaluate(lon, lat, energy)
  
      def evaluate_geom(self, geom):
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat, coords["energy"])
  
      def _update_processing(self, processing):
          if processing in [None, "None"]:
              self.processing = {"psf": 0, "edisp": 0}
          else:
              for key in processing:
                  self.processing[key] = processing[key]
  
  
  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
      """
  
      tag = "SkyModel"
  
      def __init__(
          self,
          spectral_model,
          spatial_model=None,
          temporal_model=None,
          name=None,
          processing={},
      ):
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          self.temporal_model = temporal_model
          super().__init__()
          # TODO: this hack is needed for compound models to work
          self.__dict__.pop("_parameters")
  
          self._name = make_name(name)
          self._update_processing(processing)
  
      @property
      def name(self):
          return self._name
  
      @property
      def parameters(self):
          parameters = []
  
          if self.spatial_model is not None:
              parameters.append(self.spatial_model.parameters)
  
          parameters.append(self.spectral_model.parameters)
          return Parameters.from_stack(parameters)
  
      @property
      def spatial_model(self):
          """`~gammapy.modeling.models.SpatialModel`"""
          return self._spatial_model
  
      @spatial_model.setter
      def spatial_model(self, model):
          from .spatial import SpatialModel
  
          if not (model is None or isinstance(model, SpatialModel)):
              raise TypeError(f"Invalid type: {model!r}")
  
          self._spatial_model = model
  
      @property
      def spectral_model(self):
          """`~gammapy.modeling.models.SpectralModel`"""
          return self._spectral_model
  
      @spectral_model.setter
      def spectral_model(self, model):
          from .spectral import SpectralModel
  
          if not (model is None or isinstance(model, SpectralModel)):
              raise TypeError(f"Invalid type: {model!r}")
          self._spectral_model = model
  
      @property
      def temporal_model(self):
          """`~gammapy.modeling.models.TemporalModel`"""
          return self._temporal_model
  
      @temporal_model.setter
      def temporal_model(self, model):
          from .temporal import TemporalModel
  
          if not (model is None or isinstance(model, TemporalModel)):
              raise TypeError(f"Invalid type: {model!r}")
          self._temporal_model = model
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.spatial_model.position
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return self.spatial_model.evaluation_radius
  
      @property
      def frame(self):
          return self.spatial_model.frame
  
      def __repr__(self):
          return (
              f"{self.__class__.__name__}("
              f"spatial_model={self.spatial_model!r}, "
              f"spectral_model={self.spectral_model!r})"
              f"temporal_model={self.temporal_model!r})"
          )
  
      def evaluate(self, lon, lat, energy):
          """Evaluate the model at given points.
  
          The model evaluation follows numpy broadcasting rules.
  
          Return differential surface brightness cube.
          At the moment in units: ``cm-2 s-1 TeV-1 deg-2``
  
          Parameters
          ----------
          lon, lat : `~astropy.units.Quantity`
              Spatial coordinates
          energy : `~astropy.units.Quantity`
              Energy coordinate
  
          Returns
          -------
          value : `~astropy.units.Quantity`
              Model value at the given point.
          """
          value = self.spectral_model(energy)  # pylint:disable=not-callable
          # TODO: case if self.temporal_model is not None, introduce time in arguments ?
  
          if self.spatial_model is not None:
              value = value * self.spatial_model(lon, lat)  # pylint:disable=not-callable
  
          return value
  
      def evaluate_geom(self, geom):
          """Evaluate model on `~gammapy.maps.Geom`."""
          energy = geom.get_axis_by_name("energy").center[:, np.newaxis, np.newaxis]
          value = self.spectral_model(energy)
          # TODO: case with temporal_model is not None
          if self.spatial_model is not None:
              value = value * self.spatial_model.evaluate_geom(geom.to_image())
  
          return value
  
      def copy(self, name=None, **kwargs):
          """Copy SkyModel"""
          if self.spatial_model is not None:
              spatial_model = self.spatial_model.copy()
          else:
              spatial_model = None
  
          if self.temporal_model is not None:
              temporal_model = self.temporal_model.copy()
          else:
              temporal_model = None
  
          kwargs.setdefault("name", make_name(name))
          kwargs.setdefault("spectral_model", self.spectral_model.copy())
          kwargs.setdefault("spatial_model", spatial_model)
          kwargs.setdefault("temporal_model", temporal_model)
  
          return self.__class__(**kwargs)
  
      def to_dict(self):
          """Create dict for YAML serilisation"""
          data = {}
          data["name"] = self.name
          data["type"] = self.tag
          data["spectral"] = self.spectral_model.to_dict()
  
          if self.spatial_model is not None:
              data["spatial"] = self.spatial_model.to_dict()
  
          if self.temporal_model is not None:
              data["temporal"] = self.temporal_model.to_dict()
  
          if self.processing != {"psf": 1, "edisp": 1}:
              data["processing"] = self.processing
  
          return data
  
      @classmethod
      def from_dict(cls, data):
          """Create SkyModel from dict"""
          from gammapy.modeling.models import (
              SPATIAL_MODELS,
              SPECTRAL_MODELS,
              TEMPORAL_MODELS,
          )
  
          model_class = SPECTRAL_MODELS.get_cls(data["spectral"]["type"])
          spectral_model = model_class.from_dict(data["spectral"])
  
          spatial_data = data.get("spatial")
          if spatial_data is not None:
              model_class = SPATIAL_MODELS.get_cls(spatial_data["type"])
              spatial_model = model_class.from_dict(spatial_data)
          else:
              spatial_model = None
  
          temporal_data = data.get("temporal")
          if temporal_data is not None:
              model_class = TEMPORAL_MODELS.get_cls(temporal_data["type"])
              temporal_model = model_class.from_dict(temporal_data)
          else:
              temporal_model = None
  
          return cls(
              name=data["name"],
              spatial_model=spatial_model,
              spectral_model=spectral_model,
              temporal_model=temporal_model,
              processing=data.get("processing", {}),
          )
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
  
          str_ += "\t{:26}: {}\n".format("Spectral model type", self.spectral_model.tag)
  
          if self.spatial_model is not None:
              spatial_type = self.spatial_model.tag
          else:
              spatial_type = "None"
          str_ += "\t{:26}: {}\n".format("Spatial  model type", spatial_type)
  
          if self.temporal_model is not None:
              temporal_type = self.temporal_model.tag
          else:
              temporal_type = "None"
          str_ += "\t{:26}: {}\n".format("Temporal model type", temporal_type)
  
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self.parameters)
          lines = info.split("\n")
          str_ += "\t" + "\n\t".join(lines[:-1])
  
          str_ += "\n\n"
          return str_.expandtabs(tabsize=2)
  
  
  class SkyDiffuseCube(SkyModelBase):
      """Cube sky map template model (3D).
  
      This is for a 3D map with an energy axis.
      Use `~gammapy.modeling.models.TemplateSpatialModel` for 2D maps.
  
      Parameters
      ----------
      map : `~gammapy.maps.Map`
          Map template
      norm : float
          Norm parameter (multiplied with map values)
      tilt : float
          Additional tilt in the spectrum
      reference : `~astropy.units.Quantity`
          Reference energy of the tilt.
      meta : dict, optional
          Meta information, meta['filename'] will be used for serialization
      interp_kwargs : dict
          Interpolation keyword arguments passed to `gammapy.maps.Map.interp_by_coord`.
          Default arguments are {'interp': 'linear', 'fill_value': 0}.
      """
  
      tag = "SkyDiffuseCube"
      norm = Parameter("norm", 1)
      tilt = Parameter("tilt", 0, unit="", frozen=True)
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      def __init__(
          self,
          map,
          norm=norm.quantity,
          tilt=tilt.quantity,
          reference=reference.quantity,
          meta=None,
          interp_kwargs=None,
          name=None,
          filename=None,
          processing={},
      ):
  
          self._name = make_name(name)
  
          axis = map.geom.get_axis_by_name("energy")
  
          if axis.node_type != "center":
              raise ValueError('Need a map with energy axis node_type="center"')
  
          self.map = map
          self.meta = {} if meta is None else meta
          self.filename = filename
  
          interp_kwargs = {} if interp_kwargs is None else interp_kwargs
          interp_kwargs.setdefault("interp", "linear")
          interp_kwargs.setdefault("fill_value", 0)
          self._interp_kwargs = interp_kwargs
  
          # TODO: onve we have implement a more general and better model caching
          #  remove this again
          self._cached_value = None
          self._cached_coordinates = (None, None, None)
          self._update_processing(processing)
  
          super().__init__(norm=norm, tilt=tilt, reference=reference)
  
      @property
      def name(self):
          return self._name
  
      @classmethod
      def read(cls, filename, name=None, **kwargs):
          """Read map from FITS file.
  
          The default unit used if none is found in the file is ``cm-2 s-1 MeV-1 sr-1``.
  
          Parameters
          ----------
          filename : str
              FITS image filename.
          name : str
              Name of the output model
              The default used if none is filename.
          """
          m = Map.read(filename, **kwargs)
          if m.unit == "":
              m.unit = "cm-2 s-1 MeV-1 sr-1"
          if name is None:
              name = Path(filename).stem
          return cls(m, name=name, filename=filename)
  
      def _interpolate(self, lon, lat, energy):
          coord = {
              "lon": lon.to_value("deg"),
              "lat": lat.to_value("deg"),
              "energy": energy,
          }
          return self.map.interp_by_coord(coord, **self._interp_kwargs)
  
      def evaluate(self, lon, lat, energy):
          """Evaluate model."""
          is_cached_coord = [
              _ is coord for _, coord in zip((lon, lat, energy), self._cached_coordinates)
          ]
  
          # reset cache
          if not np.all(is_cached_coord):
              self._cached_value = None
  
          if self._cached_value is None:
              self._cached_coordinates = (lon, lat, energy)
              self._cached_value = self._interpolate(lon, lat, energy)
  
          norm = self.parameters["norm"].value
  
          tilt = self.parameters["tilt"].value
          reference = self.parameters["reference"].quantity
          tilt_factor = np.power((energy / reference).to(""), -tilt)
  
          val = norm * self._cached_value * tilt_factor.value
          return u.Quantity(val, self.map.unit, copy=False)
  
      def copy(self, name=None):
          """A shallow copy"""
          new = copy.copy(self)
          new._name = make_name(name)
          return new
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.map.geom.center_skydir
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return np.max(self.map.geom.width) / 2.0
  
      @property
      def frame(self):
          return self.position.frame.name
  
      @classmethod
      def from_dict(cls, data):
          model = cls.read(data["filename"])
          model._update_from_dict(data)
          processing = data.get("processing", {})
          model._update_processing(processing)
          return model
  
      def to_dict(self):
          data = super().to_dict()
          data["name"] = self.name
          data["type"] = data.pop("type")
          data["filename"] = self.filename
  
          # Move parameters at the end
          data["parameters"] = data.pop("parameters")
          if self.processing != {"psf": 1, "edisp": 1}:
              data["processing"] = self.processing
  
          return data
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self.parameters)
          lines = info.split("\n")
          str_ += "\t" + "\n\t".join(lines[:-1])
  
          str_ += "\n\n"
          return str_.expandtabs(tabsize=2)
  
  
  class BackgroundModel(Model):
      """Background model.
  
      Create a new map by a tilt and normalization on the available map
  
      Parameters
      ----------
      map : `~gammapy.maps.Map`
          Background model map
      norm : float
          Background normalization
      tilt : float
          Additional tilt in the spectrum
      reference : `~astropy.units.Quantity`
          Reference energy of the tilt.
      """
  
      tag = "BackgroundModel"
      norm = Parameter("norm", 1, unit="", min=0)
      tilt = Parameter("tilt", 0, unit="", frozen=True)
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      def __init__(
          self,
          map,
          norm=norm.quantity,
          tilt=tilt.quantity,
          reference=reference.quantity,
          name=None,
          filename=None,
      ):
          axis = map.geom.get_axis_by_name("energy")
          if axis.node_type != "edges":
              raise ValueError('Need an integrated map, energy axis node_type="edges"')
  
          self.map = map
  
          self._name = make_name(name)
          self.filename = filename
  
          super().__init__(norm=norm, tilt=tilt, reference=reference)
  
      @property
      def name(self):
          return self._name
  
      @property
      def energy_center(self):
          """True energy axis bin centers (`~astropy.units.Quantity`)"""
          energy_axis = self.map.geom.get_axis_by_name("energy")
          energy = energy_axis.center
          return energy[:, np.newaxis, np.newaxis]
  
      def evaluate(self):
          """Evaluate background model.
  
          Returns
          -------
          background_map : `~gammapy.maps.Map`
              Background evaluated on the Map
          """
          norm = self.parameters["norm"].value
          tilt = self.parameters["tilt"].value
          reference = self.parameters["reference"].quantity
          tilt_factor = np.power((self.energy_center / reference).to(""), -tilt)
          back_values = norm * self.map.data * tilt_factor.value
          return self.map.copy(data=back_values)
  
      def to_dict(self):
          data = {}
          data["name"] = self.name
          data.update(super().to_dict())
          if self.filename is not None:
              data["filename"] = self.filename
          data["parameters"] = data.pop("parameters")
          return data
  
      @classmethod
      def from_dict(cls, data):
          if "filename" in data:
              map = Map.read(data["filename"])
          elif "map" in data:
              map = data["map"]
          else:
              raise ValueError("Requires either filename or `Map` object")
  
          model = cls(map=map, name=data["name"])
          model._update_from_dict(data)
          return model
  
      def copy(self, name=None):
          """A deep copy."""
          new = copy.deepcopy(self)
          new._name = make_name(name)
          return new
  
      def cutout(self, position, width, mode="trim", name=None):
          """Cutout background model.
  
          Parameters
          ----------
          position : `~astropy.coordinates.SkyCoord`
              Center position of the cutout region.
          width : tuple of `~astropy.coordinates.Angle`
              Angular sizes of the region in (lon, lat) in that specific order.
              If only one value is passed, a square region is extracted.
          mode : {'trim', 'partial', 'strict'}
              Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
          name : str
              Name of the returned background model.
  
          Returns
          -------
          cutout : `BackgroundModel`
              Cutout background model.
          """
          cutout_kwargs = {"position": position, "width": width, "mode": mode}
  
          bkg_map = self.map.cutout(**cutout_kwargs)
          model = self.__class__(bkg_map, name=name)
          model.parameters.update_from_dict(self.parameters.to_dict())
          return model
  
      def stack(self, other, weights=None):
          """Stack background model in place.
  
          Stacking the background model resets the current parameters values.
  
          Parameters
          ----------
          other : `BackgroundModel`
              Other background model.
          """
          bkg = self.evaluate()
          other_bkg = other.evaluate()
          bkg.stack(other_bkg, weights=weights)
          self.map = bkg
  
          # reset parameter values
          self.norm.value = 1
          self.tilt.value = 0
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self.parameters)
          lines = info.split("\n")
          str_ += "\t" + "\n\t".join(lines[:-1])
  
          str_ += "\n\n"
          return str_.expandtabs(tabsize=2)
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.map.geom.center_skydir
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return np.max(self.map.geom.width) / 2.0
  
  
  def create_fermi_isotropic_diffuse_model(filename, **kwargs):
      """Read Fermi isotropic diffuse model.
  
      See `LAT Background models <https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html>`_
  
      Parameters
      ----------
      filename : str
          filename
      kwargs : dict
          Keyword arguments forwarded to `TemplateSpectralModel`
  
      Returns
      -------
      diffuse_model : `SkyModel`
          Fermi isotropic diffuse sky model.
      """
      from .spectral import TemplateSpectralModel
      from .spatial import ConstantSpatialModel
  
      vals = np.loadtxt(make_path(filename))
      energy = u.Quantity(vals[:, 0], "MeV", copy=False)
      values = u.Quantity(vals[:, 1], "MeV-1 s-1 cm-2", copy=False)
  
      spatial_model = ConstantSpatialModel()
      spectral_model = TemplateSpectralModel(energy=energy, values=values, **kwargs)
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
          name="fermi-diffuse-iso",
      )
  
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
  from gammapy.modeling import Parameter, Parameters
  from gammapy.modeling.parameter import _get_parameters_str
  from gammapy.utils.scripts import make_name, make_path
  from .core import Model, Models
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([self, *other])
          elif isinstance(other, (SkyModel, SkyDiffuseCube)):
              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):
          return self.evaluate(lon, lat, energy)
  
      def evaluate_geom(self, geom):
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat, coords["energy"])
  
  
  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
      """
  
      tag = "SkyModel"
  
      def __init__(
          self, spectral_model, spatial_model=None, temporal_model=None, name=None
      ):
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          self.temporal_model = temporal_model
          super().__init__()
          # TODO: this hack is needed for compound models to work
          self.__dict__.pop("_parameters")
  
          self._name = make_name(name)
  
      @property
      def name(self):
          return self._name
  
      @property
      def parameters(self):
          parameters = []
  
          if self.spatial_model is not None:
              parameters.append(self.spatial_model.parameters)
  
          parameters.append(self.spectral_model.parameters)
          return Parameters.from_stack(parameters)
  
      @property
      def spatial_model(self):
          """`~gammapy.modeling.models.SpatialModel`"""
          return self._spatial_model
  
      @spatial_model.setter
      def spatial_model(self, model):
          from .spatial import SpatialModel
  
          if not (model is None or isinstance(model, SpatialModel)):
              raise TypeError(f"Invalid type: {model!r}")
  
          self._spatial_model = model
  
      @property
      def spectral_model(self):
          """`~gammapy.modeling.models.SpectralModel`"""
          return self._spectral_model
  
      @spectral_model.setter
      def spectral_model(self, model):
          from .spectral import SpectralModel
  
          if not (model is None or isinstance(model, SpectralModel)):
              raise TypeError(f"Invalid type: {model!r}")
          self._spectral_model = model
  
      @property
      def temporal_model(self):
          """`~gammapy.modeling.models.TemporalModel`"""
          return self._temporal_model
  
      @temporal_model.setter
      def temporal_model(self, model):
          from .temporal import TemporalModel
  
          if not (model is None or isinstance(model, TemporalModel)):
              raise TypeError(f"Invalid type: {model!r}")
          self._temporal_model = model
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.spatial_model.position
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return self.spatial_model.evaluation_radius
  
      @property
      def frame(self):
          return self.spatial_model.frame
  
      def __repr__(self):
          return (
              f"{self.__class__.__name__}("
              f"spatial_model={self.spatial_model!r}, "
              f"spectral_model={self.spectral_model!r})"
              f"temporal_model={self.temporal_model!r})"
          )
  
      def evaluate(self, lon, lat, energy):
          """Evaluate the model at given points.
  
          The model evaluation follows numpy broadcasting rules.
  
          Return differential surface brightness cube.
          At the moment in units: ``cm-2 s-1 TeV-1 deg-2``
  
          Parameters
          ----------
          lon, lat : `~astropy.units.Quantity`
              Spatial coordinates
          energy : `~astropy.units.Quantity`
              Energy coordinate
  
          Returns
          -------
          value : `~astropy.units.Quantity`
              Model value at the given point.
          """
          value = self.spectral_model(energy)  # pylint:disable=not-callable
          # TODO: case if self.temporal_model is not None, introduce time in arguments ?
  
          if self.spatial_model is not None:
              value = value * self.spatial_model(lon, lat)  # pylint:disable=not-callable
  
          return value
  
      def evaluate_geom(self, geom):
          """Evaluate model on `~gammapy.maps.Geom`."""
          energy = geom.get_axis_by_name("energy").center[:, np.newaxis, np.newaxis]
          value = self.spectral_model(energy)
          # TODO: case with temporal_model is not None
          if self.spatial_model is not None:
              value = value * self.spatial_model.evaluate_geom(geom.to_image())
  
          return value
  
      def copy(self, name=None, **kwargs):
          """Copy SkyModel"""
          if self.spatial_model is not None:
              spatial_model = self.spatial_model.copy()
          else:
              spatial_model = None
  
          if self.temporal_model is not None:
              temporal_model = self.temporal_model.copy()
          else:
              temporal_model = None
  
          kwargs.setdefault("name", make_name(name))
          kwargs.setdefault("spectral_model", self.spectral_model.copy())
          kwargs.setdefault("spatial_model", spatial_model)
          kwargs.setdefault("temporal_model", temporal_model)
  
          return self.__class__(**kwargs)
  
      def to_dict(self):
          """Create dict for YAML serilisation"""
          data = {}
          data["name"] = self.name
          data["type"] = self.tag
          data["spectral"] = self.spectral_model.to_dict()
  
          if self.spatial_model is not None:
              data["spatial"] = self.spatial_model.to_dict()
  
          if self.temporal_model is not None:
              data["temporal"] = self.temporal_model.to_dict()
  
          return data
  
      @classmethod
      def from_dict(cls, data):
          """Create SkyModel from dict"""
          from gammapy.modeling.models import (
              SPATIAL_MODELS,
              SPECTRAL_MODELS,
              TEMPORAL_MODELS,
          )
  
          model_class = SPECTRAL_MODELS.get_cls(data["spectral"]["type"])
          spectral_model = model_class.from_dict(data["spectral"])
  
          spatial_data = data.get("spatial")
          if spatial_data is not None:
              model_class = SPATIAL_MODELS.get_cls(spatial_data["type"])
              spatial_model = model_class.from_dict(spatial_data)
          else:
              spatial_model = None
  
          temporal_data = data.get("temporal")
          if temporal_data is not None:
              model_class = TEMPORAL_MODELS.get_cls(temporal_data["type"])
              temporal_model = model_class.from_dict(temporal_data)
          else:
              temporal_model = None
  
          return cls(
              name=data["name"],
              spatial_model=spatial_model,
              spectral_model=spectral_model,
              temporal_model=temporal_model,
          )
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
  
          str_ += "\t{:26}: {}\n".format("Spectral model type", self.spectral_model.tag)
  
          if self.spatial_model is not None:
              spatial_type = self.spatial_model.tag
          else:
              spatial_type = "None"
          str_ += "\t{:26}: {}\n".format("Spatial  model type", spatial_type)
  
          if self.temporal_model is not None:
              temporal_type = self.temporal_model.tag
          else:
              temporal_type = "None"
          str_ += "\t{:26}: {}\n".format("Temporal model type", temporal_type)
  
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self.parameters)
          lines = info.split("\n")
          str_ += "\t" + "\n\t".join(lines[:-1])
  
          str_ += "\n\n"
          return str_.expandtabs(tabsize=2)
  
  
  class SkyDiffuseCube(SkyModelBase):
      """Cube sky map template model (3D).
  
      This is for a 3D map with an energy axis.
      Use `~gammapy.modeling.models.TemplateSpatialModel` for 2D maps.
  
      Parameters
      ----------
      map : `~gammapy.maps.Map`
          Map template
      norm : float
          Norm parameter (multiplied with map values)
      tilt : float
          Additional tilt in the spectrum
      reference : `~astropy.units.Quantity`
          Reference energy of the tilt.
      meta : dict, optional
          Meta information, meta['filename'] will be used for serialization
      interp_kwargs : dict
          Interpolation keyword arguments passed to `gammapy.maps.Map.interp_by_coord`.
          Default arguments are {'interp': 'linear', 'fill_value': 0}.
      """
  
      tag = "SkyDiffuseCube"
      norm = Parameter("norm", 1)
      tilt = Parameter("tilt", 0, unit="", frozen=True)
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      def __init__(
          self,
          map,
          norm=norm.quantity,
          tilt=tilt.quantity,
          reference=reference.quantity,
          meta=None,
          interp_kwargs=None,
          name=None,
          filename=None,
      ):
  
          self._name = make_name(name)
  
          axis = map.geom.get_axis_by_name("energy")
  
          if axis.node_type != "center":
              raise ValueError('Need a map with energy axis node_type="center"')
  
          self.map = map
          self.meta = {} if meta is None else meta
          self.filename = filename
  
          interp_kwargs = {} if interp_kwargs is None else interp_kwargs
          interp_kwargs.setdefault("interp", "linear")
          interp_kwargs.setdefault("fill_value", 0)
          self._interp_kwargs = interp_kwargs
  
          # TODO: onve we have implement a more general and better model caching
          #  remove this again
          self._cached_value = None
          self._cached_coordinates = (None, None, None)
  
          super().__init__(norm=norm, tilt=tilt, reference=reference)
  
      @property
      def name(self):
          return self._name
  
      @classmethod
      def read(cls, filename, name=None, **kwargs):
          """Read map from FITS file.
  
          The default unit used if none is found in the file is ``cm-2 s-1 MeV-1 sr-1``.
  
          Parameters
          ----------
          filename : str
              FITS image filename.
          name : str
              Name of the output model
              The default used if none is filename.
          """
          m = Map.read(filename, **kwargs)
          if m.unit == "":
              m.unit = "cm-2 s-1 MeV-1 sr-1"
          if name is None:
              name = Path(filename).stem
          return cls(m, name=name, filename=filename)
  
      def _interpolate(self, lon, lat, energy):
          coord = {
              "lon": lon.to_value("deg"),
              "lat": lat.to_value("deg"),
              "energy": energy,
          }
          return self.map.interp_by_coord(coord, **self._interp_kwargs)
  
      def evaluate(self, lon, lat, energy):
          """Evaluate model."""
          is_cached_coord = [
              _ is coord for _, coord in zip((lon, lat, energy), self._cached_coordinates)
          ]
  
          # reset cache
          if not np.all(is_cached_coord):
              self._cached_value = None
  
          if self._cached_value is None:
              self._cached_coordinates = (lon, lat, energy)
              self._cached_value = self._interpolate(lon, lat, energy)
  
          norm = self.parameters["norm"].value
  
          tilt = self.parameters["tilt"].value
          reference = self.parameters["reference"].quantity
          tilt_factor = np.power((energy / reference).to(""), -tilt)
  
          val = norm * self._cached_value * tilt_factor.value
          return u.Quantity(val, self.map.unit, copy=False)
  
      def copy(self, name=None):
          """A shallow copy"""
          new = copy.copy(self)
          new._name = make_name(name)
          return new
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.map.geom.center_skydir
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return np.max(self.map.geom.width) / 2.0
  
      @property
      def frame(self):
          return self.position.frame.name
  
      @classmethod
      def from_dict(cls, data):
          model = cls.read(data["filename"])
          model._update_from_dict(data)
          return model
  
      def to_dict(self):
          data = super().to_dict()
          data["name"] = self.name
          data["type"] = data.pop("type")
          data["filename"] = self.filename
  
          # Move parameters at the end
          data["parameters"] = data.pop("parameters")
          return data
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self.parameters)
          lines = info.split("\n")
          str_ += "\t" + "\n\t".join(lines[:-1])
  
          str_ += "\n\n"
          return str_.expandtabs(tabsize=2)
  
  
  class BackgroundModel(Model):
      """Background model.
  
      Create a new map by a tilt and normalization on the available map
  
      Parameters
      ----------
      map : `~gammapy.maps.Map`
          Background model map
      norm : float
          Background normalization
      tilt : float
          Additional tilt in the spectrum
      reference : `~astropy.units.Quantity`
          Reference energy of the tilt.
      """
  
      tag = "BackgroundModel"
      norm = Parameter("norm", 1, unit="", min=0)
      tilt = Parameter("tilt", 0, unit="", frozen=True)
      reference = Parameter("reference", "1 TeV", frozen=True)
  
      def __init__(
          self,
          map,
          norm=norm.quantity,
          tilt=tilt.quantity,
          reference=reference.quantity,
          name=None,
          filename=None,
      ):
          axis = map.geom.get_axis_by_name("energy")
          if axis.node_type != "edges":
              raise ValueError('Need an integrated map, energy axis node_type="edges"')
  
          self.map = map
  
          self._name = make_name(name)
          self.filename = filename
  
          super().__init__(norm=norm, tilt=tilt, reference=reference)
  
      @property
      def name(self):
          return self._name
  
      @property
      def energy_center(self):
          """True energy axis bin centers (`~astropy.units.Quantity`)"""
          energy_axis = self.map.geom.get_axis_by_name("energy")
          energy = energy_axis.center
          return energy[:, np.newaxis, np.newaxis]
  
      def evaluate(self):
          """Evaluate background model.
  
          Returns
          -------
          background_map : `~gammapy.maps.Map`
              Background evaluated on the Map
          """
          norm = self.parameters["norm"].value
          tilt = self.parameters["tilt"].value
          reference = self.parameters["reference"].quantity
          tilt_factor = np.power((self.energy_center / reference).to(""), -tilt)
          back_values = norm * self.map.data * tilt_factor.value
          return self.map.copy(data=back_values)
  
      def to_dict(self):
          data = {}
          data["name"] = self.name
          data.update(super().to_dict())
          if self.filename is not None:
              data["filename"] = self.filename
          data["parameters"] = data.pop("parameters")
          return data
  
      @classmethod
      def from_dict(cls, data):
          if "filename" in data:
              map = Map.read(data["filename"])
          elif "map" in data:
              map = data["map"]
          else:
              raise ValueError("Requires either filename or `Map` object")
  
          model = cls(map=map, name=data["name"])
          model._update_from_dict(data)
          return model
  
      def copy(self, name=None):
          """A deep copy."""
          new = copy.deepcopy(self)
          new._name = make_name(name)
          return new
  
      def cutout(self, position, width, mode="trim", name=None):
          """Cutout background model.
  
          Parameters
          ----------
          position : `~astropy.coordinates.SkyCoord`
              Center position of the cutout region.
          width : tuple of `~astropy.coordinates.Angle`
              Angular sizes of the region in (lon, lat) in that specific order.
              If only one value is passed, a square region is extracted.
          mode : {'trim', 'partial', 'strict'}
              Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
          name : str
              Name of the returned background model.
  
          Returns
          -------
          cutout : `BackgroundModel`
              Cutout background model.
          """
          cutout_kwargs = {"position": position, "width": width, "mode": mode}
  
          bkg_map = self.map.cutout(**cutout_kwargs)
          model = self.__class__(bkg_map, name=name)
          model.parameters.update_from_dict(self.parameters.to_dict())
          return model
  
      def stack(self, other, weights=None):
          """Stack background model in place.
  
          Stacking the background model resets the current parameters values.
  
          Parameters
          ----------
          other : `BackgroundModel`
              Other background model.
          """
          bkg = self.evaluate()
          other_bkg = other.evaluate()
          bkg.stack(other_bkg, weights=weights)
          self.map = bkg
  
          # reset parameter values
          self.norm.value = 1
          self.tilt.value = 0
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self.parameters)
          lines = info.split("\n")
          str_ += "\t" + "\n\t".join(lines[:-1])
  
          str_ += "\n\n"
          return str_.expandtabs(tabsize=2)
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.map.geom.center_skydir
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return np.max(self.map.geom.width) / 2.0
  
  
  def create_fermi_isotropic_diffuse_model(filename, **kwargs):
      """Read Fermi isotropic diffuse model.
  
      See `LAT Background models <https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html>`_
  
      Parameters
      ----------
      filename : str
          filename
      kwargs : dict
          Keyword arguments forwarded to `TemplateSpectralModel`
  
      Returns
      -------
      diffuse_model : `SkyModel`
          Fermi isotropic diffuse sky model.
      """
      from .spectral import TemplateSpectralModel
      from .spatial import ConstantSpatialModel
  
      vals = np.loadtxt(make_path(filename))
      energy = u.Quantity(vals[:, 0], "MeV", copy=False)
      values = u.Quantity(vals[:, 1], "MeV-1 s-1 cm-2", copy=False)
  
      spatial_model = ConstantSpatialModel()
      spectral_model = TemplateSpectralModel(energy=energy, values=values, **kwargs)
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
          name="fermi-diffuse-iso",
      )
  
files: gammapy/modeling/models/cube.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import copy
  import logging
  import numpy as np
  import astropy.units as u
  from astropy.io import fits
  from astropy.table import Table
  from astropy.wcs import WCS
  from regions import fits_region_objects_to_table, FITSRegionParser
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.utils.fits import ebounds_to_energy_axis, energy_axis_to_ebounds
  from gammapy.utils.regions import compound_region_to_list, list_to_compound_region
  from gammapy.utils.scripts import make_path
  
  __all__ = ["CountsSpectrum"]
  
  log = logging.getLogger(__name__)
  
  
  class CountsSpectrum:
      """Generic counts spectrum.
  
      Parameters
      ----------
      energy_lo : `~astropy.units.Quantity`
          Lower bin edges of energy axis
      energy_hi : `~astropy.units.Quantity`
          Upper bin edges of energy axis
      data : `~numpy.ndarray`
          Spectrum data.
      unit : str or `~astropy.units.Unit`
          Data unit
      region : `~regions.SkyRegion`
          Region the spectrum is defined for.
      wcs : `~astropy.wcs.WCS`
          the wcs system used to perform region based event selection
  
      Examples
      --------
      .. plot::
          :include-source:
  
          from gammapy.spectrum import CountsSpectrum
          import numpy as np
          import astropy.units as u
  
          ebounds = np.logspace(0,1,11) * u.TeV
          data = np.arange(10)
          spec = CountsSpectrum(
              energy_lo=ebounds[:-1],
              energy_hi=ebounds[1:],
              data=data,
          )
          spec.plot(show_poisson_errors=True)
      """
  
      def __init__(self, energy_lo, energy_hi, data=None, unit="", region=None, wcs=None):
          e_edges = edges_from_lo_hi(energy_lo, energy_hi)
          self.energy = MapAxis.from_edges(e_edges, interp="log", name="energy")
  
          if data is None:
              data = np.zeros(self.energy.nbin)
  
          self.data = np.array(data)
          if not self.energy.nbin == self.data.size:
              raise ValueError("Incompatible data and energy axis size.")
  
          self.unit = u.Unit(unit)
          self.region = region
          self.wcs = wcs
  
      @property
      def quantity(self):
          return self.data * self.unit
  
      @quantity.setter
      def quantity(self, quantity):
          self.data = quantity.value
          self.unit = quantity.unit
  
      @staticmethod
      def read_region_table(hdu):
          """Read region table and convert it to region list."""
          region_table = Table.read(hdu)
          parser = FITSRegionParser(region_table)
          pix_region = parser.shapes.to_regions()
          wcs = WCS(region_table.meta)
          regions = []
          for reg in pix_region:
              regions.append(reg.to_sky(wcs))
          region = list_to_compound_region(regions)
          return region, wcs
  
      @classmethod
      def from_hdulist(cls, hdulist, hdu1="COUNTS", hdu2="EBOUNDS", hdu3="REGION"):
          """Read from HDU list in OGIP format."""
          table = Table.read(hdulist[hdu1])
          counts = table["COUNTS"].data
          ebounds = ebounds_to_energy_axis(hdulist[hdu2])
  
          # TODO: add region serilisation
          region = None
          wcs = None
          if hdu3 in hdulist:
              region, wcs = cls.read_region_table(hdulist[hdu3])
  
          return cls(
              data=counts,
              energy_lo=ebounds[:-1],
              energy_hi=ebounds[1:],
              region=region,
              wcs=wcs,
          )
  
      @classmethod
      def read(cls, filename, hdu1="COUNTS", hdu2="EBOUNDS", hdu3="REGION"):
          """Read from file."""
          with fits.open(make_path(filename), memmap=False) as hdulist:
              return cls.from_hdulist(hdulist, hdu1=hdu1, hdu2=hdu2, hdu3=hdu3)
  
      def to_table(self):
          """Convert to `~astropy.table.Table`.
  
          Data format specification: :ref:`gadf:ogip-pha`
          """
          channel = np.arange(self.energy.nbin, dtype=np.int16)
          counts = np.array(self.data, dtype=np.int32)
  
          names = ["CHANNEL", "COUNTS"]
          meta = {"name": "COUNTS"}
  
          return Table([channel, counts], names=names, meta=meta)
  
      def _to_region_table(self):
          """Export region to a FITS region table."""
          region_list = compound_region_to_list(self.region)
          pixel_region_list = []
          for reg in region_list:
              pixel_region_list.append(reg.to_pixel(self.wcs))
          table = fits_region_objects_to_table(pixel_region_list)
          table.meta.update(self.wcs.to_header())
          return table
  
      def to_hdulist(self, use_sherpa=False):
          """Convert to `~astropy.io.fits.HDUList`.
  
          This adds an ``EBOUNDS`` extension to the ``BinTableHDU`` produced by
          ``to_table``, in order to store the energy axis
          """
          table = self.to_table()
          name = table.meta["name"]
          hdu = fits.BinTableHDU(table, name=name)
  
          energy = self.energy.edges
  
          if use_sherpa:
              energy = energy.to("keV")
  
          ebounds = energy_axis_to_ebounds(energy)
  
          region_table = self._to_region_table()
          region_hdu = fits.BinTableHDU(region_table, name="REGION")
          return fits.HDUList([fits.PrimaryHDU(), hdu, ebounds, region_hdu])
  
      def write(self, filename, use_sherpa=False, **kwargs):
          """Write to file."""
          filename = make_path(filename)
          self.to_hdulist(use_sherpa=use_sherpa).writeto(filename, **kwargs)
  
      def fill_events(self, events):
          """Fill events (`gammapy.data.EventList`)."""
          self.fill_energy(events.energy)
  
      def fill_energy(self, energy):
          """Fill energy values (`~astropy.units.Quantity`)"""
          energy = energy.to_value(self.energy.unit)
          edges = self.energy.edges.to_value(self.energy.unit)
          self.data = np.histogram(energy, edges)[0]
  
      @property
      def total_counts(self):
          """Total number of counts."""
          return self.data.sum()
  
      def plot(
          self,
          ax=None,
          energy_unit="TeV",
          show_poisson_errors=False,
          show_energy=None,
          **kwargs,
      ):
          """Plot as data points.
  
          kwargs are forwarded to `~matplotlib.pyplot.errorbar`
  
          Parameters
          ----------
          ax : `~matplotlib.axis` (optional)
              Axis instance to be used for the plot
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          show_poisson_errors : bool, optional
              Show poisson errors on the plot
          show_energy : `~astropy.units.Quantity`, optional
              Show energy, e.g. threshold, as vertical line
  
          Returns
          -------
          ax: `~matplotlib.axis`
              Axis instance used for the plot
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          counts = self.data
          x = self.energy.center.to_value(energy_unit)
          bounds = self.energy.edges.to_value(energy_unit)
          xerr = [x - bounds[:-1], bounds[1:] - x]
          yerr = np.sqrt(counts) if show_poisson_errors else 0
          kwargs.setdefault("fmt", ".")
          ax.errorbar(x, counts, xerr=xerr, yerr=yerr, **kwargs)
          if show_energy is not None:
              ener_val = u.Quantity(show_energy).to_value(energy_unit)
              ax.vlines(ener_val, 0, 1.1 * max(self.data), linestyles="dashed")
          ax.set_xlabel(f"Energy [{energy_unit}]")
          ax.set_ylabel("Counts")
          ax.set_xscale("log")
          return ax
  
      def plot_hist(self, ax=None, energy_unit="TeV", show_energy=None, **kwargs):
          """Plot as histogram.
  
          kwargs are forwarded to `~matplotlib.pyplot.hist`
  
          Parameters
          ----------
          ax : `~matplotlib.axis` (optional)
              Axis instance to be used for the plot
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          show_energy : `~astropy.units.Quantity`, optional
              Show energy, e.g. threshold, as vertical line
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          kwargs.setdefault("lw", 2)
          kwargs.setdefault("histtype", "step")
          weights = self.data
          bins = self.energy.edges.to_value(energy_unit)
          x = self.energy.center.to_value(energy_unit)
          ax.hist(x, bins=bins, weights=weights, **kwargs)
          if show_energy is not None:
              ener_val = u.Quantity(show_energy).to_value(energy_unit)
              ax.vlines(ener_val, 0, 1.1 * max(self.data), linestyles="dashed")
          ax.set_xlabel(f"Energy [{energy_unit}]")
          ax.set_ylabel("Counts")
          ax.set_xscale("log")
          return ax
  
      def plot_region(self, ax=None, **kwargs):
          """Plot region
  
          Parameters
          ----------
          ax : `~astropy.vizualisation.WCSAxes`
              Axes to plot on.
          **kwargs : dict
              Keyword arguments forwarded to `~regions.PixelRegion.as_artist`
          """
          import matplotlib.pyplot as plt
          from matplotlib.collections import PatchCollection
  
          ax = plt.gca() or ax
          regions = compound_region_to_list(self.region)
          artists = [region.to_pixel(wcs=ax.wcs).as_artist() for region in regions]
  
          patches = PatchCollection(artists, **kwargs)
          ax.add_collection(patches)
          return ax
  
      def peek(self, figsize=(5, 10)):
          """Quick-look summary plots."""
          import matplotlib.pyplot as plt
  
          fig, ax = plt.subplots(1, 1, figsize=figsize)
          self.plot_hist(ax=ax)
          return ax
  
      def copy(self):
          """A deep copy of self."""
          return copy.deepcopy(self)
  
      def downsample(self, factor):
          """Downsample spectrum.
  
          Parameters
          ----------
          factor : int
              Downsampling factor.
  
          Returns
          -------
          spectrum : `~gammapy.spectrum.CountsSpectrum`
              Downsampled spectrum.
          """
          from gammapy.extern.skimage import block_reduce
  
          data = block_reduce(self.data, block_size=(factor,))
          energy = self.energy.downsample(factor).edges
          return self.__class__(energy_lo=energy[:-1], energy_hi=energy[1:], data=data)
  
      def energy_mask(self, emin=None, emax=None):
          """Create a mask for a given energy range.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Energy range
          """
          edges = self.energy.edges
  
          # set default values
          emin = emin if emin is not None else edges[0]
          emax = emax if emax is not None else edges[-1]
  
          return (edges[:-1] >= emin) & (edges[1:] <= emax)
  
      def _arithmetics(self, operator, other, copy):
          """Perform arithmetics on maps after checking geometry consistency."""
          if isinstance(other, CountsSpectrum):
              q = other.quantity
          else:
              q = u.Quantity(other, copy=False)
  
          out = self.copy() if copy else self
          out.quantity = operator(out.quantity, q)
          return out
  
      def __add__(self, other):
          return self._arithmetics(np.add, other, copy=True)
  
      def __iadd__(self, other):
          return self._arithmetics(np.add, other, copy=False)
  
      def __sub__(self, other):
          return self._arithmetics(np.subtract, other, copy=True)
  
      def __isub__(self, other):
          return self._arithmetics(np.subtract, other, copy=False)
  
      def __mul__(self, other):
          return self._arithmetics(np.multiply, other, copy=True)
  
      def __imul__(self, other):
          return self._arithmetics(np.multiply, other, copy=False)
  
      def __truediv__(self, other):
          return self._arithmetics(np.true_divide, other, copy=True)
  
      def __itruediv__(self, other):
          return self._arithmetics(np.true_divide, other, copy=False)
  
      def __array__(self):
          return self.data
  
  
  class SpectrumEvaluator:
      """Calculate number of predicted counts (``npred``).
  
      The true and reconstructed energy binning are inferred from the provided IRFs.
  
      Parameters
      ----------
      model : `~gammapy.modeling.models.SkyModel`
          Spectral model
      aeff : `~gammapy.irf.EffectiveAreaTable`
          EffectiveArea
      edisp : `~gammapy.irf.EDispKernel`, optional
          Energy dispersion kernel.
      livetime : `~astropy.units.Quantity`
          Observation duration (may be contained in aeff)
      e_true : `~astropy.units.Quantity`, optional
          Desired energy axis of the prediced counts vector if no IRFs are given
  
      Examples
      --------
      Calculate predicted counts in a desired reconstruced energy binning
  
      .. plot::
          :include-source:
  
          import numpy as np
          import astropy.units as u
          import matplotlib.pyplot as plt
          from gammapy.irf import EffectiveAreaTable, EDispKernel
          from gammapy.modeling.models import PowerLawSpectralModel, SkyModel
          from gammapy.spectrum import SpectrumEvaluator
  
          e_true = np.logspace(-2, 2.5, 109) * u.TeV
          e_reco = np.logspace(-2, 2, 73) * u.TeV
  
          aeff = EffectiveAreaTable.from_parametrization(energy=e_true)
          edisp = EDispKernel.from_gauss(e_true=e_true, e_reco=e_reco, sigma=0.3, bias=0)
  
          pwl = PowerLawSpectralModel(index=2.3, amplitude="2.5e-12 cm-2 s-1 TeV-1", reference="1 TeV")
          model = SkyModel(spectral_model=pwl)
  
          predictor = SpectrumEvaluator(model=model, aeff=aeff, edisp=edisp, livetime="1 hour")
          predictor.compute_npred().plot_hist()
          plt.show()
      """
  
      def __init__(self, model, aeff=None, edisp=None, livetime=None):
          self.model = model
          self.aeff = aeff
          self.edisp = edisp
          self.livetime = livetime
  
      def compute_npred(self):
          e_true = self.aeff.energy.edges
          integral_flux = self.model.spectral_model.integral(
              emin=e_true[:-1], emax=e_true[1:], intervals=True
          )
  
          true_counts = self.apply_aeff(integral_flux)
          return self.apply_edisp(true_counts)
  
      def apply_aeff(self, integral_flux):
          if self.aeff is not None:
              cts = integral_flux * self.aeff.data.data
          else:
              cts = integral_flux
  
          # Multiply with livetime if not already contained in aeff or model
          if cts.unit.is_equivalent("s-1"):
              cts *= self.livetime
  
          return cts.to("")
  
      def apply_edisp(self, true_counts):
          from . import CountsSpectrum
  
          if self.edisp is not None and self.model.processing["edisp"] == True:
              cts = self.edisp.apply(true_counts)
              e_reco = self.edisp.e_reco.edges
          else:
              cts = true_counts
              e_reco = self.aeff.energy.edges
  
          return CountsSpectrum(data=cts, energy_lo=e_reco[:-1], energy_hi=e_reco[1:])
  
files: gammapy/spectrum/core.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import copy
  import logging
  import numpy as np
  import astropy.units as u
  from astropy.io import fits
  from astropy.table import Table
  from astropy.wcs import WCS
  from regions import fits_region_objects_to_table, FITSRegionParser
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.utils.fits import ebounds_to_energy_axis, energy_axis_to_ebounds
  from gammapy.utils.regions import compound_region_to_list, list_to_compound_region
  from gammapy.utils.scripts import make_path
  
  __all__ = ["CountsSpectrum"]
  
  log = logging.getLogger(__name__)
  
  
  class CountsSpectrum:
      """Generic counts spectrum.
  
      Parameters
      ----------
      energy_lo : `~astropy.units.Quantity`
          Lower bin edges of energy axis
      energy_hi : `~astropy.units.Quantity`
          Upper bin edges of energy axis
      data : `~numpy.ndarray`
          Spectrum data.
      unit : str or `~astropy.units.Unit`
          Data unit
      region : `~regions.SkyRegion`
          Region the spectrum is defined for.
      wcs : `~astropy.wcs.WCS`
          the wcs system used to perform region based event selection
  
      Examples
      --------
      .. plot::
          :include-source:
  
          from gammapy.spectrum import CountsSpectrum
          import numpy as np
          import astropy.units as u
  
          ebounds = np.logspace(0,1,11) * u.TeV
          data = np.arange(10)
          spec = CountsSpectrum(
              energy_lo=ebounds[:-1],
              energy_hi=ebounds[1:],
              data=data,
          )
          spec.plot(show_poisson_errors=True)
      """
  
      def __init__(self, energy_lo, energy_hi, data=None, unit="", region=None, wcs=None):
          e_edges = edges_from_lo_hi(energy_lo, energy_hi)
          self.energy = MapAxis.from_edges(e_edges, interp="log", name="energy")
  
          if data is None:
              data = np.zeros(self.energy.nbin)
  
          self.data = np.array(data)
          if not self.energy.nbin == self.data.size:
              raise ValueError("Incompatible data and energy axis size.")
  
          self.unit = u.Unit(unit)
          self.region = region
          self.wcs = wcs
  
      @property
      def quantity(self):
          return self.data * self.unit
  
      @quantity.setter
      def quantity(self, quantity):
          self.data = quantity.value
          self.unit = quantity.unit
  
      @staticmethod
      def read_region_table(hdu):
          """Read region table and convert it to region list."""
          region_table = Table.read(hdu)
          parser = FITSRegionParser(region_table)
          pix_region = parser.shapes.to_regions()
          wcs = WCS(region_table.meta)
          regions = []
          for reg in pix_region:
              regions.append(reg.to_sky(wcs))
          region = list_to_compound_region(regions)
          return region, wcs
  
      @classmethod
      def from_hdulist(cls, hdulist, hdu1="COUNTS", hdu2="EBOUNDS", hdu3="REGION"):
          """Read from HDU list in OGIP format."""
          table = Table.read(hdulist[hdu1])
          counts = table["COUNTS"].data
          ebounds = ebounds_to_energy_axis(hdulist[hdu2])
  
          # TODO: add region serilisation
          region = None
          wcs = None
          if hdu3 in hdulist:
              region, wcs =cls.read_region_table(hdulist[hdu3])
  
          return cls(data=counts, energy_lo=ebounds[:-1], energy_hi=ebounds[1:], region=region, wcs=wcs)
  
      @classmethod
      def read(cls, filename, hdu1="COUNTS", hdu2="EBOUNDS", hdu3="REGION"):
          """Read from file."""
          with fits.open(make_path(filename), memmap=False) as hdulist:
              return cls.from_hdulist(hdulist, hdu1=hdu1, hdu2=hdu2, hdu3=hdu3)
  
      def to_table(self):
          """Convert to `~astropy.table.Table`.
  
          Data format specification: :ref:`gadf:ogip-pha`
          """
          channel = np.arange(self.energy.nbin, dtype=np.int16)
          counts = np.array(self.data, dtype=np.int32)
  
          names = ["CHANNEL", "COUNTS"]
          meta = {"name": "COUNTS"}
  
  
  
          return Table([channel, counts], names=names, meta=meta)
  
      def _to_region_table(self):
          """Export region to a FITS region table."""
          region_list = compound_region_to_list(self.region)
          pixel_region_list = []
          for reg in region_list:
              pixel_region_list.append(reg.to_pixel(self.wcs))
          table = fits_region_objects_to_table(pixel_region_list)
          table.meta.update(self.wcs.to_header())
          return table
  
      def to_hdulist(self, use_sherpa=False):
          """Convert to `~astropy.io.fits.HDUList`.
  
          This adds an ``EBOUNDS`` extension to the ``BinTableHDU`` produced by
          ``to_table``, in order to store the energy axis
          """
          table = self.to_table()
          name = table.meta["name"]
          hdu = fits.BinTableHDU(table, name=name)
  
          energy = self.energy.edges
  
          if use_sherpa:
              energy = energy.to("keV")
  
          ebounds = energy_axis_to_ebounds(energy)
  
          region_table = self._to_region_table()
          region_hdu = fits.BinTableHDU(region_table, name='REGION')
          return fits.HDUList([fits.PrimaryHDU(), hdu, ebounds, region_hdu])
  
      def write(self, filename, use_sherpa=False, **kwargs):
          """Write to file."""
          filename = make_path(filename)
          self.to_hdulist(use_sherpa=use_sherpa).writeto(filename, **kwargs)
  
      def fill_events(self, events):
          """Fill events (`gammapy.data.EventList`)."""
          self.fill_energy(events.energy)
  
      def fill_energy(self, energy):
          """Fill energy values (`~astropy.units.Quantity`)"""
          energy = energy.to_value(self.energy.unit)
          edges = self.energy.edges.to_value(self.energy.unit)
          self.data = np.histogram(energy, edges)[0]
  
      @property
      def total_counts(self):
          """Total number of counts."""
          return self.data.sum()
  
      def plot(
          self,
          ax=None,
          energy_unit="TeV",
          show_poisson_errors=False,
          show_energy=None,
          **kwargs,
      ):
          """Plot as data points.
  
          kwargs are forwarded to `~matplotlib.pyplot.errorbar`
  
          Parameters
          ----------
          ax : `~matplotlib.axis` (optional)
              Axis instance to be used for the plot
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          show_poisson_errors : bool, optional
              Show poisson errors on the plot
          show_energy : `~astropy.units.Quantity`, optional
              Show energy, e.g. threshold, as vertical line
  
          Returns
          -------
          ax: `~matplotlib.axis`
              Axis instance used for the plot
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          counts = self.data
          x = self.energy.center.to_value(energy_unit)
          bounds = self.energy.edges.to_value(energy_unit)
          xerr = [x - bounds[:-1], bounds[1:] - x]
          yerr = np.sqrt(counts) if show_poisson_errors else 0
          kwargs.setdefault("fmt", ".")
          ax.errorbar(x, counts, xerr=xerr, yerr=yerr, **kwargs)
          if show_energy is not None:
              ener_val = u.Quantity(show_energy).to_value(energy_unit)
              ax.vlines(ener_val, 0, 1.1 * max(self.data), linestyles="dashed")
          ax.set_xlabel(f"Energy [{energy_unit}]")
          ax.set_ylabel("Counts")
          ax.set_xscale("log")
          return ax
  
      def plot_hist(self, ax=None, energy_unit="TeV", show_energy=None, **kwargs):
          """Plot as histogram.
  
          kwargs are forwarded to `~matplotlib.pyplot.hist`
  
          Parameters
          ----------
          ax : `~matplotlib.axis` (optional)
              Axis instance to be used for the plot
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          show_energy : `~astropy.units.Quantity`, optional
              Show energy, e.g. threshold, as vertical line
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          kwargs.setdefault("lw", 2)
          kwargs.setdefault("histtype", "step")
          weights = self.data
          bins = self.energy.edges.to_value(energy_unit)
          x = self.energy.center.to_value(energy_unit)
          ax.hist(x, bins=bins, weights=weights, **kwargs)
          if show_energy is not None:
              ener_val = u.Quantity(show_energy).to_value(energy_unit)
              ax.vlines(ener_val, 0, 1.1 * max(self.data), linestyles="dashed")
          ax.set_xlabel(f"Energy [{energy_unit}]")
          ax.set_ylabel("Counts")
          ax.set_xscale("log")
          return ax
  
      def plot_region(self, ax=None, **kwargs):
          """Plot region
  
          Parameters
          ----------
          ax : `~astropy.vizualisation.WCSAxes`
              Axes to plot on.
          **kwargs : dict
              Keyword arguments forwarded to `~regions.PixelRegion.as_artist`
          """
          import matplotlib.pyplot as plt
          from matplotlib.collections import PatchCollection
  
          ax = plt.gca() or ax
          regions = compound_region_to_list(self.region)
          artists = [region.to_pixel(wcs=ax.wcs).as_artist() for region in regions]
  
          patches = PatchCollection(artists, **kwargs)
          ax.add_collection(patches)
          return ax
  
      def peek(self, figsize=(5, 10)):
          """Quick-look summary plots."""
          import matplotlib.pyplot as plt
  
          fig, ax = plt.subplots(1, 1, figsize=figsize)
          self.plot_hist(ax=ax)
          return ax
  
      def copy(self):
          """A deep copy of self."""
          return copy.deepcopy(self)
  
      def downsample(self, factor):
          """Downsample spectrum.
  
          Parameters
          ----------
          factor : int
              Downsampling factor.
  
          Returns
          -------
          spectrum : `~gammapy.spectrum.CountsSpectrum`
              Downsampled spectrum.
          """
          from gammapy.extern.skimage import block_reduce
  
          data = block_reduce(self.data, block_size=(factor,))
          energy = self.energy.downsample(factor).edges
          return self.__class__(energy_lo=energy[:-1], energy_hi=energy[1:], data=data)
  
      def energy_mask(self, emin=None, emax=None):
          """Create a mask for a given energy range.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Energy range
          """
          edges = self.energy.edges
  
          # set default values
          emin = emin if emin is not None else edges[0]
          emax = emax if emax is not None else edges[-1]
  
          return (edges[:-1] >= emin) & (edges[1:] <= emax)
  
      def _arithmetics(self, operator, other, copy):
          """Perform arithmetics on maps after checking geometry consistency."""
          if isinstance(other, CountsSpectrum):
              q = other.quantity
          else:
              q = u.Quantity(other, copy=False)
  
          out = self.copy() if copy else self
          out.quantity = operator(out.quantity, q)
          return out
  
      def __add__(self, other):
          return self._arithmetics(np.add, other, copy=True)
  
      def __iadd__(self, other):
          return self._arithmetics(np.add, other, copy=False)
  
      def __sub__(self, other):
          return self._arithmetics(np.subtract, other, copy=True)
  
      def __isub__(self, other):
          return self._arithmetics(np.subtract, other, copy=False)
  
      def __mul__(self, other):
          return self._arithmetics(np.multiply, other, copy=True)
  
      def __imul__(self, other):
          return self._arithmetics(np.multiply, other, copy=False)
  
      def __truediv__(self, other):
          return self._arithmetics(np.true_divide, other, copy=True)
  
      def __itruediv__(self, other):
          return self._arithmetics(np.true_divide, other, copy=False)
  
      def __array__(self):
          return self.data
  
  
  class SpectrumEvaluator:
      """Calculate number of predicted counts (``npred``).
  
      The true and reconstructed energy binning are inferred from the provided IRFs.
  
      Parameters
      ----------
      model : `~gammapy.modeling.models.SkyModel`
          Spectral model
      aeff : `~gammapy.irf.EffectiveAreaTable`
          EffectiveArea
      edisp : `~gammapy.irf.EDispKernel`, optional
          Energy dispersion kernel.
      livetime : `~astropy.units.Quantity`
          Observation duration (may be contained in aeff)
      e_true : `~astropy.units.Quantity`, optional
          Desired energy axis of the prediced counts vector if no IRFs are given
  
      Examples
      --------
      Calculate predicted counts in a desired reconstruced energy binning
  
      .. plot::
          :include-source:
  
          import numpy as np
          import astropy.units as u
          import matplotlib.pyplot as plt
          from gammapy.irf import EffectiveAreaTable, EDispKernel
          from gammapy.modeling.models import PowerLawSpectralModel, SkyModel
          from gammapy.spectrum import SpectrumEvaluator
  
          e_true = np.logspace(-2, 2.5, 109) * u.TeV
          e_reco = np.logspace(-2, 2, 73) * u.TeV
  
          aeff = EffectiveAreaTable.from_parametrization(energy=e_true)
          edisp = EDispKernel.from_gauss(e_true=e_true, e_reco=e_reco, sigma=0.3, bias=0)
  
          pwl = PowerLawSpectralModel(index=2.3, amplitude="2.5e-12 cm-2 s-1 TeV-1", reference="1 TeV")
          model = SkyModel(spectral_model=pwl)
  
          predictor = SpectrumEvaluator(model=model, aeff=aeff, edisp=edisp, livetime="1 hour")
          predictor.compute_npred().plot_hist()
          plt.show()
      """
  
      def __init__(self, model, aeff=None, edisp=None, livetime=None):
          self.model = model
          self.aeff = aeff
          self.edisp = edisp
          self.livetime = livetime
  
      def compute_npred(self):
          e_true = self.aeff.energy.edges
          integral_flux = self.model.spectral_model.integral(
              emin=e_true[:-1], emax=e_true[1:], intervals=True
          )
  
          true_counts = self.apply_aeff(integral_flux)
          return self.apply_edisp(true_counts)
  
      def apply_aeff(self, integral_flux):
          if self.aeff is not None:
              cts = integral_flux * self.aeff.data.data
          else:
              cts = integral_flux
  
          # Multiply with livetime if not already contained in aeff or model
          if cts.unit.is_equivalent("s-1"):
              cts *= self.livetime
  
          return cts.to("")
  
      def apply_edisp(self, true_counts):
          from . import CountsSpectrum
  
          if self.edisp is not None:
              cts = self.edisp.apply(true_counts)
              e_reco = self.edisp.e_reco.edges
          else:
              cts = true_counts
              e_reco = self.aeff.energy.edges
  
          return CountsSpectrum(data=cts, energy_lo=e_reco[:-1], energy_hi=e_reco[1:])
  
files: gammapy/spectrum/core.py