Source code for gammapy.modeling.models.temporal

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Time-dependent models."""

import logging
import numpy as np
import scipy.interpolate
from astropy import units as u
from astropy.io import fits
from astropy.table import Table
from astropy.time import Time
from gammapy.maps import MapAxis, RegionNDMap, TimeMapAxis
from astropy.utils import lazyproperty
from gammapy.modeling import Parameter
from gammapy.utils.compat import COPY_IF_NEEDED
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, time_ref_to_dict
from .core import ModelBase, _build_parameters_from_dict

__all__ = [
    "ConstantTemporalModel",
    "ExpDecayTemporalModel",
    "GaussianTemporalModel",
    "GeneralizedGaussianTemporalModel",
    "LightCurveTemplateTemporalModel",
    "LinearTemporalModel",
    "PowerLawTemporalModel",
    "SineTemporalModel",
    "TemplatePhaseCurveTemporalModel",
    "TemporalModel",
]

log = logging.getLogger(__name__)


# TODO: make this a small ABC to define a uniform interface.

[docs]
class TemporalModel(ModelBase):
    """Temporal model base class.

    Evaluates on `~astropy.time.Time` objects.
    """

    _type = "temporal"

    def __init__(self, **kwargs):
        scale = kwargs.pop("scale", "utc")

        if scale not in Time.SCALES:
            raise ValueError(
                f"{scale} is not a valid time scale. Choose from {Time.SCALES}"
            )

        self.scale = scale
        super().__init__(**kwargs)


[docs]
    def __call__(self, time, energy=None):
        """Evaluate model.

        Parameters
        ----------
        time : `~astropy.time.Time`
            Time object.
        energy : `~astropy.units.Quantity`, optional
            Energy. Default is None.

        Returns
        -------
        values : `~astropy.units.Quantity`
            Model values.
        """
        kwargs = {par.name: par.quantity for par in self.parameters}

        if energy is not None:
            kwargs["energy"] = energy

        time = Time(time, scale=self.scale).mjd * u.d
        return self.evaluate(time, **kwargs)


    @property
    def type(self):
        return self._type

    @property
    def is_energy_dependent(self):
        return False

    @property
    def reference_time(self):
        """Reference time in MJD."""
        return Time(self.t_ref.value, format="mjd", scale=self.scale)

    @reference_time.setter
    def reference_time(self, t_ref):
        """Reference time."""
        if not isinstance(t_ref, Time):
            raise TypeError(f"{t_ref} is not a {Time} object")
        self.t_ref.value = Time(t_ref, scale=self.scale).mjd


[docs]
    def to_dict(self, full_output=False):
        """Create dictionary for YAML serilisation."""
        data = super().to_dict(full_output)
        data["temporal"]["scale"] = self.scale
        return data



[docs]
    @classmethod
    def from_dict(cls, data, **kwargs):
        """Create a temporal model from a dictionary.

        Parameters
        ----------
        data : dict
            Dictionary containing the model parameters.
        **kwargs : dict
            Keyword arguments passed to `~TemporalModel.from_parameters`.
        """
        kwargs = kwargs or {}
        temporal_data = data.get("temporal", data)
        if "scale" in temporal_data:
            kwargs["scale"] = temporal_data["scale"]
        return super().from_dict(data, **kwargs)



[docs]
    @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
        return np.sum(diff).to(u.day)



[docs]
    def plot(self, time_range, ax=None, n_points=100, **kwargs):
        """
        Plot the temporal model.

        Parameters
        ----------
        time_range : `~astropy.time.Time`
            Times to plot the model.
        ax : `~matplotlib.axes.Axes`, optional
            Axis to plot on.
        n_points : int
            Number of bins to plot model. Default is 100.
        **kwargs : dict
            Keywords forwarded to `~matplotlib.pyplot.errorbar`.

        Returns
        -------
        ax : `~matplotlib.axes.Axes`, optional
            Matplotlib axes.
        """
        time_min, time_max = time_range
        time_axis = TimeMapAxis.from_time_bounds(
            time_min=time_min, time_max=time_max, nbin=n_points
        )

        m = RegionNDMap.create(region=None, axes=[time_axis])
        kwargs.setdefault("marker", "None")
        kwargs.setdefault("ls", "-")
        kwargs.setdefault("xerr", None)
        m.quantity = self(time_axis.time_mid).to(u.one)
        ax = m.plot(ax=ax, **kwargs)
        ax.set_ylabel("Norm / A.U.")
        return ax



[docs]
    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`, optional
            Time step used for sampling of the temporal model. Default is 1 s.
        random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
            Defines random number generator initialisation.
            Passed to `~gammapy.utils.random.get_random_state`.
            Default is 0.

        Returns
        -------
        time : `~astropy.units.Quantity`
            Array with times of the sampled events.
        """
        t_min = Time(t_min, scale=self.scale)
        t_max = Time(t_max, scale=self.scale)
        t_delta = u.Quantity(t_delta)
        random_state = get_random_state(random_state)

        ontime = (t_max - t_min).to("s")
        n_step = (ontime / t_delta).to_value("").item()
        t_step = ontime / n_step

        indices = np.arange(n_step + 1)
        steps = indices * t_step
        t = Time(t_min + steps, format="mjd")

        pdf = self(t)

        sampler = InverseCDFSampler(pdf=pdf, random_state=random_state)
        time_pix = sampler.sample(n_events)[0]
        time = np.interp(time_pix, indices, steps)
        return t_min + time



