Repository: gammapy/gammapy
Pull request #5423 (2024-09-03 13:35:21)

Function for the discrete cross correlation function
This PR introduces a function in `variability.py` to compute the discrete correlation function as proposed [by Edelson et al](https://ui.adsabs.harvard.edu/abs/1988ApJ...333..646E/abstract). I added the usual wrapper function to use directly `FluxPoints` objects.



Patch:
diff --git a/gammapy/estimators/utils.py b/gammapy/estimators/utils.py
index da76c2d40..61b476160 100644
--- a/gammapy/estimators/utils.py
+++ b/gammapy/estimators/utils.py
@@ -15,7 +15,12 @@ from gammapy.modeling.models import (
     PowerLawSpectralModel,
     SkyModel,
 )
-from gammapy.stats import compute_flux_doubling, compute_fpp, compute_fvar
+from gammapy.stats import (
+    compute_flux_doubling,
+    compute_fpp,
+    compute_fvar,
+    discrete_correlation,
+)
 from gammapy.stats.utils import ts_to_sigma
 from .map.core import FluxMaps
 
@@ -29,6 +34,7 @@ __all__ = [
     "compute_lightcurve_fvar",
     "compute_lightcurve_fpp",
     "compute_lightcurve_doublingtime",
+    "compute_lightcurve_discrete_correlation",
 ]
 
 
@@ -536,6 +542,88 @@ def compute_lightcurve_doublingtime(lightcurve, flux_quantity="flux"):
     return table
 
 
+def compute_lightcurve_discrete_correlation(
+    lightcurve1, lightcurve2=None, flux_quantity="flux", tau=None
+):
+    r"""Compute the discrete correlation function for two lightcurves, or the discrete autocorrelation if only one lightcurve is provided.
+    NaN values will be ignored in the computation in order to account for possible gaps in the data.
+
+    Internally calls the `~gammapy.stats.discrete_correlation` function
+
+    Parameters
+    ----------
+    lightcurve1 : `~gammapy.estimators.FluxPoints`
+        The first lightcurve object.
+    lightcurve2 : `~gammapy.estimators.FluxPoints`, optional
+        The second lightcurve object. If not provided, the autocorrelation for the first lightcurve will be computed.
+        Default is None.
+    flux_quantity : str
+        Flux quantity to use for calculation. Should be 'dnde', 'flux', 'e2dnde' or 'eflux'.
+        The choice does not affect the computation. Default is 'flux'.
+    tau : `~astropy.units.Quantity`, optional
+        Size of the bins to compute the discrete correlation.
+        If None, the bin size will be double the bins of the first lightcurve. Default is None.
+
+
+    Returns
+    -------
+    discrete_correlation_dict : dict
+        Dictionary containing:
+            "bins" : the array of discrete time bins
+            "discrete_correlation" : discrete correlation function values
+            "discrete_correlation_err" : associated error
+
+    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 = getattr(lightcurve1, flux_quantity)
+    flux_err1 = getattr(lightcurve1, flux_quantity + "_err")
+    coords1 = lightcurve1.geom.axes["time"].center
+    axis = flux1.geom.axes.index_data("time")
+
+    if tau is None:
+        tau = (coords1[-1] - coords1[0]) / (0.5 * len(coords1))
+
+    if lightcurve2:
+        flux2 = getattr(lightcurve2, flux_quantity)
+        flux_err2 = getattr(lightcurve2, flux_quantity + "_err")
+        coords2 = lightcurve2.geom.axes["time"].center
+        bins, dcf, dcf_err = discrete_correlation(
+            flux1.data,
+            flux_err1.data,
+            flux2.data,
+            flux_err2.data,
+            coords1,
+            coords2,
+            tau,
+            axis,
+        )
+
+    else:
+        bins, dcf, dcf_err = discrete_correlation(
+            flux1.data,
+            flux_err1.data,
+            flux1.data,
+            flux_err1.data,
+            coords1,
+            coords1,
+            tau,
+            axis,
+        )
+
+    discrete_correlation_dict = {
+        "bins": bins,
+        "discrete_correlation": dcf,
+        "discrete_correlation_err": dcf_err,
+    }
+
+    return discrete_correlation_dict
+
+
 def get_edges_fixed_bins(fluxpoint, group_size, axis_name="energy"):
     """Rebin the flux point to combine value adjacent bins.
 
@@ -817,7 +905,6 @@ def combine_flux_maps(
         reference_model = maps[0].reference_model
 
     if method == "gaussian_errors":
-
         means = [map_.dnde.copy() for map_ in maps]
         sigmas = [map_.dnde_err.copy() for map_ in maps]
         # compensate for the ts deviation from gaussian approximation expectation in each map
@@ -833,7 +920,6 @@ def combine_flux_maps(
         sigma = sigmas[0]
 
         for k in range(1, len(means)):
-
             mean_k = means[k].quantity.to_value(mean.unit)
             sigma_k = sigmas[k].quantity.to_value(sigma.unit)
 
diff --git a/gammapy/stats/__init__.py b/gammapy/stats/__init__.py
index 9689ddab8..b5216e560 100644
--- a/gammapy/stats/__init__.py
+++ b/gammapy/stats/__init__.py
@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Statistics."""
+
 from .counts_statistic import CashCountsStatistic, WStatCountsStatistic
 from .fit_statistics import cash, cstat, get_wstat_gof_terms, get_wstat_mu_bkg, wstat
 from .fit_statistics_cython import (
@@ -13,6 +14,7 @@ from .variability import (
     compute_flux_doubling,
     compute_fpp,
     compute_fvar,
+    discrete_correlation,
     structure_function,
 )
 
@@ -32,5 +34,6 @@ __all__ = [
     "compute_flux_doubling",
     "compute_chisq",
     "structure_function",
+    "discrete_correlation",
     "TimmerKonig_lightcurve_simulator",
 ]
diff --git a/gammapy/stats/variability.py b/gammapy/stats/variability.py
index 31b56aad5..0f6bf121e 100644
--- a/gammapy/stats/variability.py
+++ b/gammapy/stats/variability.py
@@ -11,6 +11,7 @@ __all__ = [
     "compute_flux_doubling",
     "structure_function",
     "TimmerKonig_lightcurve_simulator",
+    "discrete_correlation",
 ]
 
 
@@ -286,6 +287,93 @@ def structure_function(flux, flux_err, time, tdelta_precision=5):
     return sf, distances
 
 
+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)
+
+
 def TimmerKonig_lightcurve_simulator(
     power_spectrum,
     npoints,


Test patch:
diff --git a/gammapy/estimators/tests/test_utils.py b/gammapy/estimators/tests/test_utils.py
index cbd38d51c..f3513005b 100644
--- a/gammapy/estimators/tests/test_utils.py
+++ b/gammapy/estimators/tests/test_utils.py
@@ -12,6 +12,7 @@ from gammapy.estimators.utils import (
     compute_lightcurve_doublingtime,
     compute_lightcurve_fpp,
     compute_lightcurve_fvar,
+    compute_lightcurve_discrete_correlation,
     find_peaks,
     find_peaks_in_flux_map,
     get_rebinned_axis,
@@ -21,6 +22,64 @@ from gammapy.maps import Map, MapAxis
 from gammapy.utils.testing import assert_time_allclose, requires_data
 
 
+@pytest.fixture()
+def lc():
+    meta = dict(TIMESYS="utc", SED_TYPE="flux")
+
+    table = Table(
+        meta=meta,
+        data=[
+            Column(Time(["2010-01-01", "2010-01-03"]).mjd, "time_min"),
+            Column(Time(["2010-01-03", "2010-01-10"]).mjd, "time_max"),
+            Column([[1.0, 2.0], [1.0, 2.0]], "e_min", unit="TeV"),
+            Column([[2.0, 5.0], [2.0, 5.0]], "e_max", unit="TeV"),
+            Column([[1e-11, 4e-12], [3e-11, 7e-12]], "flux", unit="cm-2 s-1"),
+            Column(
+                [[0.1e-11, 0.4e-12], [0.3e-11, 0.7e-12]], "flux_err", unit="cm-2 s-1"
+            ),
+            Column([[np.nan, np.nan], [3.6e-11, 1e-11]], "flux_ul", unit="cm-2 s-1"),
+            Column([[False, False], [True, True]], "is_ul"),
+            Column([[True, True], [True, True]], "success"),
+        ],
+    )
+
+    return FluxPoints.from_table(table=table, format="lightcurve")
+
+
+@pytest.fixture(scope="session")
+def lc2():
+    meta = dict(TIMESYS="utc", SED_TYPE="flux")
+
+    table = Table(
+        meta=meta,
+        data=[
+            Column(Time(["2010-01-01", "2010-01-03", "2010-01-05"]).mjd, "time_min"),
+            Column(Time(["2010-01-03", "2010-01-05", "2010-01-07"]).mjd, "time_max"),
+            Column([[1.0, 2.0], [1.0, 2.0], [1.0, 2.0]], "e_min", unit="GeV"),
+            Column([[2.0, 5.0], [2.0, 5.0], [2.0, 5.0]], "e_max", unit="GeV"),
+            Column(
+                [[1.51e-7, 3.4e-8], [3.1e-7, 6.7e-8], [3.1e-7, 7.5e-8]],
+                "flux",
+                unit="m-2 s-1",
+            ),
+            Column(
+                [[0.1e-7, 0.4e-8], [0.3e-7, 0.7e-8], [0.31e-7, 0.72e-8]],
+                "flux_err",
+                unit="m-2 s-1",
+            ),
+            Column(
+                [[np.nan, np.nan], [3.6e-7, 1e-7], [3.6e-7, 1e-7]],
+                "flux_ul",
+                unit="m-2 s-1",
+            ),
+            Column([[False, False], [True, True], [True, True]], "is_ul"),
+            Column([[True, True], [True, True], [True, True]], "success"),
+        ],
+    )
+
+    return FluxPoints.from_table(table=table, format="lightcurve")
+
+
 class TestFindPeaks:
     def test_simple(self):
         """Test a simple example"""
@@ -138,33 +197,8 @@ def test_resample_energy_edges(spectrum_dataset):
     assert_allclose(np.squeeze(grouped.background)[-1], 200)
 
 
-def lc():
-    meta = dict(TIMESYS="utc", SED_TYPE="flux")
-
-    table = Table(
-        meta=meta,
-        data=[
-            Column(Time(["2010-01-01", "2010-01-03"]).mjd, "time_min"),
-            Column(Time(["2010-01-03", "2010-01-10"]).mjd, "time_max"),
-            Column([[1.0, 2.0], [1.0, 2.0]], "e_min", unit="TeV"),
-            Column([[2.0, 5.0], [2.0, 5.0]], "e_max", unit="TeV"),
-            Column([[1e-11, 4e-12], [3e-11, 7e-12]], "flux", unit="cm-2 s-1"),
-            Column(
-                [[0.1e-11, 0.4e-12], [0.3e-11, 0.7e-12]], "flux_err", unit="cm-2 s-1"
-            ),
-            Column([[np.nan, np.nan], [3.6e-11, 1e-11]], "flux_ul", unit="cm-2 s-1"),
-            Column([[False, False], [True, True]], "is_ul"),
-            Column([[True, True], [True, True]], "success"),
-        ],
-    )
-
-    return FluxPoints.from_table(table=table, format="lightcurve")
-
-
-def test_compute_lightcurve_fvar():
-    lightcurve = lc()
-
-    fvar = compute_lightcurve_fvar(lightcurve)
+def test_compute_lightcurve_fvar(lc):
+    fvar = compute_lightcurve_fvar(lc)
     ffvar = fvar["fvar"].quantity
     ffvar_err = fvar["fvar_err"].quantity
 
@@ -172,10 +206,8 @@ def test_compute_lightcurve_fvar():
     assert_allclose(ffvar_err, [[[0.0795621]], [[0.074706]]])
 
 
-def test_compute_lightcurve_fpp():
-    lightcurve = lc()
-
-    fpp = compute_lightcurve_fpp(lightcurve)
+def test_compute_lightcurve_fpp(lc):
+    fpp = compute_lightcurve_fpp(lc)
     ffpp = fpp["fpp"].quantity
     ffpp_err = fpp["fpp_err"].quantity
 
@@ -183,10 +215,8 @@ def test_compute_lightcurve_fpp():
     assert_allclose(ffpp_err, [[[0.07930673]], [[0.07397653]]])
 
 
-def test_compute_lightcurve_doublingtime():
-    lightcurve = lc()
-
-    dtime = compute_lightcurve_doublingtime(lightcurve)
+def test_compute_lightcurve_doublingtime(lc):
+    dtime = compute_lightcurve_doublingtime(lc)
     ddtime = dtime["doublingtime"].quantity
     ddtime_err = dtime["doubling_err"].quantity
     dcoord = dtime["doubling_coord"]
@@ -199,6 +229,48 @@ def test_compute_lightcurve_doublingtime():
     )
 
 
+def test_compute_dcf(lc, lc2):
+    dict = compute_lightcurve_discrete_correlation(lc, lc2, tau=3 * u.d)
+
+    assert_allclose(dict["bins"], [-388800.0, -129600.0, 129600.0, 388800.0] * u.s)
+    assert_allclose(
+        dict["discrete_correlation"],
+        [
+            [-0.760599, -1.052783],
+            [-0.760599, -0.537134],
+            [1.014132, 1.059945],
+            [-1.521198, -1.589918],
+        ],
+        rtol=1e-6,
+    )
+    assert_allclose(
+        dict["discrete_correlation_err"],
+        [[np.nan, np.nan], [np.nan, np.nan], [0.310513, 0.372241], [np.nan, np.nan]],
+        rtol=1e-6,
+    )
+
+    dict2 = compute_lightcurve_discrete_correlation(lc2, tau=3 * u.d)
+    assert_allclose(dict2["bins"], [-388800.0, -129600.0, 129600.0, 388800.0] * u.s)
+    assert_allclose(
+        dict2["discrete_correlation"],
+        [
+            [-1.11074, -1.448801],
+            [-0.277685, -0.124862],
+            [0.55537, 0.629465],
+            [-1.11074, -1.448801],
+        ],
+        rtol=1e-5,
+    )
+    assert_allclose(
+        dict2["discrete_correlation_err"],
+        [[np.nan, np.nan], [1.178118, 0.868782], [0.589059, 0.53472], [np.nan, np.nan]],
+        rtol=1e-6,
+    )
+
+    dict3 = compute_lightcurve_discrete_correlation(lc2)
+    assert_allclose(dict3["bins"], [-345600.0, -115200.0, 115200.0, 345600.0] * u.s)
+
+
 @requires_data()
 def test_get_rebinned_axis():
     lc_1d = FluxPoints.read(
diff --git a/gammapy/stats/tests/test_variability.py b/gammapy/stats/tests/test_variability.py
index b9ac407c4..cd404dbdf 100644
--- a/gammapy/stats/tests/test_variability.py
+++ b/gammapy/stats/tests/test_variability.py
@@ -13,6 +13,7 @@ from gammapy.stats.variability import (
     compute_fpp,
     compute_fvar,
     structure_function,
+    discrete_correlation,
 )
 from gammapy.utils.testing import assert_quantity_allclose
 
@@ -36,6 +37,7 @@ def lc_table():
     return table
 
 
+@pytest.fixture(scope="session")
 def lc():
     meta = dict(TIMESYS="utc", SED_TYPE="flux")
 
@@ -257,3 +259,96 @@ def test_structure_function():
         distances,
         [3600.0, 7000.0, 7200.0, 10800.0, 14200.0, 14400.0, 17800.0, 21400.0] * u.s,
     )
+
+
+def test_dcf():
+    flux = np.array(
+        [
+            [1e-11, 4e-12],
+            [3e-11, 2.5e-12],
+            [1e-11, 1e-12],
+            [0.8e-11, 0.8e-12],
+            [1e-11, 1e-12],
+        ]
+    )
+    flux2 = np.array(
+        [
+            [1e-11, 3.5e-12],
+            [4e-11, 1.5e-12],
+            [1.7e-11, 2e-12],
+            [0.4e-11, 1.8e-12],
+            [1e-11, 1.5e-12],
+            [1.1e-11, 0.7e-12],
+        ]
+    )
+
+    flux_err = np.array(
+        [
+            [0.1e-11, 0.4e-12],
+            [0.3e-11, 0.2e-12],
+            [0.1e-11, 0.1e-12],
+            [0.08e-11, 0.08e-12],
+            [0.1e-11, 0.1e-12],
+        ]
+    )
+    flux_err2 = np.array(
+        [
+            [0.2e-11, 0.24e-12],
+            [0.4e-11, 0.21e-12],
+            [0.07e-11, 0.14e-12],
+            [0.12e-11, 0.18e-12],
+            [0.1e-11, 0.13e-12],
+            [0.4e-11, 0.2e-12],
+        ]
+    )
+
+    time = (
+        np.array(
+            [6.31157019e08, 6.31160619e08, 6.31164219e08, 6.31171419e08, 6.31178419e08]
+        )
+        * u.s
+    )
+    time2 = (
+        np.array(
+            [
+                6.31156919e08,
+                6.31160519e08,
+                6.31163819e08,
+                6.31171419e08,
+                6.31178219e08,
+                6.31178819e08,
+            ]
+        )
+        * u.s
+    )
+
+    bins, dcf, dcf_err = discrete_correlation(
+        flux, flux_err, flux2, flux_err2, time, time2, tau=1.5e4 * u.s
+    )
+
+    assert_allclose(
+        dcf,
+        [
+            [-0.332361, -1.009638],
+            [-0.065639, 0.313054],
+            [0.215329, 0.133132],
+            [-0.373906, -0.56788],
+        ],
+        rtol=1e-5,
+    )
+
+    assert_allclose(
+        dcf_err,
+        [
+            [0.34881, 0.553113],
+            [0.236548, 0.174538],
+            [0.401313, 0.361757],
+            [0.820525, 1.204655],
+        ],
+        rtol=1e-5,
+    )
+
+    assert_allclose(
+        bins,
+        u.Quantity([-22500.0, -7500.0, 7500.0, 22500.0], unit="s"),
+    )