feat: Bispectrum update and expansion

[bjricketts]

Relevant Issue(s)/PR(s)

This is related to and supersedes #640. This will close #640 upon merge.

Provide an overview of the implemented solution or the fix and elaborate on the modifications.

The current bispectrum implementation is based on the 3rd order cumulant, but astronomical applications of the bispectrum generally use the Fourier method to calculate the bispectrum (Maccarone 2013, Nathan et al 2022). This PR replaces the old methodology with that most commonly used in the field as well as expands features such as:

• Aligning the new code with the standards of more updated modules such as PowerSpectrum
• The cross-bispectrum
• The dynamical bispectrum and associated plotting functions
• The dynamical cross-bispectrum
• Jellyfish plots for the auto-bispectrum
• Normalisation of the bicoherence with all 3 bicoherence normalizations: Kim-Powers (default), Hagihira, Sigl and Chamoun.
• Poisson noise bias subtraction: Wirnitzer (1985)
• Biphase error
• Calculation from all relevant Stingray data objects
• Memory options for dealing with the large data cubes that result from bispectral calculations.
• Tests for all of the above: these are exclusively additions as there was no tests for bispectrum.py previously

Note: As this is quite a large PR, I'm happy to split this into two PRs, one overhauling the Bispectrum calculation, and the other adding the cross-bispectrum and dynamical counterparts.

Fourier-decomposition bi-spectrum

The new method of calculating the bi-spectrum follows Maccarone (2013):

The light curve is cut into m non-overlapping segments. Each segment i is Fourier-transformed to give complex amplitudes Xᵢ(f). The bispectrum is the segment-averaged triple product

B(f₁, f₂) = ⟨ Xᵢ(f₁) · Xᵢ(f₂) · Xᵢ*(f₁ + f₂) ⟩_i

Each factor is a complex Fourier amplitude with its own phase, so the product carries a total phase φ(f₁) + φ(f₂) − φ(f₁+f₂). The biphase is the argument of B:

β(f₁, f₂) = arg B(f₁, f₂)

Practically speaking, we sample a 2-D grid of Fourier frequencies with _bispectrum_frequency_grid(n_bin, dt) and obtain:

• idx1, idx2 which are (nf, nf) bin indices of f_1 and f_2, from np.meshgrid(kbins, kbins, indexing="ij"), so matrix[i, j] is indexed (f₁ = freq[i], f₂ = freq[j]).
• idx3 — the bin index of f₁ + f₂, simply idx1 + idx2.

The code only keeps the positive Nyquist triangle for the auto-bi-spectrum (note that this is not the case for the cross-bi-spectrum).

The actual estimation of the bi-spectrum is done via get_flux_iterable_from_segments (like that of other Fourier objects) and then FFTs of each segment:

ft = fft(flux)                       # one FFT per segment
x1 = ft[idx1]; x2 = ft[idx2]; x3 = ft[idx3]
triple = x1 * x2 * np.conj(x3)       # Tᵢ  — the per-segment bispectrum
denom1 = |x1 * x2|²                  # for the Kim & Powers denominator
denom2 = |x3|²

Rather than storing each 2-D triple (which would cost m × nf × nf complex numbers), the loop folds every segment straight into a fixed set of running sums:

Accumulator	Quantity	Used for
bispec_sum	Σ Tᵢ (complex)	the bispectrum B = ΣTᵢ / m and biphase
denom1_sum	Σ |Xᵢ(f₁)Xᵢ(f₂)|²	bicoherence denominator
denom2_sum	Σ |Xᵢ(f₁+f₂)|²	bicoherence denominator
abs_triple_sum	Σ |Tᵢ|	Hagihira normalization
sum_sq_re, sum_sq_im	Σ (Re Tᵢ)², Σ (Im Tᵢ)²	error on B
cos_sum, sin_sum	Σ Tᵢ/|Tᵢ| (unit phasors)	circular error on the biphase