[docs]
    def integral(self, t_min, t_max, oversampling_factor=100, **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.
        oversampling_factor : int, optional
            Oversampling factor to be used for numerical integration.
            Default is 100.

        Returns
        -------
        norm : float
            Integrated flux norm on the given time intervals.
        """
        t_values, steps = np.linspace(
            t_min.mjd, t_max.mjd, oversampling_factor, retstep=True, axis=-1
        )
        times = Time(t_values, format="mjd", scale=self.scale)
        values = self(times)
        integral = np.sum(values, axis=-1) * steps
        return integral / self.time_sum(t_min, t_max).to_value("d")





[docs]
class ConstantTemporalModel(TemporalModel):
    """Constant temporal model.

    For more information see :ref:`constant-temporal-model`.
    """

    tag = ["ConstantTemporalModel", "const"]


[docs]
    @staticmethod
    def evaluate(time):
        """Evaluate at given times."""
        return np.ones(time.shape) * u.one



[docs]
    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)





[docs]
class LinearTemporalModel(TemporalModel):
    """Temporal model with a linear variation.

    For more information see :ref:`linear-temporal-model`.

    Parameters
    ----------
    alpha : float
        Constant term of the baseline flux.
        Default is 1.
    beta : `~astropy.units.Quantity`
        Time variation coefficient of the flux.
        Default is 0.
    t_ref : `~astropy.units.Quantity`
        The reference time in mjd.
        Frozen per default, at 2000-01-01.
    """

    tag = ["LinearTemporalModel", "linear"]

    alpha = Parameter("alpha", 1.0, frozen=False)
    beta = Parameter("beta", 0.0, unit="d-1", frozen=False)
    _t_ref_default = Time("2000-01-01")
    t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=True)


[docs]
    @staticmethod
    def evaluate(time, alpha, beta, t_ref):
        """Evaluate at given times."""
        return alpha + beta * (time - t_ref)



[docs]
    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
        alpha = pars["alpha"]
        beta = pars["beta"].quantity
        t_ref = self.reference_time
        value = alpha * (t_max - t_min) + beta / 2.0 * (
            (t_max - t_ref) * (t_max - t_ref) - (t_min - t_ref) * (t_min - t_ref)
        )
        return value / self.time_sum(t_min, t_max)





[docs]
class ExpDecayTemporalModel(TemporalModel):
    r"""Temporal model with an exponential decay.

    For more information see :ref:`expdecay-temporal-model`.

    Parameters
    ----------
    t0 : `~astropy.units.Quantity`
        Decay timescale. Default is 1 day.
    t_ref : `~astropy.units.Quantity`
        The reference time in mjd. Frozen per default, at 2000-01-01.
    """

    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)


[docs]
    @staticmethod
    def evaluate(time, t0, t_ref):
        """Evaluate at given times."""
        return np.exp(-(time - t_ref) / t0)



[docs]
    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 = self.reference_time
        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)





[docs]
class GaussianTemporalModel(TemporalModel):
    r"""A Gaussian temporal profile.

    For more information see :ref:`gaussian-temporal-model`.

    Parameters
    ----------
    t_ref : `~astropy.units.Quantity`
        The reference time in mjd at the peak.
        Default is 2000-01-01.
    sigma : `~astropy.units.Quantity`
        Width of the gaussian profile.
        Default is 1 day.
    """

    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)


[docs]
    @staticmethod
    def evaluate(time, t_ref, sigma):
        return np.exp(-((time - t_ref) ** 2) / (2 * sigma**2))



[docs]
    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 = self.reference_time
        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)





[docs]
class GeneralizedGaussianTemporalModel(TemporalModel):
    r"""A generalized Gaussian temporal profile.

    For more information see :ref:`generalized-gaussian-temporal-model`.

    Parameters
    ----------
    t_ref : `~astropy.units.Quantity`
        The time of the pulse's maximum intensity.
        Default is 2000-01-01.
    t_rise : `~astropy.units.Quantity`
        Rise time constant.
        Default is 1 day.
    t_decay : `~astropy.units.Quantity`
        Decay time constant.
        Default is 1 day.
    eta : `~astropy.units.Quantity`
        Inverse pulse sharpness -> higher values implies a more peaked pulse.
        Default is 1/2.

    """

    tag = ["GeneralizedGaussianTemporalModel", "gengauss"]

    _t_ref_default = Time("2000-01-01")
    t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=False)
    t_rise = Parameter("t_rise", "1d", frozen=False)
    t_decay = Parameter("t_decay", "1d", frozen=False)
    eta = Parameter("eta", 1 / 2, unit="", frozen=False)


[docs]
    @staticmethod
    def evaluate(time, t_ref, t_rise, t_decay, eta):
        val_rise = np.exp(
            -0.5
            * (np.abs(u.Quantity(time - t_ref, "d")) ** (1 / eta))
            / (t_rise ** (1 / eta))
        )
        val_decay = np.exp(
            -0.5
            * (np.abs(u.Quantity(time - t_ref, "d")) ** (1 / eta))
            / (t_decay ** (1 / eta))
        )
        val = np.where(time < t_ref, val_rise, val_decay)
        return val





[docs]
class LightCurveTemplateTemporalModel(TemporalModel):
    """Temporal light curve model.

    The lightcurve is given at specific times (and optionally energies) as a ``norm``
    It can be serialised either as an astropy table or a `~gammapy.maps.RegionNDMap`

    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, energy, norm)``
    values.

    When the temporal model is energy-dependent, the default interpolation scheme is
    linear with a log scale for the values. The interpolation method and scale values
    can be changed with the ``method`` and ``values_scale`` arguments.

    For more information see :ref:`LightCurve-temporal-model`.

    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)
    LightCurveTemplateTemporalModel model summary:
     Reference time: 59000.49919925926 MJD
     Start time: 58999.99919925926 MJD
     End time: 61862.99919925926 MJD
     Norm min: 0.01551196351647377
    Norm max: 1.0

    <BLANKLINE>

    Compute ``norm`` at a given time:

    >>> from astropy.time import Time
    >>> t = Time(59001.195, format="mjd")
    >>> light_curve.evaluate(t)
    <Quantity [0.02288737]>

    Compute mean ``norm`` in a given time interval:

    >>> import astropy.units as u
    >>> t_r = Time(59000.5, format='mjd')
    >>> t_min = t_r + [1, 4, 8] * u.d
    >>> t_max = t_r + [1.5, 6, 9] * u.d
    >>> light_curve.integral(t_min, t_max)
    <Quantity [0.00375698, 0.0143724 , 0.00688029]>
    """

    tag = ["LightCurveTemplateTemporalModel", "template"]

    _t_ref_default = Time("2000-01-01")
    t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=True)

    def __init__(self, map, t_ref=None, filename=None, method=None, values_scale=None):
        if (map.data < 0).any():
            log.warning("Map has negative values. Check and fix this!")

        self.map = map.copy()
        super().__init__()

        if t_ref:
            self.reference_time = t_ref

        self.filename = filename

        if method is None:
            method = "linear"

        if values_scale is None:
            if self.is_energy_dependent:
                values_scale = "log"
            else:
                values_scale = "lin"

        self.method = method
        self.values_scale = values_scale

    def __str__(self):
        start_time = self.t_ref.quantity + self.map.geom.axes["time"].edges[0]
        end_time = self.t_ref.quantity + self.map.geom.axes["time"].edges[-1]
        norm_min = np.nanmin(self.map.data)
        norm_max = np.nanmax(self.map.data)

        prnt = (
            f"{self.__class__.__name__} model summary:\n "
            f"Reference time: {self.t_ref.value} MJD \n "
            f"Start time: {start_time.value} MJD \n "
            f"End time: {end_time.value} MJD \n "
            f"Norm min: {norm_min} \n"
            f"Norm max: {norm_max}"
        )

        if self.is_energy_dependent:
            energy_min = self.map.geom.axes["energy"].center[0]
            energy_max = self.map.geom.axes["energy"].center[-1]
            prnt1 = f"Energy min: {energy_min} \n" f"Energy max: {energy_max} \n"
            prnt = prnt + prnt1

        return prnt

    @property
    def is_energy_dependent(self):
        """Whether the model is energy dependent."""
        return self.map.geom.has_energy_axis


[docs]
    @classmethod
    def from_table(cls, table, filename=None):
        """Create a template model from an astropy table.

        Parameters
        ----------
        table : `~astropy.table.Table`
            Table containing the template model.
        filename : str, optional
            Name of input file. Default is None.

        Returns
        -------
        model : `LightCurveTemplateTemporalModel`
            Light curve template model.
        """
        columns = [_.lower() for _ in table.colnames]
        if "time" not in columns:
            raise ValueError("A TIME column is necessary")

        t_ref = time_ref_from_dict(table.meta, scale="utc")
        nodes = table["TIME"]

        ax_unit = nodes.quantity.unit

        if not ax_unit.is_equivalent("d"):
            try:
                ax_unit = u.Unit(table.meta["TIMEUNIT"])
            except KeyError:
                raise ValueError("Time unit not found in the table")

        time_axis = MapAxis.from_nodes(nodes=nodes, name="time", unit=ax_unit)
        axes = [time_axis]
        m = RegionNDMap.create(region=None, axes=axes, data=table["NORM"])

        return cls(m, t_ref=t_ref, filename=filename)



[docs]
    @classmethod
    def read(cls, filename, format="table"):
        """Read a template model.

        Parameters
        ----------
        filename : str
            Name of file to read.
        format : {"table", "map"}
            Format of the input file.

        Returns
        -------
        model : `LightCurveTemplateTemporalModel`
            Light curve template model.
        """
        filename = str(make_path(filename))
        if format == "table":
            table = Table.read(filename)
            return cls.from_table(table, filename=filename)

        elif format == "map":
            with fits.open(filename) as hdulist:
                header = hdulist["SKYMAP_BANDS"].header
                t_ref = time_ref_from_dict(header)
                # TODO : Ugly hack to prevent creating a TimeMapAxis
                # By default, MapAxis.from_table tries to create a
                # TimeMapAxis, failing which, it creates a normal MapAxis.
                # This ugly hack forces the fail. We need a normal Axis to
                # have the evaluate method work
                hdulist["SKYMAP_BANDS"].header.pop("MJDREFI")
                m = RegionNDMap.from_hdulist(hdulist)
            return cls(m, t_ref=t_ref, filename=filename)

        else:
            raise ValueError(
                f"Not a valid format: '{format}', choose from: {'table', 'map'}"
            )



[docs]
    def to_table(self):
        """Convert model to an astropy table."""
        if self.is_energy_dependent:
            raise NotImplementedError("Not supported for energy dependent models")
        table = Table(
            data=[self.map.geom.axes["time"].center, self.map.quantity],
            names=["TIME", "NORM"],
            meta=time_ref_to_dict(self.reference_time, scale=self.scale),
        )
        return table



[docs]
    def write(self, filename, format="table", overwrite=False):
        """Write a model to disk as per the specified format.

        Parameters:
            filename : str
                Name of output file.
            format : {"table" or "map"}
                If format is "table", it is serialised as a `~astropy.table.Table`.
                If "map", then it is serialised as a `~gammapy.maps.RegionNDMap`.
                Default is "table".
            overwrite : bool, optional
                Overwrite existing file. Default is False.
        """
        if self.filename is None:
            raise IOError("Missing filename")

        if format == "table":
            table = self.to_table()
            table.write(filename, overwrite=overwrite)
        elif format == "map":
            # RegionNDMap.from_hdulist does not update the header
            hdulist = self.map.to_hdulist()
            hdulist["SKYMAP_BANDS"].header.update(
                time_ref_to_dict(self.reference_time, scale=self.scale)
            )
            hdulist.writeto(filename, overwrite=overwrite)
        else:
            raise ValueError("Not a valid format, choose from ['map', 'table']")



[docs]
    def evaluate(self, time, t_ref=None, energy=None):
        """Evaluate the model at given coordinates.

        Parameters
        ----------
        time: `~astropy.time.Time`
            Time.
        t_ref: `~gammapy.modeling.Parameter`, optional
            Reference time for the model. Default is None.
        energy: `~astropy.units.Quantity`, optional
            Energy. Default is None.

        Returns
        -------
        values : `~astropy.units.Quantity`
            Model values.
        """
        if t_ref is None:
            t_ref = self.reference_time

        t = (time - t_ref).to_value(self.map.geom.axes["time"].unit)
        coords = {"time": t}

        if self.is_energy_dependent:
            if energy is None:
                energy = self.map.geom.axes["energy"].center

            coords["energy"] = energy.reshape(-1, 1)

        val = self.map.interp_by_coord(
            coords, method=self.method, values_scale=self.values_scale
        )
        val = np.clip(val, 0, a_max=None)
        return u.Quantity(val, unit=self.map.unit, copy=COPY_IF_NEEDED)



[docs]
    def integral(self, t_min, t_max, oversampling_factor=100, **kwargs):
        if self.is_energy_dependent:
            raise NotImplementedError(
                "Integral not supported for energy dependent models"
            )

        return super().integral(t_min, t_max, oversampling_factor, **kwargs)



[docs]
    @classmethod
    def from_dict(cls, data):
        data = data["temporal"]
        filename = data["filename"]
        format = data.get("format", "table")
        return cls.read(filename, format)


    def _guess_format(self):
        if self.is_energy_dependent:
            format = "map"
        else:
            format = "table"
        log.info("Inferred format: " + format)
        return format


[docs]
    def to_dict(self, full_output=False, format=None):
        """Create dictionary for YAML serialisation."""
        data = super().to_dict(full_output)
        if format is None:
            format = self._guess_format()
        data["temporal"]["filename"] = self.filename
        data["temporal"]["format"] = format
        data["temporal"]["unit"] = str(self.map.unit)
        return data



[docs]
    def plot(self, time_range, ax=None, n_points=100, energy=None, **kwargs):
        """
        Plot the temporal model.

        Parameters
        ----------
        time_range : `~astropy.time.Time`
            Times to plot the model.
        ax : `~matplotlib.axes.Axes`, optional
            Axis to plot on. Default is None.
        n_points : int, optional
            Number of bins to plot model. Default is 100.
        energy : `~astropy.units.quantity`, optional
            Energies to compute the model at for energy dependent models. Default is None.
        **kwargs : dict
            Keywords forwarded to `~matplotlib.pyplot.errorbar`.
        Returns
        -------
        ax : `~matplotlib.axes.Axes`, optional
            Matplotlib axes.
        """
        if not self.is_energy_dependent:
            super().plot(time_range=time_range, ax=ax, n_points=n_points, **kwargs)
        else:
            time_min, time_max = Time(time_range, scale=self.scale)
            time_axis = TimeMapAxis.from_time_bounds(
                time_min=time_min, time_max=time_max, nbin=n_points
            )
            if energy is None:
                energy_axis = self.map.geom.axes["energy"]
            else:
                energy_axis = MapAxis.from_nodes(
                    nodes=energy, name="energy", interp="log"
                )

            m = RegionNDMap.create(region=None, axes=[time_axis, energy_axis])
            kwargs.setdefault("marker", "None")
            kwargs.setdefault("ls", "-")
            m.quantity = self.evaluate(
                time=time_axis.time_mid, energy=energy_axis.center
            )
            ax = m.plot(axis_name="time", ax=ax, **kwargs)
            ax.set_ylabel("Norm / A.U.")

            return ax, m





[docs]
class PowerLawTemporalModel(TemporalModel):
    """Temporal model with a Power Law decay.

    For more information see :ref:`powerlaw-temporal-model`.

    Parameters
    ----------
    alpha : float
        Decay time power. Default is 1.
    t_ref: `~astropy.units.Quantity`
        The reference time in mjd.
        Frozen by default, at 2000-01-01.
    t0: `~astropy.units.Quantity`
        The scaling time in mjd.
        Fixed by default, at 1 day.
    """

    tag = ["PowerLawTemporalModel", "powerlaw"]

    alpha = Parameter("alpha", 1.0, frozen=False)
    _t_ref_default = Time("2000-01-01")
    t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=True)
    t0 = Parameter("t0", "1 d", frozen=True)


[docs]
    @staticmethod
    def evaluate(time, alpha, t_ref, t0=1 * u.day):
        """Evaluate at given times."""
        return np.power((time - t_ref) / t0, alpha)



[docs]
    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
        alpha = pars["alpha"].quantity
        t0 = pars["t0"].quantity
        t_ref = self.reference_time
        if alpha != -1:
            value = self.evaluate(t_max, alpha + 1.0, t_ref, t0) - self.evaluate(
                t_min, alpha + 1.0, t_ref, t0
            )
            return t0 / (alpha + 1.0) * value / self.time_sum(t_min, t_max)
        else:
            value = np.log((t_max - t_ref) / (t_min - t_ref))
            return t0 * value / self.time_sum(t_min, t_max)





[docs]
class SineTemporalModel(TemporalModel):
    """Temporal model with a sinusoidal modulation.

    For more information see :ref:`sine-temporal-model`.

    Parameters
    ----------
    amp : float
        Amplitude of the sinusoidal function.
        Default is 1.
    t_ref: `~astropy.units.Quantity`
        The reference time in mjd.
        Default is 2000-01-01.
    omega: `~astropy.units.Quantity`
        Pulsation of the signal.
        Default is 1 rad/day.
    """

    tag = ["SineTemporalModel", "sinus"]

    amp = Parameter("amp", 1.0, frozen=False)
    omega = Parameter("omega", "1. rad/day", frozen=False)
    _t_ref_default = Time("2000-01-01")
    t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=False)


[docs]
    @staticmethod
    def evaluate(time, amp, omega, t_ref):
        """Evaluate at given times."""
        return 1.0 + amp * np.sin(omega * (time - t_ref))



[docs]
    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
        omega = pars["omega"].quantity.to_value("rad/day")
        amp = pars["amp"].value
        t_ref = self.reference_time

        value = (t_max - t_min).to_value(u.day) - amp / omega * (
            np.sin(omega * (t_max - t_ref).to_value(u.day))
            - np.sin(omega * (t_min - t_ref).to_value(u.day))
        )
        return value / self.time_sum(t_min, t_max).to_value(u.day)





[docs]
class TemplatePhaseCurveTemporalModel(TemporalModel):
    """Temporal phase curve model.

    A timing solution is used to compute the phase corresponding to time and
    a template phase curve is used to determine the associated ``norm``.

    The phasecurve is given as a table with columns ``phase`` 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 ``(phase, 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 'PHASE' vs 'NORM'.
    filename : str
        The name of the file containing the phase curve.
    t_ref : `~astropy.units.Quantity`
        The reference time in mjd.
        Default is 48442.5 mjd.
    phi_ref : `~astropy.units.Quantity`
        The phase at reference time.
        Default is 0.
    f0 : `~astropy.units.Quantity`
        The frequency at t_ref in s-1.
        Default is 29.946923 s-1.
    f1 : `~astropy.units.Quantity`
        The frequency derivative at t_ref in s-2.
        Default is 0 s-2.
    f2 : `~astropy.units.Quantity`
        The frequency second derivative at t_ref in s-3.
        Default is 0 s-3.
    """

    tag = ["TemplatePhaseCurveTemporalModel", "template-phase"]
    _t_ref_default = Time(48442.5, format="mjd")
    _phi_ref_default = 0
    _f0_default = 29.946923 * u.s**-1
    _f1_default = 0 * u.s**-2
    _f2_default = 0 * u.s**-3

    t_ref = Parameter("t_ref", _t_ref_default.mjd, unit="day", frozen=True)
    phi_ref = Parameter("phi_ref", _phi_ref_default, unit="", frozen=True)
    f0 = Parameter("f0", _f0_default, frozen=True)
    f1 = Parameter("f1", _f1_default, frozen=True)
    f2 = Parameter("f2", _f2_default, frozen=True)

    def __init__(self, table, filename=None, **kwargs):
        self.table = self._normalise_table(table)
        if filename is not None:
            filename = str(make_path(filename))
        self.filename = filename
        super().__init__(**kwargs)

    @staticmethod
    def _normalise_table(table):
        x = table["PHASE"].data
        y = table["NORM"].data

        interpolator = scipy.interpolate.InterpolatedUnivariateSpline(
            x, y, k=1, ext=2, bbox=[0.0, 1.0]
        )

        integral = interpolator.integral(0, 1)

        table["NORM"] *= 1 / integral
        return table


[docs]
    @classmethod
    def read(
        cls,
        path,
        t_ref=_t_ref_default.mjd * u.d,
        phi_ref=_phi_ref_default,
        f0=_f0_default,
        f1=_f1_default,
        f2=_f2_default,
    ):
        """Read phasecurve model table from FITS file.

        Beware : this does **not** read parameters.
        They will be set to defaults.

        Parameters
        ----------
        path : str or `~pathlib.Path`
            Filename with path.
        """
        filename = str(make_path(path))

        return cls(
            Table.read(filename),
            filename=filename,
            t_ref=t_ref,
            phi_ref=phi_ref,
            f0=f0,
            f1=f1,
            f2=f2,
        )


    @staticmethod
    def _time_to_phase(time, t_ref, phi_ref, f0, f1, f2):
        """Convert time to phase given timing solution parameters.

        Parameters
        ----------
        time : `~astropy.units.Quantity`
            The time at which to compute the phase.
        t_ref : `~astropy.units.Quantity`
            The reference time in mjd.
        phi_ref : `~astropy.units.Quantity`
            The phase at reference time.
            Default is 0.
        f0 : `~astropy.units.Quantity`
            The frequency at t_ref in s-1.
        f1 : `~astropy.units.Quantity`
            The frequency derivative at t_ref in s-2.
        f2 : `~astropy.units.Quantity`
            The frequency second derivative at t_ref in s-3.

        Returns
        -------
        phase : float
            Phase.
        period_number : int
            Number of period since t_ref.
        """
        delta_t = time - t_ref
        phase = (
            phi_ref + delta_t * (f0 + delta_t / 2.0 * (f1 + delta_t / 3 * f2))
        ).to_value("")

        period_number = np.floor(phase)
        phase -= period_number
        return phase, period_number


[docs]
    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):
        x = self.table["PHASE"].data
        y = self.table["NORM"].data

        return scipy.interpolate.InterpolatedUnivariateSpline(
            x, y, k=1, ext=2, bbox=[0.0, 1.0]
        )


[docs]
    def evaluate(self, time, t_ref, phi_ref, f0, f1, f2):
        phase, _ = self._time_to_phase(time, t_ref, phi_ref, f0, f1, f2)
        return self._interpolator(phase) * u.one



[docs]
    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.
        """
        kwargs = {par.name: par.quantity for par in self.parameters}
        ph_min, n_min = self._time_to_phase(t_min.mjd * u.d, **kwargs)
        ph_max, n_max = self._time_to_phase(t_max.mjd * u.d, **kwargs)

        # here we assume that the frequency does not change during the integration boundaries
        delta_t = (t_min - self.reference_time).to(u.d)
        frequency = self.f0.quantity + delta_t * (
            self.f1.quantity + delta_t * self.f2.quantity / 2
        )

        # Compute integral of one phase
        phase_integral = self._interpolator.antiderivative()(
            1
        ) - self._interpolator.antiderivative()(0)
        # Multiply by the total number of phases
        phase_integral *= n_max - n_min - 1

        # Compute integrals before first full phase and after the last full phase
        end_integral = self._interpolator.antiderivative()(
            ph_max
        ) - self._interpolator.antiderivative()(0)
        start_integral = self._interpolator.antiderivative()(
            1
        ) - self._interpolator.antiderivative()(ph_min)

        # Divide by Jacobian (here we neglect variations of frequency during the integration period)
        total = phase_integral + start_integral + end_integral
        # Normalize by total integration time
        n_period = (self.time_sum(t_min, t_max) * frequency).to("")
        if int(n_period) == 0:
            n_period = 1

        integral_norm = total / n_period

        return integral_norm



[docs]
    @classmethod
    def from_dict(cls, data):
        params = _build_parameters_from_dict(
            data["temporal"]["parameters"], cls.default_parameters
        )
        filename = data["temporal"]["filename"]
        kwargs = {par.name: par for par in params}
        return cls.read(filename, **kwargs)



[docs]
    def to_dict(self, full_output=False):
        """Create dictionary for YAML serialisation."""
        model_dict = super().to_dict()
        model_dict["temporal"]["filename"] = self.filename
        return model_dict



[docs]
    def plot_phasogram(self, ax=None, n_points=100, **kwargs):
        """
        Plot phasogram of the phase model.

        Parameters
        ----------
        ax : `~matplotlib.axes.Axes`, optional
            Axis to plot on. Default is None.
        n_points : int, optional
            Number of bins to plot model. Default is 100.
        **kwargs : dict
            Keywords forwarded to `~matplotlib.pyplot.errorbar`.

        Returns
        -------
        ax : `~matplotlib.axes.Axes`, optional
            Matplotlib axes.
        """
        phase_axis = MapAxis.from_bounds(0.0, 1, nbin=n_points, name="Phase", unit="")

        m = RegionNDMap.create(region=None, axes=[phase_axis])
        kwargs.setdefault("marker", "None")
        kwargs.setdefault("ls", "-")
        kwargs.setdefault("xerr", None)
        m.quantity = self._interpolator(phase_axis.center)
        ax = m.plot(ax=ax, **kwargs)
        ax.set_ylabel("Norm / A.U.")
        return ax



[docs]
    def sample_time(self, n_events, t_min, t_max, t_delta="1 s", random_state=0):
        """Sample arrival times of events.

        To fully cover the phase range, t_delta is the minimum between the input
        and product of the period at 0.5*(t_min + t_max) and the table bin size.

        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.
        """
        t_delta = u.Quantity(t_delta)

        # Determine period at the mid time
        t_mid = Time(t_min, scale=self.scale) + 0.5 * (t_max - t_min)
        delta_t = (t_mid - self.reference_time).to(u.d)
        frequency = self.f0.quantity + delta_t * (
            self.f1.quantity + delta_t * self.f2.quantity / 2
        )
        period = 1 / frequency

        # Take minimum time delta between user input and the period divided by the number of rows in the model table
        # this assumes that phase values are evenly spaced.
        t_delta = np.minimum(period / len(self.table), t_delta)

        return super().sample_time(n_events, t_min, t_max, t_delta, random_state)