diff --git a/AGENTS.md b/AGENTS.md index c4918858..083b7c0f 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -145,10 +145,20 @@ 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`. `FitModel` is plain data (limits, fixed flags, start values, +minimizer settings) and selects one of three minimizers through +`aare::Minimizer`: Minuit2's Migrad (default) and Fumili, or the built-in +Levenberg-Marquardt solver in `src/LevenbergMarquardt.hpp`. `src/Fit.cpp` +holds the `FitModel` method bodies, the per-thread `detail::PixelFitter` +that dispatches on the minimizer, `fit_pixel`/`fit_3d`, and the explicit +instantiations; `src/FitHelpers.hpp` holds the start-value precedence and +shape checks. `src/FitMinuit2.cpp` (with `src/Chi2.hpp` and +`src/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 three +minimizers and check that they agree; run them when changing any solver. ### Python layer diff --git a/CMakeLists.txt b/CMakeLists.txt index e75ea5d4..4f088391 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -401,6 +401,7 @@ set(SourceFiles ${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/FitMinuit2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp diff --git a/RELEASE.md b/RELEASE.md index 8e20c37d..0f8a0b87 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -29,16 +29,33 @@ - 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 Minuit2 minimizer for the fit models. ``Minimizer.Migrad`` - remains the default; ``Minimizer.Fumili`` uses the analytic model derivatives - for a Gauss-Newton style minimisation that needs far fewer function - evaluations. Select it with the ``minimizer`` constructor argument or property - in Python, or ``FitModel::SetMinimizer`` in C++. With ``compute_errors``, - Fumili reports parameter errors from its linearised covariance instead of - running Hesse. +- 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. 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. ### API Changes: +- ``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, diff --git a/benchmarks/fit_benchmark.cpp b/benchmarks/fit_benchmark.cpp index 6b537f75..df3ea412 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 @@ -127,6 +128,16 @@ 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); +} + // ---------------------------------------------------------------- // Rising S-curve (6 parameters), typical for threshold scans // ---------------------------------------------------------------- @@ -170,8 +181,88 @@ 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); +} + +// ---------------------------------------------------------------- +// 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); +} + 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_FitGausMigradHesse) ->DenseRange(0, 5) @@ -179,8 +270,18 @@ BENCHMARK(BM_FitGausMigradHesse) BENCHMARK(BM_FitGausFumiliErrors) ->DenseRange(0, 5) ->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausLMErrors)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); BENCHMARK(BM_FitScurveMigrad)->Unit(benchmark::kMicrosecond); BENCHMARK(BM_FitScurveFumili)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveLM)->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_MAIN(); diff --git a/docs/src/pyFit.rst b/docs/src/pyFit.rst index 17888299..36c59033 100644 --- a/docs/src/pyFit.rst +++ b/docs/src/pyFit.rst @@ -4,7 +4,8 @@ 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 or a built-in Levenberg-Marquardt solver. Create a model object and +call its :meth:`fit` method:: model = Gaussian(compute_errors=True) result = model.fit(x, y, y_err) @@ -30,24 +31,58 @@ The result dictionary contains ``par`` and ``chi2``. It also contains Choosing the minimizer ~~~~~~~~~~~~~~~~~~~~~~ -Fits use Minuit2's Migrad by default. Fumili is a Gauss-Newton style -minimizer that takes the gradient and a linearised Hessian from the analytic -derivatives of the model, so it usually converges in far fewer function -evaluations. Select it with the ``minimizer`` argument or property:: +Fits use Minuit2's Migrad by default. Two alternatives share the same +interface and are selected with the ``minimizer`` argument or property:: - model = Gaussian(minimizer=Minimizer.Fumili) - model.minimizer = Minimizer.Migrad + model = Gaussian(minimizer=Minimizer.LevenbergMarquardt) + model.minimizer = Minimizer.Fumili -In the bundled ``fit_benchmark`` Fumili is about 1.5 times faster than Migrad -for Gaussian peaks with little noise and about three times faster for the -six-parameter S-curves. For very noisy data the Gauss-Newton steps converge -more slowly, so compare both minimizers on representative data. +- ``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 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. + 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. +- ``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 J)^-1))``. -With ``compute_errors`` enabled, Migrad reports parameter errors from Hesse -while Fumili reports them from the covariance of its linearised Hessian, so -Fumili does not run an extra Hesse step. +All three 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 and a data cube. For very noisy data the Gauss-Newton +steps of Fumili and LevenbergMarquardt converge more slowly, so compare the +minimizers on representative data. + +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`` caps the work per pixel: Minuit2 function calls for Migrad and + Fumili, model evaluations for LevenbergMarquardt. ``0`` selects Minuit's + default budget. A pixel that does not converge within the budget returns + zeros. +- ``tolerance`` is the EDM tolerance in Minuit's convention. Migrad and + LevenbergMarquardt stop when the estimated distance to the minimum is below + ``0.002 * tolerance``, Fumili below ``1e-4 * tolerance``. +- ``strategy`` is the Minuit2 strategy (0 is fast, 1 is Minuit2's default). + LevenbergMarquardt ignores 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. + +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 + :members: Migrad, Fumili, LevenbergMarquardt diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index 34402c6f..cc7edce1 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -7,21 +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 or Fumili) is selected by the model, see - * aare::Minimizer. Parameter errors, when requested, come from MnHesse with - * Migrad and from the linearised covariance with Fumili. + * 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 4670a7cc..ca77985e 100644 --- a/include/aare/FitModel.hpp +++ b/include/aare/FitModel.hpp @@ -3,7 +3,6 @@ #include "aare/Models.hpp" #include -#include #include #include #include @@ -11,68 +10,79 @@ namespace aare { /** - * @brief Minuit2 minimizer used by FitModel. + * @brief Minimizer used by FitModel. * - * - Migrad: variable-metric minimizer driven by the analytic gradient of the - * chi-squared. Parameter errors, when requested, come from MnHesse. - * - Fumili: Gauss-Newton style minimizer for least-squares problems. It takes + * - 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, so it typically needs far fewer function evaluations than - * Migrad. Parameter errors, when requested, come from the covariance of that - * linearised Hessian without an extra Hesse step. + * 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. + * + * 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 and LevenbergMarquardt stop + * below 0.002 * tolerance, Fumili below 1e-4 * tolerance). `max_calls` + * bounds the work per pixel in the minimizer's own units (Minuit2 function + * calls, or model evaluations for LevenbergMarquardt). `strategy` only + * affects the Minuit2 minimizers. */ -enum class Minimizer { Migrad, Fumili }; +enum class Minimizer { Migrad, Fumili, LevenbergMarquardt }; +/** + * @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_; - Minimizer minimizer_; - - 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 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 Fumili covariance with Fumili. - * @param minimizer Minuit2 minimizer, see aare::Minimizer. + * 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, 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); @@ -95,16 +105,34 @@ template class FitModel { 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..1c673e99 100644 --- a/include/aare/Models.hpp +++ b/include/aare/Models.hpp @@ -39,7 +39,7 @@ inline double fast_erf(double x) { /** * @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 { @@ -50,23 +50,6 @@ struct ParamInfo { inline constexpr double no_bound = std::numeric_limits::infinity(); -/** - * @brief Compute data-range statistics used by 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); -} - // _____________________________________________________________________ // // Pol1 @@ -120,14 +103,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; - } }; // _____________________________________________________________________ @@ -220,15 +195,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); - } }; // _____________________________________________________________________ @@ -326,18 +292,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; - } }; // _____________________________________________________________________ @@ -490,16 +444,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; - } }; // _____________________________________________________________________ @@ -647,18 +591,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); - } }; // _____________________________________________________________________ @@ -828,20 +760,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); - } }; // _____________________________________________________________________ @@ -995,18 +913,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; - } }; // _____________________________________________________________________ @@ -1148,14 +1054,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/python/src/fit.hpp b/python/src/fit.hpp index 0a1bfbaf..5c28af40 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -26,7 +26,31 @@ template void bind_fit_model(py::module &m, const char *name) { aare::Minimizer>(), py::arg("strategy") = 0, py::arg("max_calls") = 100, py::arg("tolerance") = 0.5, py::arg("compute_errors") = false, - py::arg("minimizer") = aare::Minimizer::Migrad) + 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``. + max_calls : int + Work budget per pixel: Minuit2 function calls, or model + evaluations for ``Minimizer.LevenbergMarquardt``. 0 selects + Minuit's default budget. + tolerance : float + EDM tolerance in Minuit's convention. Migrad and + LevenbergMarquardt 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 and LevenbergMarquardt. + Fixed parameters and parameters ending on a limit report 0. + minimizer : Minimizer + ``Minimizer.Migrad`` (default), ``Minimizer.Fumili`` or + ``Minimizer.LevenbergMarquardt``. + )doc") .def("SetParLimits", py::overload_cast(&FM::SetParLimits), py::arg("idx"), py::arg("lo"), py::arg("hi")) @@ -92,7 +116,8 @@ template void bind_fit_model(py::module &m, const char *name) { Fit this model to 1D or 3D data using Minuit2. The minimizer is selected by the ``minimizer`` property - (``Minimizer.Migrad`` by default, or ``Minimizer.Fumili``). + (``Minimizer.Migrad`` by default, ``Minimizer.Fumili`` or + ``Minimizer.LevenbergMarquardt``). Parameters ---------- @@ -232,22 +257,67 @@ 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( - Minuit2 minimizer used by the fit models. + Minimizer used by the fit models. - ``Migrad`` is the default variable-metric minimizer. With + ``Migrad`` (default) is Minuit2's variable-metric minimizer. With ``compute_errors`` its parameter errors come from Hesse. - ``Fumili`` is a Gauss-Newton style minimizer that takes the gradient - and a linearised Hessian from the analytic derivatives of the model. - It usually needs far fewer function evaluations than Migrad. With - ``compute_errors`` its parameter errors come from the covariance of - the linearised Hessian, without an extra Hesse step. + ``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. + + 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("Fumili", aare::Minimizer::Fumili) + .value("LevenbergMarquardt", aare::Minimizer::LevenbergMarquardt); // ── Bind model classes ────────────────────────────────────────── bind_fit_model(m, "Gaussian"); @@ -267,88 +337,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. - The model's ``minimizer`` property selects Migrad (default) - or Fumili. + 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) @@ -360,15 +361,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 e1b77cdb..3a89c654 100644 --- a/python/tests/test_Fit.py +++ b/python/tests/test_Fit.py @@ -1,13 +1,41 @@ # SPDX-License-Identifier: MPL-2.0 import numpy as np +import pytest + import aare +MINIMIZERS = [ + aare.Minimizer.Migrad, + aare.Minimizer.Fumili, + aare.Minimizer.LevenbergMarquardt, +] -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 != aare.Minimizer.LevenbergMarquardt + + +@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) @@ -15,39 +43,185 @@ def test_gaussian_model_evaluates_and_fits_data(): np.testing.assert_allclose(result["par"], expected, atol=2e-3) -def test_minimizer_can_be_selected_and_fumili_fits_data(): - x = np.linspace(-5.0, 5.0, 51) - expected = np.array([20.0, 0.5, 1.2]) - +def test_minimizer_can_be_selected(): assert aare.Gaussian().minimizer == aare.Minimizer.Migrad + assert set(aare.Minimizer.__members__) == {"Migrad", "Fumili", "LevenbergMarquardt"} - model = aare.Gaussian(minimizer=aare.Minimizer.Fumili) - assert model.minimizer == aare.Minimizer.Fumili + for minimizer in MINIMIZERS: + assert aare.Gaussian(minimizer=minimizer).minimizer == minimizer - y = model(x, expected) - result = model.fit(x, y) - np.testing.assert_allclose(result["par"], expected, atol=2e-3) - - model.minimizer = aare.Minimizer.Migrad - assert model.minimizer == aare.Minimizer.Migrad + model = aare.Gaussian() + model.minimizer = aare.Minimizer.LevenbergMarquardt + assert model.minimizer == aare.Minimizer.LevenbergMarquardt -def test_fumili_and_migrad_agree_on_3d_data_with_errors(): +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) + + +@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]) +def test_fit_3d_matches_fit_1d(with_errors, 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 with_errors else None + + model = aare.Gaussian(compute_errors=with_errors, minimizer=minimizer) + cube = model.fit(x, y, y_err, n_threads=3) + + assert cube["par"].shape == (rows, cols, 3) + assert cube["chi2"].shape == (rows, cols) + 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) - migrad = aare.Gaussian(compute_errors=True) - fumili = aare.Gaussian(compute_errors=True, minimizer=aare.Minimizer.Fumili) - - with_migrad = migrad.fit(x, y, y_err, n_threads=2) - with_fumili = fumili.fit(x, y, y_err, n_threads=2) - - assert with_fumili["par"].shape == (2, 3, 3) - assert with_fumili["par_err"].shape == (2, 3, 3) - for result in (with_migrad, with_fumili): + 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(with_fumili["par_err"], with_migrad["par_err"], rtol=2e-2) + np.testing.assert_allclose(result["par_err"], reference["par_err"], rtol=2e-2) diff --git a/src/Fit.cpp b/src/Fit.cpp index 3e4a15b7..855f3486 100644 --- a/src/Fit.cpp +++ b/src/Fit.cpp @@ -1,88 +1,45 @@ // SPDX-License-Identifier: MPL-2.0 #include "aare/Fit.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 "FitHelpers.hpp" +#include "FitMinuit2.hpp" +#include "LevenbergMarquardt.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 -// (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, Minimizer minimizer) - : impl_(std::make_unique(strategy)), max_calls_(max_calls), - tolerance_(tolerance), compute_errors_(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) { - 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); - } + lower_[i] = Model::param_info[i].default_lo; + upper_[i] = Model::param_info[i].default_hi; } } -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_), minimizer_(other.minimizer_), - 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_; - minimizer_ = other.minimizer_; - user_fixed_ = other.user_fixed_; - user_start_ = other.user_start_; - } - return *this; +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 (impl_->upar.Name(i) == name) + if (name == Model::param_info[i].name) return static_cast(i); } throw std::runtime_error("FitModel: unknown parameter name '" + name + "'"); @@ -90,19 +47,24 @@ unsigned int FitModel::checked_index(const std::string &name) const { template void FitModel::SetParLimits(unsigned int idx, double lo, double hi) { - impl_->upar.SetLimits(idx, lo, 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); - impl_->upar.Fix(idx); user_fixed_[idx] = true; } template void FitModel::ReleaseParameter(unsigned int idx) { - impl_->upar.Release(idx); + check_index(idx); user_fixed_[idx] = false; } @@ -113,7 +75,8 @@ void FitModel::ReleaseParameter(const std::string &name) { template void FitModel::SetParameter(unsigned int idx, double val) { - impl_->upar.SetValue(idx, val); + check_index(idx); + value_[idx] = val; user_start_[idx] = true; } @@ -135,111 +98,81 @@ void FitModel::SetParLimits(const std::string &name, double lo, template std::string FitModel::GetParName(unsigned int idx) const { - return impl_->upar.GetName(idx); + 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.push_back(GetParName(i)); + names.emplace_back(Model::param_info[i].name); return names; } // ============================================================================ -// fit_pixel / fit_3d — Minuit2 template implementations +// fit_pixel / fit_3d // ============================================================================ -namespace { +namespace detail { /** - * @brief Copy the fitted values, optionally the errors, and the objective - * value of a minimum into the flat result layout documented in Fit.hpp. + * @brief Fits pixels with the minimizer selected by the model. + * + * Keep one instance per thread: the Levenberg-Marquardt buffers are reused + * between pixels so that a data cube is fitted without per-pixel + * allocations. The Minuit2 minimizers create their state per pixel. */ -template -NDArray pack_minimum(const ROOT::Minuit2::FunctionMinimum &min, - bool want_errors) { - const ssize_t result_size = want_errors ? (2 * npar + 1) : (npar + 1); - NDArray result({result_size}); +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) { + const auto start = start_values(model, x, y); - const auto &state = min.UserState(); - const std::vector values = state.Params(); - for (std::size_t k = 0; k < npar; ++k) - result[k] = values[k]; - - if (want_errors) { - const std::vector errors = state.Errors(); - for (std::size_t k = 0; k < npar; ++k) - result[npar + k] = errors[k]; + if (model.minimizer() == Minimizer::LevenbergMarquardt) { + const auto res = + lm_.fit(model, x, y, y_err, start, par_out, err_out); + chi2 = res.valid ? res.chi2 : 0.0; + return res.valid; + } + return fit_pixel_minuit2(model, x, y, y_err, start, par_out, err_out, + chi2); } - result[result_size - 1] = min.Fval(); - return result; -} + private: + LevenbergMarquardt lm_; +}; -} // namespace +} // namespace detail template NDArray fit_pixel(const FitModel &model, NDView x, NDView y, NDView y_err) { 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); + const auto result_size = + static_cast(want_errors ? (2 * npar + 1) : (npar + 1)); + NDArray result({result_size}, 0.0); - auto start = Model::estimate_par(x, y); + 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 - 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]); - } - - 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 chi2(x, y, y_err); - ROOT::Minuit2::MnFumiliMinimize fumili( - chi2, ROOT::Minuit2::MnUserParameterState(upar_local), - model.impl()->strategy); - const ROOT::Minuit2::FunctionMinimum min = - fumili(model.max_calls(), model.tolerance()); - - if (!min.IsValid()) - return NDArray({result_size}, 0.0); - - return pack_minimum(min, want_errors); - } - - func::Chi2Model1DGrad chi2(x, y, y_err); - 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); - - // MnHesse cannot handle a fit in which every parameter is fixed. - if (want_errors && upar_local.VariableParameters() > 0) { - ROOT::Minuit2::MnHesse hesse; - hesse(chi2, min); - } - - return pack_minimum(min, want_errors); + result(result_size - 1) = chi2; + return result; } template @@ -253,55 +186,35 @@ 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]."); + 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(); - 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) { + 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{}; - auto res = fit_pixel(model, x, values, errors); + 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, k) = res(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, k) = res(npar + k); - chi2_out(row, col) = res(2 * npar); - } else { - chi2_out(row, col) = res(npar); + err_out(row, col, static_cast(k)) = err[k]; } + chi2_out(row, col) = chi2; } } }; diff --git a/src/Fit.test.cpp b/src/Fit.test.cpp index 97243e14..fc6b7226 100644 --- a/src/Fit.test.cpp +++ b/src/Fit.test.cpp @@ -3,16 +3,32 @@ #include "aare/FitModel.hpp" #include "aare/Models.hpp" +#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; -void fill_gaussian(aare::NDArray &x, aare::NDArray &y, +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); @@ -21,231 +37,708 @@ void fill_gaussian(aare::NDArray &x, aare::NDArray &y, } } +// 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"; + } + return "unknown"; +} + +bool uses_minuit2(Minimizer m) { return m != Minimizer::LevenbergMarquardt; } + +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}) { + 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("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); +TEST_CASE("Fit unweighted and weighted noise-free Gaussian data", "[fit]") { + const auto minimizer = GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); - const aare::FitModel unweighted_model; + const auto unweighted_model = make_model(minimizer); const auto unweighted = - aare::fit_pixel(unweighted_model, x.view(), y.view()); + aare::fit_pixel(unweighted_model, d.x.view(), d.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)); + 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 aare::FitModel weighted_model(0, 100, 0.5, - true); + const auto weighted_model = + make_model(minimizer, true); const auto weighted = - aare::fit_pixel(weighted_model, x.view(), y.view(), y_err.view()); + aare::fit_pixel(weighted_model, d.x.view(), d.y.view(), d.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)); + 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("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)); - } - } -} - -TEST_CASE("Fumili 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); - - aare::FitModel unweighted_model; - REQUIRE(unweighted_model.minimizer() == aare::Minimizer::Migrad); - unweighted_model.SetMinimizer(aare::Minimizer::Fumili); - REQUIRE(unweighted_model.minimizer() == aare::Minimizer::Fumili); - - 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, aare::Minimizer::Fumili); - 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)); - for (ssize_t k = 3; k < 6; ++k) - CHECK(weighted(k) > 0.0); - CHECK(weighted(6) == Catch::Approx(0.0).margin(1e-4)); -} - -TEST_CASE("Fumili parameter errors agree with Migrad and Hesse", "[fit]") { - aare::NDArray x({n_points}); - aare::NDArray y({n_points}); - aare::NDArray y_err({n_points}, 2.0); - fill_gaussian(x, y, 120.0, 0.8, 1.3); - - const aare::FitModel migrad(0, 100, 0.5, true, - aare::Minimizer::Migrad); - const aare::FitModel fumili(0, 100, 0.5, true, - aare::Minimizer::Fumili); - - const auto with_migrad = - aare::fit_pixel(migrad, x.view(), y.view(), y_err.view()); - const auto with_fumili = - aare::fit_pixel(fumili, x.view(), y.view(), y_err.view()); - - REQUIRE(with_migrad.size() == 7); - REQUIRE(with_fumili.size() == 7); - for (ssize_t k = 0; k < 3; ++k) - CHECK(with_fumili(k) == Catch::Approx(with_migrad(k)).epsilon(1e-5)); - for (ssize_t k = 3; k < 6; ++k) - CHECK(with_fumili(k) == Catch::Approx(with_migrad(k)).epsilon(1e-2)); -} - -TEST_CASE("Fumili honours fixed parameters and limits", "[fit]") { - aare::NDArray x({n_points}); - aare::NDArray y({n_points}); - fill_gaussian(x, y, 120.0, 0.8, 1.3); - - aare::FitModel model(0, 100, 0.5, false, - aare::Minimizer::Fumili); - model.FixParameter("mu", 0.8); - model.SetParLimits("sigma", 0.5, 5.0); - - const auto result = aare::fit_pixel(model, x.view(), y.view()); - - REQUIRE(result.size() == 4); - CHECK(result(0) == Catch::Approx(120.0).epsilon(1e-5)); - CHECK(result(1) == 0.8); - CHECK(result(2) == Catch::Approx(1.3).epsilon(1e-5)); - CHECK(result(3) == Catch::Approx(0.0).margin(1e-4)); -} - -TEST_CASE("Fumili fits a rising S-curve", "[fit]") { - // f(x) = (p0 + p1*x) + 0.5*(1 + erf((x - mu)/(sqrt(2)*sigma))) * (A + C*(x - // - mu)) - const std::vector truth = {5.0, 0.1, 30.0, 4.0, 200.0, 0.5}; - aare::NDArray x({n_points}); - aare::NDArray y({n_points}); - for (ssize_t i = 0; i < n_points; ++i) { +TEST_CASE("Fit weighted Pol1 errors match the analytic formula", "[fit]") { + const auto minimizer = GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt); + 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) = aare::model::RisingScurve::eval(x(i), truth); + 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 aare::FitModel model( - 0, 200, 0.5, false, aare::Minimizer::Fumili); - const auto result = aare::fit_pixel(model, x.view(), y.view()); + const auto model = make_model(minimizer, true); + const auto res = aare::fit_pixel(model, x.view(), y.view(), y_err.view()); - REQUIRE(result.size() == 7); - for (std::size_t k = 0; k < truth.size(); ++k) - CHECK(result(k) == Catch::Approx(truth[k]).epsilon(1e-3)); + 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("Fumili 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); +TEST_CASE("Fit with a fixed parameter", "[fit]") { + const auto minimizer = GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); - aare::NDArray y({2, 2, n_points}); - aare::NDArray y_err({2, 2, n_points}, 1.0); - 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); - } + 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)); } - aare::NDArray par({2, 2, 3}); - aare::NDArray err({2, 2, 3}); - aare::NDArray chi2({2, 2}); - const aare::FitModel model(0, 100, 0.5, true, - aare::Minimizer::Fumili); + 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)); - aare::fit_3d(model, x.view(), y.view(), y_err.view(), par.view(), - err.view(), chi2.view(), 2); - - // Fumili stops once the estimated remaining chi2 decrease is below - // tolerance * 1e-4, which leaves deviations of order 1e-5 relative here. - 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-4)); - CHECK(par(row, col, 1) == Catch::Approx(-0.6).epsilon(1e-4)); - CHECK(par(row, col, 2) == Catch::Approx(0.9).epsilon(1e-4)); - for (ssize_t k = 0; k < 3; ++k) - CHECK(err(row, col, k) > 0.0); - CHECK(chi2(row, col) == Catch::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); + 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); + 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); + 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("Copying a FitModel keeps the minimizer", "[fit]") { - const aare::FitModel original(0, 100, 0.5, false, - aare::Minimizer::Fumili); + const aare::FitModel original( + 0, 100, 0.5, false, Minimizer::LevenbergMarquardt); const aare::FitModel copied(original); aare::FitModel assigned; assigned = original; - CHECK(copied.minimizer() == aare::Minimizer::Fumili); - CHECK(assigned.minimizer() == aare::Minimizer::Fumili); + 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); + 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); + 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); + 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]") { - aare::NDArray x({n_points}); - aare::NDArray y({n_points}); - fill_gaussian(x, y, 120.0, 0.8, 1.3); + const auto minimizer = GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); - for (const auto minimizer : - {aare::Minimizer::Migrad, aare::Minimizer::Fumili}) { - aare::FitModel model(0, 100, 0.5, true, - minimizer); - model.FixParameter("A", 120.0); - model.FixParameter("mu", 0.8); - model.FixParameter("sigma", 1.3); + 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, x.view(), y.view()); + 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); - CHECK(result(6) == Catch::Approx(0.0).margin(1e-12)); + 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}) { + 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)); + } + } } } diff --git a/src/FitHelpers.hpp b/src/FitHelpers.hpp new file mode 100644 index 00000000..062fa640 --- /dev/null +++ b/src/FitHelpers.hpp @@ -0,0 +1,73 @@ +// 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 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/FitMinuit2.cpp b/src/FitMinuit2.cpp new file mode 100644 index 00000000..4c07040a --- /dev/null +++ b/src/FitMinuit2.cpp @@ -0,0 +1,166 @@ +// 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. + const double span = std::max(1e6, 1e3 * std::abs(start[i])); + if (std::isfinite(lo)) + upar.SetLimits(idx, lo, lo + span); + else + upar.SetLimits(idx, 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/FitMinuit2.hpp b/src/FitMinuit2.hpp new file mode 100644 index 00000000..c0240be8 --- /dev/null +++ b/src/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/LevenbergMarquardt.hpp b/src/LevenbergMarquardt.hpp new file mode 100644 index 00000000..1c88c396 --- /dev/null +++ b/src/LevenbergMarquardt.hpp @@ -0,0 +1,419 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include "aare/FitModel.hpp" +#include "aare/NDView.hpp" + +#include +#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. + * + * - Fixed parameters are left out of the linear system. + * - 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. + * - Damping: Marquardt scaling u * diag(J^T J), starting at u = 0.1, with the + * update rule of Nielsen (Madsen, Nielsen and Tingleff, Methods for + * Non-Linear Least Squares Problems, 2004). + * - Convergence follows Minuit's convention: the fit stops when the + * estimated distance to the minimum of chi2 (EDM, computed from the + * Gauss-Newton Hessian) drops below 0.002 * tolerance, or when the step + * becomes negligible. + * - max_calls limits the number of model evaluations: the initial point, + * every trial point and every accepted point count. As in Minuit2, 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"); + + struct Result { + bool valid = false; + int calls = 0; + int iterations = 0; + double chi2 = 0.0; + }; + + /** + * @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()); + nfree_ = 0; + 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); + p_[k] = start[static_cast(k)]; + if (!model.is_user_fixed(idx)) + free_[nfree_++] = k; + } + for (int a = 0; a < nfree_; ++a) + q_[a] = p_[free_[a]]; + r_.resize(n); + r_new_.resize(n); + J_.resize(n * static_cast(npar)); + + Result res; + auto fail = [&]() { + std::fill(par_out, par_out + npar, 0.0); + if (err_out) + std::fill(err_out, err_out + npar, 0.0); + return res; + }; + + if (!evaluate(q_, r_, true, F_)) + return fail(); + res.calls = 1; + const double edm_target = 0.002 * model.tolerance(); + const int max_calls = model.max_calls() > 0 + ? static_cast(model.max_calls()) + : 200 + 100 * npar + 5 * npar * npar; + double u = 0.1; // damping relative to diag(J^T J) + double nu = 2.0; + + for (;; ++res.iterations) { + build_normal_equations(); + mark_active(); + if (edm() < edm_target || gmax_ < 1e-12) { + res.valid = true; + break; + } + if (res.calls >= max_calls) + break; + if (!solve_step(u)) { + u *= nu; + nu *= 2.0; + if (u > 1e32) + break; + continue; + } + double step2 = 0.0; + double q2 = 0.0; + 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; + const double d = q_new_[a] - q_[a]; + step2 += d * d; + q2 += q_[a] * q_[a]; + } + if (std::sqrt(step2) <= 1e-8 * (std::sqrt(q2) + 1e-8)) { + res.valid = true; + break; + } + double F_new = 0.0; + ++res.calls; + bool ok = evaluate(q_new_, r_new_, false, F_new); + if (crossed) { + // Also try the step truncated at the limit, keep the better. + double F_alt = 0.0; + ++res.calls; + if (evaluate(q_alt_, r_new_, false, F_alt) && + (!ok || F_alt < F_new)) { + q_new_ = q_alt_; + F_new = F_alt; + ok = true; + } + } + const double rho = ok ? (F_ - F_new) / predicted_ : -1.0; + if (rho > 0.0) { + q_ = q_new_; + ++res.calls; + if (!evaluate(q_, r_, true, F_)) + break; // same point as the accepted trial, cannot happen + 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 fail(); + + // p_ tracks the last evaluated point, which may be a rejected trial. + for (int a = 0; a < nfree_; ++a) + p_[free_[a]] = q_[a]; + for (int k = 0; k < npar; ++k) + par_out[k] = p_[k]; + res.chi2 = 2.0 * F_; + + if (err_out) { + std::fill(err_out, err_out + npar, 0.0); + double A[npar][npar]; + double b[npar]; + double inv[npar][npar]; + const int nact = pack(A, b, 0.0); + if (nact > 0 && invert_spd(nact, A, inv)) { + int i = 0; + for (int a = 0; a < nfree_; ++a) { + if (active_[a]) + continue; + err_out[free_[a]] = std::sqrt(std::max(inv[i][i], 0.0)); + ++i; + } + } + } + return res; + } + + private: + double weight(ssize_t i) const { + return weighted_ ? (s_[i] != 0.0 ? 1.0 / s_[i] : 0.0) : 1.0; + } + + // Residuals r = (y - f) / s and, optionally, the Jacobian dr/dp for the + // free parameters at the point q. F = chi2 / 2. Returns false when the + // model rejects the parameters. + bool evaluate(const std::array &q, + std::vector &r, bool with_jacobian, double &F) { + for (int a = 0; a < nfree_; ++a) + p_[free_[a]] = q[a]; + pvec_.assign(p_.begin(), p_.end()); + if (!Model::is_valid(pvec_)) + return false; + const ssize_t n = x_.size(); + F = 0.0; + if (with_jacobian) { + 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 a = 0; a < nfree_; ++a) + Ji[a] = -g[free_[a]] * w; + } + } else { + for (ssize_t i = 0; i < n; ++i) { + r[i] = (y_[i] - Model::eval(x_[i], pvec_)) * weight(i); + F += 0.5 * r[i] * r[i]; + } + } + return true; + } + + void build_normal_equations() { + const ssize_t n = x_.size(); + for (int a = 0; a < nfree_; ++a) { + g_[a] = 0.0; + for (int b = 0; b < nfree_; ++b) + jtj_[a][b] = 0.0; + } + for (ssize_t i = 0; i < n; ++i) { + const double *Ji = &J_[static_cast(i) * npar]; + for (int a = 0; a < nfree_; ++a) { + g_[a] += Ji[a] * r_[i]; + for (int b = a; b < nfree_; ++b) + jtj_[a][b] += Ji[a] * Ji[b]; + } + } + gmax_ = 0.0; + for (int a = 0; a < nfree_; ++a) { + for (int b = 0; b < a; ++b) + jtj_[a][b] = jtj_[b][a]; + 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. + void mark_active() { + for (int a = 0; a < nfree_; ++a) { + const int k = free_[a]; + active_[a] = (p_[k] <= lower_[k] && g_[a] > 0.0) || + (p_[k] >= upper_[k] && g_[a] < 0.0); + } + } + + // Minuit's estimated distance to the minimum of chi2, using the + // Gauss-Newton Hessian: g^T (J^T J)^-1 g with g the gradient of chi2/2. + double edm() { + double A[npar][npar]; + double b[npar]; + double z[npar]; + const int m = pack(A, b, 0.0); + if (m == 0) + return 0.0; + if (!cholesky_solve(m, A, b, z)) + return std::numeric_limits::infinity(); + double e = 0.0; + for (int i = 0; i < m; ++i) + e += b[i] * z[i]; + return e; + } + + // Copy the non-frozen block of J^T J (+ u * diag) and g into A and b. + int pack(double A[npar][npar], double b[npar], 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 (J^T J + u D) delta = -g for the non-frozen + // parameters; also the predicted decrease of F for that step. + bool solve_step(double u) { + double A[npar][npar]; + double b[npar]; + double z[npar]; + const int m = pack(A, b, u); + for (int a = 0; a < nfree_; ++a) + delta_[a] = 0.0; + predicted_ = 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]; + predicted_ += + 0.5 * delta_[a] * + (u * std::max(jtj_[a][a], 1e-300) * delta_[a] - g_[a]); + ++i; + } + return predicted_ > 0.0; + } + + static bool cholesky_solve(int m, const double A[npar][npar], + const double *b, double *z) { + double L[npar][npar]; + 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]; + } + } + } + double t[npar]; + 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]; + } + return true; + } + + static bool invert_spd(int m, const double A[npar][npar], + double inv[npar][npar]) { + double e[npar]; + double col[npar]; + 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; + } + + NDView x_, y_, s_; + bool weighted_ = false; + int nfree_ = 0; + int free_[npar] = {}; + bool active_[npar] = {}; + double lower_[npar] = {}; + double upper_[npar] = {}; + std::array p_{}; + std::array q_{}; + std::array q_new_{}; + std::array q_alt_{}; + double g_[npar] = {}; + double delta_[npar] = {}; + double jtj_[npar][npar] = {}; + double F_ = 0.0; + double gmax_ = 0.0; + double predicted_ = 0.0; + std::vector r_; + std::vector r_new_; + std::vector J_; + std::vector pvec_; +}; + +} // namespace aare::detail diff --git a/src/MinuitSteps.hpp b/src/MinuitSteps.hpp new file mode 100644 index 00000000..d5aa8bbe --- /dev/null +++ b/src/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