Source code for gammapy.modeling.models.cube

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Cube models (axes: lon, lat, energy)."""

import logging
import os
import numpy as np
import astropy.units as u
from astropy.nddata import NoOverlapError
from astropy.time import Time
from gammapy.maps import Map, MapAxis, WcsGeom
from gammapy.modeling import Covariance, Parameters
from gammapy.modeling.parameter import _get_parameters_str
from gammapy.utils.fits import LazyFitsData
from gammapy.utils.scripts import make_name, make_path
from .core import Model, ModelBase, Models
from .spatial import ConstantSpatialModel, SpatialModel
from .spectral import PowerLawNormSpectralModel, SpectralModel, TemplateSpectralModel
from .temporal import TemporalModel

log = logging.getLogger(__name__)


__all__ = [
    "create_fermi_isotropic_diffuse_model",
    "FoVBackgroundModel",
    "SkyModel",
    "TemplateNPredModel",
]


[docs]class SkyModel(ModelBase): """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", "sky-model"] _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, covariance_data=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__(covariance_data=covariance_data) @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): axis = MapAxis.from_energy_bounds( "0.1 TeV", "10 TeV", nbin=1, name="energy_true" ) geom = WcsGeom.create(skydir=self.position, npix=(2, 2), axes=[axis]) time = Time(55555, format="mjd") if self.apply_irf["exposure"]: ref_unit = u.Unit("cm-2 s-1 MeV-1") else: ref_unit = u.Unit("") obt_unit = self.spectral_model(axis.center).unit if self.spatial_model: obt_unit = obt_unit * self.spatial_model.evaluate_geom(geom).unit ref_unit = ref_unit / u.sr if self.temporal_model: if u.Quantity(self.temporal_model(time)).unit.is_equivalent( self.spectral_model(axis.center).unit ): obt_unit = ( ( self.temporal_model(time) * self.spatial_model.evaluate_geom(geom).unit ) .to(obt_unit.to_string()) .unit ) else: obt_unit = obt_unit * u.Quantity(self.temporal_model(time)).unit if not obt_unit.is_equivalent(ref_unit): raise ValueError( f"SkyModel unit {obt_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): """Spatial model as a `~gammapy.modeling.models.SpatialModel` object.""" 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): """Spectral model as a `~gammapy.modeling.models.SpectralModel` object.""" 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): """Temporal model as a `~gammapy.modeling.models.TemporalModel` object.""" 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): """Position as a `~astropy.coordinates.SkyCoord`.""" return getattr(self.spatial_model, "position", None) @property def position_lonlat(self): """Spatial model center position `(lon, lat)` in radians and frame of the model.""" return getattr(self.spatial_model, "position_lonlat", None) @property def evaluation_bin_size_min(self): """Minimal spatial bin size for spatial model evaluation.""" if ( self.spatial_model is not None and self.spatial_model.evaluation_bin_size_min is not None ): return self.spatial_model.evaluation_bin_size_min else: return None @property def evaluation_radius(self): """Evaluation radius as an `~astropy.coordinates.Angle`.""" return self.spatial_model.evaluation_radius @property def evaluation_region(self): """Evaluation region as an `~astropy.coordinates.Angle`.""" return self.spatial_model.evaluation_region @property def frame(self): return self.spatial_model.frame def __add__(self, other): if isinstance(other, (Models, list)): return Models([self, *other]) elif isinstance(other, (SkyModel, TemplateNPredModel)): return Models([self, other]) else: raise TypeError(f"Invalid type: {other!r}") def __radd__(self, model): return self.__add__(model)
[docs] def __call__(self, lon, lat, energy, time=None): return self.evaluate(lon, lat, energy, time)
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})" )
[docs] def contributes(self, mask, margin="0 deg"): """Check if a sky model contributes within a mask map. Parameters ---------- mask : `~gammapy.maps.WcsNDMap` of boolean type Map containing a boolean mask. margin : `~astropy.units.Quantity` Add a margin in degree to the source evaluation radius. Used to take into account PSF width. Returns ------- models : `DatasetModels` Selected models contributing inside the region where mask is True. """ from gammapy.datasets.evaluator import CUTOUT_MARGIN margin = u.Quantity(margin) if not mask.geom.is_image: mask = mask.reduce_over_axes(func=np.logical_or) if mask.geom.is_region and mask.geom.region is not None: if mask.geom.is_all_point_sky_regions: return True geom = mask.geom.to_wcs_geom() mask = geom.region_mask([mask.geom.region]) try: mask_cutout = mask.cutout( position=self.position, width=(2 * self.evaluation_radius + CUTOUT_MARGIN + margin), ) contributes = np.any(mask_cutout.data) except (NoOverlapError, ValueError): contributes = False return contributes
[docs] 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`, optional Time coordinate. Default is None. 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
[docs] def evaluate_geom(self, geom, gti=None): """Evaluate model on `~gammapy.maps.Geom`.""" coords = geom.get_coord(sparse=True) value = self.spectral_model(coords["energy_true"]) additional_axes = set(coords.axis_names) - { "lon", "lat", "energy_true", "time", } for axis in additional_axes: value = value * np.ones_like(coords[axis]) if self.spatial_model: value = value * self.spatial_model.evaluate_geom(geom) if self.temporal_model: value = self._compute_time_integral(value, geom, gti) return value
[docs] def integrate_geom(self, geom, gti=None, oversampling_factor=None): """Integrate model on `~gammapy.maps.Geom`. See `~gammapy.modeling.models.SpatialModel.integrate_geom` and `~gammapy.modeling.models.SpectralModel.integral`. Parameters ---------- geom : `Geom` or `~gammapy.maps.RegionGeom` Map geometry. gti : `GTI`, optional GIT table. Default is None. oversampling_factor : int, optional The oversampling factor to use for spatial integration. Default is None: the factor is estimated from the model minimal bin size. Returns ------- flux : `Map` Predicted flux map. """ energy = geom.axes["energy_true"].edges shape = len(geom.data_shape) * [ 1, ] shape[geom.axes.index_data("energy_true")] = -1 value = self.spectral_model.integral( energy[:-1], energy[1:], ).reshape(shape) if self.spatial_model: value = ( value * self.spatial_model.integrate_geom( geom, oversampling_factor=oversampling_factor ).quantity ) if self.temporal_model: value = self._compute_time_integral(value, geom, gti) value = value * np.ones(geom.data_shape) return Map.from_geom(geom=geom, data=value.value, unit=value.unit)
def _compute_time_integral(self, value, geom, gti): """Multiply input value with time integral for the given geometry and GTI.""" if "time" in geom.axes.names: if geom.axes.names[-1] != "time": raise ValueError( "Incorrect axis order. The time axis must be the last axis" ) time_axis = geom.axes["time"] temp_eval = np.ones(time_axis.nbin) for idx in range(time_axis.nbin): if gti is None: t1, t2 = time_axis.time_min[idx], time_axis.time_max[idx] else: gti_in_bin = gti.select_time( time_interval=[ time_axis.time_min[idx], time_axis.time_max[idx], ] ) t1, t2 = gti_in_bin.time_start, gti_in_bin.time_stop integral = self.temporal_model.integral(t1, t2) temp_eval[idx] = np.sum(integral) value = (value.T * temp_eval).T else: if gti is not None: integral = self.temporal_model.integral(gti.time_start, gti.time_stop) value = value * np.sum(integral) return value
[docs] def copy(self, name=None, copy_data=False, **kwargs): """Copy sky model. Parameters ---------- name : str, optional Assign a new name to the copied model. Default is None. copy_data : bool, optional Copy the data arrays attached to models. Default is False. **kwargs : dict Keyword arguments forwarded to `SkyModel`. Returns ------- model : `SkyModel` Copied sky model. """ if self.spatial_model is not None: spatial_model = self.spatial_model.copy(copy_data=copy_data) 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) kwargs.setdefault("covariance_data", self.covariance.data.copy()) return self.__class__(**kwargs)
[docs] def to_dict(self, full_output=False): """Create dictionary for YAML serilisation.""" data = {} data["name"] = self.name data["type"] = self.tag[0] 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 data.update(self.spectral_model.to_dict(full_output)) if self.spatial_model is not None: data.update(self.spatial_model.to_dict(full_output)) if self.temporal_model is not None: data.update(self.temporal_model.to_dict(full_output)) return data
[docs] @classmethod def from_dict(cls, data, **kwargs): """Create SkyModel from dictionary.""" 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({"spectral": 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": 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": 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_ = f"{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_ += "\t" + "\n\t".join(lines[:-1]) str_ += "\n\n" return str_.expandtabs(tabsize=2)
[docs] @classmethod def create(cls, spectral_model, spatial_model=None, temporal_model=None, **kwargs): """Create a model instance. Parameters ---------- spectral_model : str Tag to create spectral model. spatial_model : str, optional Tag to create spatial model. Default is None. temporal_model : str, optional Tag to create temporal model. Default is None. **kwargs : dict Keyword arguments passed to `SkyModel`. Returns ------- model : SkyModel Sky model. """ spectral_model = Model.create(spectral_model, model_type="spectral") if spatial_model: spatial_model = Model.create(spatial_model, model_type="spatial") if temporal_model: temporal_model = Model.create(temporal_model, model_type="temporal") return cls( spectral_model=spectral_model, spatial_model=spatial_model, temporal_model=temporal_model, **kwargs, )
[docs] def freeze(self, model_type=None): """Freeze parameters depending on model type. Parameters ---------- model_type : {None, "spatial", "spectral", "temporal"} Freeze all parameters or only spatial/spectral/temporal. Default is None, such that all parameters are frozen. """ if model_type is None: self.parameters.freeze_all() else: model = getattr(self, f"{model_type}_model") model.freeze()
[docs] def unfreeze(self, model_type=None): """Restore parameters frozen status to default depending on model type. Parameters ---------- model_type : {None, "spatial", "spectral", "temporal"} Restore frozen status to default for all parameters or only spatial/spectral/temporal. Default is None, such that all parameters are restored to default frozen status. """ if model_type is None: for model_type in ["spectral", "spatial", "temporal"]: self.unfreeze(model_type) else: model = getattr(self, f"{model_type}_model") if model: model.unfreeze()
[docs]class FoVBackgroundModel(ModelBase): """Field of view background model. The background model holds the correction parameters applied to the instrumental background attached to a `MapDataset` or `SpectrumDataset`. Parameters ---------- spectral_model : `~gammapy.modeling.models.SpectralModel` Normalized spectral model. Default is `~gammapy.modeling.models.PowerLawNormSpectralModel` dataset_name : str Dataset name. spatial_model : `~gammapy.modeling.models.SpatialModel` Unitless Spatial model (unit is dropped on evaluation if defined). Default is None. """ tag = ["FoVBackgroundModel", "fov-bkg"] def __init__( self, spectral_model=None, dataset_name=None, spatial_model=None, covariance_data=None, ): if dataset_name is None: raise ValueError("Dataset name a is required argument") self.datasets_names = [dataset_name] if spectral_model is None: spectral_model = PowerLawNormSpectralModel() if not spectral_model.is_norm_spectral_model: raise ValueError("A norm spectral model is required.") self._spatial_model = spatial_model self._spectral_model = spectral_model super().__init__(covariance_data=covariance_data)
[docs] @staticmethod def contributes(*args, **kwargs): """FoV background models always contribute.""" return True
@property def spectral_model(self): """Spectral norm model.""" return self._spectral_model @property def spatial_model(self): """Spatial norm model.""" return self._spatial_model @property def name(self): """Model name.""" return self.datasets_names[0] + "-bkg" @property def parameters(self): """Model parameters.""" parameters = [] parameters.append(self.spectral_model.parameters) if self.spatial_model is not None: parameters.append(self.spatial_model.parameters) return Parameters.from_stack(parameters) def __str__(self): str_ = f"{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: str_ += "\t{:26}: {}\n".format( "Spatial model type", self.spatial_model.__class__.__name__ ) str_ += "\tParameters:\n" info = _get_parameters_str(self.parameters) lines = info.split("\n") str_ += "\t" + "\n\t".join(lines[:-1]) str_ += "\n\n" return str_.expandtabs(tabsize=2)
[docs] def evaluate_geom(self, geom): """Evaluate map.""" coords = geom.get_coord(sparse=True) return self.evaluate(**coords._data)
[docs] def evaluate(self, energy, lon=None, lat=None): """Evaluate model.""" value = self.spectral_model(energy) if self.spatial_model is not None: if lon is not None and lat is not None: if self.spatial_model.is_energy_dependent: return self.spatial_model(lon, lat, energy).value * value else: return self.spatial_model(lon, lat).value * value else: raise ValueError( "lon and lat are required if a spatial model is defined" ) else: return value
[docs] def copy(self, name=None, copy_data=False, **kwargs): """Copy the `FoVBackgroundModel` instance. Parameters ---------- name : str, optional Ignored, present for API compatibility. Default is None. copy_data : bool, optional Ignored, present for API compatibility. Default is False. **kwargs : dict Keyword arguments forwarded to `FoVBackgroundModel`. Returns ------- model : `FoVBackgroundModel` Copied FoV background model. """ kwargs.setdefault("spectral_model", self.spectral_model.copy()) kwargs.setdefault("dataset_name", self.datasets_names[0]) kwargs.setdefault("covariance_data", self.covariance.data.copy()) if self.spatial_model is not None: kwargs.setdefault("spatial_model", self.spatial_model.copy()) return self.__class__(**kwargs)
[docs] def to_dict(self, full_output=False): data = {} data["type"] = self.tag[0] data["datasets_names"] = self.datasets_names data.update(self.spectral_model.to_dict(full_output=full_output)) if self.spatial_model is not None: data.update(self.spatial_model.to_dict(full_output)) return data
[docs] @classmethod def from_dict(cls, data, **kwargs): """Create model from dictionary. Parameters ---------- data : dict Data dictionary. """ from gammapy.modeling.models import ( SPATIAL_MODEL_REGISTRY, SPECTRAL_MODEL_REGISTRY, ) spectral_data = data.get("spectral") if spectral_data is not None: model_class = SPECTRAL_MODEL_REGISTRY.get_cls(spectral_data["type"]) spectral_model = model_class.from_dict({"spectral": spectral_data}) else: spectral_model = None 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": spatial_data}) else: spatial_model = None datasets_names = data.get("datasets_names") if datasets_names is None: raise ValueError("FoVBackgroundModel must define a dataset name") if len(datasets_names) > 1: raise ValueError("FoVBackgroundModel can only be assigned to one dataset") return cls( spatial_model=spatial_model, spectral_model=spectral_model, dataset_name=datasets_names[0], )
[docs] def reset_to_default(self): """Reset parameter values to default.""" values = self.spectral_model.default_parameters.value self.spectral_model.parameters.value = values
[docs] def freeze(self, model_type="spectral"): """Freeze model parameters.""" if model_type is None or model_type == "spectral": self._spectral_model.freeze()
[docs] def unfreeze(self, model_type="spectral"): """Restore parameters frozen status to default.""" if model_type is None or model_type == "spectral": self._spectral_model.unfreeze()
[docs]class TemplateNPredModel(ModelBase): """Background model. Create a new map by a tilt and normalization on the available map. Parameters ---------- map : `~gammapy.maps.Map` Background model map. spectral_model : `~gammapy.modeling.models.SpectralModel` Normalized spectral model. Default is `~gammapy.modeling.models.PowerLawNormSpectralModel`. copy_data : bool Create a deepcopy of the map data or directly use the original. Default is True. Use False to save memory in case of large maps. spatial_model : `~gammapy.modeling.models.SpatialModel` Unitless Spatial model (unit is dropped on evaluation if defined). Default is None. """ tag = "TemplateNPredModel" map = LazyFitsData(cache=True) def __init__( self, map, spectral_model=None, name=None, filename=None, datasets_names=None, copy_data=True, spatial_model=None, covariance_data=None, ): if isinstance(map, Map): axis = map.geom.axes["energy"] if axis.node_type != "edges": raise ValueError( 'Need an integrated map, energy axis node_type="edges"' ) if copy_data: self.map = map.copy() else: self.map = map self._name = make_name(name) self.filename = filename if spectral_model is None: spectral_model = PowerLawNormSpectralModel() spectral_model.tilt.frozen = True self.spatial_model = spatial_model self.spectral_model = spectral_model if isinstance(datasets_names, str): datasets_names = [datasets_names] if isinstance(datasets_names, list): if len(datasets_names) != 1: raise ValueError( "Currently background models can only be assigned to one dataset." ) self.datasets_names = datasets_names super().__init__(covariance_data=covariance_data)
[docs] def copy(self, name=None, copy_data=False, **kwargs): """Copy template npred model. Parameters ---------- name : str, optional Assign a new name to the copied model. Default is None. copy_data : bool, optional Copy the data arrays attached to models. Default is False. **kwargs : dict Keyword arguments forwarded to `TemplateNPredModel`. Returns ------- model : `TemplateNPredModel` Copied template npred model. """ name = make_name(name) kwargs.setdefault("map", self.map) kwargs.setdefault("spectral_model", self.spectral_model.copy()) kwargs.setdefault("filename", self.filename) kwargs.setdefault("datasets_names", self.datasets_names) kwargs.setdefault("covariance_data", self.covariance.data.copy()) return self.__class__(name=name, copy_data=copy_data, **kwargs)
@property def name(self): return self._name @property def energy_center(self): """True energy axis bin centers as a `~astropy.units.Quantity`.""" energy_axis = self.map.geom.axes["energy"] energy = energy_axis.center return energy[:, np.newaxis, np.newaxis] @property def spectral_model(self): """Spectral model as a `~gammapy.modeling.models.SpectralModel` object.""" 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 parameters(self): parameters = [] parameters.append(self.spectral_model.parameters) return Parameters.from_stack(parameters)
[docs] def evaluate(self): """Evaluate background model. Returns ------- background_map : `~gammapy.maps.Map` Background evaluated on the Map. """ value = self.spectral_model(self.energy_center).value back_values = self.map.data * value if self.spatial_model is not None: value = self.spatial_model.evaluate_geom(self.map.geom).value back_values *= value return self.map.copy(data=back_values)
[docs] def to_dict(self, full_output=False): data = {} data["name"] = self.name data["type"] = self.tag if self.spatial_model is not None: data["spatial"] = self.spatial_model.to_dict(full_output)["spatial"] data["spectral"] = self.spectral_model.to_dict(full_output)["spectral"] if self.filename is not None: data["filename"] = self.filename if self.datasets_names is not None: data["datasets_names"] = self.datasets_names return data
[docs] def write(self, overwrite=False): """ Write the map. Parameters ---------- overwrite: bool, optional Overwrite existing file. Default is False, which will raise a warning if the template file exists already. """ if self.filename is None: raise IOError("Missing filename") elif os.path.isfile(make_path(self.filename)) and not overwrite: log.warning("Template file already exits, and overwrite is False") else: self.map.write(self.filename, overwrite=overwrite)
[docs] @classmethod def from_dict(cls, data, **kwargs): from gammapy.modeling.models import ( SPATIAL_MODEL_REGISTRY, SPECTRAL_MODEL_REGISTRY, ) spectral_data = data.get("spectral") if spectral_data is not None: model_class = SPECTRAL_MODEL_REGISTRY.get_cls(spectral_data["type"]) spectral_model = model_class.from_dict({"spectral": spectral_data}) else: spectral_model = None 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": spatial_data}) else: spatial_model = None if "filename" in data: bkg_map = Map.read(data["filename"]) else: raise IOError("Missing filename") return cls( map=bkg_map, spatial_model=spatial_model, spectral_model=spectral_model, name=data["name"], datasets_names=data.get("datasets_names"), filename=data.get("filename"), )
[docs] def cutout(self, position, width, mode="trim", name=None): """Cutout background model. Parameters ---------- position : `~astropy.coordinates.SkyCoord` Center position of the cutout region. width : tuple of `~astropy.coordinates.Angle` Angular sizes of the region in (lon, lat) in that specific order. If only one value is passed, a square region is extracted. mode : {'trim', 'partial', 'strict'} Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`. Default is "trim". name : str, optional Name of the returned background model. Default is None. Returns ------- cutout : `TemplateNPredModel` Cutout background model. """ cutout_kwargs = {"position": position, "width": width, "mode": mode} bkg_map = self.map.cutout(**cutout_kwargs) spectral_model = self.spectral_model.copy() return self.__class__(bkg_map, spectral_model=spectral_model, name=name)
[docs] def stack(self, other, weights=None, nan_to_num=True): """Stack background model in place. Stacking the background model resets the current parameters values. Parameters ---------- other : `TemplateNPredModel` Other background model. weights : float, optional Weights. Default is None. nan_to_num: bool, optional Non-finite values are replaced by zero if True. Default is True. """ bkg = self.evaluate() if nan_to_num: bkg.data[~np.isfinite(bkg.data)] = 0 other_bkg = other.evaluate() bkg.stack(other_bkg, weights=weights, nan_to_num=nan_to_num) self.map = bkg # reset parameter values self.spectral_model.norm.value = 1 self.spectral_model.tilt.value = 0
[docs] def slice_by_energy(self, energy_min=None, energy_max=None, name=None): """Select and slice model template in energy range Parameters ---------- energy_min, energy_max : `~astropy.units.Quantity` Energy bounds of the slice. Default is None. name : str Name of the sliced model. Default is None. Returns ------- model : `TemplateNpredModel` Sliced Model. """ name = make_name(name) energy_axis = self.map._geom.axes["energy"] if energy_min is None: energy_min = energy_axis.bounds[0] if energy_max is None: energy_max = energy_axis.bounds[1] if name is None: name = self.name energy_min, energy_max = u.Quantity(energy_min), u.Quantity(energy_max) group = energy_axis.group_table(edges=[energy_min, energy_max]) is_normal = group["bin_type"] == "normal " group = group[is_normal] slices = { "energy": slice(int(group["idx_min"][0]), int(group["idx_max"][0]) + 1) } model = self.copy(name=name) model.map = model.map.slice_by_idx(slices=slices) return model
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_ += "\tParameters:\n" info = _get_parameters_str(self.parameters) lines = info.split("\n") str_ += "\t" + "\n\t".join(lines[:-1]) str_ += "\n\n" return str_.expandtabs(tabsize=2) @property def position(self): """Position as a `~astropy.coordinates.SkyCoord`.""" return self.map.geom.center_skydir @property def evaluation_radius(self): """Evaluation radius as a `~astropy.coordinates.Angle`.""" return np.max(self.map.geom.width) / 2.0
[docs] def freeze(self, model_type="spectral"): """Freeze model parameters.""" if model_type is None or model_type == "spectral": self._spectral_model.freeze()
[docs] def unfreeze(self, model_type="spectral"): """Restore parameters frozen status to default.""" if model_type is None or model_type == "spectral": self._spectral_model.unfreeze()
[docs]def create_fermi_isotropic_diffuse_model(filename, **kwargs): """Read Fermi isotropic diffuse model. See `LAT Background models <https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html>`__ # noqa: E501 Parameters ---------- filename : str Filename. kwargs : dict Keyword arguments forwarded to `TemplateSpectralModel`. Returns ------- diffuse_model : `SkyModel` Fermi isotropic diffuse sky model. """ vals = np.loadtxt(make_path(filename)) energy = u.Quantity(vals[:, 0], "MeV", copy=False) values = u.Quantity(vals[:, 1], "MeV-1 s-1 cm-2", copy=False) kwargs.setdefault("interp_kwargs", {"fill_value": None}) spatial_model = ConstantSpatialModel() spectral_model = ( TemplateSpectralModel(energy=energy, values=values, **kwargs) * PowerLawNormSpectralModel() ) return SkyModel( spatial_model=spatial_model, spectral_model=spectral_model, name="fermi-diffuse-iso", apply_irf={"psf": False, "exposure": True, "edisp": True}, )