Rewrite trigger interpolation to fix numerous bugs and improve unit test coverage
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 togsl_poly_solve_cubic
andgsl_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.
Merge request reports
Activity
added 16 commits
- b3916577 - Use gsl_test for unit tests for helpful messages
- af29febf - Have poly_der fill out new coefficients rather than modifying in place
- 1862b080 - Add code comments
- 6de8b71f - Fix handling of gsl_poly_complex_solve corner cases
- bd90f9b6 - C99 cleanup: move loop variable declarations into for statements
- 865a3b2a - Fix corrupted design matrix in XLALCreateQuadraticFitTriggerInterpolant
- 8408be9a - Fix typo in assignment of output vector
- b761e4ef - Fix conversion of extremum to bounded maximum
- b6f18d9d - Add linear fit for phase
- dfc3d034 - Add nonzero absolute test tolerances
- 950eec81 - Lanczos kernel is zero where |t| > a
- 0bfac5fc - Update Lanczos test outputs
- 18273935 - Update copyright date
- 5f2908ff - Add amplitude-phase interpolant
- 10489dae - Rewrite to avoid heap allocation
- 636c7333 - WIP: better unit test coverage
Toggle commit listadded infra_general lal + 1 deleted label
added 51 commits
-
636c7333...4af4488d - 34 commits from branch
lscsoft:master
- c092aa9a - Use gsl_test for unit tests for helpful messages
- 8ae97013 - Have poly_der fill out new coefficients rather than modifying in place
- dd6550a3 - Add code comments
- 64c2c516 - Fix handling of gsl_poly_complex_solve corner cases
- d1e71b49 - C99 cleanup: move loop variable declarations into for statements
- 21b71177 - Fix corrupted design matrix in XLALCreateQuadraticFitTriggerInterpolant
- d29f81ba - Fix typo in assignment of output vector
- 7673c659 - Fix conversion of extremum to bounded maximum
- 2098e8fb - Add linear fit for phase
- 0c095e29 - Add nonzero absolute test tolerances
- 12a7640b - Lanczos kernel is zero where |t| > a
- 11ca3a93 - Update Lanczos test outputs
- af5cb36d - Update copyright date
- 872aadae - Add amplitude-phase interpolant
- 6821b2c8 - Rewrite to avoid heap allocation
- f92fc842 - WIP: better unit test coverage
- 476482cd - Work around API difference in old versions of GSL
Toggle commit list-
636c7333...4af4488d - 34 commits from branch
added 2 commits
added 1 commit
- 16085be2 - Handle some corner cases for Lanczos interpolant
added 30 commits
-
16085be2...382fec5b - 8 commits from branch
lscsoft:master
- 91e64cda - Use gsl_test for unit tests for helpful messages
- 050e3c04 - Have poly_der fill out new coefficients rather than modifying in place
- 8228e095 - Add code comments
- ddc65827 - Fix handling of gsl_poly_complex_solve corner cases
- 8d72826f - C99 cleanup: move loop variable declarations into for statements
- c1e56c12 - Fix corrupted design matrix in XLALCreateQuadraticFitTriggerInterpolant
- 68ab5b68 - Fix typo in assignment of output vector
- 2c71f1ba - Fix conversion of extremum to bounded maximum
- 94a3a3de - Add linear fit for phase
- 230303a1 - Add nonzero absolute test tolerances
- f26a76c7 - Lanczos kernel is zero where |t| > a
- b0998b3c - Update Lanczos test outputs
- 160dccf7 - Update copyright date
- 66f4725f - Add amplitude-phase interpolant
- b0cf9f51 - Rewrite to avoid heap allocation
- ab300e38 - WIP: better unit test coverage
- 97b1eb6f - Work around API difference in old versions of GSL
- b7a98844 - Skip in SWIG bindings
- 8f0c81fb - WIP: move fast and break things
- 4d3d8e61 - Expand test coverage for CubicSpline
- b6808ce1 - Handle some corner cases for Lanczos interpolant
- 7958b525 - Add Lanczos test
Toggle commit list-
16085be2...382fec5b - 8 commits from branch
added 25 commits
-
37ea4032...8e0116a3 - 2 commits from branch
lscsoft:master
- ef81f0a7 - Use gsl_test for unit tests for helpful messages
- f39d2441 - Have poly_der fill out new coefficients rather than modifying in place
- 5cc343fe - Add code comments
- 5212ff80 - Fix handling of gsl_poly_complex_solve corner cases
- adfb3e53 - C99 cleanup: move loop variable declarations into for statements
- c14c2d5e - Fix corrupted design matrix in XLALCreateQuadraticFitTriggerInterpolant
- 03656362 - Fix typo in assignment of output vector
- 417eae8c - Fix conversion of extremum to bounded maximum
- e13f497a - Add linear fit for phase
- b6982bd6 - Add nonzero absolute test tolerances
- 48119dd2 - Lanczos kernel is zero where |t| > a
- 63033939 - Update Lanczos test outputs
- 85986725 - Update copyright date
- 436e040c - Add amplitude-phase interpolant
- 09124c20 - Rewrite to avoid heap allocation
- be80345d - WIP: better unit test coverage
- 1af32104 - Work around API difference in old versions of GSL
- 54e22db4 - Skip in SWIG bindings
- 92361618 - WIP: move fast and break things
- 7024f1b5 - Expand test coverage for CubicSpline
- 3808b9c4 - Handle some corner cases for Lanczos interpolant
- e5d9eea7 - Add Lanczos test
- e7cb26b5 - Add comment
Toggle commit list-
37ea4032...8e0116a3 - 2 commits from branch
Have you considered wrapping the time series interpolator in TimeSeriesInterp.{c,h} instead? It might not be what you need, but some patches might turn it into what you need, removing the need to maintain separate implementations of similar solutions.
Edited by Kipp Cannon@kipp.cannon, TimeSeriesInterp.h does interpolation of real time series, so it's a very different use case.
added 42 commits
-
e7cb26b5...d4e637fd - 18 commits from branch
lscsoft:master
- ae6ce4cd - Use gsl_test for unit tests for helpful messages
- 7f858762 - Have poly_der fill out new coefficients rather than modifying in place
- 2b2c3863 - Add code comments
- 96185630 - Fix handling of gsl_poly_complex_solve corner cases
- 800043bc - C99 cleanup: move loop variable declarations into for statements
- 8acc8b09 - Fix corrupted design matrix in XLALCreateQuadraticFitTriggerInterpolant
- 80302623 - Fix typo in assignment of output vector
- 083bc0c0 - Fix conversion of extremum to bounded maximum
- 715fd49c - Add linear fit for phase
- 1f0ee4e9 - Add nonzero absolute test tolerances
- 7e420586 - Lanczos kernel is zero where |t| > a
- 9f831865 - Update Lanczos test outputs
- b3e9da54 - Update copyright date
- 24fa67ac - Add amplitude-phase interpolant
- b2b0c0a9 - Rewrite to avoid heap allocation
- 49a48986 - WIP: better unit test coverage
- 52bf96ab - Work around API difference in old versions of GSL
- 106080f4 - Skip in SWIG bindings
- f455a807 - WIP: move fast and break things
- e8fc6656 - Expand test coverage for CubicSpline
- 04a49a27 - Handle some corner cases for Lanczos interpolant
- baf05407 - Add Lanczos test
- 61aeeee2 - Add comment
- 72bda5bc - Remove gsl_finite check
Toggle commit list-
e7cb26b5...d4e637fd - 18 commits from branch
added 144 commits
-
b8694fb6...17abe490 - 121 commits from branch
lscsoft:master
- e48ebb65 - Use gsl_test for unit tests for helpful messages
- 8cc4e793 - Have poly_der fill out new coefficients rather than modifying in place
- ec8b9f6b - Add code comments
- 6d46e7c2 - Fix handling of gsl_poly_complex_solve corner cases
- 5334c315 - C99 cleanup: move loop variable declarations into for statements
- 95f05ecd - Fix corrupted design matrix in XLALCreateQuadraticFitTriggerInterpolant
- 30eef0c3 - Fix typo in assignment of output vector
- 0f91c286 - Fix conversion of extremum to bounded maximum
- 15449736 - Add linear fit for phase
- 9ae02920 - Add nonzero absolute test tolerances
- 17997652 - Lanczos kernel is zero where |t| > a
- 97ca93b5 - Update Lanczos test outputs
- 9fcc69f8 - Update copyright date
- f36bb586 - Add amplitude-phase interpolant
- e66f9d58 - Rewrite to avoid heap allocation
- f958dbba - Improve unit test coverage
- 5d083f01 - Work around API difference in old versions of GSL
- f24521f6 - Skip in SWIG bindings
- b7815130 - Expand test coverage for CubicSpline
- 4f4fd7a3 - Handle some corner cases for Lanczos interpolant
- 57bec625 - Add Lanczos test
- 5a7211de - Add comment
- 9c9fc7e7 - Remove gsl_finite check
Toggle commit list-
b8694fb6...17abe490 - 121 commits from branch
@leo-singer, can you run a new pipeline for this branch?
added 66 commits
-
9c9fc7e7...307aacc6 - 43 commits from branch
lscsoft:master
- 0274a5c2 - Use gsl_test for unit tests for helpful messages
- 180dc7df - Have poly_der fill out new coefficients rather than modifying in place
- 5af2c1ec - Add code comments
- e6a3a97e - Fix handling of gsl_poly_complex_solve corner cases
- ea520d44 - C99 cleanup: move loop variable declarations into for statements
- 2a5015ab - Fix corrupted design matrix in XLALCreateQuadraticFitTriggerInterpolant
- 3211095c - Fix typo in assignment of output vector
- 8ae902c5 - Fix conversion of extremum to bounded maximum
- 842e23b1 - Add linear fit for phase
- 34cfd315 - Add nonzero absolute test tolerances
- 7b1fda03 - Lanczos kernel is zero where |t| > a
- 65d0e2c1 - Update Lanczos test outputs
- 4a0002f1 - Update copyright date
- 6f8e2fa7 - Add amplitude-phase interpolant
- 9a7104bf - Rewrite to avoid heap allocation
- f41ed91f - Improve unit test coverage
- 35f4244e - Work around API difference in old versions of GSL
- 537d2625 - Skip in SWIG bindings
- 99a8437a - Expand test coverage for CubicSpline
- 4b3f25bd - Handle some corner cases for Lanczos interpolant
- 7a912f7a - Add Lanczos test
- 1893070a - Add comment
- beffa7d6 - Remove gsl_finite check
Toggle commit list-
9c9fc7e7...307aacc6 - 43 commits from branch