Because the raw sums are kept, the bicoherence can be re-derived under any normalization after the fact with no re-FFT (recompute_bicoherence), and time-resolved variants can re-combine the sums across time bins. All the sums are additive over segments, so the aberaging is cheap to compute incrementally. To be honest, you could do this the over way as well where you recompute the FFTs, but I think its mostly a matter of taste.

After the pass, with m segments accumulated:

Bispectrum and biphase — directly Maccarone's definitions:

bispec  = bispec_sum / m
biphase = np.angle(bispec)

Bicoherence — calculated in bicoherence_from_sums, which implements the three different normalization conventions from the same sums. The Kim & Powers form Maccarone (2013) uses is |ΣTᵢ|² / (denom1_sum · denom2_sum); sigl_chamoun is its square root; hagihira is |ΣTᵢ| / Σ|Tᵢ|. Each is clipped to [0, 1]. The module defaults to the Kim & Powers normalizations.

Errors — a biphase is only meaningful with a spread, so two are produced (both zero for m = 1 by construction):

• bispec_err: the standard error of the mean of the complex triple, from the variance Σ(Re T)²/m − (Re B)² (and likewise imaginary), divided by m.
• biphase_err: a circular standard error (Fisher 1993). Averaging angles naively is wrong, so the loop sums unit phasors Tᵢ/|Tᵢ|; their mean resultant length r̄ = √(cos_sum² + sin_sum²)/m measures phase concentration, and the error is √(−2 ln r̄)/√m. Tightly coupled frequencies (r̄ → 1) get a small error; scattered phases (r̄ → 0) a large one.
  Finally the redundant region is set to NaN in every 2-D array, and the raw sums plus diagnostics (m, n, dt, df, nphots, segment_size) are packed into the result table's metadata for the class layer to copy onto the object.

Poisson-noise bias subtraction

For photon-counting data the triple product is biased: Poisson noise is white and correlated with itself across the three factors, adding a real offset that Maccarone and later Nathan et al. (2022) treat with the Wirnitzer (1985) correction. When poisson_subtract is set, each per-segment triple is de-biased before accumulation:

triple = triple - ( |x1|² + |x2|² + |x3|² - 2·Nᵢ )

where Nᵢ is the segment's photon count. Subtracting per segment (not once at the end) keeps the correction exact under the segment averaging.

Jellyfish plots

plot_jellyfish takes the per-segment diagonal triples (bispec_diagonal from save_diagonal, or the diagonal of bispec_all), forms a cumulative sum along segments, and normalizes each path by √(denom1·denom2) on the diagonal so the endpoint radius is the Sigl–Chamoun bi-coherence and its angle is the bi-phase. It draws each diagonal frequency as a path in the complex plane, with optional highlighting of the path nearest f0 (fundamental) and f0/2 (sub-harmonic), plus dashed reference circles at the requested bi-coherence levels.

General plots

As users will frequently wish to plot the biphase, bispectrum magnitude and bicoherence, I have also implemented plot functions for each of these, which return the generated axes for free-hand use (for instance if one wants to make multiple panels).

Cross-bi-spectrum

This is the generic form of Bispectrum and allows for the comparison across multiple data channels. The only real difference here is that the frequency grid must also take into account the negative frequencies as the negative frequencies no longer obey f_1<->f_2 that allows for ignoring of the negative frequencies. It is permitted to input 3 different lightcurves, but if only two are provided, it defaults to the first inputted lightcurve as the 3rd curve (and of course reduces to the auto-bi-spectrum case when only one data object is provided).

Dynamical bi-spectrum

I have never seen this used in the literature yet, but I thought I might as well do it while I was here. For the actual calculation of the bi-spectrum, see above. This is implemented as a special case of DynamicalCrossBispectrum which inherits from AveragedCrossBispectrum which inherits from CrossBispectrum. The main inclusions here are the actual calculation of the dynamical bi-spectrum itself and some associated plots for how I think it might be useful (as we're in somewhat "here be dragons" space).

Unlike the dynamical power spectrum, we can't simply take a single bin and compute the bi-spectrum as the bi-coherence is not meaningful without averaging. Therefore, I have implemented such that there is both bin sizes as well as segment sizes (where bin size > segment size). The bicoherence is not meaningful without averaging over about 10 segments, so the code warns the user of this. Otherwise, this is reusing all of the code to calculate the bispectrum. There are a number of helper methods that are essentially redos of DynamicalPowerSpectrum.

There are however 5 plotting options:

• Diagonal: plots the diagonal dynamical bicoherence b(nu, nu, t) which is closest to a dynamical power spectrum.
• Slice: plots the bicoherence at fixed f1 as a function of (f2, time).
• Frame: plots the full (f1, f2) bicoherence map at the time bin nearest t and is best used to plot just one time bin.
• Montage: plots the full bicoherence map at multiple times, inputted by the user.
• Trace: plots a 1-D representation of the biphase and bicoherence with respect to t at fixed coupled frequencies f1, f2.

Finally, there is the shift_and_add method which allows one to track relative frequencies over an observation that would otherwise be smeared out: for instance, a QPO frequency that wanders over time.

Memory options

There are two memory options when calculating the bispectrum (and other derived types): save_all and save_diagonal. save_all retains the whole data cube, while save_diagonal only saves the f1=f2 diagonal, which is useful to measure self-coupling or for producing jellyfish plots.

Any other comments?

Full disclosure: the code in this PR was written with considerable help from Claude. That being said, I've reviewed the produced code and design process carefully, and I'm very much willing to be responsible for the code its outputted. I've kept the design style to as close as other modern modules as possible and reused already existing stingray functions as much as possible.

See the new and accompanying notebook showing all the new features off at StingraySoftware/notebooks#130.

[codecov]

Codecov Report

❌ Patch coverage is 96.12403% with 10 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.05%. Comparing base (840569c) to head (3c94cb5).

Files with missing lines	Patch %	Lines
stingray/fourier.py	96.12%	10 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #982      +/-   ##
==========================================
- Coverage   96.23%   94.05%   -2.18%     
==========================================
  Files          48       48              
  Lines       10099    11039     +940     
==========================================
+ Hits         9719    10383     +664     
- Misses        380      656     +276     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[matteobachetti]

@bjricketts thanks a ton for this! We've been waiting for good bispectrum code for long. For now, I would like to keep both the new and the old codes, so that we can make a through comparison of the two and maybe take the time to deprecate the old one for a while if it just performs worse (or just keep it as an alternative algorithm if it still gives sensible results on real data). The same goes for the notebooks.

[bjricketts]

@matteobachetti I've added back the cumulant method to the code (as well as adding some more test coverage of both the new code and the old cumulant method). It's only implemented for the auto-bispectrum for now as I don't know the 3rd order cumulant method well enough to expand it to the cross-bispectrum. It might be easy, but I haven't read up on it enough to know.

The speed of the methods are essentially dependent on frequency resolution (see the figure below).

The two methods agree with each other to some stupidly small number (1e-16) so I think both methods are perfectly valid to me.

I'll add back the showcase of the cumulant method to the notebook. I'll also add a showcase of using the bispectrum on actual data (GRS1915 RXTE data). I assume I can package the data into the notebook module?

[bjricketts]

I have also updated the notebook to include an explanation/demonstration of the cumulant method as well as demonstrate it on actual data (see the same notebooks PR).

[matteobachetti]

@bjricketts thanks! For real data, I suggest to create a zendo entry or something equivalent outside the repo, that can be downloaded. This helps keeping the size of the repo clean

[bjricketts]

@matteobachetti Okay, that is all done! Let me know if you'd like other changes/additions.

[bjricketts]

Those final two commits should mean that all the workflows should pass minus the doc-links which I haven't changed and are from other notebooks unrelated to this PR.