Description

Numerical convergence testing revealed several serious bugs in TriggerInterpolation.{h,c} as well as a fundamental unsuitability of the Lanczos and real/imaginary part cubic spline interpolation schemes. This code is used by gstlal and affects localization results for offline runs but not online runs; gstlal online runs provide SNR time series to BAYESTAR and so the trigger times are ignored.

This is a rewrite to fix those bugs, improve test coverage, eliminate unnecessary heap allocations, and add one more interpolation algorithm.

Bugs fixed:

• Lanczos kernel was not zeroed outside of filter window.
• Handle corner cases of polynomials with leading zero coefficients. The polynomial that is solved in the cubic spline interplant is generally a quintic, but in degenerate cases may be of any order. Unfortunately, gsl_poly_complex_solve does not automatically remove leading zero coefficients. Since we are taking care of these special cases, also dispatch to gsl_poly_solve_cubic and gsl_poly_solve_quadratic when appropriate.
• Fix corruption of design matrix passed to gsl_multifit_linear due to improper casting of unsigned integer values.
• Fix a typo in the quadratic fit implementation that resulted in invalid data being used for the fit.

The added interpolation algorithm does cubic spline interpolation of the amplitude and phase separately (CubicSplineAmpPhase). This method exhibits less ripple in the amplitude than the Lanczos algorithm or cubic spline interpolation of the real and imaginary parts. (See, for example, ligo.skymap!155 (merged).) This new algorithm will likely be what I recommend that search pipelines use.

In addition to bug fixes and the new interpolation algorithm, I did some refactoring to improve code reuse and make better use of C99 syntax.

Also, I noticed that most of the heap allocations were pretty small. My taste has evolved since I first wrote this code, and I decided to move all of the allocations onto the stack using C99 variable-length arrays where necessary. This should make the code a bit more performant. It will also be a bit easier for client, since it will no longer require setup/teardown calls to allocate workspaces.

Along with the elimination of the setup/teardown calls, I am providing a new API. Previously, usage required three calls (example for the CubicSpline algorithm):

/* these symbols defined elsewhere */
unsigned int window;
double tmax;
double complex ymax;
const double complex *y;

CubicSplineTriggerInterpolant *interp = XLALCreateCubicSplineTriggerInterpolant(window);
XLALCOMPLEX16ApplyCubicSplineTriggerInterpolant(interp, &tmax, &ymax, y);
XLALDestroyCubicSplineTriggerInterpolant(interp);

These are still supported, but the new API consists of just one call:

XLALTriggerInterpolateCubicSpline(&tmax, &ymax, y, window);

The old API is still available in TriggerInterpolation.h, but will issue GCC deprecation warnings if used. The new API is available from TriggerInterpolate.h.

API Changes and Justification

Backwards Compatible Changes

• This change introduces no API changes
• This change adds new API calls

Backwards Incompatible Changes

• This change modifies an existing API
• This change removes an existing API

If any of the Backwards Incompatible check boxes are ticked please provide a justification why this change is necessary and why it needs to be done in a backwards incompatible way.

Review Status

This will be part of the BAYESTAR review. CC @zoheyr-doctor, @rosa.poggiani, @surabhi.sachdev.