Feature/minuit2 wrapper (#279)

## Unified Minuit2 fitting framework with FitModel API

### Models (`Models.hpp`)
Consolidate all model structs (Gaussian, RisingScurve, FallingScurve)
into a
single header. Each model provides: `eval`, `eval_and_grad`, `is_valid`,
`estimate_par`, `compute_steps`, and `param_info` metadata. No Minuit2
dependency.

### Chi2 functors (`Chi2.hpp`)
Generic `Chi2Model1DGrad` (analytic gradient) templated on the model
struct.
Replaces the separate Chi2Gaussian, Chi2GaussianGradient,
Chi2Scurves, and Chi2ScurvesGradient headers.

### FitModel (`FitModel.hpp`)
Configuration object wrapping `MnUserParameters`, strategy, tolerance,
and
user-override tracking. User constraints (fixed parameters, start
values, limits)
always take precedence over automatic data-driven estimates.

### Fit functions (`Fit.hpp`)
- `fit_pixel<Model, FCN>(model, x, y, y_err)` -> single-pixel,
self-contained
- `fit_pixel<Model, FCN>(model, upar_local, x, y, y_err)` -> pre-cloned
upar for hot loops
- `fit_3d<Model, FCN>(model, x, y, y_err, ..., n_threads)` ->
row-parallel over pixel grid

### Python bindings
- `Pol1`, `Pol2`, `Gaussian`, `RisingScurve`, `FallingScurve` model
classes with
  `FixParameter`, `SetParLimits`, `SetParameter`, and properties for
  `max_calls`, `tolerance`, `compute_errors`
- Single `fit(model, x, y, y_err, n_threads)` dispatch replacing the old
`fit_gaus_minuit`, `fit_gaus_minuit_grad`, `fit_scurve_minuit_grad`,
etc.

### Benchmarks
- Updated `fit_benchmark.cpp` (Google Benchmark) to use the new FitModel
API
- Jupyter notebooks for 1D and 3D S-curve fitting (lmfit vs Minuit2
analytic)
- ~1.8x speedup over lmfit, near-linear thread scaling up to physical
core count

---------

Co-authored-by: Erik Fröjdh <erik.frojdh@psi.ch>

python/CMakeLists.txt

@@ -26,6 +26,10 @@ set_target_properties(_aare PROPERTIES
 )
 target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
 
+target_include_directories(_aare SYSTEM PRIVATE
+    $<TARGET_PROPERTY:Minuit2::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>
+)
+
 # List of python files to be copied to the build directory
 set( PYTHON_FILES
     aare/__init__.py

python/aare/__init__.py

@@ -17,7 +17,8 @@ from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, Clu
 from .ClusterVector import ClusterVector
 from .Cluster import Cluster
 
-
+from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2
+from ._aare import fit
 from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
 from ._aare import Interpolator
 from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2

python/examples/fits.py

@@ -1,8 +1,11 @@
 # SPDX-License-Identifier: MPL-2.0
 import matplotlib.pyplot as plt
 import numpy as np
+import sys
+sys.path.insert(0, '/home/kferjaoui/sw/aare/build')
 from aare import fit_gaus, fit_pol1
-from aare import gaus, pol1
+from aare import Gaussian, fit
+from aare import pol1
 
 textpm = f"±"  #
 textmu = f"μ"  #
@@ -15,35 +18,57 @@ textsigma = f"σ"  #
 mu = np.random.uniform(1, 100)  # Mean of Gaussian
 sigma = np.random.uniform(4, 20)  # Standard deviation
 num_points = 10000  # Number of points for smooth distribution
-noise_sigma = 100
+# noise_sigma = 10
 
 # Generate Gaussian distribution
 data = np.random.normal(mu, sigma, num_points)
+counts, edges = np.histogram(data, bins=100)
 
-# Generate errors for each point
-errors = np.abs(np.random.normal(0, sigma, num_points))  # Errors with mean 0, std 0.5
+x = 0.5 * (edges[:-1] + edges[1:])       # proper bin centers
+y = counts.astype(np.float64)
+
+# Poisson noise
+yerr = np.sqrt(np.maximum(y, 1))
+# yerr = np.abs(np.random.normal(0, noise_sigma, len(x)))
 
 # Create subplot
 fig0, ax0 = plt.subplots(1, 1, num=0, figsize=(12, 8))
 
-x = np.histogram(data, bins=30)[1][:-1] + 0.05
-y = np.histogram(data, bins=30)[0]
-yerr = errors[:30]
-
-
 # Add the errors as error bars in the step plot
 ax0.errorbar(x, y, yerr=yerr, fmt=". ", capsize=5)
 ax0.grid()
 
-par, err = fit_gaus(x, y, yerr)
-print(par, err)
+# Fit with lmfit
+result_lm = fit_gaus(x, y, yerr)
+par_lm = result_lm["par"]
+err_lm = result_lm["par_err"]
+chi2_lm = result_lm["chi2"] 
+print("[lmfit] fit_gaus:            ", par_lm, err_lm, chi2_lm)
+
+# Fit with Minuit2 + analytic gradient + Hesse errors
+gaussian = Gaussian()
+gaussian.compute_errors = True
+result_m2 = gaussian.fit(x, y, yerr)
+par_m2 = result_m2['par']
+err_m2 = result_m2['par_err']
+chi2_m2 = result_m2['chi2']
+print(f"[minuit2] gaussian.fit: par={par_m2}, err={err_m2}, chi2={chi2_m2}")
 
 x = np.linspace(x[0], x[-1], 1000)
-ax0.plot(x, gaus(x, par), marker="")
-ax0.set(xlabel="x", ylabel="Counts", title=f"A0 = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
-                                           f"{textmu} = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
-                                           f"{textsigma} = {par[2]:0.2f}{textpm}{err[2]:0.2f}\n"
-                                           f"(init: {textmu}: {mu:0.2f}, {textsigma}: {sigma:0.2f})")
+ax0.plot(x, gaussian(x, par_lm), marker="", label="fit_gaus")
+ax0.plot(x, gaussian(x, par_m2), marker="", linestyle=":", label="fit_gaus_minuit_grad")
+ax0.legend()
+ax0.set(xlabel="x", ylabel="Counts", 
+    title=(
+        f"fit_gaus:       A={par_lm[0]:0.2f}{textpm}{err_lm[0]:0.2f}  "
+        f"{textmu}={par_lm[1]:0.2f}{textpm}{err_lm[1]:0.2f}  "
+        f"{textsigma}={par_lm[2]:0.2f}{textpm}{err_lm[2]:0.2f}\n"
+        f"minuit_grad: A={par_m2[0]:0.2f}{textpm}{err_m2[0]:0.2f}  "
+        f"{textmu}={par_m2[1]:0.2f}{textpm}{err_m2[1]:0.2f}  "
+        f"{textsigma}={par_m2[2]:0.2f}{textpm}{err_m2[2]:0.2f}\n"
+        f"(truth: {textmu}={mu:0.2f}, {textsigma}={sigma:0.2f})"
+    ),
+)
 fig0.tight_layout()
 
 
@@ -66,8 +91,9 @@ y_values = slope * x_values + intercept + var_points
 
 fig1, ax1 = plt.subplots(1, 1, num=1, figsize=(12, 8))
 ax1.errorbar(x_values, y_values, yerr=errors, fmt=". ", capsize=5)
-par, err = fit_pol1(x_values, y_values, errors)
-
+result_pol = fit_pol1(x_values, y_values, errors)
+par = result_pol["par"]
+err = result_pol["par_err"]
 
 x = np.linspace(np.min(x_values), np.max(x_values), 1000)
 ax1.plot(x, pol1(x, par), marker="")

python/src/fit.hpp

@@ -6,10 +6,204 @@
 #include <pybind11/stl_bind.h>
 
 #include "aare/Fit.hpp"
+#include "aare/FitModel.hpp"
+#include "aare/Models.hpp"
+#include "aare/Chi2.hpp"
 
 namespace py = pybind11;
 using namespace pybind11::literals;
 
+template <typename Model, typename FCN>
+py::object fit_dispatch(
+    const aare::FitModel<Model>& model,
+    py::array_t<double, py::array::c_style | py::array::forcecast> x,
+    py::array_t<double, py::array::c_style | py::array::forcecast> y,
+    py::object y_err_obj,
+    int n_threads);
+
+template <typename Model>
+void bind_fit_model(py::module& m, const char* name) {
+    using FM = aare::FitModel<Model>;
+    using FCN = aare::func::Chi2Model1DGrad<Model>;
+    py::class_<FM>(m, name)
+        .def(py::init<unsigned int, unsigned int, double, bool>(),
+             py::arg("strategy")       = 0,
+             py::arg("max_calls")      = 100,
+             py::arg("tolerance")      = 0.5,
+             py::arg("compute_errors") = false)
+        .def("SetParLimits", &FM::SetParLimits, py::arg("idx"), py::arg("lo"), py::arg("hi"))
+        .def("FixParameter", &FM::FixParameter, py::arg("idx"), py::arg("val"))
+        .def("ReleaseParameter", &FM::ReleaseParameter, py::arg("idx"))
+        .def("SetParameter", &FM::SetParameter, py::arg("idx"), py::arg("val"))
+        .def_property("max_calls", &FM::max_calls, &FM::SetMaxCalls)
+        .def_property("tolerance", &FM::tolerance, &FM::SetTolerance)
+        .def_property("compute_errors", &FM::compute_errors, &FM::SetComputeErrors)
+        .def("__call__",
+            [](const FM& /*self*/,
+                py::array_t<double, py::array::c_style | py::array::forcecast> x,
+                py::array_t<double, py::array::c_style | py::array::forcecast> par)
+            {
+                auto x_view = make_view_1d(x);
+                auto p_view = make_view_1d(par);
+
+                std::vector<double> pvec(p_view.begin(), p_view.end());
+
+                auto* result = new aare::NDArray<double, 1>({x_view.size()});
+                for (ssize_t i = 0; i < x_view.size(); ++i)
+                    (*result)(i) = Model::eval(x_view[i], pvec);
+
+                return return_image_data(result);
+            },
+            py::arg("x"), py::arg("par"))
+        .def("fit",
+            [](const FM& self,
+            py::array_t<double, py::array::c_style | py::array::forcecast> x,
+            py::array_t<double, py::array::c_style | py::array::forcecast> y,
+            py::object y_err_obj,
+            int n_threads) -> py::object
+            {
+                return fit_dispatch<Model, FCN>(self, x, y, y_err_obj, n_threads);
+            },
+            R"doc(
+            Fit this model to 1D or 3D data using Minuit2.
+
+            Parameters
+            ----------
+            x : array_like, shape (n_scan,)
+                Scan points.
+            y : array_like, shape (n_scan,) or (rows, cols, n_scan)
+                Measured data.
+            y_err : array_like or None
+                Per-point uncertainties. None for unweighted fit.
+            n_threads : int
+                Number of threads for 3D parallel loop.
+            )doc",
+            py::arg("x"),
+            py::arg("y"),
+            py::arg("y_err")    = py::none(),
+            py::arg("n_threads") = 4);
+}
+
+template <typename Model>
+py::dict pack_1d_result_dict(const aare::NDArray<double, 1>& result,
+                             bool compute_errors)
+{
+    constexpr std::size_t npar = Model::npar;
+
+    auto res = result.view();
+
+    auto par_out  = new NDArray<double, 1>({npar}, 0.0);
+    auto chi2_out = new NDArray<double, 1>({1},    0.0);
+
+    auto par_view  = par_out->view();
+    auto chi2_view = chi2_out->view();
+
+    for (std::size_t i = 0; i < npar; ++i) {
+        par_view(i) = res(i);
+    }
+
+    if (compute_errors) {
+        auto err_out  = new NDArray<double, 1>({npar}, 0.0);
+        auto err_view = err_out->view();
+
+        for (std::size_t i = 0; i < npar; ++i) {
+            err_view(i) = res(npar + i);
+        }
+
+        chi2_view(0) = res(2 * npar);
+
+        return py::dict(
+            "par"_a     = return_image_data(par_out),
+            "par_err"_a = return_image_data(err_out),
+            "chi2"_a    = return_image_data(chi2_out));
+    } else {
+        chi2_view(0) = res(npar);
+
+        return py::dict(
+            "par"_a  = return_image_data(par_out),
+            "chi2"_a = return_image_data(chi2_out));
+    }
+}
+
+// Helper: typed dispatch for one Model, handles 1D/3D + y_err logic
+template <typename Model, typename FCN>
+py::object fit_dispatch(
+    const aare::FitModel<Model>& model,
+    py::array_t<double, py::array::c_style | py::array::forcecast> x,
+    py::array_t<double, py::array::c_style | py::array::forcecast> y,
+    py::object y_err_obj,
+    int n_threads)
+{
+    constexpr std::size_t npar = Model::npar;
+
+    if (y.ndim() == 3) {
+        auto par_out = new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
+        auto chi2_out= new NDArray<double, 2>({y.shape(0), y.shape(1)}, 0.0);
+
+        auto x_view = make_view_1d(x);
+        auto y_view = make_view_3d(y);
+
+        if (!y_err_obj.is_none()) {
+            auto y_err = py::cast<py::array_t<double,
+                py::array::c_style | py::array::forcecast>>(y_err_obj);
+
+            if (y_err.ndim() != 3) {
+                throw std::runtime_error("For 3D input y, y_err must also be 3D.");
+            }
+            
+            auto err_out = new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
+            auto y_view_err = make_view_3d(y_err);
+
+            aare::fit_3d<Model, FCN>(model, x_view, y_view, y_view_err, 
+                            par_out->view(), err_out->view(), chi2_out->view(), n_threads);
+            
+            if (model.compute_errors()) {
+                return py::dict("par"_a     = return_image_data(par_out),
+                                "par_err"_a = return_image_data(err_out),
+                                "chi2"_a    = return_image_data(chi2_out));
+            } else {
+                delete err_out;
+                return py::dict("par"_a  = return_image_data(par_out),
+                                "chi2"_a = return_image_data(chi2_out));
+            }
+        } else {
+
+            NDView<double, 3> dummy_err{};
+            NDView<double, 3> dummy_err_out{};
+
+            aare::fit_3d<Model, FCN>(model, x_view, y_view, dummy_err, 
+                            par_out->view(), dummy_err_out, chi2_out->view(), n_threads);
+                
+            return py::dict("par"_a = return_image_data(par_out),
+                            "chi2"_a = return_image_data(chi2_out));
+        }
+    } else if (y.ndim() == 1) {
+        NDArray<double, 1> result{};
+
+        auto x_view = make_view_1d(x);
+        auto y_view = make_view_1d(y);
+
+        if (!y_err_obj.is_none()) {
+            auto y_err = py::cast<py::array_t<double,
+                py::array::c_style | py::array::forcecast>>(y_err_obj);
+
+            if (y_err.ndim() != 1) {
+                throw std::runtime_error("For 1D input y, y_err must also be 1D.");
+            }
+
+            auto y_view_err = make_view_1d(y_err);
+            result = aare::fit_pixel<Model, FCN>(model, x_view, y_view, y_view_err);
+        } else {
+            result = aare::fit_pixel<Model, FCN>(model, x_view, y_view);
+        }
+
+        return pack_1d_result_dict<Model>(result, model.compute_errors());
+
+    } else {
+        throw std::runtime_error("Data must be 1D or 3D.");
+    }
+}
+
 void define_fit_bindings(py::module &m) {
 
     // TODO! Evaluate without converting to double
@@ -121,18 +315,17 @@ void define_fit_bindings(py::module &m) {
             }
         },
         R"(
+        Fit a 1D Gaussian to data.
 
-Fit a 1D Gaussian to data.
-
-Parameters
-----------
-x : array_like
-    The x values.
-y : array_like
-    The y values.
-n_threads : int, optional
-    The number of threads to use. Default is 4.
-)",
+        Parameters
+        ----------
+        x : array_like
+            The x values.
+        y : array_like
+            The y values.
+        n_threads : int, optional
+            The number of threads to use. Default is 4.
+        )",
         py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
 
     m.def(
@@ -187,22 +380,22 @@ n_threads : int, optional
             }
         },
         R"(
