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
  
  [... the rest is omitted ...]
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
  
  [... the rest is omitted ...]
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
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/cube.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Cube models (axes: lon, lat, energy)."""
  import copy
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  from gammapy.maps import Map
  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
  
     
  
  [... the rest is omitted ...]
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)
  
  
  
  [... the rest is omitted ...]
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
  
  [... the rest is omitted ...]
files: gammapy/spectrum/core.py