Interpolation: add derivative method to SpanInterpolators by gda98re · Pull Request #7294 · sxs-collaboration/spectre

gda98re · 2026-06-09T21:25:09Z

Adds a derivative(source_points, values, target_point) method to the intrp::SpanInterpolator interface and its derived classes, returning the derivative of the interpolant at the requested point (the analogue of the existing interpolate). This is needed by downstream code that must time-differentiate an interpolated worldtube quantity.

SpanInterpolator declares the pure-virtual real overload and a virtual complex overload; the base implementation of the complex overload splits the values into real/imaginary parts and calls the real overload. Derived classes override the complex overload when a specialized version is more efficient.
LinearSpanInterpolator returns the constant slope between the two points.
CubicSpanInterpolator adds derivative_impl, the analytic derivative of the cubic Lagrange interpolant.
BarycentricRationalSpanInterpolator uses boost's barycentric_rational::prime.

Tested with an approximate-fidelity check against the analytic derivative of A cos(w t), and -- since a degree-N interpolant reproduces a degree-N polynomial exactly -- with exact-polynomial derivative checks for the linear, cubic, and barycentric interpolators.

Co-Authored-By: Claude Opus 4.8 noreply@anthropic.com

Proposed changes

Code review checklist

The code is documented and the documentation renders correctly. Run
make doc to generate the documentation locally into BUILD_DIR/docs/html.
Then open index.html.
The code follows the stylistic and code quality guidelines listed in the
code review guide.
The PR lists upgrade instructions and is labeled bugfix or
new feature if appropriate.
If a coding agent is used, have one of
"Co-Authored-By: Claude Sonnet 4.6 noreply@anthropic.com",
"Co-Authored-by: Codex noreply@openai.com", or
"Co-Authored-By: GitHub Copilot CLI noreply@microsoft.com"
as the last line of the commit, depending on the agent.

Further comments

Part of the CauchySecondOrder initial-data stack.

wthrowe · 2026-06-10T18:28:05Z

+  /// Evaluate the derivative of the interpolant of the function represented by
+  /// complex `values` at `source_points`, evaluated at the requested
+  /// `target_point`.
+  std::complex<double> derivative(


I see that this is just copying the implementation of interpolate, but don't you really want this to be virtual? This generic implementation is extremely inefficient compared to the specialized implementations in the derived classes.

Related: do we want an (additional?) function that returns the value and derivative? I could see it being cheaper to compute both simultaneously than back-to-back.

For the other PRs for the CCE ID it won't be needed as I need Du(DyJ) but not DyJ itself. But if you want I can add it

I made the complex derivative overload virtual in the base class and marked the Linear/Cubic complex overloads override. The generic split into real/imaginary implementation is now just the fallback for derived classes (like BarycentricRational) that don't specialize the complex path.

wthrowe · 2026-06-10T18:30:44Z

+                              linear_interpolator_values.size()},
+      target_point);
+  CHECK(real_linear_derivative ==
+        approx(linear_interpolator_values[1] - linear_interpolator_values[0]));


Test with a denominator that isn't 1.

test_linear_interpolator now uses points {0.5, 2.0} (spacing 1.5, nonzero offset), and the expected value/derivative are computed with the actual spacing and offset, so the interpolation denominator is not 1.

wthrowe · 2026-06-10T18:37:13Z

+    interpolator_values[i] =
+        amplitude * cos(frequency * interpolator_points[i]);
+  }
+  const double target_time = value_dist(*gen) / 100.0;


I'm confused about why this is working. The control points generated above range from 0 to a bit less than 0.01, and this generates a target point between 0.001 and 0.01, so this could generate a target outside the data range. Even beyond that, don't you have to obey required_number_of_points_before_and_after, which should restrict the valid target points to a small interval in the center?

The points are now a uniform grid point_spacing * i, and the target is drawn inside the central interval [points[n-1], points[n]) with n = required_number_of_points_before_and_after(), so the required number of points always lies on each side of the target.

wthrowe · 2026-06-10T18:46:10Z

+    return result;
+  };
+  // Unit-spaced points keep the Lagrange/barycentric formulas well conditioned
+  // so the exact result is recovered to roundoff.


This would be a much stronger test if you didn't do this. You want to make sure no points are close together, but making all the spacings exactly 1 is not necessary. Perturbing each point by a random number of magnitude below 0.1 should work fine.

Ok, the points are now i + delta_i with delta_i drawn uniformly from [-0.1, 0.1] and the points stay well separated

wthrowe · 2026-06-10T18:47:24Z

+  const auto interpolator_derivative_result = interpolator.derivative(
+      gsl::span<const double>{points.data(), points.size()},
+      gsl::span<const ValueType>{values.data(), values.size()}, target_point);
+  const Approx exact_approx = Approx::custom().epsilon(1.0e-10).scale(1.0);


You claim repeatedly above that this should be exact to roundoff.

here you mean to tighten up the epsilon to true roundoff, or rewording?

Depends on if it actually is accurate to roundoff. If it is, using a tighter tolerance would be preferable, but if it's not that won't work and the comments should be changed.

