Source code for gammapy.stats.variability

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import scipy.stats as stats
import astropy.units as u
from gammapy.utils.random import get_random_state

__all__ = [
    "compute_fvar",
    "compute_fpp",
    "compute_chisq",
    "compute_flux_doubling",
    "structure_function",
    "TimmerKonig_lightcurve_simulator",
    "discrete_correlation",
]



[docs]
def compute_fvar(flux, flux_err, axis=0):
    r"""Calculate the fractional excess variance.

    This method accesses the ``FLUX`` and ``FLUX_ERR`` columns
    from the lightcurve data.

    The fractional excess variance :math:`F_{var}`, an intrinsic
    variability estimator, is given by:

    .. math::
        F_{var} = \sqrt{ \frac{S^{2} - \bar{ \sigma^{2}}}{ \bar{x}^{2}}}

    It is the excess variance after accounting for the measurement errors
    on the light curve :math:`\sigma`. :math:`S` is the variance.

    It is important to note that the errors on the flux must be gaussian.

    Parameters
    ----------
    flux : `~astropy.units.Quantity`
        The measured fluxes.
    flux_err : `~astropy.units.Quantity`
        The error on measured fluxes.
    axis : int, optional
        Axis along which the excess variance is computed.
        Default is 0.

    Returns
    -------
    fvar, fvar_err : `~numpy.ndarray`
        Fractional excess variance.

    References
    ----------
    .. [Vaughan2003] "On characterizing the variability properties of X-ray light
       curves from active galaxies", Vaughan et al. (2003)
       https://ui.adsabs.harvard.edu/abs/2003MNRAS.345.1271V
    """
    flux_mean = np.nanmean(flux, axis=axis)
    n_points = np.count_nonzero(~np.isnan(flux), axis=axis)

    s_square = np.nansum((flux - flux_mean) ** 2, axis=axis) / (n_points - 1)
    sig_square = np.nansum(flux_err**2, axis=axis) / n_points
    fvar = np.sqrt(np.abs(s_square - sig_square)) / flux_mean

    sigxserr_a = np.sqrt(2 / n_points) * sig_square / flux_mean**2
    sigxserr_b = np.sqrt(sig_square / n_points) * (2 * fvar / flux_mean)
    sigxserr = np.sqrt(sigxserr_a**2 + sigxserr_b**2)
    fvar_err = sigxserr / (2 * fvar)

    return fvar, fvar_err




[docs]
def compute_fpp(flux, flux_err, axis=0):
    r"""Calculate the point-to-point excess variance.

    F_pp is a quantity strongly related to the fractional excess variance F_var
    implemented in `~gammapy.stats.compute_fvar`; F_pp probes the variability
    on a shorter timescale.

    For white noise, F_pp and F_var give the same value.
    However, for red noise, F_var will be larger
    than F_pp, as the variations will be larger on longer timescales.

    It is important to note that the errors on the flux must be Gaussian.

    Parameters
    ----------
    flux : `~astropy.units.Quantity`
        The measured fluxes.
    flux_err : `~astropy.units.Quantity`
        The error on measured fluxes.
    axis : int, optional
        Axis along which the excess variance is computed.
        Default is 0.

    Returns
    -------
    fpp, fpp_err : `~numpy.ndarray`
        Point-to-point excess variance.

    References
    ----------
    .. [Edelson2002] "X-Ray Spectral Variability and Rapid Variability
       of the Soft X-Ray Spectrum Seyfert 1 Galaxies
       Arakelian 564 and Ton S180", Edelson et al. (2002), equation 3,
       https://iopscience.iop.org/article/10.1086/323779
    """
    flux_mean = np.nanmean(flux, axis=axis)
    n_points = np.count_nonzero(~np.isnan(flux), axis=axis)
    flux = flux.swapaxes(0, axis).T

    s_square = np.nansum((flux[..., 1:] - flux[..., :-1]) ** 2, axis=-1) / (
        n_points.T - 1
    )
    sig_square = np.nansum(flux_err**2, axis=axis) / n_points
    fpp = np.sqrt(np.abs(s_square.T - sig_square)) / flux_mean

    sigxserr_a = np.sqrt(2 / n_points) * sig_square / flux_mean**2
    sigxserr_b = np.sqrt(sig_square / n_points) * (2 * fpp / flux_mean)
    sigxserr = np.sqrt(sigxserr_a**2 + sigxserr_b**2)
    fpp_err = sigxserr / (2 * fpp)

    return fpp, fpp_err




