diff --git a/AGENTS.md b/AGENTS.md index c4918858..5646aaae 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -15,7 +15,8 @@ footprint. - `include/aare/` contains the public C++ API and header-defined templates. - `src/` contains C++ implementations and colocated Catch2 tests named - `*.test.cpp`. + `*.test.cpp`, grouped into `src/fit/`, `src/hist/` and `src/utils/` where + a component has several private files. - `tests/` contains the C++ test executable, configuration, and shared helpers. - `python/src/` contains pybind11 bindings and module registration. - `python/aare/` contains Python facades, convenience APIs, and public exports. @@ -145,10 +146,26 @@ lives in `python/aare/transform.py` on top of C++ generators. `Fit.hpp`, `FitModel`, and `Models.hpp` (`Gaussian`, `Pol1`, `Pol2`, `RisingScurve`, `FallingScurve`, `GaussianErfcPlateau`, -`GaussianChargeSharing`, ...) wrap Minuit2. Template bodies and explicit -instantiations live in `src/Fit.cpp`; `FitModelImpl` is a pimpl so Minuit2 -headers never leak into the public API. Minuit2 is a private, -`BUILD_INTERFACE`-only dependency of `aare_core`. +`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 +(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`). +`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 +every configuration through all four minimizers and check that they +agree; run them when changing any solver. ### Python layer diff --git a/CMakeLists.txt b/CMakeLists.txt index e75ea5d4..8a1aa4eb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -400,7 +400,8 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/fit/Fit.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/fit/FitMinuit2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp @@ -461,7 +462,7 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.test.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/fit/Fit.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolation.test.cpp diff --git a/RELEASE.md b/RELEASE.md index 7fa0910e..4c8b7ffb 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -32,6 +32,48 @@ - Added string representator in python for Cluster and Eta - Added roi slice method in python for easy slicing of numpy arrays ``array[roi.slice()]``. - added context manager for ``aare.RawMasterFile`` +- Added a configurable minimizer for the fit models, selected with the + ``minimizer`` constructor argument or property in Python, or + ``FitModel::SetMinimizer`` in C++. ``Minimizer.Migrad`` remains the default. + ``Minimizer.Fumili`` is Minuit2's Gauss-Newton minimizer, which uses the + analytic model derivatives and needs far fewer function evaluations. + ``Minimizer.LevenbergMarquardt`` is a built-in, dependency-free + Levenberg-Marquardt solver with the analytic Jacobian; it reflects steps at + parameter limits and fits data cubes without per-pixel allocations. A fit + that stalls away from a minimum, or whose width collapses below the spacing + of the scan points, is reported as failed (all zeros) rather than returned + as a result, so noisy or signal-free pixels may fail. With + ``compute_errors``, Fumili and LevenbergMarquardt report parameter errors + from their linearised covariance instead of running Hesse. Pixels for which + Fumili does not reach a valid minimum, which Minuit2's implementation cannot + once a two-sided limit becomes active, are refitted with Migrad. +- ``Minimizer.LevenbergMarquardt`` is about twice as fast for every model: it + evaluates the Jacobian at the trial point instead of re-evaluating accepted + points, accumulates the normal equations with fixed-size loops, needs a + single factorisation per iteration, and the erf-based models compute one + exponential per point instead of two. Results are unchanged, except that + ``GaussianChargeSharing`` and ``GaussianChargeSharingKb`` now evaluate the + erf from the Gaussian's exponential, a rounding-level change. One iteration + costs one model evaluation of the ``max_calls`` budget, so pixels that used + 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 + 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 + linear parameters. When the basis functions are nearly collinear, as for a + polynomial on a narrow range far from x = 0, one more pass over the data + refines the linear solution to the precision of the other minimizers. + 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. - ``NDArray``/``NDView`` expressions now support scalar operands (``2 * a + b / 4``) and can be assigned to an existing ``NDArray``, reusing its buffer. @@ -48,6 +90,19 @@ ``NDView`` overloads are unchanged. - ``NDArray`` ``+ - * /`` with a scalar now returns a lazy expression instead of an ``NDArray`` and no longer converts the scalar to the element type. +- ``FitModel`` is now plain data without a pimpl (C++ users need to rebuild); + it gained ``strategy()``, ``lower_limit()``, ``upper_limit()`` and + ``value()`` accessors and lost ``impl()``. Bad parameter indices raise + ``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. +- 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`` + while the Minuit2 minimizers use a wide two-sided range as before; with + ``compute_errors``, fixed parameters and parameters ending on a limit report + an error of 0 (fixed parameters previously reported 1.0 with Minuit2). + - ``FastPedestal`` variance is now a private ``double`` intermediate. Removed the C++ ``variance()``/``variance_unchecked()`` APIs and Python ``var()``. Standard deviation is calculated before conversion to the output type, @@ -93,6 +148,16 @@ ### Bugfixes: +- Fitting a 3D data cube in Python without ``y_err`` now returns ``par_err`` + when ``compute_errors`` is set, as fitting a single pixel already did. +- Samples with a zero error are ignored by every minimizer, as documented, + also when they hold NaN. They are now removed before the start values are + estimated; previously a masked NaN at a sample used by the estimate, such + as the first or last point of ``Pol1``, made even Migrad fail. A pixel with + fewer remaining samples than free parameters is reported as failed. +- Fitting a single pixel checks that ``x``, ``y`` and ``y_err`` have the same + size and raises ``RuntimeError`` otherwise. Previously a shorter ``y`` or + ``y_err`` was read past its end. - Setting nSigma on ``ClusterFinder`` and ``ClusterFinderMT`` before the pedestal is ready no longer raises an error. The new value is used when the threshold is initialized. Previously ``ClusterFinderMT.set_nSigma()`` also diff --git a/benchmarks/fit_benchmark.cpp b/benchmarks/fit_benchmark.cpp index 626bde8d..8e522f80 100644 --- a/benchmarks/fit_benchmark.cpp +++ b/benchmarks/fit_benchmark.cpp @@ -3,6 +3,7 @@ #include "aare/FitModel.hpp" #include "aare/Models.hpp" +#include #include #include #include @@ -75,52 +76,31 @@ static void report_accuracy(benchmark::State &state, const TestCase &tc, // Benchmarks // ---------- -// Minuit2, analytic gradient (no Hesse) -static void BM_FitGausMinuitGrad(benchmark::State &state) { - const auto &tc = get_test_cases()[state.range(0)]; - auto data = generate_gaussian_data(tc); - auto xv = data.x.view(); - auto yv = data.y.view(); - - const auto model = aare::FitModel( - /*strategy = */ 0, - /*max_calls = */ 500, // increase for noisy signals - /*tolerance = */ 0.5, - /*compute_errors = */ false); - - aare::NDArray result; - for (auto _ : state) { - result = aare::fit_pixel(model, xv, yv); - benchmark::DoNotOptimize(result.data()); - } - - report_accuracy(state, tc, result); - state.SetLabel(tc.name); -} - -// Minuit2, analytic gradient + Hesse -static void BM_FitGausMinuitGradHesse(benchmark::State &state) { +static void run_gaussian_fit(benchmark::State &state, aare::Minimizer minimizer, + bool compute_errors) { const auto &tc = get_test_cases()[state.range(0)]; auto data = generate_gaussian_data(tc); auto xv = data.x.view(); auto yv = data.y.view(); auto ev = data.y_err.view(); - const auto model = aare::FitModel( - 0, 500, 0.5, true); // compute_errors = true -> Runs Hesse and provides - // errors on fitted params + const aare::FitModel model( + /*strategy = */ 0, + /*max_calls = */ 500, // increase for noisy signals + /*tolerance = */ 0.5, compute_errors, minimizer); aare::NDArray result; for (auto _ : state) { - result = aare::fit_pixel(model, xv, yv, ev); + result = compute_errors + ? aare::fit_pixel(model, xv, yv, ev) + : aare::fit_pixel(model, xv, yv); benchmark::DoNotOptimize(result.data()); } - // result has 6 elements: [A, mu, sig, err_A, err_mu, err_sig] report_accuracy(state, tc, result); - // Also report Hesse uncertainties - if (result.size() >= 6) { + // With errors the result is [A, mu, sig, err_A, err_mu, err_sig, chi2] + if (compute_errors && result.size() >= 6) { state.counters["errA"] = result(3); state.counters["errMu"] = result(4); state.counters["errSig"] = result(5); @@ -128,12 +108,372 @@ static void BM_FitGausMinuitGradHesse(benchmark::State &state) { state.SetLabel(tc.name); } -BENCHMARK(BM_FitGausMinuitGrad) +// Migrad, analytic gradient (no Hesse) +static void BM_FitGausMigrad(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Migrad, false); +} + +// Migrad, analytic gradient + Hesse for the parameter errors +static void BM_FitGausMigradHesse(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Migrad, true); +} + +// Fumili, analytic gradient and linearised Hessian +static void BM_FitGausFumili(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Fumili, false); +} + +// Fumili, parameter errors from the linearised covariance +static void BM_FitGausFumiliErrors(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Fumili, true); +} + +// Built-in Levenberg-Marquardt, analytic Jacobian +static void BM_FitGausLM(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::LevenbergMarquardt, false); +} + +// Built-in Levenberg-Marquardt, Gauss-Newton parameter errors +static void BM_FitGausLMErrors(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::LevenbergMarquardt, true); +} + +// Built-in variable projection: A solved exactly, iteration over mu and sigma +static void BM_FitGausVarPro(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::VarPro, false); +} + +static void BM_FitGausVarProErrors(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::VarPro, true); +} + +// Gaussian with sigma fixed: exercises the fixed-parameter path of the +// minimizers on the Moderate_noise case. +static void run_gaussian_fixed_sigma_fit(benchmark::State &state, + aare::Minimizer minimizer) { + const auto &tc = get_test_cases()[1]; + auto data = generate_gaussian_data(tc); + auto xv = data.x.view(); + auto yv = data.y.view(); + + aare::FitModel model(0, 500, 0.5, false, minimizer); + model.FixParameter(2, tc.true_sig); + + aare::NDArray result; + for (auto _ : state) { + result = aare::fit_pixel(model, xv, yv); + benchmark::DoNotOptimize(result.data()); + } + report_accuracy(state, tc, result); +} + +static void BM_FitGausFixedSigmaMigrad(benchmark::State &state) { + run_gaussian_fixed_sigma_fit(state, aare::Minimizer::Migrad); +} +static void BM_FitGausFixedSigmaFumili(benchmark::State &state) { + run_gaussian_fixed_sigma_fit(state, aare::Minimizer::Fumili); +} +static void BM_FitGausFixedSigmaLM(benchmark::State &state) { + run_gaussian_fixed_sigma_fit(state, aare::Minimizer::LevenbergMarquardt); +} +static void BM_FitGausFixedSigmaVarPro(benchmark::State &state) { + run_gaussian_fixed_sigma_fit(state, aare::Minimizer::VarPro); +} + +// ---------------------------------------------------------------- +// Rising S-curve (6 parameters), typical for threshold scans +// ---------------------------------------------------------------- +static constexpr ssize_t N_SCURVE_POINTS = 100; + +static void run_scurve_fit(benchmark::State &state, aare::Minimizer minimizer, + bool weighted = false) { + const std::vector truth = {5.0, 0.1, 50.0, 4.0, 200.0, 0.5}; + const double noise_sigma = 0.02 * truth[4]; + + aare::NDArray x({N_SCURVE_POINTS}); + aare::NDArray y({N_SCURVE_POINTS}); + aare::NDArray y_err({N_SCURVE_POINTS}, noise_sigma); + std::mt19937 rng(SEED); + std::normal_distribution noise(0.0, noise_sigma); + for (ssize_t i = 0; i < N_SCURVE_POINTS; ++i) { + x[i] = static_cast(i); + y[i] = aare::model::RisingScurve::eval(x[i], truth) + noise(rng); + } + auto xv = x.view(); + auto yv = y.view(); + auto ev = y_err.view(); + + const aare::FitModel model(0, 500, 0.5, false, + minimizer); + + aare::NDArray result; + for (auto _ : state) { + result = + weighted + ? aare::fit_pixel(model, xv, yv, ev) + : aare::fit_pixel(model, xv, yv); + benchmark::DoNotOptimize(result.data()); + } + + state.counters["dMu"] = result(2) - truth[2]; + state.counters["dSig"] = result(3) - truth[3]; + state.counters["dA"] = result(4) - truth[4]; + state.counters["chi2"] = result(6); +} + +static void BM_FitScurveMigrad(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::Migrad); +} + +static void BM_FitScurveFumili(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::Fumili); +} + +static void BM_FitScurveLM(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::LevenbergMarquardt); +} + +static void BM_FitScurveVarPro(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::VarPro); +} + +// Weighted fits take the per-point uncertainties into the residuals. +static void BM_FitScurveWeightedMigrad(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::Migrad, true); +} +static void BM_FitScurveWeightedFumili(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::Fumili, true); +} +static void BM_FitScurveWeightedLM(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::LevenbergMarquardt, true); +} +static void BM_FitScurveWeightedVarPro(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::VarPro, true); +} + +// ---------------------------------------------------------------- +// Pol2 (3 parameters, cheap model): exposes the solver overhead per +// iteration rather than the model evaluation. +// ---------------------------------------------------------------- +static void run_pol2_fit(benchmark::State &state, aare::Minimizer minimizer) { + constexpr ssize_t n = 51; + const std::vector truth = {12.0, -0.8, 0.05}; + + aare::NDArray x({n}); + aare::NDArray y({n}); + std::mt19937 rng(SEED); + std::normal_distribution noise(0.0, 2.0); + for (ssize_t i = 0; i < n; ++i) { + x[i] = static_cast(i); + y[i] = aare::model::Pol2::eval(x[i], truth) + noise(rng); + } + auto xv = x.view(); + auto yv = y.view(); + + const aare::FitModel model(0, 500, 0.5, false, + minimizer); + + aare::NDArray result; + for (auto _ : state) { + result = aare::fit_pixel(model, xv, yv); + benchmark::DoNotOptimize(result.data()); + } + state.counters["dp2"] = result(2) - truth[2]; +} + +static void BM_FitPol2Migrad(benchmark::State &state) { + run_pol2_fit(state, aare::Minimizer::Migrad); +} +static void BM_FitPol2Fumili(benchmark::State &state) { + run_pol2_fit(state, aare::Minimizer::Fumili); +} +static void BM_FitPol2LM(benchmark::State &state) { + run_pol2_fit(state, aare::Minimizer::LevenbergMarquardt); +} +static void BM_FitPol2VarPro(benchmark::State &state) { + run_pol2_fit(state, aare::Minimizer::VarPro); +} + +// ---------------------------------------------------------------- +// GaussianChargeSharingKb (8 parameters, two exponentials and two erfs per +// point): the most expensive model per evaluation. +// ---------------------------------------------------------------- +static void run_charge_sharing_kb_fit(benchmark::State &state, + aare::Minimizer minimizer) { + constexpr ssize_t n = 200; + const std::vector truth = {20.0, 0.05, 110.0, 6.0, + 1500.0, 0.1, 1.1, 0.1}; + + aare::NDArray x({n}); + aare::NDArray y({n}); + std::mt19937 rng(SEED); + std::normal_distribution noise(0.0, 1.0); + for (ssize_t i = 0; i < n; ++i) { + x[i] = static_cast(i); + const double clean = + aare::model::GaussianChargeSharingKb::eval(x[i], truth); + y[i] = clean + std::sqrt(std::max(clean, 1.0)) * noise(rng); + } + auto xv = x.view(); + auto yv = y.view(); + + const aare::FitModel model( + 0, 500, 0.5, false, minimizer); + + aare::NDArray result; + for (auto _ : state) { + result = aare::fit_pixel(model, + xv, yv); + benchmark::DoNotOptimize(result.data()); + } + state.counters["dMu"] = result(2) - truth[2]; + state.counters["dN"] = result(4) - truth[4]; + state.counters["chi2"] = result(8); +} + +static void BM_FitChargeSharingKbMigrad(benchmark::State &state) { + run_charge_sharing_kb_fit(state, aare::Minimizer::Migrad); +} +static void BM_FitChargeSharingKbFumili(benchmark::State &state) { + run_charge_sharing_kb_fit(state, aare::Minimizer::Fumili); +} +static void BM_FitChargeSharingKbLM(benchmark::State &state) { + run_charge_sharing_kb_fit(state, aare::Minimizer::LevenbergMarquardt); +} +static void BM_FitChargeSharingKbVarPro(benchmark::State &state) { + run_charge_sharing_kb_fit(state, aare::Minimizer::VarPro); +} + +// ---------------------------------------------------------------- +// Data cube through fit_3d on one thread: exposes the per-pixel overhead +// of each minimizer on top of the minimisation itself. +// ---------------------------------------------------------------- +static constexpr ssize_t CUBE_ROWS = 32; +static constexpr ssize_t CUBE_COLS = 32; + +static void run_cube_fit(benchmark::State &state, aare::Minimizer minimizer, + bool compute_errors) { + const auto &tc = get_test_cases()[1]; // Moderate_noise + const double x_min = tc.true_mu - 5.0 * tc.true_sig; + const double x_max = tc.true_mu + 5.0 * tc.true_sig; + const double dx = (x_max - x_min) / (N_POINTS - 1); + const double noise_sigma = tc.noise_frac * tc.true_A; + + aare::NDArray x({N_POINTS}); + aare::NDArray y({CUBE_ROWS, CUBE_COLS, N_POINTS}); + aare::NDArray y_err({CUBE_ROWS, CUBE_COLS, N_POINTS}, + noise_sigma); + std::mt19937 rng(SEED); + std::normal_distribution noise(0.0, noise_sigma); + for (ssize_t i = 0; i < N_POINTS; ++i) + x[i] = x_min + i * dx; + for (ssize_t row = 0; row < CUBE_ROWS; ++row) { + for (ssize_t col = 0; col < CUBE_COLS; ++col) { + for (ssize_t i = 0; i < N_POINTS; ++i) { + const double clean = + tc.true_A * std::exp(-std::pow(x[i] - tc.true_mu, 2) / + (2.0 * std::pow(tc.true_sig, 2))); + y(row, col, i) = clean + noise(rng); + } + } + } + + aare::NDArray par({CUBE_ROWS, CUBE_COLS, 3}); + aare::NDArray err({CUBE_ROWS, CUBE_COLS, 3}); + aare::NDArray chi2({CUBE_ROWS, CUBE_COLS}); + const aare::FitModel model( + 0, 500, 0.5, compute_errors, minimizer); + + for (auto _ : state) { + if (compute_errors) { + aare::fit_3d(model, x.view(), y.view(), y_err.view(), par.view(), + err.view(), chi2.view(), 1); + } else { + aare::fit_3d(model, x.view(), y.view(), aare::NDView{}, + par.view(), aare::NDView{}, chi2.view(), 1); + } + benchmark::DoNotOptimize(par.data()); + } + + state.SetItemsProcessed(state.iterations() * CUBE_ROWS * CUBE_COLS); + state.counters["failed"] = + static_cast(std::count(chi2.begin(), chi2.end(), 0.0)); +} + +static void BM_FitCubeMigrad(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::Migrad, false); +} +static void BM_FitCubeFumili(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::Fumili, false); +} +static void BM_FitCubeLM(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::LevenbergMarquardt, false); +} +static void BM_FitCubeMigradErrors(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::Migrad, true); +} +static void BM_FitCubeFumiliErrors(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::Fumili, true); +} +static void BM_FitCubeLMErrors(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::LevenbergMarquardt, true); +} +static void BM_FitCubeVarPro(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::VarPro, false); +} +static void BM_FitCubeVarProErrors(benchmark::State &state) { + run_cube_fit(state, aare::Minimizer::VarPro, true); +} + +BENCHMARK(BM_FitGausMigrad)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFumili)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausLM)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausVarPro)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_FitGausMigradHesse) + ->DenseRange(0, 5) + ->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFumiliErrors) + ->DenseRange(0, 5) + ->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausLMErrors)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausVarProErrors) ->DenseRange(0, 5) ->Unit(benchmark::kMicrosecond); -BENCHMARK(BM_FitGausMinuitGradHesse) - ->DenseRange(0, 5) - ->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFixedSigmaMigrad)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFixedSigmaFumili)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFixedSigmaLM)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFixedSigmaVarPro)->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_FitScurveMigrad)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveFumili)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveLM)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveVarPro)->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_FitScurveWeightedMigrad)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveWeightedFumili)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveWeightedLM)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveWeightedVarPro)->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_FitPol2Migrad)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitPol2Fumili)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitPol2LM)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitPol2VarPro)->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_FitChargeSharingKbMigrad)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitChargeSharingKbFumili)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitChargeSharingKbLM)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitChargeSharingKbVarPro)->Unit(benchmark::kMicrosecond); + +// fit_3d runs in worker threads, so measure wall time. +BENCHMARK(BM_FitCubeMigrad)->Unit(benchmark::kMillisecond)->UseRealTime(); +BENCHMARK(BM_FitCubeFumili)->Unit(benchmark::kMillisecond)->UseRealTime(); +BENCHMARK(BM_FitCubeLM)->Unit(benchmark::kMillisecond)->UseRealTime(); +BENCHMARK(BM_FitCubeMigradErrors)->Unit(benchmark::kMillisecond)->UseRealTime(); +BENCHMARK(BM_FitCubeFumiliErrors)->Unit(benchmark::kMillisecond)->UseRealTime(); +BENCHMARK(BM_FitCubeLMErrors)->Unit(benchmark::kMillisecond)->UseRealTime(); +BENCHMARK(BM_FitCubeVarPro)->Unit(benchmark::kMillisecond)->UseRealTime(); +BENCHMARK(BM_FitCubeVarProErrors)->Unit(benchmark::kMillisecond)->UseRealTime(); BENCHMARK_MAIN(); diff --git a/docs/src/pyFit.rst b/docs/src/pyFit.rst index a8375772..142b1ca6 100644 --- a/docs/src/pyFit.rst +++ b/docs/src/pyFit.rst @@ -4,11 +4,18 @@ Fitting .. py:currentmodule:: aare Aare fits one-dimensional scans and three-dimensional pixel data with -Minuit2. Create a model object and call its :meth:`fit` method:: +Minuit2's Migrad and Fumili or with the built-in Levenberg-Marquardt and +variable projection solvers. Create a model object and call its :meth:`fit` +method:: model = Gaussian(compute_errors=True) result = model.fit(x, y, y_err) +``y_err`` holds the standard deviation of each point and weights the +residuals; points with ``y_err == 0`` are left out of the fit, which skips +empty histogram bins when ``y_err = sqrt(y)``. Without ``y_err`` every point +has unit weight. + The model object is also callable, which evaluates it at the supplied points:: fitted_y = model(x, result["par"]) @@ -20,12 +27,14 @@ accepts the same model objects when a functional interface is preferred. For three-dimensional data, pass an array with shape ``(rows, columns, scan_points)`` and select the worker count with -``n_threads``:: +``n_threads`` (default 4):: result = model.fit(x, image_data, image_errors, n_threads=8) -The result dictionary contains ``par`` and ``chi2``. It also contains -``par_err`` when ``compute_errors`` is enabled. +The result dictionary contains ``par`` and ``chi2``, with shapes ``(npar,)`` +and ``(1,)`` for a scan or ``(rows, columns, npar)`` and ``(rows, columns)`` +for a cube. It also contains ``par_err`` when ``compute_errors`` is enabled. +A pixel that fails to converge returns zeros for all three. .. note:: @@ -37,4 +46,108 @@ The result dictionary contains ``par`` and ``chi2``. It also contains ``float64`` arrays to avoid it. A ``float64`` C-contiguous array is used in place without copying. +Choosing the minimizer +~~~~~~~~~~~~~~~~~~~~~~ + +.. note:: + + The built-in VarPro solver is up to ten times faster than Migrad. It and + the built-in Levenberg-Marquardt solver do not offer Minuit2's strategy + setting or Hesse errors. Try VarPro first and fall back to Migrad if you + run into problems. Minuit2's Fumili and the built-in Levenberg-Marquardt + solver are included for evaluation and may turn out to be redundant. + +Fits use Minuit2's Migrad by default. Three alternatives share the same +interface and are selected with the ``minimizer`` argument or property:: + + model = Gaussian(minimizer=Minimizer.LevenbergMarquardt) + model.minimizer = Minimizer.Fumili + +- ``Minimizer.Migrad``: Minuit2's variable-metric minimizer, driven by the + analytic gradient of the chi-squared. Parameter errors come from Hesse. +- ``Minimizer.Fumili``: Minuit2's Gauss-Newton minimizer, which controls its + steps with a trust region by default. It takes the gradient and the + linearised Hessian from the analytic derivatives of the model and typically + needs far fewer function evaluations than Migrad. The linearised Hessian at + the minimum also provides the parameter covariance, so no Hesse step runs. + + Parameter limits make Minuit2 map the bounded external parameters onto + unbounded internal ones. On an active limit the derivative of that mapping + vanishes, and Fumili's linearised Hessian, which contains only this first + derivative, loses the curvature along the bounded direction. The Hessian + can then become singular and the iteration fails to converge. Pixels for + which Fumili does not return a valid minimum are refitted with Migrad. + +- ``Minimizer.LevenbergMarquardt``: a built-in damped Gauss-Newton solver with + the analytic Jacobian of the model. It reflects steps at parameter limits + and reuses its buffers between pixels, so data cubes are fitted without + per-pixel allocations. Parameter errors are ``sqrt(diag((J^T W J)^-1))`` + with ``W`` the weights from ``y_err``. +- ``Minimizer.VarPro``: a built-in variable projection solver. + Every bundled model is linear in some of its parameters: every parameter of + a polynomial, the amplitude of a Gaussian, the amplitude and plateau of + ``GaussianErfcPlateau``, the background and amplitude of the charge-sharing + models, and the background, amplitude and charge-sharing term of the + S-curves. Each trial point solves those exactly and the Levenberg-Marquardt + iteration runs over the remaining nonlinear parameters only, so it needs + fewer evaluations than LevenbergMarquardt and is several times faster than + Migrad. Start values of the linear parameters are not used, so + ``SetParameter`` has no effect on them; ``FixParameter`` does. A pixel + whose linear solution violates a limit on a linear parameter, or that does + not converge, is refitted with LevenbergMarquardt with a fresh + ``max_calls`` budget. Parameter errors are the same Gauss-Newton estimates. + On data that barely constrains a model, VarPro can settle at a point where + an amplitude is driven to zero and the parameters attached to it lose + their meaning; such a pixel reports a valid fit with a poor ``chi2``. + +All four converge to the same minimum on well-behaved data; the ``[fit]`` +test suite checks that. The bundled ``fit_benchmark`` compares their speed on +Gaussian peaks, S-curves, polynomials, charge-sharing spectra and a data +cube. On ill-conditioned data (few points, many parameters, or parameters +that nearly trade off against each other) the minimizers stop at different +points of a flat chi-squared valley, and Migrad in particular can settle on a +parameter limit with a higher ``chi2`` than the others. Compare ``chi2`` +rather than the parameters when they disagree, fix or bound the poorly +constrained parameters, and provide start values for the rest. + +Parameters and errors +~~~~~~~~~~~~~~~~~~~~~ + +Start values are estimated from the data. ``SetParameter`` overrides the +estimate, ``FixParameter`` keeps a parameter at a given value, and +``SetParLimits`` restricts it to a range where ``-inf`` and ``inf`` leave a +side open (the Minuit2 minimizers turn an open side into a wide finite one); +free start values are clamped into their limits. The constructor arguments +control the minimizer: + +- ``max_calls`` (default 100) caps the work per pixel: Minuit2 function calls + for Migrad and Fumili, model evaluations for LevenbergMarquardt and VarPro, + which spend one evaluation per iteration (two when a step crosses a limit). + ``0`` selects Minuit's default max calls of ``200 + 100 * npar + 5 * npar**2``. + A pixel that does not converge within the budget returns zeros; the six- + and eight-parameter models on noisy data need several hundred calls, so + raise the default when such pixels fail. +- ``tolerance`` is the EDM tolerance in Minuit's convention. Migrad, + LevenbergMarquardt and VarPro stop when the estimated distance + to the minimum is below ``0.002 * tolerance``, Fumili below + ``1e-4 * tolerance``. The EDM is measured in units of ``chi2``, not + relative to the data: without ``y_err`` every sample counts with an error + of 1, so a fit of small values (amplitudes well below 1) meets the + tolerance while still far from the minimum. Pass ``y_err`` or lower + ``tolerance`` for such data. +- ``strategy`` is the Minuit2 strategy (0 is fast, 1 is Minuit2's default). + LevenbergMarquardt and VarPro ignore it. +- ``compute_errors`` adds ``par_err`` to the result. Fixed parameters and + parameters that end on a limit report an error of 0 with every minimizer. + Migrad's Hesse errors include the residual-weighted second derivatives of + the model, which the Gauss-Newton estimates of the other three omit; the + two agree when the model describes the data (``chi2`` per degree of freedom + near 1) and the Gauss-Newton errors come out smaller for a poorer fit. + +Unknown parameter names raise ``RuntimeError``, bad indices ``IndexError``, +and ``SetParLimits`` rejects a lower limit at or above the upper one. + .. autofunction:: fit + +.. autoclass:: Minimizer + :members: Migrad, Fumili, LevenbergMarquardt, VarPro diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index 1a2fa55f..cc7edce1 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -7,17 +7,24 @@ namespace aare { // --------------------------------------------------------------------------- -// Minuit2-based pixel fitting. +// Pixel fitting with the minimizer selected by the FitModel. // Template bodies and explicit instantiations live in src/Fit.cpp. // --------------------------------------------------------------------------- /** - * @brief Fit a single pixel's data using Minuit2. + * @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. * * User-precedence rules: - * - Fixed parameters: untouched (value and fixed flag preserved from model). - * - User-set start: value preserved, step size auto-filled. - * - Neither: both value and step size auto-filled from data. + * - Fixed parameters: kept at their value, excluded from the fit. + * - User-set start: value used as the starting point. + * - Neither: starting value estimated from the data and clamped + * into the parameter limits. * * @tparam Model Model struct (Gaussian, RisingScurve, …). * diff --git a/include/aare/FitModel.hpp b/include/aare/FitModel.hpp index a52ea9fa..f485849e 100644 --- a/include/aare/FitModel.hpp +++ b/include/aare/FitModel.hpp @@ -3,59 +3,102 @@ #include "aare/Models.hpp" #include -#include #include #include #include namespace aare { +/** + * @brief Minimizer used by FitModel. + * + * - Migrad: Minuit2's variable-metric minimizer, driven by the analytic + * gradient of the chi-squared. Parameter errors come from MnHesse. + * - Fumili: Minuit2's Gauss-Newton minimizer with a trust region. It takes + * the gradient and a linearised Hessian from the model's analytic + * derivatives and typically needs far fewer function evaluations than + * Migrad. Parameter errors come from the covariance of that linearised + * Hessian, without a Hesse step. Minuit2's Fumili cannot converge once a + * 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 + * 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 + * 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 + * start values of the linear parameters. Start values of free linear + * parameters are ignored. A limit on a linear parameter cannot be + * enforced by the linear solve: a pixel whose solution violates one, a + * pixel that does not converge, and every pixel of a model without the + * separable structure are fitted with LevenbergMarquardt instead. + * Parameter errors are the same Gauss-Newton estimates as with + * LevenbergMarquardt. + * + * With every minimizer, fixed parameters and parameters that end on a limit + * report an error of 0 and a failed fit returns zeros. `tolerance` is the + * EDM tolerance in Minuit's convention (Migrad, LevenbergMarquardt and + * VarPro stop below 0.002 * tolerance, Fumili below + * 1e-4 * tolerance). The EDM is in units of chi2, so an unweighted fit of + * small values meets it far from the minimum; pass errors or a lower + * tolerance for such data. `max_calls` bounds the work per pixel in the + * minimizer's own units (Minuit2 function calls, or model evaluations for + * LevenbergMarquardt and VarPro, which spend one evaluation per + * iteration and two when a step crosses a limit). `strategy` only affects + * the Minuit2 minimizers. + */ +enum class Minimizer { Migrad, Fumili, LevenbergMarquardt, VarPro }; + +/** + * @brief Fit configuration for one model: parameter limits, fixed + * 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. + */ template class FitModel { - - // Forward declaration only — full definition lives in src/FitModelImpl.hpp - // and is never installed. Downstream code cannot dereference this pointer. - struct FitModelImpl; // stores Minuit2 user-set parameters - - std::unique_ptr impl_; - unsigned int max_calls_; - double tolerance_; - bool compute_errors_; - - std::array user_fixed_{}; - std::array user_start_{}; - - /** @brief Safely resolve a parameter name to its index. */ - unsigned int checked_index(const std::string &name) const; - public: static constexpr std::size_t npar = Model::npar; /** * @brief Construct a fit model with sensible defaults. * - * @param strategy Minuit2 strategy level (0 = fast/gradient, 1 = - * default). - * @param max_calls Maximum FCN calls per pixel minimisation. - * @param tolerance Minuit2 EDM tolerance. - * @param compute_errors If true, run MnHesse after minimisation. + * @param strategy Minuit2 strategy level (0 = fast, 1 = Minuit2's + * default). Ignored by LevenbergMarquardt. + * @param max_calls Work budget per pixel: Minuit2 function calls, + * or model evaluations for LevenbergMarquardt. + * @param tolerance EDM tolerance in Minuit's convention, see + * aare::Minimizer. + * @param compute_errors If true, also report parameter errors: from + * MnHesse with Migrad, from the linearised + * covariance with Fumili and LevenbergMarquardt. + * Fixed parameters and parameters ending on a limit + * report 0. + * @param minimizer Minimizer, see aare::Minimizer. */ explicit FitModel(unsigned int strategy = 0, unsigned int max_calls = 100, - double tolerance = 0.5, bool compute_errors = false); + double tolerance = 0.5, bool compute_errors = false, + Minimizer minimizer = Minimizer::Migrad); - // Destructor must be defined in Fit.cpp where FitModelImpl is complete. - ~FitModel(); - FitModel(const FitModel &); - FitModel &operator=(const FitModel &); - FitModel(FitModel &&) noexcept = default; - FitModel &operator=(FitModel &&) noexcept = default; - - /** @brief Set lower and upper bounds for parameter idx.*/ + /** + * @brief Set lower and upper limits for parameter idx. + * + * Infinite values leave that side open (the Minuit2 minimizers turn an + * open side into a wide finite one, at least 1e6 away from the other). + * Throws std::runtime_error when lo >= hi and std::out_of_range for a + * bad index. + */ void SetParLimits(unsigned int idx, double lo, double hi); /** * @brief Fix parameter idx at value val. * - * Excluded from minimisation. Automatic estimates will not touch it. + * Excluded from minimisation. Automatic estimates will not touch it. */ void FixParameter(unsigned int idx, double val); @@ -76,16 +119,36 @@ template class FitModel { void SetMaxCalls(unsigned int n) { max_calls_ = n; } void SetTolerance(double t) { tolerance_ = t; } void SetComputeErrors(bool b) { compute_errors_ = b; } + void SetMinimizer(Minimizer m) { minimizer_ = m; } + unsigned int strategy() const { return strategy_; } unsigned int max_calls() const { return max_calls_; } double tolerance() const { return tolerance_; } bool compute_errors() const { return compute_errors_; } + Minimizer minimizer() const { return minimizer_; } bool is_user_fixed(unsigned int idx) const { return user_fixed_[idx]; } bool is_user_start(unsigned int idx) const { return user_start_[idx]; } + /** @brief Lower limit of parameter idx, -infinity when open. */ + double lower_limit(unsigned int idx) const { return lower_[idx]; } + /** @brief Upper limit of parameter idx, +infinity when open. */ + double upper_limit(unsigned int idx) const { return upper_[idx]; } + /** @brief User start or fixed value of parameter idx, 0 when unset. */ + double value(unsigned int idx) const { return value_[idx]; } - // Returns the internal Minuit2 state. FitModelImpl is an incomplete type - // here; only callers that include src/FitModelImpl.hpp can dereference it. - FitModelImpl *impl() const { return impl_.get(); } + private: + void check_index(unsigned int idx) const; + unsigned int checked_index(const std::string &name) const; + + unsigned int strategy_; + unsigned int max_calls_; + double tolerance_; + bool compute_errors_; + Minimizer minimizer_; + std::array lower_{}; + std::array upper_{}; + std::array value_{}; + std::array user_fixed_{}; + std::array user_start_{}; }; // Suppress implicit instantiation for all supported model types. diff --git a/include/aare/Models.hpp b/include/aare/Models.hpp index eee84e2d..da908f14 100644 --- a/include/aare/Models.hpp +++ b/include/aare/Models.hpp @@ -2,10 +2,15 @@ #pragma once #include "aare/NDView.hpp" +#include "aare/defs.hpp" #include #include #include +#include +#include +#include #include +#include #include namespace aare::model { @@ -13,10 +18,14 @@ namespace aare::model { inline constexpr double inv_sqrt2 = 0.70710678118654752440; inline constexpr double inv_sqrt_2pi = 0.39894228040143267794; -inline double fast_erf(double x) { - // Abramowitz–Stegun Handbook of Mathematical Functions - // erf approximation with max error ~1.5e-7, faster than std::erf. - +/** + * @brief erf(z) given exp(-z^2), from the Abramowitz-Stegun approximation + * 7.1.26 (Handbook of Mathematical Functions, max error about 1.5e-7). + * + * Models whose derivatives need exp(-z^2) anyway pass it in, so that a + * single exponential serves both the value and the gradient. + */ +ALWAYS_INLINE double fast_erf_from_exp(double z, double exp_minus_z2) { const double a1 = 0.254829592; const double a2 = -0.284496736; const double a3 = 1.421413741; @@ -24,22 +33,69 @@ inline double fast_erf(double x) { const double a5 = 1.061405429; const double p = 0.3275911; - const int sign = x < 0 ? -1 : 1; - x = std::abs(x); - - // 7.1.26 - const double t = 1.0 / (1.0 + p * x); + const double t = 1.0 / (1.0 + p * std::abs(z)); const double y = 1.0 - - (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t) * std::exp(-x * x); + (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t) * exp_minus_z2; - return sign * y; + return z < 0 ? -y : y; +} + +/** @brief erf approximation, faster than std::erf; see fast_erf_from_exp. */ +ALWAYS_INLINE double fast_erf(double x) { + return fast_erf_from_exp(x, std::exp(-x * x)); +} + +/** + * @brief exp(x) without a call into libm, for loops that should vectorise. + * + * Range reduction x = k ln2 + r with the integer k taken from the mantissa + * of x * log2(e) + 1.5 * 2^52, a degree-13 Taylor polynomial for exp(r) on + * |r| <= ln2 / 2, and 2^k assembled in the exponent bits, so that a loop + * calling it compiles to vector instructions on SSE2, AVX and NEON. The + * relative error is below 4e-16 for finite x; |x| is clamped to 700. The + * bulk basis functions of the models (basis_columns) use it; eval and + * eval_and_grad keep std::exp. It is ALWAYS_INLINE because gcc otherwise + * leaves it out of line in the large fitting translation unit, and a call + * per point defeats the vectorisation of the caller. + */ +ALWAYS_INLINE double fast_exp(double x) { + constexpr double log2e = 1.4426950408889634074; + constexpr double magic = 6755399441055744.0; // 1.5 * 2^52 + constexpr double ln2_hi = 0.693145751953125; + constexpr double ln2_lo = 1.42860682030941723212e-6; + x = std::clamp(x, -700.0, 700.0); + const double t = x * log2e + magic; + const double k = t - magic; + const double r = (x - k * ln2_hi) - k * ln2_lo; + double p = 1.0 / 6227020800.0; + p = p * r + 1.0 / 479001600.0; + p = p * r + 1.0 / 39916800.0; + p = p * r + 1.0 / 3628800.0; + p = p * r + 1.0 / 362880.0; + p = p * r + 1.0 / 40320.0; + p = p * r + 1.0 / 5040.0; + p = p * r + 1.0 / 720.0; + p = p * r + 1.0 / 120.0; + p = p * r + 1.0 / 24.0; + p = p * r + 1.0 / 6.0; + p = p * r + 0.5; + p = p * r + 1.0; + p = p * r + 1.0; + std::int64_t ti = 0; + std::int64_t mi = 0; + std::memcpy(&ti, &t, sizeof(ti)); + std::memcpy(&mi, &magic, sizeof(mi)); + const std::int64_t bits = ((ti - mi) + 1023) << 52; + double scale = 0.0; + std::memcpy(&scale, &bits, sizeof(scale)); + return p * scale; } /** * @brief Per-parameter metadata: name and optional default bounds. * - * Used by FitModel to generically initialise MnUserParameters. + * Used by FitModel for parameter names and default limits. * Unbounded directions use ±no_bound as sentinels. */ struct ParamInfo { @@ -51,22 +107,89 @@ struct ParamInfo { inline constexpr double no_bound = std::numeric_limits::infinity(); /** - * @brief Compute data-range statistics used by step-size estimators. + * @brief Detects a separable model, one that declares its linear parameters + * for Minimizer::VarPro. * - * Model-independent, called once per pixel. + * 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: + * + * 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 + * 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. */ -inline void compute_ranges(NDView x, NDView y, - double &x_range, double &y_range, - double &slope_scale) { - 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()); +template +struct is_separable : std::false_type {}; - x_range = std::max(*x_max - *x_min, 1e-9); - y_range = std::max(*y_max - *y_min, 1e-9); +template +struct is_separable> + : std::true_type {}; - slope_scale = std::max(y_range / x_range, 1e-9); +/** @brief Detects a separable model that provides basis_columns. */ +template +struct has_basis_columns : std::false_type {}; + +template +struct has_basis_columns> + : std::true_type {}; + +/** @brief The parameter indices not listed in linear, in ascending order. */ +template +constexpr std::array +nonlinear_indices(const std::array &linear) { + std::array out{}; + std::size_t m = 0; + for (std::size_t k = 0; k < Npar; ++k) { + bool is_linear = false; + for (std::size_t j = 0; j < Nlin; ++j) + is_linear = is_linear || linear[j] == k; + if (!is_linear) + out[m++] = k; + } + return out; } +/** @brief Index bookkeeping of a separable model, see is_separable. */ +template struct separable_traits { + static constexpr std::size_t npar = Model::npar; + static constexpr std::size_t nlin = Model::linear_par.size(); + static constexpr std::size_t nnl = npar - nlin; + static constexpr std::array nonlinear_par = + nonlinear_indices(Model::linear_par); + 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>; +}; + // _____________________________________________________________________ // // Pol1 @@ -105,6 +228,17 @@ struct Pol1 { g[1] = x; // df/dp1 } + /** @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 bool is_valid([[maybe_unused]] const std::vector &par) { return true; // always valid } @@ -120,14 +254,6 @@ struct Pol1 { return {intercept, slope}; } - - static void compute_steps(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; - } }; // _____________________________________________________________________ @@ -172,6 +298,18 @@ struct Pol2 { g[2] = x * x; // df/dp2 } + /** @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 bool is_valid([[maybe_unused]] const std::vector &par) { return true; // always valid } @@ -220,15 +358,6 @@ struct Pol2 { return {p0, p1, p2}; } - - static void compute_steps(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); - } }; // _____________________________________________________________________ @@ -288,6 +417,44 @@ struct Gaussian { sig; // df/dsigma = A*e*(x-mu)^2/sig^3 } + /** @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) { + const double mu = par[1]; + const double inv_sig = 1.0 / par[2]; + const double inv_sig2 = inv_sig * inv_sig; + double *e_col = phi; + double *dmu = dphi; + double *dsig = dphi + n; + for (ssize_t i = 0; i < n; ++i) { + const double dx = x[i] - mu; + const double e = fast_exp(-0.5 * dx * dx * inv_sig2); + e_col[i] = e; + dmu[i] = e * dx * inv_sig2; + dsig[i] = e * dx * dx * inv_sig2 * inv_sig; + } + } + /** @brief Reject degenerate sigma (zero width). */ static bool is_valid(const std::vector &par) { return par[2] != 0.0; @@ -326,18 +493,6 @@ struct Gaussian { return {A, mu, sig}; } - - /** - * @brief Data-driven Minuit step sizes. - */ - static void compute_steps(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; - } }; // _____________________________________________________________________ @@ -393,7 +548,7 @@ struct GaussianErfcPlateau { const double z = dx * inv_sqrt2 / sig; const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf(z)); + const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); return A * e + S * step; } @@ -409,7 +564,7 @@ struct GaussianErfcPlateau { const double z = dx * inv_sqrt2 / sig; const double e = std::exp(-z * z); - const double step = 0.5 * (1.0 - fast_erf(z)); + const double step = 0.5 * (1.0 - fast_erf_from_exp(z, e)); f = A * e + S * step; @@ -426,6 +581,63 @@ struct GaussianErfcPlateau { S * inv_sqrt_2pi * e * dx * inv_sig2; // df/dsigma } + /** @brief A and S are linear: f = A * e + S * step. */ + 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; + 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; + } + } + /** @brief Reject degenerate or negative width. */ static bool is_valid(const std::vector &par) { return par[3] > 0.0; @@ -490,16 +702,6 @@ struct GaussianErfcPlateau { return {A, S, mu, sig}; } - - static void compute_steps(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; - } }; // _____________________________________________________________________ @@ -556,7 +758,7 @@ struct GaussianChargeSharing { const double u = dx / sig; const double G = std::exp(-0.5 * u * u); - const double H = 0.5 * (1.0 - fast_erf(u * inv_sqrt2)); + const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); return p0 - p1 * x + N * (G + C * H); } @@ -574,7 +776,7 @@ struct GaussianChargeSharing { const double u = dx / sig; const double G = std::exp(-0.5 * u * u); - const double H = 0.5 * (1.0 - fast_erf(u * inv_sqrt2)); + const double H = 0.5 * (1.0 - fast_erf_from_exp(u * inv_sqrt2, G)); f = p0 - p1 * x + N * (G + C * H); @@ -596,6 +798,84 @@ struct GaussianChargeSharing { g[5] = N * H; } + /** + * @brief p0, p1 and N are linear: f = p0 * 1 + p1 * (-x) + N * (G + C H). + * C multiplies N, so it stays a nonlinear parameter. + */ + 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; + // d phi_0 and d phi_1 vanish for every nonlinear parameter. + std::fill(dphi, dphi + 2 * n, 0.0); + std::fill(dphi + 3 * n, dphi + 5 * n, 0.0); + std::fill(dphi + 6 * n, dphi + 8 * n, 0.0); + double *d_mu = dphi + 2 * n; + double *d_sig = dphi + 5 * n; + 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; + 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; + } + } + static bool is_valid(const std::vector &par) { return par[3] > 0.0; } @@ -647,18 +927,6 @@ struct GaussianChargeSharing { return {p0, p1, mu, sigma, N, C}; } - - static void compute_steps(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); - } }; // _____________________________________________________________________ @@ -730,8 +998,8 @@ struct GaussianChargeSharingKb { 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(u * inv_sqrt2)); - const double Hb = 0.5 * (1.0 - fast_erf(ub * inv_sqrt2)); + 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; @@ -761,8 +1029,8 @@ struct GaussianChargeSharingKb { 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(u * inv_sqrt2)); - const double Hb = 0.5 * (1.0 - fast_erf(ub * inv_sqrt2)); + 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; @@ -807,6 +1075,133 @@ struct GaussianChargeSharingKb { g[7] = N * kb; } + /** + * @brief p0, p1 and N are linear: f = p0 * 1 + p1 * (-x) + N * (ka + q kb). + * C, r and q multiply N or enter the peak shapes, so they stay nonlinear. + */ + 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) { + const double mu = par[2]; + const double inv_sig = 1.0 / par[3]; + 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; + double *shape = phi + 2 * n; + // d phi_0 and d phi_1 vanish for every nonlinear parameter; the + // derivative of phi_2 with respect to nonlinear parameter k is + // column k * 3 + 2. + for (int k = 0; k < 5; ++k) + std::fill(dphi + (k * 3) * n, dphi + (k * 3 + 2) * n, 0.0); + double *d_mu = dphi + 2 * n; + double *d_sig = dphi + 5 * n; + double *d_C = dphi + 8 * n; + double *d_r = dphi + 11 * n; + 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; + 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_r[i] = q * mu * dkb_dmu_b; + d_q[i] = kb; + } + } + static bool is_valid(const std::vector &par) { return par[3] > 0.0; } @@ -828,20 +1223,6 @@ struct GaussianChargeSharingKb { return {p0, p1, mu, sigma, N, C, kb_mean, kb_frac}; } - - static void compute_steps(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); - } }; // _____________________________________________________________________ @@ -908,12 +1289,12 @@ struct RisingScurve { const double dx = x - p2; const double z = dx * inv_sqrt2 / p3; - const double step = 0.5 * (1.0 + fast_erf(z)); + 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 e = std::exp(-z * z); const double dSdp2 = -inv_sqrt_2pi * e / p3; const double dSdp3 = -inv_sqrt_2pi * e * dx / (p3 * p3); @@ -925,6 +1306,79 @@ struct RisingScurve { g[5] = step * dx; // df/dp5 } + /** + * @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; + 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; + 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; @@ -995,18 +1449,6 @@ struct RisingScurve { return {p0, p1, p2, p3, p4, p5}; } - - static void compute_steps(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; - } }; // _____________________________________________________________________ @@ -1061,12 +1503,12 @@ struct FallingScurve { const double dx = x - p2; const double z = dx * inv_sqrt2 / p3; - const double step = 0.5 * (1.0 - fast_erf(z)); + 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 e = std::exp(-z * z); const double dSdp2 = +inv_sqrt_2pi * e / p3; // sign flipped vs rising const double dSdp3 = +inv_sqrt_2pi * e * dx / (p3 * p3); @@ -1078,6 +1520,80 @@ struct FallingScurve { 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; @@ -1148,14 +1664,6 @@ struct FallingScurve { return {p0, p1, p2, p3, p4, p5}; } - - static void compute_steps(const std::array &start, - double x_range, double y_range, - double slope_scale, - std::array &steps) { - RisingScurve::compute_steps(start, x_range, y_range, slope_scale, - steps); - } }; } // namespace aare::model diff --git a/include/aare/calibration.hpp b/include/aare/calibration.hpp index 28550b3a..aac09cba 100644 --- a/include/aare/calibration.hpp +++ b/include/aare/calibration.hpp @@ -11,20 +11,6 @@ namespace aare { -// Really try to convince the compile to inline this function -// TODO! Clang? -#if (defined(_MSC_VER) || defined(__INTEL_COMPILER)) -#define STRONG_INLINE __forceinline -#else -#define STRONG_INLINE inline -#endif - -#if defined(__GNUC__) -#define ALWAYS_INLINE __attribute__((always_inline)) inline -#else -#define ALWAYS_INLINE STRONG_INLINE -#endif - /** * @brief Get the gain from the raw ADC value. In Jungfrau the gain is * encoded in the left most 2 bits of the raw value. diff --git a/include/aare/defs.hpp b/include/aare/defs.hpp index b0701912..174bfb6b 100644 --- a/include/aare/defs.hpp +++ b/include/aare/defs.hpp @@ -32,6 +32,26 @@ } while (0) #endif +/** + * ALWAYS_INLINE forces a function to be inlined at every call site. gcc + * refuses to inline a plain `inline` function once its translation-unit + * growth budget (--param inline-unit-growth) is used up, which happens in + * the fitting code where every model is instantiated for every minimizer; + * a per-point helper that is then called out of line also stops the + * surrounding loop from vectorising. + */ +#if (defined(_MSC_VER) || defined(__INTEL_COMPILER)) +#define STRONG_INLINE __forceinline +#else +#define STRONG_INLINE inline +#endif + +#if defined(__GNUC__) +#define ALWAYS_INLINE __attribute__((always_inline)) inline +#else +#define ALWAYS_INLINE STRONG_INLINE +#endif + namespace aare { inline constexpr size_t bits_per_byte = 8; diff --git a/python/aare/__init__.py b/python/aare/__init__.py index b4be0755..570176f9 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -40,6 +40,7 @@ from .Cluster import Cluster from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb from ._aare import fit +from ._aare import Minimizer from ._aare import Interpolator from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2 from ._aare import reduce_to_2x2, reduce_to_3x3 diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 1e8d9c39..f26cc20c 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -22,9 +22,40 @@ fit_dispatch(const aare::FitModel &model, template void bind_fit_model(py::module &m, const char *name) { using FM = aare::FitModel; py::class_(m, name) - .def(py::init(), + .def(py::init(), py::arg("strategy") = 0, py::arg("max_calls") = 100, - py::arg("tolerance") = 0.5, py::arg("compute_errors") = false) + py::arg("tolerance") = 0.5, py::arg("compute_errors") = false, + py::arg("minimizer") = aare::Minimizer::Migrad, + R"doc( + Fit configuration for this model. + + Parameters + ---------- + strategy : int + Minuit2 strategy (0 = fast, 1 = Minuit2's default). Ignored + by ``Minimizer.LevenbergMarquardt`` and + ``Minimizer.VarPro``. + max_calls : int + Work budget per pixel: Minuit2 function calls, or model + evaluations for ``Minimizer.LevenbergMarquardt`` and + ``Minimizer.VarPro`` (one per iteration, two + when a step crosses a limit). 0 selects Minuit's default + budget. + tolerance : float + EDM tolerance in Minuit's convention. Migrad, + LevenbergMarquardt and VarPro stop below + ``0.002 * tolerance``, Fumili below ``1e-4 * tolerance``. + compute_errors : bool + Also return parameter errors: from Hesse with Migrad, from + the linearised covariance with Fumili, LevenbergMarquardt + and VarPro. Fixed parameters and parameters + ending on a limit report 0. + minimizer : Minimizer + ``Minimizer.Migrad`` (default), ``Minimizer.Fumili``, + ``Minimizer.LevenbergMarquardt`` or + ``Minimizer.VarPro``. + )doc") .def("SetParLimits", py::overload_cast(&FM::SetParLimits), py::arg("idx"), py::arg("lo"), py::arg("hi")) @@ -59,6 +90,7 @@ template void bind_fit_model(py::module &m, const char *name) { .def_property("tolerance", &FM::tolerance, &FM::SetTolerance) .def_property("compute_errors", &FM::compute_errors, &FM::SetComputeErrors) + .def_property("minimizer", &FM::minimizer, &FM::SetMinimizer) .def( "__call__", // conversion ok, we want to be able to call with any // dtype @@ -89,6 +121,11 @@ template void bind_fit_model(py::module &m, const char *name) { R"doc( Fit this model to 1D or 3D data using Minuit2. + The minimizer is selected by the ``minimizer`` property + (``Minimizer.Migrad`` by default, ``Minimizer.Fumili``, + ``Minimizer.LevenbergMarquardt`` or + ``Minimizer.VarPro``). + Parameters ---------- x : array_like, shape (n_scan,) @@ -182,8 +219,23 @@ fit_dispatch(const aare::FitModel &model, "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, @@ -217,7 +269,82 @@ fit_dispatch(const aare::FitModel &model, } } +// 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; + +#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); \ + } + + 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) + +#undef AARE_DISPATCH_MODEL + + throw std::runtime_error( + "Unknown model type. Expected Pol1, Pol2, Gaussian, " + "GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb, " + "RisingScurve or FallingScurve."); +} + void define_fit_bindings(py::module &m) { + // The enum must exist before it is used as a constructor default below. + py::enum_(m, "Minimizer", R"doc( + Minimizer used by the fit models. + + ``Migrad`` (default) is Minuit2's variable-metric minimizer. With + ``compute_errors`` its parameter errors come from Hesse. + + ``Fumili`` is Minuit2's Gauss-Newton minimizer with a trust region. + It takes the gradient and a linearised Hessian from the analytic + derivatives of the model and usually needs far fewer function + evaluations than Migrad. Its parameter errors come from the + covariance of the linearised Hessian, without an extra Hesse step. + Minuit2's Fumili cannot converge once a two-sided limit becomes + active; a pixel without a valid Fumili minimum is refitted with + Migrad. + + ``LevenbergMarquardt`` is the built-in damped Gauss-Newton solver. It + uses the analytic Jacobian, reflects steps at parameter limits and + reuses its buffers between pixels, so data cubes are fitted without + per-pixel allocations. Its parameter errors are + ``sqrt(diag((J^T J)^-1))``; ``max_calls`` counts model evaluations + and ``strategy`` is ignored. + + ``VarPro`` is the built-in variable projection solver. + Every bundled model is linear in some of its parameters (the + amplitude of a Gaussian, the baseline and amplitude of an S-curve, + every parameter of a polynomial); each trial point solves those + exactly and the Levenberg-Marquardt iteration runs over the + remaining nonlinear parameters only. It needs fewer evaluations than + ``LevenbergMarquardt`` and cannot be misled by poor start values of + the linear parameters, which it ignores. A pixel whose linear + solution violates a limit on a linear parameter, or that does not + converge, is refitted with ``LevenbergMarquardt``. Its parameter + errors are the same Gauss-Newton estimates; ``max_calls`` counts + evaluations and ``strategy`` is ignored. + + With every minimizer, fixed parameters and parameters that end on a + limit report an error of 0 and a failed fit returns zeros. + )doc") + .value("Migrad", aare::Minimizer::Migrad) + .value("Fumili", aare::Minimizer::Fumili) + .value("LevenbergMarquardt", aare::Minimizer::LevenbergMarquardt) + .value("VarPro", aare::Minimizer::VarPro); + // ── Bind model classes ────────────────────────────────────────── bind_fit_model(m, "Gaussian"); bind_fit_model(m, "GaussianErfcPlateau"); @@ -237,86 +364,19 @@ void define_fit_bindings(py::module &m) { py::array_t x, py::array_t y, py::object y_err_obj, int n_threads) -> py::object { - using namespace aare::model; - - // ── Polynomial of degree 1 ─────── - if (py::isinstance>(model_obj)) { - const auto &mdl = - model_obj.cast &>(); - return fit_dispatch(mdl, x, y, y_err_obj, n_threads); - } - - // ── Polynomial of degree 2 ─────── - if (py::isinstance>(model_obj)) { - const auto &mdl = - model_obj.cast &>(); - return fit_dispatch(mdl, x, y, y_err_obj, n_threads); - } - - // ── Gaussian ─────── - if (py::isinstance>(model_obj)) { - const auto &mdl = - model_obj.cast &>(); - return fit_dispatch(mdl, x, y, y_err_obj, n_threads); - } - - // ── GaussianErfcPlateau ─────── - if (py::isinstance>( - model_obj)) { - const auto &mdl = - model_obj - .cast &>(); - return fit_dispatch(mdl, x, y, y_err_obj, - n_threads); - } - - // ── GaussianChargeSharing ─────── - if (py::isinstance>( - model_obj)) { - const auto &mdl = - model_obj - .cast &>(); - return fit_dispatch(mdl, x, y, y_err_obj, - n_threads); - } - - // ── GaussianChargeSharingKb ─────── - if (py::isinstance>( - model_obj)) { - const auto &mdl = model_obj.cast< - const aare::FitModel &>(); - return fit_dispatch( - mdl, x, y, y_err_obj, n_threads); - } - - // ── Rising Scurve ─────── - if (py::isinstance>(model_obj)) { - const auto &mdl = - model_obj.cast &>(); - return fit_dispatch(mdl, x, y, y_err_obj, - n_threads); - } - - // ── Falling Scurve ─────── - if (py::isinstance>(model_obj)) { - const auto &mdl = - model_obj.cast &>(); - return fit_dispatch(mdl, x, y, y_err_obj, - n_threads); - } - - throw std::runtime_error( - "Unknown model type. Expected Pol1, Pol2, Gaussian, " - "RisingScurve or FallingScurve."); + return dispatch_any_model(model_obj, x, y, y_err_obj, n_threads); }, R"( - Fit a model to 1D or 3D data using Minuit2. + Fit a model to 1D or 3D data with the minimizer selected by the + model's ``minimizer`` property. Parameters ---------- - model : Pol1, Pol2, Gaussian, RisingScurve, or FallingScurve - Configured model object. User-set limits, fixed parameters, - and start values take precedence over automatic estimates. + model : object + Configured model object: Pol1, Pol2, Gaussian, GaussianErfcPlateau, + GaussianChargeSharing, GaussianChargeSharingKb, RisingScurve or + FallingScurve. User-set limits, fixed parameters, and start values + take precedence over automatic estimates. x : array_like, shape (n_scan,) Scan points (e.g. energy or threshold values). y : array_like, shape (n_scan,) or (rows, cols, n_scan) @@ -328,15 +388,12 @@ void define_fit_bindings(py::module &m) { Returns ------- - For 1D input: - numpy array of shape (2*npar+1,) if compute_errors else (npar+1,). - Layout: [params..., (errors...,) chi2]. - - For 3D input: - dict with keys: - "par" : (rows, cols, npar) fitted parameters. - "par_err" : (rows, cols, npar) parameter errors (if compute_errors). - "chi2" : (rows, cols) chi-squared per pixel. + dict + ``par`` holds the fitted parameters, shape (npar,) for 1D input or + (rows, cols, npar) for 3D input. ``chi2`` holds the chi-squared, + shape (1,) or (rows, cols). ``par_err`` holds the parameter + 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); diff --git a/python/tests/test_Fit.py b/python/tests/test_Fit.py index 3a477b27..2e876ed9 100644 --- a/python/tests/test_Fit.py +++ b/python/tests/test_Fit.py @@ -1,15 +1,401 @@ # SPDX-License-Identifier: MPL-2.0 import numpy as np +import pytest + import aare +MINIMIZERS = [ + aare.Minimizer.Migrad, + aare.Minimizer.Fumili, + aare.Minimizer.LevenbergMarquardt, + aare.Minimizer.VarPro, +] -def test_gaussian_model_evaluates_and_fits_data(): +def gaussian(x, A, mu, sigma): + return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2) + + +def uses_minuit2(minimizer): + return minimizer not in ( + aare.Minimizer.LevenbergMarquardt, + aare.Minimizer.VarPro, + ) + + +@pytest.fixture(params=MINIMIZERS, ids=lambda m: m.name) +def minimizer(request): + return request.param + + +@pytest.fixture +def gaussian_data(): + x = np.linspace(-6.0, 6.0, 61) + y = gaussian(x, 120.0, 0.8, 1.3) + y_err = np.ones_like(y) + return x, y, y_err + + +def test_gaussian_model_evaluates_and_fits_data(minimizer): x = np.linspace(-5.0, 5.0, 51) expected = np.array([20.0, 0.5, 1.2]) - model = aare.Gaussian() + model = aare.Gaussian(minimizer=minimizer) y = model(x, expected) result = model.fit(x, y) np.testing.assert_allclose(result["par"], expected, atol=2e-3) + + +def test_minimizer_can_be_selected(): + assert aare.Gaussian().minimizer == aare.Minimizer.Migrad + assert set(aare.Minimizer.__members__) == { + "Migrad", + "Fumili", + "LevenbergMarquardt", + "VarPro", + } + + for minimizer in MINIMIZERS: + assert aare.Gaussian(minimizer=minimizer).minimizer == minimizer + + model = aare.Gaussian() + model.minimizer = aare.Minimizer.LevenbergMarquardt + assert model.minimizer == aare.Minimizer.LevenbergMarquardt + + +def test_fit_returns_errors_and_chi2(gaussian_data, minimizer): + x, y, y_err = gaussian_data + model = aare.Gaussian(compute_errors=True, minimizer=minimizer) + result = model.fit(x, y, y_err) + + assert set(result) == {"par", "par_err", "chi2"} + assert result["par"].shape == (3,) + assert result["par_err"].shape == (3,) + assert result["chi2"].shape == (1,) + np.testing.assert_allclose(result["par"], [120.0, 0.8, 1.3], rtol=1e-4) + np.testing.assert_allclose( + result["par_err"], [0.361, 0.00451, 0.00451], rtol=2e-2 + ) + assert result["chi2"][0] == pytest.approx(0.0, abs=1e-3) + + without_errors = aare.fit(aare.Gaussian(minimizer=minimizer), x, y) + assert set(without_errors) == {"par", "chi2"} + + +def test_pol1_weighted_errors_match_analytic_formula(minimizer): + x = np.arange(20, dtype=float) + y = 2.0 + 0.5 * x + y_err = 1.0 + 0.1 * x + w = 1.0 / y_err**2 + S, Sx, Sxx = w.sum(), (w * x).sum(), (w * x * x).sum() + delta = S * Sxx - Sx * Sx + + result = aare.Pol1(compute_errors=True, minimizer=minimizer).fit(x, y, y_err) + + np.testing.assert_allclose(result["par"], [2.0, 0.5], rtol=1e-6) + np.testing.assert_allclose( + result["par_err"], + [np.sqrt(Sxx / delta), np.sqrt(S / delta)], + rtol=1e-6, + ) + + +def test_fixed_parameter_reports_zero_error(gaussian_data, minimizer): + x, y, y_err = gaussian_data + model = aare.Gaussian(compute_errors=True, minimizer=minimizer) + model.FixParameter("sigma", 1.5) + result = model.fit(x, y, y_err) + + assert result["par"][2] == 1.5 + assert result["par_err"][2] == 0.0 + assert result["par"][0] == pytest.approx(111.15, rel=1e-3) + assert result["chi2"][0] == pytest.approx(1684.25, rel=1e-3) + + model.ReleaseParameter("sigma") + result = model.fit(x, y, y_err) + assert result["par"][2] == pytest.approx(1.3, rel=1e-4) + + +def test_active_limit_ends_on_the_limit(gaussian_data, minimizer): + x, y, y_err = gaussian_data + model = aare.Gaussian(compute_errors=True, minimizer=minimizer) + model.SetParLimits("mu", -5.0, 0.5) + result = model.fit(x, y, y_err) + + margin = 1e-6 if uses_minuit2(minimizer) else 1e-9 + assert result["par"][1] == pytest.approx(0.5, abs=margin) + assert result["par_err"][1] == 0.0 + assert result["par"][0] == pytest.approx(116.88, rel=1e-3) + assert result["chi2"][0] == pytest.approx(4301.95, rel=1e-3) + + inactive = aare.Gaussian(minimizer=minimizer) + inactive.SetParLimits(1, -5.0, 5.0) + np.testing.assert_allclose( + inactive.fit(x, y)["par"], [120.0, 0.8, 1.3], rtol=1e-4 + ) + + +def test_user_start_value(gaussian_data, minimizer): + x, y, _ = gaussian_data + model = aare.Gaussian(minimizer=minimizer) + model.SetParameter("mu", 0.0) + np.testing.assert_allclose(model.fit(x, y)["par"], [120.0, 0.8, 1.3], rtol=1e-4) + + +def test_max_calls_failure_returns_zeros(gaussian_data, minimizer): + x, y, y_err = gaussian_data + model = aare.Gaussian(max_calls=3, compute_errors=True, minimizer=minimizer) + # start far from the minimum so that no minimizer converges in 3 calls + model.SetParameter("mu", -3.0) + model.SetParameter("sigma", 3.0) + result = model.fit(x, y, y_err) + assert np.all(result["par"] == 0.0) + assert np.all(result["par_err"] == 0.0) + assert np.all(result["chi2"] == 0.0) + + +def test_strategy_is_ignored_by_levenberg_marquardt(gaussian_data): + x, y, _ = gaussian_data + lm = aare.Minimizer.LevenbergMarquardt + a = aare.Gaussian(strategy=0, minimizer=lm).fit(x, y) + b = aare.Gaussian(strategy=1, minimizer=lm).fit(x, y) + np.testing.assert_array_equal(a["par"], b["par"]) + + +def test_invalid_parameter_access_raises(): + model = aare.Gaussian() + assert model.par_names == ["A", "mu", "sigma"] + assert model.n_par == 3 + with pytest.raises(RuntimeError): + model.SetParameter("nope", 1.0) + with pytest.raises(IndexError): + model.SetParameter(3, 1.0) + with pytest.raises(RuntimeError): + model.SetParLimits("mu", 1.0, 1.0) + + +def test_mismatched_array_sizes_raise(gaussian_data, minimizer): + x, y, y_err = gaussian_data + model = aare.Gaussian(minimizer=minimizer) + with pytest.raises(RuntimeError): + model.fit(x, y[:-1]) + with pytest.raises(RuntimeError): + model.fit(x, y, y_err[:-1]) + + +@pytest.mark.parametrize("cls", [aare.RisingScurve, aare.FallingScurve]) +def test_scurve_fit(cls, minimizer): + truth = np.array([10.0, 0.1, 50.0, 3.0, 1000.0, 2.0]) + x = np.linspace(0.0, 100.0, 100) + model = cls(max_calls=500, compute_errors=True, minimizer=minimizer) + y = model(x, truth) + result = model.fit(x, y, np.ones_like(y)) + + np.testing.assert_allclose(result["par"], truth, rtol=1e-3) + assert np.all(result["par_err"] > 0.0) + + +@pytest.mark.parametrize("with_errors", [False, True]) +@pytest.mark.parametrize("weighted", [False, True]) +def test_fit_3d_matches_fit_1d(with_errors, weighted, minimizer): + x = np.linspace(-6.0, 6.0, 61) + rows, cols = 4, 3 + y = np.empty((rows, cols, x.size)) + for r in range(rows): + for c in range(cols): + y[r, c] = gaussian(x, 80.0 + 10.0 * r, -0.6 + 0.3 * c, 0.9 + 0.1 * r) + y_err = np.ones_like(y) if weighted else None + + model = aare.Gaussian(compute_errors=with_errors, minimizer=minimizer) + cube = model.fit(x, y, y_err, n_threads=3) + + expected_keys = {"par", "par_err", "chi2"} if with_errors else {"par", "chi2"} + assert set(cube) == expected_keys + assert cube["par"].shape == (rows, cols, 3) + assert cube["chi2"].shape == (rows, cols) + if with_errors: + assert cube["par_err"].shape == (rows, cols, 3) + for r in range(rows): + for c in range(cols): + single = model.fit(x, y[r, c], None if y_err is None else y_err[r, c]) + np.testing.assert_allclose(cube["par"][r, c], single["par"], rtol=1e-10) + assert cube["chi2"][r, c] == pytest.approx(single["chi2"][0], abs=1e-10) + if with_errors: + np.testing.assert_allclose( + cube["par_err"][r, c], single["par_err"], rtol=1e-10 + ) + np.testing.assert_allclose(cube["par"][3, 2], [110.0, 0.0, 1.2], rtol=1e-4, atol=1e-6) + + +def test_all_minimizers_agree_on_3d_data_with_errors(): + x = np.linspace(-5.0, 5.0, 51) + expected = np.array([20.0, 0.5, 1.2]) + single = aare.Gaussian()(x, expected) + y = np.tile(single, (2, 3, 1)) + y_err = np.ones_like(y) + + results = { + minimizer: aare.Gaussian(compute_errors=True, minimizer=minimizer).fit( + x, y, y_err, n_threads=2 + ) + for minimizer in MINIMIZERS + } + reference = results[aare.Minimizer.Migrad] + for minimizer, result in results.items(): + assert result["par"].shape == (2, 3, 3) + assert result["par_err"].shape == (2, 3, 3) + np.testing.assert_allclose(result["par"], np.broadcast_to(expected, (2, 3, 3)), atol=2e-3) + np.testing.assert_allclose(result["chi2"], 0.0, atol=1e-3) + np.testing.assert_allclose(result["par_err"], reference["par_err"], rtol=2e-2) + + +# Regressions of the built-in solvers: fits that used to be reported as valid +# although they had stopped away from the minimum or returned garbage. + + +def test_gaussian_fit_to_noise_fails_or_stays_bounded(minimizer): + # A width collapsing between two scan points let the amplitude diverge + # (up to 1e159) in a fit that was still reported as valid. + x = np.linspace(-6.0, 6.0, 61) + y = np.random.default_rng(0).normal(0.0, 1.0, (20, 20, x.size)) + result = aare.Gaussian(minimizer=minimizer).fit(x, y, n_threads=2) + + par = result["par"] + failed = np.all(par == 0.0, axis=-1) + bounded = np.all(np.isfinite(par) & (np.abs(par) < 1e6), axis=-1) + assert np.all(failed | bounded) + assert np.all(np.isfinite(result["chi2"])) + + +def test_step_reflected_onto_start_is_not_convergence(minimizer): + # The step of +2 reflects at the limit 1 back onto the start 0; the + # truncated step must still be tried. + model = aare.Pol1(minimizer=minimizer) + model.SetParameter("p0", 0.0) + model.FixParameter("p1", 0.0) + model.SetParLimits("p0", -1.0, 1.0) + result = model.fit(np.arange(4.0), np.full(4, 2.2)) + + assert result["par"][0] == pytest.approx(1.0, abs=1e-6) + assert result["chi2"][0] == pytest.approx(5.76, rel=1e-6) + + +@pytest.mark.parametrize( + "amplitude, offset", [(1e8, 0.0), (1e10, 0.0), (120.0, 1e9)] +) +def test_large_parameter_does_not_stop_the_fit(minimizer, amplitude, offset): + # One relative step norm over all parameters let a large amplitude or + # position stop the fit at the start estimate of the width. + x = np.linspace(-6.0, 6.0, 61) + y = gaussian(x, amplitude, 0.8, 1.3) + y_err = np.sqrt(y) + 1.0 + result = aare.Gaussian(minimizer=minimizer, max_calls=500).fit( + x + offset, y, y_err + ) + + par = result["par"] + assert par[0] == pytest.approx(amplitude, rel=1e-3) + assert par[1] - offset == pytest.approx(0.8, abs=1e-3) + assert par[2] == pytest.approx(1.3, rel=1e-3) + + +def test_rejected_trial_on_limit_does_not_freeze_interior_parameter(minimizer): + # A rejected trial truncated onto the lower limit froze mu at its start. + x = np.linspace(-0.1, 0.1, 11) + mu = np.sqrt(2.0 * np.log(2.0)) + y = gaussian(x, 10.0, mu, 1.0) + model = aare.Gaussian(minimizer=minimizer) + model.FixParameter("A", 10.0) + model.FixParameter("sigma", 1.0) + model.SetParameter("mu", 3.0) + model.SetParLimits("mu", 0.0, 5.0) + result = model.fit(x, y) + + assert result["par"][1] == pytest.approx(mu, abs=1e-3) + assert result["chi2"][0] < 1e-4 + + +# Samples with a zero error are ignored, also when they hold NaN (a masked +# dead frame): a zero weight does not remove a NaN from a sum. + + +@pytest.mark.parametrize("masked", [4, 0, 9]) +def test_masked_nan_sample_is_ignored(minimizer, masked): + # The start estimate of Pol1 uses the first and the last sample. + x = np.arange(10.0) + y = 2.0 + 0.5 * x + y_err = np.ones_like(y) + y[masked] = np.nan + y_err[masked] = 0.0 + result = aare.Pol1(minimizer=minimizer, compute_errors=True).fit(x, y, y_err) + + np.testing.assert_allclose(result["par"], [2.0, 0.5], atol=1e-5) + assert np.all(result["par_err"] > 0.0) + assert result["chi2"][0] == pytest.approx(0.0, abs=1e-6) + + +def test_too_few_unmasked_samples_fail(minimizer): + x = np.arange(10.0) + y = np.full_like(x, np.nan) + y_err = np.zeros_like(x) + y[0], y_err[0] = 2.0, 1.0 # one sample left for two free parameters + result = aare.Pol1(minimizer=minimizer).fit(x, y, y_err) + assert np.all(result["par"] == 0.0) + assert np.all(result["chi2"] == 0.0) + + +def test_masked_samples_in_cube_match_1d_fits(minimizer): + x = np.linspace(-6.0, 6.0, 61) + y = np.stack([gaussian(x, 80.0 + 10.0 * c, 0.8, 1.2) for c in range(4)]) + y = y.reshape(1, 4, x.size) + y_err = np.ones_like(y) + for col, masked in [(1, 10), (2, 0), (3, 34)]: + y[0, col, masked] = np.nan + y_err[0, col, masked] = 0.0 + model = aare.Gaussian(minimizer=minimizer, compute_errors=True) + cube = model.fit(x, y, y_err, n_threads=1) + + for col in range(4): + single = model.fit(x, y[0, col], y_err[0, col]) + np.testing.assert_allclose(cube["par"][0, col], single["par"], atol=1e-12) + np.testing.assert_allclose(cube["par_err"][0, col], single["par_err"], atol=1e-12) + np.testing.assert_allclose( + cube["par"][0, col], [80.0 + 10.0 * col, 0.8, 1.2], rtol=1e-3 + ) + + +def test_parameter_whose_minimum_lies_on_its_limit_reports_no_error(minimizer): + # At a minimum that lies exactly on the limit the gradient is zero. + x = np.array([-1.0, -1.0, 1.0, 1.0]) + model = aare.Pol1(minimizer=minimizer, compute_errors=True) + model.SetParLimits("p0", 0.0, 1.0) + result = model.fit(x, 1.0 + 2.0 * x) + + assert result["par"][0] == pytest.approx(1.0, abs=1e-6) + assert result["par_err"][0] == 0.0 + assert result["par_err"][1] == pytest.approx(0.5, rel=1e-6) + + +@pytest.mark.parametrize("amplitude", [0.1, 1e-4]) +def test_fit_does_not_stop_at_start_without_a_step(minimizer, amplitude): + # The EDM tolerance is absolute in chi2 units: an unweighted fit of small + # data passes it at the start estimate (sigma 1.3617) already. + x = np.linspace(-6.0, 6.0, 61) + y = gaussian(x, amplitude, 0.8, 1.3) + result = aare.Gaussian(minimizer=minimizer).fit(x, y) + assert abs(result["par"][2] - 1.3) < 0.03 # Migrad ends at 1.3205 + + one_call = aare.Gaussian(minimizer=minimizer, max_calls=1).fit(x, y) + assert np.all(one_call["par"] == 0.0) + + +@pytest.mark.parametrize("lo, span", [(8000, 1000), (20000, 500), (60000, 200)]) +def test_polynomial_far_from_zero_keeps_its_precision(minimizer, lo, span): + # Nearly collinear columns 1, x, x^2: VarPro's normal equations lost six + # digits at x = 20000. + x = np.linspace(lo, lo + span, 40) + c = np.array([1.5 - 2e-3 * lo + 3e-7 * lo**2, 2e-3 - 6e-7 * lo, 3e-7]) + y = c[0] + c[1] * x + c[2] * x**2 + result = aare.Pol2(minimizer=minimizer).fit(x, y) + np.testing.assert_allclose(result["par"], c, rtol=1e-9) diff --git a/src/Chi2.hpp b/src/Chi2.hpp deleted file mode 100644 index 2aeb0b93..00000000 --- a/src/Chi2.hpp +++ /dev/null @@ -1,140 +0,0 @@ -// SPDX-License-Identifier: MPL-2.0 -#pragma once -#include "aare/Models.hpp" -#include "aare/NDView.hpp" -#include -#include -#include -#include -#include - -namespace aare { - -namespace func { - -/** - * @brief Generic chi-squared FCN with analytic gradient for a 1D model. - * - * @tparam Model A model struct that satisfies: - * - static constexpr std::size_t npar; - * - static double eval(double x, const std::vector& par); - * - static void eval_and_grad(double x, const std::vector& par, - * double& f, std::array& g); - * - static bool is_valid(const std::vector& par); - * - * Gradient: - * d(chi2)/dp_k = -2 * sum_i w_i * (y_i - f_i) * df_i/dp_k - * - * where w_i = 1/sigma_i^2 (weighted) or 1 (unweighted). - * - * By providing analytic gradients we avoid 2*npar extra function evaluations - * per Minuit step that would otherwise be spent on finite differences. - * - * @throws std::invalid_argument if par.size() != Model::npar. - * - * Invalid model parameters do not throw; they return a large penalty - * (and a zero gradient fallback) so the minimizer can remain in control. - */ -template -class Chi2Model1DGrad : public ROOT::Minuit2::FCNGradientBase { - public: - Chi2Model1DGrad(NDView x, NDView y) - : x_(x), y_(y), s_(), weighted_(false) {} - - Chi2Model1DGrad(NDView x, NDView y, - NDView y_err) - : x_(x), y_(y), s_(y_err), weighted_(true) {} - - ~Chi2Model1DGrad() override = default; - - double operator()(const std::vector &par) const override { - if (par.size() != Model::npar) { - throw std::invalid_argument( - "Chi2Model1DGrad: wrong parameter vector size."); - } - - if (!Model::is_valid(par)) - return 1e20; - - double chi2 = 0.0; - - if (weighted_) { - for (ssize_t i = 0; i < x_.size(); ++i) { - const double si = s_[i]; - if (si == 0.0) - continue; - - const double f_i = Model::eval(x_[i], par); - const double r_i = y_[i] - f_i; - chi2 += (r_i * r_i) / (si * si); - } - } else { - for (ssize_t i = 0; i < x_.size(); ++i) { - const double f_i = Model::eval(x_[i], par); - const double r_i = y_[i] - f_i; - chi2 += r_i * r_i; - } - } - - return chi2; - } - - std::vector - Gradient(const std::vector &par) const override { - if (par.size() != Model::npar) { - throw std::invalid_argument( - "Chi2Model1DGrad: wrong parameter vector size."); - } - - std::vector grad(Model::npar, 0.0); - if (!Model::is_valid(par)) - return grad; - - std::array df{}; - double f_i = 0.0; - - if (weighted_) { - for (ssize_t i = 0; i < x_.size(); ++i) { - const double si = s_[i]; - if (si == 0.0) - continue; - - Model::eval_and_grad(x_[i], par, f_i, df); - - const double r_i = y_[i] - f_i; - const double c = -2.0 * r_i / (si * si); - - for (std::size_t k = 0; k < Model::npar; ++k) { - grad[k] += c * df[k]; - } - } - } else { - for (ssize_t i = 0; i < x_.size(); ++i) { - Model::eval_and_grad(x_[i], par, f_i, df); - - const double r_i = y_[i] - f_i; - const double c = -2.0 * r_i; - - for (std::size_t k = 0; k < Model::npar; ++k) { - grad[k] += c * df[k]; - } - } - } - - return grad; - } - - /** @brief Error definition: 1.0 for chi-squared (delta_chi2 = 1 -> - * 1-sigma). */ - double Up() const override { return 1.0; } - - private: - NDView x_; - NDView y_; - NDView s_; - bool weighted_; -}; - -} // namespace func - -} // namespace aare \ No newline at end of file diff --git a/src/Fit.cpp b/src/Fit.cpp deleted file mode 100644 index f11aa611..00000000 --- a/src/Fit.cpp +++ /dev/null @@ -1,308 +0,0 @@ -// SPDX-License-Identifier: MPL-2.0 -#include "aare/Fit.hpp" -#include "Chi2.hpp" -#include "Minuit2/FunctionMinimum.h" -#include "Minuit2/MnHesse.h" -#include "Minuit2/MnMigrad.h" -#include "Minuit2/MnStrategy.h" -#include "Minuit2/MnUserParameters.h" -#include "aare/Models.hpp" -#include "aare/utils/par.hpp" -#include "aare/utils/task.hpp" -#include -#include -#include -#include -#include - -namespace aare { -// ============================================================================ -// FitModel — method definitions -// (constructor, destructor, copy, and all methods that touch Minuit2 state) -// ============================================================================ - -template struct FitModel::FitModelImpl { - ROOT::Minuit2::MnUserParameters upar; - ROOT::Minuit2::MnStrategy strategy; - - explicit FitModelImpl(unsigned int strategy_level) - : strategy(strategy_level) {} - - FitModelImpl(const FitModelImpl &) = default; - FitModelImpl &operator=(const FitModelImpl &) = default; -}; - -template -FitModel::FitModel(unsigned int strategy, unsigned int max_calls, - double tolerance, bool compute_errors) - : impl_(std::make_unique(strategy)), max_calls_(max_calls), - tolerance_(tolerance), compute_errors_(compute_errors) { - for (std::size_t i = 0; i < npar; ++i) { - const auto pi = Model::param_info[i]; - const bool has_lo = std::isfinite(pi.default_lo); - const bool has_hi = std::isfinite(pi.default_hi); - if (has_lo && has_hi) { - impl_->upar.Add(pi.name, 0.0, 1.0, pi.default_lo, pi.default_hi); - } else if (has_lo) { - impl_->upar.Add(pi.name, 0.0, 1.0, pi.default_lo, 1e6); - } else { - impl_->upar.Add(pi.name, 0.0, 1.0); - } - } -} - -template FitModel::~FitModel() = default; - -template -FitModel::FitModel(const FitModel &other) - : impl_(std::make_unique(*other.impl_)), - max_calls_(other.max_calls_), tolerance_(other.tolerance_), - compute_errors_(other.compute_errors_), user_fixed_(other.user_fixed_), - user_start_(other.user_start_) {} - -template -FitModel &FitModel::operator=(const FitModel &other) { - if (this != &other) { - impl_ = std::make_unique(*other.impl_); - max_calls_ = other.max_calls_; - tolerance_ = other.tolerance_; - compute_errors_ = other.compute_errors_; - user_fixed_ = other.user_fixed_; - user_start_ = other.user_start_; - } - return *this; -} - -template -unsigned int FitModel::checked_index(const std::string &name) const { - for (std::size_t i = 0; i < npar; ++i) { - if (impl_->upar.Name(i) == name) - return static_cast(i); - } - throw std::runtime_error("FitModel: unknown parameter name '" + name + "'"); -} - -template -void FitModel::SetParLimits(unsigned int idx, double lo, double hi) { - impl_->upar.SetLimits(idx, lo, hi); -} - -template -void FitModel::FixParameter(unsigned int idx, double val) { - SetParameter(idx, val); - impl_->upar.Fix(idx); - user_fixed_[idx] = true; -} - -template -void FitModel::ReleaseParameter(unsigned int idx) { - impl_->upar.Release(idx); - user_fixed_[idx] = false; -} - -template -void FitModel::ReleaseParameter(const std::string &name) { - ReleaseParameter(checked_index(name)); -} - -template -void FitModel::SetParameter(unsigned int idx, double val) { - impl_->upar.SetValue(idx, val); - user_start_[idx] = true; -} - -template -void FitModel::SetParameter(const std::string &name, double val) { - SetParameter(checked_index(name), val); -} - -template -void FitModel::FixParameter(const std::string &name, double val) { - FixParameter(checked_index(name), val); -} - -template -void FitModel::SetParLimits(const std::string &name, double lo, - double hi) { - SetParLimits(checked_index(name), lo, hi); -} - -template -std::string FitModel::GetParName(unsigned int idx) const { - return impl_->upar.GetName(idx); -} - -template -std::vector FitModel::GetParNames() const { - std::vector names; - for (std::size_t i = 0; i < npar; ++i) - names.push_back(GetParName(i)); - return names; -} - -// ============================================================================ -// fit_pixel / fit_3d — Minuit2 template implementations -// ============================================================================ - -template -NDArray fit_pixel(const FitModel &model, NDView x, - NDView y, NDView y_err) { - using FCN = func::Chi2Model1DGrad; - constexpr std::size_t npar = Model::npar; - const bool want_errors = model.compute_errors(); - const ssize_t result_size = want_errors ? (2 * npar + 1) : (npar + 1); - - auto start = Model::estimate_par(x, y); - - if (!Model::is_valid(std::vector(start.begin(), start.end()))) { - return NDArray({result_size}, 0.0); - } - - double x_range, y_range, slope_scale; - model::compute_ranges(x, y, x_range, y_range, slope_scale); - - std::array steps{}; - Model::compute_steps(start, x_range, y_range, slope_scale, steps); - - // thread-local copy of starting parameters - auto upar_local = model.impl()->upar; - - for (std::size_t i = 0; i < npar; ++i) { - if (model.is_user_fixed(i)) - continue; - if (!model.is_user_start(i)) - upar_local.SetValue(i, start[i]); - upar_local.SetError(i, steps[i]); - } - - auto chi2 = (y_err.size() > 0) ? FCN(x, y, y_err) : FCN(x, y); - - ROOT::Minuit2::MnMigrad migrad(chi2, upar_local, model.impl()->strategy); - ROOT::Minuit2::FunctionMinimum min = - migrad(model.max_calls(), model.tolerance()); - - if (!min.IsValid()) - return NDArray({result_size}, 0.0); - - if (want_errors) { - ROOT::Minuit2::MnHesse hesse; - hesse(chi2, min); - - const auto &values = min.UserState().Params(); - const auto &errors = min.UserState().Errors(); - - NDArray result({result_size}); - for (std::size_t k = 0; k < npar; ++k) { - result[k] = values[k]; - result[npar + k] = errors[k]; - } - result[2 * npar] = min.Fval(); - return result; - } - - const auto &values = min.UserState().Params(); - NDArray result({result_size}); - for (std::size_t k = 0; k < npar; ++k) - result[k] = values[k]; - result[npar] = min.Fval(); - return result; -} - -template -NDArray fit_pixel(const FitModel &model, NDView x, - NDView y) { - return fit_pixel(model, x, y, NDView{}); -} - -template -void fit_3d(const FitModel &model, NDView x, - NDView y, NDView y_err, - NDView par_out, NDView err_out, - NDView chi2_out, int n_threads) { - const std::size_t npar = Model::npar; - - if (x.size() != y.shape(2)) - throw std::runtime_error("fit_3d: x.size() must match y.shape(2)."); - - if (par_out.shape(0) != y.shape(0) || par_out.shape(1) != y.shape(1) || - par_out.shape(2) != npar) - throw std::runtime_error("par_out must have shape [rows, cols, npar]."); - - 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 has_errors = (y_err.size() > 0); - const bool want_par_errors = (err_out.size() > 0) && model.compute_errors(); - - if (has_errors) { - 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.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]."); - } - - auto process = [&](ssize_t first_row, ssize_t last_row) { - for (ssize_t row = first_row; row < last_row; row++) { - for (ssize_t col = 0; col < y.shape(1); col++) { - - NDView values(&y(row, col, 0), {y.shape(2)}); - NDView errors = - has_errors ? NDView(&y_err(row, col, 0), - {y_err.shape(2)}) - : NDView{}; - - auto res = fit_pixel(model, x, values, errors); - - for (std::size_t k = 0; k < npar; ++k) - par_out(row, col, k) = res(k); - - if (want_par_errors) { - for (std::size_t k = 0; k < npar; ++k) - err_out(row, col, k) = res(npar + k); - chi2_out(row, col) = res(2 * npar); - } else { - chi2_out(row, col) = res(npar); - } - } - } - }; - - auto tasks = split_task(0, static_cast(y.shape(0)), n_threads); - RunInParallel(process, tasks); -} - -// ============================================================================ -// Explicit instantiations for all supported model types -// ============================================================================ - -// NOLINTBEGIN -#define AARE_INSTANTIATE_FIT(Model) \ - 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); - -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) - -#undef AARE_INSTANTIATE_FIT -// NOLINTEND - -} // namespace aare diff --git a/src/Fit.test.cpp b/src/Fit.test.cpp deleted file mode 100644 index 2b0fd1bb..00000000 --- a/src/Fit.test.cpp +++ /dev/null @@ -1,81 +0,0 @@ -// SPDX-License-Identifier: MPL-2.0 -#include "aare/Fit.hpp" -#include "aare/FitModel.hpp" -#include "aare/Models.hpp" - -#include -#include -#include - -namespace { - -constexpr ssize_t n_points = 61; - -void fill_gaussian(aare::NDArray &x, aare::NDArray &y, - double amplitude, double mean, double sigma) { - for (ssize_t i = 0; i < n_points; ++i) { - x(i) = -6.0 + 0.2 * static_cast(i); - const double z = (x(i) - mean) / sigma; - y(i) = amplitude * std::exp(-0.5 * z * z); - } -} - -} // namespace - -TEST_CASE("Minuit2 fits weighted and unweighted Gaussian data", "[fit]") { - aare::NDArray x({n_points}); - aare::NDArray y({n_points}); - aare::NDArray y_err({n_points}, 1.0); - fill_gaussian(x, y, 120.0, 0.8, 1.3); - - const aare::FitModel unweighted_model; - const auto unweighted = - aare::fit_pixel(unweighted_model, x.view(), y.view()); - - REQUIRE(unweighted.size() == 4); - CHECK(unweighted(0) == Catch::Approx(120.0).epsilon(1e-5)); - CHECK(unweighted(1) == Catch::Approx(0.8).epsilon(1e-5)); - CHECK(unweighted(2) == Catch::Approx(1.3).epsilon(1e-5)); - CHECK(unweighted(3) == Catch::Approx(0.0).margin(1e-4)); - - const aare::FitModel weighted_model(0, 100, 0.5, - true); - const auto weighted = - aare::fit_pixel(weighted_model, x.view(), y.view(), y_err.view()); - - REQUIRE(weighted.size() == 7); - CHECK(weighted(0) == Catch::Approx(120.0).epsilon(1e-5)); - CHECK(weighted(1) == Catch::Approx(0.8).epsilon(1e-5)); - CHECK(weighted(2) == Catch::Approx(1.3).epsilon(1e-5)); - CHECK(weighted(6) == Catch::Approx(0.0).margin(1e-4)); -} - -TEST_CASE("Minuit2 fits a Gaussian data cube in parallel", "[fit]") { - aare::NDArray x({n_points}); - aare::NDArray values({n_points}); - fill_gaussian(x, values, 80.0, -0.6, 0.9); - - aare::NDArray y({2, 2, n_points}); - for (ssize_t row = 0; row < 2; ++row) { - for (ssize_t col = 0; col < 2; ++col) { - for (ssize_t i = 0; i < n_points; ++i) - y(row, col, i) = values(i); - } - } - - aare::NDArray par({2, 2, 3}); - aare::NDArray chi2({2, 2}); - const aare::FitModel model; - - aare::fit_3d(model, x.view(), y.view(), aare::NDView{}, - par.view(), aare::NDView{}, chi2.view(), 2); - - for (ssize_t row = 0; row < 2; ++row) { - for (ssize_t col = 0; col < 2; ++col) { - CHECK(par(row, col, 0) == Catch::Approx(80.0).epsilon(1e-5)); - CHECK(par(row, col, 1) == Catch::Approx(-0.6).epsilon(1e-5)); - CHECK(par(row, col, 2) == Catch::Approx(0.9).epsilon(1e-5)); - CHECK(chi2(row, col) == Catch::Approx(0.0).margin(1e-8)); - } - } -} diff --git a/src/fit/Chi2.hpp b/src/fit/Chi2.hpp new file mode 100644 index 00000000..27fcf1af --- /dev/null +++ b/src/fit/Chi2.hpp @@ -0,0 +1,253 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once +#include "aare/Models.hpp" +#include "aare/NDView.hpp" +#include +#include +#include +#include +#include +#include + +namespace aare { + +namespace func { + +/** + * @brief Chi-squared objective of a 1D model and its analytic derivatives. + * + * Shared by the Minuit2 FCN adapters below. + * + * chi2(par) = sum_i w_i * (y_i - f(x_i; par))^2 + * + * with w_i = 1/sigma_i^2 for weighted fits (points with sigma_i == 0 are + * skipped) and w_i = 1 otherwise. + * + * @tparam Model A model struct that satisfies: + * - static constexpr std::size_t npar; + * - static double eval(double x, const std::vector& par); + * - static void eval_and_grad(double x, const std::vector& par, + * double& f, std::array& g); + * - static bool is_valid(const std::vector& par); + * + * Derivatives, with df_i/dp_k taken from Model::eval_and_grad: + * + * d(chi2)/dp_k = -2 * sum_i w_i * (y_i - f_i) * df_i/dp_k + * d2(chi2)/dp_j dp_k ~= 2 * sum_i w_i * df_i/dp_j * df_i/dp_k + * + * The second line drops the residual-weighted second derivatives of the + * model, which is the Fumili (Gauss-Newton) approximation of the Hessian. + * + * Invalid model parameters do not throw; they return a large penalty (with + * a zero gradient and unit curvature) so the minimizer can remain in control. + * + * @throws std::invalid_argument if par.size() != Model::npar. + */ +template class Chi2Data1D { + public: + static constexpr std::size_t npar = Model::npar; + /** @brief Number of elements of the packed symmetric Hessian. */ + static constexpr std::size_t nhess = npar * (npar + 1) / 2; + static constexpr double invalid_penalty = 1e20; + + /** @brief An empty y_err view selects an unweighted fit. */ + Chi2Data1D(NDView x, NDView y, + NDView y_err = NDView{}) + : x_(x), y_(y), s_(y_err), weighted_(y_err.size() > 0) {} + + double value(const std::vector &par) const { + check_size(par); + if (!Model::is_valid(par)) + return invalid_penalty; + return weighted_ ? value_impl(par) : value_impl(par); + } + + std::vector gradient(const std::vector &par) const { + check_size(par); + std::vector grad(npar, 0.0); + if (!Model::is_valid(par)) + return grad; + if (weighted_) + gradient_impl(par, grad); + else + gradient_impl(par, grad); + return grad; + } + + /** + * @brief Compute chi2, its gradient and the Gauss-Newton Hessian in a + * single pass over the data. + * + * @param grad Receives d(chi2)/dp_k, resized to npar. + * @param hess Receives the packed symmetric Hessian in the layout used by + * ROOT::Minuit2::FumiliFCNBase: element (j, k) with j <= k is + * stored at index j + k * (k + 1) / 2. Resized to nhess. + * @return chi2 at par. + */ + double evaluate_all(const std::vector &par, + std::vector &grad, + std::vector &hess) const { + check_size(par); + grad.assign(npar, 0.0); + hess.assign(nhess, 0.0); + if (!Model::is_valid(par)) { + for (std::size_t k = 0; k < npar; ++k) + hess[k + k * (k + 1) / 2] = 1.0; + return invalid_penalty; + } + return weighted_ ? evaluate_all_impl(par, grad, hess) + : evaluate_all_impl(par, grad, hess); + } + + private: + static void check_size(const std::vector &par) { + if (par.size() != npar) { + throw std::invalid_argument( + "Chi2Data1D: wrong parameter vector size."); + } + } + + template + double value_impl(const std::vector &par) const { + double chi2 = 0.0; + for (ssize_t i = 0; i < x_.size(); ++i) { + if constexpr (Weighted) { + const double si = s_[i]; + if (si == 0.0) + continue; + const double r_i = y_[i] - Model::eval(x_[i], par); + chi2 += (r_i * r_i) / (si * si); + } else { + const double r_i = y_[i] - Model::eval(x_[i], par); + chi2 += r_i * r_i; + } + } + return chi2; + } + + template + void gradient_impl(const std::vector &par, + std::vector &grad) const { + std::array df{}; + double f_i = 0.0; + for (ssize_t i = 0; i < x_.size(); ++i) { + double c = 0.0; + if constexpr (Weighted) { + const double si = s_[i]; + if (si == 0.0) + continue; + Model::eval_and_grad(x_[i], par, f_i, df); + c = -2.0 * (y_[i] - f_i) / (si * si); + } else { + Model::eval_and_grad(x_[i], par, f_i, df); + c = -2.0 * (y_[i] - f_i); + } + for (std::size_t k = 0; k < npar; ++k) + grad[k] += c * df[k]; + } + } + + template + double evaluate_all_impl(const std::vector &par, + std::vector &grad, + std::vector &hess) const { + std::array df{}; + double f_i = 0.0; + double chi2 = 0.0; + for (ssize_t i = 0; i < x_.size(); ++i) { + double w = 1.0; + if constexpr (Weighted) { + const double si = s_[i]; + if (si == 0.0) + continue; + w = 1.0 / (si * si); + } + Model::eval_and_grad(x_[i], par, f_i, df); + const double r_i = y_[i] - f_i; + chi2 += w * r_i * r_i; + const double c = -2.0 * w * r_i; + for (std::size_t k = 0; k < npar; ++k) { + grad[k] += c * df[k]; + const double h_k = 2.0 * w * df[k]; + double *column = hess.data() + k * (k + 1) / 2; + for (std::size_t j = 0; j <= k; ++j) + column[j] += h_k * df[j]; + } + } + return chi2; + } + + NDView x_; + NDView y_; + NDView s_; + bool weighted_; +}; + +/** + * @brief Chi-squared FCN with analytic gradient, used by Migrad. + * + * Providing the analytic gradient avoids the 2*npar extra function + * evaluations per Minuit step that finite differences would need. + */ +template +class Chi2Model1DGrad : public ROOT::Minuit2::FCNGradientBase { + public: + /** @brief An empty y_err view selects an unweighted fit. */ + Chi2Model1DGrad(NDView x, NDView y, + NDView y_err = NDView{}) + : data_(x, y, y_err) {} + + ~Chi2Model1DGrad() override = default; + + double operator()(const std::vector &par) const override { + return data_.value(par); + } + + std::vector + Gradient(const std::vector &par) const override { + return data_.gradient(par); + } + + /** @brief Error definition: 1.0 for chi-squared (delta_chi2 = 1 -> + * 1-sigma). */ + double Up() const override { return 1.0; } + + private: + Chi2Data1D data_; +}; + +/** + * @brief Chi-squared FCN for Fumili. + * + * EvaluateAll fills the value, gradient and Gauss-Newton Hessian cached by + * FumiliFCNBase in one pass over the data, so Fumili never needs numerical + * derivatives. + */ +template +class Chi2Model1DFumili : public ROOT::Minuit2::FumiliFCNBase { + public: + /** @brief An empty y_err view selects an unweighted fit. */ + Chi2Model1DFumili(NDView x, NDView y, + NDView y_err = NDView{}) + : ROOT::Minuit2::FumiliFCNBase(static_cast(Model::npar)), + data_(x, y, y_err) {} + + double operator()(const std::vector &par) const override { + return data_.value(par); + } + + /** @brief Error definition: 1.0 for chi-squared (delta_chi2 = 1 -> + * 1-sigma). */ + double Up() const override { return 1.0; } + + void EvaluateAll(const std::vector &par) override { + SetFCNValue(data_.evaluate_all(par, Gradient(), Hessian())); + } + + private: + Chi2Data1D data_; +}; + +} // namespace func + +} // namespace aare diff --git a/src/fit/DampedGaussNewton.hpp b/src/fit/DampedGaussNewton.hpp new file mode 100644 index 00000000..fdf1670b --- /dev/null +++ b/src/fit/DampedGaussNewton.hpp @@ -0,0 +1,597 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include +#include +#include +#include +#include + +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; + int calls = 0; + int iterations = 0; + 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. + */ +template +bool gauss_newton_errors(const Mat &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)) + return false; + 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)); + ++i; + } + return true; +} + +/** + * @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 + * 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. + * - 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 + * points outward is frozen for that iteration. + * - Convergence in Minuit's convention: the iteration stops when the + * estimated distance to the minimum of chi2 (EDM, from the Gauss-Newton + * Hessian) drops below 0.002 * tolerance. A start point that passes + * already still takes one undamped step, see converged(). When the damped + * step becomes negligible for every parameter before that, the iteration + * has stalled and the fit is valid only if the EDM test passes at that + * point. + * - 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. + * - With polish, a converged fit takes one more undamped Gauss-Newton step + * and keeps it when it lowers chi2, see polish_step(). + * + * After fit() the final point, its Normal and the frozen parameters remain + * available for the error computation. + */ +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 + }; + + 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; + }; + + template + MinimizerResult fit(Eval &eval, const Vec &start, const Vec &lower, + const Vec &upper, const std::array &fixed, + const Options &opt) { + 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_; + free_[nfree_++] = k; + } + } + for (int a = 0; a < nfree_; ++a) { + q_[a] = p_[free_[a]]; + active_[a] = false; + } + + MinimizerResult res; + if (!evaluate(eval, q_, cur_)) + 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; + + 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; + } + } + 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. */ + const Normal &normal() const { return cur_; } + + /** + * @brief True for a parameter without an error: fixed, frozen on a limit, + * or ending exactly on one. The frozen flag alone misses the last case, + * because it needs the gradient to point outward and that gradient is + * zero when the unconstrained minimum coincides with the limit. + */ + bool skipped(int k) const { + return fixed_[k] || (free_pos_[k] >= 0 && active_[free_pos_[k]]) || + p_[k] <= lower_[k] || p_[k] >= upper_[k]; + } + + /** @brief Gauss-Newton errors at the final point, see + * gauss_newton_errors(). */ + void errors(double *err_out) const { + std::array skip; + for (int k = 0; k < NP; ++k) + skip[k] = skipped(k); + gauss_newton_errors(cur_.H, skip, err_out); + } + + private: + template + bool evaluate(Eval &eval, const Vec &q, Normal &out) { + for (int a = 0; a < nfree_; ++a) + p_[free_[a]] = q[a]; + return eval(p_, out); + } + + // Gradient and Hessian of the free parameters at the current point. + void extract_free() { + 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]; + } + gmax_ = std::max(gmax_, std::abs(g_[a])); + } + } + + // A parameter sitting on a limit with the descent direction pointing + // outward is frozen for this iteration. Position and gradient are both + // taken at the current point q_, so a parameter stays frozen at the end + // only where the constrained minimum lies on its limit. An earlier + // version tested the last evaluated point p_, a rejected trial, so that a + // trial that came off a limit unfroze the parameter: testing q_ had left + // more FallingScurve and GaussianChargeSharing fits stuck with the width + // 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(). + void mark_active() { + for (int a = 0; a < nfree_; ++a) { + const int k = free_[a]; + active_[a] = (q_[a] <= lower_[k] && g_[a] > 0.0) || + (q_[a] >= upper_[k] && g_[a] < 0.0); + } + } + + // True when the damped step is negligible for every parameter that is + // not frozen: below 1e-8 of its value plus its own scale 1/sqrt(H_aa), + // the parameter change that raises chi2 by about one. Judging each + // parameter on its own scale keeps a large parameter (an amplitude of + // 1e10, a peak position at 1e9) from hiding the steps of the others. + bool negligible_step() const { + for (int a = 0; a < nfree_; ++a) { + if (active_[a]) + continue; + const double scale = + 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; + } + return true; + } + + // 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 + // many accepted steps, each of which passes the test against its + // predecessor. + bool degenerate(const Normal &trial) const { + for (int a = 0; a < nfree_; ++a) { + 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)) + 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 + // does, so it still takes the undamped step of polish_step(), as Migrad + // always iterates at least once. Without the budget for that step the + // 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) { + const bool first = accepted_ == 0 && nfree_ > 0; + if (!first && !polish) + return true; + if (res.calls >= max_calls) + return !first; + polish_step(eval, res, reject_degenerate); + return true; + } + + // One undamped Gauss-Newton step from a point that passed the EDM test, + // kept when it lowers F. The test stops the iteration within the + // tolerance, but a solver that converges in few iterations stops while + // the damping is still noticeable, so its last accepted step fell short + // 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) { + 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_); + 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 { + double mean_diag = 0.0; + for (int a = 0; a < nfree_; ++a) + mean_diag += jtj_[a][a]; + mean_diag = std::max(mean_diag / std::max(nfree_, 1), 1e-300); + int m = 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; + } + return m; + } + + // 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)) + 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; + } + return predicted_ > 0.0; + } + + int nfree_ = 0; + int accepted_ = 0; // trial points accepted in this fit + std::array free_{}; + std::array free_pos_{}; + std::array fixed_{}; + std::array active_{}; + Vec lower_{}; + Vec upper_{}; + Vec p_{}; + Vec q_{}; + Vec q_new_{}; + Vec q_alt_{}; + 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; + double gmax_ = 0.0; + double predicted_ = 0.0; + double bound_ = 0.0; + Normal cur_; + Normal trial_; + Normal alt_; +}; + +} // namespace aare::detail diff --git a/src/fit/Fit.cpp b/src/fit/Fit.cpp new file mode 100644 index 00000000..8b34c177 --- /dev/null +++ b/src/fit/Fit.cpp @@ -0,0 +1,337 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/Fit.hpp" +#include "FitHelpers.hpp" +#include "FitMinuit2.hpp" +#include "LevenbergMarquardt.hpp" +#include "VariableProjection.hpp" +#include "aare/Models.hpp" +#include "aare/utils/par.hpp" +#include "aare/utils/task.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +namespace aare { +// ============================================================================ +// FitModel — method definitions +// ============================================================================ + +template +FitModel::FitModel(unsigned int strategy, unsigned int max_calls, + double tolerance, bool compute_errors, + Minimizer minimizer) + : strategy_(strategy), max_calls_(max_calls), tolerance_(tolerance), + compute_errors_(compute_errors), minimizer_(minimizer) { + for (std::size_t i = 0; i < npar; ++i) { + lower_[i] = Model::param_info[i].default_lo; + upper_[i] = Model::param_info[i].default_hi; + } +} + +template +void FitModel::check_index(unsigned int idx) const { + if (idx >= npar) + throw std::out_of_range("FitModel: parameter index " + + std::to_string(idx) + " out of range"); +} + +template +unsigned int FitModel::checked_index(const std::string &name) const { + for (std::size_t i = 0; i < npar; ++i) { + if (name == Model::param_info[i].name) + return static_cast(i); + } + throw std::runtime_error("FitModel: unknown parameter name '" + name + "'"); +} + +template +void FitModel::SetParLimits(unsigned int idx, double lo, double hi) { + check_index(idx); + if (!(lo < hi)) + throw std::runtime_error( + "FitModel: lower limit must be below upper limit for parameter '" + + GetParName(idx) + "'"); + lower_[idx] = lo; + upper_[idx] = hi; +} + +template +void FitModel::FixParameter(unsigned int idx, double val) { + SetParameter(idx, val); + user_fixed_[idx] = true; +} + +template +void FitModel::ReleaseParameter(unsigned int idx) { + check_index(idx); + user_fixed_[idx] = false; +} + +template +void FitModel::ReleaseParameter(const std::string &name) { + ReleaseParameter(checked_index(name)); +} + +template +void FitModel::SetParameter(unsigned int idx, double val) { + check_index(idx); + value_[idx] = val; + user_start_[idx] = true; +} + +template +void FitModel::SetParameter(const std::string &name, double val) { + SetParameter(checked_index(name), val); +} + +template +void FitModel::FixParameter(const std::string &name, double val) { + FixParameter(checked_index(name), val); +} + +template +void FitModel::SetParLimits(const std::string &name, double lo, + double hi) { + SetParLimits(checked_index(name), lo, hi); +} + +template +std::string FitModel::GetParName(unsigned int idx) const { + check_index(idx); + return Model::param_info[idx].name; +} + +template +std::vector FitModel::GetParNames() const { + std::vector names; + names.reserve(npar); + for (std::size_t i = 0; i < npar; ++i) + names.emplace_back(Model::param_info[i].name); + return names; +} + +// ============================================================================ +// fit_pixel / fit_3d +// ============================================================================ + +namespace detail { + +struct NotSeparable {}; + +/** + * @brief Fits pixels with the minimizer selected by the model. + * + * Keep one instance per thread: the buffers of the built-in minimizers are + * reused between pixels so that a data cube is fitted without per-pixel + * allocations. The Minuit2 minimizers create their state per pixel. + */ +template class PixelFitter { + public: + /** + * @brief Fit one pixel. + * + * @param par_out Receives npar fitted values; zeros on failure. + * @param err_out Receives npar errors when non-null; zeros on failure. + * @param chi2 Receives the chi-squared at the minimum; 0 on failure. + * @return true for a valid minimum. + */ + bool fit(const FitModel &model, NDView x, + NDView y, NDView y_err, double *par_out, + double *err_out, double &chi2) { + // Samples with a zero error are ignored. They are dropped here, + // before the start estimate and the minimizer: a masked sample may + // hold NaN, and a zero weight does not remove a NaN from a sum. + if (y_err.size() > 0 && + std::find(y_err.begin(), y_err.end(), 0.0) != y_err.end()) { + keep_x_.clear(); + keep_y_.clear(); + keep_err_.clear(); + for (ssize_t i = 0; i < y.size(); ++i) { + if (y_err[i] == 0.0) + continue; + keep_x_.push_back(x[i]); + keep_y_.push_back(y[i]); + keep_err_.push_back(y_err[i]); + } + const auto n = static_cast(keep_x_.size()); + x = NDView(keep_x_.data(), {n}); + y = NDView(keep_y_.data(), {n}); + y_err = NDView(keep_err_.data(), {n}); + } + // Fewer samples than free parameters leave the fit undetermined, + // and the start estimates read at least two samples. + ssize_t n_free = 0; + for (std::size_t k = 0; k < Model::npar; ++k) + n_free += model.is_user_fixed(static_cast(k)) ? 0 : 1; + if (y.size() < std::max(n_free, 2)) + 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. + 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; + } + } + return fit_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); + } + + 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])); + 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. + static bool fail(double *par_out, double *err_out, double &chi2) { + std::fill(par_out, par_out + Model::npar, 0.0); + if (err_out) + std::fill(err_out, err_out + Model::npar, 0.0); + chi2 = 0.0; + return false; + } + + // The samples of the current pixel that have a nonzero error, see fit(). + std::vector keep_x_; + std::vector keep_y_; + std::vector keep_err_; + LevenbergMarquardt lm_; + std::conditional_t::value, + VariableProjection, NotSeparable> + vp_; +}; + +} // namespace detail + +template +NDArray fit_pixel(const FitModel &model, NDView x, + NDView y, NDView y_err) { + constexpr std::size_t npar = Model::npar; + detail::check_fit_pixel_shapes(x, y, y_err); + const bool want_errors = model.compute_errors(); + const auto result_size = + static_cast(want_errors ? (2 * npar + 1) : (npar + 1)); + NDArray result({result_size}, 0.0); + + detail::PixelFitter fitter; + double chi2 = 0.0; + double *par = result.data(); + double *err = want_errors ? result.data() + npar : nullptr; + if (!fitter.fit(model, x, y, y_err, par, err, chi2)) + return result; // all zeros + + result(result_size - 1) = chi2; + return result; +} + +template +NDArray fit_pixel(const FitModel &model, NDView x, + NDView y) { + return fit_pixel(model, x, y, NDView{}); +} + +template +void fit_3d(const FitModel &model, NDView x, + NDView y, NDView y_err, + NDView par_out, NDView err_out, + NDView chi2_out, int n_threads) { + constexpr std::size_t npar = Model::npar; + detail::check_fit_3d_shapes(x, y, y_err, par_out, err_out, chi2_out); + + const bool has_errors = (y_err.size() > 0); + const bool want_par_errors = (err_out.size() > 0) && model.compute_errors(); + + auto process = [&](ssize_t first_row, ssize_t last_row) { + detail::PixelFitter fitter; // buffers reused within the task + std::array par{}; + std::array err{}; + for (ssize_t row = first_row; row < last_row; row++) { + for (ssize_t col = 0; col < y.shape(1); col++) { + NDView values(&y(row, col, 0), {y.shape(2)}); + NDView errors = + has_errors ? NDView(&y_err(row, col, 0), + {y_err.shape(2)}) + : NDView{}; + + double chi2 = 0.0; + fitter.fit(model, x, values, errors, par.data(), + want_par_errors ? err.data() : nullptr, chi2); + + for (std::size_t k = 0; k < npar; ++k) + par_out(row, col, static_cast(k)) = par[k]; + if (want_par_errors) { + for (std::size_t k = 0; k < npar; ++k) + err_out(row, col, static_cast(k)) = err[k]; + } + chi2_out(row, col) = chi2; + } + } + }; + + auto tasks = split_task(0, static_cast(y.shape(0)), n_threads); + RunInParallel(process, tasks); +} + +// ============================================================================ +// Explicit instantiations for all supported model types +// ============================================================================ + +// NOLINTBEGIN +#define AARE_INSTANTIATE_FIT(Model) \ + 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); + +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) + +#undef AARE_INSTANTIATE_FIT +// NOLINTEND + +} // namespace aare diff --git a/src/fit/Fit.test.cpp b/src/fit/Fit.test.cpp new file mode 100644 index 00000000..992500e8 --- /dev/null +++ b/src/fit/Fit.test.cpp @@ -0,0 +1,1308 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/Fit.hpp" +#include "aare/FitModel.hpp" +#include "aare/Models.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using aare::Minimizer; +using aare::NDArray; +using aare::NDView; +using Catch::Approx; + +namespace { + +constexpr ssize_t n_points = 61; + +struct Data { + NDArray x; + NDArray y; + NDArray y_err; +}; + +void fill_gaussian(NDArray &x, NDArray &y, + double amplitude, double mean, double sigma) { + for (ssize_t i = 0; i < n_points; ++i) { + x(i) = -6.0 + 0.2 * static_cast(i); + const double z = (x(i) - mean) / sigma; + y(i) = amplitude * std::exp(-0.5 * z * z); + } +} + +// Gaussian with A = 120, mu = 0.8, sigma = 1.3 and unit errors, optionally +// with deterministic pseudo-noise. +Data gaussian_data(bool add_noise = false) { + Data d{NDArray({n_points}), NDArray({n_points}), + NDArray({n_points}, 1.0)}; + fill_gaussian(d.x, d.y, 120.0, 0.8, 1.3); + if (add_noise) { + for (ssize_t i = 0; i < n_points; ++i) + d.y(i) += 3.0 * std::sin(7.3 * static_cast(i)); + } + return d; +} + +// S-curve on 100 points in [0, 100] with unit errors, optionally with +// deterministic pseudo-noise and errors of 30. +template +Data scurve_data(const std::vector &par, bool add_noise = false) { + constexpr ssize_t n = 100; + Data d{NDArray({n}), NDArray({n}), + NDArray({n}, 1.0)}; + for (ssize_t i = 0; i < n; ++i) { + d.x(i) = 100.0 * static_cast(i) / static_cast(n - 1); + d.y(i) = Model::eval(d.x(i), par); + if (add_noise) { + d.y(i) += 30.0 * std::sin(7.3 * static_cast(i)); + d.y_err(i) = 30.0; + } + } + return d; +} + +const char *name(Minimizer m) { + switch (m) { + case Minimizer::Migrad: + return "Migrad"; + case Minimizer::Fumili: + return "Fumili"; + case Minimizer::LevenbergMarquardt: + return "LevenbergMarquardt"; + case Minimizer::VarPro: + return "VarPro"; + } + return "unknown"; +} + +bool uses_minuit2(Minimizer m) { + return m != Minimizer::LevenbergMarquardt && m != Minimizer::VarPro; +} + +template +aare::FitModel make_model(Minimizer minimizer, bool errors = false, + unsigned int max_calls = 100) { + return aare::FitModel(0, max_calls, 0.5, errors, minimizer); +} + +// Minuit2 reaches an active limit through its parameter transformation; +// the built-in solver lands on it exactly. +double limit_margin(Minimizer m) { return uses_minuit2(m) ? 1e-6 : 1e-9; } + +// Relative agreement with the true parameters: the Minuit2 minimizers stop +// at Minuit's EDM tolerance, the built-in solver lands tighter. +double truth_tolerance(Minimizer m) { return uses_minuit2(m) ? 1e-3 : 1e-4; } + +// Hesse errors at a point forced away from the free minimum differ from the +// Gauss-Newton estimates of Fumili and the built-in solver. +double misfit_error_tolerance(Minimizer m) { + return uses_minuit2(m) ? 0.1 : 5e-2; +} + +template +bool on_limit(const aare::FitModel &model, unsigned int idx, + double value) { + const double margin = 1e-6 * std::max(1.0, std::abs(value)); + return std::abs(value - model.lower_limit(idx)) <= margin || + std::abs(value - model.upper_limit(idx)) <= margin; +} + +/** + * Run one configuration through every minimizer and compare Fumili and the + * built-in solver against Migrad: parameters and chi2 within tol (relative), + * unless the candidate found a clearly lower chi2. Errors are compared for + * free parameters that do not end on a limit with tolerance err_tol. + */ +template +void check_minimizers_agree(const std::string &label, Data d, bool weighted, + bool errors, unsigned int max_calls, double tol, + Setup setup, double err_tol = 0.0) { + if (err_tol == 0.0) + err_tol = 10.0 * tol; + constexpr std::size_t npar = Model::npar; + + const auto run = [&](Minimizer minimizer) { + aare::FitModel model(0, max_calls, 0.5, errors, minimizer); + setup(model); + return weighted ? aare::fit_pixel(model, d.x.view(), d.y.view(), + d.y_err.view()) + : aare::fit_pixel(model, d.x.view(), d.y.view()); + }; + aare::FitModel reference_model(0, max_calls, 0.5, errors, + Minimizer::Migrad); + setup(reference_model); + const auto ref = run(Minimizer::Migrad); + + for (const auto candidate : + {Minimizer::Fumili, Minimizer::LevenbergMarquardt, + Minimizer::VarPro}) { + INFO(label << ": " << name(candidate) << " vs Migrad"); + const auto res = run(candidate); + REQUIRE(res.size() == ref.size()); + + const auto ichi = static_cast(errors ? 2 * npar : npar); + const double chi_ref = ref(ichi); + const double chi_res = res(ichi); + if (chi_ref > 0.0 && chi_res > 0.0 && + chi_ref - chi_res > tol * std::max(chi_ref, 1.0)) { + WARN(label + ": " + name(candidate) + " found a lower chi2 (" + + std::to_string(chi_res) + " vs " + std::to_string(chi_ref) + + ")"); + continue; + } + double scale = 0.0; + for (std::size_t k = 0; k < npar; ++k) + scale = std::max(scale, std::abs(ref(static_cast(k)))); + for (std::size_t k = 0; k < npar; ++k) { + const auto i = static_cast(k); + const double denom = + std::max({std::abs(ref(i)), 1e-2 * scale, 1e-12}); + INFO("parameter " + << reference_model.GetParName(static_cast(k))); + CHECK(std::abs(res(i) - ref(i)) / denom <= tol); + } + CHECK(std::abs(chi_res - chi_ref) / std::max(std::abs(chi_ref), 1.0) <= + tol); + if (!errors) + continue; + for (std::size_t k = 0; k < npar; ++k) { + const auto idx = static_cast(k); + const auto i = static_cast(k); + if (reference_model.is_user_fixed(idx) || + on_limit(reference_model, idx, ref(i))) + continue; + const auto ie = static_cast(npar + k); + INFO("error of " << reference_model.GetParName(idx)); + CHECK(std::abs(res(ie) - ref(ie)) / + std::max(std::abs(ref(ie)), 1e-12) <= + err_tol); + } + } +} + +} // namespace + +TEST_CASE("Fit unweighted and weighted noise-free Gaussian data", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + + const auto unweighted_model = make_model(minimizer); + const auto unweighted = + aare::fit_pixel(unweighted_model, d.x.view(), d.y.view()); + + REQUIRE(unweighted.size() == 4); + CHECK(unweighted(0) == Approx(120.0).epsilon(1e-4).margin(1e-6)); + CHECK(unweighted(1) == Approx(0.8).epsilon(1e-4).margin(1e-6)); + CHECK(unweighted(2) == Approx(1.3).epsilon(1e-4).margin(1e-6)); + CHECK(unweighted(3) == Approx(0.0).margin(1e-3)); + + const auto weighted_model = + make_model(minimizer, true); + const auto weighted = + aare::fit_pixel(weighted_model, d.x.view(), d.y.view(), d.y_err.view()); + + REQUIRE(weighted.size() == 7); + CHECK(weighted(0) == Approx(120.0).epsilon(1e-4).margin(1e-6)); + CHECK(weighted(1) == Approx(0.8).epsilon(1e-4).margin(1e-6)); + CHECK(weighted(2) == Approx(1.3).epsilon(1e-4).margin(1e-6)); + // At zero residual the Hesse, Fumili and Gauss-Newton errors coincide. + CHECK(weighted(3) == Approx(0.361).epsilon(2e-2)); + CHECK(weighted(4) == Approx(0.00451).epsilon(2e-2)); + CHECK(weighted(5) == Approx(0.00451).epsilon(2e-2)); + CHECK(weighted(6) == Approx(0.0).margin(1e-3)); +} + +TEST_CASE("Fit weighted Pol1 errors match the analytic formula", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t n = 20; + NDArray x({n}); + NDArray y({n}); + NDArray y_err({n}); + double S = 0.0, Sx = 0.0, Sxx = 0.0; + for (ssize_t i = 0; i < n; ++i) { + x(i) = static_cast(i); + y(i) = 2.0 + 0.5 * x(i); + y_err(i) = 1.0 + 0.1 * static_cast(i); + const double w = 1.0 / (y_err(i) * y_err(i)); + S += w; + Sx += w * x(i); + Sxx += w * x(i) * x(i); + } + const double delta = S * Sxx - Sx * Sx; + + const auto model = make_model(minimizer, true); + const auto res = aare::fit_pixel(model, x.view(), y.view(), y_err.view()); + + REQUIRE(res.size() == 5); + CHECK(res(0) == Approx(2.0).epsilon(1e-6)); + CHECK(res(1) == Approx(0.5).epsilon(1e-6)); + CHECK(res(2) == Approx(std::sqrt(Sxx / delta)).epsilon(1e-6)); + CHECK(res(3) == Approx(std::sqrt(S / delta)).epsilon(1e-6)); + CHECK(res(4) == Approx(0.0).margin(1e-8)); +} + +TEST_CASE("Fit with a fixed parameter", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + + SECTION("fixed at a wrong value") { + auto model = make_model(minimizer, true); + model.FixParameter("sigma", 1.5); + REQUIRE(model.is_user_fixed(2)); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + REQUIRE(res.size() == 7); + CHECK(res(2) == 1.5); + CHECK(res(0) == Approx(111.15).epsilon(1e-3)); + CHECK(res(1) == Approx(0.8).epsilon(1e-3)); + CHECK(res(3) == + Approx(0.274).epsilon(misfit_error_tolerance(minimizer))); + CHECK(res(5) == 0.0); // fixed parameters report no error + CHECK(res(6) == Approx(1684.25).epsilon(1e-3)); + } + + SECTION("fixed at the true value, then released") { + auto model = make_model(minimizer, true); + model.FixParameter(2, 1.3); + auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(0) == Approx(120.0).epsilon(1e-4).margin(1e-6)); + CHECK(res(1) == Approx(0.8).epsilon(1e-4).margin(1e-6)); + CHECK(res(2) == 1.3); + CHECK(res(5) == 0.0); + CHECK(res(6) == Approx(0.0).margin(1e-3)); + + model.ReleaseParameter("sigma"); + REQUIRE_FALSE(model.is_user_fixed(2)); + res = aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(2) == Approx(1.3).epsilon(1e-4).margin(1e-6)); + CHECK(res(5) == Approx(0.00451).epsilon(2e-2)); + } +} + +TEST_CASE("Fit with parameter limits", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + const double margin = limit_margin(minimizer); + const double err_tol = misfit_error_tolerance(minimizer); + + SECTION("active upper limit on mu") { + auto model = make_model(minimizer); + model.SetParLimits("mu", -5.0, 0.5); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + CHECK(res(1) == Approx(0.5).margin(margin)); + CHECK(res(0) == Approx(116.88).epsilon(1e-3)); + CHECK(res(2) == Approx(1.335).epsilon(1e-3)); + CHECK(res(3) == Approx(4301.95).epsilon(1e-3)); + } + + SECTION("active upper limit on sigma, with errors") { + auto model = make_model(minimizer, true); + model.SetParLimits(2, 0.5, 1.0); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(2) == Approx(1.0).margin(margin)); + CHECK(res(1) == Approx(0.8).epsilon(1e-3)); + CHECK(res(3) == Approx(0.336).epsilon(err_tol)); + CHECK(res(5) == 0.0); // on a limit: no error reported + CHECK(res(6) == Approx(5550.61).epsilon(1e-3)); + } + + SECTION("fixed amplitude and an active limit on mu") { + auto model = make_model(minimizer, true); + model.FixParameter("A", 100.0); + model.SetParLimits("mu", 0.6, 0.7); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(0) == 100.0); + CHECK(res(1) == Approx(0.7).margin(margin)); + CHECK(res(2) == Approx(1.4603).epsilon(1e-3)); + CHECK(res(3) == 0.0); + CHECK(res(4) == 0.0); + CHECK(res(5) == Approx(0.00469).epsilon(err_tol)); + CHECK(res(6) == Approx(3631.01).epsilon(1e-3)); + } + + SECTION("inactive limits do not change the result") { + const auto free_model = make_model(minimizer); + const auto reference = + aare::fit_pixel(free_model, d.x.view(), d.y.view()); + + auto model = make_model(minimizer); + model.SetParLimits("mu", -5.0, 5.0); + model.SetParLimits("sigma", 0.5, + std::numeric_limits::infinity()); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + for (ssize_t k = 0; k < 3; ++k) + CHECK(res(k) == Approx(reference(k)).epsilon(1e-4).margin(1e-6)); + } +} + +TEST_CASE("Fit honours user start values", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + + SECTION("start value for mu") { + auto model = make_model(minimizer); + model.SetParameter("mu", 0.0); + REQUIRE(model.is_user_start(1)); + CHECK(model.value(1) == 0.0); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + CHECK(res(0) == Approx(120.0).epsilon(1e-4).margin(1e-6)); + CHECK(res(1) == Approx(0.8).epsilon(1e-4).margin(1e-6)); + CHECK(res(2) == Approx(1.3).epsilon(1e-4).margin(1e-6)); + } + + SECTION("negative Gaussian with all start values given") { + Data n{NDArray({n_points}), NDArray({n_points}), + NDArray({n_points}, 1.0)}; + fill_gaussian(n.x, n.y, -120.0, 0.8, 1.3); + auto model = make_model(minimizer); + model.SetParameter("A", -100.0); + model.SetParameter("mu", 0.5); + model.SetParameter("sigma", 1.0); + const auto res = aare::fit_pixel(model, n.x.view(), n.y.view()); + CHECK(res(0) == Approx(-120.0).epsilon(1e-4).margin(1e-6)); + CHECK(res(1) == Approx(0.8).epsilon(1e-4).margin(1e-6)); + CHECK(res(2) == Approx(1.3).epsilon(1e-4).margin(1e-6)); + } +} + +TEST_CASE("Fit reports failure as zeros when max_calls is exhausted", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + + // Start far from the minimum so that no minimizer converges in 3 calls. + auto model = make_model(minimizer, false, 3); + model.SetParameter("mu", -3.0); + model.SetParameter("sigma", 3.0); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + REQUIRE(res.size() == 4); + for (ssize_t k = 0; k < res.size(); ++k) + CHECK(res(k) == 0.0); + + auto model_err = make_model(minimizer, true, 3); + model_err.SetParameter("mu", -3.0); + model_err.SetParameter("sigma", 3.0); + const auto res_err = + aare::fit_pixel(model_err, d.x.view(), d.y.view(), d.y_err.view()); + REQUIRE(res_err.size() == 7); + for (ssize_t k = 0; k < res_err.size(); ++k) + CHECK(res_err(k) == 0.0); +} + +TEST_CASE("FitModel validates parameter names, indices and limits", "[fit]") { + aare::FitModel model; + CHECK(model.GetParNames() == std::vector{"A", "mu", "sigma"}); + CHECK(model.GetParName(2) == "sigma"); + CHECK(model.GetNpar() == 3); + CHECK(model.strategy() == 0); + CHECK(model.minimizer() == Minimizer::Migrad); + CHECK(model.lower_limit(2) == 1e-12); + CHECK(std::isinf(model.upper_limit(2))); + + CHECK_THROWS_AS(model.SetParameter("nope", 1.0), std::runtime_error); + CHECK_THROWS_AS(model.SetParameter(3, 1.0), std::out_of_range); + CHECK_THROWS_AS(model.GetParName(5), std::out_of_range); + CHECK_THROWS_AS(model.SetParLimits(0, 1.0, 1.0), std::runtime_error); + CHECK_THROWS_AS(model.FixParameter("nope", 1.0), std::runtime_error); + CHECK_NOTHROW(model.SetParLimits("mu", -1.0, 1.0)); + CHECK(model.lower_limit(1) == -1.0); + CHECK(model.upper_limit(1) == 1.0); + + model.SetMinimizer(Minimizer::LevenbergMarquardt); + CHECK(model.minimizer() == Minimizer::LevenbergMarquardt); +} + +TEST_CASE("fit_pixel rejects arrays of different sizes", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + const auto model = make_model(minimizer); + NDView short_y(d.y.data(), {n_points - 1}); + NDView short_err(d.y_err.data(), {n_points - 1}); + CHECK_THROWS_AS(aare::fit_pixel(model, d.x.view(), short_y), + std::runtime_error); + CHECK_THROWS_AS(aare::fit_pixel(model, d.x.view(), d.y.view(), short_err), + std::runtime_error); +} + +TEST_CASE("Copying a FitModel keeps the minimizer", "[fit]") { + const aare::FitModel original( + 0, 100, 0.5, false, Minimizer::LevenbergMarquardt); + const aare::FitModel copied(original); + aare::FitModel assigned; + assigned = original; + + CHECK(copied.minimizer() == Minimizer::LevenbergMarquardt); + CHECK(assigned.minimizer() == Minimizer::LevenbergMarquardt); +} + +TEST_CASE("Fit rising and falling S-curves", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const std::vector truth{10.0, 0.1, 50.0, 3.0, 1000.0, 2.0}; + + SECTION("rising, free fit with errors") { + auto d = scurve_data(truth); + const auto model = + make_model(minimizer, true, 500); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + REQUIRE(res.size() == 13); + for (std::size_t k = 0; k < truth.size(); ++k) { + CHECK(res(static_cast(k)) == + Approx(truth[k]).epsilon(1e-3)); + CHECK(res(static_cast(6 + k)) > 0.0); + } + CHECK(res(12) == Approx(0.0).margin(1e-3)); + } + + SECTION("rising, fixed slope and limited width") { + auto d = scurve_data(truth); + auto model = + make_model(minimizer, true, 500); + model.FixParameter("p1", 0.1); + model.SetParLimits("sigma", 1.0, 4.0); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(1) == 0.1); + CHECK(res(7) == 0.0); + CHECK(res(2) == Approx(50.0).epsilon(1e-3)); + CHECK(res(3) == Approx(3.0).epsilon(1e-3)); + CHECK(res(9) > 0.0); + } + + SECTION("rising, active limit on the threshold") { + auto d = scurve_data(truth); + auto model = + make_model(minimizer, false, 500); + model.SetParLimits("mu", 40.0, 48.0); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + CHECK(res(2) == Approx(48.0).margin(limit_margin(minimizer))); + CHECK(res(6) > 0.0); // chi2 of the constrained fit + } + + SECTION("falling") { + auto d = scurve_data(truth); + const auto model = + make_model(minimizer, false, 500); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + REQUIRE(res.size() == 7); + for (std::size_t k = 0; k < truth.size(); ++k) + CHECK(res(static_cast(k)) == + Approx(truth[k]).epsilon(1e-3)); + } +} + +TEST_CASE("Fit a Gaussian data cube in parallel", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + NDArray x({n_points}); + NDArray values({n_points}); + fill_gaussian(x, values, 80.0, -0.6, 0.9); + + NDArray y({2, 2, n_points}); + for (ssize_t row = 0; row < 2; ++row) { + for (ssize_t col = 0; col < 2; ++col) { + for (ssize_t i = 0; i < n_points; ++i) + y(row, col, i) = values(i); + } + } + + NDArray par({2, 2, 3}); + NDArray chi2({2, 2}); + const auto model = make_model(minimizer); + + aare::fit_3d(model, x.view(), y.view(), NDView{}, par.view(), + NDView{}, chi2.view(), 2); + + for (ssize_t row = 0; row < 2; ++row) { + for (ssize_t col = 0; col < 2; ++col) { + CHECK(par(row, col, 0) == Approx(80.0).epsilon(1e-4).margin(1e-6)); + CHECK(par(row, col, 1) == Approx(-0.6).epsilon(1e-4).margin(1e-6)); + CHECK(par(row, col, 2) == Approx(0.9).epsilon(1e-4).margin(1e-6)); + CHECK(chi2(row, col) == Approx(0.0).margin(1e-3)); + } + } + + CHECK_THROWS_AS(aare::fit_3d(model, NDView(x.data(), {10}), + y.view(), NDView{}, par.view(), + NDView{}, chi2.view(), 2), + std::runtime_error); +} + +TEST_CASE("Fit a data cube with errors matches per-pixel fits", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t rows = 4; + constexpr ssize_t cols = 3; + NDArray x({n_points}); + NDArray y({rows, cols, n_points}); + NDArray y_err({rows, cols, n_points}, 1.0); + NDArray tmp({n_points}); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + fill_gaussian(x, tmp, 80.0 + 10.0 * static_cast(row), + -0.6 + 0.3 * static_cast(col), + 0.9 + 0.1 * static_cast(row)); + for (ssize_t i = 0; i < n_points; ++i) + y(row, col, i) = tmp(i); + } + } + + NDArray par({rows, cols, 3}); + NDArray err({rows, cols, 3}); + NDArray chi2({rows, cols}); + const auto model = make_model(minimizer, true); + + aare::fit_3d(model, x.view(), y.view(), y_err.view(), par.view(), + err.view(), chi2.view(), 3); + + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + NDView yv(&y(row, col, 0), {n_points}); + NDView ev(&y_err(row, col, 0), {n_points}); + const auto single = aare::fit_pixel(model, x.view(), yv, ev); + // reused solver buffers must not leak state between pixels + for (ssize_t k = 0; k < 3; ++k) { + CHECK(par(row, col, k) == Approx(single(k)).margin(1e-12)); + CHECK(err(row, col, k) == Approx(single(3 + k)).margin(1e-12)); + CHECK(err(row, col, k) > 0.0); + } + CHECK(chi2(row, col) == Approx(single(6)).margin(1e-12)); + const double tol = truth_tolerance(minimizer); + CHECK(par(row, col, 0) == + Approx(80.0 + 10.0 * static_cast(row)) + .epsilon(tol) + .margin(1e-6)); + CHECK(par(row, col, 1) == + Approx(-0.6 + 0.3 * static_cast(col)) + .epsilon(tol) + .margin(1e-6)); + CHECK(par(row, col, 2) == + Approx(0.9 + 0.1 * static_cast(row)) + .epsilon(tol) + .margin(1e-6)); + } + } +} + +TEST_CASE("Fitting with every parameter fixed returns the fixed values", + "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + + auto model = make_model(minimizer, true); + model.FixParameter("A", 120.0); + model.FixParameter("mu", 0.8); + model.FixParameter("sigma", 1.3); + + const auto result = aare::fit_pixel(model, d.x.view(), d.y.view()); + + REQUIRE(result.size() == 7); + CHECK(result(0) == 120.0); + CHECK(result(1) == 0.8); + CHECK(result(2) == 1.3); + for (ssize_t k = 3; k < 6; ++k) + CHECK(result(k) == 0.0); + CHECK(result(6) == Approx(0.0).margin(1e-12)); +} + +TEST_CASE("All minimizers agree on Gaussian fits", "[fit]") { + using aare::model::Gaussian; + const auto none = [](auto &) {}; + check_minimizers_agree("exact, unweighted", gaussian_data(), + false, false, 100, 1e-4, none); + check_minimizers_agree("exact, weighted with errors", + gaussian_data(), true, true, 100, 1e-4, + none); + check_minimizers_agree( + "sigma fixed at a wrong value", gaussian_data(), true, true, 100, 1e-3, + [](auto &m) { m.FixParameter("sigma", 1.5); }, 0.2); + check_minimizers_agree( + "sigma fixed at the true value", gaussian_data(), true, true, 100, 1e-3, + [](auto &m) { m.FixParameter("sigma", 1.3); }); + check_minimizers_agree( + "active upper limit on mu", gaussian_data(), false, false, 100, 1e-3, + [](auto &m) { m.SetParLimits("mu", -5.0, 0.5); }); + check_minimizers_agree( + "inactive limits", gaussian_data(), false, false, 100, 1e-4, + [](auto &m) { m.SetParLimits("mu", -5.0, 5.0); }); + check_minimizers_agree( + "active upper limit on sigma", gaussian_data(), true, true, 100, 1e-3, + [](auto &m) { m.SetParLimits("sigma", 0.5, 1.0); }, 0.2); + check_minimizers_agree( + "user start value", gaussian_data(), false, false, 100, 1e-4, + [](auto &m) { m.SetParameter("mu", 0.0); }); + check_minimizers_agree( + "fixed amplitude and active limit on mu", gaussian_data(), true, true, + 100, 1e-3, + [](auto &m) { + m.FixParameter("A", 100.0); + m.SetParLimits("mu", 0.6, 0.7); + }, + 0.2); + check_minimizers_agree("noisy, weighted with errors", + gaussian_data(true), true, true, 500, 1e-3, + none); + check_minimizers_agree( + "noisy, fixed mu and active limit on sigma", gaussian_data(true), true, + true, 500, 1e-3, + [](auto &m) { + m.FixParameter("mu", 0.8); + m.SetParLimits("sigma", 0.5, 1.1); + }, + 0.2); +} + +TEST_CASE("All minimizers agree on S-curve fits", "[fit]") { + using aare::model::RisingScurve; + const std::vector truth{10.0, 0.1, 50.0, 3.0, 1000.0, 2.0}; + const auto none = [](auto &) {}; + check_minimizers_agree("noisy, weighted with errors", + scurve_data(truth, true), + true, true, 500, 2e-2, none); + check_minimizers_agree("noisy, fixed slope and limited width", + scurve_data(truth, true), + true, true, 500, 1e-2, [](auto &m) { + m.FixParameter("p1", 0.1); + m.SetParLimits("sigma", 1.0, 4.0); + }); + check_minimizers_agree( + "noisy, active limit on the threshold", + scurve_data(truth, true), false, false, 500, 1e-2, + [](auto &m) { m.SetParLimits("mu", 40.0, 48.0); }); + check_minimizers_agree( + "noisy, two fixed, limited amplitude, start value", + scurve_data(truth, true), true, true, 500, 1e-2, + [](auto &m) { + m.FixParameter("sigma", 3.0); + m.FixParameter("p0", 10.0); + m.SetParLimits("A", 900.0, 1100.0); + m.SetParameter("C", 0.0); + }); +} + +TEST_CASE("All minimizers agree on a noisy data cube", "[fit]") { + constexpr ssize_t rows = 3; + constexpr ssize_t cols = 2; + NDArray x({n_points}); + NDArray y({rows, cols, n_points}); + NDArray y_err({rows, cols, n_points}, 1.0); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + const double mu = -0.6 + 0.3 * static_cast(col); + const double sigma = 0.9 + 0.1 * static_cast(row); + for (ssize_t i = 0; i < n_points; ++i) { + x(i) = -6.0 + 0.2 * static_cast(i); + const double z = (x(i) - mu) / sigma; + y(row, col, i) = 80.0 * std::exp(-0.5 * z * z) + + 2.0 * std::sin(7.3 * static_cast(i)); + } + } + } + + struct Cube { + NDArray par{{rows, cols, 3}}; + NDArray err{{rows, cols, 3}}; + NDArray chi2{{rows, cols}}; + }; + const auto fit_cube = [&](Minimizer minimizer) { + Cube c; + const auto model = + make_model(minimizer, true, 500); + aare::fit_3d(model, x.view(), y.view(), y_err.view(), c.par.view(), + c.err.view(), c.chi2.view(), 2); + return c; + }; + const auto ref = fit_cube(Minimizer::Migrad); + + for (const auto candidate : + {Minimizer::Fumili, Minimizer::LevenbergMarquardt, + Minimizer::VarPro}) { + INFO(name(candidate) << " vs Migrad"); + const auto res = fit_cube(candidate); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + for (ssize_t k = 0; k < 3; ++k) { + CHECK(res.par(row, col, k) == + Approx(ref.par(row, col, k)).epsilon(1e-3)); + CHECK(res.err(row, col, k) == + Approx(ref.err(row, col, k)).epsilon(2e-2)); + } + CHECK(res.chi2(row, col) == + Approx(ref.chi2(row, col)).epsilon(1e-4)); + } + } + } +} + +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]. +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{}; + 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); + 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)); + } + CHECK(Model::eval(x, 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]; + 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 + +TEST_CASE("Separable models are consistent with eval and eval_and_grad", + "[fit]") { + using namespace aare::model; + const std::vector xs = {-4.0, -1.5, 0.3, 2.0, + 7.5, 40.0, 55.0, 120.0}; + check_separable({2.0, -0.5}, xs); + check_separable({2.0, -0.5, 0.05}, xs); + check_separable({120.0, 0.8, 1.3}, xs); + check_separable({80.0, 15.0, 2.0, 1.5}, xs); + check_separable({20.0, 0.05, 50.0, 4.0, 1500.0, 0.3}, + xs); + check_separable( + {20.0, 0.05, 50.0, 4.0, 1500.0, 0.3, 1.1, 0.1}, xs); + check_separable({10.0, 0.1, 50.0, 3.0, 1000.0, 2.0}, xs); + check_separable({10.0, 0.1, 50.0, 3.0, 1000.0, 2.0}, xs); +} + +TEST_CASE("VarPro falls back to LevenbergMarquardt for a limit " + "on a linear parameter", + "[fit]") { + using aare::model::FallingScurve; + const std::vector truth{10.0, 0.1, 50.0, 3.0, 1000.0, 2.0}; + auto d = scurve_data(truth); + + auto vp = make_model(Minimizer::VarPro, true, 500); + auto lm = + make_model(Minimizer::LevenbergMarquardt, true, 500); + // The unconstrained amplitude is 1000, so the linear solve cannot + // honour this limit and the pixel is refitted with the full solver. + vp.SetParLimits("A", 0.0, 500.0); + lm.SetParLimits("A", 0.0, 500.0); + + const auto res_vp = + aare::fit_pixel(vp, d.x.view(), d.y.view(), d.y_err.view()); + const auto res_lm = + aare::fit_pixel(lm, d.x.view(), d.y.view(), d.y_err.view()); + REQUIRE(res_vp.size() == res_lm.size()); + for (ssize_t k = 0; k < res_vp.size(); ++k) + CHECK(res_vp(k) == res_lm(k)); + CHECK(res_vp(4) == Approx(500.0).margin(1e-9)); + CHECK(res_vp(12) > 0.0); +} + +TEST_CASE("fast_exp matches std::exp", "[fit]") { + double max_rel = 0.0; + for (int i = 0; i <= 14000; ++i) { + const double x = -700.0 + 1400.0 * static_cast(i) / 14000.0; + const double reference = std::exp(x); + max_rel = + std::max(max_rel, std::abs(aare::model::fast_exp(x) - reference) / + reference); + } + CHECK(max_rel < 1e-15); + CHECK(aare::model::fast_exp(0.0) == 1.0); + CHECK(aare::model::fast_exp(-1000.0) == aare::model::fast_exp(-700.0)); +} + +TEST_CASE("Fit with a one-sided limit far from the start value", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t n = 20; + NDArray x({n}); + NDArray y({n}); + for (ssize_t i = 0; i < n; ++i) { + x(i) = static_cast(i); + y(i) = 2.0 + 0.5 * x(i); + } + const double inf = std::numeric_limits::infinity(); + const bool upper = GENERATE(true, false); + auto model = make_model(minimizer); + if (upper) + model.SetParLimits(0, -inf, 1e9); + else + model.SetParLimits(0, -1e9, inf); + + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(0) == Approx(2.0).epsilon(1e-6)); + CHECK(res(1) == Approx(0.5).epsilon(1e-6)); +} + +TEST_CASE("A parameter ending exactly on a limit reports an error of 0", + "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t n = 3; + NDArray x({n}); + NDArray y({n}, 2.0); + for (ssize_t i = 0; i < n; ++i) + x(i) = static_cast(i); + // The unconstrained minimum p0 = 2 coincides with the limit. + const bool lower = GENERATE(true, false); + auto model = make_model(minimizer, true); + model.FixParameter(1, 0.0); + if (lower) + model.SetParLimits(0, 2.0, 3.0); + else + model.SetParLimits(0, 1.0, 2.0); + + const auto res = aare::fit_pixel(model, x.view(), y.view()); + REQUIRE(res.size() == 5); + CHECK(res(0) == Approx(2.0).epsilon(1e-6)); + CHECK(res(2) == 0.0); + CHECK(res(3) == 0.0); +} + +// --------------------------------------------------------------------------- +// Regressions of the built-in solvers: fits that used to be reported as valid +// although they had stopped away from the minimum or returned garbage. +// --------------------------------------------------------------------------- + +namespace { + +// Standard normal noise from std::mt19937, whose output the standard fixes, +// and a Box-Muller transform written out here: std::normal_distribution +// differs between standard libraries. +std::vector normal_noise(std::size_t n, unsigned seed) { + std::mt19937 rng(seed); + const auto uniform = [&rng] { + return (static_cast(rng()) + 0.5) / 4294967296.0; + }; + std::vector out(n); + for (std::size_t i = 0; i < n; i += 2) { + const double r = std::sqrt(-2.0 * std::log(uniform())); + const double phi = 2.0 * std::acos(-1.0) * uniform(); + out[i] = r * std::cos(phi); + if (i + 1 < n) + out[i + 1] = r * std::sin(phi); + } + return out; +} + +} // namespace + +TEST_CASE("A Gaussian fit to pure noise fails or stays finite and bounded", + "[fit]") { + // Fitting noise, a narrow width can collapse between two scan points; + // the amplitude that fits a single point then diverges. Such a fit must + // be reported as failed, not as a valid result with |A| up to 1e159. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t rows = 20; + constexpr ssize_t cols = 20; + NDArray x({n_points}); + for (ssize_t i = 0; i < n_points; ++i) + x(i) = -6.0 + 0.2 * static_cast(i); + NDArray y({rows, cols, n_points}); + const auto noise = normal_noise(static_cast(y.size()), 1); + std::copy(noise.begin(), noise.end(), y.begin()); + + NDArray par({rows, cols, 3}); + NDArray chi2({rows, cols}); + const auto model = make_model(minimizer); + aare::fit_3d(model, x.view(), y.view(), NDView{}, par.view(), + NDView{}, chi2.view(), 1); + + int absurd = 0; + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + bool ok = std::isfinite(chi2(row, col)); + for (ssize_t k = 0; k < 3; ++k) + ok = ok && std::isfinite(par(row, col, k)) && + std::abs(par(row, col, k)) < 1e6; + absurd += ok ? 0 : 1; + } + } + CHECK(absurd == 0); +} + +TEST_CASE("A step reflected back onto the start point is not convergence", + "[fit]") { + // p0 starts at 0 with limits [-1, 1]; the damped step of +2 reflects at + // the upper limit back onto 0. The step truncated at the limit must still + // be tried: the constrained minimum is p0 = 1. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + NDArray x({4}); + NDArray y({4}, 2.2); + for (ssize_t i = 0; i < 4; ++i) + x(i) = static_cast(i); + + auto model = make_model(minimizer); + model.SetParameter("p0", 0.0); + model.FixParameter("p1", 0.0); + model.SetParLimits("p0", -1.0, 1.0); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(0) == Approx(1.0).margin(limit_margin(minimizer))); + CHECK(res(2) == Approx(5.76).epsilon(1e-6)); +} + +TEST_CASE("A large parameter does not hide the steps of the others", "[fit]") { + // The negligible-step test used one relative norm over all parameters, + // so an amplitude of 1e10 or a position of 1e9 stopped the fit at the + // start estimate of the width. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const auto [amplitude, offset] = GENERATE( + std::pair{1e8, 0.0}, std::pair{1e10, 0.0}, std::pair{120.0, 1e9}); + INFO("amplitude " << amplitude << ", offset " << offset); + + Data d{NDArray({n_points}), NDArray({n_points}), + NDArray({n_points})}; + fill_gaussian(d.x, d.y, amplitude, 0.8, 1.3); + for (ssize_t i = 0; i < n_points; ++i) { + d.x(i) += offset; + d.y_err(i) = std::sqrt(d.y(i)) + 1.0; + } + const auto model = make_model(minimizer, false, 500); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + // The EDM tolerance leaves a few 1e-4 relative, far below the errors; + // the stuck fits were 0.7 % (A = 1e8) and 5 % (start estimate) off. + CHECK(res(0) == Approx(amplitude).epsilon(1e-3)); + CHECK(res(1) - offset == Approx(0.8).margin(1e-3)); + CHECK(res(2) == Approx(1.3).epsilon(1e-3)); +} + +TEST_CASE("A rejected trial on a limit does not freeze an interior parameter", + "[fit]") { + // mu starts at 3 inside [0, 5]. A trial truncated onto the lower limit is + // rejected; testing that trial's position against the gradient of the + // accepted point froze mu, and the fit stopped at the start value. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t n = 11; + const double mu = std::sqrt(2.0 * std::log(2.0)); + NDArray x({n}); + NDArray y({n}); + for (ssize_t i = 0; i < n; ++i) { + x(i) = -0.1 + 0.02 * static_cast(i); + y(i) = 10.0 * std::exp(-0.5 * (x(i) - mu) * (x(i) - mu)); + } + auto model = make_model(minimizer); + model.FixParameter("A", 10.0); + model.FixParameter("sigma", 1.0); + model.SetParameter("mu", 3.0); + model.SetParLimits("mu", 0.0, 5.0); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(1) == Approx(mu).margin(1e-3)); + CHECK(res(3) < 1e-4); +} + +// --------------------------------------------------------------------------- +// Samples with a zero error are ignored. They may hold NaN (a masked dead +// frame), which a zero weight does not remove from a sum, so they must not +// reach the start estimate or the minimizer at all. +// --------------------------------------------------------------------------- + +TEST_CASE("Samples with a zero error are ignored even when NaN", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const double nan = std::numeric_limits::quiet_NaN(); + + constexpr ssize_t n = 10; + NDArray x({n}); + NDArray y({n}); + NDArray y_err({n}, 1.0); + for (ssize_t i = 0; i < n; ++i) { + x(i) = static_cast(i); + y(i) = 2.0 + 0.5 * x(i); + } + const auto line = make_model(minimizer, true); + + // The start estimate of Pol1 uses the first and the last sample. + const auto masked = GENERATE(as{}, 4, 0, n - 1); + SECTION("masked NaN on a line") { + INFO("masked sample " << masked); + y(masked) = nan; + y_err(masked) = 0.0; + const auto res = + aare::fit_pixel(line, x.view(), y.view(), y_err.view()); + CHECK(res(0) == Approx(2.0).epsilon(1e-5).margin(1e-6)); + CHECK(res(1) == Approx(0.5).epsilon(1e-5).margin(1e-6)); + CHECK(res(2) > 0.0); + CHECK(res(4) == Approx(0.0).margin(1e-6)); + } + + SECTION("masked NaN at the peak of a Gaussian") { + auto d = gaussian_data(); + constexpr ssize_t peak = 34; // x = 0.8 + d.y(peak) = nan; + d.y_err(peak) = 0.0; + const auto model = make_model(minimizer); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + const double tol = truth_tolerance(minimizer); + CHECK(res(0) == Approx(120.0).epsilon(tol)); + CHECK(res(1) == Approx(0.8).epsilon(tol)); + CHECK(res(2) == Approx(1.3).epsilon(tol)); + } + + SECTION("too few samples left fail cleanly") { + // One sample left for two free parameters, then none at all. + for (ssize_t i = 1; i < n; ++i) { + y(i) = nan; + y_err(i) = 0.0; + } + auto res = aare::fit_pixel(line, x.view(), y.view(), y_err.view()); + for (ssize_t k = 0; k < res.size(); ++k) + CHECK(res(k) == 0.0); + y_err(0) = 0.0; + res = aare::fit_pixel(line, x.view(), y.view(), y_err.view()); + for (ssize_t k = 0; k < res.size(); ++k) + CHECK(res(k) == 0.0); + } +} + +TEST_CASE("Masked samples in a data cube match per-pixel fits", "[fit]") { + // The kept samples of each pixel go through reused buffers, which must + // not carry anything over from one pixel to the next. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t rows = 2; + constexpr ssize_t cols = 3; + // masked sample per pixel, -1 for none + constexpr ssize_t masked[rows][cols] = {{-1, 10, 0}, + {34, -1, n_points - 1}}; + NDArray x({n_points}); + NDArray tmp({n_points}); + NDArray y({rows, cols, n_points}); + NDArray y_err({rows, cols, n_points}, 1.0); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + fill_gaussian(x, tmp, 80.0 + 10.0 * static_cast(col), 0.8, + 1.0 + 0.2 * static_cast(row)); + for (ssize_t i = 0; i < n_points; ++i) + y(row, col, i) = tmp(i); + if (masked[row][col] >= 0) { + y(row, col, masked[row][col]) = + std::numeric_limits::quiet_NaN(); + y_err(row, col, masked[row][col]) = 0.0; + } + } + } + + NDArray par({rows, cols, 3}); + NDArray err({rows, cols, 3}); + NDArray chi2({rows, cols}); + const auto model = make_model(minimizer, true); + aare::fit_3d(model, x.view(), y.view(), y_err.view(), par.view(), + err.view(), chi2.view(), 1); + + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + INFO("pixel " << row << ", " << col); + NDView yv(&y(row, col, 0), {n_points}); + NDView ev(&y_err(row, col, 0), {n_points}); + const auto single = aare::fit_pixel(model, x.view(), yv, ev); + for (ssize_t k = 0; k < 3; ++k) { + CHECK(par(row, col, k) == Approx(single(k)).margin(1e-12)); + CHECK(err(row, col, k) == Approx(single(3 + k)).margin(1e-12)); + } + CHECK(chi2(row, col) == Approx(single(6)).margin(1e-12)); + // Close to the truth within the EDM tolerance; the exact check + // is the agreement with the per-pixel fit above. + const double tol = 1e-3; + CHECK(par(row, col, 0) == + Approx(80.0 + 10.0 * static_cast(col)).epsilon(tol)); + CHECK(par(row, col, 1) == Approx(0.8).epsilon(tol)); + CHECK(par(row, col, 2) == + Approx(1.0 + 0.2 * static_cast(row)).epsilon(tol)); + } + } +} + +TEST_CASE("A parameter whose minimum lies on its limit reports no error", + "[fit]") { + // At a minimum that lies exactly on a limit the gradient is zero, so a + // test for an outward gradient does not see the parameter as bound. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + + SECTION("linear parameter") { + NDArray x({4}); + NDArray y({4}); + const double xs[4] = {-1.0, -1.0, 1.0, 1.0}; + for (ssize_t i = 0; i < 4; ++i) { + x(i) = xs[i]; + y(i) = 1.0 + 2.0 * xs[i]; + } + auto model = make_model(minimizer, true); + model.SetParLimits("p0", 0.0, 1.0); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(0) == Approx(1.0).margin(limit_margin(minimizer))); + CHECK(res(2) == 0.0); + CHECK(res(3) == Approx(0.5).epsilon(1e-6)); // 1 / sqrt(sum x^2) + } + + SECTION("nonlinear parameter") { + auto d = gaussian_data(); + auto model = make_model(minimizer, true); + model.SetParLimits("mu", -5.0, 0.8); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(1) == Approx(0.8).margin(1e-6)); + CHECK(res(4) == 0.0); + CHECK(res(3) > 0.0); + CHECK(res(5) > 0.0); + } +} + +TEST_CASE("A fit does not stop at its start point without taking a step", + "[fit]") { + // The EDM tolerance is absolute in chi2 units. An unweighted fit of data + // of small magnitude passes it at the start estimate already; the fit + // must still take a step, as Migrad always iterates at least once. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const double amplitude = GENERATE(0.1, 1e-4); + INFO("amplitude " << amplitude); + Data d{NDArray({n_points}), NDArray({n_points}), + NDArray({n_points}, 1.0)}; + fill_gaussian(d.x, d.y, amplitude, 0.8, 1.3); + + // The start estimate of sigma is 1.3617; Migrad ends at 1.3205. + const auto model = make_model(minimizer); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + CHECK(std::abs(res(2) - 1.3) < 0.03); + + // Without the budget for a single step the fit has not converged. + const auto one_call = + make_model(minimizer, false, 1); + const auto none = aare::fit_pixel(one_call, d.x.view(), d.y.view()); + for (ssize_t k = 0; k < none.size(); ++k) + CHECK(none(k) == 0.0); +} + +TEST_CASE("A polynomial far from x = 0 keeps its precision", "[fit]") { + // On a narrow range far from 0 the columns 1, x, x^2 are nearly + // collinear. VarPro solves them from the normal equations and lost six + // digits at x = 20000; from about x = 60000 it falls back to LM. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const auto [lo, span] = + GENERATE(std::pair{8000.0, 1000.0}, std::pair{20000.0, 500.0}, + std::pair{60000.0, 200.0}); + INFO("x in [" << lo << ", " << lo + span << "]"); + // 1.5 + 2e-3 (x - lo) + 3e-7 (x - lo)^2 expanded in powers of x. + const double c[3] = {1.5 - 2e-3 * lo + 3e-7 * lo * lo, 2e-3 - 6e-7 * lo, + 3e-7}; + constexpr ssize_t n = 40; + NDArray x({n}); + NDArray y({n}); + for (ssize_t i = 0; i < n; ++i) { + x(i) = lo + span * static_cast(i) / static_cast(n - 1); + y(i) = c[0] + c[1] * x(i) + c[2] * x(i) * x(i); + } + const auto model = make_model(minimizer); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + for (ssize_t k = 0; k < 3; ++k) + CHECK(res(k) == Approx(c[k]).epsilon(1e-9)); +} diff --git a/src/fit/FitHelpers.hpp b/src/fit/FitHelpers.hpp new file mode 100644 index 00000000..a82599cf --- /dev/null +++ b/src/fit/FitHelpers.hpp @@ -0,0 +1,83 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include "aare/FitModel.hpp" +#include "aare/NDView.hpp" +#include +#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. + +namespace aare::detail { + +/** + * @brief Starting values for one pixel. + * + * Model::estimate_par supplies data-driven estimates. User start values and + * fixed values take precedence, and free parameters are clamped into their + * limits. + */ +template +std::array start_values(const FitModel &model, + NDView x, + NDView y) { + auto start = Model::estimate_par(x, y); + for (std::size_t k = 0; k < Model::npar; ++k) { + const auto idx = static_cast(k); + if (model.is_user_fixed(idx) || model.is_user_start(idx)) + start[k] = model.value(idx); + if (!model.is_user_fixed(idx)) + start[k] = std::clamp(start[k], model.lower_limit(idx), + model.upper_limit(idx)); + } + return start; +} + +/** @brief Validate the array sizes passed to fit_pixel. */ +inline void check_fit_pixel_shapes(NDView x, NDView y, + NDView y_err) { + if (x.size() != y.size()) + throw std::runtime_error("fit_pixel: x.size() must match y.size()."); + if (y_err.size() > 0 && y_err.size() != y.size()) + throw std::runtime_error( + "fit_pixel: y and y_err must have identical size."); +} + +/** @brief Validate the array shapes passed to fit_3d. */ +template +void check_fit_3d_shapes(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView err_out, + NDView chi2_out) { + constexpr auto npar = static_cast(Model::npar); + + if (x.size() != y.shape(2)) + throw std::runtime_error("fit_3d: x.size() must match y.shape(2)."); + + if (par_out.shape(0) != y.shape(0) || par_out.shape(1) != y.shape(1) || + par_out.shape(2) != npar) + throw std::runtime_error("par_out must have shape [rows, cols, npar]."); + + 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) { + 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) { + 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 new file mode 100644 index 00000000..326c4b11 --- /dev/null +++ b/src/fit/FitMinuit2.cpp @@ -0,0 +1,168 @@ +// SPDX-License-Identifier: MPL-2.0 +// Minuit2 backends of the pixel fit: Migrad (MnHesse for the errors) and +// Fumili. This is the only translation unit that includes Minuit2 headers. +#include "FitMinuit2.hpp" +#include "Chi2.hpp" +#include "Minuit2/FunctionMinimum.h" +#include "Minuit2/MnFumiliMinimize.h" +#include "Minuit2/MnHesse.h" +#include "Minuit2/MnMigrad.h" +#include "Minuit2/MnStrategy.h" +#include "Minuit2/MnUserParameterState.h" +#include "Minuit2/MnUserParameters.h" +#include "MinuitSteps.hpp" +#include "aare/Models.hpp" +#include +#include +#include +#include +#include + +namespace aare::detail { + +namespace { + +// Translate the FitModel state into Minuit2 user parameters for one pixel. +template +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); + + ROOT::Minuit2::MnUserParameters upar; + for (std::size_t i = 0; i < 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); + const double hi = model.upper_limit(idx); + if (std::isfinite(lo) && std::isfinite(hi)) { + upar.SetLimits(idx, lo, hi); + } else if (std::isfinite(lo) || std::isfinite(hi)) { + // Minuit2's transformation for one-sided limits converges markedly + // slower than the two-sided one, so an open side gets a wide bound. + // It is placed relative to the start value, which then always + // lies inside the range however far it is from the finite limit. + const double span = std::max(1e6, 1e3 * std::abs(start[i])); + if (std::isfinite(lo)) + upar.SetLimits(idx, lo, std::max(start[i], lo) + span); + else + upar.SetLimits(idx, std::min(start[i], hi) - span, hi); + } + if (model.is_user_fixed(idx)) + upar.Fix(idx); + } + return upar; +} + +// Minuit2 approaches an active limit asymptotically through its parameter +// transformation, so "on a limit" is decided with a small relative margin. +template +bool on_limit(const FitModel &model, unsigned int idx, double value) { + const double margin = 1e-6 * std::max(1.0, std::abs(value)); + return std::abs(value - model.lower_limit(idx)) <= margin || + std::abs(value - model.upper_limit(idx)) <= margin; +} + +// Copy the minimum into the flat result. Fixed parameters and parameters on +// a limit report an error of 0. +template +void unpack_minimum(const FitModel &model, + const ROOT::Minuit2::FunctionMinimum &min, double *par_out, + double *err_out, double &chi2) { + constexpr std::size_t npar = Model::npar; + const auto &state = min.UserState(); + const std::vector values = state.Params(); + for (std::size_t k = 0; k < npar; ++k) + par_out[k] = values[k]; + + if (err_out) { + const std::vector errors = state.Errors(); + for (std::size_t k = 0; k < npar; ++k) { + const auto idx = static_cast(k); + const bool no_error = + model.is_user_fixed(idx) || on_limit(model, idx, values[k]); + err_out[k] = no_error ? 0.0 : errors[k]; + } + } + chi2 = min.Fval(); +} + +} // namespace + +template +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; + + const auto upar = user_parameters(model, x, y, start); + const ROOT::Minuit2::MnStrategy strategy(model.strategy()); + + if (model.minimizer() == Minimizer::Fumili) { + // Fumili takes the gradient and the Gauss-Newton Hessian from the FCN, + // so the covariance of the minimum is available without MnHesse. + func::Chi2Model1DFumili fcn(x, y, y_err); + ROOT::Minuit2::MnFumiliMinimize fumili( + fcn, ROOT::Minuit2::MnUserParameterState(upar), strategy); + const ROOT::Minuit2::FunctionMinimum min = + fumili(model.max_calls(), model.tolerance()); + if (min.IsValid()) { + unpack_minimum(model, min, par_out, err_out, chi2); + return true; + } + // Minuit2's Fumili cannot invert its linearised Hessian once a + // two-sided limit becomes active. Refit such pixels with Migrad. + } + + func::Chi2Model1DGrad fcn(x, y, y_err); + ROOT::Minuit2::MnMigrad migrad(fcn, upar, strategy); + ROOT::Minuit2::FunctionMinimum min = + migrad(model.max_calls(), model.tolerance()); + if (!min.IsValid()) + return false; + + // MnHesse cannot handle a fit in which every parameter is fixed. + if (err_out && upar.VariableParameters() > 0) { + ROOT::Minuit2::MnHesse hesse; + hesse(fcn, min); + } + unpack_minimum(model, min, par_out, err_out, chi2); + return true; +} + +// NOLINTBEGIN +#define AARE_INSTANTIATE_MINUIT2_FIT(Model) \ + 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) + +#undef AARE_INSTANTIATE_MINUIT2_FIT +// NOLINTEND + +} // namespace aare::detail diff --git a/src/fit/FitMinuit2.hpp b/src/fit/FitMinuit2.hpp new file mode 100644 index 00000000..c0240be8 --- /dev/null +++ b/src/fit/FitMinuit2.hpp @@ -0,0 +1,30 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include "aare/FitModel.hpp" +#include "aare/NDView.hpp" +#include + +// Minuit2 backends (Migrad and Fumili) of the pixel fit. FitMinuit2.cpp is +// the only translation unit that includes Minuit2 headers; it defines and +// explicitly instantiates the function below for every supported model. + +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. + * @return true when Minuit2 reports a valid minimum. + */ +template +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); + +} // namespace aare::detail diff --git a/src/fit/LevenbergMarquardt.hpp b/src/fit/LevenbergMarquardt.hpp new file mode 100644 index 00000000..8a6a7213 --- /dev/null +++ b/src/fit/LevenbergMarquardt.hpp @@ -0,0 +1,173 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include "DampedGaussNewton.hpp" +#include "aare/FitModel.hpp" +#include "aare/NDView.hpp" + +#include +#include +#include +#include +#include + +namespace aare::detail { + +/** + * @brief Levenberg-Marquardt least-squares minimiser for one FitModel. + * + * Minimises chi2 = sum_i ((y_i - f(x_i; p)) / s_i)^2 using the analytic + * gradient of the model (s_i = 1 for unweighted fits; points with s_i == 0 + * are ignored). Keep one instance per thread: the buffers are reused between + * pixels so that after the first fit no memory is allocated. + * + * Every evaluation computes the residuals r = (y - f) / s and the Jacobian + * J = dr/dp in one pass over the data and reduces them to the normal + * 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. + * + * 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 + * report 0. + */ +template class LevenbergMarquardt { + public: + static constexpr int npar = static_cast(Model::npar); + 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; + + /** + * @brief Fit one pixel. + * + * @param model Limits, fixed flags and minimiser settings. + * @param x Scan points. + * @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. + */ + Result 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()); + Vec lower{}; + Vec upper{}; + std::array fixed{}; + for (int k = 0; k < npar; ++k) { + const auto idx = static_cast(k); + lower[k] = model.lower_limit(idx); + 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; + + 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); + return res; + } + std::copy(driver_.point().begin(), driver_.point().end(), par_out); + if (err_out) + driver_.errors(err_out); + return res; + } + + private: + double weight(ssize_t i) const { + return weighted_ ? w_[static_cast(i)] : 1.0; + } + + // 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. + // 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. + bool evaluate(const Vec &p, Normal &out) { + pvec_.assign(p.begin(), p.end()); + if (!Model::is_valid(pvec_)) + return false; + const ssize_t n = x_.size(); + double F = 0.0; + std::array g{}; + double f = 0.0; + for (ssize_t i = 0; i < n; ++i) { + const double w = weight(i); + Model::eval_and_grad(x_[i], pvec_, f, g); + r_[i] = (y_[i] - f) * w; + F += 0.5 * r_[i] * r_[i]; + double *Ji = &J_[static_cast(i) * npar]; + for (int k = 0; k < npar; ++k) + Ji[k] = -g[k] * w; + } + Vec G{}; + Mat 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) { + G[a] += Ji[a] * ri; + for (int b = a; b < npar; ++b) + H[a][b] += Ji[a] * Ji[b]; + } + } + 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]; + } + return true; + } + + NDView x_, y_, s_; + bool weighted_ = false; + Driver driver_; + std::vector r_; + std::vector J_; + std::vector w_; // 1 / s_i for weighted fits + std::vector pvec_; +}; + +} // namespace aare::detail diff --git a/src/fit/MinuitSteps.hpp b/src/fit/MinuitSteps.hpp new file mode 100644 index 00000000..d5aa8bbe --- /dev/null +++ b/src/fit/MinuitSteps.hpp @@ -0,0 +1,126 @@ +// 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. + +#include "aare/Models.hpp" +#include "aare/NDView.hpp" +#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) { + 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); +} + +/** @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; + } +}; + +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); + } +}; + +} // namespace aare::minuit2 diff --git a/src/fit/VariableProjection.hpp b/src/fit/VariableProjection.hpp new file mode 100644 index 00000000..0ea0cdbc --- /dev/null +++ b/src/fit/VariableProjection.hpp @@ -0,0 +1,511 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include "DampedGaussNewton.hpp" +#include "aare/FitModel.hpp" +#include "aare/Models.hpp" +#include "aare/NDView.hpp" + +#include +#include +#include +#include +#include + +namespace aare::detail { + +/** + * @brief Variable projection least-squares minimiser for a separable + * FitModel (Golub and Pereyra, 1973). + * + * The model declares the parameters it is linear in, see + * model::is_separable. For a given point of the nonlinear parameters the + * linear ones are solved exactly from their normal equations, and the + * damped Gauss-Newton iteration of DampedGaussNewton runs over the + * nonlinear parameters only, on the projected residual + * r(alpha) = y - Phi(alpha) beta(alpha). Its Jacobian uses Kaufman's + * approximation, the derivatives of the model at fixed beta projected out + * of the span of the basis functions, and the reduced Gauss-Newton Hessian + * 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. + * + * 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. + * + * Keep one instance per thread: the buffers are reused between pixels so + * that after the first fit no memory is allocated. + * + * - Fixed nonlinear parameters are left out of the iteration. A fixed + * linear parameter keeps its value and moves to the data side of the + * linear solve. + * - Limits on nonlinear parameters are handled by the driver. A limit on a + * linear parameter cannot be enforced by the linear solve: when the + * 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. + * - 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. + */ +template class VariableProjection { + public: + static_assert(model::is_separable::value, + "VariableProjection needs a model with linear_par and " + "basis_and_grad"); + using Traits = model::separable_traits; + static constexpr int npar = static_cast(Model::npar); + static constexpr int nlin = static_cast(Traits::nlin); + static constexpr int nnl = static_cast(Traits::nnl); + static_assert(Model::npar >= 1 && Model::npar <= 16, + "VariableProjection keeps npar x npar matrices on the stack"); + 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. + static constexpr int ndrv = nnl > 0 ? nnl : 1; + static constexpr int ncol = nnl * nlin > 0 ? nnl *nlin : 1; + 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) { + x_ = x; + y_ = y; + const auto n = static_cast(x.size()); + phi_.resize(n * nlin); + dphi_.resize(n * nlin * nnl); + weighted_ = s.size() > 0; + if (weighted_) { + // 1 / s_i^2 once per pixel; points with s_i == 0 are ignored. + w2_.resize(n); + for (std::size_t i = 0; i < n; ++i) { + const double si = s[static_cast(i)]; + w2_[i] = si != 0.0 ? 1.0 / (si * si) : 0.0; + } + } + pvec_.assign(start.begin(), start.end()); + have_last_ = false; + have_ref_ = false; + nfree_lin_ = 0; + for (int j = 0; j < nlin; ++j) { + const auto k = static_cast(Traits::linear_index(j)); + lin_fixed_[j] = model.is_user_fixed(k); + if (!lin_fixed_[j]) + free_lin_[nfree_lin_++] = j; + } + Vec alpha{}; + Vec lower{}; + Vec upper{}; + std::array fixed{}; + for (int k = 0; k < nnl; ++k) { + const auto idx = + static_cast(Traits::nonlinear_par[k]); + alpha[k] = start[idx]; + lower[k] = model.lower_limit(idx); + 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; + + 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; + }; + + 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. + for (int j = 0; j < nlin; ++j) { + if (lin_fixed_[j]) + continue; + const auto k = static_cast(Traits::linear_index(j)); + const double lo = model.lower_limit(k); + 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(); + 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; + } + + private: + // F = chi2 / 2, gradient and Kaufman Hessian of the projected residual + // at the nonlinear point alpha, and the linear 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(). + bool evaluate(const Vec &alpha, Normal &out) { + for (int k = 0; k < nnl; ++k) + 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); + } + + // 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) { + const ssize_t n = x_.size(); + const auto un = static_cast(n); + std::array pc; + for (int j = 0; j < nlin; ++j) + pc[j] = phi_.data() + static_cast(j) * un; + std::array dc; + for (int c = 0; c < nnl * nlin; ++c) + dc[c] = dphi_.data() + static_cast(c) * un; + const double *yd = y_.data(); + const double *w2 = Weighted ? w2_.data() : nullptr; + + // Pass 1: normal equations M = Phi^T W Phi, t = Phi^T W y of all + // linear parameters. + Mat M{}; + Vec t{}; + for (ssize_t i = 0; i < n; ++i) { + const double wi = Weighted ? w2[i] : 1.0; + const double yi = yd[i]; + Vec f; + for (int j = 0; j < nlin; ++j) + f[j] = pc[j][i]; + for (int a = 0; a < nlin; ++a) { + const double fa = f[a] * wi; + t[a] += fa * yi; + for (int b = a; b < nlin; ++b) + M[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{}; + Vec v{}; + for (int a = 0; a < nlin; ++a) { + if (a >= nfree_lin_) { + v[a] = 0.0; + continue; + } + const int ja = free_lin_[a]; + double va = t[ja]; + for (int j = 0; j < nlin; ++j) + if (lin_fixed_[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]]; + } + // A free basis function that has all but vanished on the scan points + // (a Gaussian narrower than the point spacing, centred between two + // points) leaves its coefficient undetermined, and the solve returns + // an absurd one that fits a single point. Its squared norm is compared + // 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]; + 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])) + return false; + } + Mat L{}; + if (!cholesky_factor(nfree_lin_, A, L)) + 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]; + } + Vec beta_free{}; + cholesky_substitute(nfree_lin_, L, 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 + // recovered by one step of iterative refinement on the data (the + // corrected semi-normal equations): the correction solves the same + // equations for Phi^T W r of the residual of the first solution. A + // polynomial fitted far from x = 0 needs it; x in [20000, 20500] + // otherwise costs a Pol2 six digits. + if (refine) { + Vec beta; + for (int j = 0; j < nlin; ++j) + beta[j] = pvec_[Traits::linear_index(j)]; + Vec s{}; + for (ssize_t i = 0; i < n; ++i) { + double r = yd[i]; + for (int j = 0; j < nlin; ++j) + r -= beta[j] * pc[j][i]; + const double rw = Weighted ? r * w2[i] : r; + for (int j = 0; j < nlin; ++j) + 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; + Vec delta{}; + cholesky_substitute(nfree_lin_, L, 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) { + 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{}; + for (ssize_t i = 0; i < n; ++i) { + const double wi = Weighted ? w2[i] : 1.0; + Vec f; + for (int j = 0; j < nlin; ++j) + f[j] = pc[j][i]; + double r = yd[i]; + for (int j = 0; j < nlin; ++j) + r -= beta[j] * f[j]; + Vec u{}; + for (int k = 0; k < nnl; ++k) + for (int j = 0; j < nlin; ++j) + u[k] += beta[j] * dc[k * nlin + j][i]; + const double rw = r * wi; + F += 0.5 * r * rw; + for (int k = 0; k < nnl; ++k) { + const double uw = u[k] * wi; + b[k] += uw * r; + for (int l = k; l < nnl; ++l) + N[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; + + // 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{}; + 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]); + } + 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]; + for (int a = 0; a < nfree_lin_; ++a) + h -= qf[k][a] * z[l][a]; + out.H[k][l] = h; + } + } + return true; + } + + // 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; + } + } + // Fixed parameters and parameters that end on a limit have no error. + // The driver knows the nonlinear ones; the linear solve can land a + // linear parameter exactly on a limit without violating it. + std::array skip; + for (int k = 0; k < npar; ++k) { + const auto idx = static_cast(k); + const double v = pvec_[static_cast(k)]; + 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); + } + 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_{}; + static constexpr double vanished_basis = 1e-8; + static constexpr double collinear_basis = 1e-10; + // Pivot ratio below which the linear solution is refined, about four + // lost digits. + static constexpr double refine_basis = 1e-4; + int nfree_lin_ = 0; + std::array free_lin_{}; + std::array lin_fixed_{}; + std::vector phi_; // nlin columns of n basis values + std::vector dphi_; // nnl * nlin columns of their derivatives + std::vector w2_; // 1 / s_i^2 for weighted fits + std::vector pvec_; +}; + +} // namespace aare::detail