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, SkyModel, SkyModels
  from gammapy.modeling.parameter import _get_parameters_str
  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.SkyModels`
          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
      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 : `~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,
          background_model=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.evaluation_mode = evaluation_mode
          self.counts = counts
          self.exposure = exposure
          self.mask_fit = mask_fit
          self.psf = psf
          self.edisp = edisp
          self.background_model = background_model
          self.models = models
          self.mask_safe = mask_safe
          self.gti = gti
  
          if name is None:
              self.name = make_name()
          else:
              self.name = name
  
          # check whether a reference geom is defined
          _ = self._geom
  
      def __str__(self):
          str_ = f"{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 or self.background_model 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:
              mask = self.mask_safe.reduce_over_axes(np.logical_or).data
              if not mask.any():
                  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 or self.background_model 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)
  
          if self.background_model is not None:
              n_models += 1
  
          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)
          )
  
          components = []
  
          if self.models is not None:
              components += self.models
  
          if self.background_model is not None:
              components += [self.background_model]
  
          for idx, model in enumerate(components):
              str_ += f"\tComponent {idx}: \n"
              str_ += "\t\t{:28}: {}\n".format("Name", model.name)
              str_ += "\t\t{:28}: {}\n".format("Type", model.__class__.__name__)
  
              if isinstance(model, SkyModel):
                  str_ += "\t\t{:28}: {}\n".format(
                      "Spatial  model type", model.spatial_model.__class__.__name__
                  )
                  str_ += "\t\t{:28}: {}\n".format(
                      "Spectral model type", model.spectral_model.__class__.__name__
                  )
  
              str_ += "\t\tParameters:\n"
              info = _get_parameters_str(model.parameters)
              lines = info.split("\n")
              str_ += "\t\t" + "\n\t\t".join(lines[:-
  
  [... 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 EffectiveAreaTable, EDispKernel
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling import Dataset, Parameters
  from gammapy.modeling.models import BackgroundModel, SkyModel, SkyModels
  from gammapy.modeling.parameter import _get_parameters_str
  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_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.SkyModels`
          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
      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 : `~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,
          background_model=None,
          name="",
          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.evaluation_mode = evaluation_mode
          self.counts = counts
          self.exposure = exposure
          self.mask_fit = mask_fit
          self.psf = psf
          self.edisp = edisp
          self.background_model = background_model
          self.models = models
          self.name = name
          self.mask_safe = mask_safe
          self.gti = gti
  
          # check whether a reference geom is defined
          _ = self._geom
  
      def __str__(self):
          str_ = f"{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 or self.background_model 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:
              mask = self.mask_safe.reduce_over_axes(np.logical_or).data
              if not mask.any():
                  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 or self.background_model 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)
  
          if self.background_model is not None:
              n_models += 1
  
          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)
          )
  
          components = []
  
          if self.models is not None:
              components += self.models
  
          if self.background_model is not None:
              components += [self.background_model]
  
          for idx, model in enumerate(components):
              str_ += f"\tComponent {idx}: \n"
              str_ += "\t\t{:28}: {}\n".format("Name", model.name)
              str_ += "\t\t{:28}: {}\n".format("Type", model.__class__.__name__)
  
              if isinstance(model, SkyModel):
                  str_ += "\t\t{:28}: {}\n".format(
                      "Spatial  model type", model.spatial_model.__class__.__name__
                  )
                  str_ += "\t\t{:28}: {}\n".format(
                      "Spectral model type", model.spectral_model.__class__.__name__
                  )
  
              str_ += "\t\tParameters:\n"
              info = _get_parameters_str(model.parameters)
              lines = info.split("\n")
              str_ += "\t\t" + "\n\t\t".join(lines[:-1])
  
              str_ += "\n\n"
  
          return str_.expandtabs(tabs
  
  [... the rest is omitted ...]
files: gammapy/cube/fit.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import abc
  import collections.abc
  import copy
  from warnings import warn
  import numpy as np
  from gammapy.utils.scripts import make_name, make_path, read_yaml, write_yaml
  from gammapy.utils.table import table_from_row_data
  from ..maps import WcsNDMap
  from .parameter import Parameters
  
  __all__ = ["Dataset", "Datasets"]
  
  
  class Dataset(abc.ABC):
      """Dataset abstract base class.
  
      TODO: add tutorial how to create your own dataset types.
  
      For now, see existing examples in Gammapy how this works:
  
      - `gammapy.cube.MapDataset`
      - `gammapy.spectrum.SpectrumDataset`
      - `gammapy.spectrum.FluxPointsDataset`
      """
  
      _residuals_labels = {
          "diff": "data - model",
          "diff/model": "(data - model) / model",
          "diff/sqrt(model)": "(data - model) / sqrt(model)",
      }
  
      @property
      def mask(self):
          """Combined fit and safe mask"""
          mask_safe = (
              self.mask_safe.data
              if isinstance(self.mask_safe, WcsNDMap)
              else self.mask_safe
          )
          mask_fit = (
              self.mask_fit.data if isinstance(self.mask_fit, WcsNDMap) else self.mask_fit
          )
          if mask_safe is not None and mask_fit is not None:
              mask = mask_safe & mask_fit
          elif mask_fit is not None:
              mask = mask_fit
          elif mask_safe is not None:
              mask = mask_safe
          else:
              mask = None
          return mask
  
      def stat_sum(self):
          """Total statistic given the current model parameters."""
          stat = self.stat_array()
  
          if self.mask is not None:
              stat = stat[self.mask]
  
          return np.sum(stat, dtype=np.float64)
  
      @abc.abstractmethod
      def stat_array(self):
          """Statistic array, one value per data point."""
  
      def copy(self, name=None):
          """A deep copy."""
          new = copy.deepcopy(self)
          if name is None:
              new.name = make_name()
          else:
              new.name = name
          return new
  
      @staticmethod
      def _compute_residuals(data, model, method="diff"):
          with np.errstate(invalid="ignore"):
              if method == "diff":
                  residuals = data - model
              elif method == "diff/model":
                  residuals = (data - model) / model
              elif method == "diff/sqrt(model)":
                  residuals = (data - model) / np.sqrt(model)
              else:
                  raise AttributeError(
                      f"Invalid method: {method!r}. Choose between 'diff',"
                      " 'diff/model' and 'diff/sqrt(model)'"
                  )
          return residuals
  
  
  class Datasets(collections.abc.Sequence):
      """Dataset collection.
  
      Parameters
      ----------
      datasets : `Dataset` or list of `Dataset`
          Datasets
      """
  
      def __init__(self, datasets):
          if isinstance(datasets, Datasets):
              datasets = list(datasets)
          elif isinstance(datasets, list):
              pass
          else:
              raise TypeError(f"Invalid type: {datasets!r}")
  
          unique_names = []
          for dataset in datasets:
              while dataset.name in unique_names:
                  dataset.name = make_name()  # replace duplicate
                  warn("Dateset names must be unique, auto-replaced duplicates")
              unique_names.append(dataset.name)
  
          self._datasets = datasets
  
      @property
      def parameters(self):
          """Unique parameters (`~gammapy.modeling.Parameters`).
  
          Duplicate parameter objects have been removed.
          The order of the unique parameters remains.
          """
          parameters = Parameters.from_stack(_.parameters for _ in self)
          return parameters.unique_parameters
  
      @property
      def is_all_same_type(self):
          """Whether all contained datasets are of the same type"""
          return len(set(_.__class__ for _ in self)) == 1
  
      @property
      def is_all_same_shape(self):
          """Whether all contained datasets have the same data shape"""
          return len(set(_.data_shape for _ in self)) == 1
  
      def stat_sum(self):
          """Compute joint likelihood"""
          stat_sum = 0
          # TODO: add parallel evaluation of likelihoods
          for dataset in self:
              stat_sum += dataset.stat_sum()
          return stat_sum
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n"
          str_ += "--------\n"
  
          for idx, dataset in enumerate(self):
              str_ += f"idx={idx}, id={hex(id(dataset))!r}, name={dataset.name!r}\n"
  
          return str_
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      @classmethod
      def read(cls, filedata, filemodel):
          """De-serialize datasets from YAML and FITS files.
  
          Parameters
          ----------
          filedata : str
              filepath to yaml datasets file
          filemodel : str
              filepath to yaml models file
  
          Returns
          -------
          dataset : 'gammapy.modeling.Datasets'
              Datasets
          """
          from .serialize import dict_to_datasets
  
          components = read_yaml(make_path(filemodel))
          data_list = read_yaml(make_path(filedata))
          datasets = dict_to_datasets(data_list, components)
          return cls(datasets)
  
      def write(self, path, prefix="", overwrite=False):
          """Serialize datasets to YAML and FITS files.
  
          Parameters
          ----------
          path : `pathlib.Path`
              path to write files
          prefix : str
              common prefix of file names
          overwrite : bool
              overwrite datasets FITS files
          """
          from .serialize import datasets_to_dict
  
          path = make_path(path)
  
          datasets_dict, components_dict = datasets_to_dict(self, path, prefix, overwrite)
          write_yaml(datasets_dict, path / f"{prefix}_datasets.yaml", sort_keys=False)
          write_yaml(components_dict, path / f"{prefix}_models.yaml", sort_keys=False)
  
      def stack_reduce(self):
          """Reduce the Datasets to a unique Dataset by stacking them together.
  
          This works only if all Dataset are of the same type and if a proper
          in-place stack method exists for the Dataset type.
  
          Returns
          -------
          dataset : ~gammapy.utils.Dataset
              the stacked dataset
          """
          if not self.is_all_same_type:
              raise ValueError(
                  "Stacking impossible: all Datasets contained are not of a unique type."
              )
  
          dataset = self[0].copy()
          for ds in self[1:]:
              dataset.stack(ds)
          return dataset
  
      def info_table(self, cumulative=False):
          """Get info table for datasets.
  
          Parameters
          ----------
          cumulative : bool
              Cumulate info across all observations
  
          Returns
          -------
          info_table : `~astropy.table.Table`
              Info table.
          """
          if not self.is_all_same_type:
              raise ValueError("Info table not supported for mixed dataset type.")
  
          stacked = self[0].copy()
  
          rows = [stacked.info_dict()]
  
          for dataset in self[1:]:
              if cumulative:
                  stacked.stack(dataset)
                  row = stacked.info_dict()
              else:
                  row = dataset.info_dict()
  
              rows.append(row)
  
          return table_from_row_data(rows=rows)
  
      def __getitem__(self, val):
          if isinstance(val, (int, slice)):
              return self._datasets[val]
          elif isinstance(val, str):
              for idx, dataset in enumerate(self._datasets):
                  if val == dataset.name:
                      return self._datasets[idx]
              raise IndexError(f"No dataset: {val!r}")
          else:
              raise TypeError(f"Invalid type: {type(val)!r}")
  
      def __len__(self):
          return len(self._datasets)
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/datasets.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import abc
  import collections.abc
  import copy
  import numpy as np
  from gammapy.utils.scripts import make_path, read_yaml, write_yaml
  from gammapy.utils.table import table_from_row_data
  from ..maps import WcsNDMap
  from .parameter import Parameters
  
  __all__ = ["Dataset", "Datasets"]
  
  
  class Dataset(abc.ABC):
      """Dataset abstract base class.
  
      TODO: add tutorial how to create your own dataset types.
  
      For now, see existing examples in Gammapy how this works:
  
      - `gammapy.cube.MapDataset`
      - `gammapy.spectrum.SpectrumDataset`
      - `gammapy.spectrum.FluxPointsDataset`
      """
  
      _residuals_labels = {
          "diff": "data - model",
          "diff/model": "(data - model) / model",
          "diff/sqrt(model)": "(data - model) / sqrt(model)",
      }
  
      @property
      def mask(self):
          """Combined fit and safe mask"""
          mask_safe = (
              self.mask_safe.data
              if isinstance(self.mask_safe, WcsNDMap)
              else self.mask_safe
          )
          mask_fit = (
              self.mask_fit.data if isinstance(self.mask_fit, WcsNDMap) else self.mask_fit
          )
          if mask_safe is not None and mask_fit is not None:
              mask = mask_safe & mask_fit
          elif mask_fit is not None:
              mask = mask_fit
          elif mask_safe is not None:
              mask = mask_safe
          else:
              mask = None
          return mask
  
      def stat_sum(self):
          """Total statistic given the current model parameters."""
          stat = self.stat_array()
  
          if self.mask is not None:
              stat = stat[self.mask]
  
          return np.sum(stat, dtype=np.float64)
  
      @abc.abstractmethod
      def stat_array(self):
          """Statistic array, one value per data point."""
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      @staticmethod
      def _compute_residuals(data, model, method="diff"):
          with np.errstate(invalid="ignore"):
              if method == "diff":
                  residuals = data - model
              elif method == "diff/model":
                  residuals = (data - model) / model
              elif method == "diff/sqrt(model)":
                  residuals = (data - model) / np.sqrt(model)
              else:
                  raise AttributeError(
                      f"Invalid method: {method!r}. Choose between 'diff',"
                      " 'diff/model' and 'diff/sqrt(model)'"
                  )
          return residuals
  
  
  class Datasets(collections.abc.Sequence):
      """Dataset collection.
  
      Parameters
      ----------
      datasets : `Dataset` or list of `Dataset`
          Datasets
      """
  
      def __init__(self, datasets):
          if isinstance(datasets, Datasets):
              self._datasets = list(datasets)
          elif isinstance(datasets, list):
              self._datasets = datasets
          else:
              raise TypeError(f"Invalid type: {datasets!r}")
  
      @property
      def parameters(self):
          """Unique parameters (`~gammapy.modeling.Parameters`).
  
          Duplicate parameter objects have been removed.
          The order of the unique parameters remains.
          """
          parameters = Parameters.from_stack(_.parameters for _ in self)
          return parameters.unique_parameters
  
      @property
      def is_all_same_type(self):
          """Whether all contained datasets are of the same type"""
          return len(set(_.__class__ for _ in self)) == 1
  
      @property
      def is_all_same_shape(self):
          """Whether all contained datasets have the same data shape"""
          return len(set(_.data_shape for _ in self)) == 1
  
      def stat_sum(self):
          """Compute joint likelihood"""
          stat_sum = 0
          # TODO: add parallel evaluation of likelihoods
          for dataset in self:
              stat_sum += dataset.stat_sum()
          return stat_sum
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n"
          str_ += "--------\n"
  
          for idx, dataset in enumerate(self):
              str_ += f"idx={idx}, id={hex(id(dataset))!r}, name={dataset.name!r}\n"
  
          return str_
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      @classmethod
      def read(cls, filedata, filemodel):
          """De-serialize datasets from YAML and FITS files.
  
          Parameters
          ----------
          filedata : str
              filepath to yaml datasets file
          filemodel : str
              filepath to yaml models file
  
          Returns
          -------
          dataset : 'gammapy.modeling.Datasets'
              Datasets
          """
          from .serialize import dict_to_datasets
  
          components = read_yaml(make_path(filemodel))
          data_list = read_yaml(make_path(filedata))
          datasets = dict_to_datasets(data_list, components)
          return cls(datasets)
  
      def write(self, path, prefix="", overwrite=False):
          """Serialize datasets to YAML and FITS files.
  
          Parameters
          ----------
          path : `pathlib.Path`
              path to write files
          prefix : str
              common prefix of file names
          overwrite : bool
              overwrite datasets FITS files
          """
          from .serialize import datasets_to_dict
  
          path = make_path(path)
  
          datasets_dict, components_dict = datasets_to_dict(self, path, prefix, overwrite)
          write_yaml(datasets_dict, path / f"{prefix}_datasets.yaml", sort_keys=False)
          write_yaml(components_dict, path / f"{prefix}_models.yaml", sort_keys=False)
  
      def stack_reduce(self):
          """Reduce the Datasets to a unique Dataset by stacking them together.
  
          This works only if all Dataset are of the same type and if a proper
          in-place stack method exists for the Dataset type.
  
          Returns
          -------
          dataset : ~gammapy.utils.Dataset
              the stacked dataset
          """
          if not self.is_all_same_type:
              raise ValueError(
                  "Stacking impossible: all Datasets contained are not of a unique type."
              )
  
          dataset = self[0].copy()
          for ds in self[1:]:
              dataset.stack(ds)
          return dataset
  
      def info_table(self, cumulative=False):
          """Get info table for datasets.
  
          Parameters
          ----------
          cumulative : bool
              Cumulate info across all observations
  
          Returns
          -------
          info_table : `~astropy.table.Table`
              Info table.
          """
          if not self.is_all_same_type:
              raise ValueError("Info table not supported for mixed dataset type.")
  
          stacked = self[0].copy()
  
          rows = [stacked.info_dict()]
  
          for dataset in self[1:]:
              if cumulative:
                  stacked.stack(dataset)
                  row = stacked.info_dict()
              else:
                  row = dataset.info_dict()
  
              rows.append(row)
  
          return table_from_row_data(rows=rows)
  
      def __getitem__(self, val):
          if isinstance(val, (int, slice)):
              return self._datasets[val]
          elif isinstance(val, str):
              for idx, dataset in enumerate(self._datasets):
                  if val == dataset.name:
                      return self._datasets[idx]
              raise IndexError(f"No dataset: {val!r}")
          else:
              raise TypeError(f"Invalid type: {type(val)!r}")
  
      def __len__(self):
          return len(self._datasets)
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/datasets.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import copy
  import astropy.units as u
  from gammapy.utils.scripts import make_name
  from .parameter import Parameter, Parameters
  
  __all__ = ["Model"]
  
  
  class Model:
      """Model base class."""
  
      def __init__(self, **kwargs):
          # Copy default parameters from the class to the instance
          self._parameters = self.__class__.default_parameters.copy()
          for parameter in self._parameters:
              if parameter.name in self.__dict__:
                  raise ValueError(
                      f"Invalid parameter name: {parameter.name!r}."
                      f"Attribute exists already: {getattr(self, parameter.name)!r}"
                  )
  
              setattr(self, parameter.name, parameter)
  
          # Update parameter information from kwargs
          for name, value in kwargs.items():
              if name not in self.parameters.names:
                  raise ValueError(
                      f"Invalid argument: {name!r}. Parameter names are: {self.parameters.names}"
                  )
  
              self._parameters[name].quantity = u.Quantity(value)
  
      def __init_subclass__(cls, **kwargs):
          # Add parameters list on the model sub-class (not instances)
          cls.default_parameters = Parameters(
              [_ for _ in cls.__dict__.values() if isinstance(_, Parameter)]
          )
  
      def _init_from_parameters(self, parameters):
          """Create model from list of parameters.
  
          This should be called for models that generate
          the parameters dynamically in ``__init__``,
          like the ``NaimaSpectralModel``
          """
          # TODO: should we pass through `Parameters` here? Why?
          parameters = Parameters(parameters)
          self._parameters = parameters
          for parameter in parameters:
              setattr(self, parameter.name, parameter)
  
      @property
      def parameters(self):
          """Parameters (`~gammapy.modeling.Parameters`)"""
          return self._parameters
  
      def copy(self, name=None):
          """A deep copy."""
          new = copy.deepcopy(self)
          if name is None:
              new.name = make_name()
          else:
              new.name = name
          return new
  
      def __str__(self):
          return f"{self.__class__.__name__}\n\n{self.parameters.to_table()}"
  
      def to_dict(self):
          """Create dict for YAML serialisation"""
          return {"type": self.tag, "parameters": self.parameters.to_dict()["parameters"]}
  
      @classmethod
      def from_dict(cls, data):
          params = {
              x["name"].split("@")[0]: x["value"] * u.Unit(x["unit"])
              for x in data["parameters"]
          }
  
          # TODO: this is a special case for spatial models, maybe better move to `SpatialModel` base class
          if "frame" in data:
              params["frame"] = data["frame"]
  
          model = cls(**params)
          model._update_from_dict(data)
          return model
  
      # TODO: try to get rid of this
      def _update_from_dict(self, data):
          self._parameters.update_from_dict(data)
          for parameter in self.parameters:
              setattr(self, parameter.name, parameter)
  
      @staticmethod
      def create(tag, *args, **kwargs):
          """Create a model instance.
  
          Examples
          --------
          >>> from gammapy.modeling import Model
          >>> spectral_model = Model.create("PowerLaw2SpectralModel", amplitude="1e-10 cm-2 s-1", index=3)
          >>> type(spectral_model)
          gammapy.modeling.models.spectral.PowerLaw2SpectralModel
          """
          from .models import MODELS
  
          cls = MODELS.get_cls(tag)
          return cls(*args, **kwargs)
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/model.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import copy
  import astropy.units as u
  from .parameter import Parameter, Parameters
  
  __all__ = ["Model"]
  
  
  class Model:
      """Model base class."""
  
      def __init__(self, **kwargs):
          # Copy default parameters from the class to the instance
          self._parameters = self.__class__.default_parameters.copy()
          for parameter in self._parameters:
              if parameter.name in self.__dict__:
                  raise ValueError(
                      f"Invalid parameter name: {parameter.name!r}."
                      f"Attribute exists already: {getattr(self, parameter.name)!r}"
                  )
  
              setattr(self, parameter.name, parameter)
  
          # Update parameter information from kwargs
          for name, value in kwargs.items():
              if name not in self.parameters.names:
                  raise ValueError(
                      f"Invalid argument: {name!r}. Parameter names are: {self.parameters.names}"
                  )
  
              self._parameters[name].quantity = u.Quantity(value)
  
      def __init_subclass__(cls, **kwargs):
          # Add parameters list on the model sub-class (not instances)
          cls.default_parameters = Parameters(
              [_ for _ in cls.__dict__.values() if isinstance(_, Parameter)]
          )
  
      def _init_from_parameters(self, parameters):
          """Create model from list of parameters.
  
          This should be called for models that generate
          the parameters dynamically in ``__init__``,
          like the ``NaimaSpectralModel``
          """
          # TODO: should we pass through `Parameters` here? Why?
          parameters = Parameters(parameters)
          self._parameters = parameters
          for parameter in parameters:
              setattr(self, parameter.name, parameter)
  
      @property
      def parameters(self):
          """Parameters (`~gammapy.modeling.Parameters`)"""
          return self._parameters
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      def __str__(self):
          return f"{self.__class__.__name__}\n\n{self.parameters.to_table()}"
  
      def to_dict(self):
          """Create dict for YAML serialisation"""
          return {"type": self.tag, "parameters": self.parameters.to_dict()["parameters"]}
  
      @classmethod
      def from_dict(cls, data):
          params = {
              x["name"].split("@")[0]: x["value"] * u.Unit(x["unit"])
              for x in data["parameters"]
          }
  
          # TODO: this is a special case for spatial models, maybe better move to `SpatialModel` base class
          if "frame" in data:
              params["frame"] = data["frame"]
  
          model = cls(**params)
          model._update_from_dict(data)
          return model
  
      # TODO: try to get rid of this
      def _update_from_dict(self, data):
          self._parameters.update_from_dict(data)
          for parameter in self.parameters:
              setattr(self, parameter.name, parameter)
  
      @staticmethod
      def create(tag, *args, **kwargs):
          """Create a model instance.
  
          Examples
          --------
          >>> from gammapy.modeling import Model
          >>> spectral_model = Model.create("PowerLaw2SpectralModel", amplitude="1e-10 cm-2 s-1", index=3)
          >>> type(spectral_model)
          gammapy.modeling.models.spectral.PowerLaw2SpectralModel
          """
          from .models import MODELS
  
          cls = MODELS.get_cls(tag)
          return cls(*args, **kwargs)
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/model.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Cube models (axes: lon, lat, energy)."""
  import collections.abc
  import copy
  from pathlib import Path
  from warnings import warn
  import numpy as np
  import astropy.units as u
  import yaml
  from gammapy.maps import Map
  from gammapy.modeling import Model, Parameter, Parameters
  from gammapy.utils.scripts import make_name, make_path
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      def __add__(self, other):
          if isinstance(other, (SkyModels, list)):
              return SkyModels([self, *other])
          elif isinstance(other, (SkyModel, SkyDiffuseCube)):
              return SkyModels([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 SkyModels(collections.abc.Sequence):
      """Sky model collection.
  
      Parameters
      ----------
      skymodels : `SkyModel`, list of `SkyModel` or `SkyModels`
          Sky models
      """
  
      def __init__(self, skymodels):
          if isinstance(skymodels, SkyModels):
              models = skymodels._skymodels
          elif isinstance(skymodels, SkyModel):
              models = [skymodels]
          elif isinstance(skymodels, list):
              models = skymodels
          else:
              raise TypeError(f"Invalid type: {skymodels!r}")
  
          unique_names = []
          for model in models:
              while model.name in unique_names:
                  model.name = make_name()  # replace duplicate
                  warn("SkyModel names must be unique, auto-replaced duplicates")
              unique_names.append(model.name)
  
          self._skymodels = models
  
      @property
      def parameters(self):
          return Parameters.from_stack([_.parameters for _ in self._skymodels])
  
      @classmethod
      def read(cls, filename):
          """Read from YAML file."""
          yaml_str = Path(filename).read_text()
          return cls.from_yaml(yaml_str)
  
      @classmethod
      def from_yaml(cls, yaml_str):
          """Create from YAML string."""
          from gammapy.modeling.serialize import dict_to_models
  
          data = yaml.safe_load(yaml_str)
          skymodels = dict_to_models(data)
          return cls(skymodels)
  
      def write(self, path, overwrite=False):
          """Write to YAML file."""
          path = make_path(path)
          if path.exists() and not overwrite:
              raise IOError(f"File exists already: {path}")
          path.write_text(self.to_yaml())
  
      def to_yaml(self):
          """Convert to YAML string."""
          from gammapy.modeling.serialize import models_to_dict
  
          data = models_to_dict(self._skymodels)
          return yaml.dump(
              data, sort_keys=False, indent=4, width=80, default_flow_style=None
          )
  
      def __str__(self):
          str_ = f"{self.__class__.__name__}\n\n"
  
          for idx, skymodel in enumerate(self):
              str_ += f"Component {idx}: {skymodel}\n\n\t\n\n"
  
          return str_
  
      def __add__(self, other):
          if isinstance(other, (SkyModels, list)):
              return SkyModels([*self, *other])
          elif isinstance(other, (SkyModel, SkyDiffuseCube)):
              return SkyModels([*self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __getitem__(self, val):
          if isinstance(val, int):
              return self._skymodels[val]
          elif isinstance(val, str):
              for idx, model in enumerate(self._skymodels):
                  if val == model.name:
                      return self._skymodels[idx]
              raise IndexError(f"No model: {val!r}")
          else:
              raise TypeError(f"Invalid type: {type(val)!r}")
  
      def __len__(self):
          return len(self._skymodels)
  
  
  class SkyModel(SkyModelBase):
      """Sky model component.
  
      This model represents a factorised sky model.
      It has `~gammapy.modeling.Parameters`
      combining the spatial and spectral parameters.
  
      TODO: add possibility to have a temporal model component also.
  
      Parameters
      ----------
      spectral_model : `~gammapy.modeling.models.SpectralModel`
          Spectral model
      spatial_model : `~gammapy.modeling.models.SpatialModel`
          Spatial model (must be normalised to integrate to 1)
      name : str
          Model identifier
      """
  
      tag = "SkyModel"
  
      def __init__(self, spectral_model, spatial_model=None, name=None):
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          super().__init__()
          # TODO: this hack is needed for compound models to work
          self.__dict__.pop("_parameters")
  
          if name is None:
              self.name = make_name()
          else:
              self.name = 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 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})"
          )
  
      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
  
          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)
  
          if self.spatial_model is not None:
              value = value * self.spatial_model.evaluate_geom(geom.to_image())
  
          return value
  
      def copy(self, **kwargs):
          """Copy SkyModel"""
          if self.spatial_model is not None:
              spatial_model = self.spatial_model.copy()
          else:
              spatial_model = None
  
          kwargs.setdefault("spatial_model", spatial_model)
          kwargs.setdefault("spectral_model", self.spectral_model.copy())
          kwargs.setdefault("name", make_name())
          return self.__class__(**kwargs)
  
      def to_dict(self):
          """Create dict for YAML serilisation"""
          data = {}
          data["name"] = self.name
          data["type"] = self.tag
          data
  
  [... 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 collections.abc
  import copy
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  import yaml
  from gammapy.maps import Map
  from gammapy.modeling import Model, Parameter, Parameters
  from gammapy.utils.scripts import make_path
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      def __add__(self, other):
          if isinstance(other, (SkyModels, list)):
              return SkyModels([self, *other])
          elif isinstance(other, (SkyModel, SkyDiffuseCube)):
              return SkyModels([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 SkyModels(collections.abc.Sequence):
      """Sky model collection.
  
      Parameters
      ----------
      skymodels : `SkyModel`, list of `SkyModel` or `SkyModels`
          Sky models
      """
  
      def __init__(self, skymodels):
          if isinstance(skymodels, SkyModels):
              models = skymodels._skymodels
          elif isinstance(skymodels, SkyModel):
              models = [skymodels]
          elif isinstance(skymodels, list):
              models = skymodels
          else:
              raise TypeError(f"Invalid type: {skymodels!r}")
  
          self._skymodels = models
  
      @property
      def parameters(self):
          return Parameters.from_stack([_.parameters for _ in self._skymodels])
  
      @classmethod
      def read(cls, filename):
          """Read from YAML file."""
          yaml_str = Path(filename).read_text()
          return cls.from_yaml(yaml_str)
  
      @classmethod
      def from_yaml(cls, yaml_str):
          """Create from YAML string."""
          from gammapy.modeling.serialize import dict_to_models
  
          data = yaml.safe_load(yaml_str)
          skymodels = dict_to_models(data)
          return cls(skymodels)
  
      def write(self, path, overwrite=False):
          """Write to YAML file."""
          path = make_path(path)
          if path.exists() and not overwrite:
              raise IOError(f"File exists already: {path}")
          path.write_text(self.to_yaml())
  
      def to_yaml(self):
          """Convert to YAML string."""
          from gammapy.modeling.serialize import models_to_dict
  
          data = models_to_dict(self._skymodels)
          return yaml.dump(
              data, sort_keys=False, indent=4, width=80, default_flow_style=None
          )
  
      def __str__(self):
          str_ = f"{self.__class__.__name__}\n\n"
  
          for idx, skymodel in enumerate(self):
              str_ += f"Component {idx}: {skymodel}\n\n\t\n\n"
  
          return str_
  
      def __add__(self, other):
          if isinstance(other, (SkyModels, list)):
              return SkyModels([*self, *other])
          elif isinstance(other, (SkyModel, SkyDiffuseCube)):
              return SkyModels([*self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __getitem__(self, val):
          if isinstance(val, int):
              return self._skymodels[val]
          elif isinstance(val, str):
              for idx, model in enumerate(self._skymodels):
                  if val == model.name:
                      return self._skymodels[idx]
              raise IndexError(f"No model: {val!r}")
          else:
              raise TypeError(f"Invalid type: {type(val)!r}")
  
      def __len__(self):
          return len(self._skymodels)
  
  
  class SkyModel(SkyModelBase):
      """Sky model component.
  
      This model represents a factorised sky model.
      It has `~gammapy.modeling.Parameters`
      combining the spatial and spectral parameters.
  
      TODO: add possibility to have a temporal model component also.
  
      Parameters
      ----------
      spectral_model : `~gammapy.modeling.models.SpectralModel`
          Spectral model
      spatial_model : `~gammapy.modeling.models.SpatialModel`
          Spatial model (must be normalised to integrate to 1)
      name : str
          Model identifier
      """
  
      tag = "SkyModel"
  
      def __init__(self, spectral_model, spatial_model=None, name="source"):
          self.name = name
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          super().__init__()
          # TODO: this hack is needed for compound models to work
          self.__dict__.pop("_parameters")
  
      @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 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})"
          )
  
      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
  
          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)
  
          if self.spatial_model is not None:
              value = value * self.spatial_model.evaluate_geom(geom.to_image())
  
          return value
  
      def copy(self, **kwargs):
          """Copy SkyModel"""
          if self.spatial_model is not None:
              spatial_model = self.spatial_model.copy()
          else:
              spatial_model = None
  
          kwargs.setdefault("spatial_model", spatial_model)
          kwargs.setdefault("spectral_model", self.spectral_model.copy())
          kwargs.setdefault("name", self.name + "-copy")
          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()
  
          return data
  
      @classmethod
      def from_dict(cls, data):
          """Create SkyModel from dict"""
          from gammapy.modeling.models import SP
  
  [... the rest is omitted ...]
files: gammapy/modeling/models/cube.py
files: """Create example model YAML files programmatically.
  
  (some will be also written manually)
  """
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  from gammapy.cube import MapDataset, MapDatasetMaker
  from gammapy.data import DataStore
  from gammapy.maps import MapAxis, WcsGeom
  from gammapy.modeling import Datasets
  from gammapy.modeling.models import (
      ExpCutoffPowerLawSpectralModel,
      GaussianSpatialModel,
      PointSpatialModel,
      PowerLawSpectralModel,
      SkyDiffuseCube,
      SkyModel,
      SkyModels,
  )
  
  DATA_PATH = Path("./")
  
  
  def make_example_2():
      spatial = GaussianSpatialModel(lon_0="0 deg", lat_0="0 deg", sigma="1 deg")
      model = SkyModel(PowerLawSpectralModel(), spatial)
      models = SkyModels([model])
      models.write(DATA_PATH / "example2.yaml")
  
  
  def make_datasets_example():
      # Define which data to use and print some information
  
      energy_axis = MapAxis.from_edges(
          np.logspace(-1.0, 1.0, 4), unit="TeV", name="energy", interp="log"
      )
      geom0 = WcsGeom.create(
          skydir=(0, 0),
          binsz=0.1,
          width=(1, 1),
          frame="galactic",
          proj="CAR",
          axes=[energy_axis],
      )
      geom1 = WcsGeom.create(
          skydir=(1, 0),
          binsz=0.1,
          width=(1, 1),
          frame="galactic",
          proj="CAR",
          axes=[energy_axis],
      )
      geoms = [geom0, geom1]
  
      sources_coords = [(0, 0), (0.9, 0.1)]
      names = ["gc", "g09"]
      models = []
  
      for idx, (lon, lat) in enumerate(sources_coords):
          spatial_model = PointSpatialModel(
              lon_0=lon * u.deg, lat_0=lat * u.deg, frame="galactic"
          )
          spectral_model = ExpCutoffPowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=3e-12 * u.Unit("cm-2 s-1 TeV-1"),
              reference=1.0 * u.TeV,
              lambda_=0.1 / u.TeV,
          )
          model_ecpl = SkyModel(
              spatial_model=spatial_model, spectral_model=spectral_model, name=names[idx]
          )
          models.append(model_ecpl)
  
      # test to link a spectral parameter
      params0 = models[0].spectral_model.parameters
      params1 = models[1].spectral_model.parameters
      params0.link("reference", params1["reference"])
      # update the sky model
      models[0].parameters.link("reference", params1["reference"])
  
      obs_ids = [110380, 111140, 111159]
      data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
  
      diffuse_model = SkyDiffuseCube.read(
          "$GAMMAPY_DATA/fermi_3fhl/gll_iem_v06_cutout.fits"
      )
  
      datasets_list = []
      for idx, geom in enumerate(geoms):
          observations = data_store.get_observations(obs_ids)
  
          stacked = MapDataset.create(geom=geom)
          stacked.background_model.name = "background_irf_" + names[idx]
  
          maker = MapDatasetMaker(offset_max=4.0 * u.deg)
  
          for obs in observations:
              dataset = maker.run(stacked, obs)
              stacked.stack(dataset)
  
          stacked.psf = stacked.psf.get_psf_kernel(
              position=geom.center_skydir, geom=geom, max_radius="0.3 deg"
          )
  
          stacked.name = names[idx]
          stacked.models = models[idx] + diffuse_model
          datasets_list.append(stacked)
  
      datasets = Datasets(datasets_list)
  
      dataset0 = datasets[0]
      print("dataset0")
      print("counts sum : ", dataset0.counts.data.sum())
      print("expo sum : ", dataset0.exposure.data.sum())
      print("bkg0 sum : ", dataset0.background_model.evaluate().data.sum())
  
      datasets.write("$GAMMAPY_DATA/tests/models", prefix="gc_example_", overwrite=True)
  
  
  if __name__ == "__main__":
      make_example_2()
      make_datasets_example()
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/tests/data/make.py
files: """Create example model YAML files programmatically.
  
  (some will be also written manually)
  """
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  from gammapy.cube import MapDataset, MapDatasetMaker
  from gammapy.data import DataStore
  from gammapy.maps import MapAxis, WcsGeom
  from gammapy.modeling import Datasets
  from gammapy.modeling.models import (
      ExpCutoffPowerLawSpectralModel,
      GaussianSpatialModel,
      PointSpatialModel,
      PowerLawSpectralModel,
      SkyDiffuseCube,
      SkyModel,
      SkyModels,
  )
  
  DATA_PATH = Path("./")
  
  
  def make_example_2():
      spatial = GaussianSpatialModel(lon_0="0 deg", lat_0="0 deg", sigma="1 deg")
      model = SkyModel(spatial, PowerLawSpectralModel())
      models = SkyModels([model])
      models.write(DATA_PATH / "example2.yaml")
  
  
  def make_datasets_example():
      # Define which data to use and print some information
  
      energy_axis = MapAxis.from_edges(
          np.logspace(-1.0, 1.0, 4), unit="TeV", name="energy", interp="log"
      )
      geom0 = WcsGeom.create(
          skydir=(0, 0),
          binsz=0.1,
          width=(1, 1),
          frame="galactic",
          proj="CAR",
          axes=[energy_axis],
      )
      geom1 = WcsGeom.create(
          skydir=(1, 0),
          binsz=0.1,
          width=(1, 1),
          frame="galactic",
          proj="CAR",
          axes=[energy_axis],
      )
      geoms = [geom0, geom1]
  
      sources_coords = [(0, 0), (0.9, 0.1)]
      names = ["gc", "g09"]
      models = []
  
      for idx, (lon, lat) in enumerate(sources_coords):
          spatial_model = PointSpatialModel(
              lon_0=lon * u.deg, lat_0=lat * u.deg, frame="galactic"
          )
          spectral_model = ExpCutoffPowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=3e-12 * u.Unit("cm-2 s-1 TeV-1"),
              reference=1.0 * u.TeV,
              lambda_=0.1 / u.TeV,
          )
          model_ecpl = SkyModel(
              spatial_model=spatial_model, spectral_model=spectral_model, name=names[idx]
          )
          models.append(model_ecpl)
  
      # test to link a spectral parameter
      params0 = models[0].spectral_model.parameters
      params1 = models[1].spectral_model.parameters
      params0.link("reference", params1["reference"])
      # update the sky model
      models[0].parameters.link("reference", params1["reference"])
  
      obs_ids = [110380, 111140, 111159]
      data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
  
      diffuse_model = SkyDiffuseCube.read(
          "$GAMMAPY_DATA/fermi_3fhl/gll_iem_v06_cutout.fits"
      )
  
      datasets_list = []
      for idx, geom in enumerate(geoms):
          observations = data_store.get_observations(obs_ids)
  
          stacked = MapDataset.create(geom=geom)
          stacked.background_model.name = "background_irf_" + names[idx]
  
          maker = MapDatasetMaker(offset_max=4.0 * u.deg)
  
          for obs in observations:
              dataset = maker.run(stacked, obs)
              stacked.stack(dataset)
  
          stacked.psf = stacked.psf.get_psf_kernel(
              position=geom.center_skydir, geom=geom, max_radius="0.3 deg"
          )
  
          stacked.name = names[idx]
          stacked.models = models[idx] + diffuse_model
          datasets_list.append(stacked)
  
      datasets = Datasets(datasets_list)
  
      dataset0 = datasets[0]
      print("dataset0")
      print("counts sum : ", dataset0.counts.data.sum())
      print("expo sum : ", dataset0.exposure.data.sum())
      print("bkg0 sum : ", dataset0.background_model.evaluate().data.sum())
  
      datasets.write("$GAMMAPY_DATA/tests/models", prefix="gc_example_", overwrite=True)
  
  
  if __name__ == "__main__":
      make_example_2()
      make_datasets_example()
  
  
  [... the rest is omitted ...]
files: gammapy/modeling/tests/data/make.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from pathlib import Path
  import numpy as np
  from astropy import units as u
  from astropy.io import fits
  from astropy.table import Table
  from gammapy.data import GTI
  from gammapy.irf import EDispKernel, EffectiveAreaTable, IRFStacker
  from gammapy.modeling import Dataset, Parameters
  from gammapy.modeling.models import SkyModel, SkyModels
  from gammapy.stats import cash, significance_on_off, significance, wstat
  from gammapy.utils.fits import energy_axis_to_ebounds
  from gammapy.utils.random import get_random_state
  from gammapy.utils.scripts import make_name, make_path
  from .core import CountsSpectrum, SpectrumEvaluator
  
  __all__ = [
      "SpectrumDatasetOnOff",
      "SpectrumDataset",
      "plot_spectrum_datasets_off_regions",
  ]
  
  
  class SpectrumDataset(Dataset):
      """Spectrum dataset for likelihood fitting.
  
      The spectrum dataset bundles reduced counts data, with a spectral model,
      background model and instrument response function to compute the fit-statistic
      given the current model and data.
  
      Parameters
      ----------
      models : `~gammapy.modeling.models.SkyModels`
          Fit model
      counts : `~gammapy.spectrum.CountsSpectrum`
          Counts spectrum
      livetime : `~astropy.units.Quantity`
          Livetime
      aeff : `~gammapy.irf.EffectiveAreaTable`
          Effective area
      edisp : `~gammapy.irf.EDispKernel`
          Energy dispersion
      background : `~gammapy.spectrum.CountsSpectrum`
          Background to use for the fit.
      mask_safe : `~numpy.ndarray`
          Mask defining the safe data range.
      mask_fit : `~numpy.ndarray`
          Mask to apply to the likelihood for fitting.
      name : str
          Dataset name.
      gti : `~gammapy.data.GTI`
          GTI of the observation or union of GTI if it is a stacked observation
  
      See Also
      --------
      SpectrumDatasetOnOff, FluxPointsDataset, gammapy.cube.MapDataset
      """
  
      stat_type = "cash"
      tag = "SpectrumDataset"
  
      def __init__(
          self,
          models=None,
          counts=None,
          livetime=None,
          aeff=None,
          edisp=None,
          background=None,
          mask_safe=None,
          mask_fit=None,
          name=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.counts = counts
  
          if livetime is not None:
              livetime = u.Quantity(livetime)
  
          self.livetime = livetime
          self.mask_fit = mask_fit
          self.aeff = aeff
          self.edisp = edisp
          self.background = background
          self.models = models
          self.mask_safe = mask_safe
          self.gti = gti
  
          if name is None:
              self.name = make_name()
          else:
              self.name = name
  
      def __str__(self):
          str_ = self.__class__.__name__
          str_ += "\n\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)
  
          counts_off = np.nan
          if getattr(self, "counts_off", None) is not None:
              counts_off = np.sum(self.counts_off.data)
              str_ += "\t{:32}: {:.2f}\n\n".format("Total off counts", counts_off)
  
          background = np.nan
          if getattr(self, "background", None) is not None:
              background = np.sum(self.background.data)
              str_ += "\t{:32}: {:.2f}\n\n".format("Total background counts", background)
  
          aeff_min, aeff_max, aeff_unit = np.nan, np.nan, ""
          if self.aeff is not None:
              aeff_min = np.min(self.aeff.data.data.value[self.aeff.data.data.value > 0])
              aeff_max = np.max(self.aeff.data.data.value)
              aeff_unit = self.aeff.data.data.unit
  
          str_ += "\t{:32}: {:.2e} {}\n".format("Effective area min", aeff_min, aeff_unit)
          str_ += "\t{:32}: {:.2e} {}\n\n".format(
              "Effective area max", aeff_max, aeff_unit
          )
  
          livetime = np.nan
          if self.livetime is not None:
              livetime = self.livetime
          str_ += "\t{:32}: {:.2e}\n\n".format("Livetime", livetime)
  
          # 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)
          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.models is not None and self.counts is not None:
              stat = self.stat_sum()
          str_ += "\t{:32}: {:.2f}\n\n".format("Fit statistic value (-2 log(L))", stat)
  
          n_pars, n_free_pars = 0, 0
          if self.models is not None:
              n_pars = len(self.models.parameters)
              n_free_pars = len(self.parameters.free_parameters)
  
          str_ += "\t{:32}: {}\n".format("Number of parameters", n_pars)
          str_ += "\t{:32}: {}\n\n".format("Number of free parameters", n_free_pars)
  
          if self.models is not None:
              str_ += "\t{:32}: {}\n".format("Model type", self.models.__class__.__name__)
              info = str(self.models.parameters)
              lines = info.split("\n")
              for line in lines[2:-1]:
                  str_ += "\t" + line.replace(":", "\t:") + "\n"
  
          return str_.expandtabs(tabsize=4)
  
      @property
      def models(self):
          """Models (`gammapy.modeling.models.SkyModels`)."""
          return self._models
  
      @models.setter
      def models(self, value):
          if value is not None:
              self._models = SkyModels(value)
          else:
              self._models = None
  
          evaluators = []
          if self.models is not None:
              for model in self.models:
                  evaluator = SpectrumEvaluator(
                      model=model,
                      livetime=self.livetime,
                      aeff=self.aeff,
                      edisp=self.edisp,
                  )
                  evaluators.append(evaluator)
  
          self._evaluators = evaluators
  
      @property
      def mask_safe(self):
          if self._mask_safe is None:
              return np.ones(self.data_shape, bool)
          else:
              return self._mask_safe
  
      @mask_safe.setter
      def mask_safe(self, mask):
          self._mask_safe = mask
  
      @property
      def parameters(self):
          """List of parameters (`~gammapy.modeling.Parameters`)"""
          parameters = []
  
          for component in self.models:
              parameters.append(component.spectral_model.parameters)
  
          return Parameters.from_stack(parameters)
  
      @property
      def _energy_axis(self):
          if self.counts is not None:
              e_axis = self.counts.energy
          elif self.edisp is not None:
              e_axis = self.edisp.data.axis("e_reco")
          elif self.aeff is not None:
              # assume e_reco = e_true
              e_axis = self.aeff.data.axis("energy")
          return e_axis
  
      @
  
  [... the rest is omitted ...]
files: gammapy/spectrum/dataset.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from pathlib import Path
  import numpy as np
  from astropy import units as u
  from astropy.io import fits
  from astropy.table import Table
  from gammapy.data import GTI
  from gammapy.irf import EffectiveAreaTable, EDispKernel, IRFStacker
  from gammapy.modeling import Dataset, Parameters
  from gammapy.modeling.models import SkyModel, SkyModels
  from gammapy.stats import cash, significance_on_off, significance, wstat
  from gammapy.utils.fits import energy_axis_to_ebounds
  from gammapy.utils.random import get_random_state
  from gammapy.utils.scripts import make_path
  from .core import CountsSpectrum, SpectrumEvaluator
  
  __all__ = [
      "SpectrumDatasetOnOff",
      "SpectrumDataset",
      "plot_spectrum_datasets_off_regions",
  ]
  
  
  class SpectrumDataset(Dataset):
      """Spectrum dataset for likelihood fitting.
  
      The spectrum dataset bundles reduced counts data, with a spectral model,
      background model and instrument response function to compute the fit-statistic
      given the current model and data.
  
      Parameters
      ----------
      models : `~gammapy.modeling.models.SkyModels`
          Fit model
      counts : `~gammapy.spectrum.CountsSpectrum`
          Counts spectrum
      livetime : `~astropy.units.Quantity`
          Livetime
      aeff : `~gammapy.irf.EffectiveAreaTable`
          Effective area
      edisp : `~gammapy.irf.EDispKernel`
          Energy dispersion
      background : `~gammapy.spectrum.CountsSpectrum`
          Background to use for the fit.
      mask_safe : `~numpy.ndarray`
          Mask defining the safe data range.
      mask_fit : `~numpy.ndarray`
          Mask to apply to the likelihood for fitting.
      name : str
          Dataset name.
      gti : `~gammapy.data.GTI`
          GTI of the observation or union of GTI if it is a stacked observation
  
      See Also
      --------
      SpectrumDatasetOnOff, FluxPointsDataset, gammapy.cube.MapDataset
      """
  
      stat_type = "cash"
      tag = "SpectrumDataset"
  
      def __init__(
          self,
          models=None,
          counts=None,
          livetime=None,
          aeff=None,
          edisp=None,
          background=None,
          mask_safe=None,
          mask_fit=None,
          name="",
          gti=None,
      ):
  
          if mask_fit is not None and mask_fit.dtype != np.dtype("bool"):
              raise ValueError("mask data must have dtype bool")
  
          self.counts = counts
  
          if livetime is not None:
              livetime = u.Quantity(livetime)
  
          self.livetime = livetime
          self.mask_fit = mask_fit
          self.aeff = aeff
          self.edisp = edisp
          self.background = background
          self.models = models
          self.mask_safe = mask_safe
          self.name = name
          self.gti = gti
  
      def __str__(self):
          str_ = self.__class__.__name__
          str_ += "\n\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)
  
          counts_off = np.nan
          if getattr(self, "counts_off", None) is not None:
              counts_off = np.sum(self.counts_off.data)
              str_ += "\t{:32}: {:.2f}\n\n".format("Total off counts", counts_off)
  
          background = np.nan
          if getattr(self, "background", None) is not None:
              background = np.sum(self.background.data)
              str_ += "\t{:32}: {:.2f}\n\n".format("Total background counts", background)
  
          aeff_min, aeff_max, aeff_unit = np.nan, np.nan, ""
          if self.aeff is not None:
              aeff_min = np.min(self.aeff.data.data.value[self.aeff.data.data.value > 0])
              aeff_max = np.max(self.aeff.data.data.value)
              aeff_unit = self.aeff.data.data.unit
  
          str_ += "\t{:32}: {:.2e} {}\n".format("Effective area min", aeff_min, aeff_unit)
          str_ += "\t{:32}: {:.2e} {}\n\n".format(
              "Effective area max", aeff_max, aeff_unit
          )
  
          livetime = np.nan
          if self.livetime is not None:
              livetime = self.livetime
          str_ += "\t{:32}: {:.2e}\n\n".format("Livetime", livetime)
  
          # 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)
          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.models is not None and self.counts is not None:
              stat = self.stat_sum()
          str_ += "\t{:32}: {:.2f}\n\n".format("Fit statistic value (-2 log(L))", stat)
  
          n_pars, n_free_pars = 0, 0
          if self.models is not None:
              n_pars = len(self.models.parameters)
              n_free_pars = len(self.parameters.free_parameters)
  
          str_ += "\t{:32}: {}\n".format("Number of parameters", n_pars)
          str_ += "\t{:32}: {}\n\n".format("Number of free parameters", n_free_pars)
  
          if self.models is not None:
              str_ += "\t{:32}: {}\n".format("Model type", self.models.__class__.__name__)
              info = str(self.models.parameters)
              lines = info.split("\n")
              for line in lines[2:-1]:
                  str_ += "\t" + line.replace(":", "\t:") + "\n"
  
          return str_.expandtabs(tabsize=4)
  
      @property
      def models(self):
          """Models (`gammapy.modeling.models.SkyModels`)."""
          return self._models
  
      @models.setter
      def models(self, value):
          if value is not None:
              self._models = SkyModels(value)
          else:
              self._models = None
  
          evaluators = []
          if self.models is not None:
              for model in self.models:
                  evaluator = SpectrumEvaluator(
                      model=model,
                      livetime=self.livetime,
                      aeff=self.aeff,
                      edisp=self.edisp,
                  )
                  evaluators.append(evaluator)
  
          self._evaluators = evaluators
  
      @property
      def mask_safe(self):
          if self._mask_safe is None:
              return np.ones(self.data_shape, bool)
          else:
              return self._mask_safe
  
      @mask_safe.setter
      def mask_safe(self, mask):
          self._mask_safe = mask
  
      @property
      def parameters(self):
          """List of parameters (`~gammapy.modeling.Parameters`)"""
          parameters = []
  
          for component in self.models:
              parameters.append(component.spectral_model.parameters)
  
          return Parameters.from_stack(parameters)
  
      @property
      def _energy_axis(self):
          if self.counts is not None:
              e_axis = self.counts.energy
          elif self.edisp is not None:
              e_axis = self.edisp.data.axis("e_reco")
          elif self.aeff is not None:
              # assume e_reco = e_true
              e_axis = self.aeff.data.axis("energy")
          return e_axis
  
      @property
      def data_shape(self):
          """Shape of the counts data"""
          return (
  
  [... the rest is omitted ...]
files: gammapy/spectrum/dataset.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import logging
  import numpy as np
  from astropy import units as u
  from astropy.io.registry import IORegistryError
  from astropy.table import Table, vstack
  from gammapy.modeling import Dataset, Datasets, Fit, Parameters
  from gammapy.modeling.models import (
      PowerLawSpectralModel,
      ScaleSpectralModel,
      SkyModel,
      SkyModels,
  )
  from gammapy.utils.interpolation import interpolate_profile
  from gammapy.utils.scripts import make_name, make_path
  from gammapy.utils.table import table_from_row_data, table_standardise_units_copy
  from .dataset import SpectrumDatasetOnOff
  
  __all__ = ["FluxPoints", "FluxPointsEstimator", "FluxPointsDataset"]
  
  log = logging.getLogger(__name__)
  
  REQUIRED_COLUMNS = {
      "dnde": ["e_ref", "dnde"],
      "e2dnde": ["e_ref", "e2dnde"],
      "flux": ["e_min", "e_max", "flux"],
      "eflux": ["e_min", "e_max", "eflux"],
      # TODO: extend required columns
      "likelihood": [
          "e_min",
          "e_max",
          "e_ref",
          "ref_dnde",
          "norm",
          "norm_scan",
          "stat_scan",
      ],
  }
  
  OPTIONAL_COLUMNS = {
      "dnde": ["dnde_err", "dnde_errp", "dnde_errn", "dnde_ul", "is_ul"],
      "e2dnde": ["e2dnde_err", "e2dnde_errp", "e2dnde_errn", "e2dnde_ul", "is_ul"],
      "flux": ["flux_err", "flux_errp", "flux_errn", "flux_ul", "is_ul"],
      "eflux": ["eflux_err", "eflux_errp", "eflux_errn", "eflux_ul", "is_ul"],
  }
  
  DEFAULT_UNIT = {
      "dnde": u.Unit("cm-2 s-1 TeV-1"),
      "e2dnde": u.Unit("erg cm-2 s-1"),
      "flux": u.Unit("cm-2 s-1"),
      "eflux": u.Unit("erg cm-2 s-1"),
  }
  
  
  class FluxPoints:
      """Flux points container.
  
      The supported formats are described here: :ref:`gadf:flux-points`
  
      In summary, the following formats and minimum required columns are:
  
      * Format ``dnde``: columns ``e_ref`` and ``dnde``
      * Format ``e2dnde``: columns ``e_ref``, ``e2dnde``
      * Format ``flux``: columns ``e_min``, ``e_max``, ``flux``
      * Format ``eflux``: columns ``e_min``, ``e_max``, ``eflux``
  
      Parameters
      ----------
      table : `~astropy.table.Table`
          Table with flux point data
  
      Attributes
      ----------
      table : `~astropy.table.Table`
          Table with flux point data
  
      Examples
      --------
      The `FluxPoints` object is most easily created by reading a file with
      flux points given in one of the formats documented above::
  
          from gammapy.spectrum import FluxPoints
          filename = '$GAMMAPY_DATA/tests/spectrum/flux_points/flux_points.fits'
          flux_points = FluxPoints.read(filename)
          flux_points.plot()
  
      An instance of `FluxPoints` can also be created by passing an instance of
      `astropy.table.Table`, which contains the required columns, such as `'e_ref'`
      and `'dnde'`. The corresponding `sed_type` has to be defined in the meta data
      of the table::
  
          from astropy import units as u
          from astropy.table import Table
          from gammapy.spectrum import FluxPoints
          from gammapy.modeling.models import PowerLawSpectralModel
  
          table = Table()
          pwl = PowerLawSpectralModel()
          e_ref = np.logspace(0, 2, 7) * u.TeV
          table['e_ref'] = e_ref
          table['dnde'] = pwl(e_ref)
          table.meta['SED_TYPE'] = 'dnde'
  
          flux_points = FluxPoints(table)
          flux_points.plot()
  
      If you have flux points in a different data format, the format can be changed
      by renaming the table columns and adding meta data::
  
  
          from astropy import units as u
          from astropy.table import Table
          from gammapy.spectrum import FluxPoints
  
          table = Table.read('$GAMMAPY_DATA/tests/spectrum/flux_points/flux_points_ctb_37b.txt',
                             format='ascii.csv', delimiter=' ', comment='#')
          table.meta['SED_TYPE'] = 'dnde'
          table.rename_column('Differential_Flux', 'dnde')
          table['dnde'].unit = 'cm-2 s-1 TeV-1'
  
          table.rename_column('lower_error', 'dnde_errn')
          table['dnde_errn'].unit = 'cm-2 s-1 TeV-1'
  
          table.rename_column('upper_error', 'dnde_errp')
          table['dnde_errp'].unit = 'cm-2 s-1 TeV-1'
  
          table.rename_column('E', 'e_ref')
          table['e_ref'].unit = 'TeV'
  
          flux_points = FluxPoints(table)
          flux_points.plot()
  
      """
  
      def __init__(self, table):
          self.table = table_standardise_units_copy(table)
          # validate that the table is a valid representation
          # of the given flux point sed type
          self._validate_table(self.table, table.meta["SED_TYPE"])
  
      def __repr__(self):
          return f"{self.__class__.__name__}(sed_type={self.sed_type!r}, n_points={len(self.table)})"
  
      @property
      def table_formatted(self):
          """Return formatted version of the flux points table. Used for pretty printing"""
          table = self.table.copy()
  
          for column in table.colnames:
              if column.startswith(("dnde", "eflux", "flux", "e2dnde", "ref")):
                  table[column].format = ".3e"
              elif column.startswith(
                  ("e_min", "e_max", "e_ref", "sqrt_ts", "norm", "ts", "stat")
              ):
                  table[column].format = ".3f"
  
          return table
  
      @classmethod
      def read(cls, filename, **kwargs):
          """Read flux points.
  
          Parameters
          ----------
          filename : str
              Filename
          kwargs : dict
              Keyword arguments passed to `astropy.table.Table.read`.
          """
          filename = make_path(filename)
          try:
              table = Table.read(filename, **kwargs)
          except IORegistryError:
              kwargs.setdefault("format", "ascii.ecsv")
              table = Table.read(filename, **kwargs)
  
          if "SED_TYPE" not in table.meta.keys():
              sed_type = cls._guess_sed_type(table)
              table.meta["SED_TYPE"] = sed_type
  
          # TODO: check sign and factor 2 here
          # https://github.com/gammapy/gammapy/pull/2546#issuecomment-554274318
          # The idea below is to support the format here:
          # https://gamma-astro-data-formats.readthedocs.io/en/latest/spectra/flux_points/index.html#likelihood-columns
          # but internally to go to the uniform "stat"
  
          if "loglike" in table.colnames and "stat" not in table.colnames:
              table["stat"] = 2 * table["loglike"]
  
          if "loglike_null" in table.colnames and "stat_null" not in table.colnames:
              table["stat_null"] = 2 * table["loglike_null"]
  
          if "dloglike_scan" in table.colnames and "stat_scan" not in table.colnames:
              table["stat_scan"] = 2 * table["dloglike_scan"]
  
          return cls(table=table)
  
      def write(self, filename, **kwargs):
          """Write flux points.
  
          Parameters
          ----------
          filename : str
              Filename
          kwargs : dict
              Keyword arguments passed to `astropy.table.Table.write`.
          """
          filename = make_path(filename)
          try:
              self.table.write(filename, **kwargs)
          except IORegistryError:
              kwargs.setdefault("format", "ascii.ecsv")
              self.table.write(filename, **kwargs)
  
      @classmethod
      def stack(cls, flux_points):
          """Create flux points by stacking list of flux points.
  
          The first `FluxPoints`
  
  [... the rest is omitted ...]
files: gammapy/spectrum/flux_point.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import logging
  import numpy as np
  from astropy import units as u
  from astropy.io.registry import IORegistryError
  from astropy.table import Table, vstack
  from gammapy.modeling import Dataset, Datasets, Fit, Parameters
  from gammapy.modeling.models import (
      PowerLawSpectralModel,
      ScaleSpectralModel,
      SkyModel,
      SkyModels,
  )
  from gammapy.utils.interpolation import interpolate_profile
  from gammapy.utils.scripts import make_path
  from gammapy.utils.table import table_from_row_data, table_standardise_units_copy
  from .dataset import SpectrumDatasetOnOff
  
  __all__ = ["FluxPoints", "FluxPointsEstimator", "FluxPointsDataset"]
  
  log = logging.getLogger(__name__)
  
  REQUIRED_COLUMNS = {
      "dnde": ["e_ref", "dnde"],
      "e2dnde": ["e_ref", "e2dnde"],
      "flux": ["e_min", "e_max", "flux"],
      "eflux": ["e_min", "e_max", "eflux"],
      # TODO: extend required columns
      "likelihood": [
          "e_min",
          "e_max",
          "e_ref",
          "ref_dnde",
          "norm",
          "norm_scan",
          "stat_scan",
      ],
  }
  
  OPTIONAL_COLUMNS = {
      "dnde": ["dnde_err", "dnde_errp", "dnde_errn", "dnde_ul", "is_ul"],
      "e2dnde": ["e2dnde_err", "e2dnde_errp", "e2dnde_errn", "e2dnde_ul", "is_ul"],
      "flux": ["flux_err", "flux_errp", "flux_errn", "flux_ul", "is_ul"],
      "eflux": ["eflux_err", "eflux_errp", "eflux_errn", "eflux_ul", "is_ul"],
  }
  
  DEFAULT_UNIT = {
      "dnde": u.Unit("cm-2 s-1 TeV-1"),
      "e2dnde": u.Unit("erg cm-2 s-1"),
      "flux": u.Unit("cm-2 s-1"),
      "eflux": u.Unit("erg cm-2 s-1"),
  }
  
  
  class FluxPoints:
      """Flux points container.
  
      The supported formats are described here: :ref:`gadf:flux-points`
  
      In summary, the following formats and minimum required columns are:
  
      * Format ``dnde``: columns ``e_ref`` and ``dnde``
      * Format ``e2dnde``: columns ``e_ref``, ``e2dnde``
      * Format ``flux``: columns ``e_min``, ``e_max``, ``flux``
      * Format ``eflux``: columns ``e_min``, ``e_max``, ``eflux``
  
      Parameters
      ----------
      table : `~astropy.table.Table`
          Table with flux point data
  
      Attributes
      ----------
      table : `~astropy.table.Table`
          Table with flux point data
  
      Examples
      --------
      The `FluxPoints` object is most easily created by reading a file with
      flux points given in one of the formats documented above::
  
          from gammapy.spectrum import FluxPoints
          filename = '$GAMMAPY_DATA/tests/spectrum/flux_points/flux_points.fits'
          flux_points = FluxPoints.read(filename)
          flux_points.plot()
  
      An instance of `FluxPoints` can also be created by passing an instance of
      `astropy.table.Table`, which contains the required columns, such as `'e_ref'`
      and `'dnde'`. The corresponding `sed_type` has to be defined in the meta data
      of the table::
  
          from astropy import units as u
          from astropy.table import Table
          from gammapy.spectrum import FluxPoints
          from gammapy.modeling.models import PowerLawSpectralModel
  
          table = Table()
          pwl = PowerLawSpectralModel()
          e_ref = np.logspace(0, 2, 7) * u.TeV
          table['e_ref'] = e_ref
          table['dnde'] = pwl(e_ref)
          table.meta['SED_TYPE'] = 'dnde'
  
          flux_points = FluxPoints(table)
          flux_points.plot()
  
      If you have flux points in a different data format, the format can be changed
      by renaming the table columns and adding meta data::
  
  
          from astropy import units as u
          from astropy.table import Table
          from gammapy.spectrum import FluxPoints
  
          table = Table.read('$GAMMAPY_DATA/tests/spectrum/flux_points/flux_points_ctb_37b.txt',
                             format='ascii.csv', delimiter=' ', comment='#')
          table.meta['SED_TYPE'] = 'dnde'
          table.rename_column('Differential_Flux', 'dnde')
          table['dnde'].unit = 'cm-2 s-1 TeV-1'
  
          table.rename_column('lower_error', 'dnde_errn')
          table['dnde_errn'].unit = 'cm-2 s-1 TeV-1'
  
          table.rename_column('upper_error', 'dnde_errp')
          table['dnde_errp'].unit = 'cm-2 s-1 TeV-1'
  
          table.rename_column('E', 'e_ref')
          table['e_ref'].unit = 'TeV'
  
          flux_points = FluxPoints(table)
          flux_points.plot()
  
      """
  
      def __init__(self, table):
          self.table = table_standardise_units_copy(table)
          # validate that the table is a valid representation
          # of the given flux point sed type
          self._validate_table(self.table, table.meta["SED_TYPE"])
  
      def __repr__(self):
          return f"{self.__class__.__name__}(sed_type={self.sed_type!r}, n_points={len(self.table)})"
  
      @property
      def table_formatted(self):
          """Return formatted version of the flux points table. Used for pretty printing"""
          table = self.table.copy()
  
          for column in table.colnames:
              if column.startswith(("dnde", "eflux", "flux", "e2dnde", "ref")):
                  table[column].format = ".3e"
              elif column.startswith(
                  ("e_min", "e_max", "e_ref", "sqrt_ts", "norm", "ts", "stat")
              ):
                  table[column].format = ".3f"
  
          return table
  
      @classmethod
      def read(cls, filename, **kwargs):
          """Read flux points.
  
          Parameters
          ----------
          filename : str
              Filename
          kwargs : dict
              Keyword arguments passed to `astropy.table.Table.read`.
          """
          filename = make_path(filename)
          try:
              table = Table.read(filename, **kwargs)
          except IORegistryError:
              kwargs.setdefault("format", "ascii.ecsv")
              table = Table.read(filename, **kwargs)
  
          if "SED_TYPE" not in table.meta.keys():
              sed_type = cls._guess_sed_type(table)
              table.meta["SED_TYPE"] = sed_type
  
          # TODO: check sign and factor 2 here
          # https://github.com/gammapy/gammapy/pull/2546#issuecomment-554274318
          # The idea below is to support the format here:
          # https://gamma-astro-data-formats.readthedocs.io/en/latest/spectra/flux_points/index.html#likelihood-columns
          # but internally to go to the uniform "stat"
  
          if "loglike" in table.colnames and "stat" not in table.colnames:
              table["stat"] = 2 * table["loglike"]
  
          if "loglike_null" in table.colnames and "stat_null" not in table.colnames:
              table["stat_null"] = 2 * table["loglike_null"]
  
          if "dloglike_scan" in table.colnames and "stat_scan" not in table.colnames:
              table["stat_scan"] = 2 * table["dloglike_scan"]
  
          return cls(table=table)
  
      def write(self, filename, **kwargs):
          """Write flux points.
  
          Parameters
          ----------
          filename : str
              Filename
          kwargs : dict
              Keyword arguments passed to `astropy.table.Table.write`.
          """
          filename = make_path(filename)
          try:
              self.table.write(filename, **kwargs)
          except IORegistryError:
              kwargs.setdefault("format", "ascii.ecsv")
              self.table.write(filename, **kwargs)
  
      @classmethod
      def stack(cls, flux_points):
          """Create flux points by stacking list of flux points.
  
          The first `FluxPoints` object in the
  
  [... the rest is omitted ...]
files: gammapy/spectrum/flux_point.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Utils to create scripts and command-line tools"""
  import os.path
  from pathlib import Path
  from uuid import uuid4
  import yaml
  
  __all__ = ["read_yaml", "write_yaml", "make_path", "recursive_merge_dicts"]
  
  
  def read_yaml(filename, logger=None):
      """Read YAML file.
  
      Parameters
      ----------
      filename : `~pathlib.Path`
          Filename
      logger : `~logging.Logger`
          Logger
  
      Returns
      -------
      data : dict
          YAML file content as a dict
      """
      path = make_path(filename)
      if logger is not None:
          logger.info(f"Reading {path}")
  
      text = path.read_text()
      return yaml.safe_load(text)
  
  
  def write_yaml(dictionary, filename, logger=None, sort_keys=True):
      """Write YAML file.
  
      Parameters
      ----------
      dictionary : dict
          Python dictionary
      filename : `~pathlib.Path`
          Filename
      logger : `~logging.Logger`
          Logger
      sort_keys : bool
          Whether to sort keys.
      """
      text = yaml.safe_dump(dictionary, default_flow_style=False, sort_keys=sort_keys)
  
      path = make_path(filename)
      path.parent.mkdir(exist_ok=True)
      if logger is not None:
          logger.info(f"Writing {path}")
      path.write_text(text)
  
  
  def make_name():
      return uuid4().hex[:8]
  
  
  def make_path(path):
      """Expand environment variables on `~pathlib.Path` construction.
  
      Parameters
      ----------
      path : str, `pathlib.Path`
          path to expand
      """
      # TODO: raise error or warning if environment variables that don't resolve are used
      # e.g. "spam/$DAMN/ham" where `$DAMN` is not defined
      # Otherwise this can result in cryptic errors later on
      return Path(os.path.expandvars(path))
  
  
  def recursive_merge_dicts(a, b):
      """Recursively merge two dictionaries.
  
      Entries in b override entries in a. The built-in update function cannot be
      used for hierarchical dicts, see:
      http://stackoverflow.com/questions/3232943/update-value-of-a-nested-dictionary-of-varying-depth/3233356#3233356
  
      Parameters
      ----------
      a : dict
          dictionary to be merged
      b : dict
          dictionary to be merged
  
      Returns
      -------
      c : dict
          merged dict
  
      Examples
      --------
      >>> from gammapy.utils.scripts import recursive_merge_dicts
      >>> a = dict(a=42, b=dict(c=43, e=44))
      >>> b = dict(d=99, b=dict(c=50, g=98))
      >>> c = recursive_merge_dicts(a, b)
      >>> print(c)
      {'a': 42, 'b': {'c': 50, 'e': 44, 'g': 98}, 'd': 99}
      """
      c = a.copy()
      for k, v in b.items():
          if k in c and isinstance(c[k], dict):
              c[k] = recursive_merge_dicts(c[k], v)
          else:
              c[k] = v
      return c
  
  
  [... the rest is omitted ...]
files: gammapy/utils/scripts.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Utils to create scripts and command-line tools"""
  import os.path
  from pathlib import Path
  import yaml
  
  __all__ = ["read_yaml", "write_yaml", "make_path", "recursive_merge_dicts"]
  
  
  def read_yaml(filename, logger=None):
      """Read YAML file.
  
      Parameters
      ----------
      filename : `~pathlib.Path`
          Filename
      logger : `~logging.Logger`
          Logger
  
      Returns
      -------
      data : dict
          YAML file content as a dict
      """
      path = make_path(filename)
      if logger is not None:
          logger.info(f"Reading {path}")
  
      text = path.read_text()
      return yaml.safe_load(text)
  
  
  def write_yaml(dictionary, filename, logger=None, sort_keys=True):
      """Write YAML file.
  
      Parameters
      ----------
      dictionary : dict
          Python dictionary
      filename : `~pathlib.Path`
          Filename
      logger : `~logging.Logger`
          Logger
      sort_keys : bool
          Whether to sort keys.
      """
      text = yaml.safe_dump(dictionary, default_flow_style=False, sort_keys=sort_keys)
  
      path = make_path(filename)
      path.parent.mkdir(exist_ok=True)
      if logger is not None:
          logger.info(f"Writing {path}")
      path.write_text(text)
  
  
  def make_path(path):
      """Expand environment variables on `~pathlib.Path` construction.
  
      Parameters
      ----------
      path : str, `pathlib.Path`
          path to expand
      """
      # TODO: raise error or warning if environment variables that don't resolve are used
      # e.g. "spam/$DAMN/ham" where `$DAMN` is not defined
      # Otherwise this can result in cryptic errors later on
      return Path(os.path.expandvars(path))
  
  
  def recursive_merge_dicts(a, b):
      """Recursively merge two dictionaries.
  
      Entries in b override entries in a. The built-in update function cannot be
      used for hierarchical dicts, see:
      http://stackoverflow.com/questions/3232943/update-value-of-a-nested-dictionary-of-varying-depth/3233356#3233356
  
      Parameters
      ----------
      a : dict
          dictionary to be merged
      b : dict
          dictionary to be merged
  
      Returns
      -------
      c : dict
          merged dict
  
      Examples
      --------
      >>> from gammapy.utils.scripts import recursive_merge_dicts
      >>> a = dict(a=42, b=dict(c=43, e=44))
      >>> b = dict(d=99, b=dict(c=50, g=98))
      >>> c = recursive_merge_dicts(a, b)
      >>> print(c)
      {'a': 42, 'b': {'c': 50, 'e': 44, 'g': 98}, 'd': 99}
      """
      c = a.copy()
      for k, v in b.items():
          if k in c and isinstance(c[k], dict):
              c[k] = recursive_merge_dicts(c[k], v)
          else:
              c[k] = v
      return c
  
  
  [... the rest is omitted ...]
files: gammapy/utils/scripts.py