[docs]
def compute_chisq(flux):
    r"""Calculate the chi-square test for `LightCurve`.

    Chisquare test is a variability estimator. It computes
    deviations from the expected value here mean value.

    Parameters
    ----------
    flux : `~astropy.units.Quantity`
        The measured fluxes.

    Returns
    -------
    ChiSq, P-value : tuple of float or `~numpy.ndarray`
        Tuple of Chi-square and P-value.
    """
    yexp = np.mean(flux)
    yobs = flux.data
    chi2, pval = stats.chisquare(yobs, yexp)
    return chi2, pval




[docs]
def compute_flux_doubling(flux, flux_err, coords, axis=0):
    r"""Compute the minimum characteristic flux doubling and halving
    over a certain coordinate axis for a series of measurements.

    Computing the flux doubling can give the doubling time in a lightcurve
    displaying significant temporal variability, e.g. an AGN flare.

    The variable is computed as:

     .. math::
        doubling = min(\frac{t_(i+1)-t_i}{log_2{f_(i+1)/f_i}})

    where f_i and f_(i+1) are the fluxes measured at subsequent coordinates t_i and t_(i+1).
    The error is obtained by propagating the relative errors on the flux measures.

    Parameters
    ----------
    flux : `~astropy.units.Quantity`
        The measured fluxes.
    flux_err : `~astropy.units.Quantity`
        The error on measured fluxes.
    coords : `~astropy.units.Quantity`
        The coordinates at which the fluxes are measured.
    axis : int, optional
        Axis along which the value is computed.

    Returns
    -------
    doubling_dict : dict
        Dictionary containing the characteristic flux doubling, halving and errors,
        with coordinates at which they were found.
    """
    flux = np.atleast_2d(flux).swapaxes(0, axis).T
    flux_err = np.atleast_2d(flux_err).swapaxes(0, axis).T

    axes = np.diff(coords) / np.log2(flux[..., 1:] / flux[..., :-1])
    axes_err_1 = (
        np.diff(coords)
        * np.log(2)
        / flux[..., 1:]
        * np.log(flux[..., 1:] / flux[..., :-1]) ** 2
    )
    axes_err_2 = (
        np.diff(coords)
        * np.log(2)
        / flux[..., :-1]
        * np.log(flux[..., 1:] / flux[..., :-1]) ** 2
    )
    axes_err = np.sqrt(
        (flux_err[..., 1:] * axes_err_1) ** 2 + (flux_err[..., :-1] * axes_err_2) ** 2
    )

    imin = np.expand_dims(
        np.argmin(
            np.where(
                np.logical_and(np.isfinite(axes), axes > 0), axes, np.inf * coords.unit
            ),
            axis=-1,
        ),
        axis=-1,
    )
    imax = np.expand_dims(
        np.argmax(
            np.where(
                np.logical_and(np.isfinite(axes), axes < 0), axes, -np.inf * coords.unit
            ),
            axis=-1,
        ),
        axis=-1,
    )

    index = np.concatenate([imin, imax], axis=-1)
    coord = np.take_along_axis(coords, index.flatten(), axis=0).reshape(index.shape)

    doubling = np.take_along_axis(axes, index, axis=-1)
    doubling_err = np.take_along_axis(axes_err, index, axis=-1)

    doubling_dict = {
        "doubling": doubling.T[0],
        "doubling_err": doubling_err.T[0],
        "doubling_coord": coord.T[0],
        "halving": np.abs(doubling.T[1]),
        "halving_err": doubling_err.T[1],
        "halving_coord": coord.T[1],
    }

    return doubling_dict




[docs]
def structure_function(flux, flux_err, time, tdelta_precision=5):
    """Compute the discrete structure function for a variable source.

    Parameters
    ----------
    flux : `~astropy.units.Quantity`
        The measured fluxes.
    flux_err : `~astropy.units.Quantity`
        The error on measured fluxes.
    time : `~astropy.units.Quantity`
        The time coordinates at which the fluxes are measured.
    tdelta_precision : int, optional
        The number of decimal places to check to separate the time deltas. Default is 5.

    Returns
    -------
    sf, distances : `~numpy.ndarray`, `~astropy.units.Quantity`
        Discrete structure function and array of time distances.

    References
    ----------
    .. [Emmanoulopoulos2010] "On the use of structure functions to study blazar variability:
       caveats and problems", Emmanoulopoulos et al. (2010)
       https://academic.oup.com/mnras/article/404/2/931/968488
    """
    dist_matrix = (time[np.newaxis, :] - time[:, np.newaxis]).round(
        decimals=tdelta_precision
    )
    distances = np.unique(dist_matrix)
    distances = distances[distances > 0]
    shape = distances.shape + flux.shape[1:]
    factor = np.zeros(shape)
    norm = np.zeros(shape)

    for i, distance in enumerate(distances):
        indexes = np.array(np.where(dist_matrix == distance))
        for index in indexes.T:
            f = (flux[index[1], ...] - flux[index[0], ...]) ** 2
            w = (flux[index[1], ...] / flux_err[index[1], ...]) * (
                flux[index[0], ...] / flux_err[index[0], ...]
            )

            f = np.nan_to_num(f)
            w = np.nan_to_num(w)
            factor[i] = factor[i] + f * w
            norm[i] = norm[i] + w

    sf = factor / norm
    return sf, distances