+        Fit a 1D Gaussian to data with error estimates.
 
-Fit a 1D Gaussian to data with error estimates.
-
-Parameters
-----------
-x : array_like
-    The x values.
-y : array_like
-    The y values.
-y_err : array_like
-    The error in the y values.
-n_threads : int, optional
-    The number of threads to use. Default is 4.
-)",
+        Parameters
+        ----------
+        x : array_like
+            The x values.
+        y : array_like
+            The y values.
+        y_err : array_like
+            The error in the y values.
+        n_threads : int, optional
+            The number of threads to use. Default is 4.
+        )",
         py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
 
+
     m.def(
         "fit_pol1",
         [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
@@ -273,19 +466,19 @@ n_threads : int, optional
             }
         },
         R"(
-Fit a 1D polynomial to data with error estimates.
+        Fit a 1D polynomial to data with error estimates.
 
-Parameters
-----------
-x : array_like
-    The x values.
-y : array_like
-    The y values.
-y_err : array_like
-    The error in the y values.
-n_threads : int, optional
-    The number of threads to use. Default is 4.
-)",
+        Parameters
+        ----------
+        x : array_like
+            The x values.
+        y : array_like
+            The y values.
+        y_err : array_like
+            The error in the y values.
+        n_threads : int, optional
+            The number of threads to use. Default is 4.
+        )",
         py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
 
     //=========
