Skip to content

ENH: Native Simpson, deprecate vnl integral classes, retire netlib quadrature#6470

Open
hjmjohnson wants to merge 3 commits into
InsightSoftwareConsortium:mainfrom
hjmjohnson:deprecate-vnl-simpson-integral
Open

ENH: Native Simpson, deprecate vnl integral classes, retire netlib quadrature#6470
hjmjohnson wants to merge 3 commits into
InsightSoftwareConsortium:mainfrom
hjmjohnson:deprecate-vnl-simpson-integral

Conversation

@hjmjohnson

Copy link
Copy Markdown
Member

Reimplements vnl_simpson_integral with a native composite-Simpson loop, deprecates it and vnl_adaptsimpson_integral behind ITK's three-state legacy guard, and retires the now-unused v3p/netlib/mathews/* sources. Advances the VNL→native numerics direction (#6403, #6230); follows the deprecation pattern of #6441 and #6454.

What changed (3 commits)
Commit What
COMP: Reimplement vnl_simpson_integral Computes the composite Simpson rule (Mathews Alg. 7.2) directly in C++; drops v3p_netlib_simpru_. Loop bounds/accumulation order match the Fortran exactly → bit-identical results.
COMP: Retire unused netlib/mathews Removes simpson.*, trapezod.*, adaquad.f, the standalone netlib_integral_test, and the mangle/unmangle/prototype/CMake/README entries. trapru_/adaquad had no caller; simpru_ was the only live one.
ENH: Deprecate both integral classes __has_include(<itkConfigure.h>)-gated three-state guard (silent default / ITK_LEGACY_REMOVE warning / ITK_FUTURE_LEGACY_REMOVE error), inert for the upstream VXL build. test_integral.cxx defines ITK_LEGACY_TEST.

vnl_adaptsimpson_integral was already pure C++ (it never used the Fortran).

Verification
  • ctest -R test_integral passes against the native engine (analytic checks to 1e-6).
  • All three guard states verified (default silent / ITK_LEGACY_REMOVE warns / ITK_FUTURE_LEGACY_REMOVE errors); +ITK_LEGACY_TEST opt-out keeps the in-tree test silent.
  • pre-commit run --all-files clean.
  • Phase-5 gate: native ≤ netlib in both accuracy and time across 5 integrands (bit-identical error; faster on the oscillatory case). Full report posted as a comment below.

Compute the composite Simpson rule (Mathews Algorithm 7.2) directly in
C++ instead of delegating to the f2c-translated v3p_netlib_simpru_. The
loop bounds and accumulation order match the Fortran exactly, so results
are bit-identical; this removes the only live dependency on the vendored
netlib/mathews sources.
The Mathews netlib routines (simpru_/xsimpru_, trapru_/xtrapru_, the
orphaned adaquad.f) have no remaining caller: vnl_simpson_integral now
computes Simpson natively, vnl_adaptsimpson_integral was always pure
C++, and the trapezoid/adaptive-quadrature symbols were never wired into
any vnl class. Remove the sources, their mangle/unmangle/prototype
entries, the standalone netlib_integral_test, and the mathews entry in
the f2c regeneration recipe.
Gate both headers with the standard ITK three-state legacy guard
(silent default, ITK_LEGACY_REMOVE warning, ITK_FUTURE_LEGACY_REMOVE
error), inert unless itkConfigure.h is present so the upstream VXL build
is unaffected. test_integral defines ITK_LEGACY_TEST so it keeps
exercising the classes under the ITK_LEGACY_REMOVE build.
@hjmjohnson

Copy link
Copy Markdown
Member Author

Phase 5 performance proof — native Simpson vs netlib simpru_

Verdict: PASS the gate. The native composite-Simpson loop is an exact
arithmetic transcription of the f2c simpru_ (identical loop bounds, identical
accumulation order, identical final expression), so its error is bit-identical
to netlib in every measured case (native ≤ netlib accuracy holds with equality),
and its median time is ≤ netlib in every case — notably faster on the
oscillatory integrand. Backing the kept vnl_simpson_integral API with the
native engine and retiring the netlib/mathews sources is justified.

Two-axis benchmark (relative error + median ms, lower = better)
case err_native err_netlib ms_native ms_netlib Eigen/native ≤
x/(1+x²) [0,1], n=50 1.202e-09 1.202e-09 0.0001 0.0001 both
x/(1+x²) [0,1], n=5000 9.610e-16 9.610e-16 0.0084 0.0084 both
gaussian [-6,6], n=500 0.000e+00 0.000e+00 0.0025 0.0025 both
gaussian [-6,6], n=50000 3.012e-15 3.012e-15 0.2520 0.2519 both
sin(50x)·e⁻ˣ [0,3], n=20000 1.095e-12 1.095e-12 0.2312 0.3816 both

Error = |approx − exact| / |exact| against analytic references
(½ln2; √(2π)·erf; closed form for the oscillatory case).

Environment & reproducer
  • OS: Darwin 25.5.0 arm64 — Apple M3 Max — Apple clang 21.0.0
  • Single-threaded; median of 15 reps per case; -O2 -std=c++17
  • Both engines exercised in one process on identical integrands (fairest
    same-compiler/same-CPU/same-integrand-cost comparison); netlib arm is a
    verbatim copy of v3p/netlib/mathews/simpson.c:simpru_.
  • Reproducer: /tmp/simpson_bench/bench.cxx, built with
    c++ -O2 -std=c++17 bench.cxx -o bench, run with
    ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1 ./bench.

@github-actions github-actions Bot added type:Infrastructure Infrastructure/ecosystem related changes, such as CMake or buildbots type:Enhancement Improvement of existing methods or implementation area:ThirdParty Issues affecting the ThirdParty module labels Jun 18, 2026
@hjmjohnson hjmjohnson marked this pull request as ready for review June 18, 2026 02:58
@greptile-apps

This comment was marked as resolved.

@hjmjohnson

Copy link
Copy Markdown
Member Author

@dzenanz TLDR; These functions are not used in ITK or any of the downstream BRAINSTools/ANTs/Remote Modules/Slicer/SlicerExtentions code bases. It is deprecated immediately, and only here in the very unlikely case of a downstream user who may have needed it. NOTE: The re-implementation was a trivial rewrite in C++.

No ITK-dependent tool that uses these functions. 2 of the 3 functions were provably unreachable or not compiled.

@hjmjohnson hjmjohnson requested a review from dzenanz June 18, 2026 12:19
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

area:ThirdParty Issues affecting the ThirdParty module type:Enhancement Improvement of existing methods or implementation type:Infrastructure Infrastructure/ecosystem related changes, such as CMake or buildbots

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant