Linear interpolated FFT-based time-correlated signals

[vhaasteren]

This PR introduces the linear interpolated FFT-based rank reduced time-correlated signals. These replace the diagonal Phi matrix approach that Enterprise has traditionally used based on the Lentati et al. (2013) method.

Using it is pretty simple. Instead of:

rn_pl = powerlaw(log10_A=rn_log10_A, gamma=rn_gamma)
rn_phi = gp_signals.FourierBasisGP(spectrum=rn_pl, components=n_components, Tspan=Tspan)

one would do:

rn_pl = powerlaw(log10_A=rn_log10_A, gamma=rn_gamma)
rn_fft = gp_signals.FFTBasisGP(spectrum=rn_pl, components=n_components, oversample=3, cutoff=3, Tspan=Tspan, start_time=start_time)

so the same spectral function can be used. Free spectrum is NOT supported yet. Instead of components, we can also pass knots=, where it is understood that knots=2*n_components+1. This is because components actually means frequencies. In the time-domain, the number of knots is the number of modes+1, because we cannot just omit the DC term.

The oversample parameter determines how densely the PSD is sampled in frequencies. With oversample=1 we would use frequencies at spacing of df=1/T. With oversample=3 (the default), the frequency spacing is df=1/(3T). Note that this is a way to numerically approximate the Wiener-Khinchin integral. With oversample sufficiently large, the FFT is an excellent approximation of the analytical integral. For powerlaw signals, oversample=3 seems a very reasonable choice.

The cutoff parameter is used to specify below which frequency fcut = 1 / (cutoff*Tspan) we set the PSD equal to zero. Note that this parameterization (which is also in Discovery) is a bit ambiguous, as fcut may not correspond to an actual bin of the FFT: especially if oversample is not a high number this can cause a mismatch. In case of a mismatch, fcut is rounded up to the next oversampled-FFT frequency bin. Instead of cutoff, the parameter cutbins can also be used (this overrides cutoff). With cutbins the low-frequency cutoff is set at: fcut = cutbins / (oversample * Tspan), and its interpretation is less ambiguous: it is the number of bins of the over-sampled PSD of the FFT that is being zeroed out.

For common signals, instead of:

gw_pl = powerlaw(log10_A=gw_log10_A, gamma=gw_gamma)
orf = utils.hd_orf()
crn_phi = gp_signals.FourierBasisCommonGP(gw_pl, orf, components=20, name='gw', Tspan=Tspan)

one would do:

gw_pl = powerlaw(log10_A=gw_log10_A, gamma=gw_gamma)
orf = utils.hd_orf()
crn_fft = gp_signals.FFTBasisCommonGP(gw_pl, orf, components=20, name='gw', Tspan=Tspan, start_time=start_time)

DM-variations and Chromatic noise can be similarly set up:

nknots = 81
dm_basis = utils.create_fft_time_basis_dm(nknots=nknots)
dm_pl = powerlaw(log10_A=dm_log10_A, gamma=dm_gamma)
dm_fft = gp_signals.FFTBasisGP(dm_pl, basis=dm_basis, nknots=nknots, name='dmgp')

chrom_basis = utils.create_fft_time_basis_chromatic(nknots=nknots, idx=chrom_idx)
chrom_pl = powerlaw(log10_A=chrom_log10_A, gamma=chrom_gamma)
chrom_fft = gp_signals.FFTBasisGP(chrom_pl, basis=chrom_basis, nknots=nknots, name='chromgp')

Subtleties

Enterprise allows one to combine basis functions when they are the same. This
is especially useful when analyzing common signals which have the same basis as
a single-pulsar signal, such as one would have with red noise and a correlated
GWB. This can be done with the combine=True option in FFTBasisGP and
FFTBasisCommonGP. Default is combine=True. The subtlety is that modern PTA
datasets typically have large gaps, which causes some of the time-domain basis
functions to basically be all zeros. Therefore, some basis functions that you
would not expect to be identical will be combined.

The above is not a bug. Combining such bases and the corresponding Phi matrix
does not matter, because the basis is zero, and that part of the signal has no
bearing on the data or the model. However, when doing signal reconstruction,
such as with la_forge or utils.ConditionalGP, make sure to set combine=False.

Extra comment

An additional re-design of certain patterns is needed to accelerate TNT construction using sparse matrices for faster single pulsar noise analysis. Because of the severity of changes, that will be a different PR.

[codecov]

Codecov Report

Attention: Patch coverage is 88.75000% with 18 lines in your changes missing coverage. Please review.

Project coverage is 72.53%. Comparing base (e82c676) to head (e04d9d4).
Report is 8 commits behind head on dev.

Files with missing lines	Patch %	Lines
enterprise/signals/gp_signals.py	89.81%	11 Missing ⚠️
enterprise/signals/gp_priors.py	33.33%	4 Missing ⚠️
enterprise/signals/gp_bases.py	92.00%	2 Missing ⚠️
enterprise/signals/utils.py	94.73%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##              dev     #411      +/-   ##
==========================================
+ Coverage   71.53%   72.53%   +1.00%     
==========================================
  Files          13       13              
  Lines        3246     3401     +155     
==========================================
+ Hits         2322     2467     +145     
- Misses        924      934      +10     

Files with missing lines	Coverage Δ
enterprise/signals/signal_base.py	91.50% <100.00%> (+1.25%)	⬆️
enterprise/signals/utils.py	82.18% <94.73%> (+0.07%)	⬆️
enterprise/signals/gp_bases.py	75.72% <92.00%> (+5.21%)	⬆️
enterprise/signals/gp_priors.py	84.28% <33.33%> (-4.78%)	⬇️
enterprise/signals/gp_signals.py	71.14% <89.81%> (+4.11%)	⬆️

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a536b06...e04d9d4. Read the comment docs.

🚀 New features to boost your workflow:

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