Tightened the tolerance from 1e-10 to 1e-13 and reworded the comment to say the derivative matches the analytic value "to roundoff". Confirmed the exactness test passes at 1e-13 for the linear, cubic, and barycentric interpolators.

wthrowe · 2026-06-10T18:59:57Z

+    // The linear derivative is a midpoint finite-difference approximation;
+    // its error grows by one factor of h relative to the interpolated value.
+    const Approx derivative_approx =
+        Approx::custom().epsilon(1.0e-1).scale(1.0);


You should be able to do better than this. The linear approximation should be able to get the value to better than 1e-4 and the derivative to 1e-2. Evan a zeroth-order approximation should get the value to 1e-4 for the chosen test function.

Tightened the linear tolerances to 1e-4 (value) and 1e-2 (derivative).

wthrowe · 2026-06-10T19:04:34Z

+    test_interpolator_derivative_approximate_fidelity<DataVector>(
+        make_not_null(&gen), LinearSpanInterpolator{}, derivative_approx);
+    test_interpolator_derivative_approximate_fidelity<ComplexDataVector>(
+        make_not_null(&gen), LinearSpanInterpolator{}, derivative_approx);


Call test_interpolator_derivative_is_exact?

wthrowe · 2026-06-10T20:08:34Z

@@ -41,5 +41,23 @@ double BarycentricRationalSpanInterpolator::interpolate(
  return interpolant(target_point);


[Arbitrary line for threading purposes] Do you have some reason to believe there might be out-of-tree implementations of SpanInterpolator? If not, I really don't think we need the details of the interface changes in the release notes' upgrade instructions.

No, I've dropped the interface change details from the upgrade instructions block.

wthrowe · 2026-06-10T20:11:46Z

+  return derivative_impl(source_points, values, target_point);
+}
+
+// NOLINTNEXTLINE(readability-convert-member-functions-to-static)


If you don't end up making this a virtual function overload, why not do the clang-tidy thing instead of NOLINTing?

I removed the NOLINT

Adds a `derivative(source_points, values, target_point)` method to the `intrp::SpanInterpolator` interface and its derived classes, returning the derivative of the interpolant at the requested point (the analogue of the existing `interpolate`). This is needed by downstream code that must time-differentiate an interpolated worldtube quantity. - `SpanInterpolator` declares the pure-virtual real overload and a virtual complex overload whose base implementation splits the values into real/imaginary parts; derived classes override it when a specialized complex implementation is more efficient. - `LinearSpanInterpolator` returns the constant slope between the two points. - `CubicSpanInterpolator` adds `derivative_impl`, the analytic derivative of the cubic Lagrange interpolant. - `BarycentricRationalSpanInterpolator` uses boost's `barycentric_rational::prime`. Tested with an approximate-fidelity check against the analytic derivative of `A cos(w t)`, and -- since a degree-N interpolant reproduces a degree-N polynomial exactly -- with exact-polynomial derivative checks for the linear, cubic, and barycentric interpolators. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

wthrowe

Looks good.

In the future, please push fixup commits if you have to make significant changes after review.

github-actions · 2026-06-16T23:56:57Z

@gda98re looks like this is your first contribution to SpECTRE. Welcome! 🎉 Your contribution is much appreciated, and we invite you to add your name to the author list:

Edit the Metadata.yaml file in the repository.
Add an entry to the list in Authors.Contributors with your name, affiliation, etc. Note that the list is ordered alphabetically by last name.
Commit the change on a new branch and open a pull request with the change.

Once the pull request is merged, your name will appear on the SpECTRE DOI on Zenodo with the next public release.

gda98re force-pushed the pr-cse-2-interpolator-derivative branch from fbe5323 to f72d3c4 Compare June 10, 2026 03:24

wthrowe requested changes Jun 10, 2026

View reviewed changes

gda98re force-pushed the pr-cse-2-interpolator-derivative branch 3 times, most recently from 667aac8 to ab34d2c Compare June 16, 2026 01:08

gda98re force-pushed the pr-cse-2-interpolator-derivative branch from ab34d2c to 38246fb Compare June 16, 2026 01:34

gda98re requested a review from wthrowe June 16, 2026 16:34

wthrowe mentioned this pull request Jun 16, 2026

Failure in Unit.GrMhd.ValenciaDivClean.Sources #7298

Open

wthrowe approved these changes Jun 16, 2026

View reviewed changes

wthrowe merged commit 673e9c1 into sxs-collaboration:develop Jun 16, 2026
47 of 48 checks passed

		@@ -41,5 +41,23 @@ double BarycentricRationalSpanInterpolator::interpolate(
		return interpolant(target_point);

Conversation

gda98re commented Jun 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed changes

Code review checklist

Further comments

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gda98re Jun 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gda98re Jun 16, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gda98re Jun 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

wthrowe left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

github-actions Bot commented Jun 16, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

gda98re commented Jun 9, 2026 •

edited

Loading

gda98re Jun 10, 2026 •

edited

Loading

gda98re Jun 16, 2026 •

edited

Loading

gda98re Jun 10, 2026 •

edited

Loading