[docs]
def discrete_correlation(flux1, flux_err1, flux2, flux_err2, time1, time2, tau, axis=0):
    """Compute the discrete correlation function for a variable source.

    Parameters
    ----------
    flux1, flux_err1: `~astropy.units.Quantity`
        The first set of measured fluxes and associated error.
    flux2, flux_err2 : `~astropy.units.Quantity`
        The second set of measured fluxes and associated error.
    time1, time2 : `~astropy.units.Quantity`
        The time coordinates at which the fluxes are measured.
    tau : `~astropy.units.Quantity`
        Size of the bins to compute the discrete correlation.
    axis : int, optional
        Axis along which the correlation is computed.
        Default is 0.

    Returns
    -------
    bincenters: `~astropy.units.Quantity`
        Array of discrete time bins.
    discrete_correlation: `~numpy.ndarray`
        Array of discrete correlation function values for each bin.
    discrete_correlation_err : `~numpy.ndarray`
        Error associated to the discrete correlation values.

    References
    ----------
    .. [Edelson1988] "THE DISCRETE CORRELATION FUNCTION: A NEW METHOD FOR ANALYZING
       UNEVENLY SAMPLED VARIABILITY DATA", Edelson et al. (1988)
       https://ui.adsabs.harvard.edu/abs/1988ApJ...333..646E/abstract
    """
    flux1 = np.rollaxis(flux1, axis, 0)
    flux2 = np.rollaxis(flux2, axis, 0)

    if np.squeeze(flux1).shape[1:] != np.squeeze(flux2).shape[1:]:
        raise ValueError(
            "flux1 and flux2 must have the same squeezed shape, apart from the chosen axis."
        )

    tau = tau.to(time1.unit)
    time2 = time2.to(time1.unit)

    mean1, mean2 = np.nanmean(flux1, axis=0), np.nanmean(flux2, axis=0)
    sigma1, sigma2 = np.nanstd(flux1, axis=0), np.nanstd(flux2, axis=0)

    udcf1 = (flux1 - mean1) / np.sqrt((sigma1**2 - np.nanmean(flux_err1, axis=0) ** 2))
    udcf2 = (flux2 - mean2) / np.sqrt((sigma2**2 - np.nanmean(flux_err2, axis=0) ** 2))

    udcf = np.empty(((flux1.shape[0],) + flux2.shape))
    dist = u.Quantity(np.empty(((flux1.shape[0], flux2.shape[0]))), unit=time1.unit)

    for i, x1 in enumerate(udcf1):
        for j, x2 in enumerate(udcf2):
            udcf[i, j, ...] = x1 * x2
            dist[i, j] = time1[i] - time2[j]

    maxfactor = np.floor(np.amax(dist) / tau).value + 1
    minfactor = np.floor(np.amin(dist) / tau).value

    bins = (
        np.linspace(
            minfactor, maxfactor, int(np.abs(maxfactor) + np.abs(minfactor) + 1)
        )
        * tau
    )

    bin_indices = np.digitize(dist, bins).flatten()

    udcf = np.reshape(udcf, (udcf.shape[0] * udcf.shape[1], -1))
    discrete_correlation = np.array(
        [np.nanmean(udcf[bin_indices == i], axis=0) for i in range(1, len(bins))]
    )

    discrete_correlation_err = []
    for i in range(1, len(bins)):
        terms = (discrete_correlation[i - 1] - udcf[bin_indices == i]) ** 2
        num = np.sqrt(np.nansum(terms, axis=0))
        den = len(udcf[bin_indices == i]) - 1
        discrete_correlation_err.append(num / den)

    bincenters = (bins[1:] + bins[:-1]) / 2

    return bincenters, discrete_correlation, np.array(discrete_correlation_err)




