Interpolation: add derivative method to SpanInterpolators

src/NumericalAlgorithms/Interpolation/BarycentricRationalSpanInterpolator.cpp

@@ -41,5 +41,23 @@ double BarycentricRationalSpanInterpolator::interpolate(
   return interpolant(target_point);
 }
 
+double BarycentricRationalSpanInterpolator::derivative(
+    const gsl::span<const double>& source_points,
+    const gsl::span<const double>& values, const double target_point) const {
+  if (UNLIKELY(source_points.size() < min_order_ + 1)) {
+    ERROR("provided independent values for interpolation too small.");
+  }
+  // Boost moved barycentric_rational from boost::math to
+  // boost::math::interpolators in version 1.77.
+  // NOLINTNEXTLINE(google-build-using-namespace)
+  using namespace boost::math;
+  // NOLINTNEXTLINE(google-build-using-namespace)
+  using namespace boost::math::interpolators;
+  const barycentric_rational<double> interpolant(
+      source_points.data(), values.data(), source_points.size(),
+      std::min(source_points.size() - 1, max_order_));
+  return interpolant.prime(target_point);
+}
+
 PUP::able::PUP_ID intrp::BarycentricRationalSpanInterpolator::my_PUP_ID = 0;
 }  // namespace intrp

src/NumericalAlgorithms/Interpolation/BarycentricRationalSpanInterpolator.hpp

@@ -69,13 +69,18 @@ class BarycentricRationalSpanInterpolator : public SpanInterpolator {
     return std::make_unique<BarycentricRationalSpanInterpolator>(*this);
   }
 
-  // to get the generic overload:
+  // to get the generic overloads:
+  using SpanInterpolator::derivative;
   using SpanInterpolator::interpolate;
 
   double interpolate(const gsl::span<const double>& source_points,
                      const gsl::span<const double>& values,
                      double target_point) const override;
 
+  double derivative(const gsl::span<const double>& source_points,
+                    const gsl::span<const double>& values,
+                    double target_point) const override;
+
   size_t required_number_of_points_before_and_after() const override {
     return min_order_ / 2 + 1;
   }

src/NumericalAlgorithms/Interpolation/CubicSpanInterpolator.cpp

@@ -35,6 +35,44 @@ SPECTRE_ALWAYS_INLINE ValueType interpolate_impl(
          ((t0 - t1) * (t0 - t2) * (t1 - t2) * (t0 - t3) * (t1 - t3) *
           (t2 - t3));
 }
+
+template <typename ValueType>
+SPECTRE_ALWAYS_INLINE ValueType derivative_impl(
+    const gsl::span<const double>& source_points,
+    const gsl::span<const ValueType>& values, const double target_point) {
+  const double t0 = source_points[0];
+  const double t1 = source_points[1];
+  const double t2 = source_points[2];
+  const double t3 = source_points[3];
+
+  const auto d0 = values[0];
+  const auto d1 = values[1];
+  const auto d2 = values[2];
+  const auto d3 = values[3];
+
+  const double x_minus_t0 = target_point - t0;
+  const double x_minus_t1 = target_point - t1;
+  const double x_minus_t2 = target_point - t2;
+  const double x_minus_t3 = target_point - t3;
+
+  // Sum-of-products derivative of each Lagrange basis numerator:
+  //   d/dx [(x - a)(x - b)(x - c)] = (x-b)(x-c) + (x-a)(x-c) + (x-a)(x-b)
+  const double s0 = (x_minus_t2 * x_minus_t3) + (x_minus_t1 * x_minus_t3) +
+                    (x_minus_t1 * x_minus_t2);
+  const double s1 = (x_minus_t2 * x_minus_t3) + (x_minus_t0 * x_minus_t3) +
+                    (x_minus_t0 * x_minus_t2);
+  const double s2 = (x_minus_t1 * x_minus_t3) + (x_minus_t0 * x_minus_t3) +
+                    (x_minus_t0 * x_minus_t1);
+  const double s3 = (x_minus_t1 * x_minus_t2) + (x_minus_t0 * x_minus_t2) +
+                    (x_minus_t0 * x_minus_t1);
+
+  return (d0 * s0 * (t1 - t2) * (t1 - t3) * (t2 - t3) -
+          d1 * s1 * (t0 - t2) * (t0 - t3) * (t2 - t3) +
+          d2 * s2 * (t0 - t1) * (t0 - t3) * (t1 - t3) -
+          d3 * s3 * (t0 - t1) * (t0 - t2) * (t1 - t2)) /
+         ((t0 - t1) * (t0 - t2) * (t1 - t2) * (t0 - t3) * (t1 - t3) *
+          (t2 - t3));
+}
 }  // namespace
 
 double CubicSpanInterpolator::interpolate(
@@ -50,5 +88,18 @@ std::complex<double> CubicSpanInterpolator::interpolate(
   return interpolate_impl(source_points, values, target_point);
 }
 
+double CubicSpanInterpolator::derivative(
+    const gsl::span<const double>& source_points,
+    const gsl::span<const double>& values, double target_point) const {
+  return derivative_impl(source_points, values, target_point);
+}
+
+std::complex<double> CubicSpanInterpolator::derivative(
+    const gsl::span<const double>& source_points,
+    const gsl::span<const std::complex<double>>& values,
+    double target_point) const {
+  return derivative_impl(source_points, values, target_point);
+}
+
 PUP::able::PUP_ID intrp::CubicSpanInterpolator::my_PUP_ID = 0;
 }  // namespace intrp

