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 `__) + + v2.3 (2026-02-27) ----------------- diff --git a/stingray/fourier.py b/stingray/fourier.py index 723408de5..00de2eae3 100644 --- a/stingray/fourier.py +++ b/stingray/fourier.py @@ -23,9 +23,6 @@ rebin_data, njit, vectorize, - float64, - int64, - UniTuple, ) @@ -41,7 +38,6 @@ "unnormalize_periodograms", "bias_term", "raw_coherence", - "coherence", "intrinsic_coherence", ] @@ -558,27 +554,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 @@ -598,7 +594,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. """ @@ -636,27 +632,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 @@ -676,7 +672,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. """ @@ -706,7 +702,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 @@ -738,7 +740,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 @@ -755,7 +757,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. @@ -809,7 +811,8 @@ 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)", + "float64(complex128, float64, float64, float64, float64, float64, float64)", ], nopython=True, ) @@ -820,7 +823,7 @@ def _raw_coherence( power1_noise, power2_noise, n_ave, - intrinsic_coherence=1, + intrinsic_coherence, ): """ Raw coherence estimations from cross and power spectra. @@ -839,7 +842,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 @@ -854,20 +857,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 +876,7 @@ def raw_coherence( power1_noise, power2_noise, n_ave, - intrinsic_coherence=1, + intrinsic_coherence=1.0, return_uncertainty=False, ): r""" @@ -906,7 +905,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 @@ -925,16 +924,22 @@ 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 @@ -961,7 +966,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. @@ -973,11 +978,14 @@ def _intrinsic_coherence_uncertainties( @vectorize( - ["UniTuple(float64[:, :], 2)(float64[:, :], float64[:, :], float64, float64, int64)"], + [ + "bool(float64, float64, float64, float64, int64, float64)", + "bool(float64, float64, float64, float64, float64, 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. @@ -991,12 +999,9 @@ 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 - ---------------- - threshold : float, default 3 + threshold : float The threshold in sigma above the noise level for the powers to be considered "high". Returns @@ -1009,15 +1014,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, @@ -1049,7 +1051,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 @@ -1069,7 +1071,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 +1094,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, @@ -1127,12 +1126,12 @@ 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 ---------------- - 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. @@ -1145,12 +1144,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. + not_converged_flag : bool `np.array` + Whether the bias term adjustment procedure converged (False) or not (True). """ - 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,32 +1160,144 @@ 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 + + +@njit +def _intrinsic_coherence_array( + cross_power, + power1, + power2, + power1_noise, + power2_noise, + n_ave, + 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. Must have the same shape as `cross_power`. + power2 : float `np.array` + 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 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`. + + Other Parameters + ---------------- + 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. + + Returns + ------- + 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` + The uncertainty on the intrinsic coherence, calculated according to Vaughan & Nowak + 1997, ApJ 474, L43, eq. 8. + 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. + + """ + + coherence = np.empty(cross_power.shape, dtype=np.float64) + uncertainty = np.empty(cross_power.shape, dtype=np.float64) + 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], 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], p1_n, p2_n, n_a + ) + + return coherence, uncertainty, no_conversion_flags def intrinsic_coherence( @@ -1229,16 +1342,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 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 @@ -1246,7 +1359,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. @@ -1261,16 +1374,37 @@ def intrinsic_coherence( 1997, ApJ 474, L43, eq. 8. """ - if adjust_bias: - result = _intrinsic_coherence_with_adjusted_bias( - cross_power, power1, power2, power1_noise, power2_noise, n_ave, atol=atol - ) - else: - result = _intrinsic_coherence( - cross_power, power1, power2, power1_noise, power2_noise, n_ave - ) - coherence, uncertainty = result + 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, + adjust_bias=adjust_bias, + atol=atol, + ) + 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 " + "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 " @@ -1284,6 +1418,10 @@ def intrinsic_coherence( " well defined. " ) + if single_power: + coherence = coherence[0] + uncertainty = uncertainty[0] + if return_uncertainty: return coherence, uncertainty return coherence @@ -1309,7 +1447,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 @@ -1342,13 +1480,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 @@ -1510,33 +1648,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. @@ -1546,16 +1684,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. """ @@ -1602,7 +1740,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 @@ -1689,7 +1827,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. """ @@ -1709,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): @@ -2830,7 +2966,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 @@ -2944,7 +3080,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 """ diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py index 07f4e8503..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. @@ -300,7 +313,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..4c2f1d5a7 100644 --- a/stingray/utils.py +++ b/stingray/utils.py @@ -45,12 +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, UniTuple - from numba.core.errors import NumbaValueError, NumbaNotImplementedError, TypingError except ImportError: warnings.warn( "The recommended numba package is not installed. Some functionality might be slower." @@ -78,12 +77,18 @@ 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) +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: