From 7752c262fa265a6f51c931d2424a53306c1b44e1 Mon Sep 17 00:00:00 2001 From: Erik Frojdh Date: Wed, 30 Sep 2026 15:47:50 +0200 Subject: [PATCH] first attempt at simplifying --- AGENTS.md | 22 +- RELEASE.md | 21 +- include/aare/Fit.hpp | 15 +- include/aare/FitModel.hpp | 22 +- include/aare/Models.hpp | 1045 ++++++++++---------------------- python/src/fit.hpp | 251 +++----- src/fit/DampedGaussNewton.hpp | 568 ++++++++--------- src/fit/Fit.cpp | 94 ++- src/fit/Fit.test.cpp | 62 +- src/fit/FitHelpers.hpp | 35 +- src/fit/FitMinuit2.cpp | 34 +- src/fit/FitMinuit2.hpp | 10 +- src/fit/LevenbergMarquardt.hpp | 85 +-- src/fit/LinearAlgebra.hpp | 173 ++++++ src/fit/MinuitSteps.hpp | 200 +++--- src/fit/VariableProjection.hpp | 333 ++++------ 16 files changed, 1195 insertions(+), 1775 deletions(-) create mode 100644 src/fit/LinearAlgebra.hpp diff --git a/AGENTS.md b/AGENTS.md index 5646aaae..a8a322b0 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -149,18 +149,26 @@ lives in `python/aare/transform.py` on top of C++ generators. `GaussianChargeSharing`, ...) form the public fitting API. Each model provides `eval`, `eval_and_grad`, `is_valid`, `estimate_par`, and `param_info`; the separable models also declare their linear parameters -(`linear_par`, `basis_and_grad`, `basis_columns`). `FitModel` is plain data +(`linear_par`, `basis_columns`). A model is listed once, in +`AARE_FOR_EACH_FIT_MODEL` at the end of `Models.hpp`, which generates the +explicit instantiations and the Python bindings; it also needs its step +scales in `MinuitSteps.hpp`. `FitModel` is plain data (limits, fixed flags, start values, minimizer settings) and selects one of four minimizers through `aare::Minimizer`: Minuit2's Migrad (default) and Fumili, or the built-in Levenberg-Marquardt and variable projection solvers. The private implementation lives in `src/fit/`: `Fit.cpp` holds the `FitModel` method bodies, the per-thread `detail::PixelFitter` that -dispatches on the minimizer and falls back from VarPro to -Levenberg-Marquardt, `fit_pixel`/`fit_3d`, and the explicit -instantiations; `FitHelpers.hpp` holds the start-value precedence and -shape checks; `DampedGaussNewton.hpp` is the damped Gauss-Newton driver -shared by `LevenbergMarquardt.hpp` and `VariableProjection.hpp`, with the -small fixed-size linear algebra on `std::array` (`Vec`, `Mat`). +dispatches on the minimizer, falls back from VarPro to +Levenberg-Marquardt and zeroes the outputs of a failed fit, +`fit_pixel`/`fit_3d`, and the explicit instantiations; `FitHelpers.hpp` +holds the start-value precedence, the evaluation budget and the shape +checks; `DampedGaussNewton.hpp` is the damped Gauss-Newton driver +shared by `LevenbergMarquardt.hpp` and `VariableProjection.hpp`; +`LinearAlgebra.hpp` holds their small fixed-size linear algebra on +`std::array` (`Vec`, `Mat`, `SymMat`, `Cholesky`, +`compress`). The per-point loops of `VariableProjection::reduce` are +sensitive to the code around them, see the comment there; check instruction +counts with `perf stat` when touching it. `FitMinuit2.cpp` (with `Chi2.hpp` and `MinuitSteps.hpp`) is the only translation unit that includes Minuit2 headers; Minuit2 is a private, `BUILD_INTERFACE`-only dependency of `aare_core`. The `[fit]` tests run diff --git a/RELEASE.md b/RELEASE.md index 479c3e48..2f798d94 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -55,7 +55,7 @@ to exhaust the default budget of 100 may now converge. - Added ``Minimizer.VarPro``, a built-in variable projection solver. Every bundled model now declares the parameters it is linear in - (``linear_par`` and ``basis_and_grad`` in C++, checked against the model + (``linear_par`` and ``basis_columns`` in C++, checked against the model derivatives by the tests). Each trial point solves those exactly and the Levenberg-Marquardt iteration runs over the nonlinear parameters only, so the fit needs fewer evaluations and does not depend on start values of the @@ -65,12 +65,13 @@ Limits on linear parameters, pixels that do not converge and models without the separable structure fall back to ``LevenbergMarquardt``. Parameter errors are the same Gauss-Newton - estimates. The Gaussian, plateau, charge-sharing and S-curve models also - provide their basis functions for all scan points at once - (``basis_columns``), a loop that compiles to vector instructions with the - new ``model::fast_exp``; ``eval`` and ``eval_and_grad``, and with them the - other minimizers, are unchanged. ``Minimizer.Migrad`` remains the - default. + estimates. ``basis_columns`` provides the basis functions for all scan + points at once, a loop that compiles to vector instructions with the new + ``model::fast_exp``. ``GaussianErfcPlateau``, ``GaussianChargeSharing`` and + ``GaussianChargeSharingKb`` now compute their peak and erfc terms in one + shared function (``model::peak_terms``), a rounding-level change of their + values and derivatives for every minimizer. ``Minimizer.Migrad`` remains + the default. - ``NDArray``/``NDView`` expressions now support scalar operands (``2 * a + b / 4``) and can be assigned to an existing ``NDArray``, reusing its buffer. @@ -93,6 +94,12 @@ ``std::out_of_range`` (``IndexError`` in Python) and ``SetParLimits`` rejects ``lo >= hi``. The models in ``Models.hpp`` no longer provide the Minuit specific ``compute_steps`` and ``compute_ranges`` helpers. + ``model::RisingScurve`` and ``model::FallingScurve`` now derive from one + template, ``model::Scurve<+1>`` and ``model::Scurve<-1>``. The list of fit + models lives in ``AARE_FOR_EACH_FIT_MODEL`` in ``Models.hpp``, from which + the explicit instantiations and the Python bindings are generated. +- C++ ``fit_3d`` accepts an empty ``err_out`` for a weighted fit too; it + used to require one whenever ``y_err`` was passed. - Fit behaviour common to all minimizers: user start and fixed values are applied before the validity check and free start values are clamped into their limits; one-sided limits are open-ended for ``LevenbergMarquardt`` diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index cc7edce1..d94d321a 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -8,17 +8,16 @@ namespace aare { // --------------------------------------------------------------------------- // Pixel fitting with the minimizer selected by the FitModel. -// Template bodies and explicit instantiations live in src/Fit.cpp. +// Template bodies and explicit instantiations live in src/fit/Fit.cpp. // --------------------------------------------------------------------------- /** * @brief Fit a single pixel's data. * - * The minimizer (Migrad, Fumili or LevenbergMarquardt) is selected by the - * model, see aare::Minimizer. Parameter errors, when requested, come from - * MnHesse with Migrad and from the linearised covariance otherwise; fixed - * parameters and parameters ending on a limit report 0. A failed fit - * returns zeros. + * The minimizer is selected by the model, see aare::Minimizer. Parameter + * errors, when requested, come from MnHesse with Migrad and from the linearised + * covariance otherwise; fixed parameters and parameters ending on a limit + * report 0. A failed fit returns zeros. * * User-precedence rules: * - Fixed parameters: kept at their value, excluded from the fit. @@ -57,8 +56,8 @@ NDArray fit_pixel(const FitModel &model, NDView x, * @param y_err Uncertainties, same shape as y, or empty for unweighted * fits. * @param par_out Output parameters, shape `(rows, cols, npar)`. - * @param err_out Output parameter errors, shape `(rows, cols, npar)`, if - * used. + * @param err_out Output parameter errors, shape `(rows, cols, npar)`, or + * empty when they are not wanted. Filled only with compute_errors. * @param chi2_out Output chi-squared / objective values, shape `(rows, * cols)`. * @param n_threads Number of threads used to split rows. diff --git a/include/aare/FitModel.hpp b/include/aare/FitModel.hpp index f485849e..18943424 100644 --- a/include/aare/FitModel.hpp +++ b/include/aare/FitModel.hpp @@ -22,12 +22,12 @@ namespace aare { * two-sided limit becomes active; a pixel without a valid Fumili minimum * is refitted with Migrad. * - LevenbergMarquardt: the built-in damped Gauss-Newton solver - * (src/LevenbergMarquardt.hpp). It uses the analytic Jacobian, reflects + * (src/fit/LevenbergMarquardt.hpp). It uses the analytic Jacobian, reflects * steps at parameter limits and reuses its buffers between pixels, so a * data cube is fitted without per-pixel allocations. Parameter errors are * sqrt(diag((J^T J)^-1)) over the free parameters. * - VarPro: the built-in variable projection solver - * (src/VariableProjection.hpp) for models that declare their linear + * (src/fit/VariableProjection.hpp) for models that declare their linear * parameters (model::is_separable). Every trial point solves the linear * parameters exactly and the Levenberg-Marquardt iteration runs over the * nonlinear ones only, which needs fewer evaluations and avoids poor @@ -58,8 +58,8 @@ enum class Minimizer { Migrad, Fumili, LevenbergMarquardt, VarPro }; * parameters, user start values and minimizer settings. * * The object holds plain data and can be shared read-only between threads. - * Method bodies live in src/Fit.cpp, which explicitly instantiates the class - * for every supported model. + * Method bodies live in src/fit/Fit.cpp, which explicitly instantiates the + * class for every supported model. */ template class FitModel { public: @@ -152,14 +152,10 @@ template class FitModel { }; // Suppress implicit instantiation for all supported model types. -// Definitions live in src/Fit.cpp. -extern template class FitModel; -extern template class FitModel; -extern template class FitModel; -extern template class FitModel; -extern template class FitModel; -extern template class FitModel; -extern template class FitModel; -extern template class FitModel; +// Definitions live in src/fit/Fit.cpp. +#define AARE_DECLARE_FIT_MODEL(Model) \ + extern template class FitModel; +AARE_FOR_EACH_FIT_MODEL(AARE_DECLARE_FIT_MODEL) +#undef AARE_DECLARE_FIT_MODEL } // namespace aare diff --git a/include/aare/Models.hpp b/include/aare/Models.hpp index da908f14..0d0dcd6b 100644 --- a/include/aare/Models.hpp +++ b/include/aare/Models.hpp @@ -111,53 +111,38 @@ inline constexpr double no_bound = std::numeric_limits::infinity(); * for Minimizer::VarPro. * * A model whose function is linear in some of its parameters lists them in - * ascending order in `linear_par` and provides - * - * static void basis_and_grad(double x, const std::vector &par, - * std::array &phi, - * std::array, nnl> &dphi) - * - * with nlin = linear_par.size() and nnl = npar - nlin, such that - * - * f(x; par) = sum_j par[linear_par[j]] * phi[j] - * - * and dphi[k][j] = d phi[j] / d par[nonlinear[k]], where nonlinear lists the - * remaining parameter indices in ascending order (see separable_traits). - * basis_and_grad ignores the linear entries of par. The [fit] tests check - * phi and dphi against eval and eval_and_grad for every model. - * - * A model may additionally provide the same quantities for all scan points - * at once, see has_basis_columns: + * ascending order in `linear_par` and provides its basis functions and their + * derivatives for all scan points at once: * * static void basis_columns(const double *__restrict x, ssize_t n, * const std::vector &par, * double *__restrict phi, * double *__restrict dphi) * - * writes phi[j * n + i] = phi_j(x_i) and dphi[(k * nlin + j) * n + i] = - * d phi_j / d par[nonlinear[k]] (x_i), one contiguous column per function. - * Written as a loop free of reductions with the parameters hoisted and - * fast_exp instead of std::exp, it compiles to vector instructions, which - * is where VarPro spends most of its time. The __restrict qualifiers are + * With nlin = linear_par.size() and nnl = npar - nlin it writes the nlin + * columns phi[j * n + i] = phi_j(x_i) of the basis functions, such that + * + * f(x; par) = sum_j par[linear_par[j]] * phi_j(x) + * + * and the nnl * nlin columns dphi[(k * nlin + j) * n + i] = + * d phi_j / d par[nonlinear[k]] (x_i), where nonlinear lists the remaining + * parameter indices in ascending order (see separable_traits). It ignores + * the linear entries of par. The [fit] tests check the columns against eval + * and eval_and_grad for every model. + * + * VarPro spends most of its time here, so the function is written to compile + * to vector instructions: one loop free of reductions, with the parameters + * hoisted and fast_exp instead of std::exp. The __restrict qualifiers are * needed: without them gcc has to prove that the output columns do not * overlap, exceeds its budget of run-time alias checks once a loop writes - * more than three columns, and quietly emits scalar code. The tests check - * it against basis_and_grad. + * more than three columns, and quietly emits scalar code. */ template struct is_separable : std::false_type {}; template struct is_separable> - : std::true_type {}; - -/** @brief Detects a separable model that provides basis_columns. */ -template -struct has_basis_columns : std::false_type {}; - -template -struct has_basis_columns> + decltype(&Model::basis_columns)>> : std::true_type {}; /** @brief The parameter indices not listed in linear, in ascending order. */ @@ -186,10 +171,107 @@ template struct separable_traits { static constexpr std::size_t linear_index(int j) { return Model::linear_par[static_cast(j)]; } - using Basis = std::array; - using BasisGrad = std::array, nnl>; }; +/** + * @brief Features of a peak on a background that may be higher on its left, + * as the start estimates of the peak models need them. + */ +struct PeakEstimate { + double mu; // x at the maximum of y + double y_max; // maximum of y + double left_mean; // mean of the first ~10% of y + double right_mean; // mean of the last ~10% of y + double sigma; // width of the peak +}; + +/** + * @brief Estimate the features of a peak. Assumes x is sorted ascending. + * + * sigma comes from the half width on the right side of the peak, above the + * right background, so that a plateau on the left does not bias it. It falls + * back to 10% of the x range. + */ +inline PeakEstimate estimate_peak(NDView x, NDView y) { + const ssize_t n = y.size(); + + const auto max_it = std::max_element(y.begin(), y.end()); + const ssize_t i_max = std::distance(y.begin(), max_it); + const double y_max = *max_it; + const double mu = x[i_max]; + + const double x_range = std::max(x[n - 1] - x[0], 1e-9); + + const ssize_t tail = std::min(std::max(n / 10, 2), n); + + double left_mean = 0.0; + double right_mean = 0.0; + + for (ssize_t i = 0; i < tail; ++i) + left_mean += y[i]; + left_mean /= tail; + + for (ssize_t i = n - tail; i < n; ++i) + right_mean += y[i]; + right_mean /= tail; + + const double half = right_mean + 0.5 * (y_max - right_mean); + + double sigma = 0.1 * x_range; + for (ssize_t i = i_max; i < n; ++i) { + if (y[i] <= half) { + const double half_width = std::abs(x[i] - mu); + if (half_width > 1e-12) + sigma = half_width / 1.1774100225154747; // sqrt(2*ln(2)) + break; + } + } + sigma = std::max(sigma, 1e-6); + + return {mu, y_max, left_mean, right_mean, sigma}; +} + +/** + * @brief A Gaussian peak G and the erfc step H below it at one point, with + * their derivatives with respect to the peak position and width: + * + * G = exp(-u^2 / 2), H = erfc(u / sqrt(2)) / 2, u = (x - mu) / sigma + * + * dG/dmu = G * u / sigma + * dG/dsigma = G * u^2 / sigma + * dH/dmu = G / (sqrt(2*pi) * sigma) + * dH/dsigma = G * u / (sqrt(2*pi) * sigma) + */ +struct PeakTerms { + double G; + double H; + double dG_dmu; + double dG_dsig; + double dH_dmu; + double dH_dsig; +}; + +/** + * @brief The terms of a peak from its Gaussian G = exp(-u^2 / 2) at + * u = (x - mu) / sigma, shared by the peak models. + * + * The caller evaluates the exponential: with std::exp at a single point, and + * with fast_exp in a loop that should vectorise. A model with two peaks + * evaluates both exponentials before it asks for the terms, which then stay + * in registers rather than being spilled around the second call of std::exp. + * The function is inlined, so the terms a caller does not use cost nothing. + */ +ALWAYS_INLINE PeakTerms peak_terms(double G, double u, double inv_sig) { + PeakTerms t; + t.G = G; + t.H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); + t.dG_dmu = G * u * inv_sig; + t.dG_dsig = t.dG_dmu * u; + t.dH_dmu = inv_sqrt_2pi * G * inv_sig; + t.dH_dsig = t.dH_dmu * u; + return t; +} + // _____________________________________________________________________ // // Pol1 @@ -231,12 +313,16 @@ struct Pol1 { /** @brief Both parameters are linear: f = p0 * 1 + p1 * x. */ static constexpr std::array linear_par = {0, 1}; - static void basis_and_grad( - double x, [[maybe_unused]] const std::vector &par, - std::array &phi, - [[maybe_unused]] std::array, 0> &dphi) { - phi[0] = 1.0; - phi[1] = x; + static void basis_columns(const double *__restrict x, ssize_t n, + [[maybe_unused]] const std::vector &par, + double *__restrict phi, + [[maybe_unused]] double *__restrict dphi) { + double *one = phi; + double *x_col = phi + n; + for (ssize_t i = 0; i < n; ++i) { + one[i] = 1.0; + x_col[i] = x[i]; + } } static bool is_valid([[maybe_unused]] const std::vector &par) { @@ -301,13 +387,18 @@ struct Pol2 { /** @brief All parameters are linear: f = p0 * 1 + p1 * x + p2 * x^2. */ static constexpr std::array linear_par = {0, 1, 2}; - static void basis_and_grad( - double x, [[maybe_unused]] const std::vector &par, - std::array &phi, - [[maybe_unused]] std::array, 0> &dphi) { - phi[0] = 1.0; - phi[1] = x; - phi[2] = x * x; + static void basis_columns(const double *__restrict x, ssize_t n, + [[maybe_unused]] const std::vector &par, + double *__restrict phi, + [[maybe_unused]] double *__restrict dphi) { + double *one = phi; + double *x_col = phi + n; + double *x2_col = phi + 2 * n; + for (ssize_t i = 0; i < n; ++i) { + one[i] = 1.0; + x_col[i] = x[i]; + x2_col[i] = x[i] * x[i]; + } } static bool is_valid([[maybe_unused]] const std::vector &par) { @@ -420,23 +511,6 @@ struct Gaussian { /** @brief A is linear: f = A * exp(-(x - mu)^2 / (2 sigma^2)). */ static constexpr std::array linear_par = {0}; - static void basis_and_grad(double x, const std::vector &par, - std::array &phi, - std::array, 2> &dphi) { - const double mu = par[1]; - const double sig = par[2]; - - const double dx = x - mu; - const double inv_sig = 1.0 / sig; - const double inv_sig2 = inv_sig * inv_sig; - const double e = std::exp(-0.5 * dx * dx * inv_sig2); - - phi[0] = e; - dphi[0][0] = e * dx * inv_sig2; // d/dmu - dphi[1][0] = e * dx * dx * inv_sig2 * inv_sig; // d/dsigma - } - - /** @brief basis_and_grad for all points, one column per function. */ static void basis_columns(const double *__restrict x, ssize_t n, const std::vector &par, double *__restrict phi, double *__restrict dphi) { @@ -517,16 +591,11 @@ struct Gaussian { * par[2] = mu (shared Gaussian center / step inflection point) * par[3] = sigma (shared Gaussian width / step width, must be > 0) * - * Analytic partial derivatives: - * Let dx = x - mu - * z = dx / (sqrt(2)*sigma) - * E = exp(-z^2) - * H = 0.5 * (1 - erf(z)) - * - * df/dA = E + * Analytic partial derivatives, with the peak G and the step H of PeakTerms: + * df/dA = G * df/dS = H - * df/dmu = A*E*dx/sigma^2 + S*E/(sqrt(2*pi)*sigma) - * df/dsigma = A*E*dx^2/sigma^3 + S*E*dx/(sqrt(2*pi)*sigma^2) + * df/dmu = A * dG/dmu + S * dH/dmu + * df/dsigma = A * dG/dsigma + S * dH/dsigma */ struct GaussianErfcPlateau { static constexpr std::size_t npar = 4; @@ -544,13 +613,10 @@ struct GaussianErfcPlateau { const double mu = par[2]; const double sig = par[3]; - const double dx = x - mu; - const double z = dx * inv_sqrt2 / sig; - - const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); - - return A * e + S * step; + const double inv_sig = 1.0 / sig; + const double u = (x - mu) * inv_sig; + const PeakTerms t = peak_terms(std::exp(-0.5 * u * u), u, inv_sig); + return A * t.G + S * t.H; } static void eval_and_grad(double x, const std::vector &par, @@ -560,81 +626,41 @@ struct GaussianErfcPlateau { const double mu = par[2]; const double sig = par[3]; - const double dx = x - mu; - const double z = dx * inv_sqrt2 / sig; - - const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); - - f = A * e + S * step; - const double inv_sig = 1.0 / sig; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; + const double u = (x - mu) * inv_sig; + const PeakTerms t = peak_terms(std::exp(-0.5 * u * u), u, inv_sig); - g[0] = e; // df/dA - g[1] = step; // df/dS + f = A * t.G + S * t.H; - g[2] = A * e * dx * inv_sig2 + S * inv_sqrt_2pi * e * inv_sig; // df/dmu - - g[3] = A * e * dx * dx * inv_sig3 + - S * inv_sqrt_2pi * e * dx * inv_sig2; // df/dsigma + g[0] = t.G; // df/dA + g[1] = t.H; // df/dS + g[2] = A * t.dG_dmu + S * t.dH_dmu; // df/dmu + g[3] = A * t.dG_dsig + S * t.dH_dsig; // df/dsigma } - /** @brief A and S are linear: f = A * e + S * step. */ + /** @brief A and S are linear: f = A * G + S * H. */ static constexpr std::array linear_par = {0, 1}; - static void basis_and_grad(double x, const std::vector &par, - std::array &phi, - std::array, 2> &dphi) { - const double mu = par[2]; - const double sig = par[3]; - - const double dx = x - mu; - const double inv_sig = 1.0 / sig; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; - const double z = dx * inv_sqrt2 * inv_sig; - - const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); - - phi[0] = e; - phi[1] = step; - dphi[0][0] = e * dx * inv_sig2; // de/dmu - dphi[0][1] = inv_sqrt_2pi * e * inv_sig; // dstep/dmu - dphi[1][0] = e * dx * dx * inv_sig3; // de/dsigma - dphi[1][1] = inv_sqrt_2pi * e * dx * inv_sig2; // dstep/dsigma - } - - /** @brief basis_and_grad for all points, one column per function. */ static void basis_columns(const double *__restrict x, ssize_t n, const std::vector &par, double *__restrict phi, double *__restrict dphi) { const double mu = par[2]; const double inv_sig = 1.0 / par[3]; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; - const double cz = inv_sqrt2 * inv_sig; - const double ce = inv_sqrt_2pi * inv_sig; - const double ce2 = inv_sqrt_2pi * inv_sig2; - double *e_col = phi; - double *step_col = phi + n; - double *de_dmu = dphi; - double *dstep_dmu = dphi + n; - double *de_dsig = dphi + 2 * n; - double *dstep_dsig = dphi + 3 * n; + double *G = phi; + double *H = phi + n; + double *dG_dmu = dphi; + double *dH_dmu = dphi + n; + double *dG_dsig = dphi + 2 * n; + double *dH_dsig = dphi + 3 * n; for (ssize_t i = 0; i < n; ++i) { - const double dx = x[i] - mu; - const double z = dx * cz; - const double e = fast_exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); - e_col[i] = e; - step_col[i] = step; - de_dmu[i] = e * dx * inv_sig2; - dstep_dmu[i] = ce * e; - de_dsig[i] = e * dx * dx * inv_sig3; - dstep_dsig[i] = ce2 * e * dx; + const double u = (x[i] - mu) * inv_sig; + const PeakTerms t = peak_terms(fast_exp(-0.5 * u * u), u, inv_sig); + G[i] = t.G; + H[i] = t.H; + dG_dmu[i] = t.dG_dmu; + dH_dmu[i] = t.dH_dmu; + dG_dsig[i] = t.dG_dsig; + dH_dsig[i] = t.dH_dsig; } } @@ -651,56 +677,15 @@ struct GaussianErfcPlateau { * S ~ average of first 10% of y values * mu ~ x at maximum y * A ~ y_max - 0.5*S - * sigma ~ right-side half-width of the peak, fallback to 10% x-range + * sigma ~ right-side half-width of the peak, see estimate_peak */ static std::array estimate_par(NDView x, NDView y) { - const ssize_t n = y.size(); - - const auto max_it = std::max_element(y.begin(), y.end()); - const ssize_t i_max = std::distance(y.begin(), max_it); - - const double y_max = *max_it; - - const double x_range = std::max(x[n - 1] - x[0], 1e-9); - - const ssize_t tail = std::min(std::max(n / 10, 2), n); - - double left_mean = 0.0; - double right_mean = 0.0; - - for (ssize_t i = 0; i < tail; ++i) - left_mean += y[i]; - left_mean /= tail; - - for (ssize_t i = n - tail; i < n; ++i) - right_mean += y[i]; - right_mean /= tail; - - const double S = left_mean; - const double mu = x[i_max]; - + const PeakEstimate peak = estimate_peak(x, y); + const double S = peak.left_mean; // At x = mu, the erfc step contributes roughly S/2. - const double A = y_max - 0.5 * S; - - // Estimate sigma from the right half-width of the peak. - // This avoids the left plateau biasing the FWHM estimate. - const double half = right_mean + 0.5 * (y_max - right_mean); - - double sig = 0.1 * x_range; - for (ssize_t i = i_max; i < n; ++i) { - if (y[i] <= half) { - const double half_width = std::abs(x[i] - mu); - if (half_width > 1e-12) { - sig = half_width / 1.1774100225154747; // sqrt(2*ln(2)) - } - break; - } - } - - sig = std::max(sig, 1e-6); - - return {A, S, mu, sig}; + const double A = peak.y_max - 0.5 * S; + return {A, S, peak.mu, peak.sigma}; } }; @@ -754,13 +739,10 @@ struct GaussianChargeSharing { const double N = par[4]; const double C = par[5]; - const double dx = x - mu; - const double u = dx / sig; - - const double G = std::exp(-0.5 * u * u); - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - - return p0 - p1 * x + N * (G + C * H); + const double inv_sig = 1.0 / sig; + const double u = (x - mu) * inv_sig; + const PeakTerms t = peak_terms(std::exp(-0.5 * u * u), u, inv_sig); + return p0 - p1 * x + N * (t.G + C * t.H); } static void eval_and_grad(double x, const std::vector &par, @@ -772,30 +754,18 @@ struct GaussianChargeSharing { const double N = par[4]; const double C = par[5]; - const double dx = x - mu; - const double u = dx / sig; - - const double G = std::exp(-0.5 * u * u); - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - - f = p0 - p1 * x + N * (G + C * H); - const double inv_sig = 1.0 / sig; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; + const double u = (x - mu) * inv_sig; + const PeakTerms t = peak_terms(std::exp(-0.5 * u * u), u, inv_sig); - const double dG_dmu = G * dx * inv_sig2; - const double dG_dsig = G * dx * dx * inv_sig3; - - const double dH_dmu = inv_sqrt_2pi * G * inv_sig; - const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2; + f = p0 - p1 * x + N * (t.G + C * t.H); g[0] = 1.0; g[1] = -x; - g[2] = N * (dG_dmu + C * dH_dmu); - g[3] = N * (dG_dsig + C * dH_dsig); - g[4] = G + C * H; - g[5] = N * H; + g[2] = N * (t.dG_dmu + C * t.dH_dmu); + g[3] = N * (t.dG_dsig + C * t.dH_dsig); + g[4] = t.G + C * t.H; + g[5] = N * t.H; } /** @@ -804,49 +774,12 @@ struct GaussianChargeSharing { */ static constexpr std::array linear_par = {0, 1, 4}; - static void basis_and_grad(double x, const std::vector &par, - std::array &phi, - std::array, 3> &dphi) { - const double mu = par[2]; - const double sig = par[3]; - const double C = par[5]; - - const double inv_sig = 1.0 / sig; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; - - const double dx = x - mu; - const double u = dx * inv_sig; - - const double G = std::exp(-0.5 * u * u); - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - - const double dG_dmu = G * dx * inv_sig2; - const double dG_dsig = G * dx * dx * inv_sig3; - const double dH_dmu = inv_sqrt_2pi * G * inv_sig; - const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2; - - phi[0] = 1.0; - phi[1] = -x; - phi[2] = G + C * H; - for (int k = 0; k < 3; ++k) { - dphi[k][0] = 0.0; - dphi[k][1] = 0.0; - } - dphi[0][2] = dG_dmu + C * dH_dmu; // d/dmu - dphi[1][2] = dG_dsig + C * dH_dsig; // d/dsigma - dphi[2][2] = H; // d/dC - } - - /** @brief basis_and_grad for all points, one column per function. */ static void basis_columns(const double *__restrict x, ssize_t n, const std::vector &par, double *__restrict phi, double *__restrict dphi) { const double mu = par[2]; const double inv_sig = 1.0 / par[3]; const double C = par[5]; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; double *one = phi; double *minus_x = phi + n; double *shape = phi + 2 * n; @@ -859,20 +792,14 @@ struct GaussianChargeSharing { double *d_C = dphi + 8 * n; for (ssize_t i = 0; i < n; ++i) { const double xi = x[i]; - const double dx = xi - mu; - const double u = dx * inv_sig; - const double G = fast_exp(-0.5 * u * u); - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - const double dG_dmu = G * dx * inv_sig2; - const double dG_dsig = G * dx * dx * inv_sig3; - const double dH_dmu = inv_sqrt_2pi * G * inv_sig; - const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2; + const double u = (xi - mu) * inv_sig; + const PeakTerms t = peak_terms(fast_exp(-0.5 * u * u), u, inv_sig); one[i] = 1.0; minus_x[i] = -xi; - shape[i] = G + C * H; - d_mu[i] = dG_dmu + C * dH_dmu; - d_sig[i] = dG_dsig + C * dH_dsig; - d_C[i] = H; + shape[i] = t.G + C * t.H; + d_mu[i] = t.dG_dmu + C * t.dH_dmu; + d_sig[i] = t.dG_dsig + C * t.dH_dsig; + d_C[i] = t.H; } } @@ -882,50 +809,13 @@ struct GaussianChargeSharing { static std::array estimate_par(NDView x, NDView y) { - const ssize_t n = y.size(); - - const auto max_it = std::max_element(y.begin(), y.end()); - const ssize_t i_max = std::distance(y.begin(), max_it); - - const double x_range = std::max(x[n - 1] - x[0], 1e-9); - - const ssize_t tail = std::min(std::max(n / 10, 2), n); - - double left_mean = 0.0; - double right_mean = 0.0; - - for (ssize_t i = 0; i < tail; ++i) - left_mean += y[i]; - left_mean /= tail; - - for (ssize_t i = n - tail; i < n; ++i) - right_mean += y[i]; - right_mean /= tail; - - const double p0 = right_mean; + const PeakEstimate peak = estimate_peak(x, y); + const double p0 = peak.right_mean; const double p1 = 0.0; - const double mu = x[i_max]; - - const double N = std::max(*max_it - right_mean, 1e-9); - + const double N = std::max(peak.y_max - peak.right_mean, 1e-9); // Left plateau excess is roughly N*C. - const double C = (left_mean - right_mean) / N; - - double sigma = 0.1 * x_range; - const double half = right_mean + 0.5 * ((*max_it) - right_mean); - - for (ssize_t i = i_max; i < n; ++i) { - if (y[i] <= half) { - const double half_width = std::abs(x[i] - mu); - if (half_width > 1e-12) - sigma = half_width / 1.1774100225154747; // sqrt(2*ln(2)) - break; - } - } - - sigma = std::max(sigma, 1e-6); - - return {p0, p1, mu, sigma, N, C}; + const double C = (peak.left_mean - peak.right_mean) / N; + return {p0, p1, peak.mu, peak.sigma, N, C}; } }; @@ -987,22 +877,16 @@ struct GaussianChargeSharingKb { const double r = par[6]; const double q = par[7]; - const double mu_b = r * mu; - - const double dx = x - mu; - const double dxb = x - mu_b; - - const double u = dx / sig; - const double ub = dxb / sig; - + const double inv_sig = 1.0 / sig; + const double u = (x - mu) * inv_sig; + const double ub = (x - r * mu) * inv_sig; const double G = std::exp(-0.5 * u * u); const double Gb = std::exp(-0.5 * ub * ub); + const PeakTerms a = peak_terms(G, u, inv_sig); + const PeakTerms b = peak_terms(Gb, ub, inv_sig); - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - const double Hb = 0.5 * (1.0 - fast_erf_from_exp(ub * inv_sqrt2, Gb)); - - const double ka = G + C * H; - const double kb = Gb + C * Hb; + const double ka = a.G + C * a.H; + const double kb = b.G + C * b.H; return p0 - p1 * x + N * (ka + q * kb); } @@ -1018,44 +902,24 @@ struct GaussianChargeSharingKb { const double r = par[6]; const double q = par[7]; - const double mu_b = r * mu; - - const double dx = x - mu; - const double dxb = x - mu_b; - - const double u = dx / sig; - const double ub = dxb / sig; - + const double inv_sig = 1.0 / sig; + const double u = (x - mu) * inv_sig; + const double ub = (x - r * mu) * inv_sig; const double G = std::exp(-0.5 * u * u); const double Gb = std::exp(-0.5 * ub * ub); + const PeakTerms a = peak_terms(G, u, inv_sig); + const PeakTerms b = peak_terms(Gb, ub, inv_sig); - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - const double Hb = 0.5 * (1.0 - fast_erf_from_exp(ub * inv_sqrt2, Gb)); - - const double ka = G + C * H; - const double kb = Gb + C * Hb; + const double ka = a.G + C * a.H; + const double kb = b.G + C * b.H; f = p0 - p1 * x + N * (ka + q * kb); - const double inv_sig = 1.0 / sig; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; + const double dka_dmu = a.dG_dmu + C * a.dH_dmu; + const double dka_dsig = a.dG_dsig + C * a.dH_dsig; - const double dG_dmu = G * dx * inv_sig2; - const double dG_dsig = G * dx * dx * inv_sig3; - const double dH_dmu = inv_sqrt_2pi * G * inv_sig; - const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2; - - const double dGb_dmu_b = Gb * dxb * inv_sig2; - const double dGb_dsig = Gb * dxb * dxb * inv_sig3; - const double dHb_dmu_b = inv_sqrt_2pi * Gb * inv_sig; - const double dHb_dsig = inv_sqrt_2pi * Gb * dxb * inv_sig2; - - const double dka_dmu = dG_dmu + C * dH_dmu; - const double dka_dsig = dG_dsig + C * dH_dsig; - - const double dkb_dmu_b = dGb_dmu_b + C * dHb_dmu_b; - const double dkb_dsig = dGb_dsig + C * dHb_dsig; + const double dkb_dmu_b = b.dG_dmu + C * b.dH_dmu; + const double dkb_dsig = b.dG_dsig + C * b.dH_dsig; g[0] = 1.0; g[1] = -x; @@ -1067,7 +931,7 @@ struct GaussianChargeSharingKb { g[4] = ka + q * kb; - g[5] = N * (H + q * Hb); + g[5] = N * (a.H + q * b.H); // mu_b = r * mu, so dmu_b/dr = mu g[6] = N * q * mu * dkb_dmu_b; @@ -1081,67 +945,6 @@ struct GaussianChargeSharingKb { */ static constexpr std::array linear_par = {0, 1, 4}; - static void basis_and_grad(double x, const std::vector &par, - std::array &phi, - std::array, 5> &dphi) { - const double mu = par[2]; - const double sig = par[3]; - const double C = par[5]; - const double r = par[6]; - const double q = par[7]; - - const double inv_sig = 1.0 / sig; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; - - const double mu_b = r * mu; - - const double dx = x - mu; - const double dxb = x - mu_b; - - const double u = dx * inv_sig; - const double ub = dxb * inv_sig; - - const double G = std::exp(-0.5 * u * u); - const double Gb = std::exp(-0.5 * ub * ub); - - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - const double Hb = 0.5 * (1.0 - fast_erf_from_exp(ub * inv_sqrt2, Gb)); - - const double ka = G + C * H; - const double kb = Gb + C * Hb; - - const double dG_dmu = G * dx * inv_sig2; - const double dG_dsig = G * dx * dx * inv_sig3; - const double dH_dmu = inv_sqrt_2pi * G * inv_sig; - const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2; - - const double dGb_dmu_b = Gb * dxb * inv_sig2; - const double dGb_dsig = Gb * dxb * dxb * inv_sig3; - const double dHb_dmu_b = inv_sqrt_2pi * Gb * inv_sig; - const double dHb_dsig = inv_sqrt_2pi * Gb * dxb * inv_sig2; - - const double dka_dmu = dG_dmu + C * dH_dmu; - const double dka_dsig = dG_dsig + C * dH_dsig; - - const double dkb_dmu_b = dGb_dmu_b + C * dHb_dmu_b; - const double dkb_dsig = dGb_dsig + C * dHb_dsig; - - phi[0] = 1.0; - phi[1] = -x; - phi[2] = ka + q * kb; - for (int k = 0; k < 5; ++k) { - dphi[k][0] = 0.0; - dphi[k][1] = 0.0; - } - dphi[0][2] = dka_dmu + q * r * dkb_dmu_b; // d/dmu, mu_b = r * mu - dphi[1][2] = dka_dsig + q * dkb_dsig; // d/dsigma - dphi[2][2] = H + q * Hb; // d/dC - dphi[3][2] = q * mu * dkb_dmu_b; // d/dr - dphi[4][2] = kb; // d/dq - } - - /** @brief basis_and_grad for all points, one column per function. */ static void basis_columns(const double *__restrict x, ssize_t n, const std::vector &par, double *__restrict phi, double *__restrict dphi) { @@ -1150,8 +953,6 @@ struct GaussianChargeSharingKb { const double C = par[5]; const double r = par[6]; const double q = par[7]; - const double inv_sig2 = inv_sig * inv_sig; - const double inv_sig3 = inv_sig2 * inv_sig; const double mu_b = r * mu; double *one = phi; double *minus_x = phi + n; @@ -1168,35 +969,23 @@ struct GaussianChargeSharingKb { double *d_q = dphi + 14 * n; for (ssize_t i = 0; i < n; ++i) { const double xi = x[i]; - const double dx = xi - mu; - const double dxb = xi - mu_b; - const double u = dx * inv_sig; - const double ub = dxb * inv_sig; - const double G = fast_exp(-0.5 * u * u); - const double Gb = fast_exp(-0.5 * ub * ub); - const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); - const double Hb = - 0.5 * (1.0 - fast_erf_from_exp(ub * inv_sqrt2, Gb)); - const double ka = G + C * H; - const double kb = Gb + C * Hb; - const double dG_dmu = G * dx * inv_sig2; - const double dG_dsig = G * dx * dx * inv_sig3; - const double dH_dmu = inv_sqrt_2pi * G * inv_sig; - const double dH_dsig = inv_sqrt_2pi * G * dx * inv_sig2; - const double dGb_dmu_b = Gb * dxb * inv_sig2; - const double dGb_dsig = Gb * dxb * dxb * inv_sig3; - const double dHb_dmu_b = inv_sqrt_2pi * Gb * inv_sig; - const double dHb_dsig = inv_sqrt_2pi * Gb * dxb * inv_sig2; - const double dka_dmu = dG_dmu + C * dH_dmu; - const double dka_dsig = dG_dsig + C * dH_dsig; - const double dkb_dmu_b = dGb_dmu_b + C * dHb_dmu_b; - const double dkb_dsig = dGb_dsig + C * dHb_dsig; + const double u = (xi - mu) * inv_sig; + const double ub = (xi - mu_b) * inv_sig; + const PeakTerms a = peak_terms(fast_exp(-0.5 * u * u), u, inv_sig); + const PeakTerms b = + peak_terms(fast_exp(-0.5 * ub * ub), ub, inv_sig); + const double ka = a.G + C * a.H; + const double kb = b.G + C * b.H; + const double dka_dmu = a.dG_dmu + C * a.dH_dmu; + const double dka_dsig = a.dG_dsig + C * a.dH_dsig; + const double dkb_dmu_b = b.dG_dmu + C * b.dH_dmu; + const double dkb_dsig = b.dG_dsig + C * b.dH_dsig; one[i] = 1.0; minus_x[i] = -xi; shape[i] = ka + q * kb; d_mu[i] = dka_dmu + q * r * dkb_dmu_b; d_sig[i] = dka_dsig + q * dkb_dsig; - d_C[i] = H + q * Hb; + d_C[i] = a.H + q * b.H; d_r[i] = q * mu * dkb_dmu_b; d_q[i] = kb; } @@ -1227,15 +1016,15 @@ struct GaussianChargeSharingKb { // _____________________________________________________________________ // -// RisingScurve +// RisingScurve and FallingScurve // _____________________________________________________________________ /** - * @brief Rising S-curve (error-function step) with linear baseline and - * post-step slope. + * @brief S-curve (error-function step) with linear baseline and post-step + * slope. Direction = +1 rises with x, Direction = -1 falls. * - * f(x) = (p0 + p1*x) + 0.5*(1 + erf((x - mu) / (sqrt(2)*sigma))) * (A + C*(x - - * mu))) + * f(x) = (p0 + p1*x) + * + 0.5*(1 ± erf((x - mu) / (sqrt(2)*sigma))) * (A + C*(x - mu)) * * Parameters: * par[0] = p0 (baseline offset) @@ -1244,8 +1033,15 @@ struct GaussianChargeSharingKb { * par[3] = sigma (transition width, must be > 0) * par[4] = A (step amplitude) * par[5] = C (post-step slope) + * + * With the step S(x) = 0.5*(1 ± erf(z)), z = (x - mu) / (sqrt(2)*sigma): + * dS/dmu = -(±1/sqrt(2*pi)) * exp(-z^2) / sigma + * dS/dsigma = -(±1/sqrt(2*pi)) * exp(-z^2) * (x - mu) / sigma^2 */ -struct RisingScurve { +template struct Scurve { + static_assert(Direction == 1 || Direction == -1, + "an S-curve either rises (+1) or falls (-1)"); + static constexpr std::size_t npar = 6; static constexpr std::array param_info = {{ @@ -1266,18 +1062,12 @@ struct RisingScurve { const double p5 = par[5]; const double dx = x - p2; - const double step = 0.5 * (1.0 + fast_erf(dx * inv_sqrt2 / p3)); + const double step = 0.5 * (1.0 + sign * fast_erf(dx * inv_sqrt2 / p3)); return (p0 + p1 * x) + step * (p4 + p5 * dx); } - /** - * @brief Evaluate function value and partial derivatives in a single pass. - * - * Uses: - * S(x) = 0.5*(1 + erf(z)), z = (x-p2) / (sqrt(2)*p3) - * dS/dp2 = -(1/sqrt(2*pi)) * exp(-z^2) / p3 - * dS/dp3 = -(1/sqrt(2*pi)) * exp(-z^2) * (x-p2) / p3^2 - */ + /** @brief Evaluate function value and partial derivatives in a single + * pass. */ static void eval_and_grad(double x, const std::vector &par, double &f, std::array &g) { const double p0 = par[0]; @@ -1290,13 +1080,13 @@ struct RisingScurve { const double dx = x - p2; const double z = dx * inv_sqrt2 / p3; const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 + fast_erf_from_exp(z, e)); + const double step = 0.5 * (1.0 + sign * fast_erf_from_exp(z, e)); const double amp = p4 + p5 * dx; f = (p0 + p1 * x) + step * amp; - const double dSdp2 = -inv_sqrt_2pi * e / p3; - const double dSdp3 = -inv_sqrt_2pi * e * dx / (p3 * p3); + const double dSdp2 = dstep * e / p3; + const double dSdp3 = dstep * e * dx / (p3 * p3); g[0] = 1.0; // df/dp0 g[1] = x; // df/dp1 @@ -1312,43 +1102,13 @@ struct RisingScurve { */ static constexpr std::array linear_par = {0, 1, 4, 5}; - static void basis_and_grad(double x, const std::vector &par, - std::array &phi, - std::array, 2> &dphi) { - const double p2 = par[2]; - const double p3 = par[3]; - - const double inv_p3 = 1.0 / p3; - const double dx = x - p2; - const double z = dx * inv_sqrt2 * inv_p3; - const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 + fast_erf_from_exp(z, e)); - - const double dSdp2 = -inv_sqrt_2pi * e * inv_p3; - const double dSdp3 = dSdp2 * dx * inv_p3; - - phi[0] = 1.0; - phi[1] = x; - phi[2] = step; - phi[3] = step * dx; - dphi[0][0] = 0.0; - dphi[0][1] = 0.0; - dphi[0][2] = dSdp2; // dS/dmu - dphi[0][3] = dSdp2 * dx - step; // d(S dx)/dmu - dphi[1][0] = 0.0; - dphi[1][1] = 0.0; - dphi[1][2] = dSdp3; // dS/dsigma - dphi[1][3] = dSdp3 * dx; // d(S dx)/dsigma - } - - /** @brief basis_and_grad for all points, one column per function. */ static void basis_columns(const double *__restrict x, ssize_t n, const std::vector &par, double *__restrict phi, double *__restrict dphi) { const double p2 = par[2]; const double inv_p3 = 1.0 / par[3]; const double cz = inv_sqrt2 * inv_p3; - const double ce = -inv_sqrt_2pi * inv_p3; + const double ce = dstep * inv_p3; double *one = phi; double *x_col = phi + n; double *step_col = phi + 2 * n; @@ -1365,7 +1125,7 @@ struct RisingScurve { const double dx = xi - p2; const double z = dx * cz; const double e = fast_exp(-z * z); - const double step = 0.5 * (1.0 + fast_erf_from_exp(z, e)); + const double step = 0.5 * (1.0 + sign * fast_erf_from_exp(z, e)); const double dSdp2 = ce * e; const double dSdp3 = dSdp2 * dx * inv_p3; one[i] = 1.0; @@ -1384,286 +1144,95 @@ struct RisingScurve { return par[3] != 0.0; } - /** @brief Data-driven initial parameter estimates for a rising S-curve. */ + /** + * @brief Data-driven initial parameter estimates. + * + * The baseline is a line through the ~10% of points before the step + * turns on (the first points of a rising curve, the last of a falling + * one), the plateau the mean of the ~10% at the other end. + */ static std::array estimate_par(NDView x, NDView y) { + constexpr bool rising = Direction > 0; const ssize_t n = y.size(); + const ssize_t n_tail = std::max(n / 10, 2); + const ssize_t base_begin = rising ? 0 : n - n_tail; + const ssize_t plateau_begin = rising ? n - n_tail : 0; - // baseline: average of first ~10% of points (before turn-on) - ssize_t n_base = std::max(n / 10, 2); double sum_y = 0, sum_xy = 0, sum_x = 0, sum_x2 = 0; - for (ssize_t i = 0; i < n_base; ++i) { + for (ssize_t i = base_begin; i < base_begin + n_tail; ++i) { sum_y += y[i]; sum_x += x[i]; sum_xy += x[i] * y[i]; sum_x2 += x[i] * x[i]; } - double denom = n_base * sum_x2 - sum_x * sum_x; + double denom = n_tail * sum_x2 - sum_x * sum_x; double p1 = (std::abs(denom) > 1e-30) - ? (n_base * sum_xy - sum_x * sum_y) / denom + ? (n_tail * sum_xy - sum_x * sum_y) / denom : 0.0; - double p0 = (sum_y - p1 * sum_x) / n_base; + double p0 = (sum_y - p1 * sum_x) / n_tail; - // plateau: average of last ~10% double plateau = 0; - ssize_t n_plat = std::max(n / 10, 2); - for (ssize_t i = n - n_plat; i < n; ++i) + for (ssize_t i = plateau_begin; i < plateau_begin + n_tail; ++i) plateau += y[i]; - plateau /= n_plat; + plateau /= n_tail; // amplitude: plateau minus baseline at midpoint double x_mid = 0.5 * (x[0] + x[n - 1]); double baseline_at_mid = p0 + p1 * x_mid; double p4 = plateau - baseline_at_mid; - // threshold: x where y first crosses 50% between baseline and plateau - double y_half = baseline_at_mid + 0.5 * p4; - double p2 = x_mid; // fallback - for (ssize_t i = 0; i < n; ++i) { - if (y[i] >= y_half) { - p2 = x[i]; - break; + // x of the first point past the given fraction of the step + const auto crossing = [&](double fraction, double fallback) { + const double level = baseline_at_mid + fraction * p4; + for (ssize_t i = 0; i < n; ++i) { + if (rising ? y[i] >= level : y[i] <= level) + return x[i]; } - } + return fallback; + }; - // sigma: estimate from transition width (10%-90% rise) - double y_10 = baseline_at_mid + 0.1 * p4; - double y_90 = baseline_at_mid + 0.9 * p4; - double x_10 = x[0], x_90 = x[n - 1]; - for (ssize_t i = 0; i < n; ++i) { - if (y[i] >= y_10) { - x_10 = x[i]; - break; - } - } - for (ssize_t i = 0; i < n; ++i) { - if (y[i] >= y_90) { - x_90 = x[i]; - break; - } - } - // for a Gaussian CDF: 10%-90% width = 2 * 1.2816 * sigma - double p3 = std::max((x_90 - x_10) / 2.5631, 1.0); + // threshold: halfway between baseline and plateau + double p2 = crossing(0.5, x_mid); + + // sigma: from the transition width between 10% and 90% of the step; + // for a Gaussian CDF that width is 2 * 1.2816 * sigma + double x_enter = crossing(rising ? 0.1 : 0.9, x[0]); + double x_leave = crossing(rising ? 0.9 : 0.1, x[n - 1]); + double p3 = std::max((x_leave - x_enter) / 2.5631, 1.0); double p5 = 0.0; // assume flat gain, let optimizer find the slope return {p0, p1, p2, p3, p4, p5}; } + + private: + static constexpr double sign = Direction; + // dS/dmu = dstep * exp(-z^2) / sigma + static constexpr double dstep = -sign * inv_sqrt_2pi; }; -// _____________________________________________________________________ -// -// FallingScurve -// _____________________________________________________________________ +/** @brief S-curve that rises with x, see Scurve. */ +struct RisingScurve : Scurve<+1> {}; -/** - * @brief Falling S-curve (complementary error-function step) with linear - * baseline and post-step slope. - * - * f(x) = (p0 + p1*x) + 0.5*(1 - erf((x - mu) / (sqrt(2)*sigma))) * (A + C*(x - - * mu)) - * - * Parameters are identical to RisingScurve. The only difference is the - * sign of the erf term, which flips the step direction (and the signs of - * dS/dp2 and dS/dp3 in the gradient). - */ -struct FallingScurve { - static constexpr std::size_t npar = 6; - - static constexpr std::array param_info = {{ - {"p0", -no_bound, no_bound}, - {"p1", -no_bound, no_bound}, - {"mu", -no_bound, no_bound}, - {"sigma", 1e-12, no_bound}, - {"A", -no_bound, no_bound}, - {"C", -no_bound, no_bound}, - }}; - - static double eval(double x, const std::vector &par) { - const double p0 = par[0]; - const double p1 = par[1]; - const double p2 = par[2]; - const double p3 = par[3]; - const double p4 = par[4]; - const double p5 = par[5]; - - const double dx = x - p2; - const double step = 0.5 * (1.0 - fast_erf(dx * inv_sqrt2 / p3)); - return (p0 + p1 * x) + step * (p4 + p5 * dx); - } - - static void eval_and_grad(double x, const std::vector &par, - double &f, std::array &g) { - const double p0 = par[0]; - const double p1 = par[1]; - const double p2 = par[2]; - const double p3 = par[3]; - const double p4 = par[4]; - const double p5 = par[5]; - - const double dx = x - p2; - const double z = dx * inv_sqrt2 / p3; - const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); - const double amp = p4 + p5 * dx; - - f = (p0 + p1 * x) + step * amp; - - const double dSdp2 = +inv_sqrt_2pi * e / p3; // sign flipped vs rising - const double dSdp3 = +inv_sqrt_2pi * e * dx / (p3 * p3); - - g[0] = 1.0; - g[1] = x; - g[2] = dSdp2 * amp - step * p5; - g[3] = dSdp3 * amp; - g[4] = step; - g[5] = step * dx; - } - - /** - * @brief p0, p1, A and C are linear: - * f = p0 * 1 + p1 * x + A * S + C * S * (x - mu). - */ - static constexpr std::array linear_par = {0, 1, 4, 5}; - - static void basis_and_grad(double x, const std::vector &par, - std::array &phi, - std::array, 2> &dphi) { - const double p2 = par[2]; - const double p3 = par[3]; - - const double inv_p3 = 1.0 / p3; - const double dx = x - p2; - const double z = dx * inv_sqrt2 * inv_p3; - const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); - - const double dSdp2 = - inv_sqrt_2pi * e * inv_p3; // sign flipped vs rising - const double dSdp3 = dSdp2 * dx * inv_p3; - - phi[0] = 1.0; - phi[1] = x; - phi[2] = step; - phi[3] = step * dx; - dphi[0][0] = 0.0; - dphi[0][1] = 0.0; - dphi[0][2] = dSdp2; // dS/dmu - dphi[0][3] = dSdp2 * dx - step; // d(S dx)/dmu - dphi[1][0] = 0.0; - dphi[1][1] = 0.0; - dphi[1][2] = dSdp3; // dS/dsigma - dphi[1][3] = dSdp3 * dx; // d(S dx)/dsigma - } - - /** @brief basis_and_grad for all points, one column per function. */ - static void basis_columns(const double *__restrict x, ssize_t n, - const std::vector &par, - double *__restrict phi, double *__restrict dphi) { - const double p2 = par[2]; - const double inv_p3 = 1.0 / par[3]; - const double cz = inv_sqrt2 * inv_p3; - const double ce = inv_sqrt_2pi * inv_p3; // sign flipped vs rising - double *one = phi; - double *x_col = phi + n; - double *step_col = phi + 2 * n; - double *step_dx = phi + 3 * n; - // d phi_0 and d phi_1 vanish for mu and sigma. - std::fill(dphi, dphi + 2 * n, 0.0); - std::fill(dphi + 4 * n, dphi + 6 * n, 0.0); - double *dmu_step = dphi + 2 * n; - double *dmu_step_dx = dphi + 3 * n; - double *dsig_step = dphi + 6 * n; - double *dsig_step_dx = dphi + 7 * n; - for (ssize_t i = 0; i < n; ++i) { - const double xi = x[i]; - const double dx = xi - p2; - const double z = dx * cz; - const double e = fast_exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); - const double dSdp2 = ce * e; - const double dSdp3 = dSdp2 * dx * inv_p3; - one[i] = 1.0; - x_col[i] = xi; - step_col[i] = step; - step_dx[i] = step * dx; - dmu_step[i] = dSdp2; - dmu_step_dx[i] = dSdp2 * dx - step; - dsig_step[i] = dSdp3; - dsig_step_dx[i] = dSdp3 * dx; - } - } - - /** @brief Reject degenerate width (zero transition width). */ - static bool is_valid(const std::vector &par) { - return par[3] != 0.0; - } - - /** @brief Data-driven initial parameter estimates for a falling S-curve. */ - static std::array estimate_par(NDView x, - NDView y) { - const ssize_t n = y.size(); - - // baseline: last ~10% of points (after turn-off) - ssize_t n_base = std::max(n / 10, 2); - double sum_y = 0, sum_xy = 0, sum_x = 0, sum_x2 = 0; - for (ssize_t i = n - n_base; i < n; ++i) { - sum_y += y[i]; - sum_x += x[i]; - sum_xy += x[i] * y[i]; - sum_x2 += x[i] * x[i]; - } - double denom = n_base * sum_x2 - sum_x * sum_x; - double p1 = (std::abs(denom) > 1e-30) - ? (n_base * sum_xy - sum_x * sum_y) / denom - : 0.0; - double p0 = (sum_y - p1 * sum_x) / n_base; - - // plateau: average of first ~10% - double plateau = 0; - ssize_t n_plat = std::max(n / 10, 2); - for (ssize_t i = 0; i < n_plat; ++i) - plateau += y[i]; - plateau /= n_plat; - - // amplitude: plateau minus baseline at midpoint - double x_mid = 0.5 * (x[0] + x[n - 1]); - double baseline_at_mid = p0 + p1 * x_mid; - double p4 = plateau - baseline_at_mid; - - // threshold: x where y first drops below 50% - double y_half = baseline_at_mid + 0.5 * p4; - double p2 = x_mid; // fallback - for (ssize_t i = 0; i < n; ++i) { - if (y[i] <= y_half) { - p2 = x[i]; - break; - } - } - - // sigma: estimate from transition width (90%-10% fall) - double y_90 = baseline_at_mid + 0.9 * p4; - double y_10 = baseline_at_mid + 0.1 * p4; - double x_90 = x[0], x_10 = x[n - 1]; - for (ssize_t i = 0; i < n; ++i) { - if (y[i] <= y_90) { - x_90 = x[i]; - break; - } - } - for (ssize_t i = 0; i < n; ++i) { - if (y[i] <= y_10) { - x_10 = x[i]; - break; - } - } - // same CDF relationship: 10%-90% width = 2 * 1.2816 * sigma - double p3 = std::max((x_10 - x_90) / 2.5631, 1.0); - - double p5 = 0.0; - - return {p0, p1, p2, p3, p4, p5}; - } -}; +/** @brief S-curve that falls with x, see Scurve. */ +struct FallingScurve : Scurve<-1> {}; } // namespace aare::model + +/** + * @brief The fit models: expands to X(Model) for every model in aare::model. + * + * The explicit instantiations of FitModel and the fit functions and the + * Python bindings are generated from this list, so a new model is added + * here once. It also needs its Minuit step scales in src/fit/MinuitSteps.hpp. + */ +#define AARE_FOR_EACH_FIT_MODEL(X) \ + X(Pol1) \ + X(Pol2) \ + X(Gaussian) \ + X(GaussianErfcPlateau) \ + X(GaussianChargeSharing) \ + X(GaussianChargeSharingKb) \ + X(RisingScurve) \ + X(FallingScurve) diff --git a/python/src/fit.hpp b/python/src/fit.hpp index f26cc20c..0e11b2c1 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -1,6 +1,8 @@ // SPDX-License-Identifier: MPL-2.0 +#include #include #include +#include #include #include #include @@ -12,12 +14,13 @@ namespace py = pybind11; using namespace pybind11::literals; +// Input arrays are converted, so that the fit can be called with any dtype. +using DoubleArray = + py::array_t; + template -py::object -fit_dispatch(const aare::FitModel &model, - py::array_t x, - py::array_t y, - py::object y_err_obj, int n_threads); +py::dict fit_dispatch(const aare::FitModel &model, DoubleArray x, + DoubleArray y, py::object y_err_obj, int n_threads); template void bind_fit_model(py::module &m, const char *name) { using FM = aare::FitModel; @@ -92,12 +95,8 @@ template void bind_fit_model(py::module &m, const char *name) { &FM::SetComputeErrors) .def_property("minimizer", &FM::minimizer, &FM::SetMinimizer) .def( - "__call__", // conversion ok, we want to be able to call with any - // dtype - [](const FM & /*self*/, - py::array_t x, - py::array_t - par) { + "__call__", + [](const FM & /*self*/, DoubleArray x, DoubleArray par) { auto x_view = make_view_1d(x); auto p_view = make_view_1d(par); @@ -111,15 +110,13 @@ template void bind_fit_model(py::module &m, const char *name) { }, py::arg("x"), py::arg("par")) .def( - "fit", // conversion ok - [](const FM &self, - py::array_t x, - py::array_t y, - py::object y_err_obj, int n_threads) -> py::object { + "fit", + [](const FM &self, DoubleArray x, DoubleArray y, + py::object y_err_obj, int n_threads) { return fit_dispatch(self, x, y, y_err_obj, n_threads); }, R"doc( - Fit this model to 1D or 3D data using Minuit2. + Fit this model to 1D or 3D data. The minimizer is selected by the ``minimizer`` property (``Minimizer.Migrad`` by default, ``Minimizer.Fumili``, @@ -141,163 +138,75 @@ template void bind_fit_model(py::module &m, const char *name) { py::arg("n_threads") = 4); } +// Fit one model to 1D or 3D data and pack the result: "par", "par_err" when +// the model computes errors, and "chi2". template -py::dict pack_1d_result_dict(const aare::NDArray &result, - bool compute_errors) { - constexpr std::size_t npar = Model::npar; - - auto res = result.view(); - - auto par_out = new NDArray({npar}, 0.0); - auto chi2_out = new NDArray({1}, 0.0); - - auto par_view = par_out->view(); - auto chi2_view = chi2_out->view(); - - for (std::size_t i = 0; i < npar; ++i) { - par_view(i) = res(i); - } - - if (compute_errors) { - auto err_out = new NDArray({npar}, 0.0); - auto err_view = err_out->view(); - - for (std::size_t i = 0; i < npar; ++i) { - err_view(i) = res(npar + i); - } - - chi2_view(0) = res(2 * npar); - - return py::dict("par"_a = return_image_data(par_out), - "par_err"_a = return_image_data(err_out), - "chi2"_a = return_image_data(chi2_out)); - } else { - chi2_view(0) = res(npar); - - return py::dict("par"_a = return_image_data(par_out), - "chi2"_a = return_image_data(chi2_out)); - } -} - -// Helper: typed dispatch for one Model, handles 1D/3D + y_err logic -template -py::object -fit_dispatch(const aare::FitModel &model, - py::array_t x, - py::array_t y, - py::object y_err_obj, int n_threads) { - constexpr std::size_t npar = Model::npar; +py::dict fit_dispatch(const aare::FitModel &model, DoubleArray x, + DoubleArray y, py::object y_err_obj, int n_threads) { + constexpr auto npar = static_cast(Model::npar); + const bool want_errors = model.compute_errors(); + const bool weighted = !y_err_obj.is_none(); + DoubleArray y_err; + if (weighted) + y_err = py::cast(y_err_obj); + auto x_view = make_view_1d(x); + py::dict result; if (y.ndim() == 3) { - auto par_out = - new NDArray({y.shape(0), y.shape(1), npar}, 0.0); - auto chi2_out = new NDArray({y.shape(0), y.shape(1)}, 0.0); + const ssize_t rows = y.shape(0); + const ssize_t cols = y.shape(1); + using Cube = NDArray; + using Image = NDArray; + const std::array par_shape{rows, cols, npar}; + auto par = std::make_unique(par_shape, 0.0); + auto par_err = + want_errors ? std::make_unique(par_shape, 0.0) : nullptr; + auto chi2 = + std::make_unique(std::array{rows, cols}, 0.0); - auto x_view = make_view_1d(x); - auto y_view = make_view_3d(y); + aare::fit_3d( + model, x_view, make_view_3d(y), + weighted ? make_view_3d(y_err) : NDView{}, par->view(), + par_err ? par_err->view() : NDView{}, chi2->view(), + n_threads); - if (!y_err_obj.is_none()) { - auto y_err = py::cast< - py::array_t>( - y_err_obj); - - auto err_out = - new NDArray({y.shape(0), y.shape(1), npar}, 0.0); - auto y_view_err = make_view_3d(y_err); - - aare::fit_3d(model, x_view, y_view, y_view_err, - par_out->view(), err_out->view(), - chi2_out->view(), n_threads); - - if (model.compute_errors()) { - return py::dict("par"_a = return_image_data(par_out), - "par_err"_a = return_image_data(err_out), - "chi2"_a = return_image_data(chi2_out)); - } else { - delete err_out; - return py::dict("par"_a = return_image_data(par_out), - "chi2"_a = return_image_data(chi2_out)); - } - } else { - // Unweighted fit; parameter errors are still available when - // requested, as in the 1D path. - NDView dummy_err{}; - - if (model.compute_errors()) { - auto err_out = - new NDArray({y.shape(0), y.shape(1), npar}, 0.0); - - aare::fit_3d(model, x_view, y_view, dummy_err, - par_out->view(), err_out->view(), - chi2_out->view(), n_threads); - - return py::dict("par"_a = return_image_data(par_out), - "par_err"_a = return_image_data(err_out), - "chi2"_a = return_image_data(chi2_out)); - } - - NDView dummy_err_out{}; - - aare::fit_3d(model, x_view, y_view, dummy_err, - par_out->view(), dummy_err_out, - chi2_out->view(), n_threads); - - return py::dict("par"_a = return_image_data(par_out), - "chi2"_a = return_image_data(chi2_out)); - } + result["par"] = return_image_data(par.release()); + if (par_err) + result["par_err"] = return_image_data(par_err.release()); + result["chi2"] = return_image_data(chi2.release()); } else if (y.ndim() == 1) { - NDArray result{}; - - auto x_view = make_view_1d(x); - auto y_view = make_view_1d(y); - - if (!y_err_obj.is_none()) { - auto y_err = py::cast< - py::array_t>( - y_err_obj); - - auto y_view_err = make_view_1d(y_err); - result = aare::fit_pixel(model, x_view, y_view, y_view_err); - } else { - result = aare::fit_pixel(model, x_view, y_view); - } - - return pack_1d_result_dict(result, model.compute_errors()); + // [par..., (err...,) chi2], see fit_pixel + const auto flat = aare::fit_pixel( + model, x_view, make_view_1d(y), + weighted ? make_view_1d(y_err) : NDView{}); + result["par"] = py::array_t(npar, flat.data()); + if (want_errors) + result["par_err"] = py::array_t(npar, flat.data() + npar); + result["chi2"] = py::array_t(1, flat.data() + flat.size() - 1); } else { throw std::runtime_error("Data must be 1D or 3D."); } + return result; } // Resolve the model type behind a Python object and run the fit. -py::object dispatch_any_model( - py::object model_obj, - py::array_t x, - py::array_t y, - py::object y_err_obj, int n_threads) { - using namespace aare::model; - +py::dict dispatch_any_model(py::object model_obj, DoubleArray x, DoubleArray y, + py::object y_err_obj, int n_threads) { #define AARE_DISPATCH_MODEL(Model) \ - if (py::isinstance>(model_obj)) { \ - const auto &mdl = model_obj.cast &>(); \ - return fit_dispatch(mdl, x, y, y_err_obj, n_threads); \ + if (py::isinstance>(model_obj)) { \ + const auto &mdl = \ + model_obj.cast &>(); \ + return fit_dispatch(mdl, x, y, y_err_obj, n_threads); \ } - - AARE_DISPATCH_MODEL(Pol1) - AARE_DISPATCH_MODEL(Pol2) - AARE_DISPATCH_MODEL(Gaussian) - AARE_DISPATCH_MODEL(GaussianErfcPlateau) - AARE_DISPATCH_MODEL(GaussianChargeSharing) - AARE_DISPATCH_MODEL(GaussianChargeSharingKb) - AARE_DISPATCH_MODEL(RisingScurve) - AARE_DISPATCH_MODEL(FallingScurve) - + AARE_FOR_EACH_FIT_MODEL(AARE_DISPATCH_MODEL) #undef AARE_DISPATCH_MODEL +#define AARE_MODEL_NAME(Model) " " #Model throw std::runtime_error( - "Unknown model type. Expected Pol1, Pol2, Gaussian, " - "GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb, " - "RisingScurve or FallingScurve."); + "Unknown model type. Expected one of:" AARE_FOR_EACH_FIT_MODEL( + AARE_MODEL_NAME)); +#undef AARE_MODEL_NAME } void define_fit_bindings(py::module &m) { @@ -346,27 +255,13 @@ void define_fit_bindings(py::module &m) { .value("VarPro", aare::Minimizer::VarPro); // ── Bind model classes ────────────────────────────────────────── - bind_fit_model(m, "Gaussian"); - bind_fit_model(m, "GaussianErfcPlateau"); - bind_fit_model(m, - "GaussianChargeSharing"); - bind_fit_model( - m, "GaussianChargeSharingKb"); - bind_fit_model(m, "RisingScurve"); - bind_fit_model(m, "FallingScurve"); - bind_fit_model(m, "Pol1"); - bind_fit_model(m, "Pol2"); +#define AARE_BIND_FIT_MODEL(Model) \ + bind_fit_model(m, #Model); + AARE_FOR_EACH_FIT_MODEL(AARE_BIND_FIT_MODEL) +#undef AARE_BIND_FIT_MODEL - m.def( - "fit", - [](py::object model_obj, // conversion ok, we want to be able to call - // with any dtype - py::array_t x, - py::array_t y, - py::object y_err_obj, int n_threads) -> py::object { - return dispatch_any_model(model_obj, x, y, y_err_obj, n_threads); - }, - R"( + m.def("fit", &dispatch_any_model, + R"( Fit a model to 1D or 3D data with the minimizer selected by the model's ``minimizer`` property. @@ -395,6 +290,6 @@ void define_fit_bindings(py::module &m) { errors when ``compute_errors`` is enabled; fixed parameters and parameters on a limit report 0. A failed fit returns zeros. )", - py::arg("model"), py::arg("x"), py::arg("y"), - py::arg("y_err") = py::none(), py::arg("n_threads") = 4); + py::arg("model"), py::arg("x"), py::arg("y"), + py::arg("y_err") = py::none(), py::arg("n_threads") = 4); } diff --git a/src/fit/DampedGaussNewton.hpp b/src/fit/DampedGaussNewton.hpp index fdf1670b..f5b7c01b 100644 --- a/src/fit/DampedGaussNewton.hpp +++ b/src/fit/DampedGaussNewton.hpp @@ -1,6 +1,8 @@ // SPDX-License-Identifier: MPL-2.0 #pragma once +#include "LinearAlgebra.hpp" + #include #include #include @@ -9,10 +11,6 @@ namespace aare::detail { -/** @brief Fixed-size vector and square matrix of the solvers, on the stack. */ -template using Vec = std::array; -template using Mat = std::array, N>; - /** @brief Outcome of one of the built-in minimizers. */ struct MinimizerResult { bool valid = false; @@ -21,133 +19,53 @@ struct MinimizerResult { double chi2 = 0.0; }; -/** - * @brief Cholesky factorisation of the leading m x m block of the symmetric - * positive-definite matrix A. Writes the lower triangle of L with - * A = L L^T. Returns false when the block is not positive definite. - * - * The helpers below clamp m to N. The callers never exceed it, but with the - * bound visible gcc's range analysis no longer reports subscripts of the - * one-parameter instantiations as out of range. - */ -template bool cholesky_factor(int m, const Mat &A, Mat &L) { - m = std::min(m, N); - for (int i = 0; i < m; ++i) { - for (int j = 0; j <= i; ++j) { - double sum = A[i][j]; - for (int k = 0; k < j; ++k) - sum -= L[i][k] * L[j][k]; - if (i == j) { - if (!(sum > 0.0)) - return false; - L[i][i] = std::sqrt(sum); - } else { - L[i][j] = sum / L[j][j]; - } - } - } - return true; -} - -/** @brief Solves L L^T z = b for the leading m x m block of the factor L. */ -template -void cholesky_substitute(int m, const Mat &L, const Vec &b, Vec &z) { - m = std::min(m, N); - Vec t{}; - for (int i = 0; i < m; ++i) { - double sum = b[i]; - for (int k = 0; k < i; ++k) - sum -= L[i][k] * t[k]; - t[i] = sum / L[i][i]; - } - for (int i = m - 1; i >= 0; --i) { - double sum = t[i]; - for (int k = i + 1; k < m; ++k) - sum -= L[k][i] * z[k]; - z[i] = sum / L[i][i]; - } -} - -/** @brief Solves A z = b for the leading m x m block of the SPD matrix A. */ -template -bool cholesky_solve(int m, const Mat &A, const Vec &b, Vec &z) { - Mat L{}; - if (!cholesky_factor(m, A, L)) - return false; - cholesky_substitute(m, L, b, z); - return true; -} - -/** @brief Inverse of the leading m x m block of the SPD matrix A. */ -template bool invert_spd(int m, const Mat &A, Mat &inv) { - m = std::min(m, N); - Vec e{}; - Vec col{}; - for (int j = 0; j < m; ++j) { - for (int i = 0; i < m; ++i) - e[i] = (i == j) ? 1.0 : 0.0; - if (!cholesky_solve(m, A, e, col)) - return false; - for (int i = 0; i < m; ++i) - inv[i][j] = col[i]; - } - return true; -} - /** * @brief Gauss-Newton parameter errors sqrt(diag(H^-1)) over the parameters - * that are not skipped, where H is the upper triangle of the Gauss-Newton - * Hessian of chi2 / 2 over all N parameters. Skipped parameters (fixed, or - * frozen on a limit) report 0. Returns false, leaving all errors 0, when the - * block cannot be inverted. + * that are not skipped, where H is the Gauss-Newton Hessian of chi2 / 2 over + * all N parameters. Skipped parameters (fixed, or frozen on a limit) report + * 0. Returns false, leaving all errors 0, when the block cannot be inverted. */ template -bool gauss_newton_errors(const Mat &H, const std::array &skip, +bool gauss_newton_errors(const SymMat &H, const std::array &skip, double *err) { std::fill(err, err + N, 0.0); - Mat A{}; - Mat inv{}; - int m = 0; - for (int a = 0; a < N; ++a) { - if (skip[a]) - continue; - int mb = 0; - for (int c = 0; c < N; ++c) { - if (skip[c]) - continue; - A[m][mb++] = a <= c ? H[a][c] : H[c][a]; - } - ++m; - } - if (m == 0 || !invert_spd(m, A, inv)) + SymMat A; + Cholesky cholesky; + const int m = compress(N, H, skip, A); + if (m == 0 || !cholesky.factor(m, A)) return false; + Vec variance{}; + cholesky.inverse_diag(variance); int i = 0; for (int a = 0; a < N; ++a) { if (skip[a]) continue; - err[a] = std::sqrt(std::max(inv[i][i], 0.0)); + err[a] = std::sqrt(std::max(variance[i], 0.0)); ++i; } return true; } +/** @brief Payload of an evaluator that keeps nothing with a point. */ +struct NoPayload {}; + /** * @brief Damped Gauss-Newton iteration shared by the built-in minimizers. * * The problem is described by an evaluator called as * `eval(const Vec &p, Normal &out)`: at the point p (all NP parameters, * fixed ones at their values) it fills out.F = chi2 / 2, the gradient - * out.g = dF/dp and the upper triangle of the Gauss-Newton Hessian out.H, - * both over all NP parameters, and returns false when the point is not - * admissible. LevenbergMarquardt evaluates residuals and the Jacobian of - * the full model, VariableProjection the projected problem of the nonlinear + * out.g = dF/dp and the Gauss-Newton Hessian out.H, both over all NP + * parameters, and returns false when the point is not + * admissible. It may store a Payload of its own in out.aux, which stays with + * the point. LevenbergMarquardt evaluates residuals and the Jacobian of the + * full model, VariableProjection the projected problem of the nonlinear * parameters. The driver extracts the free parameters and implements * - * - Marquardt damping u * diag(H), starting at u = damping0, with the - * update rule of Nielsen (Madsen, Nielsen and Tingleff, Methods for - * Non-Linear Least Squares Problems, 2004). The damping also bounds the - * accuracy of the last accepted step, so a solver that converges in few - * iterations should start with a small damping. + * - Marquardt damping u * diag(H) with the update rule of Nielsen (Madsen, + * Nielsen and Tingleff, Methods for Non-Linear Least Squares Problems, + * 2004). The start is cautious, u = 0.1: estimated start values can be far + * off, and the first steps are then wild. * - Limits: a trial step that leaves the allowed box is reflected at the * limit. The step truncated at the limit is evaluated as well and the * better of the two is kept. A parameter sitting on a limit whose gradient @@ -162,34 +80,32 @@ bool gauss_newton_errors(const Mat &H, const std::array &skip, * - max_calls bounds the number of evaluations: the initial point and every * trial point count, so an iteration costs one evaluation, or two when a * limit is crossed. - * - With reject_degenerate, a trial point at which the model has lost its - * sensitivity to a free parameter (its Hessian diagonal collapsed by more - * than eight orders of magnitude relative to the current point or to the - * start point) is rejected like an inadmissible point, so the damping - * grows and a shorter step is tried instead. + * - A trial point at which the model has lost its sensitivity to a free + * parameter is rejected like an inadmissible point, so the damping grows + * and a shorter step is tried instead, see degenerate(). * - With polish, a converged fit takes one more undamped Gauss-Newton step - * and keeps it when it lowers chi2, see polish_step(). + * and keeps it when it lowers chi2, see polish_step(). The damping bounds + * the accuracy of the last accepted step, so a solver that converges in + * few iterations needs it. * * After fit() the final point, its Normal and the frozen parameters remain * available for the error computation. */ -template class DampedGaussNewton { +template class DampedGaussNewton { public: static_assert(NP >= 1 && NP <= 16, "DampedGaussNewton keeps NP x NP matrices on the stack"); struct Normal { - double F = 0.0; // chi2 / 2 - Vec g{}; // dF/dp - Mat H{}; // Gauss-Newton Hessian, upper triangle - std::array aux{}; // evaluator payload, kept with the point + double F = 0.0; // chi2 / 2 + Vec g{}; // dF/dp + SymMat H; // Gauss-Newton Hessian + Payload aux{}; // evaluator payload, kept with the point }; struct Options { double tolerance = 0.5; // EDM tolerance, Minuit's convention int max_calls = 100; // budget of evaluations - double damping0 = 0.1; // initial Marquardt damping - bool reject_degenerate = false; bool polish = false; }; @@ -197,13 +113,14 @@ template class DampedGaussNewton { MinimizerResult fit(Eval &eval, const Vec &start, const Vec &lower, const Vec &upper, const std::array &fixed, const Options &opt) { + opt_ = opt; + lower_ = lower; + upper_ = upper; + fixed_ = fixed; + p_ = start; nfree_ = 0; accepted_ = 0; for (int k = 0; k < NP; ++k) { - lower_[k] = lower[k]; - upper_[k] = upper[k]; - p_[k] = start[k]; - fixed_[k] = fixed[k]; free_pos_[k] = -1; if (!fixed[k]) { free_pos_[k] = nfree_; @@ -220,126 +137,24 @@ template class DampedGaussNewton { return res; res.calls = 1; for (int a = 0; a < nfree_; ++a) - start_diag_[a] = cur_.H[free_[a]][free_[a]]; - const double edm_target = 0.002 * opt.tolerance; - const int max_calls = opt.max_calls; - const bool reject_degenerate = opt.reject_degenerate; - double u = opt.damping0; // damping relative to diag(H) - double nu = 2.0; + start_diag_[a] = cur_.H(free_[a], free_[a]); - for (;; ++res.iterations) { - extract_free(); - mark_active(); - if (gmax_ < 1e-12) { - edm(); // the undamped step for converged() - res.valid = - converged(eval, res, max_calls, reject_degenerate, false); - break; - } - // Solve for the damped step first. The solve yields - // g^T (H + u D)^-1 g, a lower bound on the EDM, so the undamped - // factorisation behind edm() is only needed when that bound - // leaves room for convergence. - const bool have_step = solve_step(u); - if (!have_step || bound_ < edm_margin * edm_target) { - const double e = edm(); - if (e < edm_target) { - // A point already deep inside the tolerance needs no - // polishing. - res.valid = - converged(eval, res, max_calls, reject_degenerate, - opt.polish && e >= 0.01 * edm_target); - break; - } - } - if (res.calls >= max_calls) - break; - if (!have_step) { - u *= nu; - nu *= 2.0; - if (u > 1e32) - break; - continue; - } - // A stalled iteration is not a minimum: repeated rejections - // shrink the damped step until it is negligible anywhere. The - // test looks at the step before it is reflected or truncated at - // a limit, so a step that a limit folds back onto the current - // point still gets its truncated alternative evaluated. - if (negligible_step()) { - res.valid = edm() < edm_target; - break; - } - bool moved = false; - bool crossed = false; - for (int a = 0; a < nfree_; ++a) { - const int k = free_[a]; - const double v = q_[a] + delta_[a]; - const double truncated = std::clamp(v, lower_[k], upper_[k]); - double reflected = v; - if (truncated != v) { - crossed = true; - if (v < lower_[k]) - reflected = 2.0 * lower_[k] - v; - if (v > upper_[k]) - reflected = 2.0 * upper_[k] - v; - reflected = std::clamp(reflected, lower_[k], upper_[k]); - } - q_new_[a] = reflected; - q_alt_[a] = truncated; - moved = moved || reflected != q_[a]; - } - // A reflection that lands on the current point needs no - // evaluation; the truncated step below is the only candidate. - bool ok = false; - if (moved) { - ++res.calls; - ok = evaluate(eval, q_new_, trial_) && - !(reject_degenerate && degenerate(trial_)); - } - if (crossed && res.calls < max_calls) { - // Also try the step truncated at the limit, keep the better. - ++res.calls; - const bool ok_alt = evaluate(eval, q_alt_, alt_) && - !(reject_degenerate && degenerate(alt_)); - if (ok_alt && (!ok || alt_.F < trial_.F)) { - q_new_ = q_alt_; - std::swap(trial_, alt_); - ok = true; - } - } - const double rho = ok ? (cur_.F - trial_.F) / predicted_ : -1.0; - if (rho > 0.0) { - // The trial becomes the current point. - q_ = q_new_; - std::swap(cur_, trial_); - for (int a = 0; a < nfree_; ++a) - p_[free_[a]] = q_[a]; - ++accepted_; - const double t = 2.0 * rho - 1.0; - u *= std::max(1.0 / 3.0, 1.0 - t * t * t); - nu = 2.0; - } else { - u *= nu; - nu *= 2.0; - if (u > 1e32) - break; - } + // With every parameter fixed the start point is the result. + res.valid = nfree_ == 0 || iterate(eval, res); + if (res.valid) { + // p_ holds the last evaluated point, which may be a rejected + // trial. + for (int a = 0; a < nfree_; ++a) + p_[free_[a]] = q_[a]; + res.chi2 = 2.0 * cur_.F; } - if (!res.valid) - return res; - - // p_ tracks the last evaluated point, which may be a rejected trial. - for (int a = 0; a < nfree_; ++a) - p_[free_[a]] = q_[a]; - res.chi2 = 2.0 * cur_.F; return res; } /** @brief The final point, all NP parameters. */ const Vec &point() const { return p_; } - /** @brief F, gradient and Hessian at the final point. */ + /** @brief F, gradient, Hessian and payload at the final point. */ const Normal &normal() const { return cur_; } /** @@ -363,6 +178,60 @@ template class DampedGaussNewton { } private: + // The iteration from the evaluated start point. Returns whether it + // converged. + template bool iterate(Eval &eval, MinimizerResult &res) { + const double edm_target = 0.002 * opt_.tolerance; + double u = initial_damping; // damping relative to diag(H) + double nu = 2.0; + for (;; ++res.iterations) { + extract_free(); + mark_active(); + if (gmax_ < 1e-12) { + edm(); // the undamped step for converged() + return converged(eval, res, false); + } + // Solve for the damped step first. The solve yields + // g^T (H + u D)^-1 g, a lower bound on the EDM, so the undamped + // factorisation behind edm() is only needed when that bound + // leaves room for convergence. + const bool have_step = solve_step(u); + if (!have_step || bound_ < edm_margin * edm_target) { + const double e = edm(); + // A point already deep inside the tolerance needs no + // polishing. + if (e < edm_target) + return converged(eval, res, + opt_.polish && e >= 0.01 * edm_target); + } + if (res.calls >= opt_.max_calls) + return false; + double rho = -1.0; // gain ratio; without a step, grow the damping + if (have_step) { + // A stalled iteration is not a minimum: repeated rejections + // shrink the damped step until it is negligible anywhere. + // The test looks at the step before it is reflected or + // truncated at a limit, so a step that a limit folds back + // onto the current point still gets its truncated + // alternative evaluated. + if (negligible_step()) + return edm() < edm_target; + rho = try_step(eval, res); + } + if (rho > 0.0) { + accept_trial(); + const double t = 2.0 * rho - 1.0; + u *= std::max(1.0 / 3.0, 1.0 - t * t * t); + nu = 2.0; + } else { + u *= nu; + nu *= 2.0; + if (u > 1e32) + return false; + } + } + } + template bool evaluate(Eval &eval, const Vec &q, Normal &out) { for (int a = 0; a < nfree_; ++a) @@ -370,18 +239,66 @@ template class DampedGaussNewton { return eval(p_, out); } + // Spends one evaluation on the trial point q. False when the point is + // inadmissible or degenerate. + template + bool try_point(Eval &eval, MinimizerResult &res, const Vec &q, + Normal &out) { + ++res.calls; + return evaluate(eval, q, out) && !degenerate(out); + } + + // Evaluates the damped step from the current point, reflected at the + // limits it crosses, and the step truncated at those limits as well; + // the better of the two ends up in q_new_ and trial_. Returns its gain + // ratio, the decrease of F over the predicted one, or a negative value + // when neither point is admissible. + template double try_step(Eval &eval, MinimizerResult &res) { + bool moved = false; + bool crossed = false; + for (int a = 0; a < nfree_; ++a) { + const int k = free_[a]; + const double v = q_[a] + delta_[a]; + const double truncated = std::clamp(v, lower_[k], upper_[k]); + double reflected = v; + if (truncated != v) { + crossed = true; + if (v < lower_[k]) + reflected = 2.0 * lower_[k] - v; + if (v > upper_[k]) + reflected = 2.0 * upper_[k] - v; + reflected = std::clamp(reflected, lower_[k], upper_[k]); + } + q_new_[a] = reflected; + q_alt_[a] = truncated; + moved = moved || reflected != q_[a]; + } + // A reflection that lands on the current point needs no evaluation; + // the truncated step is then the only candidate. + bool ok = moved && try_point(eval, res, q_new_, trial_); + if (crossed && res.calls < opt_.max_calls && + try_point(eval, res, q_alt_, alt_) && (!ok || alt_.F < trial_.F)) { + q_new_ = q_alt_; + std::swap(trial_, alt_); + ok = true; + } + return ok ? (cur_.F - trial_.F) / predicted_ : -1.0; + } + + // The trial becomes the current point. + void accept_trial() { + q_ = q_new_; + std::swap(cur_, trial_); + ++accepted_; + } + // Gradient and Hessian of the free parameters at the current point. void extract_free() { + compress(NP, cur_.H, fixed_, jtj_); + compress(NP, cur_.g, fixed_, g_); gmax_ = 0.0; - for (int a = 0; a < nfree_; ++a) { - const int ka = free_[a]; - g_[a] = cur_.g[ka]; - for (int b = 0; b < nfree_; ++b) { - const int kb = free_[b]; - jtj_[a][b] = ka <= kb ? cur_.H[ka][kb] : cur_.H[kb][ka]; - } + for (int a = 0; a < nfree_; ++a) gmax_ = std::max(gmax_, std::abs(g_[a])); - } } // A parameter sitting on a limit with the descent direction pointing @@ -394,8 +311,7 @@ template class DampedGaussNewton { // on its lower limit and a higher chi2. But pairing a trial position with // the current gradient could freeze an interior parameter and stop the // fit there. The width no longer reaches its limit that way, because - // LevenbergMarquardt now rejects degenerate trial points as - // VariableProjection does, see degenerate(). + // degenerate trial points are rejected, see degenerate(). void mark_active() { for (int a = 0; a < nfree_; ++a) { const int k = free_[a]; @@ -414,7 +330,7 @@ template class DampedGaussNewton { if (active_[a]) continue; const double scale = - std::abs(q_[a]) + 1.0 / std::sqrt(std::max(jtj_[a][a], 1e-300)); + std::abs(q_[a]) + 1.0 / std::sqrt(std::max(jtj_(a, a), 1e-300)); if (!(std::abs(delta_[a]) <= 1e-8 * scale)) return false; } @@ -422,10 +338,12 @@ template class DampedGaussNewton { } // A trial point where the Hessian diagonal of a free parameter collapsed - // relative to the current point or to the start point: the model no - // longer responds to that parameter there, so the iteration could stop - // at a spurious stationary point (a step function with the width driven - // to zero, for example). The start point catches a collapse spread over + // by more than eight orders of magnitude relative to the current point + // or to the start point: the model no longer responds to that parameter + // there, so the iteration could stop at a spurious stationary point. A + // width driven onto its lower limit, for example, turns a model into a + // step function whose chi2 no longer responds to the width and jumps + // with the position. The start point catches a collapse spread over // many accepted steps, each of which passes the test against its // predecessor. bool degenerate(const Normal &trial) const { @@ -433,39 +351,13 @@ template class DampedGaussNewton { if (active_[a]) continue; const int k = free_[a]; - const double ref = std::max(cur_.H[k][k], start_diag_[a]); - if (!(trial.H[k][k] >= 1e-8 * ref)) + const double ref = std::max(cur_.H(k, k), start_diag_[a]); + if (!(trial.H(k, k) >= 1e-8 * ref)) return true; } return false; } - // Minuit's estimated distance to the minimum of chi2, using the - // Gauss-Newton Hessian: g^T H^-1 g with g the gradient of chi2/2. - // Also keeps the undamped Gauss-Newton step for polish_step(). - double edm() { - Mat A{}; - Vec b{}; - Vec z{}; - const int m = pack(A, b, 0.0); - for (int a = 0; a < nfree_; ++a) - gn_step_[a] = 0.0; - if (m == 0) - return 0.0; - if (!cholesky_solve(m, A, b, z)) - return std::numeric_limits::infinity(); - double e = 0.0; - int i = 0; - for (int a = 0; a < nfree_; ++a) { - if (active_[a]) - continue; - gn_step_[a] = -z[i]; - e += b[i] * z[i]; - ++i; - } - return e; - } - // The convergence test passed. A point that passes it before any step // was accepted may merely lie within the tolerance of chi2's scale, as // the start estimate of an unweighted fit of data of small magnitude @@ -474,14 +366,13 @@ template class DampedGaussNewton { // fit has not converged. With polish, a converged fit takes the same // step when the budget allows. Returns whether the fit is valid. template - bool converged(Eval &eval, MinimizerResult &res, int max_calls, - bool reject_degenerate, bool polish) { + bool converged(Eval &eval, MinimizerResult &res, bool polish) { const bool first = accepted_ == 0 && nfree_ > 0; if (!first && !polish) return true; - if (res.calls >= max_calls) + if (res.calls >= opt_.max_calls) return !first; - polish_step(eval, res, reject_degenerate); + polish_step(eval, res); return true; } @@ -492,80 +383,98 @@ template class DampedGaussNewton { // of the minimum by about the damping factor. The extra evaluation lands // at the minimum to the accuracy of a full Gauss-Newton step. template - void polish_step(Eval &eval, MinimizerResult &res, bool reject_degenerate) { + void polish_step(Eval &eval, MinimizerResult &res) { for (int a = 0; a < nfree_; ++a) { const int k = free_[a]; q_new_[a] = std::clamp(q_[a] + gn_step_[a], lower_[k], upper_[k]); } - ++res.calls; - const bool ok = evaluate(eval, q_new_, trial_) && - !(reject_degenerate && degenerate(trial_)) && - trial_.F < cur_.F; - if (ok) { - q_ = q_new_; - std::swap(cur_, trial_); + if (try_point(eval, res, q_new_, trial_) && trial_.F < cur_.F) { + accept_trial(); extract_free(); - } - for (int a = 0; a < nfree_; ++a) - p_[free_[a]] = q_[a]; - if (ok) mark_active(); + } } // Copy the non-frozen block of H (+ u * diag) and g into A and b. - int pack(Mat &A, Vec &b, double u) const { + int pack(SymMat &A, Vec &b, double u) const { + const int m = compress(nfree_, jtj_, active_, A); + compress(nfree_, g_, active_, b); double mean_diag = 0.0; for (int a = 0; a < nfree_; ++a) - mean_diag += jtj_[a][a]; + mean_diag += jtj_(a, a); mean_diag = std::max(mean_diag / std::max(nfree_, 1), 1e-300); - int m = 0; + int i = 0; for (int a = 0; a < nfree_; ++a) { if (active_[a]) continue; - int mb = 0; - for (int c = 0; c < nfree_; ++c) { - if (active_[c]) - continue; - A[m][mb++] = jtj_[a][c]; - } - A[m][m] += u * std::max(jtj_[a][a], 1e-10 * mean_diag); - b[m] = g_[a]; - ++m; + A(i, i) += u * std::max(jtj_(a, a), 1e-10 * mean_diag); + ++i; } return m; } + // Solves (H + u D) z = g for the non-frozen parameters. Writes + // step = -z, 0 for a frozen parameter, and gz = g^T z. Returns the + // number of parameters solved for, or -1 when the matrix is not + // positive definite. + int solve(double u, Vec &step, double &gz) const { + SymMat A; + Vec b{}; + Vec z{}; + Cholesky cholesky; + const int m = pack(A, b, u); + for (int a = 0; a < nfree_; ++a) + step[a] = 0.0; + gz = 0.0; + if (m > 0 && !cholesky.factor(m, A)) + return -1; + cholesky.solve(b, z); + int i = 0; + for (int a = 0; a < nfree_; ++a) { + if (active_[a]) + continue; + step[a] = -z[i]; + gz += b[i] * z[i]; + ++i; + } + return m; + } + + // Minuit's estimated distance to the minimum of chi2, using the + // Gauss-Newton Hessian: g^T H^-1 g with g the gradient of chi2/2. + // Also keeps the undamped Gauss-Newton step for polish_step(). + double edm() { + double e = 0.0; + if (solve(0.0, gn_step_, e) < 0) + return std::numeric_limits::infinity(); + return e; + } + // Damped Gauss-Newton step (H + u D) delta = -g for the non-frozen // parameters, the predicted decrease of F for that step, and // bound_ = g^T (H + u D)^-1 g, which is at most the EDM because u D is // positive semi-definite. bool solve_step(double u) { - Mat A{}; - Vec b{}; - Vec z{}; - const int m = pack(A, b, u); - for (int a = 0; a < nfree_; ++a) - delta_[a] = 0.0; predicted_ = 0.0; - bound_ = 0.0; - if (m == 0) - return true; - if (!cholesky_solve(m, A, b, z)) + const int m = solve(u, delta_, bound_); + if (m < 0) return false; - int i = 0; for (int a = 0; a < nfree_; ++a) { if (active_[a]) continue; - delta_[a] = -z[i]; - bound_ += b[i] * z[i]; predicted_ += 0.5 * delta_[a] * - (u * std::max(jtj_[a][a], 1e-300) * delta_[a] - g_[a]); - ++i; + (u * std::max(jtj_(a, a), 1e-300) * delta_[a] - g_[a]); } - return predicted_ > 0.0; + return m == 0 || predicted_ > 0.0; } + static constexpr double initial_damping = 0.1; + // The two solves that bound the EDM differ by rounding only when the + // damping is negligible; the margin keeps the stop decision exact. + static constexpr double edm_margin = 1.0 + 1e-6; + + Options opt_; int nfree_ = 0; int accepted_ = 0; // trial points accepted in this fit std::array free_{}; @@ -574,18 +483,15 @@ template class DampedGaussNewton { std::array active_{}; Vec lower_{}; Vec upper_{}; - Vec p_{}; - Vec q_{}; - Vec q_new_{}; - Vec q_alt_{}; + Vec p_{}; // all parameters at the last evaluated point + Vec q_{}; // free parameters at the current point + Vec q_new_{}; // trial point + Vec q_alt_{}; // trial point truncated at the limits Vec g_{}; Vec delta_{}; Vec gn_step_{}; Vec start_diag_{}; // Hessian diagonal at the start point - Mat jtj_{}; - // The two solves that bound the EDM differ by rounding only when the - // damping is negligible; the margin keeps the stop decision exact. - static constexpr double edm_margin = 1.0 + 1e-6; + SymMat jtj_; // Hessian of the free parameters double gmax_ = 0.0; double predicted_ = 0.0; double bound_ = 0.0; diff --git a/src/fit/Fit.cpp b/src/fit/Fit.cpp index 8b34c177..5bccf2b4 100644 --- a/src/fit/Fit.cpp +++ b/src/fit/Fit.cpp @@ -172,50 +172,47 @@ template class PixelFitter { return fail(par_out, err_out, chi2); const auto start = start_values(model, x, y); - - if (model.minimizer() == Minimizer::VarPro) { - // A model without the separable structure, a pixel that does - // not converge and a pixel whose linear solution violates a - // limit fall back to the full solver. + const Minimizer minimizer = model.minimizer(); + bool valid = false; + if (minimizer == Minimizer::Migrad || minimizer == Minimizer::Fumili) { + valid = fit_pixel_minuit2(model, x, y, y_err, start, par_out, + err_out, chi2); + } else { + // VarPro falls back to the full solver for a model without the + // separable structure, a pixel that does not converge and a + // pixel whose linear solution violates a limit. if constexpr (model::is_separable::value) { - const auto res = - vp_.fit(model, x, y, y_err, start, par_out, err_out); - if (res.valid && finite(par_out, err_out, res.chi2)) { - chi2 = res.chi2; - return true; - } + valid = + minimizer == Minimizer::VarPro && + run(vp_, model, x, y, y_err, start, par_out, err_out, chi2); } - return fit_lm(model, x, y, y_err, start, par_out, err_out, chi2); + valid = valid || + run(lm_, model, x, y, y_err, start, par_out, err_out, chi2); } - if (model.minimizer() == Minimizer::LevenbergMarquardt) - return fit_lm(model, x, y, y_err, start, par_out, err_out, chi2); - return fit_pixel_minuit2(model, x, y, y_err, start, par_out, err_out, - chi2); + return valid || fail(par_out, err_out, chi2); } private: - // The built-in solvers do not guard every operation against NaN or - // infinity, so a result is checked before it is reported as a valid fit. - static bool finite(const double *par, const double *err, double chi2) { - bool ok = std::isfinite(chi2); - for (std::size_t k = 0; k < Model::npar; ++k) - ok = ok && std::isfinite(par[k]) && (!err || std::isfinite(err[k])); + // Runs a built-in solver. They do not guard every operation against NaN + // or infinity, so a result only counts as a valid fit when it is finite. + template + static bool run(Solver &solver, const FitModel &model, + NDView x, NDView y, + NDView y_err, + const std::array &start, + double *par_out, double *err_out, double &chi2) { + const auto res = + solver.fit(model, x, y, y_err, start, par_out, err_out); + bool ok = res.valid && std::isfinite(res.chi2); + for (std::size_t k = 0; ok && k < Model::npar; ++k) + ok = std::isfinite(par_out[k]) && + (!err_out || std::isfinite(err_out[k])); + if (ok) + chi2 = res.chi2; return ok; } - bool fit_lm(const FitModel &model, NDView x, - NDView y, NDView y_err, - const std::array &start, double *par_out, - double *err_out, double &chi2) { - const auto res = lm_.fit(model, x, y, y_err, start, par_out, err_out); - if (res.valid && finite(par_out, err_out, res.chi2)) { - chi2 = res.chi2; - return true; - } - return fail(par_out, err_out, chi2); - } - - // A failed fit reports zeros. + // A failed fit reports zeros, whichever minimizer it came from. static bool fail(double *par_out, double *err_out, double &chi2) { std::fill(par_out, par_out + Model::npar, 0.0); if (err_out) @@ -311,25 +308,18 @@ void fit_3d(const FitModel &model, NDView x, // NOLINTBEGIN #define AARE_INSTANTIATE_FIT(Model) \ - template class FitModel; \ - template NDArray fit_pixel( \ - const FitModel &, NDView, NDView, \ + template class FitModel; \ + template NDArray fit_pixel( \ + const FitModel &, NDView, NDView, \ NDView); \ - template NDArray fit_pixel( \ - const FitModel &, NDView, NDView); \ - template void fit_3d(const FitModel &, NDView, \ - NDView, NDView, \ - NDView, NDView, \ - NDView, int); + template NDArray fit_pixel( \ + const FitModel &, NDView, NDView); \ + template void fit_3d( \ + const FitModel &, NDView, NDView, \ + NDView, NDView, NDView, \ + NDView, int); -AARE_INSTANTIATE_FIT(model::Gaussian) -AARE_INSTANTIATE_FIT(model::GaussianErfcPlateau) -AARE_INSTANTIATE_FIT(model::GaussianChargeSharing) -AARE_INSTANTIATE_FIT(model::GaussianChargeSharingKb) -AARE_INSTANTIATE_FIT(model::Pol1) -AARE_INSTANTIATE_FIT(model::Pol2) -AARE_INSTANTIATE_FIT(model::RisingScurve) -AARE_INSTANTIATE_FIT(model::FallingScurve) +AARE_FOR_EACH_FIT_MODEL(AARE_INSTANTIATE_FIT) #undef AARE_INSTANTIATE_FIT // NOLINTEND diff --git a/src/fit/Fit.test.cpp b/src/fit/Fit.test.cpp index 992500e8..69c1dceb 100644 --- a/src/fit/Fit.test.cpp +++ b/src/fit/Fit.test.cpp @@ -620,6 +620,18 @@ TEST_CASE("Fit a data cube with errors matches per-pixel fits", "[fit]") { .margin(1e-6)); } } + + // The error output is optional for a weighted fit as well, but has to + // match the data when it is passed. + NDArray par_only({rows, cols, 3}); + aare::fit_3d(model, x.view(), y.view(), y_err.view(), par_only.view(), + NDView{}, chi2.view(), 3); + for (ssize_t i = 0; i < par.size(); ++i) + CHECK(par_only[i] == Approx(par[i]).margin(1e-9)); + NDArray short_err({rows, cols, 2}); + CHECK_THROWS_AS(aare::fit_3d(model, x.view(), y.view(), y_err.view(), + par.view(), short_err.view(), chi2.view(), 3), + std::runtime_error); } TEST_CASE("Fitting with every parameter fixed returns the fixed values", @@ -778,57 +790,43 @@ TEST_CASE("All minimizers agree on a noisy data cube", "[fit]") { namespace { // The separable structure of a model must reproduce eval and eval_and_grad: -// f = sum_j par[linear_par[j]] * phi[j], df/dpar[linear_par[j]] = phi[j] and -// df/dpar[nonlinear[k]] = sum_j par[linear_par[j]] * dphi[k][j]. +// with the columns phi_j and dphi_kj of basis_columns, +// f = sum_j par[linear_par[j]] * phi_j, df/dpar[linear_par[j]] = phi_j and +// df/dpar[nonlinear[k]] = sum_j par[linear_par[j]] * dphi_kj. template void check_separable(const std::vector &par, const std::vector &xs) { using Traits = aare::model::separable_traits; static_assert(aare::model::is_separable::value); - typename Traits::Basis phi{}; - typename Traits::BasisGrad dphi{}; + const std::size_t n = xs.size(); + std::vector phi(Traits::nlin * xs.size()); + std::vector dphi(Traits::nnl * Traits::nlin * xs.size()); + Model::basis_columns(xs.data(), static_cast(xs.size()), par, + phi.data(), dphi.data()); std::array g{}; double f = 0.0; - for (const double x : xs) { - INFO("x = " << x); - Model::basis_and_grad(x, par, phi, dphi); - Model::eval_and_grad(x, par, f, g); + for (std::size_t i = 0; i < n; ++i) { + INFO("x = " << xs[i]); + Model::eval_and_grad(xs[i], par, f, g); double sum = 0.0; for (std::size_t j = 0; j < Traits::nlin; ++j) { const auto idx = Model::linear_par[j]; - sum += par[idx] * phi[j]; - CHECK(g[idx] == Approx(phi[j]).epsilon(1e-12).margin(1e-14)); + sum += par[idx] * phi[j * n + i]; + CHECK(g[idx] == + Approx(phi[j * n + i]).epsilon(1e-12).margin(1e-14)); } - CHECK(Model::eval(x, par) == Approx(sum).epsilon(1e-12).margin(1e-14)); + CHECK(Model::eval(xs[i], par) == + Approx(sum).epsilon(1e-12).margin(1e-14)); CHECK(f == Approx(sum).epsilon(1e-12).margin(1e-14)); for (std::size_t k = 0; k < Traits::nnl; ++k) { double d = 0.0; for (std::size_t j = 0; j < Traits::nlin; ++j) - d += par[Model::linear_par[j]] * dphi[k][j]; + d += par[Model::linear_par[j]] * + dphi[(k * Traits::nlin + j) * n + i]; CHECK(g[Traits::nonlinear_par[k]] == Approx(d).epsilon(1e-12).margin(1e-14)); } } - // The bulk basis pass must agree with the per-point one. - if constexpr (aare::model::has_basis_columns::value) { - const auto n = static_cast(xs.size()); - std::vector cols(Traits::nlin * xs.size()); - std::vector dcols(Traits::nnl * Traits::nlin * xs.size()); - Model::basis_columns(xs.data(), n, par, cols.data(), dcols.data()); - for (ssize_t i = 0; i < n; ++i) { - INFO("x = " << xs[static_cast(i)]); - Model::basis_and_grad(xs[static_cast(i)], par, phi, - dphi); - for (std::size_t j = 0; j < Traits::nlin; ++j) - CHECK(cols[j * xs.size() + static_cast(i)] == - Approx(phi[j]).epsilon(1e-12).margin(1e-14)); - for (std::size_t k = 0; k < Traits::nnl; ++k) - for (std::size_t j = 0; j < Traits::nlin; ++j) - CHECK(dcols[(k * Traits::nlin + j) * xs.size() + - static_cast(i)] == - Approx(dphi[k][j]).epsilon(1e-12).margin(1e-14)); - } - } } } // namespace diff --git a/src/fit/FitHelpers.hpp b/src/fit/FitHelpers.hpp index a82599cf..98963715 100644 --- a/src/fit/FitHelpers.hpp +++ b/src/fit/FitHelpers.hpp @@ -7,8 +7,8 @@ #include #include -// Helpers shared by the built-in solver and the optional Minuit2 backend so -// that both apply the same start-value precedence and shape checks. +// Helpers shared by the built-in solvers and the Minuit2 backend so that all +// apply the same start-value precedence, budget and shape checks. namespace aare::detail { @@ -35,6 +35,16 @@ std::array start_values(const FitModel &model, return start; } +/** + * @brief Evaluation budget of the built-in minimizers. As in Minuit2, + * max_calls = 0 selects the default 200 + 100 * npar + 5 * npar^2. + */ +template int evaluation_budget(const FitModel &model) { + constexpr int npar = static_cast(Model::npar); + return model.max_calls() > 0 ? static_cast(model.max_calls()) + : 200 + 100 * npar + 5 * npar * npar; +} + /** @brief Validate the array sizes passed to fit_pixel. */ inline void check_fit_pixel_shapes(NDView x, NDView y, NDView y_err) { @@ -63,21 +73,16 @@ void check_fit_3d_shapes(NDView x, NDView y, if (chi2_out.shape(0) != y.shape(0) || chi2_out.shape(1) != y.shape(1)) throw std::runtime_error("chi2_out must have shape [rows, cols]."); - const bool err_out_matches = err_out.shape(0) == y.shape(0) && - err_out.shape(1) == y.shape(1) && - err_out.shape(2) == npar; + if (y_err.size() > 0 && + (y.shape(0) != y_err.shape(0) || y.shape(1) != y_err.shape(1) || + y.shape(2) != y_err.shape(2))) + throw std::runtime_error( + "fit_3d: y and y_err must have identical shape."); - if (y_err.size() > 0) { - if (y.shape(0) != y_err.shape(0) || y.shape(1) != y_err.shape(1) || - y.shape(2) != y_err.shape(2)) - throw std::runtime_error( - "fit_3d: y and y_err must have identical shape."); - if (!err_out_matches) - throw std::runtime_error( - "err_out must have shape [rows, cols, npar]."); - } else if (err_out.size() > 0 && !err_out_matches) { + if (err_out.size() > 0 && + (err_out.shape(0) != y.shape(0) || err_out.shape(1) != y.shape(1) || + err_out.shape(2) != npar)) throw std::runtime_error("err_out must have shape [rows, cols, npar]."); - } } } // namespace aare::detail diff --git a/src/fit/FitMinuit2.cpp b/src/fit/FitMinuit2.cpp index 326c4b11..3291cd0d 100644 --- a/src/fit/FitMinuit2.cpp +++ b/src/fit/FitMinuit2.cpp @@ -28,17 +28,10 @@ ROOT::Minuit2::MnUserParameters user_parameters(const FitModel &model, NDView x, NDView y, const std::array &start) { - constexpr std::size_t npar = Model::npar; - - double x_range = 0.0; - double y_range = 0.0; - double slope_scale = 0.0; - minuit2::compute_ranges(x, y, x_range, y_range, slope_scale); - std::array steps{}; - minuit2::Steps::compute(start, x_range, y_range, slope_scale, steps); + const auto steps = minuit2::steps(start, x, y); ROOT::Minuit2::MnUserParameters upar; - for (std::size_t i = 0; i < npar; ++i) { + for (std::size_t i = 0; i < Model::npar; ++i) { const auto idx = static_cast(i); upar.Add(Model::param_info[i].name, start[i], steps[i]); const double lo = model.lower_limit(idx); @@ -102,12 +95,6 @@ bool fit_pixel_minuit2(const FitModel &model, NDView x, NDView y, NDView y_err, const std::array &start, double *par_out, double *err_out, double &chi2) { - constexpr std::size_t npar = Model::npar; - std::fill(par_out, par_out + npar, 0.0); - if (err_out) - std::fill(err_out, err_out + npar, 0.0); - chi2 = 0.0; - if (!Model::is_valid(std::vector(start.begin(), start.end()))) return false; @@ -148,19 +135,12 @@ bool fit_pixel_minuit2(const FitModel &model, NDView x, // NOLINTBEGIN #define AARE_INSTANTIATE_MINUIT2_FIT(Model) \ - template bool fit_pixel_minuit2( \ - const FitModel &, NDView, NDView, \ - NDView, const std::array &, double *, \ - double *, double &); + template bool fit_pixel_minuit2( \ + const FitModel &, NDView, NDView, \ + NDView, const std::array &, \ + double *, double *, double &); -AARE_INSTANTIATE_MINUIT2_FIT(model::Gaussian) -AARE_INSTANTIATE_MINUIT2_FIT(model::GaussianErfcPlateau) -AARE_INSTANTIATE_MINUIT2_FIT(model::GaussianChargeSharing) -AARE_INSTANTIATE_MINUIT2_FIT(model::GaussianChargeSharingKb) -AARE_INSTANTIATE_MINUIT2_FIT(model::Pol1) -AARE_INSTANTIATE_MINUIT2_FIT(model::Pol2) -AARE_INSTANTIATE_MINUIT2_FIT(model::RisingScurve) -AARE_INSTANTIATE_MINUIT2_FIT(model::FallingScurve) +AARE_FOR_EACH_FIT_MODEL(AARE_INSTANTIATE_MINUIT2_FIT) #undef AARE_INSTANTIATE_MINUIT2_FIT // NOLINTEND diff --git a/src/fit/FitMinuit2.hpp b/src/fit/FitMinuit2.hpp index c0240be8..9ee78d5a 100644 --- a/src/fit/FitMinuit2.hpp +++ b/src/fit/FitMinuit2.hpp @@ -15,10 +15,12 @@ namespace aare::detail { * @brief Fit one pixel with the Minuit2 minimizer selected by the model. * * @param start Starting values, see start_values() in FitHelpers.hpp. - * @param par_out Receives npar fitted values; zeros on failure. - * @param err_out Receives npar errors when non-null; zeros on failure. - * Fixed parameters and parameters on a limit report 0. - * @param chi2 Receives the chi-squared at the minimum; 0 on failure. + * @param par_out Receives npar fitted values. + * @param err_out Receives npar errors when non-null. Fixed parameters and + * parameters on a limit report 0. + * @param chi2 Receives the chi-squared at the minimum. + * + * The outputs are untouched when the fit fails. * @return true when Minuit2 reports a valid minimum. */ template diff --git a/src/fit/LevenbergMarquardt.hpp b/src/fit/LevenbergMarquardt.hpp index 8a6a7213..a8e9bf2f 100644 --- a/src/fit/LevenbergMarquardt.hpp +++ b/src/fit/LevenbergMarquardt.hpp @@ -2,6 +2,7 @@ #pragma once #include "DampedGaussNewton.hpp" +#include "FitHelpers.hpp" #include "aare/FitModel.hpp" #include "aare/NDView.hpp" @@ -26,9 +27,8 @@ namespace aare::detail { * equations g = J^T r and J^T J, which DampedGaussNewton iterates on: it * handles fixed parameters, limits, damping and the convergence test, see * there. Rejected trials are rare with its damping rule, and a rejected one - * only wastes a single pass. Trial points at which the model loses its - * sensitivity to a free parameter are rejected. As in Minuit2, max_calls = 0 - * selects the default budget 200 + 100 * npar + 5 * npar^2. + * only wastes a single pass. max_calls counts evaluations, see + * evaluation_budget(). * * Errors are Gauss-Newton estimates, sqrt(diag((J^T J)^-1)), over the free * parameters that do not sit on a limit. Fixed and limit-bound parameters @@ -40,7 +40,6 @@ template class LevenbergMarquardt { static_assert(Model::npar >= 1 && Model::npar <= 16, "LevenbergMarquardt keeps npar x npar matrices on the stack"); - using Result = MinimizerResult; using Driver = DampedGaussNewton; using Normal = typename Driver::Normal; @@ -52,18 +51,28 @@ template class LevenbergMarquardt { * @param y Measured values. * @param s Per-point uncertainties, empty view for an unweighted fit. * @param start Starting values, see start_values() in FitHelpers.hpp. - * @param par_out Receives npar fitted values; zeros on failure. - * @param err_out Receives npar errors when non-null; zeros on failure. + * @param par_out Receives npar fitted values; untouched on failure. + * @param err_out Receives npar errors when non-null; untouched on + * failure. */ - Result fit(const FitModel &model, NDView x, - NDView y, NDView s, - const std::array &start, double *par_out, - double *err_out) { + MinimizerResult fit(const FitModel &model, NDView x, + NDView y, NDView s, + const std::array &start, + double *par_out, double *err_out) { x_ = x; y_ = y; - s_ = s; weighted_ = s.size() > 0; const auto n = static_cast(x.size()); + r_.resize(n); + J_.resize(n * static_cast(npar)); + if (weighted_) { + // 1 / s_i once per pixel rather than once per evaluation. + w_.resize(n); + for (std::size_t i = 0; i < n; ++i) { + const double si = s[static_cast(i)]; + w_[i] = si != 0.0 ? 1.0 / si : 0.0; + } + } Vec lower{}; Vec upper{}; std::array fixed{}; @@ -73,41 +82,16 @@ template class LevenbergMarquardt { upper[k] = model.upper_limit(idx); fixed[k] = model.is_user_fixed(idx); } - r_.resize(n); - J_.resize(n * static_cast(npar)); - if (weighted_) { - // 1 / s_i once per pixel rather than once per evaluation. - w_.resize(n); - for (std::size_t i = 0; i < n; ++i) - w_[i] = s_[static_cast(i)] != 0.0 - ? 1.0 / s_[static_cast(i)] - : 0.0; - } - const int max_calls = model.max_calls() > 0 - ? static_cast(model.max_calls()) - : 200 + 100 * npar + 5 * npar * npar; + typename Driver::Options opt; + opt.tolerance = model.tolerance(); + opt.max_calls = evaluation_budget(model); auto eval = [this](const Vec &p, Normal &out) { return evaluate(p, out); }; - typename Driver::Options opt; - opt.tolerance = model.tolerance(); - opt.max_calls = max_calls; - // Start cautiously: the estimated start values of the linear - // parameters can be far off, and the first steps are then wild. - opt.damping0 = 0.1; - // A width driven onto its lower limit turns a model into a step - // function whose chi2 no longer responds to the width and jumps - // with the position; the iteration then stalls there. Such trial - // points are rejected, as in VariableProjection. - opt.reject_degenerate = true; - const Result res = driver_.fit(eval, start, lower, upper, fixed, opt); - if (!res.valid) { - std::fill(par_out, par_out + npar, 0.0); - if (err_out) - std::fill(err_out, err_out + npar, 0.0); + const auto res = driver_.fit(eval, start, lower, upper, fixed, opt); + if (!res.valid) return res; - } std::copy(driver_.point().begin(), driver_.point().end(), par_out); if (err_out) driver_.errors(err_out); @@ -120,7 +104,7 @@ template class LevenbergMarquardt { } // Residuals r = (y - f) / s and the Jacobian J = dr/dp at the point p, - // reduced to F = chi2 / 2, g = J^T r and the upper triangle of J^T J. + // reduced to F = chi2 / 2, g = J^T r and J^T J. // The sums run over all parameters with fixed-size loops, so the // compiler can unroll them and keep the accumulators in registers. // Returns false when the model rejects the parameters. @@ -142,26 +126,21 @@ template class LevenbergMarquardt { Ji[k] = -g[k] * w; } Vec G{}; - Mat H{}; + SymMat H; for (ssize_t i = 0; i < n; ++i) { const double *Ji = &J_[static_cast(i) * npar]; const double ri = r_[i]; - for (int a = 0; a < npar; ++a) { + for (int a = 0; a < npar; ++a) G[a] += Ji[a] * ri; - for (int b = a; b < npar; ++b) - H[a][b] += Ji[a] * Ji[b]; - } + H.rank1_update(Ji); } out.F = F; - for (int a = 0; a < npar; ++a) { - out.g[a] = G[a]; - for (int b = a; b < npar; ++b) - out.H[a][b] = H[a][b]; - } + out.g = G; + out.H = H; return true; } - NDView x_, y_, s_; + NDView x_, y_; bool weighted_ = false; Driver driver_; std::vector r_; diff --git a/src/fit/LinearAlgebra.hpp b/src/fit/LinearAlgebra.hpp new file mode 100644 index 00000000..962fa7d3 --- /dev/null +++ b/src/fit/LinearAlgebra.hpp @@ -0,0 +1,173 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +// Small fixed-size linear algebra of the built-in minimizers. Everything +// lives on the stack and the sizes are compile-time constants, so the loops +// unroll. Function templates cannot deduce N from the std::array aliases +// (int against std::size_t), hence the explicit compress(...). + +#include "aare/defs.hpp" +#include +#include +#include + +namespace aare::detail { + +/** @brief Fixed-size vector. */ +template using Vec = std::array; + +/** @brief Fixed-size R x C matrix, square by default. */ +template using Mat = std::array, R>; + +/** + * @brief Symmetric N x N matrix. Only the upper triangle is stored and + * updated; (i, j) and (j, i) name the same element. + */ +template struct SymMat { + Mat upper{}; // [i][j] with i <= j; the rest is never read + + double &operator()(int i, int j) { + return upper[std::min(i, j)][std::max(i, j)]; + } + double operator()(int i, int j) const { + return upper[std::min(i, j)][std::max(i, j)]; + } + + /** + * @brief Adds w * x x^T, for x indexable with N elements. It is inlined + * because it runs once per data point. + */ + template + ALWAYS_INLINE void rank1_update(const X &x, double w = 1.0) { + for (int a = 0; a < N; ++a) { + const double xa = x[a] * w; + for (int b = a; b < N; ++b) + upper[a][b] += xa * x[b]; + } + } +}; + +/** + * @brief Cholesky factorisation A = L L^T of the leading m x m block of a + * symmetric positive-definite matrix, for any number of solves with it. + * + * The methods clamp m to N. The callers never exceed it, but with the bound + * visible gcc's range analysis no longer reports subscripts of the + * one-parameter instantiations as out of range. + */ +template class Cholesky { + public: + /** @brief Factorises; false when the block is not positive definite. */ + bool factor(int m, const SymMat &A) { + m = std::min(m, N); + m_ = 0; + for (int i = 0; i < m; ++i) { + diag_[i] = A.upper[i][i]; + for (int j = 0; j <= i; ++j) { + double sum = A.upper[j][i]; + for (int k = 0; k < j; ++k) + sum -= L_[i][k] * L_[j][k]; + if (i == j) { + if (!(sum > 0.0)) + return false; + L_[i][i] = std::sqrt(sum); + } else { + L_[i][j] = sum / L_[j][j]; + } + } + } + m_ = m; + return true; + } + + /** @brief Solves A z = b for the first m entries of b and z. */ + void solve(const Vec &b, Vec &z) const { + const int m = std::min(m_, N); + Vec t{}; + for (int i = 0; i < m; ++i) { + double sum = b[i]; + for (int k = 0; k < i; ++k) + sum -= L_[i][k] * t[k]; + t[i] = sum / L_[i][i]; + } + for (int i = m - 1; i >= 0; --i) { + double sum = t[i]; + for (int k = i + 1; k < m; ++k) + sum -= L_[k][i] * z[k]; + z[i] = sum / L_[i][i]; + } + } + + /** + * @brief True when every pivot satisfies L_ii^2 >= ratio * A_ii. A pivot + * far below its diagonal element costs a solution about + * log10(A_ii / L_ii^2) digits. + */ + bool pivots_at_least(double ratio) const { + const int m = std::min(m_, N); + for (int i = 0; i < m; ++i) { + if (!(L_[i][i] * L_[i][i] >= ratio * diag_[i])) + return false; + } + return true; + } + + /** @brief The first m diagonal elements of A^-1. */ + void inverse_diag(Vec &d) const { + const int m = std::min(m_, N); + Vec e{}; + Vec col{}; + for (int j = 0; j < m; ++j) { + for (int i = 0; i < m; ++i) + e[i] = (i == j) ? 1.0 : 0.0; + solve(e, col); + d[j] = col[j]; + } + } + + private: + int m_ = 0; + Mat L_{}; // lower triangle + Vec diag_{}; // diagonal of A +}; + +/** + * @brief Packs the rows and columns of H that drop does not mark, among the + * first n, into the leading block of A. Returns the size of that block. + */ +template +int compress(int n, const SymMat &H, const std::array &drop, + SymMat &A) { + n = std::min(n, N); + int m = 0; + for (int a = 0; a < n; ++a) { + if (drop[a]) + continue; + int mc = m; + for (int c = a; c < n; ++c) { + if (drop[c]) + continue; + A.upper[m][mc++] = H.upper[a][c]; + } + ++m; + } + return m; +} + +/** + * @brief Packs the entries of g that drop does not mark, among the first n, + * into the leading entries of b. Returns their number. + */ +template +int compress(int n, const Vec &g, const std::array &drop, + Vec &b) { + n = std::min(n, N); + int m = 0; + for (int a = 0; a < n; ++a) { + if (!drop[a]) + b[m++] = g[a]; + } + return m; +} + +} // namespace aare::detail diff --git a/src/fit/MinuitSteps.hpp b/src/fit/MinuitSteps.hpp index d5aa8bbe..4b694e4c 100644 --- a/src/fit/MinuitSteps.hpp +++ b/src/fit/MinuitSteps.hpp @@ -1,126 +1,116 @@ // SPDX-License-Identifier: MPL-2.0 #pragma once -// Data-driven Minuit step sizes for the optional Minuit2 backend. These -// estimators used to live in aare/Models.hpp; the built-in solver does not -// need them. +// Data-driven Minuit step sizes for the Minuit2 backend. The built-in +// solvers do not need them. #include "aare/Models.hpp" #include "aare/NDView.hpp" #include #include #include +#include namespace aare::minuit2 { -/** - * @brief Compute data-range statistics used by the step-size estimators. - * - * Model-independent, called once per pixel. - */ -inline void compute_ranges(NDView x, NDView y, - double &x_range, double &y_range, - double &slope_scale) { +/** @brief What a parameter measures, which sets the size of its first step. */ +enum class Scale { + Height, // a y value: offset or amplitude + Slope, // dy/dx + Curvature, // d2y/dx2 + Position, // an x value: position or width + Fraction, // dimensionless, relative to another parameter + Ratio, // dimensionless and close to one +}; + +/** @brief The scale of every parameter of a fit model. */ +template struct Scales; + +template <> struct Scales { + static constexpr std::array value = {Scale::Height, Scale::Slope}; +}; + +template <> struct Scales { + static constexpr std::array value = {Scale::Height, Scale::Slope, + Scale::Curvature}; +}; + +template <> struct Scales { + static constexpr std::array value = { + Scale::Height, Scale::Position, Scale::Position}; +}; + +template <> struct Scales { + static constexpr std::array value = { + Scale::Height, Scale::Height, Scale::Position, Scale::Position}; +}; + +template <> struct Scales { + static constexpr std::array value = { + Scale::Height, Scale::Slope, Scale::Position, + Scale::Position, Scale::Height, Scale::Fraction}; +}; + +template <> struct Scales { + static constexpr std::array value = { + Scale::Height, Scale::Slope, Scale::Position, Scale::Position, + Scale::Height, Scale::Fraction, Scale::Ratio, Scale::Fraction}; +}; + +template <> struct Scales { + static constexpr std::array value = { + Scale::Height, Scale::Slope, Scale::Position, + Scale::Position, Scale::Height, Scale::Slope}; +}; + +template <> +struct Scales : Scales {}; + +/** @brief Extent of one pixel's data, which the steps are scaled to. */ +struct Ranges { + double x; + double y; + double slope; // y / x +}; + +inline Ranges compute_ranges(NDView x, NDView y) { const auto [x_min, x_max] = std::minmax_element(x.begin(), x.end()); const auto [y_min, y_max] = std::minmax_element(y.begin(), y.end()); - x_range = std::max(*x_max - *x_min, 1e-9); - y_range = std::max(*y_max - *y_min, 1e-9); - - slope_scale = std::max(y_range / x_range, 1e-9); + const double x_range = std::max(*x_max - *x_min, 1e-9); + const double y_range = std::max(*y_max - *y_min, 1e-9); + return {x_range, y_range, std::max(y_range / x_range, 1e-9)}; } -/** @brief Per-model Minuit step sizes, specialised for every fit model. */ -template struct Steps; - -template <> struct Steps { - static void compute(const std::array &start, - [[maybe_unused]] double x_range, double y_range, - double slope_scale, std::array &steps) { - steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range); - steps[1] = 0.1 * slope_scale; +/** @brief Step of a parameter of the given scale that starts at start. */ +inline double step(Scale scale, double start, const Ranges &range) { + switch (scale) { + case Scale::Height: + return std::max(0.1 * std::abs(start), 0.1 * range.y); + case Scale::Slope: + return 0.1 * range.slope; + case Scale::Curvature: + return 0.1 * range.slope / std::max(range.x, 1e-12); + case Scale::Position: + return 0.05 * range.x; + case Scale::Fraction: + return std::max(0.1 * std::abs(start), 0.01); + case Scale::Ratio: + return std::max(0.01 * std::abs(start), 1e-3); } -}; + return 0.0; +} -template <> struct Steps { - static void compute(const std::array &start, double x_range, - double y_range, double slope_scale, - std::array &steps) { - steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range); - steps[1] = 0.1 * slope_scale; - steps[2] = 0.1 * slope_scale / std::max(x_range, 1e-12); - } -}; - -template <> struct Steps { - static void compute(const std::array &start, double x_range, - double y_range, double /*slope_scale*/, - std::array &steps) { - steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range); - steps[1] = 0.05 * x_range; - steps[2] = 0.05 * x_range; - } -}; - -template <> struct Steps { - static void compute(const std::array &start, double x_range, - double y_range, double /*slope_scale*/, - std::array &steps) { - steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range); - steps[1] = std::max(0.1 * std::abs(start[1]), 0.1 * y_range); - steps[2] = 0.05 * x_range; - steps[3] = 0.05 * x_range; - } -}; - -template <> struct Steps { - static void compute(const std::array &start, double x_range, - double y_range, double slope_scale, - std::array &steps) { - steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range); - steps[1] = 0.1 * slope_scale; - steps[2] = 0.05 * x_range; - steps[3] = 0.05 * x_range; - steps[4] = std::max(0.1 * std::abs(start[4]), 0.1 * y_range); - steps[5] = std::max(0.1 * std::abs(start[5]), 0.01); - } -}; - -template <> struct Steps { - static void compute(const std::array &start, double x_range, - double y_range, double slope_scale, - std::array &steps) { - steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range); - steps[1] = 0.1 * slope_scale; - steps[2] = 0.05 * x_range; - steps[3] = 0.05 * x_range; - steps[4] = std::max(0.1 * std::abs(start[4]), 0.1 * y_range); - steps[5] = std::max(0.1 * std::abs(start[5]), 0.01); - steps[6] = std::max(0.01 * std::abs(start[6]), 1e-3); - steps[7] = std::max(0.1 * std::abs(start[7]), 0.01); - } -}; - -template <> struct Steps { - static void compute(const std::array &start, double x_range, - double y_range, double slope_scale, - std::array &steps) { - steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range); - steps[1] = 0.1 * slope_scale; - steps[2] = 0.05 * x_range; - steps[3] = 0.05 * x_range; - steps[4] = std::max(0.1 * std::abs(start[4]), 0.1 * y_range); - steps[5] = 0.1 * slope_scale; - } -}; - -template <> struct Steps { - static void compute(const std::array &start, double x_range, - double y_range, double slope_scale, - std::array &steps) { - Steps::compute(start, x_range, y_range, - slope_scale, steps); - } -}; +/** @brief Minuit step sizes of all parameters for one pixel. */ +template +std::array +steps(const std::array &start, NDView x, + NDView y) { + const Ranges range = compute_ranges(x, y); + std::array out{}; + for (std::size_t k = 0; k < Model::npar; ++k) + out[k] = step(Scales::value[k], start[k], range); + return out; +} } // namespace aare::minuit2 diff --git a/src/fit/VariableProjection.hpp b/src/fit/VariableProjection.hpp index 0ea0cdbc..5b34c1eb 100644 --- a/src/fit/VariableProjection.hpp +++ b/src/fit/VariableProjection.hpp @@ -2,6 +2,7 @@ #pragma once #include "DampedGaussNewton.hpp" +#include "FitHelpers.hpp" #include "aare/FitModel.hpp" #include "aare/Models.hpp" #include "aare/NDView.hpp" @@ -29,15 +30,15 @@ namespace aare::detail { * is the Schur complement of the full one. Compared with * LevenbergMarquardt the iteration needs fewer evaluations, cannot be misled * by poor start values of the linear parameters, and each evaluation - * accumulates far fewer sums. The driver rejects trial points at which the - * model loses its sensitivity to a nonlinear parameter and polishes the - * converged point with one undamped step. + * accumulates far fewer sums. The few iterations stop while the damping + * still shortens the steps, so the driver polishes the converged point with + * one undamped step. * * One evaluation is a pass over the data that writes the basis functions - * and their derivatives as columns (the model's vectorisable basis_columns - * when it has one, otherwise basis_and_grad per point), a reduction to the - * normal equations of the linear parameters, the small solve, and a second - * reduction over the columns for the projected residual. + * and their derivatives as columns (the model's vectorisable + * basis_columns), a reduction to the normal equations of the linear + * parameters, the small solve, and a second reduction over the columns for + * the projected residual. * * Keep one instance per thread: the buffers are reused between pixels so * that after the first fit no memory is allocated. @@ -50,8 +51,8 @@ namespace aare::detail { * solution violates one, the fit reports failure so that the caller can * fall back to LevenbergMarquardt. * - Start values of free linear parameters are ignored. - * - max_calls counts evaluations of the projected problem, as in - * LevenbergMarquardt; 0 selects 200 + 100 * npar + 5 * npar^2. + * - max_calls counts evaluations of the projected problem, see + * evaluation_budget(). * - Errors are the Gauss-Newton estimates of the full model at the * solution, sqrt(diag((J^T J)^-1)) over the free parameters that do not * sit on a limit, the same as LevenbergMarquardt reports. @@ -60,7 +61,7 @@ template class VariableProjection { public: static_assert(model::is_separable::value, "VariableProjection needs a model with linear_par and " - "basis_and_grad"); + "basis_columns"); using Traits = model::separable_traits; static constexpr int npar = static_cast(Model::npar); static constexpr int nlin = static_cast(Traits::nlin); @@ -70,21 +71,31 @@ template class VariableProjection { static_assert(nlin >= 1, "a separable model has at least one linear " "parameter"); - using Result = MinimizerResult; // The driver iterates over the nonlinear parameters. A fully linear - // model has none and is solved without iterating. + // model has none; it gets a single fixed one, which the driver solves + // with the one evaluation of its start point. static constexpr int ndrv = nnl > 0 ? nnl : 1; static constexpr int ncol = nnl * nlin > 0 ? nnl *nlin : 1; - using Driver = DampedGaussNewton; + + // Kept with every evaluated point: the linear solution beta and the + // blocks of the Gauss-Newton Hessian [[M, Q^T], [Q, N]] of the full + // model there, which errors() needs at the final point. + struct Solution { + Vec beta{}; + SymMat M; // Phi^T W Phi + Mat Q{}; // U^T W Phi + SymMat N; // U^T W U + }; + using Driver = DampedGaussNewton; using Normal = typename Driver::Normal; /** * @brief Fit one pixel, see LevenbergMarquardt::fit for the arguments. */ - Result fit(const FitModel &model, NDView x, - NDView y, NDView s, - const std::array &start, double *par_out, - double *err_out) { + MinimizerResult fit(const FitModel &model, NDView x, + NDView y, NDView s, + const std::array &start, + double *par_out, double *err_out) { x_ = x; y_ = y; const auto n = static_cast(x.size()); @@ -100,7 +111,6 @@ template class VariableProjection { } } pvec_.assign(start.begin(), start.end()); - have_last_ = false; have_ref_ = false; nfree_lin_ = 0; for (int j = 0; j < nlin; ++j) { @@ -112,7 +122,8 @@ template class VariableProjection { Vec alpha{}; Vec lower{}; Vec upper{}; - std::array fixed{}; + std::array fixed; + fixed.fill(true); for (int k = 0; k < nnl; ++k) { const auto idx = static_cast(Traits::nonlinear_par[k]); @@ -121,58 +132,37 @@ template class VariableProjection { upper[k] = model.upper_limit(idx); fixed[k] = model.is_user_fixed(idx); } - const int max_calls = model.max_calls() > 0 - ? static_cast(model.max_calls()) - : 200 + 100 * npar + 5 * npar * npar; + typename Driver::Options opt; + opt.tolerance = model.tolerance(); + opt.max_calls = evaluation_budget(model); + opt.polish = true; - Result res; - auto fail = [&]() { - res.valid = false; - res.chi2 = 0.0; - std::fill(par_out, par_out + npar, 0.0); - if (err_out) - std::fill(err_out, err_out + npar, 0.0); - return res; + auto eval = [this](const Vec &a, Normal &out) { + return evaluate(a, out); }; + MinimizerResult res = + driver_.fit(eval, alpha, lower, upper, fixed, opt); + if (!res.valid) + return res; + const Solution &solution = driver_.normal().aux; + for (int k = 0; k < nnl; ++k) + pvec_[Traits::nonlinear_par[k]] = driver_.point()[k]; + for (int j = 0; j < nlin; ++j) + pvec_[Traits::linear_index(j)] = solution.beta[j]; + res.valid = within_linear_limits(model); + if (!res.valid) + return res; + std::copy(pvec_.begin(), pvec_.end(), par_out); + if (err_out) + errors(model, solution, err_out); + return res; + } - if constexpr (nnl == 0) { - Normal nrm; - if (!evaluate(alpha, nrm)) - return fail(); - res.calls = 1; - res.valid = true; - res.chi2 = 2.0 * nrm.F; - for (int j = 0; j < nlin; ++j) - pvec_[Traits::linear_index(j)] = - nrm.aux[static_cast(j)]; - } else { - auto eval = [this](const Vec &a, Normal &out) { - return evaluate(a, out); - }; - typename Driver::Options opt; - opt.tolerance = model.tolerance(); - opt.max_calls = max_calls; - // The same cautious start as LevenbergMarquardt: a large first - // step can still drive a width to zero. Rejecting degenerate - // trial points guards against that, and the polish step makes - // up for the few iterations, which stop while the damping still - // shortens the steps. - opt.damping0 = 0.1; - opt.reject_degenerate = true; - opt.polish = true; - res = driver_.fit(eval, alpha, lower, upper, fixed, opt); - if (!res.valid) - return fail(); - for (int k = 0; k < nnl; ++k) - pvec_[Traits::nonlinear_par[k]] = driver_.point()[k]; - for (int j = 0; j < nlin; ++j) - pvec_[Traits::linear_index(j)] = - driver_.normal().aux[static_cast(j)]; - } - - // The linear solve knows nothing about limits. A solution beyond a - // limit fails; one within rounding of a limit is moved onto it, so - // that it counts as on the limit whichever way the rounding went. + private: + // The linear solve knows nothing about limits. A solution beyond a + // limit fails; one within rounding of a limit is moved onto it, so that + // it counts as on the limit whichever way the rounding went. + bool within_linear_limits(const FitModel &model) { for (int j = 0; j < nlin; ++j) { if (lin_fixed_[j]) continue; @@ -181,21 +171,17 @@ template class VariableProjection { const double hi = model.upper_limit(k); const double tol = 1e-12 * std::max(1.0, std::abs(pvec_[k])); if (pvec_[k] < lo - tol || pvec_[k] > hi + tol) - return fail(); + return false; if (std::abs(pvec_[k] - lo) <= tol) pvec_[k] = lo; else if (std::abs(pvec_[k] - hi) <= tol) pvec_[k] = hi; } - std::copy(pvec_.begin(), pvec_.end(), par_out); - if (err_out) - errors(model, err_out); - return res; + return true; } - private: // F = chi2 / 2, gradient and Kaufman Hessian of the projected residual - // at the nonlinear point alpha, and the linear solution in out.aux. + // at the nonlinear point alpha, and its Solution in out.aux. // Returns false when the model rejects the point, or when the basis // functions are numerically linearly dependent there or one of them has // vanished on the scan points, see reduce(). @@ -204,40 +190,26 @@ template class VariableProjection { pvec_[Traits::nonlinear_par[k]] = alpha[k]; if (!Model::is_valid(pvec_)) return false; - const ssize_t n = x_.size(); - // Basis functions and their derivatives, one column per function: // phi_[j * n + i] and dphi_[(k * nlin + j) * n + i]. - if constexpr (model::has_basis_columns::value) { - Model::basis_columns(x_.data(), n, pvec_, phi_.data(), - dphi_.data()); - } else { - typename Traits::Basis phi{}; - typename Traits::BasisGrad dphi{}; - double *pc = phi_.data(); - double *dc = dphi_.data(); - const auto un = static_cast(n); - for (ssize_t i = 0; i < n; ++i) { - Model::basis_and_grad(x_[i], pvec_, phi, dphi); - const auto ui = static_cast(i); - for (int j = 0; j < nlin; ++j) - pc[static_cast(j) * un + ui] = - phi[static_cast(j)]; - for (int k = 0; k < nnl; ++k) - for (int j = 0; j < nlin; ++j) - dc[static_cast(k * nlin + j) * un + ui] = - dphi[static_cast(k)] - [static_cast(j)]; - } - } - return weighted_ ? reduce(alpha, out) : reduce(alpha, out); + Model::basis_columns(x_.data(), x_.size(), pvec_, phi_.data(), + dphi_.data()); + return weighted_ ? reduce(out) : reduce(out); } // The reductions over the columns: normal equations and solve of the // linear parameters, then the projected residual. All loops over the // linear and nonlinear parameters have compile-time bounds; the free // linear block is packed out of the full sums afterwards. - template bool reduce(const Vec &alpha, Normal &out) { + // + // The two passes over the data are written out, and the factorisation + // is a member rather than a local, because the code gcc generates for + // the second pass depends on the rest of this function: with a local + // Cholesky, with SymMat::rank1_update in the first pass or with + // compress() for the block A, the weighted S-curve fit executes 6 to + // 10 % more instructions. Check BM_FitScurveWeightedVarPro with + // perf stat when changing it. + template bool reduce(Normal &out) { const ssize_t n = x_.size(); const auto un = static_cast(n); std::array pc; @@ -251,7 +223,7 @@ template class VariableProjection { // Pass 1: normal equations M = Phi^T W Phi, t = Phi^T W y of all // linear parameters. - Mat M{}; + SymMat M; Vec t{}; for (ssize_t i = 0; i < n; ++i) { const double wi = Weighted ? w2[i] : 1.0; @@ -263,23 +235,15 @@ template class VariableProjection { const double fa = f[a] * wi; t[a] += fa * yi; for (int b = a; b < nlin; ++b) - M[a][b] += fa * f[b]; + M.upper[a][b] += fa * f[b]; } } // Free linear parameters: their block of M, with the fixed ones - // moved to the right-hand side, v_a = t_a - sum_fixed c_j M_aj. - // The lower triangle of M is filled in so that the run-time indices - // below need no branch: gcc's range analysis derives an impossible - // subscript from the branch not taken and warns about it. The loops - // run over the compile-time bound and write zeros past nfree_lin_ - // for the same reason, gcc otherwise reports v and qf as possibly - // uninitialised when they are passed to cholesky_substitute. Both - // are false positives; neither change affects the results. - for (int a = 1; a < nlin; ++a) - for (int b = 0; b < a; ++b) - M[a][b] = M[b][a]; - Mat A{}; + // moved to the right-hand side, v_a = t_a - sum_fixed c_j M_aj. The + // loop runs over the compile-time bound and writes zeros past + // nfree_lin_; gcc otherwise reports v as possibly uninitialised. + SymMat A; Vec v{}; for (int a = 0; a < nlin; ++a) { if (a >= nfree_lin_) { @@ -290,10 +254,10 @@ template class VariableProjection { double va = t[ja]; for (int j = 0; j < nlin; ++j) if (lin_fixed_[j]) - va -= pvec_[Traits::linear_index(j)] * M[ja][j]; + va -= pvec_[Traits::linear_index(j)] * M(ja, j); v[a] = va; - for (int b = 0; b < nfree_lin_; ++b) - A[a][b] = M[ja][free_lin_[b]]; + for (int b = a; b < nfree_lin_; ++b) + A.upper[a][b] = M(ja, free_lin_[b]); } // A free basis function that has all but vanished on the scan points // (a Gaussian narrower than the point spacing, centred between two @@ -302,28 +266,24 @@ template class VariableProjection { // with that at the pixel's first evaluation, the start point. if (!have_ref_) { for (int j = 0; j < nlin; ++j) - ref_norm2_[j] = M[j][j]; + ref_norm2_[j] = M(j, j); have_ref_ = true; } for (int a = 0; a < nfree_lin_; ++a) { const int j = free_lin_[a]; - if (!(M[j][j] >= vanished_basis * ref_norm2_[j])) + if (!(M(j, j) >= vanished_basis * ref_norm2_[j])) return false; } - Mat L{}; - if (!cholesky_factor(nfree_lin_, A, L)) + if (!cholesky_.factor(nfree_lin_, A)) return false; // Nearly collinear basis functions: a pivot of the factorisation far // below its diagonal element means that the linear solution loses // about log10(diagonal / pivot) digits, beyond which it is noise. - bool refine = false; - for (int a = 0; a < nfree_lin_; ++a) { - if (!(L[a][a] * L[a][a] >= collinear_basis * A[a][a])) - return false; - refine = refine || L[a][a] * L[a][a] < refine_basis * A[a][a]; - } + if (!cholesky_.pivots_at_least(collinear_basis)) + return false; + const bool refine = !cholesky_.pivots_at_least(refine_basis); Vec beta_free{}; - cholesky_substitute(nfree_lin_, L, v, beta_free); + cholesky_.solve(v, beta_free); for (int a = 0; a < nfree_lin_; ++a) pvec_[Traits::linear_index(free_lin_[a])] = beta_free[a]; // Short of that, the digits lost to the normal equations are @@ -346,25 +306,22 @@ template class VariableProjection { s[j] += pc[j][i] * rw; } Vec sf{}; - for (int a = 0; a < nlin; ++a) - sf[a] = a < nfree_lin_ ? s[free_lin_[a]] : 0.0; + compress(nlin, s, lin_fixed_, sf); Vec delta{}; - cholesky_substitute(nfree_lin_, L, sf, delta); + cholesky_.solve(sf, delta); for (int a = 0; a < nfree_lin_; ++a) pvec_[Traits::linear_index(free_lin_[a])] += delta[a]; } Vec beta; - for (int j = 0; j < nlin; ++j) { + for (int j = 0; j < nlin; ++j) beta[j] = pvec_[Traits::linear_index(j)]; - out.aux[static_cast(j)] = beta[j]; - } // Pass 2: residual r = y - Phi beta and u_k = df/dalpha_k at fixed // beta; F, b = U^T W r, N = U^T W U and Q = Phi^T W U. double F = 0.0; Vec b{}; - Mat N{}; - std::array, ndrv> Q{}; + SymMat N; + Mat Q{}; for (ssize_t i = 0; i < n; ++i) { const double wi = Weighted ? w2[i] : 1.0; Vec f; @@ -383,44 +340,33 @@ template class VariableProjection { const double uw = u[k] * wi; b[k] += uw * r; for (int l = k; l < nnl; ++l) - N[k][l] += uw * u[l]; + N.upper[k][l] += uw * u[l]; for (int j = 0; j < nlin; ++j) Q[k][j] += uw * f[j]; } } - // The full Gauss-Newton Hessian of the model at this point is - // [[M, Q^T], [Q, N]]; kept for the errors of the final point. - for (int k = 0; k < nnl; ++k) - last_alpha_[k] = alpha[k]; - for (int a = 0; a < nlin; ++a) - for (int c = a; c < nlin; ++c) - last_M_[a][c] = M[a][c]; - for (int k = 0; k < nnl; ++k) { - for (int j = 0; j < nlin; ++j) - last_Q_[k][j] = Q[k][j]; - for (int l = k; l < nnl; ++l) - last_N_[k][l] = N[k][l]; - } - have_last_ = true; + out.aux.beta = beta; + out.aux.M = M; + out.aux.Q = Q; + out.aux.N = N; // Kaufman's Jacobian: U projected out of the span of the free basis, // H = N - Q_free^T A^-1 Q_free. - std::array, ndrv> qf{}; - std::array, ndrv> z{}; + Mat qf{}; + Mat z{}; for (int k = 0; k < nnl; ++k) { - for (int a = 0; a < nlin; ++a) - qf[k][a] = a < nfree_lin_ ? Q[k][free_lin_[a]] : 0.0; - cholesky_substitute(nfree_lin_, L, qf[k], z[k]); + compress(nlin, Q[k], lin_fixed_, qf[k]); + cholesky_.solve(qf[k], z[k]); } out.F = F; for (int k = 0; k < nnl; ++k) { out.g[k] = -b[k]; for (int l = k; l < nnl; ++l) { - double h = N[k][l]; + double h = N(k, l); for (int a = 0; a < nfree_lin_; ++a) h -= qf[k][a] * z[l][a]; - out.H[k][l] = h; + out.H(k, l) = h; } } return true; @@ -428,41 +374,25 @@ template class VariableProjection { // Gauss-Newton errors from the Jacobian of the full model at the // solution, the estimate LevenbergMarquardt reports as well. The - // Jacobian's normal matrix is assembled from the sums of the last - // evaluation when that was the final point, which it is unless the last - // trial was rejected; otherwise one pass over the data computes it. - void errors(const FitModel &model, double *err_out) { - Mat H{}; - bool at_solution = have_last_; - for (int k = 0; k < nnl; ++k) - at_solution = at_solution && - pvec_[Traits::nonlinear_par[k]] == last_alpha_[k]; - if (at_solution) { - for (int a = 0; a < nlin; ++a) - for (int c = a; c < nlin; ++c) - H[Traits::linear_index(a)][Traits::linear_index(c)] = - last_M_[a][c]; - for (int k = 0; k < nnl; ++k) { - const auto m = static_cast(Traits::nonlinear_par[k]); - for (int j = 0; j < nlin; ++j) { - const auto i = static_cast(Traits::linear_index(j)); - H[std::min(i, m)][std::max(i, m)] = last_Q_[k][j]; - } - for (int l = k; l < nnl; ++l) - H[m][Traits::nonlinear_par[l]] = last_N_[k][l]; - } - } else { - const ssize_t n = x_.size(); - std::array g{}; - double f = 0.0; - for (ssize_t i = 0; i < n; ++i) { - Model::eval_and_grad(x_[i], pvec_, f, g); - const double wi = - weighted_ ? w2_[static_cast(i)] : 1.0; - for (int a = 0; a < npar; ++a) - for (int c = a; c < npar; ++c) - H[a][c] += g[a] * g[c] * wi; - } + // Jacobian's normal matrix is assembled from the sums that the + // evaluation of the final point left in its Solution. + void errors(const FitModel &model, const Solution &solution, + double *err_out) const { + SymMat H; + const auto lin = [](int j) { + return static_cast(Traits::linear_index(j)); + }; + const auto nonlin = [](int k) { + return static_cast(Traits::nonlinear_par[k]); + }; + for (int a = 0; a < nlin; ++a) + for (int c = a; c < nlin; ++c) + H(lin(a), lin(c)) = solution.M(a, c); + for (int k = 0; k < nnl; ++k) { + for (int j = 0; j < nlin; ++j) + H(lin(j), nonlin(k)) = solution.Q[k][j]; + for (int l = k; l < nnl; ++l) + H(nonlin(k), nonlin(l)) = solution.N(k, l); } // Fixed parameters and parameters that end on a limit have no error. // The driver knows the nonlinear ones; the linear solve can land a @@ -474,23 +404,15 @@ template class VariableProjection { skip[k] = model.is_user_fixed(idx) || v <= model.lower_limit(idx) || v >= model.upper_limit(idx); } - if constexpr (nnl > 0) { - for (int k = 0; k < nnl; ++k) - skip[Traits::nonlinear_par[k]] = - skip[Traits::nonlinear_par[k]] || driver_.skipped(k); - } + for (int k = 0; k < nnl; ++k) + skip[Traits::nonlinear_par[k]] = + skip[Traits::nonlinear_par[k]] || driver_.skipped(k); gauss_newton_errors(H, skip, err_out); } NDView x_, y_; Driver driver_; bool weighted_ = false; - // Sums of the last evaluation, see errors(). - bool have_last_ = false; - Vec last_alpha_{}; - Mat last_M_{}; - std::array, ndrv> last_Q_{}; - Mat last_N_{}; // Squared norms of the basis functions at the start point, see reduce(). bool have_ref_ = false; Vec ref_norm2_{}; @@ -499,6 +421,7 @@ template class VariableProjection { // Pivot ratio below which the linear solution is refined, about four // lost digits. static constexpr double refine_basis = 1e-4; + Cholesky cholesky_; // of the free linear block, see reduce() int nfree_lin_ = 0; std::array free_lin_{}; std::array lin_fixed_{};