files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import collections.abc
  import copy
  from os.path import split
  import numpy as np
  import astropy.units as u
  from astropy.table import Table
  import yaml
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.scripts import make_name, make_path
  
  
  def _set_link(shared_register, model):
      for param in model.parameters:
          name = param.name
          link_label = param._link_label_io
          if link_label is not None:
              if link_label in shared_register:
                  new_param = shared_register[link_label]
                  setattr(model, name, new_param)
              else:
                  shared_register[link_label] = param
      return shared_register
  
  
  __all__ = ["Model", "Models", "ProperModels"]
  
  
  class Model:
      """Model base class."""
  
      def __init__(self, **kwargs):
          # Copy default parameters from the class to the instance
          default_parameters = self.default_parameters.copy()
  
          for par in default_parameters:
              value = kwargs.get(par.name, par)
  
              if not isinstance(value, Parameter):
                  par.quantity = u.Quantity(value)
              else:
                  par = value
  
              setattr(self, par.name, par)
  
          self._covariance = Covariance(self.parameters)
  
      def __init_subclass__(cls, **kwargs):
          # Add parameters list on the model sub-class (not instances)
          cls.default_parameters = Parameters(
              [_ for _ in cls.__dict__.values() if isinstance(_, Parameter)]
          )
  
      @classmethod
      def from_parameters(cls, parameters, **kwargs):
          """Create model from parameter list
  
          Parameters
          ----------
          parameters : `Parameters`
              Parameters for init
  
          Returns
          -------
          model : `Model`
              Model instance
          """
          for par in parameters:
              kwargs[par.name] = par
          return cls(**kwargs)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance(self.parameters)
  
      @property
      def covariance(self):
          self._check_covariance()
          for par in self.parameters:
              pars = Parameters([par])
              error = np.nan_to_num(par.error ** 2, nan=1)
              covar = Covariance(pars, data=[[error]])
              self._covariance.set_subcovariance(covar)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for par in self.parameters:
              pars = Parameters([par])
              variance = self._covariance.get_subcovariance(pars)
              par.error = np.sqrt(variance)
  
      @property
      def parameters(self):
          """Parameters (`~gammapy.modeling.Parameters`)"""
          return Parameters(
              [getattr(self, name) for name in self.default_parameters.names]
          )
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      def to_dict(self, full_output=True):
          """Create dict for YAML serialisation"""
          tag = self.tag[0] if isinstance(self.tag, list) else self.tag
          params = self.parameters.to_dict()
  
          if full_output is False:
              base = self.__class__
              names = self.parameters.names
              for k, name in enumerate(names):
                  init = base.__dict__[name].to_dict()
                  for item in ["min", "max", "frozen", "error"]:
                      if params[k][item] == init[item] or np.isnan(init[item]):
                          del params[k][item]
          return {"type": tag, "parameters": params}
  
      @classmethod
      def from_dict(cls, data):
          kwargs = {}
          parameters = Parameters.from_dict(data["parameters"])
  
          # TODO: this is a special case for spatial models, maybe better move to `SpatialModel` base class
          if "frame" in data:
              kwargs["frame"] = data["frame"]
  
          return cls.from_parameters(parameters, **kwargs)
  
      @staticmethod
      def create(tag, model_type=None, *args, **kwargs):
          """Create a model instance.
  
          Examples
          --------
          >>> from gammapy.modeling.models import Model
          >>> spectral_model = Model.create("pl-2", model_type="spectral", amplitude="1e-10 cm-2 s-1", index=3)
          >>> type(spectral_model)
          gammapy.modeling.models.spectral.PowerLaw2SpectralModel
          """
          from . import (
              MODEL_REGISTRY,
              SPATIAL_MODEL_REGISTRY,
              SPECTRAL_MODEL_REGISTRY,
              TEMPORAL_MODEL_REGISTRY,
          )
  
          if model_type is None:
              cls = MODEL_REGISTRY.get_cls(tag)
          else:
              registry = {
                  "spatial": SPATIAL_MODEL_REGISTRY,
                  "spectral": SPECTRAL_MODEL_REGISTRY,
                  "temporal": TEMPORAL_MODEL_REGISTRY,
              }
              cls = registry[model_type].get_cls(tag)
          return cls(*args, **kwargs)
  
      def __str__(self):
          string = f"{self.__class__.__name__}\n"
          if len(self.parameters) > 0:
              string += f"\n{self.parameters.to_table()}"
          return string
  
  
  class Models(collections.abc.MutableSequence):
      """Sky model collection.
  
      Parameters
      ----------
      models : `SkyModel`, list of `SkyModel` or `Models`
          Sky models
      """
  
      def __init__(self, models=None):
          if models is None:
              models = []
  
          if isinstance(models, Models):
              models = models._models
          elif isinstance(models, Model):
              models = [models]
          elif not isinstance(models, list):
              raise TypeError(f"Invalid type: {models!r}")
  
          unique_names = []
          for model in models:
              if model.name in unique_names:
                  raise (ValueError("Model names must be unique"))
              unique_names.append(model.name)
  
          self._models = models
          self._covar_file = None
          self._covariance = Covariance(self.parameters)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in self._models]
              )
  
      @property
      def covariance(self):
          self._check_covariance()
  
          for model in self._models:
              self._covariance.set_subcovariance(model.covariance)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for model in self._models:
              subcovar = self._covariance.get_subcovariance(model.covariance.parameters)
              model.covariance = subcovar
  
      @property
      def parameters(self):
          return Parameters.from_stack([_.parameters for _ in self._models])
  
      @property
      def parameters_unique_names(self):
          """List of unique parameter names as model_name.par_type.par_name"""
          names = []
          for model in self:
              for par in model.parameters:
                  components = [model.name, par.type, par.name]
                  name = ".".join(components)
                  names.append(name)
  
          return names
  
      @property
      def names(self):
          return [m.name for m in self._models]
  
      @classmethod
      def read(cls, filename):
          """Read from YAML file."""
          yaml_str = make_path(filename).read_text()
          path, filename = split(filename)
          return cls.from_yaml(yaml_str, path=path)
  
      @classmethod
      def from_yaml(cls, yaml_str, path=""):
          """Create from YAML string."""
          data = yaml.safe_load(yaml_str)
          return cls.from_dict(data, path=path)
  
      @classmethod
      def from_dict(cls, data, path=""):
          """Create from dict."""
          from . import MODEL_REGISTRY, SkyModel
  
          models = []
  
          for component in data["components"]:
              model_cls = MODEL_REGISTRY.get_cls(component["type"])
              model = model_cls.from_dict(component)
              models.append(model)
  
          models = cls(models)
  
          if "covariance" in data:
              filename = data["covariance"]
              path = make_path(path)
              if not (path / filename).exists():
                  path, filename = split(filename)
  
              models.read_covariance(path, filename, format="ascii.fixed_width")
  
          shared_register = {}
          for model in models:
              if isinstance(model, SkyModel):
                  submodels = [
                      model.spectral_model,
                      model.spatial_model,
                      model.temporal_model,
                  ]
                  for submodel in submodels:
                      if submodel is not None:
                          shared_register = _set_link(shared_register, submodel)
              else:
                  shared_register = _set_link(shared_register, model)
          return models
  
      def write(self, path, overwrite=False, full_output=True, write_covariance=True):
          """Write to YAML file.
  
          Parameters
          ----------
          path : `pathlib.Path`
              path to write files
          overwrite : bool
              overwrite files
          write_covariance : bool
              save covariance or not
          """
  
          base_path, _ = split(path)
          path = make_path(path)
          base_path = make_path(base_path)
  
          if path.exists() and not overwrite:
              raise IOError(f"File exists already: {path}")
  
          if (
              write_covariance
              and self.covariance is not None
              and len(self.parameters) != 0
          ):
              filecovar = path.stem + "_covariance.dat"
              kwargs = dict(
                  format="ascii.fixed_width", delimiter="|", overwrite=overwrite
              )
              self.write_covariance(base_path / filecovar, **kwargs)
              self._covar_file = filecovar
  
          path.write_text(self.to_yaml())
  
      def to_yaml(self, full_output=True):
          """Convert to YAML string."""
          data = self.to_dict(full_output)
          return yaml.dump(
              data, sort_keys=False, indent=4, width=80, default_flow_style=False
          )
  
      def to_dict(self, full_output=True):
          """Convert to dict."""
          # update linked parameters labels
          params_list = []
          params_shared = []
          for param in self.parameters:
              if param not in params_list:
                  params_list.append(param)
                  params_list.append(param)
              elif param not in params_shared:
                  params_shared.append(param)
          for param in params_shared:
              param._link_label_io = param.name + "@" + make_name()
  
          models_data = []
          for model in self._models:
              model_data = model.to_dict(full_output)
              models_data.append(model_data)
          if self._covar_file is not None:
              return {
                  "components": models_data,
                  "covariance": str(self._covar_file),
              }
          else:
              return {"components": models_data}
  
      def read_covariance(self, path, filename="_covariance.dat", **kwargs):
          """Read covariance data from file
  
          Parameters
          ----------
          filename : str
              Filename
          **kwargs : dict
              Keyword arguments passed to `~astropy.table.Table.read`
  
          """
          path = make_path(path)
          filepath = str(path / filename)
          t = Table.read(filepath, **kwargs)
          t.remove_column("Parameters")
          arr = np.array(t)
          data = arr.view(float).reshape(arr.shape + (-1,))
          self.covariance = data
          self._covar_file = filename
  
      def write_covariance(self, filename, **kwargs):
          """Write covariance to file
  
          Parameters
          ----------
          filename : str
              Filename
          **kwargs : dict
              Keyword arguments passed to `~astropy.table.Table.write`
  
          """
          names = self.parameters_unique_names
          table = Table()
          table["Parameters"] = names
  
          for idx, name in enumerate(names):
              values = self.covariance.data[idx]
              table[name] = values
  
          table.write(make_path(filename), **kwargs)
  
      def __str__(self):
          str_ = f"{self.__class__.__name__}\n\n"
  
          for idx, model in enumerate(self):
              str_ += f"Component {idx}: "
              str_ += str(model)
  
          return str_.expandtabs(tabsize=2)
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([*self, *other])
          elif isinstance(other, Model):
              if other.name in self.names:
                  raise (ValueError("Model names must be unique"))
              return Models([*self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __getitem__(self, key):
          return self._models[self.index(key)]
  
      def __delitem__(self, key):
          del self._models[self.index(key)]
  
      def __setitem__(self, key, model):
          from gammapy.modeling.models import SkyModel, BackgroundModel
  
          if isinstance(model, (SkyModel, BackgroundModel)):
              self._models[self.index(key)] = model
          else:
              raise TypeError(f"Invalid type: {model!r}")
  
      def insert(self, idx, model):
          if model.name in self.names:
              raise (ValueError("Model names must be unique"))
  
          self._models.insert(idx, model)
  
      def index(self, key):
          if isinstance(key, (int, slice)):
              return key
          elif isinstance(key, str):
              return self.names.index(key)
          elif isinstance(key, Model):
              return self._models.index(key)
          else:
              raise TypeError(f"Invalid type: {type(key)!r}")
  
      def __len__(self):
          return len(self._models)
  
      def _ipython_key_completions_(self):
          return self.names
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
  
  class ProperModels(Models):
      """ Proper Models of a Dataset or Datasets."""
  
      def __init__(
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/core.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import collections.abc
  import copy
  from os.path import split
  import numpy as np
  import astropy.units as u
  from astropy.table import Table
  import yaml
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.scripts import make_name, make_path
  
  
  def _set_link(shared_register, model):
      for param in model.parameters:
          name = param.name
          link_label = param._link_label_io
          if link_label is not None:
              if link_label in shared_register:
                  new_param = shared_register[link_label]
                  setattr(model, name, new_param)
              else:
                  shared_register[link_label] = param
      return shared_register
  
  
  __all__ = ["Model", "Models", "ProperModels"]
  
  
  class Model:
      """Model base class."""
  
      def __init__(self, **kwargs):
          # Copy default parameters from the class to the instance
          default_parameters = self.default_parameters.copy()
  
          for par in default_parameters:
              value = kwargs.get(par.name, par)
  
              if not isinstance(value, Parameter):
                  par.quantity = u.Quantity(value)
              else:
                  par = value
  
              setattr(self, par.name, par)
  
          self._covariance = Covariance(self.parameters)
  
      def __init_subclass__(cls, **kwargs):
          # Add parameters list on the model sub-class (not instances)
          cls.default_parameters = Parameters(
              [_ for _ in cls.__dict__.values() if isinstance(_, Parameter)]
          )
  
      @classmethod
      def from_parameters(cls, parameters, **kwargs):
          """Create model from parameter list
  
          Parameters
          ----------
          parameters : `Parameters`
              Parameters for init
  
          Returns
          -------
          model : `Model`
              Model instance
          """
          for par in parameters:
              kwargs[par.name] = par
          return cls(**kwargs)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance(self.parameters)
  
      @property
      def covariance(self):
          self._check_covariance()
          for par in self.parameters:
              pars = Parameters([par])
              error = np.nan_to_num(par.error ** 2, nan=1)
              covar = Covariance(pars, data=[[error]])
              self._covariance.set_subcovariance(covar)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for par in self.parameters:
              pars = Parameters([par])
              variance = self._covariance.get_subcovariance(pars)
              par.error = np.sqrt(variance)
  
      @property
      def parameters(self):
          """Parameters (`~gammapy.modeling.Parameters`)"""
          return Parameters(
              [getattr(self, name) for name in self.default_parameters.names]
          )
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
      def to_dict(self):
          """Create dict for YAML serialisation"""
          tag = self.tag[0] if isinstance(self.tag, list) else self.tag
          return {"type": tag, "parameters": self.parameters.to_dict()}
  
      @classmethod
      def from_dict(cls, data):
          kwargs = {}
          parameters = Parameters.from_dict(data["parameters"])
  
          # TODO: this is a special case for spatial models, maybe better move to `SpatialModel` base class
          if "frame" in data:
              kwargs["frame"] = data["frame"]
  
          return cls.from_parameters(parameters, **kwargs)
  
      @staticmethod
      def create(tag, model_type=None, *args, **kwargs):
          """Create a model instance.
  
          Examples
          --------
          >>> from gammapy.modeling.models import Model
          >>> spectral_model = Model.create("pl-2", model_type="spectral", amplitude="1e-10 cm-2 s-1", index=3)
          >>> type(spectral_model)
          gammapy.modeling.models.spectral.PowerLaw2SpectralModel
          """
          from . import (
              MODEL_REGISTRY,
              SPATIAL_MODEL_REGISTRY,
              SPECTRAL_MODEL_REGISTRY,
              TEMPORAL_MODEL_REGISTRY,
          )
  
          if model_type is None:
              cls = MODEL_REGISTRY.get_cls(tag)
          else:
              registry = {
                  "spatial": SPATIAL_MODEL_REGISTRY,
                  "spectral": SPECTRAL_MODEL_REGISTRY,
                  "temporal": TEMPORAL_MODEL_REGISTRY,
              }
              cls = registry[model_type].get_cls(tag)
          return cls(*args, **kwargs)
  
      def __str__(self):
          string = f"{self.__class__.__name__}\n"
          if len(self.parameters) > 0:
              string += f"\n{self.parameters.to_table()}"
          return string
  
  
  class Models(collections.abc.MutableSequence):
      """Sky model collection.
  
      Parameters
      ----------
      models : `SkyModel`, list of `SkyModel` or `Models`
          Sky models
      """
  
      def __init__(self, models=None):
          if models is None:
              models = []
  
          if isinstance(models, Models):
              models = models._models
          elif isinstance(models, Model):
              models = [models]
          elif not isinstance(models, list):
              raise TypeError(f"Invalid type: {models!r}")
  
          unique_names = []
          for model in models:
              if model.name in unique_names:
                  raise (ValueError("Model names must be unique"))
              unique_names.append(model.name)
  
          self._models = models
          self._covar_file = None
          self._covariance = Covariance(self.parameters)
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in self._models]
              )
  
      @property
      def covariance(self):
          self._check_covariance()
  
          for model in self._models:
              self._covariance.set_subcovariance(model.covariance)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for model in self._models:
              subcovar = self._covariance.get_subcovariance(model.covariance.parameters)
              model.covariance = subcovar
  
      @property
      def parameters(self):
          return Parameters.from_stack([_.parameters for _ in self._models])
  
      @property
      def parameters_unique_names(self):
          """List of unique parameter names as model_name.par_type.par_name"""
          names = []
          for model in self:
              for par in model.parameters:
                  components = [model.name, par.type, par.name]
                  name = ".".join(components)
                  names.append(name)
  
          return names
  
      @property
      def names(self):
          return [m.name for m in self._models]
  
      @classmethod
      def read(cls, filename):
          """Read from YAML file."""
          yaml_str = make_path(filename).read_text()
          path, filename = split(filename)
          return cls.from_yaml(yaml_str, path=path)
  
      @classmethod
      def from_yaml(cls, yaml_str, path=""):
          """Create from YAML string."""
          data = yaml.safe_load(yaml_str)
          return cls.from_dict(data, path=path)
  
      @classmethod
      def from_dict(cls, data, path=""):
          """Create from dict."""
          from . import MODEL_REGISTRY, SkyModel
  
          models = []
  
          for component in data["components"]:
              model_cls = MODEL_REGISTRY.get_cls(component["type"])
              model = model_cls.from_dict(component)
              models.append(model)
  
          models = cls(models)
  
          if "covariance" in data:
              filename = data["covariance"]
              path = make_path(path)
              if not (path / filename).exists():
                  path, filename = split(filename)
  
              models.read_covariance(path, filename, format="ascii.fixed_width")
  
          shared_register = {}
          for model in models:
              if isinstance(model, SkyModel):
                  submodels = [
                      model.spectral_model,
                      model.spatial_model,
                      model.temporal_model,
                  ]
                  for submodel in submodels:
                      if submodel is not None:
                          shared_register = _set_link(shared_register, submodel)
              else:
                  shared_register = _set_link(shared_register, model)
          return models
  
      def write(self, path, overwrite=False, write_covariance=True):
          """Write to YAML file.
  
          Parameters
          ----------
          path : `pathlib.Path`
              path to write files
          overwrite : bool
              overwrite files
          write_covariance : bool
              save covariance or not
          """
  
          base_path, _ = split(path)
          path = make_path(path)
          base_path = make_path(base_path)
  
          if path.exists() and not overwrite:
              raise IOError(f"File exists already: {path}")
  
          if (
              write_covariance
              and self.covariance is not None
              and len(self.parameters) != 0
          ):
              filecovar = path.stem + "_covariance.dat"
              kwargs = dict(
                  format="ascii.fixed_width", delimiter="|", overwrite=overwrite
              )
              self.write_covariance(base_path / filecovar, **kwargs)
              self._covar_file = filecovar
  
          path.write_text(self.to_yaml())
  
      def to_yaml(self):
          """Convert to YAML string."""
          data = self.to_dict()
          return yaml.dump(
              data, sort_keys=False, indent=4, width=80, default_flow_style=False
          )
  
      def to_dict(self):
          """Convert to dict."""
          # update linked parameters labels
          params_list = []
          params_shared = []
          for param in self.parameters:
              if param not in params_list:
                  params_list.append(param)
                  params_list.append(param)
              elif param not in params_shared:
                  params_shared.append(param)
          for param in params_shared:
              param._link_label_io = param.name + "@" + make_name()
  
          models_data = []
          for model in self._models:
              model_data = model.to_dict()
              models_data.append(model_data)
          if self._covar_file is not None:
              return {
                  "components": models_data,
                  "covariance": str(self._covar_file),
              }
          else:
              return {"components": models_data}
  
      def read_covariance(self, path, filename="_covariance.dat", **kwargs):
          """Read covariance data from file
  
          Parameters
          ----------
          filename : str
              Filename
          **kwargs : dict
              Keyword arguments passed to `~astropy.table.Table.read`
  
          """
          path = make_path(path)
          filepath = str(path / filename)
          t = Table.read(filepath, **kwargs)
          t.remove_column("Parameters")
          arr = np.array(t)
          data = arr.view(float).reshape(arr.shape + (-1,))
          self.covariance = data
          self._covar_file = filename
  
      def write_covariance(self, filename, **kwargs):
          """Write covariance to file
  
          Parameters
          ----------
          filename : str
              Filename
          **kwargs : dict
              Keyword arguments passed to `~astropy.table.Table.write`
  
          """
          names = self.parameters_unique_names
          table = Table()
          table["Parameters"] = names
  
          for idx, name in enumerate(names):
              values = self.covariance.data[idx]
              table[name] = values
  
          table.write(make_path(filename), **kwargs)
  
      def __str__(self):
          str_ = f"{self.__class__.__name__}\n\n"
  
          for idx, model in enumerate(self):
              str_ += f"Component {idx}: "
              str_ += str(model)
  
          return str_.expandtabs(tabsize=2)
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([*self, *other])
          elif isinstance(other, Model):
              if other.name in self.names:
                  raise (ValueError("Model names must be unique"))
              return Models([*self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __getitem__(self, key):
          return self._models[self.index(key)]
  
      def __delitem__(self, key):
          del self._models[self.index(key)]
  
      def __setitem__(self, key, model):
          from gammapy.modeling.models import SkyModel, BackgroundModel
  
          if isinstance(model, (SkyModel, BackgroundModel)):
              self._models[self.index(key)] = model
          else:
              raise TypeError(f"Invalid type: {model!r}")
  
      def insert(self, idx, model):
          if model.name in self.names:
              raise (ValueError("Model names must be unique"))
  
          self._models.insert(idx, model)
  
      def index(self, key):
          if isinstance(key, (int, slice)):
              return key
          elif isinstance(key, str):
              return self.names.index(key)
          elif isinstance(key, Model):
              return self._models.index(key)
          else:
              raise TypeError(f"Invalid type: {type(key)!r}")
  
      def __len__(self):
          return len(self._models)
  
      def _ipython_key_completions_(self):
          return self.names
  
      def copy(self):
          """A deep copy."""
          return copy.deepcopy(self)
  
  
  class ProperModels(Models):
      """ Proper Models of a Dataset or Datasets."""
  
      def __init__(self, parent):
          from gammapy.datasets import Dataset, Datasets
  
          if isinstance(parent, Dataset):
              self._datasets = [parent]
              self._is_dataset = True
          elif isinstance(parent, Datasets):
              self._datasets = parent._datasets
              self._is_dataset = False
          else:
              raise TypeError(f"Invalid type: {type(parent)!r}")
  
          unique_models = []
          for d in self._datasets:
              if d._models is not None:
                  for model in d._
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/core.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Cube models (axes: lon, lat, energy)."""
  import copy
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.modeling.parameter import _get_parameters_str
  from gammapy.utils.scripts import make_name, make_path
  from gammapy.utils.fits import LazyFitsData, HDULocation
  from .core import Model, Models
  from .spatial import SpatialModel
  from .spectral import SpectralModel, PowerLawNormSpectralModel
  from .temporal import TemporalModel
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([self, *other])
          elif isinstance(other, (SkyModel, BackgroundModel)):
              return Models([self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __call__(self, lon, lat, energy, time=None):
          return self.evaluate(lon, lat, energy, time)
  
      def evaluate_geom(self, geom, gti=None):
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat, coords["energy_true"])
  
  
  class SkyModel(SkyModelBase):
      """Sky model component.
  
      This model represents a factorised sky model.
      It has `~gammapy.modeling.Parameters`
      combining the spatial and spectral parameters.
  
      Parameters
      ----------
      spectral_model : `~gammapy.modeling.models.SpectralModel`
          Spectral model
      spatial_model : `~gammapy.modeling.models.SpatialModel`
          Spatial model (must be normalised to integrate to 1)
      temporal_model : `~gammapy.modeling.models.temporalModel`
          Temporal model
      name : str
          Model identifier
      apply_irf : dict
          Dictionary declaring which IRFs should be applied to this model. Options
          are {"exposure": True, "psf": True, "edisp": True}
      datasets_names : list of str
          Which datasets this model is applied to.
      """
  
      tag = "SkyModel"
      _apply_irf_default = {"exposure": True, "psf": True, "edisp": True}
  
      def __init__(
          self,
          spectral_model,
          spatial_model=None,
          temporal_model=None,
          name=None,
          apply_irf=None,
          datasets_names=None,
      ):
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          self.temporal_model = temporal_model
          self._name = make_name(name)
  
          if apply_irf is None:
              apply_irf = self._apply_irf_default.copy()
  
          self.apply_irf = apply_irf
          self.datasets_names = datasets_names
          self._check_unit()
          super().__init__()
  
      @property
      def _models(self):
          models = self.spectral_model, self.spatial_model, self.temporal_model
          return [model for model in models if model is not None]
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in self._models],
              )
  
      def _check_unit(self):
          from astropy.time import Time
          from gammapy.data.gti import GTI
  
          # evaluate over a test geom to check output unit
          # TODO simpler way to test this ?
          axis = MapAxis.from_edges(np.logspace(-1, 1, 3), unit=u.TeV, name="energy_true")
          geom = WcsGeom.create(skydir=(0, 0), npix=(2, 2), frame="galactic", axes=[axis])
          t_ref = Time(55555, format="mjd")
          gti = GTI.create([1, 5] * u.day, [2, 6] * u.day, reference_time=t_ref)
          value = self.evaluate_geom(geom, gti)
          if self.spatial_model is not None:
              ref_unit = "cm-2 s-1 MeV-1 sr-1"
          else:
              ref_unit = "cm-2 s-1 MeV-1"
          if not value.unit.is_equivalent(ref_unit):
              raise ValueError(
                  f"SkyModel unit {value.unit} is not equivalent to {ref_unit}"
              )
  
      @property
      def covariance(self):
          self._check_covariance()
  
          for model in self._models:
              self._covariance.set_subcovariance(model.covariance)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for model in self._models:
              subcovar = self._covariance.get_subcovariance(model.covariance.parameters)
              model.covariance = subcovar
  
      @property
      def name(self):
          return self._name
  
      @property
      def parameters(self):
          parameters = []
  
          parameters.append(self.spectral_model.parameters)
  
          if self.spatial_model is not None:
              parameters.append(self.spatial_model.parameters)
  
          if self.temporal_model is not None:
              parameters.append(self.temporal_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):
          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):
          if not (model is None or isinstance(model, SpectralModel)):
              raise TypeError(f"Invalid type: {model!r}")
          self._spectral_model = model
  
      @property
      def temporal_model(self):
          """`~gammapy.modeling.models.TemporalModel`"""
          return self._temporal_model
  
      @temporal_model.setter
      def temporal_model(self, model):
          if not (model is None or isinstance(model, TemporalModel)):
              raise TypeError(f"Invalid type: {model!r}")
  
          self._temporal_model = model
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.spatial_model.position
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return self.spatial_model.evaluation_radius
  
      @property
      def frame(self):
          return self.spatial_model.frame
  
      def __repr__(self):
          return (
              f"{self.__class__.__name__}("
              f"spatial_model={self.spatial_model!r}, "
              f"spectral_model={self.spectral_model!r})"
              f"temporal_model={self.temporal_model!r})"
          )
  
      def evaluate(self, lon, lat, energy, time=None):
          """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
          time: `~astropy.time.Time`
              Time co-ordinate
  
          Returns
          -------
          value : `~astropy.units.Quantity`
              Model value at the given point.
          """
          value = self.spectral_model(energy)  # pylint:disable=not-callable
          # TODO: case if self.temporal_model is not None, introduce time in arguments ?
  
          if self.spatial_model is not None:
              if self.spatial_model.is_energy_dependent:
                  spatial = self.spatial_model(lon, lat, energy)
              else:
                  spatial = self.spatial_model(lon, lat)
  
              value = value * spatial  # pylint:disable=not-callable
  
          if (self.temporal_model is not None) and (time is not None):
              value = value * self.temporal_model(time)
  
          return value
  
      def evaluate_geom(self, geom, gti=None):
          """Evaluate model on `~gammapy.maps.Geom`."""
          energy = geom.axes["energy_true"].center[:, np.newaxis, np.newaxis]
          value = self.spectral_model(energy)
  
          if self.spatial_model:
              value = value * self.spatial_model.evaluate_geom(geom)
  
          if self.temporal_model:
              integral = self.temporal_model.integral(gti.time_start, gti.time_stop)
              value = value * np.sum(integral)
  
          return value
  
      def integrate_geom(self, geom, gti=None):
          """Integrate model on `~gammapy.maps.Geom`.
  
          Parameters
          ----------
          geom : `Geom`
              Map geometry
          gti : `GTI`
              GIT table
  
          Returns
          -------
          flux : `Map`
              Predicted flux map
          """
          energy = geom.axes["energy_true"].edges
          value = self.spectral_model.integral(
              energy[:-1], energy[1:], intervals=True
          ).reshape((-1, 1, 1))
  
          if self.spatial_model and not isinstance(geom, RegionGeom):
              # TODO: integrate spatial model over region to correct for
              #  containment
              value = value * self.spatial_model.integrate_geom(geom).quantity
  
          if self.temporal_model:
              integral = self.temporal_model.integral(gti.time_start, gti.time_stop)
              value = value * np.sum(integral)
  
          return Map.from_geom(geom=geom, data=value.value, unit=value.unit)
  
      def copy(self, name=None, **kwargs):
          """Copy SkyModel"""
          if self.spatial_model is not None:
              spatial_model = self.spatial_model.copy()
          else:
              spatial_model = None
  
          if self.temporal_model is not None:
              temporal_model = self.temporal_model.copy()
          else:
              temporal_model = None
  
          kwargs.setdefault("name", make_name(name))
          kwargs.setdefault("spectral_model", self.spectral_model.copy())
          kwargs.setdefault("spatial_model", spatial_model)
          kwargs.setdefault("temporal_model", temporal_model)
          kwargs.setdefault("apply_irf", self.apply_irf.copy())
          kwargs.setdefault("datasets_names", self.datasets_names)
  
          return self.__class__(**kwargs)
  
      def to_dict(self, full_output=True):
          """Create dict for YAML serilisation"""
          data = {}
          data["name"] = self.name
          data["type"] = self.tag
          data["spectral"] = self.spectral_model.to_dict(full_output)
  
          if self.spatial_model is not None:
              data["spatial"] = self.spatial_model.to_dict(full_output)
  
          if self.temporal_model is not None:
              data["temporal"] = self.temporal_model.to_dict(full_output)
  
          if self.apply_irf != self._apply_irf_default:
              data["apply_irf"] = self.apply_irf
  
          if self.datasets_names is not None:
              data["datasets_names"] = self.datasets_names
  
          return data
  
      @classmethod
      def from_dict(cls, data):
          """Create SkyModel from dict"""
          from gammapy.modeling.models import (
              SPATIAL_MODEL_REGISTRY,
              SPECTRAL_MODEL_REGISTRY,
              TEMPORAL_MODEL_REGISTRY,
          )
  
          model_class = SPECTRAL_MODEL_REGISTRY.get_cls(data["spectral"]["type"])
          spectral_model = model_class.from_dict(data["spectral"])
  
          spatial_data = data.get("spatial")
  
          if spatial_data is not None:
              model_class = SPATIAL_MODEL_REGISTRY.get_cls(spatial_data["type"])
              spatial_model = model_class.from_dict(spatial_data)
          else:
              spatial_model = None
  
          temporal_data = data.get("temporal")
  
          if temporal_data is not None:
              model_class = TEMPORAL_MODEL_REGISTRY.get_cls(temporal_data["type"])
              temporal_model = model_class.from_dict(temporal_data)
          else:
              temporal_model = None
  
          return cls(
              name=data["name"],
              spatial_model=spatial_model,
              spectral_model=spectral_model,
              temporal_model=temporal_model,
              apply_irf=data.get("apply_irf", cls._apply_irf_default),
              datasets_names=data.get("datasets_names"),
          )
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
  
          str_ += "\t{:26}: {}\n".format("Datasets names", self.datasets_names)
  
          str_ += "\t{:26}: {}\n".format(
              "Spectral model type", self.spectral_model.__class__.__name__
          )
  
          if self.spatial_model is not None:
              spatial_type = self.spatial_model.__class__.__name__
          else:
              spatial_type = ""
          str_ += "\t{:26}: {}\n".format("Spatial  model type", spatial_type)
  
          if self.temporal_model is not None:
              temporal_type = self.temporal_model.__class__.__name__
          else:
              temporal_type = ""
          str_ += "\t{:26}: {}\n".format("Temporal model type", temporal_type)
  
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/cube.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Cube models (axes: lon, lat, energy)."""
  import copy
  from pathlib import Path
  import numpy as np
  import astropy.units as u
  from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.modeling.parameter import _get_parameters_str
  from gammapy.utils.scripts import make_name, make_path
  from gammapy.utils.fits import LazyFitsData, HDULocation
  from .core import Model, Models
  from .spatial import SpatialModel
  from .spectral import SpectralModel, PowerLawNormSpectralModel
  from .temporal import TemporalModel
  
  
  class SkyModelBase(Model):
      """Sky model base class"""
  
      def __add__(self, other):
          if isinstance(other, (Models, list)):
              return Models([self, *other])
          elif isinstance(other, (SkyModel, BackgroundModel)):
              return Models([self, other])
          else:
              raise TypeError(f"Invalid type: {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __call__(self, lon, lat, energy, time=None):
          return self.evaluate(lon, lat, energy, time)
  
      def evaluate_geom(self, geom, gti=None):
          coords = geom.get_coord(frame=self.frame)
          return self(coords.lon, coords.lat, coords["energy_true"])
  
  
  class SkyModel(SkyModelBase):
      """Sky model component.
  
      This model represents a factorised sky model.
      It has `~gammapy.modeling.Parameters`
      combining the spatial and spectral parameters.
  
      Parameters
      ----------
      spectral_model : `~gammapy.modeling.models.SpectralModel`
          Spectral model
      spatial_model : `~gammapy.modeling.models.SpatialModel`
          Spatial model (must be normalised to integrate to 1)
      temporal_model : `~gammapy.modeling.models.temporalModel`
          Temporal model
      name : str
          Model identifier
      apply_irf : dict
          Dictionary declaring which IRFs should be applied to this model. Options
          are {"exposure": True, "psf": True, "edisp": True}
      datasets_names : list of str
          Which datasets this model is applied to.
      """
  
      tag = "SkyModel"
      _apply_irf_default = {"exposure": True, "psf": True, "edisp": True}
  
      def __init__(
          self,
          spectral_model,
          spatial_model=None,
          temporal_model=None,
          name=None,
          apply_irf=None,
          datasets_names=None,
      ):
          self.spatial_model = spatial_model
          self.spectral_model = spectral_model
          self.temporal_model = temporal_model
          self._name = make_name(name)
  
          if apply_irf is None:
              apply_irf = self._apply_irf_default.copy()
  
          self.apply_irf = apply_irf
          self.datasets_names = datasets_names
          self._check_unit()
          super().__init__()
  
      @property
      def _models(self):
          models = self.spectral_model, self.spatial_model, self.temporal_model
          return [model for model in models if model is not None]
  
      def _check_covariance(self):
          if not self.parameters == self._covariance.parameters:
              self._covariance = Covariance.from_stack(
                  [model.covariance for model in self._models],
              )
  
      def _check_unit(self):
          from astropy.time import Time
          from gammapy.data.gti import GTI
  
          # evaluate over a test geom to check output unit
          # TODO simpler way to test this ?
          axis = MapAxis.from_edges(np.logspace(-1, 1, 3), unit=u.TeV, name="energy_true")
          geom = WcsGeom.create(skydir=(0, 0), npix=(2, 2), frame="galactic", axes=[axis])
          t_ref = Time(55555, format="mjd")
          gti = GTI.create([1, 5] * u.day, [2, 6] * u.day, reference_time=t_ref)
          value = self.evaluate_geom(geom, gti)
          if self.spatial_model is not None:
              ref_unit = "cm-2 s-1 MeV-1 sr-1"
          else:
              ref_unit = "cm-2 s-1 MeV-1"
          if not value.unit.is_equivalent(ref_unit):
              raise ValueError(
                  f"SkyModel unit {value.unit} is not equivalent to {ref_unit}"
              )
  
      @property
      def covariance(self):
          self._check_covariance()
  
          for model in self._models:
              self._covariance.set_subcovariance(model.covariance)
  
          return self._covariance
  
      @covariance.setter
      def covariance(self, covariance):
          self._check_covariance()
          self._covariance.data = covariance
  
          for model in self._models:
              subcovar = self._covariance.get_subcovariance(model.covariance.parameters)
              model.covariance = subcovar
  
      @property
      def name(self):
          return self._name
  
      @property
      def parameters(self):
          parameters = []
  
          parameters.append(self.spectral_model.parameters)
  
          if self.spatial_model is not None:
              parameters.append(self.spatial_model.parameters)
  
          if self.temporal_model is not None:
              parameters.append(self.temporal_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):
          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):
          if not (model is None or isinstance(model, SpectralModel)):
              raise TypeError(f"Invalid type: {model!r}")
          self._spectral_model = model
  
      @property
      def temporal_model(self):
          """`~gammapy.modeling.models.TemporalModel`"""
          return self._temporal_model
  
      @temporal_model.setter
      def temporal_model(self, model):
          if not (model is None or isinstance(model, TemporalModel)):
              raise TypeError(f"Invalid type: {model!r}")
  
          self._temporal_model = model
  
      @property
      def position(self):
          """`~astropy.coordinates.SkyCoord`"""
          return self.spatial_model.position
  
      @property
      def evaluation_radius(self):
          """`~astropy.coordinates.Angle`"""
          return self.spatial_model.evaluation_radius
  
      @property
      def frame(self):
          return self.spatial_model.frame
  
      def __repr__(self):
          return (
              f"{self.__class__.__name__}("
              f"spatial_model={self.spatial_model!r}, "
              f"spectral_model={self.spectral_model!r})"
              f"temporal_model={self.temporal_model!r})"
          )
  
      def evaluate(self, lon, lat, energy, time=None):
          """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
          time: `~astropy.time.Time`
              Time co-ordinate
  
          Returns
          -------
          value : `~astropy.units.Quantity`
              Model value at the given point.
          """
          value = self.spectral_model(energy)  # pylint:disable=not-callable
          # TODO: case if self.temporal_model is not None, introduce time in arguments ?
  
          if self.spatial_model is not None:
              if self.spatial_model.is_energy_dependent:
                  spatial = self.spatial_model(lon, lat, energy)
              else:
                  spatial = self.spatial_model(lon, lat)
  
              value = value * spatial  # pylint:disable=not-callable
  
          if (self.temporal_model is not None) and (time is not None):
              value = value * self.temporal_model(time)
  
          return value
  
      def evaluate_geom(self, geom, gti=None):
          """Evaluate model on `~gammapy.maps.Geom`."""
          energy = geom.axes["energy_true"].center[:, np.newaxis, np.newaxis]
          value = self.spectral_model(energy)
  
          if self.spatial_model:
              value = value * self.spatial_model.evaluate_geom(geom)
  
          if self.temporal_model:
              integral = self.temporal_model.integral(gti.time_start, gti.time_stop)
              value = value * np.sum(integral)
  
          return value
  
      def integrate_geom(self, geom, gti=None):
          """Integrate model on `~gammapy.maps.Geom`.
  
          Parameters
          ----------
          geom : `Geom`
              Map geometry
          gti : `GTI`
              GIT table
  
          Returns
          -------
          flux : `Map`
              Predicted flux map
          """
          energy = geom.axes["energy_true"].edges
          value = self.spectral_model.integral(
              energy[:-1], energy[1:], intervals=True
          ).reshape((-1, 1, 1))
  
          if self.spatial_model and not isinstance(geom, RegionGeom):
              # TODO: integrate spatial model over region to correct for
              #  containment
              value = value * self.spatial_model.integrate_geom(geom).quantity
  
          if self.temporal_model:
              integral = self.temporal_model.integral(gti.time_start, gti.time_stop)
              value = value * np.sum(integral)
  
          return Map.from_geom(geom=geom, data=value.value, unit=value.unit)
  
      def copy(self, name=None, **kwargs):
          """Copy SkyModel"""
          if self.spatial_model is not None:
              spatial_model = self.spatial_model.copy()
          else:
              spatial_model = None
  
          if self.temporal_model is not None:
              temporal_model = self.temporal_model.copy()
          else:
              temporal_model = None
  
          kwargs.setdefault("name", make_name(name))
          kwargs.setdefault("spectral_model", self.spectral_model.copy())
          kwargs.setdefault("spatial_model", spatial_model)
          kwargs.setdefault("temporal_model", temporal_model)
          kwargs.setdefault("apply_irf", self.apply_irf.copy())
          kwargs.setdefault("datasets_names", self.datasets_names)
  
          return self.__class__(**kwargs)
  
      def to_dict(self):
          """Create dict for YAML serilisation"""
          data = {}
          data["name"] = self.name
          data["type"] = self.tag
          data["spectral"] = self.spectral_model.to_dict()
  
          if self.spatial_model is not None:
              data["spatial"] = self.spatial_model.to_dict()
  
          if self.temporal_model is not None:
              data["temporal"] = self.temporal_model.to_dict()
  
          if self.apply_irf != self._apply_irf_default:
              data["apply_irf"] = self.apply_irf
  
          if self.datasets_names is not None:
              data["datasets_names"] = self.datasets_names
  
          return data
  
      @classmethod
      def from_dict(cls, data):
          """Create SkyModel from dict"""
          from gammapy.modeling.models import (
              SPATIAL_MODEL_REGISTRY,
              SPECTRAL_MODEL_REGISTRY,
              TEMPORAL_MODEL_REGISTRY,
          )
  
          model_class = SPECTRAL_MODEL_REGISTRY.get_cls(data["spectral"]["type"])
          spectral_model = model_class.from_dict(data["spectral"])
  
          spatial_data = data.get("spatial")
  
          if spatial_data is not None:
              model_class = SPATIAL_MODEL_REGISTRY.get_cls(spatial_data["type"])
              spatial_model = model_class.from_dict(spatial_data)
          else:
              spatial_model = None
  
          temporal_data = data.get("temporal")
  
          if temporal_data is not None:
              model_class = TEMPORAL_MODEL_REGISTRY.get_cls(temporal_data["type"])
              temporal_model = model_class.from_dict(temporal_data)
          else:
              temporal_model = None
  
          return cls(
              name=data["name"],
              spatial_model=spatial_model,
              spectral_model=spectral_model,
              temporal_model=temporal_model,
              apply_irf=data.get("apply_irf", cls._apply_irf_default),
              datasets_names=data.get("datasets_names"),
          )
  
      def __str__(self):
          str_ = self.__class__.__name__ + "\n\n"
          str_ += "\t{:26}: {}\n".format("Name", self.name)
  
          str_ += "\t{:26}: {}\n".format("Datasets names", self.datasets_names)
  
          str_ += "\t{:26}: {}\n".format(
              "Spectral model type", self.spectral_model.__class__.__name__
          )
  
          if self.spatial_model is not None:
              spatial_type = self.spatial_model.__class__.__name__
          else:
              spatial_type = ""
          str_ += "\t{:26}: {}\n".format("Spatial  model type", spatial_type)
  
          if self.temporal_model is not None:
              temporal_type = self.temporal_model.__class__.__name__
          else:
              temporal_type = ""
          str_ += "\t{:26}: {}\n".format("Temporal model type", temporal_type)
  
          str_ += "\tParameters:\n"
          info = _get_parameters_str(self.parameters)
          lines = info.split("\n")
          str
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/cube.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spatial models."""
  import logging
  import numpy as np
  import scipy.integrate
  import scipy.special
  import astropy.units as u
  from astropy.utils import lazyproperty
  from astropy.coordinates import Angle, SkyCoord
  from astropy.coordinates.angle_utilities import angular_separation, position_angle
  from regions import (
      CircleAnnulusSkyRegion,
      EllipseSkyRegion,
      PointSkyRegion,
      PolygonSkyRegion,
  )
  from gammapy.maps import Map, WcsGeom
  from gammapy.modeling import Parameter
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  log = logging.getLogger(__name__)
  
  
  def compute_sigma_eff(lon_0, lat_0, lon, lat, phi, major_axis, e):
      """Effective radius, used for the evaluation of elongated models"""
      phi_0 = position_angle(lon_0, lat_0, lon, lat)
      d_phi = phi - phi_0
      minor_axis = Angle(major_axis * np.sqrt(1 - e ** 2))
  
      a2 = (major_axis * np.sin(d_phi)) ** 2
      b2 = (minor_axis * np.cos(d_phi)) ** 2
      denominator = np.sqrt(a2 + b2)
      sigma_eff = major_axis * minor_axis / denominator
      return minor_axis, sigma_eff
  
  
  class SpatialModel(Model):
      """Spatial model base class."""
  
      _type = "spatial"
  
      def __init__(self, **kwargs):
          frame = kwargs.pop("frame", "icrs")
          super().__init__(**kwargs)
          if not hasattr(self, "frame"):
              self.frame = frame
  
      def __call__(self, lon, lat, energy=None):
          """Call evaluate method"""
          kwargs = {par.name: par.quantity for par in self.parameters}
  
          if energy is None and self.is_energy_dependent:
              raise ValueError("Missing energy value for evaluation")
  
          if energy is not None:
              kwargs["energy"] = energy
  
          return self.evaluate(lon, lat, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      # TODO: make this a hard-coded class attribute?
      @lazyproperty
      def is_energy_dependent(self):
          varnames = self.evaluate.__code__.co_varnames
          return "energy" in varnames
  
      @property
      def position(self):
          """Spatial model center position"""
          lon = self.lon_0.quantity
          lat = self.lat_0.quantity
          return SkyCoord(lon, lat, frame=self.frame)
  
      @position.setter
      def position(self, skycoord):
          """Spatial model center position"""
          coord = skycoord.transform_to(self.frame)
          self.lon_0.quantity = coord.data.lon
          self.lat_0.quantity = coord.data.lat
  
      # TODO: get rid of this!
      _phi_0 = 0.0
  
      @property
      def phi_0(self):
          return self._phi_0
  
      @phi_0.setter
      def phi_0(self, phi_0=0.0):
          self._phi_0 = phi_0
  
      @property
      def position_error(self):
          """Get 95% containment position error as (`~regions.EllipseSkyRegion`)"""
          if self.covariance is None:
              return EllipseSkyRegion(
                  center=self.position,
                  height=np.nan * u.deg,
                  width=np.nan * u.deg,
                  angle=np.nan * u.deg,
              )
          pars = self.parameters
          sub_covar = self.covariance.get_subcovariance(["lon_0", "lat_0"]).data.copy()
          cos_lat = np.cos(self.lat_0.quantity.to_value("rad"))
          sub_covar[0, 0] *= cos_lat ** 2.0
          sub_covar[0, 1] *= cos_lat
          sub_covar[1, 0] *= cos_lat
          eig_vals, eig_vecs = np.linalg.eig(sub_covar)
          lon_err, lat_err = np.sqrt(eig_vals)
          y_vec = eig_vecs[:, 0]
          phi = (np.arctan2(y_vec[1], y_vec[0]) * u.rad).to("deg") + self.phi_0
          err = np.sort([lon_err, lat_err])
          scale_r95 = Gauss2DPDF().containment_radius(0.95)
          err *= scale_r95
          if err[1] == lon_err * scale_r95:
              phi += 90 * u.deg
              height = 2 * err[1] * pars["lon_0"].unit
              width = 2 * err[0] * pars["lat_0"].unit
          else:
              height = 2 * err[1] * pars["lat_0"].unit
              width = 2 * err[0] * pars["lon_0"].unit
          return EllipseSkyRegion(
              center=self.position, height=height, width=width, angle=phi
          )
  
      def evaluate_geom(self, geom):
          coords = geom.to_image().get_coord(frame=self.frame)
  
          if self.is_energy_dependent:
              energy = geom.axes["energy_true"].center
              return self(coords.lon, coords.lat, energy[:, np.newaxis, np.newaxis])
          else:
              return self(coords.lon, coords.lat)
  
      def integrate_geom(self, geom):
          """Integrate model on `~gammapy.maps.Geom`."""
          values = self.evaluate_geom(geom)
          data = values * geom.solid_angle()
          return Map.from_geom(geom=geom, data=data.value, unit=data.unit)
  
      def to_dict(self, full_output=True):
          """Create dict for YAML serilisation"""
          data = super().to_dict(full_output)
          data["frame"] = self.frame
          data["parameters"] = data.pop("parameters")
          return data
  
      def _get_plot_map(self, geom):
          if self.evaluation_radius is None and geom is None:
              raise ValueError(
                  f"{self.__class__.__name__} requires geom to be defined for plotting."
              )
  
          if geom is None:
              width = 2 * max(self.evaluation_radius, 0.1 * u.deg)
              geom = WcsGeom.create(
                  skydir=self.position, frame=self.frame, width=width, binsz=0.02
              )
          data = self.evaluate_geom(geom)
          return Map.from_geom(geom, data=data.value, unit=data.unit)
  
      def plot(self, ax=None, geom=None, **kwargs):
          """Plot spatial model.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          geom : `~gammapy.maps.WcsGeom`, optional
              Geom to use for plotting.
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          m = self._get_plot_map(geom)
          if not m.geom.is_flat:
              raise TypeError("Use .plot_interactive() for Map dimension > 2")
          _, ax, _ = m.plot(ax=ax, **kwargs)
          return ax
  
      def plot_interative(self, ax=None, geom=None, **kwargs):
          """Plot spatial model.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          geom : `~gammapy.maps.WcsGeom`, optional
              Geom to use for plotting.
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
  
          m = self._get_plot_map(geom)
          if m.geom.is_image:
              raise TypeError("Use .plot() for 2D Maps")
          m.plot_interactive(ax=ax, **kwargs)
  
      def plot_error(self, ax=None, **kwargs):
          """Plot position error
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          # plot center position
          lon, lat = self.lon_0.value, self.lat_0.value
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("marker", "x")
          kwargs.setdefault("color", "red")
          kwargs.setdefault("label", "position")
  
          ax.scatter(lon, lat, transform=ax.get_transform(self.frame), **kwargs)
  
          # plot position error
          if not np.all(self.covariance.data == 0):
              region = self.position_error.to_pixel(ax.wcs)
              artist = region.as_artist(facecolor="none", edgecolor=kwargs["color"])
              ax.add_artist(artist)
  
          return ax
  
  
  class PointSpatialModel(SpatialModel):
      r"""Point Source.
  
      For more information see :ref:`point-spatial-model`.
  
      Parameters
      ----------
      lon_0, lat_0 : `~astropy.coordinates.Angle`
          Center position
      frame : {"icrs", "galactic"}
          Center position coordinate frame
      """
  
      tag = ["PointSpatialModel", "point"]
      lon_0 = Parameter("lon_0", "0 deg")
      lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
  
      @property
      def evaluation_radius(self):
          """Evaluation radius (`~astropy.coordinates.Angle`).
  
          Set as zero degrees.
          """
          return 0 * u.deg
  
      @staticmethod
      def _grid_weights(x, y, x0, y0):
          """Compute 4-pixel weights such that centroid is preserved."""
          dx = np.abs(x - x0)
          dx = np.where(dx < 1, 1 - dx, 0)
  
          dy = np.abs(y - y0)
          dy = np.where(dy < 1, 1 - dy, 0)
  
          return dx * dy
  
      def evaluate_geom(self, geom):
          """Evaluate model on `~gammapy.maps.Geom`."""
          values = self.integrate_geom(geom).data
          return values / geom.solid_angle()
  
      def integrate_geom(self, geom):
          """Integrate model on `~gammapy.maps.Geom`
  
          Parameters
          ----------
          geom : `Geom`
              Map geometry
  
          Returns
          -------
          flux : `Map`
              Predicted flux map
          """
          x, y = geom.get_pix()[0:2]
          x0, y0 = self.position.to_pixel(geom.wcs)
          data = self._grid_weights(x, y, x0, y0)
          return Map.from_geom(geom=geom, data=data, unit="")
  
      def to_region(self, **kwargs):
          """Model outline (`~regions.PointSkyRegion`)."""
          return PointSkyRegion(center=self.position, **kwargs)
  
  
  class GaussianSpatialModel(SpatialModel):
      r"""Two-dimensional Gaussian model.
  
      For more information see :ref:`gaussian-spatial-model`.
  
      Parameters
      ----------
      lon_0, lat_0 : `~astropy.coordinates.Angle`
          Center position
      sigma : `~astropy.coordinates.Angle`
          Length of the major semiaxis of the Gaussian, in angular units.
      e : `float`
          Eccentricity of the Gaussian (:math:`0< e< 1`).
      phi : `~astropy.coordinates.Angle`
          Rotation angle :math:`\phi`: of the major semiaxis.
          Increases counter-clockwise from the North direction.
      frame : {"icrs", "galactic"}
          Center position coordinate frame
      """
  
      tag = ["GaussianSpatialModel", "gauss"]
  
      lon_0 = Parameter("lon_0", "0 deg")
      lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
      sigma = Parameter("sigma", "1 deg", min=0)
      e = Parameter("e", 0, min=0, max=1, frozen=True)
      phi = Parameter("phi", "0 deg", frozen=True)
  
      @property
      def evaluation_radius(self):
          r"""Evaluation radius (`~astropy.coordinates.Angle`).
  
          Set as :math:`5\sigma`.
          """
          return 5 * self.parameters["sigma"].quantity
  
      @staticmethod
      def evaluate(lon, lat, lon_0, lat_0, sigma, e, phi):
          """Evaluate model."""
          sep = angular_separation(lon, lat, lon_0, lat_0)
  
          if e == 0:
              a = 1.0 - np.cos(sigma)
              norm = (1 / (4 * np.pi * a * (1.0 - np.exp(-1.0 / a)))).value
          else:
              minor_axis, sigma_eff = compute_sigma_eff(
                  lon_0, lat_0, lon, lat, phi, sigma, e
              )
              a = 1.0 - np.cos(sigma_eff)
              norm = (1 / (2 * np.pi * sigma * minor_axis)).to_value("sr-1")
  
          exponent = -0.5 * ((1 - np.cos(sep)) / a)
          return u.Quantity(norm * np.exp(exponent).value, "sr-1", copy=False)
  
      def to_region(self, **kwargs):
          """Model outline (`
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/spatial.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spatial models."""
  import logging
  import numpy as np
  import scipy.integrate
  import scipy.special
  import astropy.units as u
  from astropy.utils import lazyproperty
  from astropy.coordinates import Angle, SkyCoord
  from astropy.coordinates.angle_utilities import angular_separation, position_angle
  from regions import (
      CircleAnnulusSkyRegion,
      EllipseSkyRegion,
      PointSkyRegion,
      PolygonSkyRegion,
  )
  from gammapy.maps import Map, WcsGeom
  from gammapy.modeling import Parameter
  from gammapy.utils.gauss import Gauss2DPDF
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  log = logging.getLogger(__name__)
  
  
  def compute_sigma_eff(lon_0, lat_0, lon, lat, phi, major_axis, e):
      """Effective radius, used for the evaluation of elongated models"""
      phi_0 = position_angle(lon_0, lat_0, lon, lat)
      d_phi = phi - phi_0
      minor_axis = Angle(major_axis * np.sqrt(1 - e ** 2))
  
      a2 = (major_axis * np.sin(d_phi)) ** 2
      b2 = (minor_axis * np.cos(d_phi)) ** 2
      denominator = np.sqrt(a2 + b2)
      sigma_eff = major_axis * minor_axis / denominator
      return minor_axis, sigma_eff
  
  
  class SpatialModel(Model):
      """Spatial model base class."""
      _type = "spatial"
  
      def __init__(self, **kwargs):
          frame = kwargs.pop("frame", "icrs")
          super().__init__(**kwargs)
          if not hasattr(self, "frame"):
              self.frame = frame
  
      def __call__(self, lon, lat, energy=None):
          """Call evaluate method"""
          kwargs = {par.name: par.quantity for par in self.parameters}
  
          if energy is None and self.is_energy_dependent:
              raise ValueError("Missing energy value for evaluation")
  
          if energy is not None:
              kwargs["energy"] = energy
  
          return self.evaluate(lon, lat, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      # TODO: make this a hard-coded class attribute?
      @lazyproperty
      def is_energy_dependent(self):
          varnames = self.evaluate.__code__.co_varnames
          return "energy" in varnames
  
      @property
      def position(self):
          """Spatial model center position"""
          lon = self.lon_0.quantity
          lat = self.lat_0.quantity
          return SkyCoord(lon, lat, frame=self.frame)
  
      @position.setter
      def position(self, skycoord):
          """Spatial model center position"""
          coord = skycoord.transform_to(self.frame)
          self.lon_0.quantity = coord.data.lon
          self.lat_0.quantity = coord.data.lat
  
      # TODO: get rid of this!
      _phi_0 = 0.0
  
      @property
      def phi_0(self):
          return self._phi_0
  
      @phi_0.setter
      def phi_0(self, phi_0=0.0):
          self._phi_0 = phi_0
  
      @property
      def position_error(self):
          """Get 95% containment position error as (`~regions.EllipseSkyRegion`)"""
          if self.covariance is None:
              return EllipseSkyRegion(
                  center=self.position,
                  height=np.nan * u.deg,
                  width=np.nan * u.deg,
                  angle=np.nan * u.deg,
              )
          pars = self.parameters
          sub_covar = self.covariance.get_subcovariance(["lon_0", "lat_0"]).data.copy()
          cos_lat = np.cos(self.lat_0.quantity.to_value("rad"))
          sub_covar[0, 0] *= cos_lat ** 2.0
          sub_covar[0, 1] *= cos_lat
          sub_covar[1, 0] *= cos_lat
          eig_vals, eig_vecs = np.linalg.eig(sub_covar)
          lon_err, lat_err = np.sqrt(eig_vals)
          y_vec = eig_vecs[:, 0]
          phi = (np.arctan2(y_vec[1], y_vec[0]) * u.rad).to("deg") + self.phi_0
          err = np.sort([lon_err, lat_err])
          scale_r95 = Gauss2DPDF().containment_radius(0.95)
          err *= scale_r95
          if err[1] == lon_err * scale_r95:
              phi += 90 * u.deg
              height = 2 * err[1] * pars["lon_0"].unit
              width = 2 * err[0] * pars["lat_0"].unit
          else:
              height = 2 * err[1] * pars["lat_0"].unit
              width = 2 * err[0] * pars["lon_0"].unit
          return EllipseSkyRegion(
              center=self.position, height=height, width=width, angle=phi
          )
  
      def evaluate_geom(self, geom):
          coords = geom.to_image().get_coord(frame=self.frame)
  
          if self.is_energy_dependent:
              energy = geom.axes["energy_true"].center
              return self(coords.lon, coords.lat, energy[:, np.newaxis, np.newaxis])
          else:
              return self(coords.lon, coords.lat)
  
      def integrate_geom(self, geom):
          """Integrate model on `~gammapy.maps.Geom`."""
          values = self.evaluate_geom(geom)
          data = values * geom.solid_angle()
          return Map.from_geom(geom=geom, data=data.value, unit=data.unit)
  
      def to_dict(self):
          """Create dict for YAML serilisation"""
          data = super().to_dict()
          data["frame"] = self.frame
          data["parameters"] = data.pop("parameters")
          return data
  
      def _get_plot_map(self, geom):
          if self.evaluation_radius is None and geom is None:
              raise ValueError(
                  f"{self.__class__.__name__} requires geom to be defined for plotting."
              )
  
          if geom is None:
              width = 2 * max(self.evaluation_radius, 0.1 * u.deg)
              geom = WcsGeom.create(
                  skydir=self.position, frame=self.frame, width=width, binsz=0.02
              )
          data = self.evaluate_geom(geom)
          return Map.from_geom(geom, data=data.value, unit=data.unit)
  
      def plot(self, ax=None, geom=None, **kwargs):
          """Plot spatial model.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          geom : `~gammapy.maps.WcsGeom`, optional
              Geom to use for plotting.
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          m = self._get_plot_map(geom)
          if not m.geom.is_flat:
              raise TypeError("Use .plot_interactive() for Map dimension > 2")
          _, ax, _ = m.plot(ax=ax, **kwargs)
          return ax
  
      def plot_interative(self, ax=None, geom=None, **kwargs):
          """Plot spatial model.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          geom : `~gammapy.maps.WcsGeom`, optional
              Geom to use for plotting.
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
  
          m = self._get_plot_map(geom)
          if m.geom.is_image:
              raise TypeError("Use .plot() for 2D Maps")
          m.plot_interactive(ax=ax, **kwargs)
  
      def plot_error(self, ax=None, **kwargs):
          """Plot position error
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          **kwargs : dict
              Keyword arguments passed to `~gammapy.maps.WcsMap.plot()`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          # plot center position
          lon, lat = self.lon_0.value, self.lat_0.value
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("marker", "x")
          kwargs.setdefault("color", "red")
          kwargs.setdefault("label", "position")
  
          ax.scatter(lon, lat, transform=ax.get_transform(self.frame), **kwargs)
  
          # plot position error
          if not np.all(self.covariance.data == 0):
              region = self.position_error.to_pixel(ax.wcs)
              artist = region.as_artist(facecolor="none", edgecolor=kwargs["color"])
              ax.add_artist(artist)
  
          return ax
  
  
  class PointSpatialModel(SpatialModel):
      r"""Point Source.
  
      For more information see :ref:`point-spatial-model`.
  
      Parameters
      ----------
      lon_0, lat_0 : `~astropy.coordinates.Angle`
          Center position
      frame : {"icrs", "galactic"}
          Center position coordinate frame
      """
  
      tag = ["PointSpatialModel", "point"]
      lon_0 = Parameter("lon_0", "0 deg")
      lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
  
      @property
      def evaluation_radius(self):
          """Evaluation radius (`~astropy.coordinates.Angle`).
  
          Set as zero degrees.
          """
          return 0 * u.deg
  
      @staticmethod
      def _grid_weights(x, y, x0, y0):
          """Compute 4-pixel weights such that centroid is preserved."""
          dx = np.abs(x - x0)
          dx = np.where(dx < 1, 1 - dx, 0)
  
          dy = np.abs(y - y0)
          dy = np.where(dy < 1, 1 - dy, 0)
  
          return dx * dy
  
      def evaluate_geom(self, geom):
          """Evaluate model on `~gammapy.maps.Geom`."""
          values = self.integrate_geom(geom).data
          return values / geom.solid_angle()
  
      def integrate_geom(self, geom):
          """Integrate model on `~gammapy.maps.Geom`
  
          Parameters
          ----------
          geom : `Geom`
              Map geometry
  
          Returns
          -------
          flux : `Map`
              Predicted flux map
          """
          x, y = geom.get_pix()[0:2]
          x0, y0 = self.position.to_pixel(geom.wcs)
          data = self._grid_weights(x, y, x0, y0)
          return Map.from_geom(geom=geom, data=data, unit="")
  
      def to_region(self, **kwargs):
          """Model outline (`~regions.PointSkyRegion`)."""
          return PointSkyRegion(center=self.position, **kwargs)
  
  
  class GaussianSpatialModel(SpatialModel):
      r"""Two-dimensional Gaussian model.
  
      For more information see :ref:`gaussian-spatial-model`.
  
      Parameters
      ----------
      lon_0, lat_0 : `~astropy.coordinates.Angle`
          Center position
      sigma : `~astropy.coordinates.Angle`
          Length of the major semiaxis of the Gaussian, in angular units.
      e : `float`
          Eccentricity of the Gaussian (:math:`0< e< 1`).
      phi : `~astropy.coordinates.Angle`
          Rotation angle :math:`\phi`: of the major semiaxis.
          Increases counter-clockwise from the North direction.
      frame : {"icrs", "galactic"}
          Center position coordinate frame
      """
  
      tag = ["GaussianSpatialModel", "gauss"]
  
      lon_0 = Parameter("lon_0", "0 deg")
      lat_0 = Parameter("lat_0", "0 deg", min=-90, max=90)
      sigma = Parameter("sigma", "1 deg", min=0)
      e = Parameter("e", 0, min=0, max=1, frozen=True)
      phi = Parameter("phi", "0 deg", frozen=True)
  
      @property
      def evaluation_radius(self):
          r"""Evaluation radius (`~astropy.coordinates.Angle`).
  
          Set as :math:`5\sigma`.
          """
          return 5 * self.parameters["sigma"].quantity
  
      @staticmethod
      def evaluate(lon, lat, lon_0, lat_0, sigma, e, phi):
          """Evaluate model."""
          sep = angular_separation(lon, lat, lon_0, lat_0)
  
          if e == 0:
              a = 1.0 - np.cos(sigma)
              norm = (1 / (4 * np.pi * a * (1.0 - np.exp(-1.0 / a)))).value
          else:
              minor_axis, sigma_eff = compute_sigma_eff(
                  lon_0, lat_0, lon, lat, phi, sigma, e
              )
              a = 1.0 - np.cos(sigma_eff)
              norm = (1 / (2 * np.pi * sigma * minor_axis)).to_value("sr-1")
  
          exponent = -0.5 * ((1 - np.cos(sep)) / a)
          return u.Quantity(norm * np.exp(exponent).value, "sr-1", copy=False)
  
      def to_region(self, **kwargs):
          """Model outline (`~regions.EllipseSkyRegion
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/spatial.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spectral models for Gammapy."""
  import operator
  import numpy as np
  import scipy.optimize
  import scipy.special
  import astropy.units as u
  from astropy import constants as const
  from astropy.table import Table
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.integrate import trapz_loglog
  from gammapy.utils.interpolation import ScaledRegularGridInterpolator
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  
  def integrate_spectrum(func, emin, emax, ndecade=100, intervals=False):
      """Integrate 1d function using the log-log trapezoidal rule.
  
      If scalar values for xmin and xmax are passed an oversampled grid is generated using the
      ``ndecade`` keyword argument. If xmin and xmax arrays are passed, no
      oversampling is performed and the integral is computed in the provided
      grid.
  
      Parameters
      ----------
      func : callable
          Function to integrate.
      emin : `~astropy.units.Quantity`
          Integration range minimum
      emax : `~astropy.units.Quantity`
          Integration range minimum
      ndecade : int, optional
          Number of grid points per decade used for the integration.
          Default : 100.
      intervals : bool, optional
          Return integrals in the grid not the sum, default: False
      """
      if emin.isscalar and emax.isscalar:
          energies = MapAxis.from_energy_bounds(
              emin=emin, emax=emax, nbin=ndecade, per_decade=True
          ).edges
      else:
          energies = edges_from_lo_hi(emin, emax)
  
      values = func(energies)
  
      integral = trapz_loglog(values, energies)
  
      if intervals:
          return integral
  
      return integral.sum()
  
  
  class SpectralModel(Model):
      """Spectral model base class."""
  
      _type = "spectral"
  
      def __call__(self, energy):
          kwargs = {par.name: par.quantity for par in self.parameters}
          kwargs = self._convert_evaluate_unit(kwargs, energy)
          return self.evaluate(energy, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def _convert_evaluate_unit(kwargs_ref, energy):
          kwargs = {}
          for name, quantity in kwargs_ref.items():
              if quantity.unit.physical_type == "energy":
                  quantity = quantity.to(energy.unit)
              kwargs[name] = quantity
          return kwargs
  
      def __add__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.add)
  
      def __mul__(self, other):
          if isinstance(other, SpectralModel):
              return CompoundSpectralModel(self, other, operator.mul)
          else:
              raise TypeError(f"Multiplication invalid for type {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __sub__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.sub)
  
      def __rsub__(self, model):
          return self.__sub__(model)
  
      def _evaluate_gradient(self, energy, eps):
          n = len(self.parameters)
          f = self(energy)
          shape = (n, len(np.atleast_1d(energy)))
          df_dp = np.zeros(shape)
  
          for idx, parameter in enumerate(self.parameters):
              if parameter.frozen or eps[idx] == 0:
                  continue
  
              parameter.value += eps[idx]
              df = self(energy) - f
              df_dp[idx] = df.value / eps[idx]
  
              # Reset model to original parameter
              parameter.value -= eps[idx]
  
          return df_dp
  
      def evaluate_error(self, energy, epsilon=1e-4):
          """Evaluate spectral model with error propagation.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to evaluate
          epsilon : float
              Step size of the gradient evaluation. Given as a
              fraction of the parameter error.
  
          Returns
          -------
          dnde, dnde_error : tuple of `~astropy.units.Quantity`
              Tuple of flux and flux error.
          """
          p_cov = self.covariance
          eps = np.sqrt(np.diag(p_cov)) * epsilon
  
          df_dp = self._evaluate_gradient(energy, eps)
          f_cov = df_dp.T @ p_cov @ df_dp
          f_err = np.sqrt(np.diagonal(f_cov))
  
          q = self(energy)
          return u.Quantity([q.value, f_err], unit=q.unit)
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate spectral model numerically if no analytical solution defined.
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE
  
          If array input for ``emin`` and ``emax`` is given you have to set
          ``intervals=True`` if you want the integral in each energy bin.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to :func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          if hasattr(self, "evaluate_integral"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_integral(emin, emax, **kwargs)
          else:
              return integrate_spectrum(self, emin, emax, **kwargs)
  
      def integral_error(self, emin, emax):
          """Evaluate the error of the integral flux of a given spectrum in
          a given energy range.
  
          Parameters
          ----------
          emin, emax :  `~astropy.units.Quantity`
              Lower and upper bound of integration range.
  
          Returns
          -------
          flux, flux_err : tuple of `~astropy.units.Quantity`
              Integral flux and flux error betwen emin and emax.
          """
          energy = np.sqrt(emin * emax)
          flux = self.integral(emin, emax, intervals=True)
          dnde, dnde_err = self.evaluate_error(energy, epsilon=1e-4)
          flux_err = flux * dnde_err / dnde
          return u.Quantity([flux.value, flux_err.value], unit=flux.unit)
  
      def energy_flux(self, emin, emax, **kwargs):
          r"""Compute energy flux in given energy range.
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          def f(x):
              return x * self(x)
  
          if hasattr(self, "evaluate_energy_flux"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_energy_flux(emin, emax, **kwargs)
          else:
              return integrate_spectrum(f, emin, emax, **kwargs)
  
      def plot(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model curve.
  
          kwargs are forwarded to `matplotlib.pyplot.plot`
  
          By default a log-log scaling of the axes is used, if you want to change
          the y axis scaling to linear you can use::
  
              from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel
              from astropy import units as u
  
              pwl = ExpCutoffPowerLawSpectralModel()
              ax = pwl.plot(energy_range=(0.1, 100) * u.TeV)
              ax.set_yscale('linear')
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          # evaluate model
          flux = self(energy).to(flux_unit)
  
          y = self._plot_scale_flux(energy, flux, energy_power)
  
          ax.plot(energy.value, y.value, **kwargs)
  
          self._plot_format_ax(ax, energy, y, energy_power)
          return ax
  
      def plot_error(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model error band.
  
          .. note::
  
              This method calls ``ax.set_yscale("log", nonposy='clip')`` and
              ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation.
              The additional argument ``nonposx='clip'`` avoids artefacts in the plot,
              when the error band extends to negative values (see also
              https://github.com/matplotlib/matplotlib/issues/8623).
  
              When you call ``plt.loglog()`` or ``plt.semilogy()`` explicitely in your
              plotting code and the error band extends to negative values, it is not
              shown correctly. To circumvent this issue also use
              ``plt.loglog(nonposx='clip', nonposy='clip')``
              or ``plt.semilogy(nonposy='clip')``.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
          **kwargs : dict
              Keyword arguments forwarded to `matplotlib.pyplot.fill_between`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("facecolor", "black")
          kwargs.setdefault("alpha", 0.2)
          kwargs.setdefault("linewidth", 0)
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          flux, flux_err = self.evaluate_error(energy).to(flux_unit)
  
          y_lo = self._plot_scale_flux(energy, flux - flux_err, energy_power)
          y_hi = self._plot_scale_flux(energy, flux + flux_err, energy_power)
  
          where = (energy >= energy_range[0]) & (energy <= energy_range[1])
          ax.fill_between(energy.value, y_lo.value, y_hi.value, where=where, **kwargs)
  
          self._plot_format_ax(ax, energy, y_lo, energy_power)
          return ax
  
      @staticmethod
      def _plot_format_ax(ax, energy, y, energy_power):
          ax.set_xlabel(f"Energy [{energy.unit}]")
          if energy_power > 0:
              ax.set_ylabel(f"E{energy_power} * Flux [{y.unit}]")
          else:
              ax.set_ylabel(f"Flux [{y.unit}]")
  
          ax.set_xscale("log", nonposx="clip")
          ax.set_yscale("log", nonposy="clip")
  
      @staticmethod
      def _plot_scale_flux(energy, flux, energy_power):
          try:
              eunit = [_ for _ in flux.unit.bases if _.physical_type == "energy"][0]
          except IndexError:
              eunit = energy.unit
          y = flux * np.power(energy, energy_power)
          return y.to(flux.unit * eunit ** energy_power)
  
      def spectral_index(self, energy, epsilon=1e-5):
          """Compute spectral index at given energy.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to estimate the index
          epsilon : float
              Fractional energy increment to use for determining the spectral index.
  
          Returns
          -------
          index : float
              Estimated spectral index.
          """
          f1 = self(energy)
          f2 = self(energy * (1 + epsilon))
          return np.log(f1 / f2) / np.log(1 + epsilon)
  
      def
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/spectral.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Spectral models for Gammapy."""
  import operator
  import numpy as np
  import scipy.optimize
  import scipy.special
  import astropy.units as u
  from astropy import constants as const
  from astropy.table import Table
  from gammapy.maps import MapAxis
  from gammapy.maps.utils import edges_from_lo_hi
  from gammapy.modeling import Covariance, Parameter, Parameters
  from gammapy.utils.integrate import trapz_loglog
  from gammapy.utils.interpolation import ScaledRegularGridInterpolator
  from gammapy.utils.scripts import make_path
  from .core import Model
  
  
  def integrate_spectrum(func, emin, emax, ndecade=100, intervals=False):
      """Integrate 1d function using the log-log trapezoidal rule.
  
      If scalar values for xmin and xmax are passed an oversampled grid is generated using the
      ``ndecade`` keyword argument. If xmin and xmax arrays are passed, no
      oversampling is performed and the integral is computed in the provided
      grid.
  
      Parameters
      ----------
      func : callable
          Function to integrate.
      emin : `~astropy.units.Quantity`
          Integration range minimum
      emax : `~astropy.units.Quantity`
          Integration range minimum
      ndecade : int, optional
          Number of grid points per decade used for the integration.
          Default : 100.
      intervals : bool, optional
          Return integrals in the grid not the sum, default: False
      """
      if emin.isscalar and emax.isscalar:
          energies = MapAxis.from_energy_bounds(
              emin=emin, emax=emax, nbin=ndecade, per_decade=True
          ).edges
      else:
          energies = edges_from_lo_hi(emin, emax)
  
      values = func(energies)
  
      integral = trapz_loglog(values, energies)
  
      if intervals:
          return integral
  
      return integral.sum()
  
  
  class SpectralModel(Model):
      """Spectral model base class."""
      _type = "spectral"
  
      def __call__(self, energy):
          kwargs = {par.name: par.quantity for par in self.parameters}
          kwargs = self._convert_evaluate_unit(kwargs, energy)
          return self.evaluate(energy, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def _convert_evaluate_unit(kwargs_ref, energy):
          kwargs = {}
          for name, quantity in kwargs_ref.items():
              if quantity.unit.physical_type == "energy":
                  quantity = quantity.to(energy.unit)
              kwargs[name] = quantity
          return kwargs
  
      def __add__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.add)
  
      def __mul__(self, other):
          if isinstance(other, SpectralModel):
              return CompoundSpectralModel(self, other, operator.mul)
          else:
              raise TypeError(f"Multiplication invalid for type {other!r}")
  
      def __radd__(self, model):
          return self.__add__(model)
  
      def __sub__(self, model):
          if not isinstance(model, SpectralModel):
              model = ConstantSpectralModel(const=model)
          return CompoundSpectralModel(self, model, operator.sub)
  
      def __rsub__(self, model):
          return self.__sub__(model)
  
      def _evaluate_gradient(self, energy, eps):
          n = len(self.parameters)
          f = self(energy)
          shape = (n, len(np.atleast_1d(energy)))
          df_dp = np.zeros(shape)
  
          for idx, parameter in enumerate(self.parameters):
              if parameter.frozen or eps[idx] == 0:
                  continue
  
              parameter.value += eps[idx]
              df = self(energy) - f
              df_dp[idx] = df.value / eps[idx]
  
              # Reset model to original parameter
              parameter.value -= eps[idx]
  
          return df_dp
  
      def evaluate_error(self, energy, epsilon=1e-4):
          """Evaluate spectral model with error propagation.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to evaluate
          epsilon : float
              Step size of the gradient evaluation. Given as a
              fraction of the parameter error.
  
          Returns
          -------
          dnde, dnde_error : tuple of `~astropy.units.Quantity`
              Tuple of flux and flux error.
          """
          p_cov = self.covariance
          eps = np.sqrt(np.diag(p_cov)) * epsilon
  
          df_dp = self._evaluate_gradient(energy, eps)
          f_cov = df_dp.T @ p_cov @ df_dp
          f_err = np.sqrt(np.diagonal(f_cov))
  
          q = self(energy)
          return u.Quantity([q.value, f_err], unit=q.unit)
  
      def integral(self, emin, emax, **kwargs):
          r"""Integrate spectral model numerically if no analytical solution defined.
  
          .. math::
              F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E) dE
  
          If array input for ``emin`` and ``emax`` is given you have to set
          ``intervals=True`` if you want the integral in each energy bin.
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to :func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          if hasattr(self, "evaluate_integral"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_integral(emin, emax, **kwargs)
          else:
              return integrate_spectrum(self, emin, emax, **kwargs)
  
      def integral_error(self, emin, emax):
          """Evaluate the error of the integral flux of a given spectrum in
          a given energy range.
  
          Parameters
          ----------
          emin, emax :  `~astropy.units.Quantity`
              Lower and upper bound of integration range.
  
          Returns
          -------
          flux, flux_err : tuple of `~astropy.units.Quantity`
              Integral flux and flux error betwen emin and emax.
          """
          energy = np.sqrt(emin * emax)
          flux = self.integral(emin, emax, intervals=True)
          dnde, dnde_err = self.evaluate_error(energy, epsilon=1e-4)
          flux_err = flux * dnde_err / dnde
          return u.Quantity([flux.value, flux_err.value], unit=flux.unit)
  
      def energy_flux(self, emin, emax, **kwargs):
          r"""Compute energy flux in given energy range.
  
          .. math::
              G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} E \phi(E) dE
  
          Parameters
          ----------
          emin, emax : `~astropy.units.Quantity`
              Lower and upper bound of integration range.
          **kwargs : dict
              Keyword arguments passed to func:`~gammapy.utils.integrate.integrate_spectrum`
          """
  
          def f(x):
              return x * self(x)
  
          if hasattr(self, "evaluate_energy_flux"):
              kwargs = {par.name: par.quantity for par in self.parameters}
              kwargs = self._convert_evaluate_unit(kwargs, emin)
              return self.evaluate_energy_flux(emin, emax, **kwargs)
          else:
              return integrate_spectrum(f, emin, emax, **kwargs)
  
      def plot(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model curve.
  
          kwargs are forwarded to `matplotlib.pyplot.plot`
  
          By default a log-log scaling of the axes is used, if you want to change
          the y axis scaling to linear you can use::
  
              from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel
              from astropy import units as u
  
              pwl = ExpCutoffPowerLawSpectralModel()
              ax = pwl.plot(energy_range=(0.1, 100) * u.TeV)
              ax.set_yscale('linear')
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          # evaluate model
          flux = self(energy).to(flux_unit)
  
          y = self._plot_scale_flux(energy, flux, energy_power)
  
          ax.plot(energy.value, y.value, **kwargs)
  
          self._plot_format_ax(ax, energy, y, energy_power)
          return ax
  
      def plot_error(
          self,
          energy_range,
          ax=None,
          energy_unit="TeV",
          flux_unit="cm-2 s-1 TeV-1",
          energy_power=0,
          n_points=100,
          **kwargs,
      ):
          """Plot spectral model error band.
  
          .. note::
  
              This method calls ``ax.set_yscale("log", nonposy='clip')`` and
              ``ax.set_xscale("log", nonposx='clip')`` to create a log-log representation.
              The additional argument ``nonposx='clip'`` avoids artefacts in the plot,
              when the error band extends to negative values (see also
              https://github.com/matplotlib/matplotlib/issues/8623).
  
              When you call ``plt.loglog()`` or ``plt.semilogy()`` explicitely in your
              plotting code and the error band extends to negative values, it is not
              shown correctly. To circumvent this issue also use
              ``plt.loglog(nonposx='clip', nonposy='clip')``
              or ``plt.semilogy(nonposy='clip')``.
  
          Parameters
          ----------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          energy_range : `~astropy.units.Quantity`
              Plot range
          energy_unit : str, `~astropy.units.Unit`, optional
              Unit of the energy axis
          flux_unit : str, `~astropy.units.Unit`, optional
              Unit of the flux axis
          energy_power : int, optional
              Power of energy to multiply flux axis with
          n_points : int, optional
              Number of evaluation nodes
          **kwargs : dict
              Keyword arguments forwarded to `matplotlib.pyplot.fill_between`
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              Axis
          """
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
  
          kwargs.setdefault("facecolor", "black")
          kwargs.setdefault("alpha", 0.2)
          kwargs.setdefault("linewidth", 0)
  
          emin, emax = energy_range
          energy = MapAxis.from_energy_bounds(emin, emax, n_points, energy_unit).edges
  
          flux, flux_err = self.evaluate_error(energy).to(flux_unit)
  
          y_lo = self._plot_scale_flux(energy, flux - flux_err, energy_power)
          y_hi = self._plot_scale_flux(energy, flux + flux_err, energy_power)
  
          where = (energy >= energy_range[0]) & (energy <= energy_range[1])
          ax.fill_between(energy.value, y_lo.value, y_hi.value, where=where, **kwargs)
  
          self._plot_format_ax(ax, energy, y_lo, energy_power)
          return ax
  
      @staticmethod
      def _plot_format_ax(ax, energy, y, energy_power):
          ax.set_xlabel(f"Energy [{energy.unit}]")
          if energy_power > 0:
              ax.set_ylabel(f"E{energy_power} * Flux [{y.unit}]")
          else:
              ax.set_ylabel(f"Flux [{y.unit}]")
  
          ax.set_xscale("log", nonposx="clip")
          ax.set_yscale("log", nonposy="clip")
  
      @staticmethod
      def _plot_scale_flux(energy, flux, energy_power):
          try:
              eunit = [_ for _ in flux.unit.bases if _.physical_type == "energy"][0]
          except IndexError:
              eunit = energy.unit
          y = flux * np.power(energy, energy_power)
          return y.to(flux.unit * eunit ** energy_power)
  
      def spectral_index(self, energy, epsilon=1e-5):
          """Compute spectral index at given energy.
  
          Parameters
          ----------
          energy : `~astropy.units.Quantity`
              Energy at which to estimate the index
          epsilon : float
              Fractional energy increment to use for determining the spectral index.
  
          Returns
          -------
          index : float
              Estimated spectral index.
          """
          f1 = self(energy)
          f2 = self(energy * (1 + epsilon))
          return np.log(f1 / f2) / np.log(1 + epsilon)
  
      def inverse
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/spectral.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Time-dependent models."""
  import numpy as np
  import scipy.interpolate
  from astropy import units as u
  from astropy.table import Table
  from astropy.time import Time
  from astropy.utils import lazyproperty
  from gammapy.modeling import Parameter
  from gammapy.utils.random import InverseCDFSampler, get_random_state
  from gammapy.utils.scripts import make_path
  from gammapy.utils.time import time_ref_from_dict
  from .core import Model
  
  
  # TODO: make this a small ABC to define a uniform interface.
  class TemporalModel(Model):
      """Temporal model base class.
      evaluates on  astropy.time.Time objects"""
  
      _type = "temporal"
  
      def __call__(self, time):
          """Evaluate model
  
          Parameters
          ----------
          time : `~astropy.time.Time`
              Time object
          """
          kwargs = {par.name: par.quantity for par in self.parameters}
          time = u.Quantity(time.mjd, "day")
          return self.evaluate(time, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def time_sum(t_min, t_max):
          """
          Total time between t_min and t_max
  
          Parameters
          ----------
          t_min, t_max: `~astropy.time.Time`
              Lower and upper bound of integration range
  
          Returns
          -------
          time_sum : `~astropy.time.TimeDelta`
              Summed time in the intervals.
  
          """
          diff = t_max - t_min
          # TODO: this is a work-around for https://github.com/astropy/astropy/issues/10501
          return u.Quantity(np.sum(diff.to_value("day")), "day")
  
      def plot(self, time_range, ax=None):
          """
          Plot Temporal Model.
  
          Parameters
          ----------
          time_range : `~astropy.time.Time`
              times to plot the model
          ax : `~matplotlib.axes.Axes`, optional
              axis
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              axis
          """
  
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          t_min, t_max = time_range
          n_value = 100
          delta = t_max - t_min
          times = t_min + delta * np.linspace(0, 1, n_value)
          val = self(times)
          ax.plot(times.mjd, val)
          return ax
  
  
  class ConstantTemporalModel(TemporalModel):
      """Constant temporal model."""
  
      tag = ["ConstantTemporalModel", "const"]
  
      @staticmethod
      def evaluate(time):
          """Evaluate at given times."""
          return np.ones(time.shape)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : `~astropy.units.Quantity`
              Integrated flux norm on the given time intervals
          """
          return (t_max - t_min) / self.time_sum(t_min, t_max)
  
      @staticmethod
      def sample_time(n_events, t_min, t_max, random_state=0):
          """Sample arrival times of events.
  
          Parameters
          ----------
          n_events : int
              Number of events to sample.
          t_min : `~astropy.time.Time`
              Start time of the sampling.
          t_max : `~astropy.time.Time`
              Stop time of the sampling.
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
              Defines random number generator initialisation.
              Passed to `~gammapy.utils.random.get_random_state`.
  
          Returns
          -------
          time : `~astropy.units.Quantity`
              Array with times of the sampled events.
          """
          random_state = get_random_state(random_state)
  
          t_min = Time(t_min)
          t_max = Time(t_max)
  
          t_stop = (t_max - t_min).sec
  
          time_delta = random_state.uniform(high=t_stop, size=n_events) * u.s
  
          return t_min + time_delta
  
  
  class ExpDecayTemporalModel(TemporalModel):
      r"""Temporal model with an exponential decay.
  
      .. math::
              F(t) = exp(t - t_ref)/t0
  
      Parameters
      ----------
      t0 : `~astropy.units.Quantity`
          Decay time scale
      t_ref: `~astropy.units.Quantity`
          The reference time in mjd
      """
  
      tag = ["ExpDecayTemporalModel", "exp-decay"]
  
      t0 = Parameter("t0", "1 d", frozen=False)
      _t_ref_default = Time("2000-01-01")
      t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=True)
  
      @staticmethod
      def evaluate(time, t0, t_ref):
          """Evaluate at given times"""
          return np.exp(-(time - t_ref) / t0)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : float
              Integrated flux norm on the given time intervals
          """
          pars = self.parameters
          t0 = pars["t0"].quantity
          t_ref = Time(pars["t_ref"].quantity, format="mjd")
          value = self.evaluate(t_max, t0, t_ref) - self.evaluate(t_min, t0, t_ref)
          return -t0 * value / self.time_sum(t_min, t_max)
  
  
  class GaussianTemporalModel(TemporalModel):
      r"""A Gaussian temporal profile
  
      ..math::
              F(t) = exp( -0.5 * \frac{ (t - t_{ref})^2 } { \sigma^2 })
  
      Parameters
      ----------
      t_ref: `~astropy.units.Quantity`
          The reference time in mjd at the peak.
      sigma : `~astropy.units.Quantity`
          Width of the gaussian profile.
      """
  
      tag = ["GaussianTemporalModel", "gauss"]
  
      _t_ref_default = Time("2000-01-01")
      t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=False)
      sigma = Parameter("sigma", "1 d", frozen=False)
  
      @staticmethod
      def evaluate(time, t_ref, sigma):
          return np.exp(-((time - t_ref) ** 2) / (2 * sigma ** 2))
  
      def integral(self, t_min, t_max, **kwargs):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : float
              Integrated flux norm on the given time intervals
          """
          pars = self.parameters
          sigma = pars["sigma"].quantity
          t_ref = Time(pars["t_ref"].quantity, format="mjd")
          norm = np.sqrt(np.pi / 2) * sigma
  
          u_min = (t_min - t_ref) / (np.sqrt(2) * sigma)
          u_max = (t_max - t_ref) / (np.sqrt(2) * sigma)
  
          integral = norm * (scipy.special.erf(u_max) - scipy.special.erf(u_min))
          return integral / self.time_sum(t_min, t_max)
  
  
  class LightCurveTemplateTemporalModel(TemporalModel):
      """Temporal light curve model.
  
      The lightcurve is given as a table with columns ``time`` and ``norm``.
  
      The ``norm`` is supposed to be a unit-less multiplicative factor in the model,
      to be multiplied with a spectral model.
  
      The model does linear interpolation for times between the given ``(time, norm)`` values.
  
      The implementation currently uses `scipy.interpolate. InterpolatedUnivariateSpline`,
      using degree ``k=1`` to get linear interpolation.
      This class also contains an ``integral`` method, making the computation of
      mean fluxes for a given time interval a one-liner.
  
      Parameters
      ----------
      table : `~astropy.table.Table`
          A table with 'TIME' vs 'NORM'
  
      Examples
      --------
      Read an example light curve object:
  
      >>> from gammapy.modeling.models import LightCurveTemplateTemporalModel
      >>> path = '$GAMMAPY_DATA/tests/models/light_curve/lightcrv_PKSB1222+216.fits'
      >>> light_curve = LightCurveTemplateTemporalModel.read(path)
  
      Show basic information about the lightcurve:
  
      >>> print(light_curve)
      LightCurve model summary:
      Start time: 59000.5 MJD
      End time: 61862.5 MJD
      Norm min: 0.01551196351647377
      Norm max: 1.0
  
      Compute ``norm`` at a given time:
  
      >>> light_curve.evaluate(46300)
      0.49059393580053845
  
      Compute mean ``norm`` in a given time interval:
  
      >>> light_curve.mean_norm_in_time_interval(46300, 46301)
      """
  
      tag = ["LightCurveTemplateTemporalModel", "template"]
  
      def __init__(self, table, filename=None):
          self.table = table
          if filename is not None:
              filename = str(make_path(filename))
          self.filename = filename
          super().__init__()
  
      def __str__(self):
          norm = self.table["NORM"]
          return (
              f"{self.__class__.__name__} model summary:\n"
              f"Start time: {self._time[0].mjd} MJD\n"
              f"End time: {self._time[-1].mjd} MJD\n"
              f"Norm min: {norm.min()}\n"
              f"Norm max: {norm.max()}\n"
          )
  
      @classmethod
      def read(cls, path):
          """Read lightcurve model table from FITS file.
  
          TODO: This doesn't read the XML part of the model yet.
          """
          filename = str(make_path(path))
          return cls(Table.read(filename), filename=filename)
  
      def write(self, path=None, overwrite=False):
          if path is None:
              path = self.filename
          if path is None:
              raise ValueError(f"filename is required for {self.tag}")
          else:
              self.filename = str(make_path(path))
              self.table.write(self.filename, overwrite=overwrite)
  
      @lazyproperty
      def _interpolator(self, ext=0):
          x = self._time.value
          y = self.table["NORM"].data
          return scipy.interpolate.InterpolatedUnivariateSpline(x, y, k=1, ext=ext)
  
      @lazyproperty
      def _time_ref(self):
          return time_ref_from_dict(self.table.meta)
  
      @lazyproperty
      def _time(self):
          return self._time_ref + self.table["TIME"].data * getattr(
              u, self.table.meta["TIMEUNIT"]
          )
  
      def evaluate(self, time, ext=0):
          """Evaluate for a given time.
  
          Parameters
          ----------
          time : array_like
              Time since the ``reference`` time.
          ext : int or str, optional, default: 0
              Parameter passed to ~scipy.interpolate.InterpolatedUnivariateSpline
              Controls the extrapolation mode for GTIs outside the range
              0 or "extrapolate", return the extrapolated value.
              1 or "zeros", return 0
              2 or "raise", raise a ValueError
              3 or "const", return the boundary value.
  
  
          Returns
          -------
          norm : array_like
              Norm at the given times.
          """
          return self._interpolator(time, ext=ext)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
          Returns
          -------
          norm: The model integrated flux
          """
  
          n1 = self._interpolator.antiderivative()(t_max.mjd)
          n2 = self._interpolator.antiderivative()(t_min.mjd)
          return u.Quantity(n1 - n2, "day") / self.time_sum(t_min, t_max)
  
      def sample_time(self, n_events, t_min, t_max, t_delta="1 s", random_state=0):
          """Sample arrival times of events.
  
          Parameters
          ----------
          n_events : int
              Number of events to sample.
          t_min : `~astropy.time.Time`
              Start time of the sampling.
          t_max : `~astropy.time.Time`
              Stop time of the sampling.
          t_delta : `~astropy.units.Quantity`
              Time step used for sampling of the temporal model.
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
              Defines random number generator initialisation.
              Passed to `~gammapy.utils.random.get_random_state`.
  
          Returns
          -------
          time : `~astropy.units.Quantity`
              Array with times of the sampled events.
          """
          time_unit = getattr(u, self.table.meta["TIMEUNIT"])
  
          t_min = Time(t_min)
          t_max = Time(t_max)
          t_delta = u.Quantity(t_delta)
          random_state = get_random_state(random_state)
  
          ontime = u.Quantity((t_max - t_min).sec, "s")
          t_stop = ontime.to_value(time_unit)
  
          # TODO: the separate time unit handling is unfortunate, but the quantity support for np.arange and np.interp
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/temporal.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """Time-dependent models."""
  import numpy as np
  import scipy.interpolate
  from astropy import units as u
  from astropy.table import Table
  from astropy.time import Time
  from astropy.utils import lazyproperty
  from gammapy.modeling import Parameter
  from gammapy.utils.random import InverseCDFSampler, get_random_state
  from gammapy.utils.scripts import make_path
  from gammapy.utils.time import time_ref_from_dict
  from .core import Model
  
  
  # TODO: make this a small ABC to define a uniform interface.
  class TemporalModel(Model):
      """Temporal model base class.
      evaluates on  astropy.time.Time objects"""
      _type = "temporal"
  
      def __call__(self, time):
          """Evaluate model
  
          Parameters
          ----------
          time : `~astropy.time.Time`
              Time object
          """
          kwargs = {par.name: par.quantity for par in self.parameters}
          time = u.Quantity(time.mjd, "day")
          return self.evaluate(time, **kwargs)
  
      @property
      def type(self):
          return self._type
  
      @staticmethod
      def time_sum(t_min, t_max):
          """
          Total time between t_min and t_max
  
          Parameters
          ----------
          t_min, t_max: `~astropy.time.Time`
              Lower and upper bound of integration range
  
          Returns
          -------
          time_sum : `~astropy.time.TimeDelta`
              Summed time in the intervals.
  
          """
          diff = t_max - t_min
          # TODO: this is a work-around for https://github.com/astropy/astropy/issues/10501
          return u.Quantity(np.sum(diff.to_value("day")), "day")
  
      def plot(self, time_range, ax=None):
          """
          Plot Temporal Model.
  
          Parameters
          ----------
          time_range : `~astropy.time.Time`
              times to plot the model
          ax : `~matplotlib.axes.Axes`, optional
              axis
  
          Returns
          -------
          ax : `~matplotlib.axes.Axes`, optional
              axis
          """
  
  
          import matplotlib.pyplot as plt
  
          ax = plt.gca() if ax is None else ax
          t_min, t_max = time_range
          n_value = 100
          delta = (t_max - t_min)
          times = t_min + delta * np.linspace(0, 1, n_value)
          val = self(times)
          ax.plot(times.mjd, val)
          return ax
  
  class ConstantTemporalModel(TemporalModel):
      """Constant temporal model."""
  
      tag = ["ConstantTemporalModel", "const"]
  
      @staticmethod
      def evaluate(time):
          """Evaluate at given times."""
          return np.ones(time.shape)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : `~astropy.units.Quantity`
              Integrated flux norm on the given time intervals
          """
          return (t_max - t_min) / self.time_sum(t_min, t_max)
  
      @staticmethod
      def sample_time(n_events, t_min, t_max, random_state=0):
          """Sample arrival times of events.
  
          Parameters
          ----------
          n_events : int
              Number of events to sample.
          t_min : `~astropy.time.Time`
              Start time of the sampling.
          t_max : `~astropy.time.Time`
              Stop time of the sampling.
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
              Defines random number generator initialisation.
              Passed to `~gammapy.utils.random.get_random_state`.
  
          Returns
          -------
          time : `~astropy.units.Quantity`
              Array with times of the sampled events.
          """
          random_state = get_random_state(random_state)
  
          t_min = Time(t_min)
          t_max = Time(t_max)
  
          t_stop = (t_max - t_min).sec
  
          time_delta = random_state.uniform(high=t_stop, size=n_events) * u.s
  
          return t_min + time_delta
  
  
  class ExpDecayTemporalModel(TemporalModel):
      r"""Temporal model with an exponential decay.
  
      .. math::
              F(t) = exp(t - t_ref)/t0
  
      Parameters
      ----------
      t0 : `~astropy.units.Quantity`
          Decay time scale
      t_ref: `~astropy.units.Quantity`
          The reference time in mjd
      """
  
      tag = ["ExpDecayTemporalModel", "exp-decay"]
  
      t0 = Parameter("t0", "1 d", frozen=False)
      _t_ref_default = Time("2000-01-01")
      t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=True)
  
      @staticmethod
      def evaluate(time, t0, t_ref):
          """Evaluate at given times"""
          return np.exp(-(time - t_ref) / t0)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : float
              Integrated flux norm on the given time intervals
          """
          pars = self.parameters
          t0 = pars["t0"].quantity
          t_ref = Time(pars["t_ref"].quantity, format="mjd")
          value = self.evaluate(t_max, t0, t_ref) - self.evaluate(t_min, t0, t_ref)
          return -t0 * value / self.time_sum(t_min, t_max)
  
  
  
  
  class GaussianTemporalModel(TemporalModel):
      r"""A Gaussian temporal profile
  
      ..math::
              F(t) = exp( -0.5 * \frac{ (t - t_{ref})^2 } { \sigma^2 })
  
      Parameters
      ----------
      t_ref: `~astropy.units.Quantity`
          The reference time in mjd at the peak.
      sigma : `~astropy.units.Quantity`
          Width of the gaussian profile.
      """
  
      tag = ["GaussianTemporalModel", "gauss"]
  
      _t_ref_default = Time("2000-01-01")
      t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=False)
      sigma = Parameter("sigma", "1 d", frozen=False)
  
      @staticmethod
      def evaluate(time, t_ref, sigma):
          return np.exp(-((time - t_ref) ** 2) / (2 * sigma ** 2))
  
      def integral(self, t_min, t_max, **kwargs):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
  
          Returns
          -------
          norm : float
              Integrated flux norm on the given time intervals
          """
          pars = self.parameters
          sigma = pars["sigma"].quantity
          t_ref = Time(pars["t_ref"].quantity, format="mjd")
          norm = np.sqrt(np.pi / 2) * sigma
  
          u_min = (t_min - t_ref) / (np.sqrt(2) * sigma)
          u_max = (t_max - t_ref) / (np.sqrt(2) * sigma)
  
          integral = norm * (scipy.special.erf(u_max) - scipy.special.erf(u_min))
          return integral / self.time_sum(t_min, t_max)
  
  
  class LightCurveTemplateTemporalModel(TemporalModel):
      """Temporal light curve model.
  
      The lightcurve is given as a table with columns ``time`` and ``norm``.
  
      The ``norm`` is supposed to be a unit-less multiplicative factor in the model,
      to be multiplied with a spectral model.
  
      The model does linear interpolation for times between the given ``(time, norm)`` values.
  
      The implementation currently uses `scipy.interpolate. InterpolatedUnivariateSpline`,
      using degree ``k=1`` to get linear interpolation.
      This class also contains an ``integral`` method, making the computation of
      mean fluxes for a given time interval a one-liner.
  
      Parameters
      ----------
      table : `~astropy.table.Table`
          A table with 'TIME' vs 'NORM'
  
      Examples
      --------
      Read an example light curve object:
  
      >>> from gammapy.modeling.models import LightCurveTemplateTemporalModel
      >>> path = '$GAMMAPY_DATA/tests/models/light_curve/lightcrv_PKSB1222+216.fits'
      >>> light_curve = LightCurveTemplateTemporalModel.read(path)
  
      Show basic information about the lightcurve:
  
      >>> print(light_curve)
      LightCurve model summary:
      Start time: 59000.5 MJD
      End time: 61862.5 MJD
      Norm min: 0.01551196351647377
      Norm max: 1.0
  
      Compute ``norm`` at a given time:
  
      >>> light_curve.evaluate(46300)
      0.49059393580053845
  
      Compute mean ``norm`` in a given time interval:
  
      >>> light_curve.mean_norm_in_time_interval(46300, 46301)
      """
  
      tag = ["LightCurveTemplateTemporalModel", "template"]
  
      def __init__(self, table, filename=None):
          self.table = table
          if filename is not None:
              filename = str(make_path(filename))
          self.filename = filename
          super().__init__()
  
      def __str__(self):
          norm = self.table["NORM"]
          return (
              f"{self.__class__.__name__} model summary:\n"
              f"Start time: {self._time[0].mjd} MJD\n"
              f"End time: {self._time[-1].mjd} MJD\n"
              f"Norm min: {norm.min()}\n"
              f"Norm max: {norm.max()}\n"
          )
  
      @classmethod
      def read(cls, path):
          """Read lightcurve model table from FITS file.
  
          TODO: This doesn't read the XML part of the model yet.
          """
          filename = str(make_path(path))
          return cls(Table.read(filename), filename=filename)
  
      def write(self, path=None, overwrite=False):
          if path is None:
              path = self.filename
          if path is None:
              raise ValueError(f"filename is required for {self.tag}")
          else:
              self.filename = str(make_path(path))
              self.table.write(self.filename, overwrite=overwrite)
  
      @lazyproperty
      def _interpolator(self, ext=0):
          x = self._time.value
          y = self.table["NORM"].data
          return scipy.interpolate.InterpolatedUnivariateSpline(x, y, k=1, ext=ext)
  
      @lazyproperty
      def _time_ref(self):
          return time_ref_from_dict(self.table.meta)
  
      @lazyproperty
      def _time(self):
          return self._time_ref + self.table["TIME"].data * getattr(
              u, self.table.meta["TIMEUNIT"]
          )
  
      def evaluate(self, time, ext=0):
          """Evaluate for a given time.
  
          Parameters
          ----------
          time : array_like
              Time since the ``reference`` time.
          ext : int or str, optional, default: 0
              Parameter passed to ~scipy.interpolate.InterpolatedUnivariateSpline
              Controls the extrapolation mode for GTIs outside the range
              0 or "extrapolate", return the extrapolated value.
              1 or "zeros", return 0
              2 or "raise", raise a ValueError
              3 or "const", return the boundary value.
  
  
          Returns
          -------
          norm : array_like
              Norm at the given times.
          """
          return self._interpolator(time, ext=ext)
  
      def integral(self, t_min, t_max):
          """Evaluate the integrated flux within the given time intervals
  
          Parameters
          ----------
          t_min: `~astropy.time.Time`
              Start times of observation
          t_max: `~astropy.time.Time`
              Stop times of observation
          Returns
          -------
          norm: The model integrated flux
          """
  
          n1 = self._interpolator.antiderivative()(t_max.mjd)
          n2 = self._interpolator.antiderivative()(t_min.mjd)
          return u.Quantity(n1 - n2, "day") / self.time_sum(t_min, t_max)
  
      def sample_time(self, n_events, t_min, t_max, t_delta="1 s", random_state=0):
          """Sample arrival times of events.
  
          Parameters
          ----------
          n_events : int
              Number of events to sample.
          t_min : `~astropy.time.Time`
              Start time of the sampling.
          t_max : `~astropy.time.Time`
              Stop time of the sampling.
          t_delta : `~astropy.units.Quantity`
              Time step used for sampling of the temporal model.
          random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
              Defines random number generator initialisation.
              Passed to `~gammapy.utils.random.get_random_state`.
  
          Returns
          -------
          time : `~astropy.units.Quantity`
              Array with times of the sampled events.
          """
          time_unit = getattr(u, self.table.meta["TIMEUNIT"])
  
          t_min = Time(t_min)
          t_max = Time(t_max)
          t_delta = u.Quantity(t_delta)
          random_state = get_random_state(random_state)
  
          ontime = u.Quantity((t_max - t_min).sec, "s")
          t_stop = ontime.to_value(time_unit)
  
          # TODO: the separate time unit handling is unfortunate, but the quantity support for np.arange and
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/temporal.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.modeling.models import Model, Parameter, Parameters
  from gammapy.datasets import Datasets
  from gammapy.utils.testing import requires_data
  
  
  class MyModel(Model):
      """Simple model example"""
  
      x = Parameter("x", 1, "cm")
      y = Parameter("y", 2)
  
  
  class CoModel(Model):
      """Compound model example"""
  
      norm = Parameter("norm", 42, "cm")
  
      def __init__(self, m1, m2, norm=norm.quantity):
          self.m1 = m1
          self.m2 = m2
          super().__init__(norm=norm)
  
      @property
      def parameters(self):
          return Parameters([self.norm]) + self.m1.parameters + self.m2.parameters
  
  
  class WrapperModel(Model):
      """Wrapper compound model.
  
      Dynamically generated parameters in `__init__`,
      and a parameter name conflict with the wrapped
      model, both have a parameter called "y".
      """
  
      def __init__(self, m1, a=1, y=99):
          self.m1 = m1
          a = Parameter("a", a)
          y = Parameter("y", y)
          self.default_parameters = Parameters([a, y])
          super().__init__(a=a, y=y)
  
      @property
      def parameters(self):
          return Parameters([self.a, self.y]) + self.m1.parameters
  
  
  def test_model_class():
      assert isinstance(MyModel.parameters, property)
      assert MyModel.x.name == "x"
      assert MyModel.default_parameters["x"] is MyModel.x
  
  
  def test_model_class_par_init():
      x = Parameter("x", 4, "cm")
      y = Parameter("y", 10)
  
      model = MyModel(x=x, y=y)
  
      assert x is model.x
      assert y is model.y
  
  
  def test_model_init():
      m = MyModel()
      assert m.x.name == "x"
      assert m.x.value == 1
      assert m.x is m.parameters[0]
      assert m.y is m.parameters[1]
      assert m.parameters is not MyModel.default_parameters
  
      m = MyModel(x="99 cm")
      assert m.x.value == 99
      assert m.y.value == 2
  
      # Currently we always convert to the default unit of a parameter
      # TODO: discuss if this is the behaviour we want, or if we instead
      # should change to the user-set unit, as long as it's compatible
      m = MyModel(x=99 * u.m)
      assert_allclose(m.x.value, 9900)
      assert m.x.unit == "cm"
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99)
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99 * u.s)
  
  
  def test_wrapper_model():
      outer = MyModel()
      m = WrapperModel(outer)
  
      assert isinstance(m.a, Parameter)
      assert m.y.value == 99
  
      assert m.parameters.names == ["a", "y", "x", "y"]
  
  
  def test_model_parameter():
      m = MyModel(x="99 cm")
      assert isinstance(m.x, Parameter)
      assert m.x.value == 99
      assert m.x.unit == "cm"
  
      with pytest.raises(TypeError):
          m.x = 99
  
      with pytest.raises(TypeError):
          m.x = 99 * u.cm
  
  
  # TODO: implement parameter linking. Not working ATM!
  @pytest.mark.xfail()
  def test_model_parameter_link():
      # Assigning a parameter should create a link
      m = MyModel()
      par = MyModel.x.copy()
      m.x = par
      assert isinstance(m.x, Parameter)
      assert m.x is par
      # model.parameters should be in sync with attributes
      assert m.x is m.parameters["x"]
  
  
  def test_model_copy():
      m = MyModel()
  
      m2 = m.copy()
  
      # Models should be independent
      assert m.parameters is not m2.parameters
      assert m.parameters[0] is not m2.parameters[0]
  
  
  def test_model_create():
      spectral_model = Model.create(
          "pl-2", model_type="spectral", amplitude="1e-10 cm-2 s-1", index=3
      )
      assert "PowerLaw2SpectralModel" in spectral_model.tag
      assert_allclose(spectral_model.index.value, 3)
  
  
  def test_compound_model():
      m1 = MyModel()
      m2 = MyModel(x=10 * u.cm, y=20)
      m = CoModel(m1, m2)
      assert len(m.parameters) == 5
      assert m.parameters.names == ["norm", "x", "y", "x", "y"]
      assert_allclose(m.parameters.values, [42, 1, 2, 10, 20])
  
  
  def test_parameter_link_init():
      m1 = MyModel()
      m2 = MyModel(y=m1.y)
  
      assert m1.y is m2.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  def test_parameter_link():
      m1 = MyModel()
      m2 = MyModel()
  
      m2.y = m1.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  @requires_data()
  def test_models_management(tmp_path):
      path = "$GAMMAPY_DATA/tests/models"
      filedata = "gc_example_datasets.yaml"
      filemodel = "gc_example_models.yaml"
  
      datasets = Datasets.read(path, filedata, filemodel)
  
      model1 = datasets.models[0].copy(name="model1", datasets_names=None)
      model2 = datasets.models[0].copy(name="model2", datasets_names=[datasets[1].name])
      model3 = datasets.models[0].copy(name="model3", datasets_names=[datasets[0].name])
  
      model1b = datasets.models[0].copy(name="model1", datasets_names=None)
      model1b.spectral_model.amplitude.value *= 2
  
      names0 = datasets[0].models.names
      names1 = datasets[1].models.names
  
      datasets[0].models.append(model1)
      _ = datasets[0].models + model2
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[0].models["model1"].datasets_names is None
      assert datasets[0].models["model2"].datasets_names == [
          datasets[1].name,
          datasets[0].name,
      ]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      # TODO consistency check at datasets level ?
      # or force same Models for each dataset._models on datasets init ?
      # here we have the right behavior: model1 and model2 are also added to dataset1
      # because serialization create a global model object shared by all datasets
      # if that was not the case we could have inconsistancies
      # such as model1.datasets_names == None added only to dataset1
      # user can still create such inconsistancies if they define datasets
      # with diferent Models objects for each dataset.
  
      del datasets[0].models["model1"]
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.remove(model2)
      assert datasets[0].models.names == names0
  
      datasets.models.append(model2)
      assert model2 in datasets.models
      assert model2 in datasets[1].models
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.extend([model1, model3])
      assert datasets[0].models.names == names0 + ["model2", "model1", "model3"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      assert datasets[0].models.names == names0
      assert datasets[1].models.names == names1
      datasets.models.extend([model1, model2, model3])
      assert datasets[0].models.names == names0 + ["model1", "model2", "model3"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      _ = datasets.models + [model1, model2]
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      datasets[0].models["model2"] = model3
      assert datasets[0].models.names == names0 + ["model1", "model3"]
      assert datasets[1].models.names == names1 + ["model1"]
  
      datasets.models.remove(model1)
      datasets[0].models = model1
      _ = datasets.models  # auto-update models
  
      npred1 = datasets[0].npred().data.sum()
      datasets.models.remove(model1)
      npred0 = datasets[0].npred().data.sum()
      datasets.models.append(model1b)
      npred1b = datasets[0].npred().data.sum()
      assert npred1b != npred1
      assert npred1b != npred0
      assert_allclose(npred1b, 5157.137554, rtol=1e-5)
  
      datasets.models.remove(model1b)
      _ = datasets.models  # auto-update models
      newmodels = [datasets.models[0].copy() for k in range(48)]
      datasets.models.extend(newmodels)
  
      datasets[0].use_cache = False
      nocache = datasets[0].npred().data.sum()
      datasets[0].use_cache = True
      assert_allclose(datasets[0].npred().data.sum(), nocache)
  
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/tests/test_core.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.modeling.models import Model, Parameter, Parameters
  from gammapy.datasets import Datasets
  from gammapy.utils.testing import requires_data
  
  
  class MyModel(Model):
      """Simple model example"""
      x = Parameter("x", 1, "cm")
      y = Parameter("y", 2)
  
  
  class CoModel(Model):
      """Compound model example"""
      norm = Parameter("norm", 42, "cm")
  
      def __init__(self, m1, m2, norm=norm.quantity):
          self.m1 = m1
          self.m2 = m2
          super().__init__(norm=norm)
  
      @property
      def parameters(self):
          return Parameters([self.norm]) + self.m1.parameters + self.m2.parameters
  
  
  class WrapperModel(Model):
      """Wrapper compound model.
  
      Dynamically generated parameters in `__init__`,
      and a parameter name conflict with the wrapped
      model, both have a parameter called "y".
      """
  
      def __init__(self, m1, a=1, y=99):
          self.m1 = m1
          a = Parameter("a", a)
          y = Parameter("y", y)
          self.default_parameters = Parameters([a, y])
          super().__init__(a=a, y=y)
  
      @property
      def parameters(self):
          return Parameters([self.a, self.y]) + self.m1.parameters
  
  
  def test_model_class():
      assert isinstance(MyModel.parameters, property)
      assert MyModel.x.name == "x"
      assert MyModel.default_parameters["x"] is MyModel.x
  
  
  def test_model_class_par_init():
      x = Parameter("x", 4, "cm")
      y = Parameter("y", 10)
  
      model = MyModel(x=x, y=y)
  
      assert x is model.x
      assert y is model.y
  
  
  def test_model_init():
      m = MyModel()
      assert m.x.name == "x"
      assert m.x.value == 1
      assert m.x is m.parameters[0]
      assert m.y is m.parameters[1]
      assert m.parameters is not MyModel.default_parameters
  
      m = MyModel(x="99 cm")
      assert m.x.value == 99
      assert m.y.value == 2
  
      # Currently we always convert to the default unit of a parameter
      # TODO: discuss if this is the behaviour we want, or if we instead
      # should change to the user-set unit, as long as it's compatible
      m = MyModel(x=99 * u.m)
      assert_allclose(m.x.value, 9900)
      assert m.x.unit == "cm"
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99)
  
      with pytest.raises(u.UnitConversionError):
          MyModel(x=99 * u.s)
  
  
  def test_wrapper_model():
      outer = MyModel()
      m = WrapperModel(outer)
  
      assert isinstance(m.a, Parameter)
      assert m.y.value == 99
  
      assert m.parameters.names == ["a", "y", "x", "y"]
  
  
  def test_model_parameter():
      m = MyModel(x="99 cm")
      assert isinstance(m.x, Parameter)
      assert m.x.value == 99
      assert m.x.unit == "cm"
  
      with pytest.raises(TypeError):
          m.x = 99
  
      with pytest.raises(TypeError):
          m.x = 99 * u.cm
  
  
  # TODO: implement parameter linking. Not working ATM!
  @pytest.mark.xfail()
  def test_model_parameter_link():
      # Assigning a parameter should create a link
      m = MyModel()
      par = MyModel.x.copy()
      m.x = par
      assert isinstance(m.x, Parameter)
      assert m.x is par
      # model.parameters should be in sync with attributes
      assert m.x is m.parameters["x"]
  
  
  def test_model_copy():
      m = MyModel()
  
      m2 = m.copy()
  
      # Models should be independent
      assert m.parameters is not m2.parameters
      assert m.parameters[0] is not m2.parameters[0]
  
  
  def test_model_create():
      spectral_model = Model.create(
          "pl-2", model_type="spectral", amplitude="1e-10 cm-2 s-1", index=3
      )
      assert "PowerLaw2SpectralModel" in spectral_model.tag
      assert_allclose(spectral_model.index.value, 3)
  
  
  def test_compound_model():
      m1 = MyModel()
      m2 = MyModel(x=10 * u.cm, y=20)
      m = CoModel(m1, m2)
      assert len(m.parameters) == 5
      assert m.parameters.names == ["norm", "x", "y", "x", "y"]
      assert_allclose(m.parameters.values, [42, 1, 2, 10, 20])
  
  
  def test_parameter_link_init():
      m1 = MyModel()
      m2 = MyModel(y=m1.y)
  
      assert m1.y is m2.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  def test_parameter_link():
      m1 = MyModel()
      m2 = MyModel()
  
      m2.y = m1.y
  
      m1.y.value = 100
      assert_allclose(m2.y.value, 100)
  
  
  @requires_data()
  def test_models_management(tmp_path):
      path = "$GAMMAPY_DATA/tests/models"
      filedata = "gc_example_datasets.yaml"
      filemodel = "gc_example_models.yaml"
  
      datasets = Datasets.read(path, filedata, filemodel)
  
      model1 = datasets.models[0].copy(name="model1", datasets_names=None)
      model2 = datasets.models[0].copy(name="model2", datasets_names=[datasets[1].name])
      model3 = datasets.models[0].copy(name="model3", datasets_names=[datasets[0].name])
  
      model1b = datasets.models[0].copy(name="model1", datasets_names=None)
      model1b.spectral_model.amplitude.value *= 2
  
      names0 = datasets[0].models.names
      names1 = datasets[1].models.names
  
      datasets[0].models.append(model1)
      _ = datasets[0].models + model2
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[0].models["model1"].datasets_names is None
      assert datasets[0].models["model2"].datasets_names == [
          datasets[1].name,
          datasets[0].name,
      ]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      # TODO consistency check at datasets level ?
      # or force same Models for each dataset._models on datasets init ?
      # here we have the right behavior: model1 and model2 are also added to dataset1
      # because serialization create a global model object shared by all datasets
      # if that was not the case we could have inconsistancies
      # such as model1.datasets_names == None added only to dataset1
      # user can still create such inconsistancies if they define datasets
      # with diferent Models objects for each dataset.
  
      del datasets[0].models["model1"]
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.remove(model2)
      assert datasets[0].models.names == names0
  
      datasets.models.append(model2)
      assert model2 in datasets.models
      assert model2 in datasets[1].models
      assert datasets[0].models.names == names0 + ["model2"]
  
      datasets[0].models.extend([model1, model3])
      assert datasets[0].models.names == names0 + ["model2", "model1", "model3"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      assert datasets[0].models.names == names0
      assert datasets[1].models.names == names1
      datasets.models.extend([model1, model2, model3])
      assert datasets[0].models.names == names0 + ["model1", "model2", "model3"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      for m in [model1, model2, model3]:
          datasets.models.remove(m)
      _ = datasets.models + [model1, model2]
      assert datasets[0].models.names == names0 + ["model1", "model2"]
      assert datasets[1].models.names == names1 + ["model1", "model2"]
  
      datasets[0].models["model2"] = model3
      assert datasets[0].models.names == names0 + ["model1", "model3"]
      assert datasets[1].models.names == names1 + ["model1"]
  
      datasets.models.remove(model1)
      datasets[0].models = model1
      _ = datasets.models  # auto-update models
  
      npred1 = datasets[0].npred().data.sum()
      datasets.models.remove(model1)
      npred0 = datasets[0].npred().data.sum()
      datasets.models.append(model1b)
      npred1b = datasets[0].npred().data.sum()
      assert npred1b != npred1
      assert npred1b != npred0
      assert_allclose(npred1b, 5157.137554, rtol=1e-5)
  
      datasets.models.remove(model1b)
      _ = datasets.models  # auto-update models
      newmodels = [datasets.models[0].copy() for k in range(48)]
      datasets.models.extend(newmodels)
  
      datasets[0].use_cache = False
      nocache = datasets[0].npred().data.sum()
      datasets[0].use_cache = True
      assert_allclose(datasets[0].npred().data.sum(), nocache)
  
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/tests/test_core.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.time import Time
  from gammapy.data.gti import GTI
  from gammapy.datasets.map import MapEvaluator
  from gammapy.irf import EDispKernel, PSFKernel
  from gammapy.maps import Map, MapAxis, WcsGeom
  from gammapy.modeling.models import (
      BackgroundModel,
      ConstantSpectralModel,
      ConstantTemporalModel,
      GaussianSpatialModel,
      Models,
      PointSpatialModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      TemplateSpatialModel,
      CompoundSpectralModel,
      SkyModel,
      create_fermi_isotropic_diffuse_model,
  )
  from gammapy.utils.testing import requires_data
  
  
  @pytest.fixture(scope="session")
  def sky_model():
      spatial_model = GaussianSpatialModel(
          lon_0="3 deg", lat_0="4 deg", sigma="3 deg", frame="galactic"
      )
      spectral_model = PowerLawSpectralModel(
          index=2, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      temporal_model = ConstantTemporalModel()
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
          temporal_model=temporal_model,
          name="source-1",
      )
  
  
  @pytest.fixture(scope="session")
  def gti():
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      t_ref = Time(55555, format="mjd")
      gti = GTI.create(start, stop, reference_time=t_ref)
      return gti
  
  
  @pytest.fixture(scope="session")
  def diffuse_model():
      axis = MapAxis.from_nodes([0.1, 100], name="energy_true", unit="TeV", interp="log")
      m = Map.create(
          npix=(4, 3), binsz=2, axes=[axis], unit="cm-2 s-1 MeV-1 sr-1", frame="galactic"
      )
      m.data += 42
      spatial_model = TemplateSpatialModel(m, normalize=False)
      return SkyModel(PowerLawNormSpectralModel(), spatial_model)
  
  
  @pytest.fixture(scope="session")
  def geom():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 3), unit=u.TeV, name="energy")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def geom_true():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 4), unit=u.TeV, name="energy_true")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def exposure(geom_true):
      m = Map.from_geom(geom_true)
      m.quantity = np.ones(geom_true.data_shape) * u.Quantity("100 m2 s")
      m.data[1] *= 10
      return m
  
  
  @pytest.fixture(scope="session")
  def background(geom):
      m = Map.from_geom(geom)
      m.quantity = np.ones(geom.data_shape) * 1e-7
      return m
  
  
  @pytest.fixture(scope="session")
  def edisp(geom, geom_true):
      e_reco = geom.axes["energy"].edges
      e_true = geom_true.axes["energy_true"].edges
      return EDispKernel.from_diagonal_response(e_true=e_true, e_reco=e_reco)
  
  
  @pytest.fixture(scope="session")
  def psf(geom_true):
      sigma = 0.5 * u.deg
      return PSFKernel.from_gauss(geom_true, sigma)
  
  
  @pytest.fixture(scope="session")
  def evaluator(sky_model, exposure, psf, edisp, gti):
      return MapEvaluator(sky_model, exposure, psf=psf, edisp=edisp, gti=gti)
  
  
  @pytest.fixture(scope="session")
  def diffuse_evaluator(diffuse_model, exposure, psf, edisp):
      return MapEvaluator(diffuse_model, exposure, psf=psf, edisp=edisp)
  
  
  @pytest.fixture(scope="session")
  def sky_models(sky_model):
      sky_model_2 = sky_model.copy(name="source-2")
      sky_model_3 = sky_model.copy(name="source-3")
      return Models([sky_model_2, sky_model_3])
  
  
  @pytest.fixture(scope="session")
  def sky_models_2(sky_model):
      sky_model_4 = sky_model.copy(name="source-4")
      sky_model_5 = sky_model.copy(name="source-5")
      return Models([sky_model_4, sky_model_5])
  
  
  def test_sky_model_init():
      with pytest.raises(TypeError):
          spatial_model = GaussianSpatialModel()
          SkyModel(spectral_model=1234, spatial_model=spatial_model)
  
      with pytest.raises(TypeError):
          SkyModel(spectral_model=PowerLawSpectralModel(), spatial_model=1234)
  
  
  def test_sky_model_spatial_none_io(tmpdir):
      pwl = PowerLawSpectralModel()
      model = SkyModel(spectral_model=pwl, name="test")
      models = Models([model])
  
      filename = tmpdir / "test-models-none.yaml"
      models.write(filename)
  
      models = Models.read(filename)
  
      assert models["test"].spatial_model is None
  
  
  def test_sky_model_spatial_none_evaluate(geom_true, gti):
      pwl = PowerLawSpectralModel()
      model = SkyModel(spectral_model=pwl, name="test")
  
      data = model.evaluate_geom(geom_true, gti).to_value("cm-2 s-1 TeV-1")
  
      assert data.shape == (3, 1, 1)
      assert_allclose(data[0], 1.256774e-11, rtol=1e-6)
  
  
  def test_skymodel_addition(sky_model, sky_models, sky_models_2, diffuse_model):
      models = sky_model + sky_model.copy()
      assert isinstance(models, Models)
      assert len(models) == 2
  
      models = sky_model + sky_models
      assert isinstance(models, Models)
      assert len(models) == 3
  
      models = sky_models + sky_model
      assert isinstance(models, Models)
      assert len(models) == 3
  
      models = sky_models + diffuse_model
      assert isinstance(models, Models)
      assert len(models) == 3
  
      models = sky_models + sky_models_2
      assert isinstance(models, Models)
      assert len(models) == 4
  
      models = sky_model + sky_models
      assert isinstance(models, Models)
      assert len(models) == 3
  
  
  def test_background_model(background):
      bkg1 = BackgroundModel(background)
      bkg1.spectral_model.norm.value = 2.0
      npred1 = bkg1.evaluate()
      assert_allclose(npred1.data[0][0][0], background.data[0][0][0] * 2.0, rtol=1e-3)
      assert_allclose(npred1.data.sum(), background.data.sum() * 2.0, rtol=1e-3)
  
      bkg2 = BackgroundModel(background)
      bkg2.spectral_model.norm.value = 2.0
      bkg2.spectral_model.tilt.value = 0.2
      bkg2.spectral_model.reference.quantity = "1000 GeV"
  
      npred2 = bkg2.evaluate()
      assert_allclose(npred2.data[0][0][0], 2.254e-07, rtol=1e-3)
      assert_allclose(npred2.data.sum(), 7.352e-06, rtol=1e-3)
  
  
  def test_background_model_io(tmpdir, background):
      filename = str(tmpdir / "test-bkg-file.fits")
      bkg = BackgroundModel(background, filename=filename)
      bkg.spectral_model.norm.value = 2.0
      bkg.map.write(filename, overwrite=True)
      bkg_dict = bkg.to_dict()
      bkg_read = bkg.from_dict(bkg_dict)
  
      assert_allclose(
          bkg_read.evaluate().data.sum(), background.data.sum() * 2.0, rtol=1e-3
      )
      assert bkg_read.filename == filename
  
  
  class TestSkyModels:
      @staticmethod
      def test_parameters(sky_models):
          parnames = [
              "index",
              "amplitude",
              "reference",
              "lon_0",
              "lat_0",
              "sigma",
              "e",
              "phi",
          ] * 2
          assert sky_models.parameters.names == parnames
  
          # Check that model parameters are references to the parts
          p1 = sky_models.parameters["lon_0"]
          p2 = sky_models[0].parameters["lon_0"]
          assert p1 is p2
  
      @staticmethod
      def test_str(sky_models):
          assert "Component 0" in str(sky_models)
          assert "Component 1" in str(sky_models)
  
      @staticmethod
      def test_get_item(sky_models):
          model = sky_models["source-2"]
          assert model.name == "source-2"
  
          model = sky_models["source-3"]
          assert model.name == "source-3"
  
          with pytest.raises(ValueError):
              sky_models["spam"]
  
      @staticmethod
      def test_names(sky_models):
          assert sky_models.names == ["source-2", "source-3"]
  
  
  @requires_data()
  def test_models_mutation(sky_model, sky_models, sky_models_2):
      mods = sky_models
  
      mods.insert(0, sky_model)
      assert mods.names == ["source-1", "source-2", "source-3"]
  
      mods.extend(sky_models_2)
      assert mods.names == ["source-1", "source-2", "source-3", "source-4", "source-5"]
  
      mod3 = mods[3]
      mods.remove(mods[3])
      assert mods.names == ["source-1", "source-2", "source-3", "source-5"]
      mods.append(mod3)
      assert mods.names == ["source-1", "source-2", "source-3", "source-5", "source-4"]
      mods.pop(3)
      assert mods.names == ["source-1", "source-2", "source-3", "source-4"]
  
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods.append(sky_model)
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods.insert(0, sky_model)
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods.extend(sky_models_2)
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods = sky_models + sky_models_2
  
  
  class TestSkyModel:
      @staticmethod
      def test_repr(sky_model):
          assert "SkyModel" in repr(sky_model)
  
      @staticmethod
      def test_str(sky_model):
          assert "SkyModel" in str(sky_model)
  
      @staticmethod
      def test_parameters(sky_model):
          # Check that model parameters are references to the spatial and spectral parts
          p1 = sky_model.parameters["lon_0"]
          p2 = sky_model.spatial_model.parameters["lon_0"]
          assert p1 is p2
  
          p1 = sky_model.parameters["amplitude"]
          p2 = sky_model.spectral_model.parameters["amplitude"]
          assert p1 is p2
  
      @staticmethod
      def test_evaluate_scalar(sky_model):
          lon = 3 * u.deg
          lat = 4 * u.deg
          energy = 1 * u.TeV
  
          q = sky_model.evaluate(lon, lat, energy)
  
          assert q.unit == "cm-2 s-1 TeV-1 sr-1"
          assert np.isscalar(q.value)
          assert_allclose(q.to_value("cm-2 s-1 TeV-1 deg-2"), 1.76879232e-13)
  
      @staticmethod
      def test_evaluate_array(sky_model):
          lon = 3 * u.deg * np.ones(shape=(3, 4))
          lat = 4 * u.deg * np.ones(shape=(3, 4))
          energy = [1, 1, 1, 1, 1] * u.TeV
  
          q = sky_model.evaluate(lon, lat, energy[:, np.newaxis, np.newaxis])
  
          assert q.shape == (5, 3, 4)
          assert_allclose(q.to_value("cm-2 s-1 TeV-1 deg-2"), 1.76879232e-13)
  
      @staticmethod
      def test_processing(sky_model):
          assert sky_model.apply_irf == {"exposure": True, "psf": True, "edisp": True}
          out = sky_model.to_dict()
          assert "apply_irf" not in out
  
          sky_model.apply_irf["edisp"] = False
          out = sky_model.to_dict()
          assert out["apply_irf"] == {"exposure": True, "psf": True, "edisp": False}
          sky_model.apply_irf["edisp"] = True
  
  
  class Test_Template_with_cube:
      @staticmethod
      def test_evaluate_scalar(diffuse_model):
          # Check pixel inside map
          val = diffuse_model.evaluate(0
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/tests/test_cube.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.time import Time
  from gammapy.data.gti import GTI
  from gammapy.datasets.map import MapEvaluator
  from gammapy.irf import EDispKernel, PSFKernel
  from gammapy.maps import Map, MapAxis, WcsGeom
  from gammapy.modeling.models import (
      BackgroundModel,
      ConstantSpectralModel,
      ConstantTemporalModel,
      GaussianSpatialModel,
      Models,
      PointSpatialModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      TemplateSpatialModel,
      CompoundSpectralModel,
      SkyModel,
      create_fermi_isotropic_diffuse_model,
  )
  from gammapy.utils.testing import requires_data
  
  
  @pytest.fixture(scope="session")
  def sky_model():
      spatial_model = GaussianSpatialModel(
          lon_0="3 deg", lat_0="4 deg", sigma="3 deg", frame="galactic"
      )
      spectral_model = PowerLawSpectralModel(
          index=2, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      temporal_model = ConstantTemporalModel()
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
          temporal_model=temporal_model,
          name="source-1",
      )
  
  
  @pytest.fixture(scope="session")
  def gti():
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      t_ref = Time(55555, format="mjd")
      gti = GTI.create(start, stop, reference_time=t_ref)
      return gti
  
  
  @pytest.fixture(scope="session")
  def diffuse_model():
      axis = MapAxis.from_nodes([0.1, 100], name="energy_true", unit="TeV", interp="log")
      m = Map.create(
          npix=(4, 3), binsz=2, axes=[axis], unit="cm-2 s-1 MeV-1 sr-1", frame="galactic"
      )
      m.data += 42
      spatial_model = TemplateSpatialModel(m, normalize=False)
      return SkyModel(
          PowerLawNormSpectralModel(),
          spatial_model
      )
  
  
  @pytest.fixture(scope="session")
  def geom():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 3), unit=u.TeV, name="energy")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def geom_true():
      axis = MapAxis.from_edges(np.logspace(-1, 1, 4), unit=u.TeV, name="energy_true")
      return WcsGeom.create(skydir=(0, 0), npix=(5, 4), frame="galactic", axes=[axis])
  
  
  @pytest.fixture(scope="session")
  def exposure(geom_true):
      m = Map.from_geom(geom_true)
      m.quantity = np.ones(geom_true.data_shape) * u.Quantity("100 m2 s")
      m.data[1] *= 10
      return m
  
  
  @pytest.fixture(scope="session")
  def background(geom):
      m = Map.from_geom(geom)
      m.quantity = np.ones(geom.data_shape) * 1e-7
      return m
  
  
  @pytest.fixture(scope="session")
  def edisp(geom, geom_true):
      e_reco = geom.axes["energy"].edges
      e_true = geom_true.axes["energy_true"].edges
      return EDispKernel.from_diagonal_response(e_true=e_true, e_reco=e_reco)
  
  
  @pytest.fixture(scope="session")
  def psf(geom_true):
      sigma = 0.5 * u.deg
      return PSFKernel.from_gauss(geom_true, sigma)
  
  
  @pytest.fixture(scope="session")
  def evaluator(sky_model, exposure, psf, edisp, gti):
      return MapEvaluator(sky_model, exposure, psf=psf, edisp=edisp, gti=gti)
  
  
  @pytest.fixture(scope="session")
  def diffuse_evaluator(diffuse_model, exposure, psf, edisp):
      return MapEvaluator(diffuse_model, exposure, psf=psf, edisp=edisp)
  
  
  @pytest.fixture(scope="session")
  def sky_models(sky_model):
      sky_model_2 = sky_model.copy(name="source-2")
      sky_model_3 = sky_model.copy(name="source-3")
      return Models([sky_model_2, sky_model_3])
  
  
  @pytest.fixture(scope="session")
  def sky_models_2(sky_model):
      sky_model_4 = sky_model.copy(name="source-4")
      sky_model_5 = sky_model.copy(name="source-5")
      return Models([sky_model_4, sky_model_5])
  
  
  def test_sky_model_init():
      with pytest.raises(TypeError):
          spatial_model = GaussianSpatialModel()
          SkyModel(spectral_model=1234, spatial_model=spatial_model)
  
      with pytest.raises(TypeError):
          SkyModel(spectral_model=PowerLawSpectralModel(), spatial_model=1234)
  
  
  def test_sky_model_spatial_none_io(tmpdir):
      pwl = PowerLawSpectralModel()
      model = SkyModel(spectral_model=pwl, name="test")
      models = Models([model])
  
      filename = tmpdir / "test-models-none.yaml"
      models.write(filename)
  
      models = Models.read(filename)
  
      assert models["test"].spatial_model is None
  
  
  def test_sky_model_spatial_none_evaluate(geom_true, gti):
      pwl = PowerLawSpectralModel()
      model = SkyModel(spectral_model=pwl, name="test")
  
      data = model.evaluate_geom(geom_true, gti).to_value("cm-2 s-1 TeV-1")
  
      assert data.shape == (3, 1, 1)
      assert_allclose(data[0], 1.256774e-11, rtol=1e-6)
  
  
  def test_skymodel_addition(sky_model, sky_models, sky_models_2, diffuse_model):
      models = sky_model + sky_model.copy()
      assert isinstance(models, Models)
      assert len(models) == 2
  
      models = sky_model + sky_models
      assert isinstance(models, Models)
      assert len(models) == 3
  
      models = sky_models + sky_model
      assert isinstance(models, Models)
      assert len(models) == 3
  
      models = sky_models + diffuse_model
      assert isinstance(models, Models)
      assert len(models) == 3
  
      models = sky_models + sky_models_2
      assert isinstance(models, Models)
      assert len(models) == 4
  
      models = sky_model + sky_models
      assert isinstance(models, Models)
      assert len(models) == 3
  
  
  def test_background_model(background):
      bkg1 = BackgroundModel(background)
      bkg1.spectral_model.norm.value = 2.0
      npred1 = bkg1.evaluate()
      assert_allclose(npred1.data[0][0][0], background.data[0][0][0] * 2.0, rtol=1e-3)
      assert_allclose(npred1.data.sum(), background.data.sum() * 2.0, rtol=1e-3)
  
      bkg2 = BackgroundModel(background)
      bkg2.spectral_model.norm.value = 2.0
      bkg2.spectral_model.tilt.value = 0.2
      bkg2.spectral_model.reference.quantity = "1000 GeV"
  
      npred2 = bkg2.evaluate()
      assert_allclose(npred2.data[0][0][0], 2.254e-07, rtol=1e-3)
      assert_allclose(npred2.data.sum(), 7.352e-06, rtol=1e-3)
  
  
  def test_background_model_io(tmpdir, background):
      filename = str(tmpdir / "test-bkg-file.fits")
      bkg = BackgroundModel(background, filename=filename)
      bkg.spectral_model.norm.value = 2.0
      bkg.map.write(filename, overwrite=True)
      bkg_dict = bkg.to_dict()
      bkg_read = bkg.from_dict(bkg_dict)
  
      assert_allclose(
          bkg_read.evaluate().data.sum(), background.data.sum() * 2.0, rtol=1e-3
      )
      assert bkg_read.filename == filename
  
  
  class TestSkyModels:
      @staticmethod
      def test_parameters(sky_models):
          parnames = [
              "index",
              "amplitude",
              "reference",
              "lon_0",
              "lat_0",
              "sigma",
              "e",
              "phi",
          ] * 2
          assert sky_models.parameters.names == parnames
  
          # Check that model parameters are references to the parts
          p1 = sky_models.parameters["lon_0"]
          p2 = sky_models[0].parameters["lon_0"]
          assert p1 is p2
  
      @staticmethod
      def test_str(sky_models):
          assert "Component 0" in str(sky_models)
          assert "Component 1" in str(sky_models)
  
      @staticmethod
      def test_get_item(sky_models):
          model = sky_models["source-2"]
          assert model.name == "source-2"
  
          model = sky_models["source-3"]
          assert model.name == "source-3"
  
          with pytest.raises(ValueError):
              sky_models["spam"]
  
      @staticmethod
      def test_names(sky_models):
          assert sky_models.names == ["source-2", "source-3"]
  
  
  @requires_data()
  def test_models_mutation(sky_model, sky_models, sky_models_2):
      mods = sky_models
  
      mods.insert(0, sky_model)
      assert mods.names == ["source-1", "source-2", "source-3"]
  
      mods.extend(sky_models_2)
      assert mods.names == ["source-1", "source-2", "source-3", "source-4", "source-5"]
  
      mod3 = mods[3]
      mods.remove(mods[3])
      assert mods.names == ["source-1", "source-2", "source-3", "source-5"]
      mods.append(mod3)
      assert mods.names == ["source-1", "source-2", "source-3", "source-5", "source-4"]
      mods.pop(3)
      assert mods.names == ["source-1", "source-2", "source-3", "source-4"]
  
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods.append(sky_model)
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods.insert(0, sky_model)
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods.extend(sky_models_2)
      with pytest.raises(ValueError, match="Model names must be unique"):
          mods = sky_models + sky_models_2
  
  
  class TestSkyModel:
      @staticmethod
      def test_repr(sky_model):
          assert "SkyModel" in repr(sky_model)
  
      @staticmethod
      def test_str(sky_model):
          assert "SkyModel" in str(sky_model)
  
      @staticmethod
      def test_parameters(sky_model):
          # Check that model parameters are references to the spatial and spectral parts
          p1 = sky_model.parameters["lon_0"]
          p2 = sky_model.spatial_model.parameters["lon_0"]
          assert p1 is p2
  
          p1 = sky_model.parameters["amplitude"]
          p2 = sky_model.spectral_model.parameters["amplitude"]
          assert p1 is p2
  
      @staticmethod
      def test_evaluate_scalar(sky_model):
          lon = 3 * u.deg
          lat = 4 * u.deg
          energy = 1 * u.TeV
  
          q = sky_model.evaluate(lon, lat, energy)
  
          assert q.unit == "cm-2 s-1 TeV-1 sr-1"
          assert np.isscalar(q.value)
          assert_allclose(q.to_value("cm-2 s-1 TeV-1 deg-2"), 1.76879232e-13)
  
      @staticmethod
      def test_evaluate_array(sky_model):
          lon = 3 * u.deg * np.ones(shape=(3, 4))
          lat = 4 * u.deg * np.ones(shape=(3, 4))
          energy = [1, 1, 1, 1, 1] * u.TeV
  
          q = sky_model.evaluate(lon, lat, energy[:, np.newaxis, np.newaxis])
  
          assert q.shape == (5, 3, 4)
          assert_allclose(q.to_value("cm-2 s-1 TeV-1 deg-2"), 1.76879232e-13)
  
      @staticmethod
      def test_processing(sky_model):
          assert sky_model.apply_irf == {"exposure": True, "psf": True, "edisp": True}
          out = sky_model.to_dict()
          assert "apply_irf" not in out
  
          sky_model.apply_irf["edisp"] = False
          out = sky_model.to_dict()
          assert out["apply_irf"] == {"exposure": True, "psf": True, "edisp": False}
          sky_model.apply_irf["edisp"] = True
  
  
  class Test_Template_with_cube:
      @staticmethod
      def test_evaluate_scalar(diffuse_model):
          # Check pixel inside map
          val = diffuse
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/tests/test_cube.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.table import Table
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling.models import (
      MODEL_REGISTRY,
      PowerLawSpectralModel,
      AbsorbedSpectralModel,
      Absorption,
      BackgroundModel,
      Model,
      Models,
  )
  from gammapy.utils.scripts import read_yaml, write_yaml
  from gammapy.utils.testing import requires_data
  
  
  @requires_data()
  def test_dict_to_skymodels():
      filename = get_pkg_data_filename("data/examples.yaml")
      models_data = read_yaml(filename)
      models = Models.from_dict(models_data)
  
      assert len(models) == 5
  
      model0 = models[0]
      assert isinstance(model0, BackgroundModel)
      assert model0.name == "background_irf"
  
      model0 = models[1]
      assert "ExpCutoffPowerLawSpectralModel" in model0.spectral_model.tag
      assert "PointSpatialModel" in model0.spatial_model.tag
  
      pars0 = model0.parameters
      assert pars0["index"].value == 2.1
      assert pars0["index"].unit == ""
      assert np.isnan(pars0["index"].max)
      assert np.isnan(pars0["index"].min)
      assert not pars0["index"].frozen
  
      assert pars0["lon_0"].value == -0.5
      assert pars0["lon_0"].unit == "deg"
      assert pars0["lon_0"].max == 180.0
      assert pars0["lon_0"].min == -180.0
      assert pars0["lon_0"].frozen
  
      assert pars0["lat_0"].value == -0.0005
      assert pars0["lat_0"].unit == "deg"
      assert pars0["lat_0"].max == 90.0
      assert pars0["lat_0"].min == -90.0
      assert pars0["lat_0"].frozen
  
      assert pars0["lambda_"].value == 0.006
      assert pars0["lambda_"].unit == "TeV-1"
      assert np.isnan(pars0["lambda_"].min)
      assert np.isnan(pars0["lambda_"].max)
  
      model1 = models[2]
      assert "pl" in model1.spectral_model.tag
      assert "PowerLawSpectralModel" in model1.spectral_model.tag
      assert "DiskSpatialModel" in model1.spatial_model.tag
      assert "disk" in model1.spatial_model.tag
      assert "LightCurveTemplateTemporalModel" in model1.temporal_model.tag
  
      pars1 = model1.parameters
      assert pars1["index"].value == 2.2
      assert pars1["index"].unit == ""
      assert pars1["lat_0"].scale == 1.0
      assert pars1["lat_0"].factor == pars1["lat_0"].value
  
      assert np.isnan(pars1["index"].max)
      assert np.isnan(pars1["index"].min)
  
      assert pars1["r_0"].unit == "deg"
  
      model2 = models[3]
      assert_allclose(model2.spectral_model.energy.data, [34.171, 44.333, 57.517])
      assert model2.spectral_model.energy.unit == "MeV"
      assert_allclose(
          model2.spectral_model.values.data, [2.52894e-06, 1.2486e-06, 6.14648e-06]
      )
      assert model2.spectral_model.values.unit == "1 / (cm2 MeV s)"
  
      assert "TemplateSpectralModel" in model2.spectral_model.tag
      assert "TemplateSpatialModel" in model2.spatial_model.tag
  
      assert not model2.spatial_model.normalize
  
  
  @requires_data()
  def test_sky_models_io(tmp_path):
      # TODO: maybe change to a test case where we create a model programatically?
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      models.covariance = np.eye(len(models.parameters))
      models.write(tmp_path / "tmp.yaml")
      models = Models.read(tmp_path / "tmp.yaml")
      assert models._covar_file == "tmp_covariance.dat"
      assert_allclose(models.covariance.data, np.eye(len(models.parameters)))
      assert_allclose(models.parameters["lat_0"].min, -90.0)
  
      # TODO: not sure if we should just round-trip, or if we should
      # check YAML file content (e.g. against a ref file in the repo)
      # or check serialised dict content
  
  
  @requires_data()
  def test_absorption_io(tmp_path):
      dominguez = Absorption.read_builtin("dominguez")
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("pl", "spectral"),
          absorption=dominguez,
          redshift=0.5,
      )
      assert len(model.parameters) == 5
  
      model_dict = model.to_dict()
      parnames = [_["name"] for _ in model_dict["parameters"]]
      assert parnames == ["redshift", "alpha_norm"]
  
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert new_model.redshift.value == 0.5
      assert new_model.alpha_norm.name == "alpha_norm"
      assert new_model.alpha_norm.value == 1
      assert "PowerLawSpectralModel" in new_model.spectral_model.tag
      assert_allclose(new_model.absorption.energy, dominguez.energy)
      assert_allclose(new_model.absorption.param, dominguez.param)
      assert len(new_model.parameters) == 5
  
      test_absorption = Absorption(
          u.Quantity(range(3), "keV"),
          u.Quantity(range(2), ""),
          u.Quantity(np.ones((2, 3)), ""),
      )
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
          absorption=test_absorption,
          redshift=0.5,
      )
      model_dict = model.to_dict()
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert_allclose(new_model.absorption.energy, test_absorption.energy)
      assert_allclose(new_model.absorption.param, test_absorption.param)
  
      write_yaml(model_dict, tmp_path / "tmp.yaml")
      read_yaml(tmp_path / "tmp.yaml")
  
  
  def make_all_models():
      """Make an instance of each model, for testing."""
      yield Model.create("ConstantSpatialModel", "spatial")
      map_constantmodel = Map.create(npix=(10, 20), unit="sr-1")
      yield Model.create("TemplateSpatialModel", "spatial", map=map_constantmodel)
      yield Model.create(
          "DiskSpatialModel", "spatial", lon_0="1 deg", lat_0="2 deg", r_0="3 deg"
      )
      yield Model.create("gauss", "spatial", lon_0="1 deg", lat_0="2 deg", sigma="3 deg")
      yield Model.create("PointSpatialModel", "spatial", lon_0="1 deg", lat_0="2 deg")
      yield Model.create(
          "ShellSpatialModel",
          "spatial",
          lon_0="1 deg",
          lat_0="2 deg",
          radius="3 deg",
          width="4 deg",
      )
      yield Model.create("ConstantSpectralModel", "spectral", const="99 cm-2 s-1 TeV-1")
      yield Model.create(
          "CompoundSpectralModel",
          "spectral",
          model1=Model.create("PowerLawSpectralModel", "spectral"),
          model2=Model.create("PowerLawSpectralModel", "spectral"),
          operator=np.add,
      )
      yield Model.create("PowerLawSpectralModel", "spectral")
      yield Model.create("PowerLawNormSpectralModel", "spectral")
      yield Model.create("PowerLaw2SpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLawSpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLawNormSpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLaw3FGLSpectralModel", "spectral")
      yield Model.create("SuperExpCutoffPowerLaw3FGLSpectralModel", "spectral")
      yield Model.create("SuperExpCutoffPowerLaw4FGLSpectralModel", "spectral")
      yield Model.create("LogParabolaSpectralModel", "spectral")
      yield Model.create("LogParabolaNormSpectralModel", "spectral")
      yield Model.create(
          "TemplateSpectralModel", "spectral", energy=[1, 2] * u.cm, values=[3, 4] * u.cm
      )  # TODO: add unit validation?
      yield Model.create("GaussianSpectralModel", "spectral")
      # TODO: yield Model.create("AbsorbedSpectralModel")
      # TODO: yield Model.create("NaimaSpectralModel")
      # TODO: yield Model.create("ScaleSpectralModel")
      yield Model.create("ConstantTemporalModel", "temporal")
      yield Model.create("LightCurveTemplateTemporalModel", "temporal", Table())
      yield Model.create(
          "SkyModel",
          spatial_model=Model.create("ConstantSpatialModel", "spatial"),
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
      )
      m1 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_nodes([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("TemplateSpatialModel", "spatial", map=m1)
      m2 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_edges([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("BackgroundModel", map=m2)
  
  
  @pytest.mark.parametrize("model_class", MODEL_REGISTRY)
  def test_all_model_classes(model_class):
      if isinstance(model_class.tag, list):
          assert model_class.tag[0] == model_class.__name__
      else:
          assert model_class.tag == model_class.__name__
  
  
  @pytest.mark.parametrize("model", make_all_models())
  def test_all_model_instances(model):
      tag = model.tag[0] if isinstance(model.tag, list) else model.tag
      assert tag == model.__class__.__name__
  
  
  @requires_data()
  def test_missing_parameters():
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      assert models["source1"].spatial_model.e in models.parameters
      assert len(models["source1"].spatial_model.parameters) == 6
  
  
  def test_simplified_output():
      model = PowerLawSpectralModel()
      full = model.to_dict()
      simplified = model.to_dict(full_output=False)
      for k, name in enumerate(model.parameters.names):
          for item in ["min", "max", "frozen", "error"]:
              assert item in full["parameters"][k]
              assert item not in simplified["parameters"][k]
  
  
  def test_registries_print():
      print(MODEL_REGISTRY)
  
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/tests/test_io.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.table import Table
  from astropy.utils.data import get_pkg_data_filename
  from gammapy.maps import Map, MapAxis
  from gammapy.modeling.models import (
      MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      BackgroundModel,
      Model,
      Models,
  )
  from gammapy.utils.scripts import read_yaml, write_yaml
  from gammapy.utils.testing import requires_data
  
  
  @requires_data()
  def test_dict_to_skymodels():
      filename = get_pkg_data_filename("data/examples.yaml")
      models_data = read_yaml(filename)
      models = Models.from_dict(models_data)
  
      assert len(models) == 5
  
      model0 = models[0]
      assert isinstance(model0, BackgroundModel)
      assert model0.name == "background_irf"
  
      model0 = models[1]
      assert "ExpCutoffPowerLawSpectralModel" in model0.spectral_model.tag
      assert "PointSpatialModel" in model0.spatial_model.tag
  
      pars0 = model0.parameters
      assert pars0["index"].value == 2.1
      assert pars0["index"].unit == ""
      assert np.isnan(pars0["index"].max)
      assert np.isnan(pars0["index"].min)
      assert not pars0["index"].frozen
  
      assert pars0["lon_0"].value == -0.5
      assert pars0["lon_0"].unit == "deg"
      assert pars0["lon_0"].max == 180.0
      assert pars0["lon_0"].min == -180.0
      assert pars0["lon_0"].frozen
  
      assert pars0["lat_0"].value == -0.0005
      assert pars0["lat_0"].unit == "deg"
      assert pars0["lat_0"].max == 90.0
      assert pars0["lat_0"].min == -90.0
      assert pars0["lat_0"].frozen
  
      assert pars0["lambda_"].value == 0.006
      assert pars0["lambda_"].unit == "TeV-1"
      assert np.isnan(pars0["lambda_"].min)
      assert np.isnan(pars0["lambda_"].max)
  
      model1 = models[2]
      assert "pl" in model1.spectral_model.tag
      assert "PowerLawSpectralModel" in model1.spectral_model.tag
      assert "DiskSpatialModel" in model1.spatial_model.tag
      assert "disk" in model1.spatial_model.tag
      assert "LightCurveTemplateTemporalModel" in model1.temporal_model.tag
  
      pars1 = model1.parameters
      assert pars1["index"].value == 2.2
      assert pars1["index"].unit == ""
      assert pars1["lat_0"].scale == 1.0
      assert pars1["lat_0"].factor == pars1["lat_0"].value
  
      assert np.isnan(pars1["index"].max)
      assert np.isnan(pars1["index"].min)
  
      assert pars1["r_0"].unit == "deg"
  
      model2 = models[3]
      assert_allclose(model2.spectral_model.energy.data, [34.171, 44.333, 57.517])
      assert model2.spectral_model.energy.unit == "MeV"
      assert_allclose(
          model2.spectral_model.values.data, [2.52894e-06, 1.2486e-06, 6.14648e-06]
      )
      assert model2.spectral_model.values.unit == "1 / (cm2 MeV s)"
  
      assert "TemplateSpectralModel" in model2.spectral_model.tag
      assert "TemplateSpatialModel" in model2.spatial_model.tag
  
      assert not model2.spatial_model.normalize
  
  
  @requires_data()
  def test_sky_models_io(tmp_path):
      # TODO: maybe change to a test case where we create a model programatically?
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      models.covariance = np.eye(len(models.parameters))
      models.write(tmp_path / "tmp.yaml")
      models = Models.read(tmp_path / "tmp.yaml")
      assert models._covar_file == "tmp_covariance.dat"
      assert_allclose(models.covariance.data, np.eye(len(models.parameters)))
      assert_allclose(models.parameters["lat_0"].min, -90.0)
  
      # TODO: not sure if we should just round-trip, or if we should
      # check YAML file content (e.g. against a ref file in the repo)
      # or check serialised dict content
  
  
  @requires_data()
  def test_absorption_io(tmp_path):
      dominguez = Absorption.read_builtin("dominguez")
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("pl", "spectral"),
          absorption=dominguez,
          redshift=0.5,
      )
      assert len(model.parameters) == 5
  
      model_dict = model.to_dict()
      parnames = [_["name"] for _ in model_dict["parameters"]]
      assert parnames == ["redshift", "alpha_norm"]
  
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert new_model.redshift.value == 0.5
      assert new_model.alpha_norm.name == "alpha_norm"
      assert new_model.alpha_norm.value == 1
      assert "PowerLawSpectralModel" in new_model.spectral_model.tag
      assert_allclose(new_model.absorption.energy, dominguez.energy)
      assert_allclose(new_model.absorption.param, dominguez.param)
      assert len(new_model.parameters) == 5
  
      test_absorption = Absorption(
          u.Quantity(range(3), "keV"),
          u.Quantity(range(2), ""),
          u.Quantity(np.ones((2, 3)), ""),
      )
      model = AbsorbedSpectralModel(
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
          absorption=test_absorption,
          redshift=0.5,
      )
      model_dict = model.to_dict()
      new_model = AbsorbedSpectralModel.from_dict(model_dict)
  
      assert_allclose(new_model.absorption.energy, test_absorption.energy)
      assert_allclose(new_model.absorption.param, test_absorption.param)
  
      write_yaml(model_dict, tmp_path / "tmp.yaml")
      read_yaml(tmp_path / "tmp.yaml")
  
  
  def make_all_models():
      """Make an instance of each model, for testing."""
      yield Model.create("ConstantSpatialModel", "spatial")
      map_constantmodel = Map.create(npix=(10, 20), unit="sr-1")
      yield Model.create("TemplateSpatialModel", "spatial", map=map_constantmodel)
      yield Model.create(
          "DiskSpatialModel", "spatial", lon_0="1 deg", lat_0="2 deg", r_0="3 deg"
      )
      yield Model.create("gauss", "spatial", lon_0="1 deg", lat_0="2 deg", sigma="3 deg")
      yield Model.create("PointSpatialModel", "spatial", lon_0="1 deg", lat_0="2 deg")
      yield Model.create(
          "ShellSpatialModel",
          "spatial",
          lon_0="1 deg",
          lat_0="2 deg",
          radius="3 deg",
          width="4 deg",
      )
      yield Model.create("ConstantSpectralModel", "spectral", const="99 cm-2 s-1 TeV-1")
      yield Model.create(
          "CompoundSpectralModel",
          "spectral",
          model1=Model.create("PowerLawSpectralModel", "spectral"),
          model2=Model.create("PowerLawSpectralModel", "spectral"),
          operator=np.add,
      )
      yield Model.create("PowerLawSpectralModel", "spectral")
      yield Model.create("PowerLawNormSpectralModel", "spectral")
      yield Model.create("PowerLaw2SpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLawSpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLawNormSpectralModel", "spectral")
      yield Model.create("ExpCutoffPowerLaw3FGLSpectralModel", "spectral")
      yield Model.create("SuperExpCutoffPowerLaw3FGLSpectralModel", "spectral")
      yield Model.create("SuperExpCutoffPowerLaw4FGLSpectralModel", "spectral")
      yield Model.create("LogParabolaSpectralModel", "spectral")
      yield Model.create("LogParabolaNormSpectralModel", "spectral")
      yield Model.create(
          "TemplateSpectralModel", "spectral", energy=[1, 2] * u.cm, values=[3, 4] * u.cm
      )  # TODO: add unit validation?
      yield Model.create("GaussianSpectralModel", "spectral")
      # TODO: yield Model.create("AbsorbedSpectralModel")
      # TODO: yield Model.create("NaimaSpectralModel")
      # TODO: yield Model.create("ScaleSpectralModel")
      yield Model.create("ConstantTemporalModel", "temporal")
      yield Model.create("LightCurveTemplateTemporalModel", "temporal", Table())
      yield Model.create(
          "SkyModel",
          spatial_model=Model.create("ConstantSpatialModel", "spatial"),
          spectral_model=Model.create("PowerLawSpectralModel", "spectral"),
      )
      m1 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_nodes([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("TemplateSpatialModel", "spatial", map=m1)
      m2 = Map.create(
          npix=(10, 20, 30), axes=[MapAxis.from_edges([1, 2] * u.TeV, name="energy")]
      )
      yield Model.create("BackgroundModel", map=m2)
  
  
  @pytest.mark.parametrize("model_class", MODEL_REGISTRY)
  def test_all_model_classes(model_class):
      if isinstance(model_class.tag, list):
          assert model_class.tag[0] == model_class.__name__
      else:
          assert model_class.tag == model_class.__name__
  
  
  @pytest.mark.parametrize("model", make_all_models())
  def test_all_model_instances(model):
      tag = model.tag[0] if isinstance(model.tag, list) else model.tag
      assert tag == model.__class__.__name__
  
  
  @requires_data()
  def test_missing_parameters():
      filename = get_pkg_data_filename("data/examples.yaml")
      models = Models.read(filename)
      assert models["source1"].spatial_model.e in models.parameters
      assert len(models["source1"].spatial_model.parameters) == 6
  
  
  def test_registries_print():
      print(MODEL_REGISTRY)
  
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/tests/test_io.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  import operator
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.maps import MapAxis
  from gammapy.modeling.models import (
      SPECTRAL_MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      ConstantSpectralModel,
      CompoundSpectralModel,
      ExpCutoffPowerLaw3FGLSpectralModel,
      ExpCutoffPowerLawSpectralModel,
      ExpCutoffPowerLawNormSpectralModel,
      GaussianSpectralModel,
      LogParabolaSpectralModel,
      LogParabolaNormSpectralModel,
      NaimaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      BrokenPowerLawSpectralModel,
      SmoothBrokenPowerLawSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
      TemplateSpectralModel,
  )
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      mpl_plot_check,
      requires_data,
      requires_dependency,
  )
  
  
  def table_model():
      energy = MapAxis.from_energy_bounds(0.1 * u.TeV, 100 * u.TeV, 1000).center
  
      model = PowerLawSpectralModel(
          index=2.3, amplitude="4 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      dnde = model(energy)
  
      return TemplateSpectralModel(energy, dnde)
  
  
  TEST_MODELS = [
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.6, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV cm-2 s-1"),
      ),
      dict(
          name="norm-powerlaw",
          model=PowerLawNormSpectralModel(
              tilt=2 * u.Unit(""), norm=4.0 * u.Unit(""), reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, ""),
          integral_1_10TeV=u.Quantity(3.6, "TeV"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV2"),
      ),
      dict(
          name="powerlaw2",
          model=PowerLaw2SpectralModel(
              amplitude=u.Quantity(2.9227116204223784, "cm-2 s-1"),
              index=2.3 * u.Unit(""),
              emin=1 * u.TeV,
              emax=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.6 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.765838739678921, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV cm-2 s-1"),
          e_peak=4 * u.TeV,
      ),
      dict(
          name="norm-ecpl",
          model=ExpCutoffPowerLawNormSpectralModel(
              index=1.6 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, ""),
          integral_1_10TeV=u.Quantity(3.765838739678921, "TeV"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV2"),
      ),
      dict(
          name="ecpl_3fgl",
          model=ExpCutoffPowerLaw3FGLSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              ecut=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.7349563611124971, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.6034046173089, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(5.340285560055799, "TeV cm-2 s-1"),
      ),
      dict(
          name="plsec_4fgl",
          model=SuperExpCutoffPowerLaw4FGLSpectralModel(
              index_1=1.5,
              index_2=2,
              amplitude=1 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              expfactor=1e-2,
          ),
          val_at_2TeV=u.Quantity(0.3431043087721737, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(1.2125247, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.38072082, "TeV cm-2 s-1"),
      ),
      dict(
          name="logpar",
          model=LogParabolaSpectralModel(
              alpha=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.255689748270628, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV cm-2 s-1"),
          e_peak=0.74082 * u.TeV,
      ),
      dict(
          name="norm-logpar",
          model=LogParabolaNormSpectralModel(
              alpha=2.3 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, ""),
          integral_1_10TeV=u.Quantity(2.255689748270628, "TeV"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV2"),
      ),
      dict(
          name="logpar10",
          model=LogParabolaSpectralModel.from_log10(
              alpha=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=1.151292546497023 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.255689748270628, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV cm-2 s-1"),
          e_peak=0.74082 * u.TeV,
      ),
      dict(
          name="constant",
          model=ConstantSpectralModel(const=4 / u.cm ** 2 / u.s / u.TeV),
          val_at_2TeV=u.Quantity(4, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(35.9999999999999, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(198.00000000000006, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw_index1",
          model=PowerLawSpectralModel(
              index=1 * u.Unit(""),
              amplitude=2 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(4.605170185, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(18.0, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl_2",
          model=ExpCutoffPowerLawSpectralModel(
              index=2.0 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.81873075, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.83075297, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.41406327, "TeV cm-2 s-1"),
          e_peak=np.nan * u.TeV,
      ),
      dict(
          name="GaussianSpectralModel",
          model=GaussianSpectralModel(
              norm=4 / u.cm ** 2 / u.s, mean=2 * u.TeV, sigma=0.2 * u.TeV
          ),
          val_at_2TeV=u.Quantity(7.978845608028654, "cm-2 s-1 TeV-1"),
          val_at_3TeV=u.Quantity(2.973439029468601e-05, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.9999988533937123, "cm-2 s-1"),
          integral_infinity=u.Quantity(4, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.999998896163037, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.8 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
              alpha=0.8,
          ),
          val_at_2TeV=u.Quantity(0.871694294554192, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.026342, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.38652453, "TeV cm-2 s-1"),
          e_peak=1.7677669529663684 * u.TeV,
      ),
      dict(
          name="bpl",
          model=BrokenPowerLawSpectralModel(
              index1=1.5 * u.Unit(""),
              index2=2.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.125, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(0.45649740094103286, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(0.9669999668731384, "TeV cm-2 s-1"),
      ),
      dict(
          name="sbpl",
          model=SmoothBrokenPowerLawSpectralModel(
              index1=1
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/tests/test_spectral.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  import operator
  from numpy.testing import assert_allclose
  import astropy.units as u
  from gammapy.maps import MapAxis
  from gammapy.modeling.models import (
      SPECTRAL_MODEL_REGISTRY,
      AbsorbedSpectralModel,
      Absorption,
      ConstantSpectralModel,
      CompoundSpectralModel,
      ExpCutoffPowerLaw3FGLSpectralModel,
      ExpCutoffPowerLawSpectralModel,
      ExpCutoffPowerLawNormSpectralModel,
      GaussianSpectralModel,
      LogParabolaSpectralModel,
      LogParabolaNormSpectralModel,
      NaimaSpectralModel,
      PowerLaw2SpectralModel,
      PowerLawSpectralModel,
      PowerLawNormSpectralModel,
      BrokenPowerLawSpectralModel,
      SmoothBrokenPowerLawSpectralModel,
      SuperExpCutoffPowerLaw4FGLSpectralModel,
      TemplateSpectralModel,
  )
  from gammapy.utils.testing import (
      assert_quantity_allclose,
      mpl_plot_check,
      requires_data,
      requires_dependency,
  )
  
  
  def table_model():
      energy = MapAxis.from_energy_bounds(0.1 * u.TeV, 100 * u.TeV, 1000).center
  
      model = PowerLawSpectralModel(
          index=2.3, amplitude="4 cm-2 s-1 TeV-1", reference="1 TeV"
      )
      dnde = model(energy)
  
      return TemplateSpectralModel(energy, dnde)
  
  
  TEST_MODELS = [
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw",
          model=PowerLawSpectralModel(
              index=2 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.6, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV cm-2 s-1"),
      ),
      dict(
          name="norm-powerlaw",
          model=PowerLawNormSpectralModel(
              tilt=2 * u.Unit(""), norm=4.0 * u.Unit(""), reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, ""),
          integral_1_10TeV=u.Quantity(3.6, "TeV"),
          eflux_1_10TeV=u.Quantity(9.210340371976184, "TeV2"),
      ),
      dict(
          name="powerlaw2",
          model=PowerLaw2SpectralModel(
              amplitude=u.Quantity(2.9227116204223784, "cm-2 s-1"),
              index=2.3 * u.Unit(""),
              emin=1 * u.TeV,
              emax=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(4 * 2.0 ** (-2.3), "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.9227116204223784, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.650836884969039, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.6 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.765838739678921, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV cm-2 s-1"),
          e_peak=4 * u.TeV,
      ),
      dict(
          name="norm-ecpl",
          model=ExpCutoffPowerLawNormSpectralModel(
              index=1.6 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.080321705479446, ""),
          integral_1_10TeV=u.Quantity(3.765838739678921, "TeV"),
          eflux_1_10TeV=u.Quantity(9.901735870666526, "TeV2"),
      ),
      dict(
          name="ecpl_3fgl",
          model=ExpCutoffPowerLaw3FGLSpectralModel(
              index=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              ecut=10 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.7349563611124971, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.6034046173089, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(5.340285560055799, "TeV cm-2 s-1"),
      ),
      dict(
          name="plsec_4fgl",
          model=SuperExpCutoffPowerLaw4FGLSpectralModel(
              index_1=1.5,
              index_2=2,
              amplitude=1 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              expfactor=1e-2,
          ),
          val_at_2TeV=u.Quantity(0.3431043087721737, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(1.2125247, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.38072082, "TeV cm-2 s-1"),
      ),
      dict(
          name="logpar",
          model=LogParabolaSpectralModel(
              alpha=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.255689748270628, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV cm-2 s-1"),
          e_peak=0.74082 * u.TeV,
      ),
      dict(
          name="norm-logpar",
          model=LogParabolaNormSpectralModel(
              alpha=2.3 * u.Unit(""),
              norm=4 * u.Unit(""),
              reference=1 * u.TeV,
              beta=0.5 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, ""),
          integral_1_10TeV=u.Quantity(2.255689748270628, "TeV"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV2"),
      ),
      dict(
          name="logpar10",
          model=LogParabolaSpectralModel.from_log10(
              alpha=2.3 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              beta=1.151292546497023 * u.Unit(""),
          ),
          val_at_2TeV=u.Quantity(0.6387956571420305, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.255689748270628, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(3.9586515834989267, "TeV cm-2 s-1"),
          e_peak=0.74082 * u.TeV,
      ),
      dict(
          name="constant",
          model=ConstantSpectralModel(const=4 / u.cm ** 2 / u.s / u.TeV),
          val_at_2TeV=u.Quantity(4, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(35.9999999999999, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(198.00000000000006, "TeV cm-2 s-1"),
      ),
      dict(
          name="powerlaw_index1",
          model=PowerLawSpectralModel(
              index=1 * u.Unit(""),
              amplitude=2 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(1.0, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(4.605170185, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(18.0, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl_2",
          model=ExpCutoffPowerLawSpectralModel(
              index=2.0 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.81873075, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(2.83075297, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(6.41406327, "TeV cm-2 s-1"),
          e_peak=np.nan * u.TeV,
      ),
      dict(
          name="GaussianSpectralModel",
          model=GaussianSpectralModel(
              norm=4 / u.cm ** 2 / u.s, mean=2 * u.TeV, sigma=0.2 * u.TeV
          ),
          val_at_2TeV=u.Quantity(7.978845608028654, "cm-2 s-1 TeV-1"),
          val_at_3TeV=u.Quantity(2.973439029468601e-05, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.9999988533937123, "cm-2 s-1"),
          integral_infinity=u.Quantity(4, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.999998896163037, "TeV cm-2 s-1"),
      ),
      dict(
          name="ecpl",
          model=ExpCutoffPowerLawSpectralModel(
              index=1.8 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              reference=1 * u.TeV,
              lambda_=0.1 / u.TeV,
              alpha=0.8,
          ),
          val_at_2TeV=u.Quantity(0.871694294554192, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(3.026342, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(7.38652453, "TeV cm-2 s-1"),
          e_peak=1.7677669529663684 * u.TeV,
      ),
      dict(
          name="bpl",
          model=BrokenPowerLawSpectralModel(
              index1=1.5 * u.Unit(""),
              index2=2.5 * u.Unit(""),
              amplitude=4 / u.cm ** 2 / u.s / u.TeV,
              ebreak=0.5 * u.TeV,
          ),
          val_at_2TeV=u.Quantity(0.125, "cm-2 s-1 TeV-1"),
          integral_1_10TeV=u.Quantity(0.45649740094103286, "cm-2 s-1"),
          eflux_1_10TeV=u.Quantity(0.9669999668731384, "TeV cm-2 s-1"),
      ),
      dict(
          name="sbpl",
          model=SmoothBrokenPowerLawSpectralModel(
              index1=1
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/tests/test_spectral.py
files.new_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  from astropy import units as u
  from astropy.table import Table
  from astropy.time import Time
  from gammapy.data.gti import GTI
  from gammapy.modeling.models import (
      ConstantTemporalModel,
      ExpDecayTemporalModel,
      GaussianTemporalModel,
      LightCurveTemplateTemporalModel,
      PowerLawSpectralModel,
      SkyModel,
  )
  from gammapy.utils.scripts import make_path
  from gammapy.utils.testing import requires_data, requires_dependency, mpl_plot_check
  
  
  # TODO: add light-curve test case from scratch
  # only use the FITS one for I/O (or not at all)
  @pytest.fixture(scope="session")
  def light_curve():
      path = "$GAMMAPY_DATA/tests/models/light_curve/lightcrv_PKSB1222+216.fits"
      return LightCurveTemplateTemporalModel.read(path)
  
  
  @requires_data()
  def test_light_curve_str(light_curve):
      ss = str(light_curve)
      assert "LightCurveTemplateTemporalModel" in ss
  
  
  @requires_data()
  def test_light_curve_evaluate(light_curve):
      t = Time(59500, format="mjd")
      val = light_curve(t)
      assert_allclose(val, 0.015512, rtol=1e-5)
  
      t = Time(46300, format="mjd")
      val = light_curve.evaluate(t.mjd, ext=3)
      assert_allclose(val, 0.01551196, rtol=1e-5)
  
  
  def rate(x, c="1e4 s"):
      c = u.Quantity(c)
      return np.exp(-x / c)
  
  
  def ph_curve(x, amplitude=0.5, x0=0.01):
      return 100.0 + amplitude * np.sin(2 * np.pi * (x - x0) / 1.0)
  
  
  def test_time_sampling(tmp_path):
      time = np.arange(0, 10, 0.06) * u.hour
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = rate(time)
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model = LightCurveTemplateTemporalModel(table)
  
      filename = str(make_path(tmp_path / "tmp.fits"))
      temporal_model.write(path=filename)
      model_read = temporal_model.read(filename)
      assert temporal_model.filename == filename
      assert model_read.filename == filename
      assert_allclose(model_read.table["TIME"].quantity.value, time.value)
  
      t_ref = "2010-01-01T00:00:00"
      t_min = "2010-01-01T00:00:00"
      t_max = "2010-01-01T08:00:00"
  
      sampler = temporal_model.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0, t_delta="10 min"
      )
  
      sampler = u.Quantity((sampler - Time(t_ref)).sec, "s")
      assert len(sampler) == 2
      assert_allclose(sampler.value, [12661.65802564, 7826.92991], rtol=1e-5)
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = np.ones(len(time))
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model_uniform = LightCurveTemplateTemporalModel(table)
  
      sampler_uniform = temporal_model_uniform.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0, t_delta="10 min"
      )
      sampler_uniform = u.Quantity((sampler_uniform - Time(t_ref)).sec, "s")
  
      assert len(sampler_uniform) == 2
      assert_allclose(sampler_uniform.value, [1261.65802564, 6026.9299098], rtol=1e-5)
  
  
  def test_lightcurve_temporal_model_integral():
      time = np.arange(0, 10, 0.06) * u.hour
      table = Table()
      table["TIME"] = time
      table["NORM"] = np.ones(len(time))
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model = LightCurveTemplateTemporalModel(table)
  
      start = [1, 3, 5] * u.hour
      stop = [2, 3.5, 6] * u.hour
      gti = GTI.create(start, stop, reference_time=Time("2010-01-01T00:00:00"))
  
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 1.0, rtol=1e-5)
  
  
  def test_constant_temporal_model_sample():
      temporal_model = ConstantTemporalModel()
  
      t_ref = "2010-01-01T00:00:00"
      t_min = "2010-01-01T00:00:00"
      t_max = "2010-01-01T08:00:00"
  
      sampler = temporal_model.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0
      )
  
      sampler = u.Quantity((sampler - Time(t_ref)).sec, "s")
  
      assert len(sampler) == 2
      assert_allclose(sampler.value, [15805.82891311, 20597.45375153], rtol=1e-5)
  
  
  def test_constant_temporal_model_evaluate():
      temporal_model = ConstantTemporalModel()
      t = Time(46300, format="mjd")
      val = temporal_model(t)
      assert_allclose(val, 1.0, rtol=1e-5)
  
  
  def test_constant_temporal_model_integral():
      temporal_model = ConstantTemporalModel()
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      gti = GTI.create(start, stop)
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 1.0, rtol=1e-5)
  
  
  def test_exponential_temporal_model_evaluate():
      t = Time(46301, format="mjd")
      t_ref = 46300 * u.d
      t0 = 2.0 * u.d
      temporal_model = ExpDecayTemporalModel(t_ref=t_ref, t0=t0)
      val = temporal_model(t)
      assert_allclose(val, 0.6065306597126334, rtol=1e-5)
  
  
  def test_exponential_temporal_model_integral():
      t_ref = Time(55555, format="mjd")
  
      temporal_model = ExpDecayTemporalModel(t_ref=t_ref.mjd * u.d)
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      gti = GTI.create(start, stop, reference_time=t_ref)
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 0.102557, rtol=1e-5)
  
  
  def test_gaussian_temporal_model_evaluate():
      t = Time(46301, format="mjd")
      t_ref = 46300 * u.d
      sigma = 2.0 * u.d
      temporal_model = GaussianTemporalModel(t_ref=t_ref, sigma=sigma)
      val = temporal_model(t)
      assert_allclose(val, 0.882497, rtol=1e-5)
  
  
  def test_gaussian_temporal_model_integral():
      temporal_model = GaussianTemporalModel(t_ref=50003 * u.d, sigma="2.0 day")
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      t_ref = Time(50000, format="mjd")
      gti = GTI.create(start, stop, reference_time=t_ref)
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 0.682679, rtol=1e-5)
  
  
  @requires_data()
  def test_to_dict(light_curve):
  
      out = light_curve.to_dict()
      assert out["type"] == "LightCurveTemplateTemporalModel"
      assert "lightcrv_PKSB1222+216.fits" in out["filename"]
  
  
  @requires_data()
  def test_with_skymodel(light_curve):
  
      sky_model = SkyModel(spectral_model=PowerLawSpectralModel())
      out = sky_model.to_dict()
      assert "temporal" not in out
  
      sky_model = SkyModel(
          spectral_model=PowerLawSpectralModel(), temporal_model=light_curve
      )
      assert "LightCurveTemplateTemporalModel" in sky_model.temporal_model.tag
  
      out = sky_model.to_dict()
      assert "temporal" in out
  
  
  @requires_dependency("matplotlib")
  def test_plot_constant_model():
      time_range = [Time.now(), Time.now() + 1 * u.d]
      constant_model = ConstantTemporalModel(const=1)
      with mpl_plot_check():
          constant_model.plot(time_range)
  
  
  [... the rest is omitted ...]
files.new_path: gammapy/modeling/models/tests/test_temporal.py
files.old_contents: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  from astropy import units as u
  from astropy.table import Table
  from astropy.time import Time
  from gammapy.data.gti import GTI
  from gammapy.modeling.models import (
      ConstantTemporalModel,
      ExpDecayTemporalModel,
      GaussianTemporalModel,
      LightCurveTemplateTemporalModel,
      PowerLawSpectralModel,
      SkyModel,
  )
  from gammapy.utils.scripts import make_path
  from gammapy.utils.testing import requires_data, requires_dependency, mpl_plot_check
  
  
  # TODO: add light-curve test case from scratch
  # only use the FITS one for I/O (or not at all)
  @pytest.fixture(scope="session")
  def light_curve():
      path = "$GAMMAPY_DATA/tests/models/light_curve/lightcrv_PKSB1222+216.fits"
      return LightCurveTemplateTemporalModel.read(path)
  
  
  @requires_data()
  def test_light_curve_str(light_curve):
      ss = str(light_curve)
      assert "LightCurveTemplateTemporalModel" in ss
  
  
  @requires_data()
  def test_light_curve_evaluate(light_curve):
      t = Time(59500, format="mjd")
      val = light_curve(t)
      assert_allclose(val, 0.015512, rtol=1e-5)
  
      t = Time(46300, format="mjd")
      val = light_curve.evaluate(t.mjd, ext=3)
      assert_allclose(val, 0.01551196, rtol=1e-5)
  
  
  def rate(x, c="1e4 s"):
      c = u.Quantity(c)
      return np.exp(-x / c)
  
  
  def ph_curve(x, amplitude=0.5, x0=0.01):
      return 100.0 + amplitude * np.sin(2 * np.pi * (x - x0) / 1.0)
  
  
  def test_time_sampling(tmp_path):
      time = np.arange(0, 10, 0.06) * u.hour
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = rate(time)
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model = LightCurveTemplateTemporalModel(table)
  
      filename = str(make_path(tmp_path / "tmp.fits"))
      temporal_model.write(path=filename)
      model_read = temporal_model.read(filename)
      assert temporal_model.filename == filename
      assert model_read.filename == filename
      assert_allclose(model_read.table["TIME"].quantity.value, time.value)
  
      t_ref = "2010-01-01T00:00:00"
      t_min = "2010-01-01T00:00:00"
      t_max = "2010-01-01T08:00:00"
  
      sampler = temporal_model.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0, t_delta="10 min"
      )
  
      sampler = u.Quantity((sampler - Time(t_ref)).sec, "s")
      assert len(sampler) == 2
      assert_allclose(sampler.value, [12661.65802564, 7826.92991], rtol=1e-5)
  
      table = Table()
      table["TIME"] = time
      table["NORM"] = np.ones(len(time))
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model_uniform = LightCurveTemplateTemporalModel(table)
  
      sampler_uniform = temporal_model_uniform.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0, t_delta="10 min"
      )
      sampler_uniform = u.Quantity((sampler_uniform - Time(t_ref)).sec, "s")
  
      assert len(sampler_uniform) == 2
      assert_allclose(sampler_uniform.value, [1261.65802564, 6026.9299098], rtol=1e-5)
  
  
  def test_lightcurve_temporal_model_integral():
      time = np.arange(0, 10, 0.06) * u.hour
      table = Table()
      table["TIME"] = time
      table["NORM"] = np.ones(len(time))
      table.meta = dict(MJDREFI=55197.0, MJDREFF=0, TIMEUNIT="hour")
      temporal_model = LightCurveTemplateTemporalModel(table)
  
      start = [1, 3, 5] * u.hour
      stop = [2, 3.5, 6] * u.hour
      gti = GTI.create(start, stop, reference_time=Time("2010-01-01T00:00:00"))
  
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 1.0, rtol=1e-5)
  
  
  def test_constant_temporal_model_sample():
      temporal_model = ConstantTemporalModel()
  
      t_ref = "2010-01-01T00:00:00"
      t_min = "2010-01-01T00:00:00"
      t_max = "2010-01-01T08:00:00"
  
      sampler = temporal_model.sample_time(
          n_events=2, t_min=t_min, t_max=t_max, random_state=0
      )
  
      sampler = u.Quantity((sampler - Time(t_ref)).sec, "s")
  
      assert len(sampler) == 2
      assert_allclose(sampler.value, [15805.82891311, 20597.45375153], rtol=1e-5)
  
  
  def test_constant_temporal_model_evaluate():
      temporal_model = ConstantTemporalModel()
      t = Time(46300, format="mjd")
      val = temporal_model(t)
      assert_allclose(val, 1.0, rtol=1e-5)
  
  
  def test_constant_temporal_model_integral():
      temporal_model = ConstantTemporalModel()
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      gti = GTI.create(start, stop)
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 1.0, rtol=1e-5)
  
  
  def test_exponential_temporal_model_evaluate():
      t = Time(46301, format="mjd")
      t_ref = 46300 * u.d
      t0 = 2.0 * u.d
      temporal_model = ExpDecayTemporalModel(t_ref=t_ref, t0=t0)
      val = temporal_model(t)
      assert_allclose(val, 0.6065306597126334, rtol=1e-5)
  
  
  def test_exponential_temporal_model_integral():
      t_ref = Time(55555, format="mjd")
  
      temporal_model = ExpDecayTemporalModel(t_ref=t_ref.mjd * u.d)
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      gti = GTI.create(start, stop, reference_time=t_ref)
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 0.102557, rtol=1e-5)
  
  
  def test_gaussian_temporal_model_evaluate():
      t = Time(46301, format="mjd")
      t_ref = 46300 * u.d
      sigma = 2.0 * u.d
      temporal_model = GaussianTemporalModel(t_ref=t_ref, sigma=sigma)
      val = temporal_model(t)
      assert_allclose(val, 0.882497, rtol=1e-5)
  
  
  def test_gaussian_temporal_model_integral():
      temporal_model = GaussianTemporalModel(t_ref=50003 * u.d, sigma="2.0 day")
      start = [1, 3, 5] * u.day
      stop = [2, 3.5, 6] * u.day
      t_ref = Time(50000, format="mjd")
      gti = GTI.create(start, stop, reference_time=t_ref)
      val = temporal_model.integral(gti.time_start, gti.time_stop)
      assert len(val) == 3
      assert_allclose(np.sum(val), 0.682679, rtol=1e-5)
  
  
  @requires_data()
  def test_to_dict(light_curve):
  
      out = light_curve.to_dict()
      assert out["type"] == "LightCurveTemplateTemporalModel"
      assert "lightcrv_PKSB1222+216.fits" in out["filename"]
  
  
  @requires_data()
  def test_with_skymodel(light_curve):
  
      sky_model = SkyModel(spectral_model=PowerLawSpectralModel())
      out = sky_model.to_dict()
      assert "temporal" not in out
  
      sky_model = SkyModel(
          spectral_model=PowerLawSpectralModel(), temporal_model=light_curve
      )
      assert "LightCurveTemplateTemporalModel" in sky_model.temporal_model.tag
  
      out = sky_model.to_dict()
      assert "temporal" in out
  
  
  @requires_dependency("matplotlib")
  def test_plot_constant_model():
      time_range = [Time.now(), Time.now() + 1 * u.d]
      constant_model = ConstantTemporalModel(const=1)
      with mpl_plot_check():
          constant_model.plot(time_range)
  
  
  
  [... the rest is omitted ...]
files.old_path: gammapy/modeling/models/tests/test_temporal.py