diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 0d8b872f7..5acef542f 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,21 @@ Changelog ========= +v1.1 (unreleased) +----------------- +Bug fixes +^^^^^^^^^ +- IMPORTANT: Fixed sign of time lags, which were calculated using the interest band as the reference. +- Fixed an issue when the fractional exposure in FITS light curves is slightly >1 (as sometimes happens in NICER data) + +New +^^^ +- Implemented the ``bexvar`` variability estimation method for light curves. + +Improvements +^^^^^^^^^^^^ +- A less confusing default value of segment_size in Z searches + v1.0 (2022-03-29) --------------------- TL,DR: these things will break your code with v1.0: diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py index 147f7a511..4d636f3b7 100644 --- a/stingray/crossspectrum.py +++ b/stingray/crossspectrum.py @@ -503,9 +503,9 @@ class Crossspectrum(StingrayObject): ---------------- gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...] Good Time intervals. Defaults to the common GTIs from the two input - objects. Could throw errors if these GTIs have overlaps with the input - `Lightcurve` GTIs! If you're getting errors regarding your GTIs, don't - use this and only give GTIs to the `Lightcurve` objects before making + objects. Could throw errors if these GTIs have overlaps with the input + `Lightcurve` GTIs! If you're getting errors regarding your GTIs, don't + use this and only give GTIs to the `Lightcurve` objects before making the cross spectrum. lc1: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects @@ -928,7 +928,7 @@ def _fourier_cross(self, lc1, lc2, fullspec=False): fourier_2 = fft(lc2.counts) # do Fourier transform 2 freqs = fftfreq(lc1.n, lc1.dt) - cross = np.multiply(fourier_1, np.conj(fourier_2)) + cross = np.multiply(fourier_2, np.conj(fourier_1)) if fullspec is True: return freqs, cross @@ -1722,13 +1722,13 @@ class AveragedCrossspectrum(Crossspectrum): Parameters ---------- data1: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects OR :class:`stingray.EventList` object - A light curve from which to compute the cross spectrum. In some cases, - this would be the light curve of the wavelength/energy/frequency band + A light curve from which to compute the cross spectrum. In some cases, + this would be the light curve of the wavelength/energy/frequency band of interest. data2: :class:`stingray.Lightcurve`object OR iterable of :class:`stingray.Lightcurve` objects OR :class:`stingray.EventList` object - A second light curve to use in the cross spectrum. In some cases, this - would be the wavelength/energy/frequency reference band to compare the + A second light curve to use in the cross spectrum. In some cases, this + would be the wavelength/energy/frequency reference band to compare the band of interest with. segment_size: float @@ -1745,9 +1745,9 @@ class AveragedCrossspectrum(Crossspectrum): ---------------- gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...] Good Time intervals. Defaults to the common GTIs from the two input - objects. Could throw errors if these GTIs have overlaps with the - input object GTIs! If you're getting errors regarding your GTIs, - don't use this and only give GTIs to the input objects before + objects. Could throw errors if these GTIs have overlaps with the + input object GTIs! If you're getting errors regarding your GTIs, + don't use this and only give GTIs to the input objects before making the cross spectrum. dt : float @@ -1794,16 +1794,16 @@ class AveragedCrossspectrum(Crossspectrum): after the average. legacy: bool - Use the legacy machinery of `AveragedCrossspectrum`. This might be - useful to compare with old results, and is also needed to use light - curve lists as an input, to conserve the spectra of each segment, or + Use the legacy machinery of `AveragedCrossspectrum`. This might be + useful to compare with old results, and is also needed to use light + curve lists as an input, to conserve the spectra of each segment, or to use the large_data option. gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...] Good Time intervals. Defaults to the common GTIs from the two input - objects. Could throw errors if these GTIs have overlaps with the - input object GTIs! If you're getting errors regarding your GTIs, - don't use this and only give GTIs to the input objects before + objects. Could throw errors if these GTIs have overlaps with the + input object GTIs! If you're getting errors regarding your GTIs, + don't use this and only give GTIs to the input objects before making the cross spectrum. Attributes @@ -1818,8 +1818,8 @@ class AveragedCrossspectrum(Crossspectrum): The uncertainties of ``power``. An approximation for each bin given by ``power_err= power/sqrt(m)``. Where ``m`` is the number of power averaged in each bin (by frequency - binning, or averaging power spectra of segments of a light curve). - Note that for a single realization (``m=1``) the error is equal to the + binning, or averaging power spectra of segments of a light curve). + Note that for a single realization (``m=1``) the error is equal to the power. df: float diff --git a/stingray/fourier.py b/stingray/fourier.py index 0413a6336..358ae1573 100644 --- a/stingray/fourier.py +++ b/stingray/fourier.py @@ -526,7 +526,7 @@ def bias_term(power1, power2, power1_noise, power2_noise, n_ave, bias : float `np.array`, same shape as ``power1`` and ``power2`` The bias term """ - if n_ave > 500: + if (isinstance(n_ave, Iterable) and np.all(n_ave > 500)) or (not isinstance(n_ave, Iterable) and n_ave > 500): return 0. * power1 bsq = power1 * power2 - intrinsic_coherence * (power1 - power1_noise) * \ (power2 - power2_noise) @@ -571,6 +571,7 @@ def raw_coherence(cross_power, power1, power2, power1_noise, power2_noise, if isinstance(num, Iterable): num[num < 0] = (cross_power * np.conj(cross_power)).real[num < 0] elif num < 0: + warnings.warn("Negative numerator in raw_coherence calculation. Setting bias term to 0") num = (cross_power * np.conj(cross_power)).real den = power1 * power2 return num / den @@ -663,7 +664,7 @@ def estimate_intrinsic_coherence(cross_power, power1, power2, power1_noise, def error_on_averaged_cross_spectrum(cross_power, seg_power, ref_power, n_ave, seg_power_noise, ref_power_noise, - common_ref="False"): + common_ref=False): """ Error on cross spectral quantities, From Ingram 2019. @@ -702,7 +703,7 @@ def error_on_averaged_cross_spectrum(cross_power, seg_power, ref_power, n_ave, Error on the modulus of the cross spectrum """ - if n_ave < 30: + if (not isinstance(n_ave, Iterable) and n_ave < 30) or (isinstance(n_ave, Iterable) and np.any(n_ave) < 30): warnings.warn("n_ave is below 30. Please note that the error bars " "on the quantities derived from the cross spectrum " "are only reliable for a large number of averaged " @@ -1160,7 +1161,7 @@ def avg_cs_from_iterables_quick( ft2 = fft(flux2) # Calculate the unnormalized cross spectrum - unnorm_power = ft1 * ft2.conj() + unnorm_power = ft1.conj() * ft2 # Accumulate the sum to calculate the total mean of the lc sum_of_photons1 += n_ph1 @@ -1358,7 +1359,7 @@ def local_show_progress(a): n_ph = np.sqrt(n_ph1 * n_ph2) # Calculate the unnormalized cross spectrum - unnorm_power = ft1 * ft2.conj() + unnorm_power = ft1.conj() * ft2 unnorm_pd1 = unnorm_pd2 = 0 # If requested, calculate the auxiliary PDSs diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py index 4c86d5104..eec68c958 100644 --- a/stingray/tests/test_crossspectrum.py +++ b/stingray/tests/test_crossspectrum.py @@ -1219,6 +1219,9 @@ def test_timelag(self): test_lc1.counts, err_dist=test_lc1.err_dist, dt=dt) + # The second light curve is delayed by two bins. + # The time lag should be -2 * dt, because this will + # become the reference band in AveragedCrossspectrum test_lc2 = Lightcurve(test_lc1.time, np.array(np.roll(test_lc1.counts, 2)), err_dist=test_lc1.err_dist, @@ -1230,7 +1233,9 @@ def test_timelag(self): time_lag, time_lag_err = cs.time_lag() - assert np.all(np.abs(time_lag[:6] - 0.1) < 3 * time_lag_err[:6]) + # The actual measured time lag will be half that for AveragedCrosspectrum + measured_lag = -dt + assert np.all(np.abs(time_lag[:6] - measured_lag) < 3 * time_lag_err[:6]) def test_errorbars_legacy(self): time = np.arange(10000) * 0.1 diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py index bb3c22925..cabb124b9 100644 --- a/stingray/tests/test_fourier.py +++ b/stingray/tests/test_fourier.py @@ -64,7 +64,7 @@ def setup_class(cls): good = (freq > 0) & (freq < 0.1) ft1, ft2 = ft1[good], ft2[good] cls.cross = normalize_periodograms( - ft1 * ft2.conj(), dt, cls.N, mean, norm="abs", power_type="all") + ft1.conj() * ft2, dt, cls.N, mean, norm="abs", power_type="all") cls.pds1 = normalize_periodograms( ft1 * ft1.conj(), dt, cls.N, mean, norm="abs", power_type="real") cls.pds2 = normalize_periodograms( diff --git a/stingray/tests/test_varenergyspectrum.py b/stingray/tests/test_varenergyspectrum.py index a0a78f212..3a79c1e04 100644 --- a/stingray/tests/test_varenergyspectrum.py +++ b/stingray/tests/test_varenergyspectrum.py @@ -358,35 +358,36 @@ def setup_class(cls): data = np.load(os.path.join(datadir, "sample_variable_lc.npy")) flux = data times = np.arange(data.size) * dt - maxfreq = 0.25 / cls.time_lag + maxfreq = 0.15 / cls.time_lag roll_amount = int(cls.time_lag // dt) good = slice(roll_amount, roll_amount + int(200 // dt)) + # When rolling, a positive delay is introduced rolled_flux = np.array(np.roll(flux, roll_amount)) times, flux, rolled_flux = times[good], flux[good], rolled_flux[good] length = times[-1] - times[0] - test_lc1 = Lightcurve(times, flux, err_dist="gauss", dt=dt, skip_checks=True) - test_lc2 = Lightcurve(test_lc1.time, rolled_flux, err_dist=test_lc1.err_dist, dt=dt) - test_ev1, test_ev2 = EventList(), EventList() - test_ev1.simulate_times(test_lc1) - test_ev2.simulate_times(test_lc2) + test_ref = Lightcurve(times, flux, err_dist="gauss", dt=dt, skip_checks=True) + test_sub = Lightcurve(test_ref.time, rolled_flux, err_dist=test_ref.err_dist, dt=dt) + test_ref_ev, test_sub_ev = EventList(), EventList() + test_ref_ev.simulate_times(test_ref) + test_sub_ev.simulate_times(test_sub) - test_ev1.energy = np.random.uniform(0.3, 9, len(test_ev1.time)) - test_ev2.energy = np.random.uniform(9, 12, len(test_ev2.time)) + test_sub_ev.energy = np.random.uniform(0.3, 9, len(test_sub_ev.time)) + test_ref_ev.energy = np.random.uniform(9, 12, len(test_ref_ev.time)) cls.lag = LagEnergySpectrum( - test_ev1, - freq_interval=[0, maxfreq], + test_sub_ev, + freq_interval=[maxfreq / 2, maxfreq], energy_spec=(0.3, 9, 1, "lin"), ref_band=[9, 12], bin_time=dt / 2, segment_size=length, - events2=test_ev2, + events2=test_ref_ev, ) # Make single event list - test_ev = test_ev1.join(test_ev2) + test_ev = test_sub_ev.join(test_ref_ev) cls.lag_same = LagEnergySpectrum( test_ev, diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py index e5060a56f..ee288df90 100644 --- a/stingray/varenergyspectrum.py +++ b/stingray/varenergyspectrum.py @@ -875,9 +875,13 @@ def _spectrum_function(self): Cmean, mean_sub_power, mean_ref_power, m_tot, sub_power_noise, ref_power_noise, common_ref=common_ref ) - lag = np.mean((np.angle(cross[good]) / (2 * np.pi * freq[good]))) + # The frequency of these lags is measured from the *weighted* mean of the frequencies + # in the cross spectrum. The weight is just the absolute value of the CS + csabs = np.abs(cross[good]) + fmean = np.sum(freq[good] * csabs) / np.sum(csabs) + lag = np.angle(Cmean) / (2 * np.pi * fmean) - lag_e = phi_e / (2 * np.pi * f) + lag_e = phi_e / (2 * np.pi * fmean) self.spectrum[i] = lag self.spectrum_error[i] = lag_e