src/NumericalAlgorithms/Interpolation/CubicSpanInterpolator.hpp

@@ -54,6 +54,15 @@ class CubicSpanInterpolator : public SpanInterpolator {
       const gsl::span<const std::complex<double>>& values,
       double target_point) const;
 
+  double derivative(const gsl::span<const double>& source_points,
+                    const gsl::span<const double>& values,
+                    double target_point) const override;
+
+  std::complex<double> derivative(
+      const gsl::span<const double>& source_points,
+      const gsl::span<const std::complex<double>>& values,
+      double target_point) const override;
+
   size_t required_number_of_points_before_and_after() const override {
     return 2;
   }

src/NumericalAlgorithms/Interpolation/LinearSpanInterpolator.cpp

@@ -20,6 +20,13 @@ SPECTRE_ALWAYS_INLINE ValueType interpolate_impl(
                          (source_points[1] - source_points[0]) *
                          (target_point - source_points[0]);
 }
+
+template <typename ValueType>
+SPECTRE_ALWAYS_INLINE ValueType
+derivative_impl(const gsl::span<const double>& source_points,
+                const gsl::span<const ValueType>& values) {
+  return (values[1] - values[0]) / (source_points[1] - source_points[0]);
+}
 }  // namespace
 
 double LinearSpanInterpolator::interpolate(
@@ -35,5 +42,19 @@ std::complex<double> LinearSpanInterpolator::interpolate(
   return interpolate_impl(source_points, values, target_point);
 }
 
+double LinearSpanInterpolator::derivative(
+    const gsl::span<const double>& source_points,
+    const gsl::span<const double>& values,
+    const double /*target_point*/) const {
+  return derivative_impl(source_points, values);
+}
+
+std::complex<double> LinearSpanInterpolator::derivative(
+    const gsl::span<const double>& source_points,
+    const gsl::span<const std::complex<double>>& values,
+    const double /*target_point*/) const {
+  return derivative_impl(source_points, values);
+}
+
 PUP::able::PUP_ID intrp::LinearSpanInterpolator::my_PUP_ID = 0;
 }  // namespace intrp

src/NumericalAlgorithms/Interpolation/LinearSpanInterpolator.hpp

@@ -50,6 +50,15 @@ class LinearSpanInterpolator : public SpanInterpolator {
       const gsl::span<const std::complex<double>>& values,
       double target_point) const;
 
+  double derivative(const gsl::span<const double>& source_points,
+                    const gsl::span<const double>& values,
+                    double target_point) const override;
+
+  std::complex<double> derivative(
+      const gsl::span<const double>& source_points,
+      const gsl::span<const std::complex<double>>& values,
+      double target_point) const override;
+
   size_t required_number_of_points_before_and_after() const override {
     return 1;
   }

src/NumericalAlgorithms/Interpolation/SpanInterpolator.cpp

@@ -31,4 +31,24 @@ std::complex<double> SpanInterpolator::interpolate(
                   gsl::span<const double>(imag_part.data(), imag_part.size()),
                   target_point)};
 }
+
+std::complex<double> SpanInterpolator::derivative(
+    const gsl::span<const double>& source_points,
+    const gsl::span<const std::complex<double>>& values,
+    const double target_point) const {
+  // the operation below to get the real and imag parts does not alter the
+  // contents of the span, so the const-cast is safe.
+  const ComplexDataVector view{
+      const_cast<std::complex<double>*>(values.data()),  // NOLINT
+      values.size()};
+  const DataVector real_part = real(view);
+  const DataVector imag_part = imag(view);
+  return std::complex<double>{
+      derivative(source_points,
+                 gsl::span<const double>(real_part.data(), real_part.size()),
+                 target_point),
+      derivative(source_points,
+                 gsl::span<const double>(imag_part.data(), imag_part.size()),
+                 target_point)};
+}
 }  // namespace intrp

src/NumericalAlgorithms/Interpolation/SpanInterpolator.hpp

@@ -55,6 +55,23 @@ class SpanInterpolator : public PUP::able {
       const gsl::span<const std::complex<double>>& values,
       double target_point) const;
 
+  /// Evaluate the derivative of the interpolant of the function represented by
+  /// `values` at `source_points`, evaluated at the requested `target_point`.
+  virtual double derivative(const gsl::span<const double>& source_points,
+                            const gsl::span<const double>& values,
+                            double target_point) const = 0;
+
+  /// Evaluate the derivative of the interpolant of the function represented by
+  /// complex `values` at `source_points`, evaluated at the requested
+  /// `target_point`. This generic implementation calls the real version on the
+  /// real and imaginary parts separately; derived classes should override it
+  /// when a specialized complex implementation that avoids the split is more
+  /// efficient.
+  virtual std::complex<double> derivative(
+      const gsl::span<const double>& source_points,
+      const gsl::span<const std::complex<double>>& values,
+      double target_point) const;
+
   /// The number of domain points that should be both before and after the
   /// requested target point for best interpolation. For instance, for a linear
   /// interpolator, this function would return `1` to request that the target is