Repository: gammapy/gammapy
Pull request #5269 (2024-05-17 09:22:30)

Add leakage protection factor in the Timmer Konig algorithm
This PR adds a small modification to the TimmerKonig algorithm. Now there is the possibility of using an additional (optional) input which is a multiplicative factor of the time series length `npoints` to avoid spectral leakage. The final time series of length `npoints` is extracted randomly at the end from the longer time series of length `npoints * leakage_protection`.
A specific test was added.



Patch:
diff --git a/gammapy/stats/variability.py b/gammapy/stats/variability.py
index 88f02c380..c3377f699 100644
--- a/gammapy/stats/variability.py
+++ b/gammapy/stats/variability.py
@@ -289,6 +289,7 @@ def TimmerKonig_lightcurve_simulator(
     power_spectrum,
     npoints,
     spacing,
+    leakage_protection=1.0,
     random_state="random-seed",
     power_spectrum_params=None,
 ):
@@ -303,6 +304,8 @@ def TimmerKonig_lightcurve_simulator(
         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.
+    leakage_protection : float, optional
+        Factor by which to multiply the length of the time series to avoid red noise leakage. Default is 1.0.
     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 "random-seed".
@@ -346,7 +349,16 @@ def TimmerKonig_lightcurve_simulator(
 
     random_state = get_random_state(random_state)
 
-    frequencies = np.fft.fftfreq(npoints, spacing.value)
+    if leakage_protection < 1.0:
+        raise ValueError("The leakage protection factor must be at least 1.0.")
+
+    npoints_extended = int(npoints * leakage_protection)
+    if leakage_protection > 1.0 and npoints_extended == npoints:
+        raise Warning(
+            "The extended time series is the same length as the final desired time series."
+            "To have an effective protection from noise leakage and avoid aliasing, increase the leakage_protection option."
+        )
+    frequencies = np.fft.fftfreq(npoints_extended, 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]))
@@ -380,6 +392,10 @@ def TimmerKonig_lightcurve_simulator(
     fourier_coeffs = np.insert(fourier_coeffs, 0, 0)
     time_series = np.fft.ifft(fourier_coeffs).real
 
+    if npoints_extended > npoints:
+        extract = random_state.randint(1, npoints_extended - npoints)
+        time_series = np.take(time_series, range(extract, extract + npoints))
+
     tmax = time_series.max()
     if tmax < np.abs(time_series.min()):
         time_series = -time_series


Test patch:
diff --git a/gammapy/stats/tests/test_variability.py b/gammapy/stats/tests/test_variability.py
index 41912a4fe..bee57c1f5 100644
--- a/gammapy/stats/tests/test_variability.py
+++ b/gammapy/stats/tests/test_variability.py
@@ -164,11 +164,28 @@ def test_tk():
         temp, 15, 1 * u.h, power_spectrum_params=params
     )
 
+    time_series4, time_axis4 = TimmerKonig_lightcurve_simulator(
+        lambda x: x ** (-3), 20, 1 * u.s, leakage_protection=100.0
+    )
+
     assert len(time_series) == 20
     assert isinstance(time_axis, u.Quantity)
     assert time_axis.unit == u.s
     assert len(time_series2) == 21
     assert len(time_series3) == 15
+    assert len(time_series4) == 20
+
+
+def test_tk_badleakage():
+    with pytest.raises(ValueError):
+        TimmerKonig_lightcurve_simulator(
+            lambda x: x ** (-3), 20, 1 * u.s, leakage_protection=0.1
+        )
+
+    with pytest.raises(Warning):
+        TimmerKonig_lightcurve_simulator(
+            lambda x: x ** (-3), 20, 1 * u.s, leakage_protection=1.001
+        )
 
 
 def test_structure_function():