[docs]
def TimmerKonig_lightcurve_simulator(
    power_spectrum,
    npoints,
    spacing,
    nchunks=10,
    random_state="random-seed",
    power_spectrum_params=None,
    mean=0.0,
    std=1.0,
    poisson=False,
):
    """Implementation of the Timmer-Koenig algorithm to simulate a time series from a power spectrum.

    Parameters
    ----------
    power_spectrum : function
        Power spectrum used to generate the time series. It is expected to be
        a function mapping the input frequencies to the periodogram.
    npoints : int
        Number of points in the output time series.
    spacing : `~astropy.units.Quantity`
        Sample spacing, inverse of the sampling rate. The units are inherited by the resulting time axis.
    nchunks : int, optional
        Factor by which to multiply the length of the time series to avoid red noise leakage. Default is 10.
    random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional
        Defines random number generator initialisation.
        Passed to `~gammapy.utils.random.get_random_state`. Default is "random-seed".
    power_spectrum_params : dict, optional
        Dictionary of parameters to be provided to the power spectrum function.
    mean : float, `~astropy.units.Quantity`, optional
        Desired mean of the final series. Default is 0.
    std : float, `~astropy.units.Quantity`, optional
        Desired standard deviation of the final series. Default is 1.
    poisson : bool, optional
        Whether to apply poissonian noise to the final time series. Default is False.


    Returns
    -------
    time_series : `~numpy.ndarray`
        Simulated time series.
    time_axis : `~astropy.units.Quantity`
        Time axis of the series in the same units as 'spacing'. It will be defined with length 'npoints', from 0 to
        'npoints'*'spacing'.

    Examples
    --------
    To pass the function to be used in the simlation one can use either the 'lambda' keyword or an extended definition.
    Parameters of the function can be passed using the 'power_spectrum_params' keyword.
    For example, these are three ways to pass a power law (red noise) with index 2:

    >>> from gammapy.stats import TimmerKonig_lightcurve_simulator
    >>> import astropy.units as u
    >>> def powerlaw(x):
    ...     return x**(-2)
    >>> def powerlaw_with_parameters(x, i):
    ...     return x**(-i)
    >>> ts, ta = TimmerKonig_lightcurve_simulator(lambda x: x**(-2), 20, 1*u.h)
    >>> ts2, ta2 = TimmerKonig_lightcurve_simulator(powerlaw, 20, 1*u.h)
    >>> ts3, ta3 = TimmerKonig_lightcurve_simulator(powerlaw_with_parameters,
    ...                                            20, 1*u.h, power_spectrum_params={"i":2})

    References
    ----------
    .. [Timmer1995] "On generating power law noise", J. Timmer and M, Konig, section 3
       https://ui.adsabs.harvard.edu/abs/1995A%26A...300..707T/abstract
    """
    if not callable(power_spectrum):
        raise ValueError(
            "The power spectrum has to be provided as a callable function."
        )

    if not isinstance(npoints * nchunks, int):
        raise TypeError("npoints and nchunks must be integers")

    if poisson:
        if isinstance(mean, u.Quantity):
            wmean = mean.value * spacing.value
        else:
            wmean = mean * spacing.value
        if wmean < 1.0:
            raise Warning(
                "Poisson noise was requested but the target mean is too low - resulting counts will likely be 0."
            )

    random_state = get_random_state(random_state)

    npoints_ext = npoints * nchunks

    frequencies = np.fft.fftfreq(npoints_ext, spacing.value)

    # To obtain real data only the positive or negative part of the frequency is necessary.
    real_frequencies = np.sort(np.abs(frequencies[frequencies < 0]))

    if power_spectrum_params:
        periodogram = power_spectrum(real_frequencies, **power_spectrum_params)
    else:
        periodogram = power_spectrum(real_frequencies)

    real_part = random_state.normal(0, 1, len(periodogram) - 1)
    imaginary_part = random_state.normal(0, 1, len(periodogram) - 1)

    # Nyquist frequency component handling
    if npoints_ext % 2 == 0:
        idx0 = -2
        random_factor = random_state.normal(0, 1)
    else:
        idx0 = -1
        random_factor = random_state.normal(0, 1) + 1j * random_state.normal(0, 1)

    fourier_coeffs = np.concatenate(
        [
            np.sqrt(0.5 * periodogram[:-1]) * (real_part + 1j * imaginary_part),
            np.sqrt(0.5 * periodogram[-1:]) * random_factor,
        ]
    )
    fourier_coeffs = np.concatenate(
        [fourier_coeffs, np.conjugate(fourier_coeffs[idx0::-1])]
    )

    fourier_coeffs = np.insert(fourier_coeffs, 0, 0)
    time_series = np.fft.ifft(fourier_coeffs).real

    ndiv = npoints_ext // (2 * nchunks)
    setstart = npoints_ext // 2 - ndiv
    setend = npoints_ext // 2 + ndiv
    if npoints % 2 != 0:
        setend = setend + 1
    time_series = time_series[setstart:setend]

    time_series = (time_series - time_series.mean()) / time_series.std()
    time_series = time_series * std + mean

    if poisson:
        time_series = (
            random_state.poisson(
                np.where(time_series >= 0, time_series, 0) * spacing.value
            )
            / spacing.value
        )

    time_axis = np.linspace(0, npoints * spacing.value, npoints) * spacing.unit

    return time_series, time_axis