NumericalDiff: Implement Ridders' method libeigen/eigen!3385
diff --git a/contrib/Eigen/NumericalDiff b/contrib/Eigen/NumericalDiff index ac98ce3..81fec45 100644 --- a/contrib/Eigen/NumericalDiff +++ b/contrib/Eigen/NumericalDiff
@@ -17,6 +17,7 @@ /** * \defgroup NumericalDiff_Module Numerical differentiation module + * \brief Finite-difference Jacobians and differentiation weights. * * \code * #include <contrib/Eigen/NumericalDiff> @@ -28,19 +29,13 @@ * is a different method and has its own module in Eigen : \ref * AutoDiff_Module. * - * NumericalDiff implements the "Forward" and "Central" schemes. Stencil - * computes finite-difference weights for arbitrary fixed-size grids and - * derivative orders; the caller supplies the points and chooses the - * step size. The "Forward" and "Central" schemes are basic methods, and - * there exist some more elaborated way of computing such approximates. - * They are implemented using both proprietary and free software, and - * usually requires linking to an external library. It is very easy for - * you to write a functor using such software, and the purpose is quite - * orthogonal to what we want to achieve with Eigen. - * - * This is why we will not provide wrappers for every great numerical - * differentiation software that exist, but should rather stick with those - * basic ones, that still are useful for testing. + * This module provides the following classes: + * - \ref NumericalDiff extends a functor with Jacobian computation using + * \c Forward or \c Central finite differences at a fixed step. + * - \ref Ridders stores a functor and computes its Jacobian using adaptive + * central differences and Richardson extrapolation, with optional error estimates. + * - \ref Stencil computes finite-difference weights for arbitrary fixed-size + * grids and derivative orders, using points and a step size chosen by the caller. * * Also, the \ref NonLinearOptimization_Module needs this in order to * provide full features compatibility with the original (c)minpack @@ -56,6 +51,7 @@ // IWYU pragma: begin_exports #include "src/NumericalDiff/NumericalDiff.h" #include "src/NumericalDiff/Stencil.h" +#include "src/NumericalDiff/Ridders.h" // IWYU pragma: end_exports #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
diff --git a/contrib/Eigen/src/NumericalDiff/Ridders.h b/contrib/Eigen/src/NumericalDiff/Ridders.h new file mode 100644 index 0000000..0c9042e --- /dev/null +++ b/contrib/Eigen/src/NumericalDiff/Ridders.h
@@ -0,0 +1,258 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#ifndef EIGEN_NUMERICALDIFF_RIDDERS_H +#define EIGEN_NUMERICALDIFF_RIDDERS_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +namespace internal { + +template <typename F, typename = void> +struct ridders_solver_traits {}; + +template <typename F> +struct ridders_solver_traits<F, void_t<typename F::QRSolver>> { + using QRSolver = typename F::QRSolver; +}; + +} // namespace internal + +/** \class Ridders + * \ingroup NumericalDiff_Module + * \brief Adaptive numerical Jacobians using Ridders' extrapolation method. + * + * Stores a functor by value and computes its Jacobian from central differences at successively + * smaller steps. Richardson extrapolation in the squared step size reduces truncation error; + * the Jacobian column with the smallest estimated error is retained. The stopping test uses + * the maximum absolute componentwise difference over the entire column. + * + * The functor supplies \c Scalar, \c InputType, \c ValueType, \c JacobianType, + * \c InputsAtCompileTime, \c ValuesAtCompileTime, \c values(), and \c operator(), + * as for NumericalDiff. The scalar type must be real and non-integer. + * The functor's \c inputs() method is required only when calling inputs(). + * + * \warning The functor must be smooth and defined throughout the initial step interval. + * + * \code + * Eigen::Ridders<MyFunctor> derivative(functor); + * derivative.df(x, jacobian); // Returns the number of function evaluations. + * derivative.df(x, jacobian, error); // Also returns an absolute error estimate per input coordinate. + * \endcode + * + * \tparam Functor_ The vector-valued functor to differentiate. + */ +template <typename Functor_> +class Ridders : public internal::ridders_solver_traits<Functor_> { + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + + using Functor = Functor_; + using Scalar = typename Functor::Scalar; + using InputType = typename Functor::InputType; + using ValueType = typename Functor::ValueType; + using JacobianType = typename Functor::JacobianType; + + static constexpr int InputsAtCompileTime = Functor::InputsAtCompileTime; + static constexpr int ValuesAtCompileTime = Functor::ValuesAtCompileTime; + + /** \brief Parameters controlling the initial interval and the maximum refinement cost. */ + struct Parameters { + /** \brief Initial step scale, giving \f$h_j = \mathrm{initialStep}\,\max(|x_j|, 1)\f$. + * Must be positive and finite. */ + Scalar initialStep = Scalar(0.01); + /** \brief Step reduction factor applied at each refinement. + * The factor must be greater than one and its square must be finite. + * Values close to one amplify rounding errors in extrapolation and may repeat sampling points. + * Refinement stops when both sampling points repeat; this does not detect all ill-conditioning. */ + Scalar stepReduction = Scalar(2); + /** \brief Maximum number of central differences per input coordinate. + * Must be at least two and at most half the maximum representable Index. */ + int maxSteps = 10; + }; + + Ridders() = default; + + /** \brief Constructs from a copy of a functor and differentiation parameters. + * \param functor The functor to differentiate. + * \param parameters The differentiation parameters. */ + explicit Ridders(const Functor& functor, const Parameters& parameters = Parameters()) : m_functor(functor) { + setParameters(parameters); + } + + /** \brief Sets the differentiation parameters. + * \param parameters The differentiation parameters. */ + void setParameters(const Parameters& parameters) { + eigen_assert(parameters.initialStep > Scalar(0) && (numext::isfinite)(parameters.initialStep)); + eigen_assert(parameters.stepReduction > Scalar(1) && + (numext::isfinite)(parameters.stepReduction * parameters.stepReduction)); + eigen_assert(parameters.maxSteps >= 2); + eigen_assert(parameters.maxSteps <= NumTraits<Index>::highest() / 2); + m_parameters = parameters; + } + + /** \returns the differentiation parameters. */ + const Parameters& parameters() const { return m_parameters; } + /** \returns the stored functor. */ + Functor& functor() { return m_functor; } + /** \returns the stored functor. */ + const Functor& functor() const { return m_functor; } + + /** \brief Evaluates the stored functor. + * \param[in] x The evaluation point. + * \param[out] values The function values. + * \returns the value returned by the functor. */ + int operator()(const InputType& x, ValueType& values) const { return m_functor(x, values); } + /** \returns the stored functor's number of inputs. */ + int inputs() const { return m_functor.inputs(); } + /** \returns the stored functor's number of values. */ + int values() const { return m_functor.values(); } + + /** + * \brief Computes the Jacobian using adaptive central differences. + * + * The caller sizes \a jac to values() by x.size(). At most \c 2*maxSteps evaluations + * are made per coordinate. Difference quotients divide by the separation of the representable + * evaluation points. Refinement stops when successive tableau diagonals differ by at least + * twice the best estimated error, when evaluation points coincide or repeat the previous pair, + * or when a non-finite value is encountered. The best finite estimate is retained; + * if none exists, the column contains NaNs. + * + * \param[in] x The evaluation point. + * \param[out] jac The Jacobian at \a x. + * \returns the number of functor evaluations. + * \sa Parameters + */ + int df(const InputType& x, JacobianType& jac) const { return dfImpl(x, jac, nullptr); } + + /** + * \brief Computes the Jacobian and an absolute error estimate per input coordinate. + * + * The estimate is the maximum componentwise difference from the two parent tableau entries; + * it is a heuristic, not an error bound. An infinite estimate means no extrapolation was available. + * + * \param[in] x The evaluation point. + * \param[out] jac The Jacobian at \a x, sized by the caller to values() by x.size(). + * \param[out] error The error estimates, resized to x.size(). + * \returns the number of functor evaluations. + * \sa Parameters + */ + int df(const InputType& x, JacobianType& jac, InputType& error) const { + error.resize(x.size()); + return dfImpl(x, jac, &error); + } + + private: + Functor m_functor; + Parameters m_parameters; + + int dfImpl(const InputType& input, JacobianType& jac, InputType* error) const { + EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex && !NumTraits<Scalar>::IsInteger, + NUMERICAL_DIFF_REQUIRES_A_REAL_FLOATING_POINT_SCALAR); + // C. J. F. Ridders, "Accurate computation of F'(x) and F'(x)F''(x)", + // Advances in Engineering Software 4(2), 75-76 (1982). Extrapolate central differences in h^2. + using Tableau = Matrix<Scalar, Functor::ValuesAtCompileTime, Dynamic>; + const Index values = m_functor.values(); + eigen_assert(jac.rows() == values && jac.cols() == input.size()); + if (values == 0 || input.size() == 0) { + if (error) error->setZero(); + return 0; + } + + Tableau workspace(values, 2 * Index(m_parameters.maxSteps)); + auto first = workspace.leftCols(m_parameters.maxSteps); + auto second = workspace.rightCols(m_parameters.maxSteps); + auto current = &first; + auto previous = &second; + + // Functors may require ValueType& outputs, which cannot bind to column views. + ValueType plus(values), minus(values); + InputType x = input; + const Scalar reductionSquared = m_parameters.stepReduction * m_parameters.stepReduction; + int evaluations = 0; + + for (Index j = 0; j < input.size(); ++j) { + Scalar h = m_parameters.initialStep * numext::maxi(numext::abs(input[j]), Scalar(1)); + Scalar previousPlus = input[j], previousMinus = input[j]; + Scalar bestError = NumTraits<Scalar>::infinity(); + jac.col(j).setConstant(NumTraits<Scalar>::quiet_NaN()); + for (int step = 0; step < m_parameters.maxSteps; ++step) { + Scalar xPlus = input[j] + h; + Scalar xMinus = input[j] - h; + internal::numerical_diff_barrier(xPlus); + internal::numerical_diff_barrier(xMinus); + if (!(numext::isfinite)(xPlus) || !(numext::isfinite)(xMinus)) break; + if (step > 0 && xPlus == previousPlus && xMinus == previousMinus) break; + + const Scalar separation = xPlus - xMinus; + if (!(separation > Scalar(0)) || !(numext::isfinite)(separation)) break; + previousPlus = xPlus; + previousMinus = xMinus; + + x[j] = xPlus; + m_functor(x, plus); + x[j] = xMinus; + m_functor(x, minus); + x[j] = input[j]; + evaluations += 2; + + current->col(0) = (plus - minus) / separation; + if (!current->col(0).allFinite()) { + if (separation <= Scalar(1)) break; + // Divide before subtracting to avoid overflow, preserving already finite quotients. + auto quotient = current->col(0).array(); + quotient = quotient.isFiniteTyped().select(quotient, plus.array() / separation - minus.array() / separation); + if (!current->col(0).allFinite()) break; + } + Index bestOrder = step == 0 ? 0 : -1; + + Scalar factor = reductionSquared; + bool finite = true; + for (int order = 1; order <= step; ++order) { + // A(k,i) = (t^(2k) * A(k-1,i) - A(k-1,i-1)) / (t^(2k) - 1). + // Avoid forming `t^(2k) * A(k-1,i)` directly, since it can overflow. + current->col(order) = + current->col(order - 1) + (current->col(order - 1) - previous->col(order - 1)) / (factor - Scalar(1)); + if (!current->col(order).allFinite()) { + finite = false; + break; + } + + const Scalar estimate = + numext::maxi((current->col(order) - current->col(order - 1)).template lpNorm<Infinity>(), + (current->col(order) - previous->col(order - 1)).template lpNorm<Infinity>()); + if ((numext::isfinite)(estimate) && estimate <= bestError) { + bestError = estimate; + bestOrder = order; + } + + factor *= reductionSquared; + } + if (bestOrder >= 0) jac.col(j) = current->col(bestOrder); + if (!finite) break; + if (step > 0 && + (current->col(step) - previous->col(step - 1)).template lpNorm<Infinity>() >= Scalar(2) * bestError) + break; + + std::swap(current, previous); + h /= m_parameters.stepReduction; + } + if (error) (*error)[j] = bestError; + } + return evaluations; + } +}; + +#if EIGEN_COMP_CXXVER < 17 +template <typename Functor> +constexpr int Ridders<Functor>::InputsAtCompileTime; +template <typename Functor> +constexpr int Ridders<Functor>::ValuesAtCompileTime; +#endif + +} // namespace Eigen + +#endif // EIGEN_NUMERICALDIFF_RIDDERS_H
diff --git a/contrib/benchmarks/CMakeLists.txt b/contrib/benchmarks/CMakeLists.txt index 2bc57a4..d037aff 100644 --- a/contrib/benchmarks/CMakeLists.txt +++ b/contrib/benchmarks/CMakeLists.txt
@@ -36,6 +36,7 @@ add_subdirectory(Splines) add_subdirectory(Polynomials) add_subdirectory(KroneckerProduct) +add_subdirectory(NumericalDiff) add_subdirectory(StructuredMatrices) # GPU benchmarks have their own CMake project (needs CUDAToolkit).
diff --git a/contrib/benchmarks/NumericalDiff/CMakeLists.txt b/contrib/benchmarks/NumericalDiff/CMakeLists.txt new file mode 100644 index 0000000..6ae891e --- /dev/null +++ b/contrib/benchmarks/NumericalDiff/CMakeLists.txt
@@ -0,0 +1,4 @@ +# SPDX-FileCopyrightText: The Eigen Authors +# SPDX-License-Identifier: MPL-2.0 + +eigen_add_benchmark(bench_ridders bench_ridders.cpp)
diff --git a/contrib/benchmarks/NumericalDiff/bench_ridders.cpp b/contrib/benchmarks/NumericalDiff/bench_ridders.cpp new file mode 100644 index 0000000..85a51cd --- /dev/null +++ b/contrib/benchmarks/NumericalDiff/bench_ridders.cpp
@@ -0,0 +1,81 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <contrib/Eigen/NumericalDiff> +#include <algorithm> +#include <cmath> + +template <typename Scalar_, int Options> +struct CubicFunctor { + using Scalar = Scalar_; + using InputType = Eigen::Matrix<Scalar, 2, 1>; + using ValueType = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>; + using JacobianType = Eigen::Matrix<Scalar, Eigen::Dynamic, 2, Options>; + + static constexpr int InputsAtCompileTime = 2; + static constexpr int ValuesAtCompileTime = Eigen::Dynamic; + + Eigen::Index count; + + Eigen::Index values() const { return count; } + + int operator()(const InputType& x, ValueType& y) const { + for (Eigen::Index i = 0; i < count; ++i) y[i] = Scalar(i + 1) * x[0] * x[0] * x[0] + x[1] * x[1] * x[1]; + return 0; + } +}; + +template <typename Derivative> +static void measureJacobian(benchmark::State& state, Derivative& derivative, typename Derivative::Scalar tolerance) { + using Functor = typename Derivative::Functor; + using Scalar = typename Functor::Scalar; + const Eigen::Index count = derivative.values(); + typename Functor::InputType x; + x << Scalar(0.5), Scalar(-1); + typename Functor::JacobianType jac(count, 2); + const int evaluations = derivative.df(x, jac); + long double relativeError = 0; + for (Eigen::Index i = 0; i < count; ++i) { + const Scalar expected = Scalar(i + 1) * Scalar(0.75); + relativeError = (std::max)(relativeError, std::abs(static_cast<long double>(jac(i, 0)) / (0.75L * (i + 1)) - 1)); + relativeError = (std::max)(relativeError, std::abs(static_cast<long double>(jac(i, 1)) / 3 - 1)); + const Scalar bound = tolerance * Eigen::numext::maxi(expected, Scalar(3)); + if (!(Eigen::numext::abs(jac(i, 0) - expected) <= bound && Eigen::numext::abs(jac(i, 1) - Scalar(3)) <= bound)) { + state.SkipWithError("Incorrect cubic Jacobian"); + return; + } + } + state.counters["evaluations"] = evaluations; + state.counters["max_relative_error"] = static_cast<double>(relativeError); + + for (auto _ : state) { + benchmark::DoNotOptimize(derivative.df(x, jac)); + benchmark::DoNotOptimize(jac.data()); + benchmark::ClobberMemory(); + } +} + +template <typename Scalar, int Options> +static void BM_Ridders(benchmark::State& state) { + using Functor = CubicFunctor<Scalar, Options>; + Eigen::Ridders<Functor> derivative(Functor{state.range(0)}); + typename Eigen::Ridders<Functor>::Parameters parameters; + parameters.maxSteps = static_cast<int>(state.range(1)); + derivative.setParameters(parameters); + measureJacobian(state, derivative, Scalar(512) * Eigen::NumTraits<Scalar>::epsilon()); +} + +template <typename Scalar, int Options> +static void BM_Central(benchmark::State& state) { + using Functor = CubicFunctor<Scalar, Options>; + Eigen::NumericalDiff<Functor, Eigen::Central> derivative(Functor{state.range(0)}); + measureJacobian(state, derivative, Scalar(16) * std::sqrt(Eigen::NumTraits<Scalar>::epsilon())); +} + +BENCHMARK_TEMPLATE(BM_Ridders, float, Eigen::ColMajor)->ArgsProduct({{1, 16, 128, 1024}, {2, 10}}); +BENCHMARK_TEMPLATE(BM_Ridders, double, Eigen::ColMajor)->ArgsProduct({{1, 16, 128, 1024}, {2, 10}}); +BENCHMARK_TEMPLATE(BM_Ridders, double, Eigen::RowMajor)->ArgsProduct({{1, 16, 128, 1024}, {2, 10}}); +BENCHMARK_TEMPLATE(BM_Central, float, Eigen::ColMajor)->Arg(1)->Arg(16)->Arg(128)->Arg(1024); +BENCHMARK_TEMPLATE(BM_Central, double, Eigen::ColMajor)->Arg(1)->Arg(16)->Arg(128)->Arg(1024); +BENCHMARK_TEMPLATE(BM_Central, double, Eigen::RowMajor)->Arg(1)->Arg(16)->Arg(128)->Arg(1024);
diff --git a/contrib/test/CMakeLists.txt b/contrib/test/CMakeLists.txt index 3186340..1e5db5d 100644 --- a/contrib/test/CMakeLists.txt +++ b/contrib/test/CMakeLists.txt
@@ -39,6 +39,7 @@ ei_add_test(NumericalDiff) ei_add_test(stencil) +ei_add_test(ridders) ei_add_test(structured_dpr1) ei_add_test(autodiff_scalar) ei_add_test(autodiff)
diff --git a/contrib/test/ridders.cpp b/contrib/test/ridders.cpp new file mode 100644 index 0000000..df350d5 --- /dev/null +++ b/contrib/test/ridders.cpp
@@ -0,0 +1,552 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include "main.h" +#include <contrib/Eigen/NumericalDiff> +#include <contrib/Eigen/LevenbergMarquardt> + +template <typename Scalar_, int Inputs = Dynamic, int Values = Dynamic, int Options = ColMajor> +struct RiddersFunctor { + using Scalar = Scalar_; + using InputType = Matrix<Scalar, Inputs, 1>; + using ValueType = Matrix<Scalar, Values, 1>; + using JacobianType = Matrix<Scalar, Values, Inputs, Options>; + + static constexpr int InputsAtCompileTime = Inputs; + static constexpr int ValuesAtCompileTime = Values; + + mutable int evaluations = 0; + + int values() const { return 4; } + + int operator()(const InputType& x, ValueType& y) const { + ++evaluations; + y[0] = numext::exp(x[0]); + y[1] = numext::sin(Scalar(50) * x[0]); + y[2] = x[0] * x[1]; + y[3] = numext::exp(x[1]); + return 0; + } +}; + +template <typename Scalar, int Inputs, int Values, int Options> +void test_ridders_jacobian() { + using F = RiddersFunctor<Scalar, Inputs, Values, Options>; + + Ridders<F> diff; + typename F::InputType x(2), error(2); + typename F::JacobianType jac(4, 2), expected(4, 2), withoutError(4, 2); + + x << Scalar(1), Scalar(-2); + const typename F::InputType original = x; + + expected.row(0) << Scalar(std::exp(1.0L)), Scalar(0); + expected.row(1) << Scalar(50 * std::cos(50.0L)), Scalar(0); + expected.row(2) << x[1], x[0]; + expected.row(3) << Scalar(0), Scalar(std::exp(-2.0L)); + + const int count = diff.df(x, jac, error); + VERIFY_IS_EQUAL(count, diff.functor().evaluations); + VERIFY(count >= 8 && count <= 40 && count % 2 == 0); + + VERIFY(jac.allFinite()); + VERIFY(error.allFinite()); + VERIFY((error.array() >= Scalar(0)).all()); + + // The initial stencil spans half a radian in the oscillatory component; extrapolation + // must remove its truncation error to reach a roundoff-scale bound. + VERIFY((jac - expected).template lpNorm<Infinity>() <= + Scalar(256) * NumTraits<Scalar>::epsilon() * expected.template lpNorm<Infinity>()); + + VERIFY_IS_EQUAL(x, original); + VERIFY_IS_EQUAL(diff.df(x, withoutError), count); + VERIFY_IS_EQUAL(jac, withoutError); +} + +struct RiddersScalarFunctor { + using Scalar = double; + using InputType = Matrix<double, 1, 1>; + using ValueType = InputType; + using JacobianType = InputType; + + static constexpr int InputsAtCompileTime = 1; + static constexpr int ValuesAtCompileTime = 1; + + enum Function { Cubic, LargeCubic, Quintic, SquareRoot, NoisyExp, Identity, Invalid, InvalidRefinement }; + Function function = Cubic; + mutable int evaluations = 0; + + int values() const { return 1; } + + int operator()(const InputType& x, ValueType& y) const { + ++evaluations; + switch (function) { + case Cubic: + y[0] = x[0] * x[0] * x[0]; + break; + + case LargeCubic: + y[0] = 2e307 * x[0] * x[0] * x[0]; + break; + + case Quintic: + y[0] = x[0] * x[0] * x[0] * x[0] * x[0]; + break; + + case SquareRoot: + y[0] = std::sqrt(x[0]); + break; + + case NoisyExp: + y[0] = double(float(std::exp(x[0]))); + break; + + case Identity: + y[0] = x[0]; + break; + + case Invalid: + y[0] = NumTraits<double>::quiet_NaN(); + break; + + case InvalidRefinement: + y[0] = evaluations <= 2 ? x[0] : NumTraits<double>::quiet_NaN(); + break; + } + return 0; + } +}; + +void test_ridders_parameters() { + using DiffType = Ridders<RiddersScalarFunctor>; + + DiffType diff; + DiffType::Parameters parameters; + parameters.initialStep = 0.25; + parameters.stepReduction = 4; + parameters.maxSteps = 2; + diff.setParameters(parameters); + + RiddersScalarFunctor::InputType x, error; + RiddersScalarFunctor::JacobianType jac; + x[0] = 2; + VERIFY_IS_EQUAL(diff.df(x, jac, error), 4); + VERIFY_IS_EQUAL(jac[0], 12.0); + + // The cubic's central differences are 12 + h^2, with h = 1/2 and 1/8. + VERIFY_IS_EQUAL(error[0], 0.25); + + parameters.maxSteps = 10; + diff.setParameters(parameters); + VERIFY(diff.df(x, jac, error) < 20); + VERIFY_IS_EQUAL(jac[0], 12.0); + VERIFY_IS_EQUAL(error[0], 0.0); + + parameters.initialStep = 0; + VERIFY_RAISES_ASSERT(diff.setParameters(parameters)); + + parameters.initialStep = 0.01; + parameters.stepReduction = 1; + VERIFY_RAISES_ASSERT(diff.setParameters(parameters)); + + parameters.stepReduction = 2; + parameters.maxSteps = 1; + VERIFY_RAISES_ASSERT(diff.setParameters(parameters)); + + parameters.maxSteps = 10; + parameters.initialStep = NumTraits<double>::infinity(); + VERIFY_RAISES_ASSERT(diff.setParameters(parameters)); + + parameters.initialStep = 0.01; + parameters.stepReduction = NumTraits<double>::infinity(); + VERIFY_RAISES_ASSERT(diff.setParameters(parameters)); + + parameters = DiffType::Parameters(); + if (NumTraits<Index>::highest() / 2 < NumTraits<int>::highest()) { + parameters.maxSteps = static_cast<int>(NumTraits<Index>::highest() / 2) + 1; + VERIFY_RAISES_ASSERT(diff.setParameters(parameters)); + --parameters.maxSteps; + } else { + parameters.maxSteps = NumTraits<int>::highest(); + } + diff.setParameters(parameters); +} + +void test_ridders_repeated_stencil() { + using DiffType = Ridders<RiddersScalarFunctor>; + + for (int multiplier : {1, 16}) { + DiffType diff; + DiffType::Parameters parameters; + parameters.stepReduction = 1 + multiplier * NumTraits<double>::epsilon(); + diff.setParameters(parameters); + + RiddersScalarFunctor::InputType x, error; + RiddersScalarFunctor::JacobianType jac; + x[0] = 1; + VERIFY_IS_EQUAL(diff.df(x, jac, error), 2); + VERIFY_IS_EQUAL(diff.functor().evaluations, 2); + VERIFY((numext::isfinite)(jac[0])); + VERIFY(std::abs(jac[0] - 3.0001) <= 32 * NumTraits<double>::epsilon()); + VERIFY((numext::isinf)(error[0]) && error[0] > 0); + } +} + +void test_ridders_dimensions() { + Ridders<RiddersFunctor<double>> diff; + VectorXd x = VectorXd::Ones(2), error; + MatrixXd jac(3, 2); + VERIFY_RAISES_ASSERT(diff.df(x, jac)); + + jac.resize(4, 1); + VERIFY_RAISES_ASSERT(diff.df(x, jac, error)); + + jac.resize(4, 0); + VERIFY_RAISES_ASSERT(diff.df(x, jac)); + VERIFY_IS_EQUAL(diff.functor().evaluations, 0); +} + +template <typename Scalar_> +struct RiddersAffineFunctor { + using Scalar = Scalar_; + using InputType = Matrix<Scalar, 1, 1>; + using ValueType = Matrix<Scalar, Dynamic, 1>; + using JacobianType = Matrix<Scalar, Dynamic, 1>; + + static constexpr int InputsAtCompileTime = 1; + static constexpr int ValuesAtCompileTime = Dynamic; + + ValueType slopes, offsets; + Scalar center = Scalar(100); + Scalar inputScale = Scalar(1); + mutable int evaluations = 0; + int finiteEvaluations = 20; + + Index values() const { return slopes.size(); } + + int operator()(const InputType& x, ValueType& y) const { + ++evaluations; + y = slopes * ((x[0] - center) * inputScale) + offsets; + if (evaluations > finiteEvaluations) y.setConstant(NumTraits<Scalar>::quiet_NaN()); + return 0; + } +}; + +template <typename Scalar> +void test_ridders_central_overflow() { + using Functor = RiddersAffineFunctor<Scalar>; + using DiffType = Ridders<Functor>; + + const Scalar large = NumTraits<Scalar>::highest() * Scalar(0.75); + const Scalar tiny = (std::numeric_limits<Scalar>::denorm_min)(); + + for (Index count : {1, 7, 17}) { + Functor functor; + functor.slopes.resize(count); + functor.offsets.setZero(count); + + for (Index i = 0; i < count; ++i) { + functor.slopes[i] = i % 3 == 0 ? large : (i % 3 == 1 ? Scalar(-8) : (std::numeric_limits<Scalar>::min)()); + if (i % 6 == 3) functor.slopes[i] = -large; + } + + DiffType diff(functor); + typename Functor::InputType x, error; + typename Functor::JacobianType jac(count, 1); + x[0] = functor.center; + VERIFY_IS_EQUAL(diff.df(x, jac, error), 4); + VERIFY(jac.allFinite()); + for (Index i = 0; i < count; ++i) VERIFY_IS_EQUAL(jac[i], functor.slopes[i]); + VERIFY_IS_EQUAL(error[0], Scalar(0)); + + // A failed refinement must retain finite entries, including subnormals, from the initial stencil. + diff.functor().evaluations = 0; + diff.functor().finiteEvaluations = 2; + for (Index i = 2; i < count; i += 3) diff.functor().slopes[i] = tiny; + VERIFY_IS_EQUAL(diff.df(x, jac, error), 4); + VERIFY_IS_EQUAL(jac, diff.functor().slopes); + VERIFY((numext::isinf)(error[0]) && error[0] > Scalar(0)); + } + + Functor functor; + functor.slopes.setZero(1); + functor.offsets.setConstant(1, large); + + DiffType diff(functor); + typename Functor::InputType x, error; + typename Functor::JacobianType jac(1, 1); + x[0] = Scalar(100); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 4); + VERIFY_IS_EQUAL(jac[0], Scalar(0)); + + x[0] = Scalar(0); + diff.functor().center = Scalar(0); + + typename DiffType::Parameters parameters = diff.parameters(); + parameters.initialStep = (std::numeric_limits<Scalar>::min)(); + diff.setParameters(parameters); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 4); + VERIFY_IS_EQUAL(jac[0], Scalar(0)); + + // Division first cannot repair a genuinely overflowing derivative or nonfinite samples. + diff.functor().slopes[0] = NumTraits<Scalar>::highest(); + diff.functor().offsets[0] = Scalar(0); + diff.functor().inputScale = Scalar(4); + parameters.initialStep = Scalar(0.25); + diff.setParameters(parameters); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 2); + VERIFY((numext::isnan)(jac[0])); + VERIFY((numext::isinf)(error[0]) && error[0] > Scalar(0)); + + diff.functor().slopes[0] = NumTraits<Scalar>::infinity(); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 2); + VERIFY((numext::isnan)(jac[0])); + VERIFY((numext::isinf)(error[0]) && error[0] > Scalar(0)); + + parameters.initialStep = Scalar(2); + diff.setParameters(parameters); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 2); + VERIFY((numext::isnan)(jac[0])); +} + +void test_ridders_domain_and_noise() { + using DiffType = Ridders<RiddersScalarFunctor>; + + DiffType diff; + RiddersScalarFunctor::InputType x, error; + RiddersScalarFunctor::JacobianType jac; + + diff.functor().function = RiddersScalarFunctor::SquareRoot; + x[0] = 0.005; + + DiffType::Parameters parameters; + parameters.initialStep = 0.01 * x[0]; + diff.setParameters(parameters); + VERIFY(diff.df(x, jac, error) <= 20); + + const double expected = double(0.5L / std::sqrt(static_cast<long double>(x[0]))); + VERIFY(std::abs(jac[0] - expected) <= 256 * NumTraits<double>::epsilon() * expected); + VERIFY((numext::isfinite)(error[0])); + + diff.functor().function = RiddersScalarFunctor::NoisyExp; + x[0] = 1; + parameters.initialStep = 0.01; + diff.setParameters(parameters); + VERIFY(diff.df(x, jac, error) < 20); + VERIFY(std::abs(jac[0] - std::exp(1.0)) <= 256 * double(NumTraits<float>::epsilon()) * std::exp(1.0)); + VERIFY((numext::isfinite)(error[0])); +} + +void test_ridders_richardson_overflow() { + Ridders<RiddersScalarFunctor> diff; + diff.functor().function = RiddersScalarFunctor::LargeCubic; + RiddersScalarFunctor::InputType x, error; + RiddersScalarFunctor::JacobianType jac; + x[0] = 1; + + const int evaluations = diff.df(x, jac, error); + VERIFY_IS_EQUAL(evaluations, diff.functor().evaluations); + VERIFY(evaluations > 4 && evaluations <= 20); + + VERIFY((numext::isfinite)(jac[0])); + VERIFY((numext::isfinite)(error[0])); + VERIFY(error[0] >= 0); + + // The cubic's truncation error cancels on extrapolation; divide before comparing near the range limit. + VERIFY(std::abs(jac[0] / 6e307 - 1) <= 256 * NumTraits<double>::epsilon()); +} + +void test_ridders_stencil() { + using DiffType = Ridders<RiddersScalarFunctor>; + + DiffType diff; + diff.functor().function = RiddersScalarFunctor::Quintic; + + DiffType::Parameters parameters; + parameters.initialStep = 0.25; + parameters.maxSteps = 3; + diff.setParameters(parameters); + + RiddersScalarFunctor::InputType x, error, value; + RiddersScalarFunctor::JacobianType jac; + x[0] = 1; + VERIFY_IS_EQUAL(diff.df(x, jac, error), 6); + VERIFY_IS_EQUAL(jac[0], 5.0); + + // At depth three the extrapolation is the degree-five-exact stencil on these six nodes. + const Stencil<1, double, 6> stencil({-0.25, 0.25, -0.125, 0.125, -0.0625, 0.0625}); + double expected = 0; + for (Index i = 0; i < 6; ++i) { + RiddersScalarFunctor::InputType point; + point[0] = x[0] + stencil.point(i); + diff(point, value); + expected += stencil.weight(i) * value[0]; + } + + VERIFY(std::abs(jac[0] - expected) <= 32 * NumTraits<double>::epsilon() * std::abs(jac[0])); +} + +// Tests construction from a functor without a default constructor and optional inputs() forwarding. +struct OwnedRiddersFunctor : RiddersScalarFunctor { + explicit OwnedRiddersFunctor(Function selected) { function = selected; } + + int inputs() const { return 1; } +}; + +void test_ridders_composition() { + using Derivative = Ridders<OwnedRiddersFunctor>; + + static_assert(!std::is_base_of<OwnedRiddersFunctor, Derivative>::value, "Ridders owns its functor by composition"); + + OwnedRiddersFunctor functor(RiddersScalarFunctor::Cubic); + Derivative::Parameters parameters; + parameters.initialStep = 0.25; + parameters.maxSteps = 2; + Derivative derivative(functor, parameters); + VERIFY_IS_EQUAL(derivative.inputs(), 1); + VERIFY_IS_EQUAL(derivative.values(), 1); + VERIFY_IS_EQUAL(derivative.parameters().maxSteps, 2); + + const int* inputs = &Derivative::InputsAtCompileTime; + const int* values = &Derivative::ValuesAtCompileTime; + VERIFY_IS_EQUAL(*inputs, 1); + VERIFY_IS_EQUAL(*values, 1); + + OwnedRiddersFunctor::InputType x, y, error; + OwnedRiddersFunctor::JacobianType jac; + x[0] = 2; + VERIFY_IS_EQUAL(derivative(x, y), 0); + VERIFY_IS_EQUAL(y[0], 8.0); + + VERIFY_IS_EQUAL(derivative.df(x, jac, error), 4); + VERIFY_IS_EQUAL(jac[0], 12.0); + VERIFY_IS_EQUAL(functor.evaluations, 0); + + const Derivative& constDerivative = derivative; + VERIFY_IS_EQUAL(constDerivative.functor().evaluations, 5); +} + +struct LinearRiddersFunctor : RiddersFunctor<double, 7, 7> { + int values() const { return 7; } + + int operator()(const InputType& x, ValueType& y) const { + y = x * 8; + return 0; + } +}; + +void test_ridders_step_alignment() { + Ridders<LinearRiddersFunctor> derivative; + LinearRiddersFunctor::InputType x; + x << 0.0, 0.082, 1.13, -2350.7, 7.7e8 + 0.3, -1.0 - std::ldexp(1.0, -27), -1.0; + LinearRiddersFunctor::JacobianType jac; + + derivative.df(x, jac); + VERIFY_IS_EQUAL(jac, LinearRiddersFunctor::JacobianType::Identity() * 8); +} + +void test_ridders_fallback() { + using DiffType = Ridders<RiddersScalarFunctor>; + + DiffType diff; + RiddersScalarFunctor::InputType x, error; + RiddersScalarFunctor::JacobianType jac; + x[0] = 1; + + diff.functor().function = RiddersScalarFunctor::Invalid; + VERIFY_IS_EQUAL(diff.df(x, jac, error), 2); + VERIFY((numext::isnan)(jac[0])); + VERIFY((numext::isinf)(error[0]) && error[0] > 0); + + diff.functor().function = RiddersScalarFunctor::InvalidRefinement; + diff.functor().evaluations = 0; + VERIFY_IS_EQUAL(diff.df(x, jac, error), 4); + VERIFY_IS_EQUAL(jac[0], 1.0); + VERIFY((numext::isinf)(error[0])); + + diff.functor().function = RiddersScalarFunctor::Identity; + + DiffType::Parameters parameters; + parameters.initialStep = NumTraits<double>::epsilon() / 16; + diff.setParameters(parameters); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 0); + VERIFY((numext::isnan)(jac[0])); + VERIFY((numext::isinf)(error[0])); + + parameters.initialStep = NumTraits<double>::epsilon(); + diff.setParameters(parameters); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 4); + VERIFY_IS_EQUAL(jac[0], 1.0); + VERIFY_IS_EQUAL(error[0], 0.0); +} + +struct RiddersLevenbergMarquardtFunctor : DenseFunctor<double> { + RiddersLevenbergMarquardtFunctor() : DenseFunctor<double>(1, 1) {} + + int operator()(const InputType& x, ValueType& y) const { + y[0] = x[0] * x[0] - 2; + return 0; + } +}; + +void test_ridders_levenberg_marquardt() { + using Functor = RiddersLevenbergMarquardtFunctor; + using Derivative = Ridders<Functor>; + + static_assert(std::is_same<typename Derivative::QRSolver, typename Functor::QRSolver>::value, + "Ridders must preserve the functor's QRSolver type"); + + Derivative derivative; + LevenbergMarquardt<Derivative> solver(derivative); + VectorXd x = VectorXd::Ones(1); + + const LevenbergMarquardtSpace::Status status = solver.minimize(x); + VERIFY(status >= LevenbergMarquardtSpace::RelativeReductionTooSmall && + status <= LevenbergMarquardtSpace::CosinusTooSmall); + VERIFY(std::abs(x[0] - std::sqrt(2.0)) <= 256 * NumTraits<double>::epsilon()); +} + +struct EmptyRiddersFunctor : RiddersFunctor<double> { + int values() const { return 0; } +}; + +void test_ridders_empty() { + Ridders<EmptyRiddersFunctor> diff; + VectorXd x = VectorXd::Ones(2), error; + MatrixXd jac(0, 2); + VERIFY_IS_EQUAL(diff.df(x, jac, error), 0); + VERIFY_IS_EQUAL(error, VectorXd::Zero(2)); + VERIFY_IS_EQUAL(diff.functor().evaluations, 0); + + Ridders<RiddersFunctor<double>> noInputs; + x.resize(0); + jac.resize(4, 0); + VERIFY_IS_EQUAL(noInputs.df(x, jac, error), 0); + VERIFY_IS_EQUAL(error.size(), 0); + VERIFY_IS_EQUAL(noInputs.functor().evaluations, 0); +} + +EIGEN_DECLARE_TEST(ridders) { + CALL_SUBTEST((test_ridders_jacobian<float, 2, 4, ColMajor>())); + CALL_SUBTEST((test_ridders_jacobian<double, 2, 4, RowMajor>())); + CALL_SUBTEST((test_ridders_jacobian<double, Dynamic, Dynamic, ColMajor>())); + CALL_SUBTEST((test_ridders_jacobian<long double, Dynamic, Dynamic, RowMajor>())); + + CALL_SUBTEST(test_ridders_parameters()); + CALL_SUBTEST(test_ridders_repeated_stencil()); + CALL_SUBTEST(test_ridders_dimensions()); + + CALL_SUBTEST(test_ridders_central_overflow<float>()); + CALL_SUBTEST(test_ridders_central_overflow<double>()); + + CALL_SUBTEST(test_ridders_domain_and_noise()); + CALL_SUBTEST(test_ridders_richardson_overflow()); + CALL_SUBTEST(test_ridders_stencil()); + + CALL_SUBTEST(test_ridders_composition()); + CALL_SUBTEST(test_ridders_step_alignment()); + CALL_SUBTEST(test_ridders_fallback()); + + CALL_SUBTEST(test_ridders_levenberg_marquardt()); + CALL_SUBTEST(test_ridders_empty()); +}