@@ -359,21 +552,22 @@ n_threads : int, optional
             }
         },
         R"(
-Fit a 1D polynomial to data with error estimates.
+        Fit a 1D polynomial to data with error estimates.
 
-Parameters
-----------
-x : array_like
-    The x values.
-y : array_like
-    The y values.
-y_err : array_like
-    The error in the y values.
-n_threads : int, optional
-    The number of threads to use. Default is 4.
-)",
+        Parameters
+        ----------
+        x : array_like
+            The x values.
+        y : array_like
+            The y values.
+        y_err : array_like
+            The error in the y values.
+        n_threads : int, optional
+            The number of threads to use. Default is 4.
+        )",
         py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
 
+
     m.def(
         "fit_scurve2",
         [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
@@ -444,18 +638,105 @@ n_threads : int, optional
             }
         },
         R"(
-Fit a 1D polynomial to data with error estimates.
+        Fit a 1D polynomial to data with error estimates.
 
-Parameters
-----------
-x : array_like
-    The x values.
-y : array_like
-    The y values.
-y_err : array_like
-    The error in the y values.
-n_threads : int, optional
-    The number of threads to use. Default is 4.
-)",
+        Parameters
+        ----------
+        x : array_like
+            The x values.
+        y : array_like
+            The y values.
+        y_err : array_like
+            The error in the y values.
+        n_threads : int, optional
+            The number of threads to use. Default is 4.
+        )",
         py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
