From d5a50d3c8a95163379ce68ab010fae90d77741e4 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Sat, 28 Feb 2026 23:55:50 +0100
Subject: [PATCH 1/9] Make coherence work again with Numba

---
 stingray/fourier.py            | 272 +++++++++++++++++++++++++++------
 stingray/tests/test_fourier.py |   2 +-
 stingray/utils.py              |   5 +-
 3 files changed, 233 insertions(+), 46 deletions(-)

diff --git a/stingray/fourier.py b/stingray/fourier.py
index 723408de5..3542e90be 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -25,7 +25,6 @@
     vectorize,
     float64,
     int64,
-    UniTuple,
 )
 
 
@@ -808,9 +807,7 @@ def _apply_low_lim_to_coherence_uncertainty(coherence, uncertainty, min_uncertai
 
 
 @vectorize(
-    [
-        "UniTuple(float64[:, :], 2)(float64[:, :], float64[:, :], float64[:, :], float64, float64, int64, float64)"
-    ],
+    ["float64(complex128, float64, float64, float64, float64, int64, float64)"],
     nopython=True,
 )
 def _raw_coherence(
@@ -820,7 +817,7 @@ def _raw_coherence(
     power1_noise,
     power2_noise,
     n_ave,
-    intrinsic_coherence=1,
+    intrinsic_coherence,
 ):
     """
     Raw coherence estimations from cross and power spectra.
@@ -854,20 +851,16 @@ def _raw_coherence(
     coherence : float `np.array`
         The raw coherence values at all frequencies.
     """
-    bsq = bias_term(
-        power1, power2, power1_noise, power2_noise, n_ave, intrinsic_coherence=intrinsic_coherence
-    )
+    bsq = _bias_term(power1, power2, power1_noise, power2_noise, n_ave, intrinsic_coherence)
     num = (cross_power * np.conj(cross_power)).real - bsq
     if num < 0:
         num = 0
+
     den = power1 * power2
 
     coherence = num / den
-    min_uncertainty = 1 / n_ave
-    uncertainty = (2**0.5 * coherence * (1 - coherence)) / (np.sqrt(coherence) * n_ave**0.5)
-    uncertainty = _apply_low_lim_to_coherence_uncertainty(coherence, uncertainty, min_uncertainty)
 
-    return coherence, uncertainty
+    return coherence
 
 
 def raw_coherence(
@@ -877,7 +870,7 @@ def raw_coherence(
     power1_noise,
     power2_noise,
     n_ave,
-    intrinsic_coherence=1,
+    intrinsic_coherence=1.0,
     return_uncertainty=False,
 ):
     r"""
@@ -925,16 +918,21 @@ def raw_coherence(
     ----------
     .. [#] https://doi.org/10.1093/mnras/stz2409
     """
-    coherence, uncertainty = _raw_coherence(
+    coherence = _raw_coherence(
         cross_power,
         power1,
         power2,
         power1_noise,
         power2_noise,
         n_ave,
-        intrinsic_coherence=intrinsic_coherence,
+        intrinsic_coherence,
     )
     if return_uncertainty:
+        min_uncertainty = 1 / n_ave
+        uncertainty = (2**0.5 * coherence * (1 - coherence)) / (np.sqrt(coherence) * n_ave**0.5)
+        uncertainty = _apply_low_lim_to_coherence_uncertainty(
+            coherence, uncertainty, min_uncertainty
+        )
         return coherence, uncertainty
     return coherence
 
@@ -973,11 +971,11 @@ def _intrinsic_coherence_uncertainties(
 
 
 @vectorize(
-    ["UniTuple(float64[:, :], 2)(float64[:, :], float64[:, :], float64, float64, int64)"],
+    ["bool(float64, float64, float64, float64, int64, float64)"],
     nopython=True,
 )
 def check_powers_for_intrinsic_coherence(
-    power1, power2, power1_noise, power2_noise, n_ave, threshold=3
+    power1, power2, power1_noise, power2_noise, n_ave, threshold
 ):
     """Check if the powers are above the threshold for the intrinsic coherence to be well defined.
 
@@ -1009,15 +1007,12 @@ def check_powers_for_intrinsic_coherence(
     """
     # Consider low power when the power is less than 3 sigma above the noise level.
     # This is somewhat arbitrary, better criteria suggestions are welcome.
-    low_power1 = (power1 - power1_noise) <= 3 * power1_noise / np.sqrt(n_ave)
-    low_power2 = (power2 - power2_noise) <= 3 * power2_noise / np.sqrt(n_ave)
+    low_power1 = (power1 - power1_noise) <= threshold * power1_noise / np.sqrt(n_ave)
+    low_power2 = (power2 - power2_noise) <= threshold * power2_noise / np.sqrt(n_ave)
     return low_power1 or low_power2
 
 
-@vectorize(
-    ["UniTuple(float64[:, :], 2)(float64, float64, float64, float64, float64, int64, float64)"],
-    nopython=True,
-)
+@njit
 def _intrinsic_coherence(
     cross_power,
     power1,
@@ -1069,7 +1064,7 @@ def _intrinsic_coherence(
         1997, ApJ 474, L43, eq. 8.
     """
 
-    if check_powers_for_intrinsic_coherence(power1, power2, power1_noise, power2_noise, n_ave):
+    if check_powers_for_intrinsic_coherence(power1, power2, power1_noise, power2_noise, n_ave, 3):
         # If the powers are too close to the noise level, the coherence is not well defined.
         return np.nan, np.nan
 
@@ -1092,10 +1087,7 @@ def _intrinsic_coherence(
     return coherence, uncertainty
 
 
-@vectorize(
-    ["UniTuple(float64[:, :], 2)(float64, float64, float64, float64, float64, int64)"],
-    nopython=True,
-)
+@njit
 def _intrinsic_coherence_with_adjusted_bias(
     cross_power,
     power1,
@@ -1145,12 +1137,14 @@ def _intrinsic_coherence_with_adjusted_bias(
     uncertainty : float `np.array`, optional
         The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak
         1997, ApJ 474, L43, eq. 8.
+    non_converged_flag : bool `np.array`
+        Whether the bias term adjustment procedure converged (True) or not (False).
 
     """
 
-    if check_powers_for_intrinsic_coherence(power1, power2, power1_noise, power2_noise, n_ave):
+    if check_powers_for_intrinsic_coherence(power1, power2, power1_noise, power2_noise, n_ave, 3):
         # If the powers are too close to the noise level, the coherence is not well defined.
-        return np.nan, np.nan
+        return np.nan, np.nan, False
 
     # Consider low power when the power is less than 3 sigma above the noise level.
     # This is somewhat arbitrary, better criteria suggestions are welcome.
@@ -1159,35 +1153,125 @@ def _intrinsic_coherence_with_adjusted_bias(
 
     current_coherence = 1.0
 
+    converged = False
     for _ in range(40):
         bsq = _bias_term(power1, power2, power1_noise, power2_noise, n_ave, current_coherence)
         num = (cross_power * np.conj(cross_power)).real - bsq
         if num <= 0:
-            coherence = 0
             # if the current coherence still drops, the bias term increases,
             # so at this point it is useless to keep iterating.
-            return 0, np.nan
+            return 0.0, np.nan, False
 
         den = power1_sub * power2_sub
 
         coherence = num / den
         if np.abs(current_coherence - coherence) < atol:
+            converged = True
             break
         current_coherence = coherence
-    else:
-        warnings.warn(
-            "The iterative procedure to adjust the bias term did not converge after 40 iterations. "
-            "Consider rebinning the spectra to increase the signal-to-noise ratio, or to use a "
-            "more permissive tolerance (``atol`` parameter)."
-        )
 
     uncertainty = _intrinsic_coherence_uncertainties(
         bsq, num, power1_sub, power2_sub, power1_noise, power2_noise, coherence, n_ave
     )
-    return coherence, uncertainty
+    return coherence, uncertainty, not converged
 
 
-def intrinsic_coherence(
+@njit
+def intrinsic_coherence_array(
+    cross_power,
+    power1,
+    power2,
+    power1_noise,
+    power2_noise,
+    n_ave,
+    return_uncertainty=False,
+    adjust_bias=False,
+    atol=0.01,
+):
+    r"""
+    Intrinsic coherence estimations from cross and power spectra.
+
+    Vaughan & Nowak 1997, ApJ 474, L43
+
+    .. math::
+
+        \tilde{\gamma}^2(f) = \frac{|\langle \tilde{C}(f) \rangle|^2 - \tilde{b}^2}
+        {(\langle \tilde{P}_1(f) \rangle - \tilde{P}_{1, \rm noise})
+        (\langle \tilde{P}_2(f) \rangle - \tilde{P}_{2, \rm noise})}
+
+    where :math:`\tilde{b}^2` is the bias term that accounts for the contribution of Poisson
+    noise to the cross spectrum (see :func:`stingray.fourier.bias_term`), and tilde generally
+    indicates noisy measurements of the cross spectrum :math:`C` and the power spectra :math:`P_n`.
+    The terms :math:`\tilde{P}_{\rm n, \rm noise}` are the estimates of the contribution of
+    Poisson noise to the power spectra.
+
+    The bias term depends on the intrinsic coherence itself, so it can be calculated iteratively
+    following the procedure from sec. 5 of Ingram 2019, MNRAS 489, 3927, which adjusts the bias
+    term until convergence is reached. This is done if ``adjust_bias`` is set to True. It is
+    typically a very small correction.
+
+    For errors, assumes **high powers, high coherence**. See eq. 8 of the paper.
+    Powers below 3 sigma above the noise level (P_noise / sqrt(MW)) are considered too low,
+    and the coherence will be set to NaN.
+
+    TODO: implement a more general treatment of the uncertainty.
+
+    Parameters
+    ----------
+    cross_power : complex `np.array`
+        cross spectrum
+    power1 : float `np.array`
+        sub-band periodogram
+    power2 : float `np.array`
+        reference-band periodogram
+    power1_noise : float
+        Poisson noise level of the sub-band periodogram
+    power2_noise : float
+        Poisson noise level of the reference-band periodogram
+    n_ave : int
+        number of intervals that have been averaged to obtain the input spectra
+
+    Other Parameters
+    ----------------
+    return_uncertainty : bool, default False
+        Whether to return the uncertainty on the coherence, calculated according to
+        Vaughan & Nowak 1997, ApJ 474, L43, eq. 8.
+    atol: float, default 0.01
+        The absolute tolerance for the convergence of the iterative procedure to adjust
+        the bias term.Only relevant if ``adjust_bias`` is True.
+
+    Returns
+    -------
+    coherence : float `np.array` or tuple of two float `np.array`
+        The intrinsic coherence values at all frequencies. It is 0 if the numerator
+        of the coherence calculation is negative, and NaN if the powers are too low compared
+        to the noise level.
+    uncertainty : float `np.array`, optional
+        The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak
+        1997, ApJ 474, L43, eq. 8.
+
+    """
+
+    coherence = np.empty(cross_power.shape, dtype=np.float64)
+    uncertainty = np.empty(cross_power.shape, dtype=np.float64)
+    flags = np.zeros(cross_power.shape, dtype=np.bool_)
+
+    for i in range(cross_power.size):
+        if adjust_bias:
+            coherence[i], uncertainty[i], flags[i] = _intrinsic_coherence_with_adjusted_bias(
+                cross_power[i], power1[i], power2[i], power1_noise, power2_noise, n_ave, atol=atol
+            )
+
+        else:
+            coherence[i], uncertainty[i] = _intrinsic_coherence(
+                cross_power[i], power1[i], power2[i], power1_noise, power2_noise, n_ave
+            )
+
+    return coherence, uncertainty, flags
+
+
+@njit
+def intrinsic_coherence_single(
     cross_power,
     power1,
     power2,
@@ -1261,16 +1345,118 @@ def intrinsic_coherence(
         1997, ApJ 474, L43, eq. 8.
 
     """
+
     if adjust_bias:
-        result = _intrinsic_coherence_with_adjusted_bias(
+        coherence, uncertainty, flags = _intrinsic_coherence_with_adjusted_bias(
             cross_power, power1, power2, power1_noise, power2_noise, n_ave, atol=atol
         )
     else:
-        result = _intrinsic_coherence(
+        coherence, uncertainty = _intrinsic_coherence(
             cross_power, power1, power2, power1_noise, power2_noise, n_ave
         )
-    coherence, uncertainty = result
+        flags = False
+
+    return coherence, uncertainty, flags
+
+
+def intrinsic_coherence(
+    cross_power,
+    power1,
+    power2,
+    power1_noise,
+    power2_noise,
+    n_ave,
+    return_uncertainty=False,
+    adjust_bias=False,
+    atol=0.01,
+):
+    r"""
+    Intrinsic coherence estimations from cross and power spectra.
+
+    Vaughan & Nowak 1997, ApJ 474, L43
+
+    .. math::
+
+        \tilde{\gamma}^2(f) = \frac{|\langle \tilde{C}(f) \rangle|^2 - \tilde{b}^2}
+        {(\langle \tilde{P}_1(f) \rangle - \tilde{P}_{1, \rm noise})
+        (\langle \tilde{P}_2(f) \rangle - \tilde{P}_{2, \rm noise})}
+
+    where :math:`\tilde{b}^2` is the bias term that accounts for the contribution of Poisson
+    noise to the cross spectrum (see :func:`stingray.fourier.bias_term`), and tilde generally
+    indicates noisy measurements of the cross spectrum :math:`C` and the power spectra :math:`P_n`.
+    The terms :math:`\tilde{P}_{\rm n, \rm noise}` are the estimates of the contribution of
+    Poisson noise to the power spectra.
+
+    The bias term depends on the intrinsic coherence itself, so it can be calculated iteratively
+    following the procedure from sec. 5 of Ingram 2019, MNRAS 489, 3927, which adjusts the bias
+    term until convergence is reached. This is done if ``adjust_bias`` is set to True. It is
+    typically a very small correction.
+
+    For errors, assumes **high powers, high coherence**. See eq. 8 of the paper.
+    Powers below 3 sigma above the noise level (P_noise / sqrt(MW)) are considered too low,
+    and the coherence will be set to NaN.
+
+    TODO: implement a more general treatment of the uncertainty.
+
+    Parameters
+    ----------
+    cross_power : complex `np.array`
+        cross spectrum
+    power1 : float `np.array`
+        sub-band periodogram
+    power2 : float `np.array`
+        reference-band periodogram
+    power1_noise : float
+        Poisson noise level of the sub-band periodogram
+    power2_noise : float
+        Poisson noise level of the reference-band periodogram
+    n_ave : int
+        number of intervals that have been averaged to obtain the input spectra
+
+    Other Parameters
+    ----------------
+    return_uncertainty : bool, default False
+        Whether to return the uncertainty on the coherence, calculated according to
+        Vaughan & Nowak 1997, ApJ 474, L43, eq. 8.
+    atol: float, default 0.01
+        The absolute tolerance for the convergence of the iterative procedure to adjust
+        the bias term.Only relevant if ``adjust_bias`` is True.
+
+    Returns
+    -------
+    coherence : float `np.array` or tuple of two float `np.array`
+        The intrinsic coherence values at all frequencies. It is 0 if the numerator
+        of the coherence calculation is negative, and NaN if the powers are too low compared
+        to the noise level.
+    uncertainty : float `np.array`, optional
+        The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak
+        1997, ApJ 474, L43, eq. 8.
+
+    """
+
+    if isinstance(cross_power, Iterable):
+        func = intrinsic_coherence_array
+    else:
+        func = intrinsic_coherence_single
+    coherence, uncertainty, flags = func(
+        cross_power,
+        power1,
+        power2,
+        power1_noise,
+        power2_noise,
+        n_ave,
+        return_uncertainty=True,
+        adjust_bias=adjust_bias,
+        atol=atol,
+    )
+    print(coherence, uncertainty, flags)
 
+    if np.any(flags):
+        warnings.warn(
+            "The iterative procedure to adjust the bias term did not converge after 40 "
+            "iterations. Consider rebinning the spectra to increase the signal-to-noise "
+            "ratio, or to use a more permissive tolerance (``atol`` parameter)."
+        )
     if np.any(np.isnan(coherence)):
         warnings.warn(
             "NaN values detected in intrinsic_coherence calculation. This happens when the powers "
diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py
index 07f4e8503..96c89a58a 100644
--- a/stingray/tests/test_fourier.py
+++ b/stingray/tests/test_fourier.py
@@ -300,7 +300,7 @@ def test_check_powers_for_intrinsic_coherence(self):
         pow2_noise = 5
         n_ave = np.array([1, 10, 10])
         # Only one power is below the threshold, due to the low number of averaged powers.
-        res = check_powers_for_intrinsic_coherence(pow1, pow2, pow1_noise, pow2_noise, n_ave)
+        res = check_powers_for_intrinsic_coherence(pow1, pow2, pow1_noise, pow2_noise, n_ave, 3.0)
         assert np.all(res == np.array([True, False, False]))
 
 
diff --git a/stingray/utils.py b/stingray/utils.py
index 113e2bd74..8c2f0174d 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -49,7 +49,8 @@
     from numba import jit
 
     HAS_NUMBA = True
-    from numba import njit, prange, vectorize, float32, float64, int32, int64, UniTuple
+    from numba import njit, prange, vectorize, float32, float64, int32, int64
+
     from numba.core.errors import NumbaValueError, NumbaNotImplementedError, TypingError
 except ImportError:
     warnings.warn(
@@ -78,7 +79,7 @@ def decorator(func, *a, **kw):
     def generic(*args, **kwargs):
         return None
 
-    float32 = float64 = int32 = int64 = UniTuple = generic
+    float32 = float64 = int32 = int64 = generic
 
     def prange(x):
         return range(x)

From 85165cf18050db62afd2b7d86bf8b09a69192fb0 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Sun, 1 Mar 2026 19:13:47 +0100
Subject: [PATCH 2/9] Fixes to intrinsic coherence functions

---
 stingray/fourier.py            | 189 ++++++++++++---------------------
 stingray/tests/test_fourier.py |  13 +++
 2 files changed, 78 insertions(+), 124 deletions(-)

diff --git a/stingray/fourier.py b/stingray/fourier.py
index 3542e90be..6f6a74ea0 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -1184,7 +1184,6 @@ def intrinsic_coherence_array(
     power1_noise,
     power2_noise,
     n_ave,
-    return_uncertainty=False,
     adjust_bias=False,
     atol=0.01,
 ):
@@ -1219,144 +1218,75 @@ def intrinsic_coherence_array(
     Parameters
     ----------
     cross_power : complex `np.array`
-        cross spectrum
+        Cross spectrum
     power1 : float `np.array`
-        sub-band periodogram
+        Sub-band periodogram. Must have the same shape as `cross_power`.
     power2 : float `np.array`
-        reference-band periodogram
-    power1_noise : float
-        Poisson noise level of the sub-band periodogram
-    power2_noise : float
-        Poisson noise level of the reference-band periodogram
-    n_ave : int
-        number of intervals that have been averaged to obtain the input spectra
+        Reference-band periodogram. Must have the same shape as `cross_power`.
+    power1_noise : float `np.array`
+        Poisson noise level of the sub-band periodogram. Can have a single value
+        or the same shape as `cross_power`.
+    power2_noise : float `np.array`
+        Poisson noise level of the reference-band periodogram. Can have a single value
+        or the same shape as `cross_power`.
+    n_ave : int `np.array`
+        number of intervals that have been averaged to obtain the input spectra. Can have
+        a single value or the same shape as `cross_power`.
 
     Other Parameters
     ----------------
-    return_uncertainty : bool, default False
-        Whether to return the uncertainty on the coherence, calculated according to
-        Vaughan & Nowak 1997, ApJ 474, L43, eq. 8.
     atol: float, default 0.01
         The absolute tolerance for the convergence of the iterative procedure to adjust
         the bias term.Only relevant if ``adjust_bias`` is True.
 
     Returns
     -------
-    coherence : float `np.array` or tuple of two float `np.array`
+    coherence : float `np.array`
         The intrinsic coherence values at all frequencies. It is 0 if the numerator
         of the coherence calculation is negative, and NaN if the powers are too low compared
         to the noise level.
-    uncertainty : float `np.array`, optional
+    uncertainty : float `np.array`
         The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak
         1997, ApJ 474, L43, eq. 8.
+    flags : bool `np.array`
+        Whether the bias term adjustment procedure converged (False) or not (True). Only
+        relevant if ``adjust_bias`` is True.
 
     """
 
     coherence = np.empty(cross_power.shape, dtype=np.float64)
     uncertainty = np.empty(cross_power.shape, dtype=np.float64)
-    flags = np.zeros(cross_power.shape, dtype=np.bool_)
+    no_conversion_flags = []
 
     for i in range(cross_power.size):
+        if np.size(power1_noise) == np.size(cross_power):
+            p1_n = float(power1_noise[i])
+        else:
+            p1_n = float(power1_noise[0])
+
+        if np.size(power2_noise) == np.size(cross_power):
+            p2_n = float(power2_noise[i])
+        else:
+            p2_n = float(power2_noise[0])
+
+        if np.size(n_ave) == np.size(cross_power):
+            n_a = int(n_ave[i])
+        else:
+            n_a = int(n_ave[0])
+
         if adjust_bias:
-            coherence[i], uncertainty[i], flags[i] = _intrinsic_coherence_with_adjusted_bias(
-                cross_power[i], power1[i], power2[i], power1_noise, power2_noise, n_ave, atol=atol
+            coherence[i], uncertainty[i], flag = _intrinsic_coherence_with_adjusted_bias(
+                cross_power[i], power1[i], power2[i], p1_n, p2_n, n_a, atol=atol
             )
+            if flag:
+                no_conversion_flags.append(i)
 
         else:
             coherence[i], uncertainty[i] = _intrinsic_coherence(
-                cross_power[i], power1[i], power2[i], power1_noise, power2_noise, n_ave
+                cross_power[i], power1[i], power2[i], p1_n, p2_n, n_a
             )
 
-    return coherence, uncertainty, flags
-
-
-@njit
-def intrinsic_coherence_single(
-    cross_power,
-    power1,
-    power2,
-    power1_noise,
-    power2_noise,
-    n_ave,
-    return_uncertainty=False,
-    adjust_bias=False,
-    atol=0.01,
-):
-    r"""
-    Intrinsic coherence estimations from cross and power spectra.
-
-    Vaughan & Nowak 1997, ApJ 474, L43
-
-    .. math::
-
-        \tilde{\gamma}^2(f) = \frac{|\langle \tilde{C}(f) \rangle|^2 - \tilde{b}^2}
-        {(\langle \tilde{P}_1(f) \rangle - \tilde{P}_{1, \rm noise})
-        (\langle \tilde{P}_2(f) \rangle - \tilde{P}_{2, \rm noise})}
-
-    where :math:`\tilde{b}^2` is the bias term that accounts for the contribution of Poisson
-    noise to the cross spectrum (see :func:`stingray.fourier.bias_term`), and tilde generally
-    indicates noisy measurements of the cross spectrum :math:`C` and the power spectra :math:`P_n`.
-    The terms :math:`\tilde{P}_{\rm n, \rm noise}` are the estimates of the contribution of
-    Poisson noise to the power spectra.
-
-    The bias term depends on the intrinsic coherence itself, so it can be calculated iteratively
-    following the procedure from sec. 5 of Ingram 2019, MNRAS 489, 3927, which adjusts the bias
-    term until convergence is reached. This is done if ``adjust_bias`` is set to True. It is
-    typically a very small correction.
-
-    For errors, assumes **high powers, high coherence**. See eq. 8 of the paper.
-    Powers below 3 sigma above the noise level (P_noise / sqrt(MW)) are considered too low,
-    and the coherence will be set to NaN.
-
-    TODO: implement a more general treatment of the uncertainty.
-
-    Parameters
-    ----------
-    cross_power : complex `np.array`
-        cross spectrum
-    power1 : float `np.array`
-        sub-band periodogram
-    power2 : float `np.array`
-        reference-band periodogram
-    power1_noise : float
-        Poisson noise level of the sub-band periodogram
-    power2_noise : float
-        Poisson noise level of the reference-band periodogram
-    n_ave : int
-        number of intervals that have been averaged to obtain the input spectra
-
-    Other Parameters
-    ----------------
-    return_uncertainty : bool, default False
-        Whether to return the uncertainty on the coherence, calculated according to
-        Vaughan & Nowak 1997, ApJ 474, L43, eq. 8.
-    atol: float, default 0.01
-        The absolute tolerance for the convergence of the iterative procedure to adjust
-        the bias term.Only relevant if ``adjust_bias`` is True.
-
-    Returns
-    -------
-    coherence : float `np.array` or tuple of two float `np.array`
-        The intrinsic coherence values at all frequencies. It is 0 if the numerator
-        of the coherence calculation is negative, and NaN if the powers are too low compared
-        to the noise level.
-    uncertainty : float `np.array`, optional
-        The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak
-        1997, ApJ 474, L43, eq. 8.
-
-    """
-
-    if adjust_bias:
-        coherence, uncertainty, flags = _intrinsic_coherence_with_adjusted_bias(
-            cross_power, power1, power2, power1_noise, power2_noise, n_ave, atol=atol
-        )
-    else:
-        coherence, uncertainty = _intrinsic_coherence(
-            cross_power, power1, power2, power1_noise, power2_noise, n_ave
-        )
-        flags = False
-
-    return coherence, uncertainty, flags
+    return coherence, uncertainty, no_conversion_flags
 
 
 def intrinsic_coherence(
@@ -1401,16 +1331,16 @@ def intrinsic_coherence(
     Parameters
     ----------
     cross_power : complex `np.array`
-        cross spectrum
+        Cross spectrum
     power1 : float `np.array`
-        sub-band periodogram
+        Sub-band periodogram. Must have the same shape as `cross_power`.
     power2 : float `np.array`
-        reference-band periodogram
-    power1_noise : float
+        Reference-band periodogram. Must have the same shape as `cross_power`.
+    power1_noise : float or float `np.array` of the same shape as `cross_power`
         Poisson noise level of the sub-band periodogram
-    power2_noise : float
+    power2_noise : float or float `np.array` of the same shape as `cross_power`
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or int `np.array` of the same shape as `cross_power`
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -1433,25 +1363,32 @@ def intrinsic_coherence(
         1997, ApJ 474, L43, eq. 8.
 
     """
-
-    if isinstance(cross_power, Iterable):
-        func = intrinsic_coherence_array
-    else:
-        func = intrinsic_coherence_single
-    coherence, uncertainty, flags = func(
+    single_power = False
+    if not isinstance(cross_power, np.ndarray):
+        cross_power = np.asanyarray([cross_power])
+        power1 = np.asanyarray([power1])
+        power2 = np.asanyarray([power2])
+        single_power = True
+
+    if not isinstance(power1_noise, Iterable):
+        power1_noise = np.asanyarray([power1_noise])
+    if not isinstance(power2_noise, Iterable):
+        power2_noise = np.asanyarray([power2_noise])
+    if not isinstance(n_ave, Iterable):
+        n_ave = np.asanyarray([n_ave])
+
+    coherence, uncertainty, flagged = intrinsic_coherence_array(
         cross_power,
         power1,
         power2,
         power1_noise,
         power2_noise,
         n_ave,
-        return_uncertainty=True,
         adjust_bias=adjust_bias,
         atol=atol,
     )
-    print(coherence, uncertainty, flags)
 
-    if np.any(flags):
+    if len(flagged) > 0:
         warnings.warn(
             "The iterative procedure to adjust the bias term did not converge after 40 "
             "iterations. Consider rebinning the spectra to increase the signal-to-noise "
@@ -1470,6 +1407,10 @@ def intrinsic_coherence(
             " well defined. "
         )
 
+    if single_power:
+        coherence = coherence[0]
+        uncertainty = uncertainty[0]
+
     if return_uncertainty:
         return coherence, uncertainty
     return coherence
diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py
index 96c89a58a..13d9c8c9b 100644
--- a/stingray/tests/test_fourier.py
+++ b/stingray/tests/test_fourier.py
@@ -224,8 +224,21 @@ def test_intrinsic_low_coh_iteration(self):
         assert np.allclose(coh, 0.22, atol=0.01)
 
         coh = intrinsic_coherence(C[0], P1[0], P2[0], 2, 2, 2, adjust_bias=True)
+        assert np.shape(coh) == ()
         assert np.isclose(coh, 0.22, atol=0.01)
 
+    def test_intrinsic_coh_all_arrays(self):
+        # bsq = (8 * 8 - 6 * 6) / 2. = 14 for coherence = 1, 36 for coherence = 0
+        # C**2 = 36, so that C**2 - bsq is positive, but coherence is less than 1.
+        # This will trigger the iteration.
+
+        C, P1, P2 = np.array([6, 6]), np.array([8, 8]), np.array([8, 8])
+        P1n = np.array([2, 2])
+        P2n = np.array([2, 2])
+        n_ave = np.array([2, 2])
+        coh = intrinsic_coherence(C, P1, P2, P1n, P2n, n_ave, adjust_bias=True)
+        assert np.allclose(coh, 0.22, atol=0.01)
+
     def test_intrinsic_low_coh_many_iteration(self):
         # bsq = (8 * 8 - 6 * 6) / 2. = 14 for coherence = 1, 36 for coherence = 0
         # C**2 > 36, so that C**2 - bsq is positive, but coherence is less than 1.

From ef8431adb4edc7b22c8a5c70edc373726761cc14 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Sun, 1 Mar 2026 19:17:43 +0100
Subject: [PATCH 3/9] Separate simple numba import from subimports (in order to
 detect API changes from missing the entire library)

---
 stingray/utils.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/stingray/utils.py b/stingray/utils.py
index 8c2f0174d..4c2f1d5a7 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -45,13 +45,11 @@
 
 # If numba is installed, import jit. Otherwise, define an empty decorator with
 # the same name.
+
 try:
     from numba import jit
 
     HAS_NUMBA = True
-    from numba import njit, prange, vectorize, float32, float64, int32, int64
-
-    from numba.core.errors import NumbaValueError, NumbaNotImplementedError, TypingError
 except ImportError:
     warnings.warn(
         "The recommended numba package is not installed. Some functionality might be slower."
@@ -85,6 +83,12 @@ def prange(x):
         return range(x)
 
 
+if HAS_NUMBA:
+    from numba import njit, prange, vectorize, float32, float64, int32, int64
+
+    from numba.core.errors import NumbaValueError, NumbaNotImplementedError, TypingError
+
+
 try:
     from tqdm import tqdm as show_progress
 except ImportError:

From cd520b7a0ccfff49a937e2ca5a936ce92bbd7a0b Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Sun, 1 Mar 2026 19:37:31 +0100
Subject: [PATCH 4/9] Add changelog

---
 docs/changes/970.bugfix.rst | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 docs/changes/970.bugfix.rst

diff --git a/docs/changes/970.bugfix.rst b/docs/changes/970.bugfix.rst
new file mode 100644
index 000000000..5e9cd44a9
--- /dev/null
+++ b/docs/changes/970.bugfix.rst
@@ -0,0 +1 @@
+Fix bug in intrinsic coherence functions that made Numba unavailable

From 686b8b94ac040972b525cb6576af867613e19fb3 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Sun, 1 Mar 2026 19:47:33 +0100
Subject: [PATCH 5/9] Be kind with int values recorded as floats

---
 stingray/fourier.py | 58 ++++++++++++++++++++++++++++++---------------
 1 file changed, 39 insertions(+), 19 deletions(-)

diff --git a/stingray/fourier.py b/stingray/fourier.py
index 6f6a74ea0..a0132a0d8 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -23,8 +23,6 @@
     rebin_data,
     njit,
     vectorize,
-    float64,
-    int64,
 )
 
 
@@ -705,7 +703,13 @@ def unnormalize_periodograms(
     raise ValueError("Unrecognized power type")
 
 
-@vectorize([float64(float64, float64, float64, float64, int64, float64)], nopython=True)
+@vectorize(
+    [
+        "float64(float64, float64, float64, float64, int64, float64)",
+        "float64(float64, float64, float64, float64, float64, float64)",
+    ],
+    nopython=True,
+)
 def _bias_term(power1, power2, power1_noise, power2_noise, n_ave, input_intrinsic_coherence):
     if n_ave > 500:
         return 0.0
@@ -737,7 +741,7 @@ def bias_term(power1, power2, power1_noise, power2_noise, n_ave, intrinsic_coher
         Poisson noise level of the sub-band periodogram
     power2_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -754,7 +758,7 @@ def bias_term(power1, power2, power1_noise, power2_noise, n_ave, intrinsic_coher
     return _bias_term(power1, power2, power1_noise, power2_noise, n_ave, intrinsic_coherence)
 
 
-@vectorize([float64(float64, float64, float64)], nopython=True)
+@vectorize(["float64(float64, float64, float64)"], nopython=True)
 def _apply_low_lim_to_coherence_uncertainty(coherence, uncertainty, min_uncertainty):
     """
     Apply a low limit to the uncertainty on the coherence, to avoid zero or negative uncertainties.
@@ -807,7 +811,10 @@ def _apply_low_lim_to_coherence_uncertainty(coherence, uncertainty, min_uncertai
 
 
 @vectorize(
-    ["float64(complex128, float64, float64, float64, float64, int64, float64)"],
+    [
+        "float64(complex128, float64, float64, float64, float64, int64, float64)",
+        "float64(complex128, float64, float64, float64, float64, float64, float64)",
+    ],
     nopython=True,
 )
 def _raw_coherence(
@@ -836,7 +843,7 @@ def _raw_coherence(
         Poisson noise level of the sub-band periodogram
     power2_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -899,7 +906,7 @@ def raw_coherence(
         Poisson noise level of the sub-band periodogram
     power2_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -918,6 +925,16 @@ def raw_coherence(
     ----------
     .. [#] https://doi.org/10.1093/mnras/stz2409
     """
+    print(
+        cross_power,
+        power1,
+        power2,
+        power1_noise,
+        power2_noise,
+        n_ave,
+        intrinsic_coherence,
+    )
+
     coherence = _raw_coherence(
         cross_power,
         power1,
@@ -959,7 +976,7 @@ def _intrinsic_coherence_uncertainties(
         The Poisson noise level of the second band periodogram
     coherence : float
         The intrinsic coherence calculated according to the _intrinsic_coherence function.
-    n_ave : int
+    n_ave : int or float
         The number of intervals that have been averaged to obtain the input spectra
     """
     # Terms from VN97, eq. 8, for the uncertainty on the coherence.
@@ -971,7 +988,10 @@ def _intrinsic_coherence_uncertainties(
 
 
 @vectorize(
-    ["bool(float64, float64, float64, float64, int64, float64)"],
+    [
+        "bool(float64, float64, float64, float64, int64, float64)",
+        "bool(float64, float64, float64, float64, float64, float64)",
+    ],
     nopython=True,
 )
 def check_powers_for_intrinsic_coherence(
@@ -989,7 +1009,7 @@ def check_powers_for_intrinsic_coherence(
         Poisson noise level of the sub-band periodogram
     power2_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -1044,7 +1064,7 @@ def _intrinsic_coherence(
         Poisson noise level of the sub-band periodogram
     power2_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -1119,7 +1139,7 @@ def _intrinsic_coherence_with_adjusted_bias(
         Poisson noise level of the sub-band periodogram
     power2_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -1229,7 +1249,7 @@ def intrinsic_coherence_array(
     power2_noise : float `np.array`
         Poisson noise level of the reference-band periodogram. Can have a single value
         or the same shape as `cross_power`.
-    n_ave : int `np.array`
+    n_ave : int or float `np.array`
         number of intervals that have been averaged to obtain the input spectra. Can have
         a single value or the same shape as `cross_power`.
 
@@ -1340,7 +1360,7 @@ def intrinsic_coherence(
         Poisson noise level of the sub-band periodogram
     power2_noise : float or float `np.array` of the same shape as `cross_power`
         Poisson noise level of the reference-band periodogram
-    n_ave : int or int `np.array` of the same shape as `cross_power`
+    n_ave : int or float or int `np.array` of the same shape as `cross_power`
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -1436,7 +1456,7 @@ def estimate_intrinsic_coherence(cross_power, power1, power2, power1_noise, powe
         Poisson noise level of the sub-band periodogram
     power2_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Returns
@@ -1729,7 +1749,7 @@ def error_on_averaged_cross_spectrum(
         Poisson noise level of the sub-band periodogram
     ref_power_noise : float
         Poisson noise level of the reference-band periodogram
-    n_ave : int
+    n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
 
     Other Parameters
@@ -2957,7 +2977,7 @@ def avg_pds_from_timeseries(
         The normalized periodogram powers
     n_bin : int
         the number of bins in the light curves used in each segment
-    n_ave : int
+    n_ave : int or float
         the number of averaged periodograms
     mean : float
         the mean flux
@@ -3071,7 +3091,7 @@ def avg_cs_from_timeseries(
         The normalized periodogram powers
     n_bin : int
         the number of bins in the light curves used in each segment
-    n_ave : int
+    n_ave : int or float
         the number of averaged periodograms
     """
 

From 14aef0e31031836a6ab34116ac2f5d1c7510f886 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Sun, 1 Mar 2026 19:48:37 +0100
Subject: [PATCH 6/9] Cleanup

---
 stingray/fourier.py | 23 +++++++----------------
 1 file changed, 7 insertions(+), 16 deletions(-)

diff --git a/stingray/fourier.py b/stingray/fourier.py
index a0132a0d8..918f6d779 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -925,15 +925,6 @@ def raw_coherence(
     ----------
     .. [#] https://doi.org/10.1093/mnras/stz2409
     """
-    print(
-        cross_power,
-        power1,
-        power2,
-        power1_noise,
-        power2_noise,
-        n_ave,
-        intrinsic_coherence,
-    )
 
     coherence = _raw_coherence(
         cross_power,
@@ -1011,11 +1002,11 @@ def check_powers_for_intrinsic_coherence(
         Poisson noise level of the reference-band periodogram
     n_ave : int or float
         number of intervals that have been averaged to obtain the input spectra
+    threshold : float
+        The threshold in sigma above the noise level for the powers to be considered "high".
 
     Other Parameters
     ----------------
-    threshold : float, default 3
-        The threshold in sigma above the noise level for the powers to be considered "high".
 
     Returns
     -------
@@ -1157,8 +1148,8 @@ def _intrinsic_coherence_with_adjusted_bias(
     uncertainty : float `np.array`, optional
         The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak
         1997, ApJ 474, L43, eq. 8.
-    non_converged_flag : bool `np.array`
-        Whether the bias term adjustment procedure converged (True) or not (False).
+    not_converged_flag : bool `np.array`
+        Whether the bias term adjustment procedure converged (False) or not (True).
 
     """
 
@@ -1197,7 +1188,7 @@ def _intrinsic_coherence_with_adjusted_bias(
 
 
 @njit
-def intrinsic_coherence_array(
+def _intrinsic_coherence_array(
     cross_power,
     power1,
     power2,
@@ -1257,7 +1248,7 @@ def intrinsic_coherence_array(
     ----------------
     atol: float, default 0.01
         The absolute tolerance for the convergence of the iterative procedure to adjust
-        the bias term.Only relevant if ``adjust_bias`` is True.
+        the bias term. Only relevant if ``adjust_bias`` is True.
 
     Returns
     -------
@@ -1397,7 +1388,7 @@ def intrinsic_coherence(
     if not isinstance(n_ave, Iterable):
         n_ave = np.asanyarray([n_ave])
 
-    coherence, uncertainty, flagged = intrinsic_coherence_array(
+    coherence, uncertainty, flagged = _intrinsic_coherence_array(
         cross_power,
         power1,
         power2,

From 81afd8b84048cafebbf3332653f403dac412d208 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Sun, 1 Mar 2026 21:38:20 +0100
Subject: [PATCH 7/9] Fix docstrings [docs only]

---
 stingray/fourier.py | 79 +++++++++++++++++++++++----------------------
 1 file changed, 40 insertions(+), 39 deletions(-)

diff --git a/stingray/fourier.py b/stingray/fourier.py
index 918f6d779..622b9f9d2 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -555,27 +555,27 @@ def normalize_periodograms(
 
     Parameters
     ----------
-    unnorm_power: numpy.ndarray
+    unnorm_power : numpy.ndarray
         The unnormalized cross spectrum.
 
-    dt: float
+    dt : float
         The sampling time of the light curve
 
-    n_bin: int
+    n_bin : int
         The number of bins in the light curve
 
     Other parameters
     ----------------
-    mean_flux: float
+    mean_flux : float
         The mean of the light curve used to calculate the powers
         (If a cross spectrum, the geometrical mean of the light
         curves in the two channels). Only relevant for "frac" normalization
 
-    n_ph: int or float
+    n_ph : int or float
         The number of counts in the light curve used to calculate
         the unnormalized periodogram. Only relevant for Leahy normalization.
 
-    variance: float
+    variance : float
         The average variance of the measurements in light curve (if a cross
         spectrum,  the geometrical mean of the variances in the two channels).
         **NOT** the variance of the light curve, but of each flux measurement
@@ -595,7 +595,7 @@ def normalize_periodograms(
 
     Returns
     -------
-    power: numpy.nd.array
+    power : numpy.nd.array
         The normalized co-spectrum (real part of the cross spectrum). For
         'none' normalization, imaginary part is returned as well.
     """
@@ -633,27 +633,27 @@ def unnormalize_periodograms(
 
     Parameters
     ----------
-    norm_power: numpy.ndarray
+    norm_power : numpy.ndarray
         The normalized cross-spectrum or poisson noise
 
-    dt: float
+    dt : float
         The sampling time of the light curve
 
-    n_bin: int
+    n_bin : int
         The number of bins in the light curve
 
     Other parameters
     ----------------
-    mean_flux: float
+    mean_flux : float
         The mean of the light curve used to calculate the powers
         (If a cross spectrum, the geometrical mean of the light
         curves in the two channels). Only relevant for "frac" normalization
 
-    n_ph: int or float
+    n_ph : int or float
         The number of counts in the light curve used to calculate
         the unnormalized periodogram. Only relevant for Leahy normalization.
 
-    variance: float
+    variance : float
         The average variance of the measurements in light curve (if a cross
         spectrum,  the geometrical mean of the variances in the two channels).
         **NOT** the variance of the light curve, but of each flux measurement
@@ -673,7 +673,7 @@ def unnormalize_periodograms(
 
     Returns
     -------
-    power: numpy.nd.array
+    power : numpy.nd.array
         The normalized co-spectrum (real part of the cross spectrum). For
         'none' normalization, imaginary part is returned as well.
     """
@@ -1005,9 +1005,6 @@ def check_powers_for_intrinsic_coherence(
     threshold : float
         The threshold in sigma above the noise level for the powers to be considered "high".
 
-    Other Parameters
-    ----------------
-
     Returns
     -------
     low_power1 : bool `np.array`
@@ -1135,7 +1132,7 @@ def _intrinsic_coherence_with_adjusted_bias(
 
     Other Parameters
     ----------------
-    atol: float, default 0.01
+    atol : float, default 0.01
         The absolute tolerance for the convergence of the iterative procedure to adjust
         the bias term.
 
@@ -1246,7 +1243,11 @@ def _intrinsic_coherence_array(
 
     Other Parameters
     ----------------
-    atol: float, default 0.01
+    adjust_bias : bool, default False
+        Whether to adjust the bias term iteratively following the procedure from sec. 5 of
+        Ingram 2019, MNRAS 489, 3927. It is typically a very small correction, but it can be
+        relevant for low coherence values and/or low number of averaged intervals
+    atol : float, default 0.01
         The absolute tolerance for the convergence of the iterative procedure to adjust
         the bias term. Only relevant if ``adjust_bias`` is True.
 
@@ -1259,9 +1260,9 @@ def _intrinsic_coherence_array(
     uncertainty : float `np.array`
         The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak
         1997, ApJ 474, L43, eq. 8.
-    flags : bool `np.array`
-        Whether the bias term adjustment procedure converged (False) or not (True). Only
-        relevant if ``adjust_bias`` is True.
+    no_conversion_flags : `list` of int
+        List of indices of ``coherence`` corresponding to frequencies where the bias term
+        adjustment procedure did not converge. Only relevant if ``adjust_bias`` is True.
 
     """
 
@@ -1359,7 +1360,7 @@ def intrinsic_coherence(
     return_uncertainty : bool, default False
         Whether to return the uncertainty on the coherence, calculated according to
         Vaughan & Nowak 1997, ApJ 474, L43, eq. 8.
-    atol: float, default 0.01
+    atol : float, default 0.01
         The absolute tolerance for the convergence of the iterative procedure to adjust
         the bias term.Only relevant if ``adjust_bias`` is True.
 
@@ -1480,13 +1481,13 @@ def get_rms_from_rms_norm_periodogram(power_sqrms, poisson_noise_sqrms, df, M, l
 
     Parameters
     ----------
-    power_sqrms: array-like
+    power_sqrms : array-like
         Powers, in units of fractional rms ($(rms/mean)^2 Hz{-1}$)
-    poisson_noise_sqrms: float
+    poisson_noise_sqrms : float
         Poisson noise level, in units of fractional rms ($(rms/mean)^2 Hz{-1}$
-    df: float or ``np.array``, same dimension of ``power_sqrms``
+    df : float or ``np.array``, same dimension of ``power_sqrms``
         The frequency resolution of each power
-    M: int or ``np.array``, same dimension of ``power_sqrms``
+    M : int or ``np.array``, same dimension of ``power_sqrms``
         The number of powers averaged to obtain each value of power.
 
     Other Parameters
@@ -1648,33 +1649,33 @@ def rms_calculation(
 
     Parameters
     ----------
-    unnorm_powers: array of float
+    unnorm_powers : array of float
         unnormalised power or cross spectrum, the array has already been
         filtered for the given frequency range
 
-    min_freq: float
+    min_freq : float
         The lower frequency bound for the calculation (from the freq grid).
 
-    max_freq: float
+    max_freq : float
         The upper frequency bound for the calculation (from the freq grid).
 
-    nphots: float
+    nphots : float
         Number of photons for the full power or cross spectrum
 
-    T: float
+    T : float
         Time length of the light curve
 
-    M_freq: scalar or array of float
+    M_freq : scalar or array of float
         If scalar, it is the number of segments in the AveragedCrossspectrum
         If array, it is the number of segments times the rebinning sample
         in the given frequency range.
 
-    K_freq: scalar or array of float
+    K_freq : scalar or array of float
         If scalar, the power or cross spectrum is not rebinned (K_freq = 1)
         If array,  the power or cross spectrum is rebinned and it is the
         rebinned sample in the given frequency range.
 
-    freq_bins: integer
+    freq_bins : integer
         if the cross or power spectrum is rebinned freq_bins = 1,
         if it NOT rebinned freq_bins is the number of frequency bins
         in the given frequency range.
@@ -1684,16 +1685,16 @@ def rms_calculation(
 
     Other parameters
     ----------------
-    deadtime: float
+    deadtime : float
         Deadtime of the instrument
 
     Returns
     -------
-    rms: float
+    rms : float
         The fractional rms amplitude contained between ``min_freq`` and
         ``max_freq``.
 
-    rms_err: float
+    rms_err : float
         The error on the fractional rms amplitude.
 
     """
@@ -1827,7 +1828,7 @@ def cross_to_covariance(cross_power, ref_power, ref_power_noise, delta_nu):
 
     Returns
     -------
-    covariance: complex `np.array`
+    covariance : complex `np.array`
         The cross spectrum, normalized as a covariance.
 
     """

From 52a8aaecaf77261dbf31836ed55ccbbc1ca077aa Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Mon, 2 Mar 2026 18:39:14 +0100
Subject: [PATCH 8/9] Edit changelog

---
 CHANGELOG.rst               | 9 +++++++++
 docs/changes/970.bugfix.rst | 1 -
 2 files changed, 9 insertions(+), 1 deletion(-)
 delete mode 100644 docs/changes/970.bugfix.rst

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 7cf874c83..85da99db8 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,3 +1,12 @@
+v2.3.1 (2026-03-02)
+-------------------
+
+Bug Fixes
+^^^^^^^^^
+
+- Fix bug in intrinsic coherence functions that made Numba unavailable (`#970 <https://github.com/StingraySoftware/stingray/pull/970>`__)
+
+
 v2.3 (2026-02-27)
 -----------------
 
diff --git a/docs/changes/970.bugfix.rst b/docs/changes/970.bugfix.rst
deleted file mode 100644
index 5e9cd44a9..000000000
--- a/docs/changes/970.bugfix.rst
+++ /dev/null
@@ -1 +0,0 @@
-Fix bug in intrinsic coherence functions that made Numba unavailable

From c5ec171c03510cf4fb87cf6a3960a36e2f9efa44 Mon Sep 17 00:00:00 2001
From: Matteo Bachetti <matteo@matteobachetti.it>
Date: Mon, 2 Mar 2026 18:40:57 +0100
Subject: [PATCH 9/9] Cleanup

---
 stingray/fourier.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/stingray/fourier.py b/stingray/fourier.py
index 622b9f9d2..00de2eae3 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -38,7 +38,6 @@
     "unnormalize_periodograms",
     "bias_term",
     "raw_coherence",
-    "coherence",
     "intrinsic_coherence",
 ]
 
@@ -1848,8 +1847,6 @@ def _which_segment_idx_fun(binned=False, dt=None):
     """
     # Make function interface equal (fluxes gets ignored)
     if not binned:
-        fun = generate_indices_of_segment_boundaries_unbinned
-
         # Define a new function, make sure that, by default, the sort check
         # is disabled.
         def fun(*args, **kwargs):