diff --git a/RELEASE.md b/RELEASE.md index 0a897d00..8e20c37d 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -29,6 +29,13 @@ - 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. ### API Changes: diff --git a/benchmarks/fit_benchmark.cpp b/benchmarks/fit_benchmark.cpp index 626bde8d..6b537f75 100644 --- a/benchmarks/fit_benchmark.cpp +++ b/benchmarks/fit_benchmark.cpp @@ -75,52 +75,31 @@ static void report_accuracy(benchmark::State &state, const TestCase &tc, // Benchmarks // ---------- -// Minuit2, analytic gradient (no Hesse) -static void BM_FitGausMinuitGrad(benchmark::State &state) { - const auto &tc = get_test_cases()[state.range(0)]; - auto data = generate_gaussian_data(tc); - auto xv = data.x.view(); - auto yv = data.y.view(); - - const auto model = aare::FitModel( - /*strategy = */ 0, - /*max_calls = */ 500, // increase for noisy signals - /*tolerance = */ 0.5, - /*compute_errors = */ false); - - aare::NDArray result; - for (auto _ : state) { - result = aare::fit_pixel(model, xv, yv); - benchmark::DoNotOptimize(result.data()); - } - - report_accuracy(state, tc, result); - state.SetLabel(tc.name); -} - -// Minuit2, analytic gradient + Hesse -static void BM_FitGausMinuitGradHesse(benchmark::State &state) { +static void run_gaussian_fit(benchmark::State &state, aare::Minimizer minimizer, + bool compute_errors) { const auto &tc = get_test_cases()[state.range(0)]; auto data = generate_gaussian_data(tc); auto xv = data.x.view(); auto yv = data.y.view(); auto ev = data.y_err.view(); - const auto model = aare::FitModel( - 0, 500, 0.5, true); // compute_errors = true -> Runs Hesse and provides - // errors on fitted params + const aare::FitModel model( + /*strategy = */ 0, + /*max_calls = */ 500, // increase for noisy signals + /*tolerance = */ 0.5, compute_errors, minimizer); aare::NDArray result; for (auto _ : state) { - result = aare::fit_pixel(model, xv, yv, ev); + result = compute_errors + ? aare::fit_pixel(model, xv, yv, ev) + : aare::fit_pixel(model, xv, yv); benchmark::DoNotOptimize(result.data()); } - // result has 6 elements: [A, mu, sig, err_A, err_mu, err_sig] report_accuracy(state, tc, result); - // Also report Hesse uncertainties - if (result.size() >= 6) { + // With errors the result is [A, mu, sig, err_A, err_mu, err_sig, chi2] + if (compute_errors && result.size() >= 6) { state.counters["errA"] = result(3); state.counters["errMu"] = result(4); state.counters["errSig"] = result(5); @@ -128,12 +107,80 @@ static void BM_FitGausMinuitGradHesse(benchmark::State &state) { state.SetLabel(tc.name); } -BENCHMARK(BM_FitGausMinuitGrad) +// Migrad, analytic gradient (no Hesse) +static void BM_FitGausMigrad(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Migrad, false); +} + +// Migrad, analytic gradient + Hesse for the parameter errors +static void BM_FitGausMigradHesse(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Migrad, true); +} + +// Fumili, analytic gradient and linearised Hessian +static void BM_FitGausFumili(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Fumili, false); +} + +// Fumili, parameter errors from the linearised covariance +static void BM_FitGausFumiliErrors(benchmark::State &state) { + run_gaussian_fit(state, aare::Minimizer::Fumili, true); +} + +// ---------------------------------------------------------------- +// Rising S-curve (6 parameters), typical for threshold scans +// ---------------------------------------------------------------- +static constexpr ssize_t N_SCURVE_POINTS = 100; + +static void run_scurve_fit(benchmark::State &state, aare::Minimizer minimizer) { + const std::vector truth = {5.0, 0.1, 50.0, 4.0, 200.0, 0.5}; + const double noise_sigma = 0.02 * truth[4]; + + aare::NDArray x({N_SCURVE_POINTS}); + aare::NDArray y({N_SCURVE_POINTS}); + std::mt19937 rng(SEED); + std::normal_distribution noise(0.0, noise_sigma); + for (ssize_t i = 0; i < N_SCURVE_POINTS; ++i) { + x[i] = static_cast(i); + y[i] = aare::model::RisingScurve::eval(x[i], truth) + noise(rng); + } + auto xv = x.view(); + auto yv = y.view(); + + const aare::FitModel model(0, 500, 0.5, false, + minimizer); + + aare::NDArray result; + for (auto _ : state) { + result = aare::fit_pixel(model, xv, yv); + benchmark::DoNotOptimize(result.data()); + } + + state.counters["dMu"] = result(2) - truth[2]; + state.counters["dSig"] = result(3) - truth[3]; + state.counters["dA"] = result(4) - truth[4]; + state.counters["chi2"] = result(6); +} + +static void BM_FitScurveMigrad(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::Migrad); +} + +static void BM_FitScurveFumili(benchmark::State &state) { + run_scurve_fit(state, aare::Minimizer::Fumili); +} + +BENCHMARK(BM_FitGausMigrad)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFumili)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond); + +BENCHMARK(BM_FitGausMigradHesse) + ->DenseRange(0, 5) + ->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitGausFumiliErrors) ->DenseRange(0, 5) ->Unit(benchmark::kMicrosecond); -BENCHMARK(BM_FitGausMinuitGradHesse) - ->DenseRange(0, 5) - ->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveMigrad)->Unit(benchmark::kMicrosecond); +BENCHMARK(BM_FitScurveFumili)->Unit(benchmark::kMicrosecond); BENCHMARK_MAIN(); diff --git a/docs/src/pyFit.rst b/docs/src/pyFit.rst index 257f7414..17888299 100644 --- a/docs/src/pyFit.rst +++ b/docs/src/pyFit.rst @@ -27,4 +27,27 @@ For three-dimensional data, pass an array with shape The result dictionary contains ``par`` and ``chi2``. It also contains ``par_err`` when ``compute_errors`` is enabled. +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:: + + model = Gaussian(minimizer=Minimizer.Fumili) + model.minimizer = Minimizer.Migrad + +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. + +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. + .. autofunction:: fit + +.. autoclass:: Minimizer + :members: Migrad, Fumili diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index 1a2fa55f..34402c6f 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -14,6 +14,10 @@ namespace aare { /** * @brief Fit a single pixel's data using Minuit2. * + * 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. + * * User-precedence rules: * - Fixed parameters: untouched (value and fixed flag preserved from model). * - User-set start: value preserved, step size auto-filled. diff --git a/include/aare/FitModel.hpp b/include/aare/FitModel.hpp index a52ea9fa..4670a7cc 100644 --- a/include/aare/FitModel.hpp +++ b/include/aare/FitModel.hpp @@ -10,6 +10,19 @@ namespace aare { +/** + * @brief Minuit2 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 + * 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. + */ +enum class Minimizer { Migrad, Fumili }; + template class FitModel { // Forward declaration only — full definition lives in src/FitModelImpl.hpp @@ -20,6 +33,7 @@ template class FitModel { unsigned int max_calls_; double tolerance_; bool compute_errors_; + Minimizer minimizer_; std::array user_fixed_{}; std::array user_start_{}; @@ -37,10 +51,13 @@ template class FitModel { * default). * @param max_calls Maximum FCN calls per pixel minimisation. * @param tolerance Minuit2 EDM tolerance. - * @param compute_errors If true, run MnHesse after minimisation. + * @param 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. */ explicit FitModel(unsigned int strategy = 0, unsigned int max_calls = 100, - double tolerance = 0.5, bool compute_errors = false); + double tolerance = 0.5, bool compute_errors = false, + Minimizer minimizer = Minimizer::Migrad); // Destructor must be defined in Fit.cpp where FitModelImpl is complete. ~FitModel(); @@ -76,10 +93,12 @@ template class FitModel { void SetMaxCalls(unsigned int n) { max_calls_ = n; } void SetTolerance(double t) { tolerance_ = t; } void SetComputeErrors(bool b) { compute_errors_ = b; } + void SetMinimizer(Minimizer m) { minimizer_ = m; } unsigned int 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]; } diff --git a/python/aare/__init__.py b/python/aare/__init__.py index b4be0755..570176f9 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -40,6 +40,7 @@ from .Cluster import Cluster from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb from ._aare import fit +from ._aare import Minimizer from ._aare import Interpolator from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2 from ._aare import reduce_to_2x2, reduce_to_3x3 diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 82bff325..0a1bfbaf 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -22,9 +22,11 @@ fit_dispatch(const aare::FitModel &model, template void bind_fit_model(py::module &m, const char *name) { using FM = aare::FitModel; py::class_(m, name) - .def(py::init(), + .def(py::init(), py::arg("strategy") = 0, py::arg("max_calls") = 100, - py::arg("tolerance") = 0.5, py::arg("compute_errors") = false) + py::arg("tolerance") = 0.5, py::arg("compute_errors") = false, + py::arg("minimizer") = aare::Minimizer::Migrad) .def("SetParLimits", py::overload_cast(&FM::SetParLimits), py::arg("idx"), py::arg("lo"), py::arg("hi")) @@ -59,6 +61,7 @@ template void bind_fit_model(py::module &m, const char *name) { .def_property("tolerance", &FM::tolerance, &FM::SetTolerance) .def_property("compute_errors", &FM::compute_errors, &FM::SetComputeErrors) + .def_property("minimizer", &FM::minimizer, &FM::SetMinimizer) .def( "__call__", [](const FM & /*self*/, @@ -88,6 +91,9 @@ template void bind_fit_model(py::module &m, const char *name) { R"doc( Fit this model to 1D or 3D data using Minuit2. + The minimizer is selected by the ``minimizer`` property + (``Minimizer.Migrad`` by default, or ``Minimizer.Fumili``). + Parameters ---------- x : array_like, shape (n_scan,) @@ -227,6 +233,22 @@ fit_dispatch(const aare::FitModel &model, } 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. + + ``Migrad`` is the default 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. + )doc") + .value("Migrad", aare::Minimizer::Migrad) + .value("Fumili", aare::Minimizer::Fumili); + // ── Bind model classes ────────────────────────────────────────── bind_fit_model(m, "Gaussian"); bind_fit_model(m, "GaussianErfcPlateau"); @@ -325,6 +347,8 @@ void define_fit_bindings(py::module &m) { 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. 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) diff --git a/python/tests/test_Fit.py b/python/tests/test_Fit.py index 3a477b27..e1b77cdb 100644 --- a/python/tests/test_Fit.py +++ b/python/tests/test_Fit.py @@ -13,3 +13,41 @@ def test_gaussian_model_evaluates_and_fits_data(): result = model.fit(x, y) 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]) + + assert aare.Gaussian().minimizer == aare.Minimizer.Migrad + + model = aare.Gaussian(minimizer=aare.Minimizer.Fumili) + assert model.minimizer == aare.Minimizer.Fumili + + 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 + + +def test_fumili_and_migrad_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): + 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) diff --git a/src/Chi2.hpp b/src/Chi2.hpp index 2aeb0b93..27fcf1af 100644 --- a/src/Chi2.hpp +++ b/src/Chi2.hpp @@ -3,8 +3,9 @@ #include "aare/Models.hpp" #include "aare/NDView.hpp" #include -#include -#include +#include +#include +#include #include #include @@ -13,7 +14,14 @@ namespace aare { namespace func { /** - * @brief Generic chi-squared FCN with analytic gradient for a 1D model. + * @brief Chi-squared objective of a 1D model and its analytic derivatives. + * + * Shared by the Minuit2 FCN adapters below. + * + * chi2(par) = sum_i w_i * (y_i - f(x_i; par))^2 + * + * with w_i = 1/sigma_i^2 for weighted fits (points with sigma_i == 0 are + * skipped) and w_i = 1 otherwise. * * @tparam Model A model struct that satisfies: * - static constexpr std::size_t npar; @@ -22,106 +30,182 @@ namespace func { * double& f, std::array& g); * - static bool is_valid(const std::vector& par); * - * Gradient: - * d(chi2)/dp_k = -2 * sum_i w_i * (y_i - f_i) * df_i/dp_k + * Derivatives, with df_i/dp_k taken from Model::eval_and_grad: * - * where w_i = 1/sigma_i^2 (weighted) or 1 (unweighted). + * d(chi2)/dp_k = -2 * sum_i w_i * (y_i - f_i) * df_i/dp_k + * d2(chi2)/dp_j dp_k ~= 2 * sum_i w_i * df_i/dp_j * df_i/dp_k * - * By providing analytic gradients we avoid 2*npar extra function evaluations - * per Minuit step that would otherwise be spent on finite differences. + * The second line drops the residual-weighted second derivatives of the + * model, which is the Fumili (Gauss-Newton) approximation of the Hessian. + * + * Invalid model parameters do not throw; they return a large penalty (with + * a zero gradient and unit curvature) so the minimizer can remain in control. * * @throws std::invalid_argument if par.size() != Model::npar. + */ +template class Chi2Data1D { + public: + static constexpr std::size_t npar = Model::npar; + /** @brief Number of elements of the packed symmetric Hessian. */ + static constexpr std::size_t nhess = npar * (npar + 1) / 2; + static constexpr double invalid_penalty = 1e20; + + /** @brief An empty y_err view selects an unweighted fit. */ + Chi2Data1D(NDView x, NDView y, + NDView y_err = NDView{}) + : x_(x), y_(y), s_(y_err), weighted_(y_err.size() > 0) {} + + double value(const std::vector &par) const { + check_size(par); + if (!Model::is_valid(par)) + return invalid_penalty; + return weighted_ ? value_impl(par) : value_impl(par); + } + + std::vector gradient(const std::vector &par) const { + check_size(par); + std::vector grad(npar, 0.0); + if (!Model::is_valid(par)) + return grad; + if (weighted_) + gradient_impl(par, grad); + else + gradient_impl(par, grad); + return grad; + } + + /** + * @brief Compute chi2, its gradient and the Gauss-Newton Hessian in a + * single pass over the data. + * + * @param grad Receives d(chi2)/dp_k, resized to npar. + * @param hess Receives the packed symmetric Hessian in the layout used by + * ROOT::Minuit2::FumiliFCNBase: element (j, k) with j <= k is + * stored at index j + k * (k + 1) / 2. Resized to nhess. + * @return chi2 at par. + */ + double evaluate_all(const std::vector &par, + std::vector &grad, + std::vector &hess) const { + check_size(par); + grad.assign(npar, 0.0); + hess.assign(nhess, 0.0); + if (!Model::is_valid(par)) { + for (std::size_t k = 0; k < npar; ++k) + hess[k + k * (k + 1) / 2] = 1.0; + return invalid_penalty; + } + return weighted_ ? evaluate_all_impl(par, grad, hess) + : evaluate_all_impl(par, grad, hess); + } + + private: + static void check_size(const std::vector &par) { + if (par.size() != npar) { + throw std::invalid_argument( + "Chi2Data1D: wrong parameter vector size."); + } + } + + template + double value_impl(const std::vector &par) const { + double chi2 = 0.0; + for (ssize_t i = 0; i < x_.size(); ++i) { + if constexpr (Weighted) { + const double si = s_[i]; + if (si == 0.0) + continue; + const double r_i = y_[i] - Model::eval(x_[i], par); + chi2 += (r_i * r_i) / (si * si); + } else { + const double r_i = y_[i] - Model::eval(x_[i], par); + chi2 += r_i * r_i; + } + } + return chi2; + } + + template + void gradient_impl(const std::vector &par, + std::vector &grad) const { + std::array df{}; + double f_i = 0.0; + for (ssize_t i = 0; i < x_.size(); ++i) { + double c = 0.0; + if constexpr (Weighted) { + const double si = s_[i]; + if (si == 0.0) + continue; + Model::eval_and_grad(x_[i], par, f_i, df); + c = -2.0 * (y_[i] - f_i) / (si * si); + } else { + Model::eval_and_grad(x_[i], par, f_i, df); + c = -2.0 * (y_[i] - f_i); + } + for (std::size_t k = 0; k < npar; ++k) + grad[k] += c * df[k]; + } + } + + template + double evaluate_all_impl(const std::vector &par, + std::vector &grad, + std::vector &hess) const { + std::array df{}; + double f_i = 0.0; + double chi2 = 0.0; + for (ssize_t i = 0; i < x_.size(); ++i) { + double w = 1.0; + if constexpr (Weighted) { + const double si = s_[i]; + if (si == 0.0) + continue; + w = 1.0 / (si * si); + } + Model::eval_and_grad(x_[i], par, f_i, df); + const double r_i = y_[i] - f_i; + chi2 += w * r_i * r_i; + const double c = -2.0 * w * r_i; + for (std::size_t k = 0; k < npar; ++k) { + grad[k] += c * df[k]; + const double h_k = 2.0 * w * df[k]; + double *column = hess.data() + k * (k + 1) / 2; + for (std::size_t j = 0; j <= k; ++j) + column[j] += h_k * df[j]; + } + } + return chi2; + } + + NDView x_; + NDView y_; + NDView s_; + bool weighted_; +}; + +/** + * @brief Chi-squared FCN with analytic gradient, used by Migrad. * - * Invalid model parameters do not throw; they return a large penalty - * (and a zero gradient fallback) so the minimizer can remain in control. + * Providing the analytic gradient avoids the 2*npar extra function + * evaluations per Minuit step that finite differences would need. */ template class Chi2Model1DGrad : public ROOT::Minuit2::FCNGradientBase { public: - Chi2Model1DGrad(NDView x, NDView y) - : x_(x), y_(y), s_(), weighted_(false) {} - + /** @brief An empty y_err view selects an unweighted fit. */ Chi2Model1DGrad(NDView x, NDView y, - NDView y_err) - : x_(x), y_(y), s_(y_err), weighted_(true) {} + NDView y_err = NDView{}) + : data_(x, y, y_err) {} ~Chi2Model1DGrad() override = default; double operator()(const std::vector &par) const override { - if (par.size() != Model::npar) { - throw std::invalid_argument( - "Chi2Model1DGrad: wrong parameter vector size."); - } - - if (!Model::is_valid(par)) - return 1e20; - - double chi2 = 0.0; - - if (weighted_) { - for (ssize_t i = 0; i < x_.size(); ++i) { - const double si = s_[i]; - if (si == 0.0) - continue; - - const double f_i = Model::eval(x_[i], par); - const double r_i = y_[i] - f_i; - chi2 += (r_i * r_i) / (si * si); - } - } else { - for (ssize_t i = 0; i < x_.size(); ++i) { - const double f_i = Model::eval(x_[i], par); - const double r_i = y_[i] - f_i; - chi2 += r_i * r_i; - } - } - - return chi2; + return data_.value(par); } std::vector Gradient(const std::vector &par) const override { - if (par.size() != Model::npar) { - throw std::invalid_argument( - "Chi2Model1DGrad: wrong parameter vector size."); - } - - std::vector grad(Model::npar, 0.0); - if (!Model::is_valid(par)) - return grad; - - std::array df{}; - double f_i = 0.0; - - if (weighted_) { - for (ssize_t i = 0; i < x_.size(); ++i) { - const double si = s_[i]; - if (si == 0.0) - continue; - - Model::eval_and_grad(x_[i], par, f_i, df); - - const double r_i = y_[i] - f_i; - const double c = -2.0 * r_i / (si * si); - - for (std::size_t k = 0; k < Model::npar; ++k) { - grad[k] += c * df[k]; - } - } - } else { - for (ssize_t i = 0; i < x_.size(); ++i) { - Model::eval_and_grad(x_[i], par, f_i, df); - - const double r_i = y_[i] - f_i; - const double c = -2.0 * r_i; - - for (std::size_t k = 0; k < Model::npar; ++k) { - grad[k] += c * df[k]; - } - } - } - - return grad; + return data_.gradient(par); } /** @brief Error definition: 1.0 for chi-squared (delta_chi2 = 1 -> @@ -129,12 +213,41 @@ class Chi2Model1DGrad : public ROOT::Minuit2::FCNGradientBase { double Up() const override { return 1.0; } private: - NDView x_; - NDView y_; - NDView s_; - bool weighted_; + Chi2Data1D data_; +}; + +/** + * @brief Chi-squared FCN for Fumili. + * + * EvaluateAll fills the value, gradient and Gauss-Newton Hessian cached by + * FumiliFCNBase in one pass over the data, so Fumili never needs numerical + * derivatives. + */ +template +class Chi2Model1DFumili : public ROOT::Minuit2::FumiliFCNBase { + public: + /** @brief An empty y_err view selects an unweighted fit. */ + Chi2Model1DFumili(NDView x, NDView y, + NDView y_err = NDView{}) + : ROOT::Minuit2::FumiliFCNBase(static_cast(Model::npar)), + data_(x, y, y_err) {} + + double operator()(const std::vector &par) const override { + return data_.value(par); + } + + /** @brief Error definition: 1.0 for chi-squared (delta_chi2 = 1 -> + * 1-sigma). */ + double Up() const override { return 1.0; } + + void EvaluateAll(const std::vector &par) override { + SetFCNValue(data_.evaluate_all(par, Gradient(), Hessian())); + } + + private: + Chi2Data1D data_; }; } // namespace func -} // namespace aare \ No newline at end of file +} // namespace aare diff --git a/src/Fit.cpp b/src/Fit.cpp index f11aa611..3e4a15b7 100644 --- a/src/Fit.cpp +++ b/src/Fit.cpp @@ -2,9 +2,11 @@ #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 "aare/Models.hpp" #include "aare/utils/par.hpp" @@ -14,6 +16,7 @@ #include #include #include +#include namespace aare { // ============================================================================ @@ -34,9 +37,11 @@ template struct FitModel::FitModelImpl { template FitModel::FitModel(unsigned int strategy, unsigned int max_calls, - double tolerance, bool compute_errors) + double tolerance, bool compute_errors, + Minimizer minimizer) : impl_(std::make_unique(strategy)), max_calls_(max_calls), - tolerance_(tolerance), compute_errors_(compute_errors) { + 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); @@ -57,8 +62,8 @@ template FitModel::FitModel(const FitModel &other) : impl_(std::make_unique(*other.impl_)), max_calls_(other.max_calls_), tolerance_(other.tolerance_), - compute_errors_(other.compute_errors_), user_fixed_(other.user_fixed_), - user_start_(other.user_start_) {} + 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) { @@ -67,6 +72,7 @@ FitModel &FitModel::operator=(const FitModel &other) { 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_; } @@ -144,10 +150,38 @@ std::vector FitModel::GetParNames() const { // fit_pixel / fit_3d — Minuit2 template implementations // ============================================================================ +namespace { + +/** + * @brief Copy the fitted values, optionally the errors, and the objective + * value of a minimum into the flat result layout documented in Fit.hpp. + */ +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}); + + 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]; + } + + result[result_size - 1] = min.Fval(); + return result; +} + +} // namespace + template NDArray fit_pixel(const FitModel &model, NDView x, NDView y, NDView y_err) { - using FCN = func::Chi2Model1DGrad; constexpr std::size_t npar = Model::npar; const bool want_errors = model.compute_errors(); const ssize_t result_size = want_errors ? (2 * npar + 1) : (npar + 1); @@ -175,8 +209,23 @@ NDArray fit_pixel(const FitModel &model, NDView x, upar_local.SetError(i, steps[i]); } - auto chi2 = (y_err.size() > 0) ? FCN(x, y, y_err) : FCN(x, y); + 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()); @@ -184,28 +233,13 @@ NDArray fit_pixel(const FitModel &model, NDView x, if (!min.IsValid()) return NDArray({result_size}, 0.0); - if (want_errors) { + // 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); - - const auto &values = min.UserState().Params(); - const auto &errors = min.UserState().Errors(); - - NDArray result({result_size}); - for (std::size_t k = 0; k < npar; ++k) { - result[k] = values[k]; - result[npar + k] = errors[k]; - } - result[2 * npar] = min.Fval(); - return result; } - const auto &values = min.UserState().Params(); - NDArray result({result_size}); - for (std::size_t k = 0; k < npar; ++k) - result[k] = values[k]; - result[npar] = min.Fval(); - return result; + return pack_minimum(min, want_errors); } template diff --git a/src/Fit.test.cpp b/src/Fit.test.cpp index 2b0fd1bb..97243e14 100644 --- a/src/Fit.test.cpp +++ b/src/Fit.test.cpp @@ -6,6 +6,7 @@ #include #include #include +#include namespace { @@ -79,3 +80,172 @@ TEST_CASE("Minuit2 fits a Gaussian data cube in parallel", "[fit]") { } } } + +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) { + x(i) = static_cast(i); + y(i) = aare::model::RisingScurve::eval(x(i), truth); + } + + const aare::FitModel model( + 0, 200, 0.5, false, aare::Minimizer::Fumili); + const auto result = aare::fit_pixel(model, x.view(), y.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)); +} + +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); + + 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); + } + } + + 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); + + 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)); + } + } +} + +TEST_CASE("Copying a FitModel keeps the minimizer", "[fit]") { + const aare::FitModel original(0, 100, 0.5, false, + aare::Minimizer::Fumili); + const aare::FitModel copied(original); + aare::FitModel assigned; + assigned = original; + + CHECK(copied.minimizer() == aare::Minimizer::Fumili); + CHECK(assigned.minimizer() == aare::Minimizer::Fumili); +} + +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); + + 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); + + const auto result = aare::fit_pixel(model, x.view(), 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)); + } +}