+
+
+    // ── Bind model classes ──────────────────────────────────────────
+    bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
+    bind_fit_model<aare::model::RisingScurve>(m, "RisingScurve");
+    bind_fit_model<aare::model::FallingScurve>(m, "FallingScurve");
+    bind_fit_model<aare::model::Pol1>(m, "Pol1");
+    bind_fit_model<aare::model::Pol2>(m, "Pol2");
+
+    m.def("fit",
+        [](py::object model_obj,
+           py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           py::object y_err_obj,
+           int n_threads) -> py::object
+        {
+            using namespace aare::model;
+            using namespace aare::func;
+            
+            // ── Polynomial of degree 1 ───────
+            if(py::isinstance< aare::FitModel<Pol1> >(model_obj)) {
+                const auto& mdl = model_obj.cast< const aare::FitModel<Pol1>& >();
+                return fit_dispatch<Pol1, Chi2Pol1>(mdl, x, y, y_err_obj, n_threads); 
+            }
+
+            // ── Polynomial of degree 2 ───────
+            if(py::isinstance< aare::FitModel<Pol2> >(model_obj)) {
+                const auto& mdl = model_obj.cast< const aare::FitModel<Pol2>& >();
+                return fit_dispatch<Pol2, Chi2Pol2>(mdl, x, y, y_err_obj, n_threads); 
+            }
+            // ── Gaussian ───────
+            if(py::isinstance< aare::FitModel<Gaussian> >(model_obj)) {
+                const auto& mdl = model_obj.cast< const aare::FitModel<Gaussian>& >();
+                return fit_dispatch<Gaussian, Chi2Gaussian>(mdl, x, y, y_err_obj, n_threads); 
+            }
+
+            // ── Rising Scurve ───────
+            if(py::isinstance< aare::FitModel<RisingScurve> >(model_obj)) {
+                const auto& mdl = model_obj.cast< const aare::FitModel<RisingScurve>& >();
+                return fit_dispatch<RisingScurve, Chi2RisingScurve>(mdl, x, y, y_err_obj, n_threads); 
+            }
+
+            // ── Falling Scurve ───────
+            if(py::isinstance< aare::FitModel<FallingScurve> >(model_obj)) {
+                const auto& mdl = model_obj.cast< const aare::FitModel<FallingScurve>& >();
+                return fit_dispatch<FallingScurve, Chi2FallingScurve>(mdl, x, y, y_err_obj, n_threads); 
+            }
+
+            throw std::runtime_error(
+                "Unknown model type. Expected Pol1, Pol2, Gaussian, RisingScurve or FallingScurve."
+            );
+        },
+        R"(
+        Fit a model to 1D or 3D data using Minuit2.
+ 
+        Parameters
+        ----------
+        model : Pol1, Pol2, Gaussian, RisingScurve, or FallingScurve
+            Configured model object.  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)
+            Measured data.  1D for a single pixel, 3D for a detector image.
+        y_err : array_like or None
+            Per-point uncertainties.  Same shape as y.  None → unweighted fit.
+        n_threads : int
+            Number of threads for the 3D parallel loop.
+ 
+        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.
+        )",
+        py::arg("model"),
+        py::arg("x"),
+        py::arg("y"), 
+        py::arg("y_err") = py::none(), 
+        py::arg("n_threads") = 4
+    );
 }