diff --git a/AGENTS.md b/AGENTS.md index c4918858..e5396d4a 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -48,7 +48,9 @@ cmake --build build -j4 Useful optional settings include `AARE_DOCS`, `AARE_BENCHMARKS`, `AARE_ASAN`, `AARE_WARNINGS_AS_ERRORS`, `AARE_VERBOSE` (raises the compile-time log -level), and `AARE_TUNE_LOCAL` (`-march=native`, not portable). Reconfigure an +level), `AARE_TUNE_LOCAL` (`-march=native`, not portable), and `AARE_MINUIT2` +(builds the optional Minuit2 fitting backend for validation; adds an LGPL +dependency and is not used for released packages). Reconfigure an existing build directory instead of creating alternate in-tree build layouts unless isolation is needed; `build/CMakeCache.txt` records the options and interpreter it was configured with. @@ -145,10 +147,19 @@ lives in `python/aare/transform.py` on top of C++ generators. `Fit.hpp`, `FitModel`, and `Models.hpp` (`Gaussian`, `Pol1`, `Pol2`, `RisingScurve`, `FallingScurve`, `GaussianErfcPlateau`, -`GaussianChargeSharing`, ...) wrap Minuit2. Template bodies and explicit -instantiations live in `src/Fit.cpp`; `FitModelImpl` is a pimpl so Minuit2 -headers never leak into the public API. Minuit2 is a private, -`BUILD_INTERFACE`-only dependency of `aare_core`. +`GaussianChargeSharing`, ...) form the public fitting API. Each model +provides `eval`, `eval_and_grad`, `is_valid`, `estimate_par`, and +`param_info`. The minimiser is the dependency-free Levenberg-Marquardt solver +in `src/LevenbergMarquardt.hpp` (Gauss-Newton errors, reflected steps at +limits, Minuit-style EDM tolerance); `src/FitHelpers.hpp` holds the start +value precedence and shape checks; template bodies and explicit +instantiations live in `src/Fit.cpp`. `FitModel` is plain data with the same +public methods as before. The earlier Minuit2 backend survives as an opt-in +validation path: `-DAARE_MINUIT2=ON` compiles `src/FitMinuit2.cpp` (with +`src/Chi2.hpp` and `src/MinuitSteps.hpp`), exposes +`aare::minuit2::fit_pixel/fit_3d` from `aare/FitMinuit2.hpp` and +`aare.experimental.fit_minuit2`, and adds the `[minuit2]` cross-check test. +Run the `[fit]` tests with that option on when changing the solver. ### Python layer diff --git a/CMakeLists.txt b/CMakeLists.txt index e75ea5d4..45dbdf67 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -66,11 +66,18 @@ option(AARE_FETCH_PYBIND11 "Use FetchContent to download pybind11" ON) option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON) option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON) option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON) -option(AARE_FETCH_MINUIT2 "Use FetchContent to download Minuit2" ON) +option(AARE_FETCH_MINUIT2 + "Use FetchContent to download Minuit2 (only with AARE_MINUIT2=ON)" ON) option(AARE_WARNINGS_AS_ERRORS "Treat warnings as errors during compilation" OFF) +# Optional Minuit2 fitting backend, kept for side-by-side validation of the +# built-in Levenberg-Marquardt solver. Adds an LGPL-2.1-or-later dependency and +# is meant for in-tree builds only. +option(AARE_MINUIT2 + "Build the optional Minuit2 fitting backend (validation only)" OFF) + # Convenience option to use system libraries only (no FetchContent) option(AARE_SYSTEM_LIBRARIES "Use system libraries" OFF) if(AARE_SYSTEM_LIBRARIES) @@ -90,8 +97,8 @@ if(AARE_SYSTEM_LIBRARIES) set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE) - # Still fetch Minuit2 when setting AARE_SYSTEM_LIBRARIES since it is not - # available on conda-forge + # AARE_FETCH_MINUIT2 is left untouched: Minuit2 is only used with + # AARE_MINUIT2=ON and is not available on conda-forge endif() if(AARE_BENCHMARKS) @@ -100,47 +107,49 @@ endif() set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -if(AARE_FETCH_MINUIT2) - # We are building Minuit2 from sources. +if(AARE_MINUIT2) + if(AARE_FETCH_MINUIT2) + # We are building Minuit2 from sources. - FetchContent_Declare( - Minuit2 - GIT_REPOSITORY https://github.com/GooFit/Minuit2.git - GIT_TAG v6-40-02 - UPDATE_DISCONNECTED 1) - # Disable Minuit2 extras we don't need - set(minuit2_mpi - OFF - CACHE BOOL "") - set(minuit2_omp - OFF - CACHE BOOL "") - set(BUILD_TESTING - OFF - CACHE BOOL "") - set(MINUIT2_INSTALL - OFF - CACHE BOOL "") - FetchContent_MakeAvailable(Minuit2) + FetchContent_Declare( + Minuit2 + GIT_REPOSITORY https://github.com/GooFit/Minuit2.git + GIT_TAG v6-40-02 + UPDATE_DISCONNECTED 1) + # Disable Minuit2 extras we don't need + set(minuit2_mpi + OFF + CACHE BOOL "") + set(minuit2_omp + OFF + CACHE BOOL "") + set(BUILD_TESTING + OFF + CACHE BOOL "") + set(MINUIT2_INSTALL + OFF + CACHE BOOL "") + FetchContent_MakeAvailable(Minuit2) - set_property(TARGET Minuit2Math PROPERTY POSITION_INDEPENDENT_CODE ON) - set_property(TARGET Minuit2 PROPERTY POSITION_INDEPENDENT_CODE ON) + set_property(TARGET Minuit2Math PROPERTY POSITION_INDEPENDENT_CODE ON) + set_property(TARGET Minuit2 PROPERTY POSITION_INDEPENDENT_CODE ON) - add_library(aare::Minuit2 ALIAS Minuit2) -else() - find_package(Minuit2 CONFIG QUIET) - if(Minuit2_FOUND) - add_library(aare::Minuit2 ALIAS Minuit2::Minuit2) + add_library(aare::Minuit2 ALIAS Minuit2) else() - message( - STATUS - "Stand alone Minuit2 not found. Trying to find it as a part of ROOT") - find_package( - ROOT - COMPONENTS Minuit2 - REQUIRED) - add_library(aare::Minuit2 ALIAS ROOT::Minuit2) - message(STATUS "Found Minuit2 as a part of ROOT: ${ROOT_DIR}") + find_package(Minuit2 CONFIG QUIET) + if(Minuit2_FOUND) + add_library(aare::Minuit2 ALIAS Minuit2::Minuit2) + else() + message( + STATUS + "Stand alone Minuit2 not found. Trying to find it as a part of ROOT") + find_package( + ROOT + COMPONENTS Minuit2 + REQUIRED) + add_library(aare::Minuit2 ALIAS ROOT::Minuit2) + message(STATUS "Found Minuit2 as a part of ROOT: ${ROOT_DIR}") + endif() endif() endif() @@ -358,6 +367,7 @@ set(PUBLICHEADERS include/aare/File.hpp include/aare/Fit.hpp include/aare/FitModel.hpp + include/aare/FitMinuit2.hpp include/aare/Models.hpp include/aare/FileInterface.hpp include/aare/FilePtr.hpp @@ -427,14 +437,20 @@ target_link_libraries( aare_core PUBLIC fmt::fmt nlohmann_json::nlohmann_json ${STD_FS_LIB} # from # helpers.cmake - PRIVATE aare_compiler_flags Threads::Threads $) - -target_include_directories( - aare_core SYSTEM - PRIVATE $) + PRIVATE aare_compiler_flags Threads::Threads) set_property(TARGET aare_core PROPERTY POSITION_INDEPENDENT_CODE ON) +if(AARE_MINUIT2) + target_sources(aare_core + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src/FitMinuit2.cpp) + target_link_libraries(aare_core PRIVATE $) + target_include_directories( + aare_core SYSTEM + PRIVATE $) + target_compile_definitions(aare_core PUBLIC AARE_MINUIT2) +endif() + if(AARE_TESTS) target_compile_definitions(aare_core PRIVATE AARE_TESTS) endif() @@ -489,6 +505,9 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/ROI.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/to_string.test.cpp) + if(AARE_MINUIT2) + list(APPEND TestSources ${CMAKE_CURRENT_SOURCE_DIR}/src/FitMinuit2.test.cpp) + endif() target_sources(tests PRIVATE ${TestSources}) endif() @@ -546,12 +565,38 @@ else() set(CLANG_TIDY_COMMAND "clang-tidy") endif() +# Files that cannot be linted with the default configuration +set(CLANG_TIDY_EXCLUDES + -not + -name + "CircularFifo.hpp" + -not + -name + "ProducerConsumerQueue.hpp" + -not + -name + "VariableSizeClusterFinder.hpp") +if(NOT AARE_MINUIT2) + # The Minuit2 backend sources need the Minuit2 headers + list( + APPEND + CLANG_TIDY_EXCLUDES + -not + -name + "FitMinuit2.cpp" + -not + -name + "Chi2.hpp" + -not + -name + "MinuitSteps.hpp") +endif() + add_custom_target( clang-tidy COMMAND find -path "./src/*" -not -path "./src/python/*" -not -name "*.test.cpp" - -not -name "CircularFifo.hpp" -not -name "ProducerConsumerQueue.hpp" -not - -name "VariableSizeClusterFinder.hpp" | xargs -I {} -n 1 -P 10 bash -c + ${CLANG_TIDY_EXCLUDES} | xargs -I {} -n 1 -P 10 bash -c "${CLANG_TIDY_COMMAND} --config-file=.clang-tidy -p build {}" WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} COMMENT "linting with clang-tidy" diff --git a/README.md b/README.md index 7a9f1c8a..b4de18ea 100644 --- a/README.md +++ b/README.md @@ -19,10 +19,11 @@ Two files contain code from other projects and keep their original licenses: The full license texts are in the `LICENSES` directory. -The Python wheels and conda packages also contain Minuit2 -(LGPL-2.1-or-later), {fmt} and nlohmann/json (MIT), and pybind11 -(BSD-3-Clause), which CMake downloads at build time. See -`THIRD-PARTY-NOTICES.txt` for their notices. +The Python wheels and conda packages also contain {fmt} and nlohmann/json +(MIT) and pybind11 (BSD-3-Clause), which CMake downloads at build time. See +`THIRD-PARTY-NOTICES.txt` for their notices. Building from source with +`-DAARE_MINUIT2=ON` adds the optional Minuit2 fitting backend +(LGPL-2.1-or-later) for validation; released packages do not include it. ## Build and install diff --git a/RELEASE.md b/RELEASE.md index 0a897d00..03f6e800 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -4,6 +4,27 @@ ### New Features: +- Replaced the Minuit2 based fitting backend with a built-in, dependency-free + Levenberg-Marquardt solver that uses the analytic gradients of the models. + The C++ and Python fitting API is unchanged. On 512x1024x200 data cubes the + new solver is 2 to 2.5 times faster, and about 6 times faster with + ``compute_errors``. Parameter errors are Gauss-Newton estimates, + ``sqrt(diag((J^T J)^-1))``: they agree with the previous Hesse errors at a + good minimum and differ by 5 to 15 percent when a parameter is fixed or + limited to a wrong value. Fixed parameters and parameters ending on a limit + report an error of 0 (fixed parameters previously reported 1.0). One-sided + limits are open-ended instead of capped at 1e6, user start and fixed values + are applied before the validity check, ``strategy`` is accepted and + ignored, ``tolerance`` keeps its EDM meaning and ``max_calls`` counts model + evaluations. ``SetParLimits`` rejects ``lo >= hi`` and bad parameter + indices raise ``IndexError``. The public models lost their Minuit specific + ``compute_steps`` and ``compute_ranges`` helpers, and ``FitModel`` no longer + has a pimpl; C++ users need to rebuild. +- The previous Minuit2 implementation remains available for validation when + building from source with ``-DAARE_MINUIT2=ON``: ``aare/FitMinuit2.hpp`` + provides ``aare::minuit2::fit_pixel`` and ``fit_3d``, and Python gains + ``aare.experimental.fit_minuit2``. Released wheels and conda packages are + built without it and no longer contain LGPL-2.1 code. - Added the Python ``FrameDiscardPolicy`` enum with ``NoDiscard``, ``Discard``, and ``DiscardPartial``, enabling access to ``RawMasterFile.frame_discard_policy``. - Added the Python ``Pedestal`` factory with ``dtype`` selection, matching @@ -94,13 +115,11 @@ and the Python wheel declare ``MPL-2.0 AND Apache-2.0 AND MIT`` with all three license files. Building the wheel now requires scikit-build-core 0.11 or newer. -- The wheel and conda package now ship ``THIRD-PARTY-NOTICES.txt`` and the - LGPL 2.1 text for the libraries compiled into the extension: Minuit2 - (LGPL-2.1-or-later), {fmt} and nlohmann/json (MIT), and pybind11 - (BSD-3-Clause). Their license expressions include these licenses. - Minuit2 is now fetched at tag ``v6-40-02`` instead of ``master``, and - libzmq at ``v4.3.5`` (MPL-2.0) instead of ``v4.3.4`` (LGPL-3.0 with a - static-linking exception). +- The wheel and conda package now ship ``THIRD-PARTY-NOTICES.txt`` for the + libraries compiled into the extension: {fmt} and nlohmann/json (MIT) and + pybind11 (BSD-3-Clause). Their license expressions include these licenses. + libzmq is fetched at ``v4.3.5`` (MPL-2.0) instead of ``v4.3.4`` (LGPL-3.0 + with a static-linking exception). - ``RawFile`` and ``File`` reject raw files with frame padding disabled unless the frame discard policy is ``discardpartial``. The constructor reports the master path before opening data subfiles. Legacy ``.raw`` master files now diff --git a/THIRD-PARTY-NOTICES.txt b/THIRD-PARTY-NOTICES.txt index 3d07d9c7..e06a2514 100644 --- a/THIRD-PARTY-NOTICES.txt +++ b/THIRD-PARTY-NOTICES.txt @@ -5,44 +5,15 @@ The aare Python wheels and conda packages compile the following libraries into the aare extension module. They are not part of the aare source repository; CMake downloads them at build time. Their license notices follow below. - Minuit2 v6-40-02 LGPL-2.1-or-later https://github.com/GooFit/Minuit2 {fmt} 12.1.0 MIT https://github.com/fmtlib/fmt nlohmann/json 3.11.3 MIT https://github.com/nlohmann/json pybind11 3.0.3 BSD-3-Clause https://github.com/pybind/pybind11 -Minuit2 is linked statically. To use aare with a modified Minuit2, build aare -from source (https://github.com/slsdetectorgroup/aare) with --DAARE_FETCH_MINUIT2=OFF and make your Minuit2 build findable by CMake's -find_package(Minuit2). - ------------------------------------------------------------------------------- -Minuit2 ------------------------------------------------------------------------------- - -Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT - -This library is free software; you can redistribute it and/or -modify it under the terms of the GNU Lesser General Public -License as published by the Free Software Foundation; either -version 2.1 of the License, or (at your option) any later version. - -This library is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -Lesser General Public License for more details. - -You should have received a copy of the GNU Lesser General Public -License along with this library; if not, write to the Free Software -Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -See LGPL2_1.txt next to this file for the text of the LGPL version 2.1. -For questions regarding the ROOT license, please contact rootdev@cern.ch. - -The precompiled binary distributions provided as a service to our community -may feature packages released with licenses different from LGPL. - -In aare packages the text of the LGPL version 2.1 is in -LICENSES/LGPL-2.1-or-later.txt. +Building aare from source with -DAARE_MINUIT2=ON additionally links the +optional Minuit2 fitting backend (LGPL-2.1-or-later, +https://github.com/GooFit/Minuit2). That option exists for validating the +built-in solver; the published wheels and conda packages are built without it +and contain no Minuit2 code. ------------------------------------------------------------------------------ {fmt} diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index e874a003..36b73a28 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -34,8 +34,5 @@ set_target_properties( add_executable(fit_benchmark fit_benchmark.cpp) target_link_libraries(fit_benchmark PRIVATE benchmark::benchmark aare_core aare_compiler_flags) -target_include_directories( - fit_benchmark SYSTEM - PRIVATE $) set_target_properties(fit_benchmark PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}) diff --git a/benchmarks/fit_benchmark.cpp b/benchmarks/fit_benchmark.cpp index 626bde8d..06e67897 100644 --- a/benchmarks/fit_benchmark.cpp +++ b/benchmarks/fit_benchmark.cpp @@ -2,6 +2,9 @@ #include "aare/Fit.hpp" #include "aare/FitModel.hpp" #include "aare/Models.hpp" +#ifdef AARE_MINUIT2 +#include "aare/FitMinuit2.hpp" +#endif #include #include @@ -71,12 +74,31 @@ static void report_accuracy(benchmark::State &state, const TestCase &tc, state.counters["dSig"] = result(2) - tc.true_sig; } +// Backend tags: the built-in solver, and Minuit2 when compiled in. +struct Builtin { + template + static aare::NDArray fit(Args &&...args) { + return aare::fit_pixel( + std::forward(args)...); + } +}; +#ifdef AARE_MINUIT2 +struct Minuit2 { + template + static aare::NDArray fit(Args &&...args) { + return aare::minuit2::fit_pixel( + std::forward(args)...); + } +}; +#endif + // ---------- // Benchmarks // ---------- -// Minuit2, analytic gradient (no Hesse) -static void BM_FitGausMinuitGrad(benchmark::State &state) { +// Gaussian, unweighted, no errors +template +static void BM_FitGaussian(benchmark::State &state) { const auto &tc = get_test_cases()[state.range(0)]; auto data = generate_gaussian_data(tc); auto xv = data.x.view(); @@ -90,7 +112,7 @@ static void BM_FitGausMinuitGrad(benchmark::State &state) { aare::NDArray result; for (auto _ : state) { - result = aare::fit_pixel(model, xv, yv); + result = Backend::fit(model, xv, yv); benchmark::DoNotOptimize(result.data()); } @@ -98,28 +120,25 @@ static void BM_FitGausMinuitGrad(benchmark::State &state) { state.SetLabel(tc.name); } -// Minuit2, analytic gradient + Hesse -static void BM_FitGausMinuitGradHesse(benchmark::State &state) { +// Gaussian, weighted, with parameter errors +template +static void BM_FitGaussianErrors(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(); 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 auto model = aare::FitModel(0, 500, 0.5, true); aare::NDArray result; for (auto _ : state) { - result = aare::fit_pixel(model, xv, yv, ev); + result = Backend::fit(model, xv, yv, ev); benchmark::DoNotOptimize(result.data()); } - // result has 6 elements: [A, mu, sig, err_A, err_mu, err_sig] + // result has 7 elements: [A, mu, sig, err_A, err_mu, err_sig, chi2] report_accuracy(state, tc, result); - - // Also report Hesse uncertainties if (result.size() >= 6) { state.counters["errA"] = result(3); state.counters["errMu"] = result(4); @@ -128,12 +147,22 @@ static void BM_FitGausMinuitGradHesse(benchmark::State &state) { state.SetLabel(tc.name); } -BENCHMARK(BM_FitGausMinuitGrad) +BENCHMARK_TEMPLATE(BM_FitGaussian, Builtin) ->DenseRange(0, 5) ->Unit(benchmark::kMicrosecond); -BENCHMARK(BM_FitGausMinuitGradHesse) +BENCHMARK_TEMPLATE(BM_FitGaussianErrors, Builtin) ->DenseRange(0, 5) ->Unit(benchmark::kMicrosecond); +#ifdef AARE_MINUIT2 +BENCHMARK_TEMPLATE(BM_FitGaussian, Minuit2) + ->DenseRange(0, 5) + ->Unit(benchmark::kMicrosecond); + +BENCHMARK_TEMPLATE(BM_FitGaussianErrors, Minuit2) + ->DenseRange(0, 5) + ->Unit(benchmark::kMicrosecond); +#endif + BENCHMARK_MAIN(); diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index d2d39b53..1fa55b59 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -50,11 +50,10 @@ test: - python -m pytest python/tests about: - license: MPL-2.0 AND Apache-2.0 AND MIT AND BSD-3-Clause AND LGPL-2.1-or-later + license: MPL-2.0 AND Apache-2.0 AND MIT AND BSD-3-Clause license_file: - LICENSES/MPL-2.0.txt - LICENSES/Apache-2.0.txt - LICENSES/MIT.txt - - LICENSES/LGPL-2.1-or-later.txt - THIRD-PARTY-NOTICES.txt summary: Data analysis library for hybrid pixel detectors from PSI diff --git a/docs/src/Installation.rst b/docs/src/Installation.rst index 2cd8c681..5a24a8ea 100644 --- a/docs/src/Installation.rst +++ b/docs/src/Installation.rst @@ -100,6 +100,15 @@ For detailed options see the CMakeLists.txt file in the root directory of the pr Use system libraries instead of using FetchContent to pull in dependencies. Default option is off. +**AARE_MINUIT2 "Build the optional Minuit2 fitting backend (validation only)" OFF** + +Also build the previous Minuit2 based fitting backend next to the built-in +Levenberg-Marquardt solver, for side-by-side validation +(``aare::minuit2::fit_pixel``, ``aare.experimental.fit_minuit2``). Adds an +LGPL-2.1-or-later dependency, fetched when ``AARE_FETCH_MINUIT2`` is on +(the default) or found with ``find_package(Minuit2)`` or through ROOT. +Released packages are built without it. Default option is off. + **AARE_PYTHON_BINDINGS "Build python bindings" OFF** Build the Python bindings. Default option is off. diff --git a/docs/src/Requirements.rst b/docs/src/Requirements.rst index 7c44dbdc..e621797a 100644 --- a/docs/src/Requirements.rst +++ b/docs/src/Requirements.rst @@ -14,11 +14,13 @@ Requirements To simplify deployment we build and statically link a few libraries. - fmt -- Minuit2 - nlohmann_json - pybind11 - ZeroMQ +Minuit2 is only needed for the optional fitting backend enabled with +``-DAARE_MINUIT2=ON`` (LGPL-2.1-or-later); it is not part of the default build. + **Extra dependencies for building documentation** - Sphinx diff --git a/docs/src/pyFit.rst b/docs/src/pyFit.rst index 257f7414..5d7f742c 100644 --- a/docs/src/pyFit.rst +++ b/docs/src/pyFit.rst @@ -3,8 +3,9 @@ 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:: +Aare fits one-dimensional scans and three-dimensional pixel data with a +built-in Levenberg-Marquardt solver that uses the analytic gradients of the +models. Create a model object and call its :meth:`fit` method:: model = Gaussian(compute_errors=True) result = model.fit(x, y, y_err) @@ -27,4 +28,30 @@ 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. +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 constructor arguments control the minimiser: + +- ``max_calls`` caps the number of model evaluations per pixel. A pixel that + does not converge within the budget returns zeros. +- ``tolerance`` is the EDM tolerance in Minuit's convention: the fit stops + when the estimated distance to the minimum is below ``0.002 * tolerance``. +- ``strategy`` is accepted for compatibility and ignored. +- ``compute_errors`` returns Gauss-Newton parameter errors, + ``sqrt(diag((J^T J)^-1))``. Fixed parameters and parameters that end on a + limit report an error of 0. + +Validation backend +^^^^^^^^^^^^^^^^^^ + +When aare is built from source with ``-DAARE_MINUIT2=ON``, the previous +Minuit2 (Migrad and Hesse) implementation is available as +``aare.experimental.fit_minuit2(model, x, y, y_err=None, n_threads=4)`` with +the same arguments and result as :func:`fit`. It exists to compare the two +solvers and is not part of the released packages. + .. autofunction:: fit diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index 1a2fa55f..f9a4fcd8 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -7,12 +7,13 @@ namespace aare { // --------------------------------------------------------------------------- -// Minuit2-based pixel fitting. +// Levenberg-Marquardt pixel fitting. // 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 with the built-in Levenberg-Marquardt + * solver. * * User-precedence rules: * - Fixed parameters: untouched (value and fixed flag preserved from model). diff --git a/include/aare/FitMinuit2.hpp b/include/aare/FitMinuit2.hpp new file mode 100644 index 00000000..29d3475a --- /dev/null +++ b/include/aare/FitMinuit2.hpp @@ -0,0 +1,33 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include "aare/FitModel.hpp" +#include "aare/NDArray.hpp" + +namespace aare::minuit2 { + +// --------------------------------------------------------------------------- +// Optional Minuit2 (Migrad + Hesse) fitting backend, kept for side-by-side +// validation of the built-in Levenberg-Marquardt solver. +// +// The definitions are only compiled into aare_core when CMake is configured +// with -DAARE_MINUIT2=ON (compile definition AARE_MINUIT2); without it these +// functions fail to link. Arguments, result layout and failure behaviour are +// the same as for aare::fit_pixel and aare::fit_3d, see aare/Fit.hpp. +// --------------------------------------------------------------------------- + +template +NDArray fit_pixel(const FitModel &model, NDView x, + NDView y, NDView y_err); + +template +NDArray fit_pixel(const FitModel &model, NDView x, + NDView y); + +template +void fit_3d(const FitModel &model, NDView x, + NDView y, NDView y_err, + NDView par_out, NDView err_out, + NDView chi2_out, int n_threads); + +} // namespace aare::minuit2 diff --git a/include/aare/FitModel.hpp b/include/aare/FitModel.hpp index a52ea9fa..2950ca0b 100644 --- a/include/aare/FitModel.hpp +++ b/include/aare/FitModel.hpp @@ -3,59 +3,54 @@ #include "aare/Models.hpp" #include -#include #include #include #include namespace aare { +/** + * @brief Fit configuration for one model: parameter limits, fixed + * parameters, user start values and minimiser settings. + * + * The object holds plain data and can be shared read-only between threads. + * Method bodies live in src/Fit.cpp, which explicitly instantiates the class + * for every supported model. + */ template class FitModel { - - // Forward declaration only — full definition lives in src/FitModelImpl.hpp - // and is never installed. Downstream code cannot dereference this pointer. - struct FitModelImpl; // stores Minuit2 user-set parameters - - std::unique_ptr impl_; - unsigned int max_calls_; - double tolerance_; - bool compute_errors_; - - std::array user_fixed_{}; - std::array user_start_{}; - - /** @brief Safely resolve a parameter name to its index. */ - unsigned int checked_index(const std::string &name) const; - public: static constexpr std::size_t npar = Model::npar; /** * @brief Construct a fit model with sensible defaults. * - * @param strategy Minuit2 strategy level (0 = fast/gradient, 1 = - * default). - * @param max_calls Maximum FCN calls per pixel minimisation. - * @param tolerance Minuit2 EDM tolerance. - * @param compute_errors If true, run MnHesse after minimisation. + * @param strategy Accepted for compatibility with the earlier + * Minuit2 based implementation. Ignored by the + * built-in solver and forwarded to the optional + * Minuit2 validation backend. + * @param max_calls Maximum number of model evaluations per pixel. + * @param tolerance EDM tolerance in Minuit's convention: the fit + * stops when the estimated distance to the + * minimum drops below 0.002 * tolerance. + * @param compute_errors If true, return Gauss-Newton parameter errors, + * sqrt(diag((J^T J)^-1)). Fixed parameters and + * parameters ending on a limit report 0. */ explicit FitModel(unsigned int strategy = 0, unsigned int max_calls = 100, double tolerance = 0.5, bool compute_errors = false); - // 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. 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); @@ -77,15 +72,32 @@ template class FitModel { void SetTolerance(double t) { tolerance_ = t; } void SetComputeErrors(bool b) { compute_errors_ = b; } + 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_; } 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_; + 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/pyproject.toml b/pyproject.toml index 6054356c..9f735557 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,7 @@ dependencies = [ "numpy", "matplotlib", ] -license = "MPL-2.0 AND Apache-2.0 AND MIT AND BSD-3-Clause AND LGPL-2.1-or-later" +license = "MPL-2.0 AND Apache-2.0 AND MIT AND BSD-3-Clause" license-files = ["LICENSES/*", "THIRD-PARTY-NOTICES.txt"] diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 7f4db4c0..01052b80 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -27,10 +27,6 @@ set_target_properties(_aare PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}) target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags) -target_include_directories( - _aare SYSTEM - PRIVATE $) - # List of python files to be copied to the build directory set(PYTHON_FILES aare/__init__.py diff --git a/python/aare/experimental.py b/python/aare/experimental.py index 2a5210d9..a3ef22a2 100644 --- a/python/aare/experimental.py +++ b/python/aare/experimental.py @@ -4,3 +4,10 @@ from ._aare.experimental import MultiThreadedFileReader __all__ = ["MultiThreadedFileReader"] + +try: + from ._aare.experimental import fit_minuit2 + + __all__.append("fit_minuit2") +except ImportError: # built without AARE_MINUIT2 + pass diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 82bff325..30da7d91 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -8,11 +8,41 @@ #include "aare/Fit.hpp" #include "aare/FitModel.hpp" #include "aare/Models.hpp" +#ifdef AARE_MINUIT2 +#include "aare/FitMinuit2.hpp" +#endif namespace py = pybind11; using namespace pybind11::literals; -template +// Backend tags select which fit functions serve a Python call. The built-in +// Levenberg-Marquardt solver is always available; the Minuit2 backend only +// when aare was configured with AARE_MINUIT2=ON. +struct LmBackend { + template + static aare::NDArray fit_pixel(Args &&...args) { + return aare::fit_pixel(std::forward(args)...); + } + template + static void fit_3d(Args &&...args) { + aare::fit_3d(std::forward(args)...); + } +}; + +#ifdef AARE_MINUIT2 +struct Minuit2Backend { + template + static aare::NDArray fit_pixel(Args &&...args) { + return aare::minuit2::fit_pixel(std::forward(args)...); + } + template + static void fit_3d(Args &&...args) { + aare::minuit2::fit_3d(std::forward(args)...); + } +}; +#endif + +template py::object fit_dispatch(const aare::FitModel &model, py::array_t x, @@ -24,7 +54,23 @@ template void bind_fit_model(py::module &m, const char *name) { py::class_(m, name) .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, + R"doc( + Fit configuration for this model. + + Parameters + ---------- + strategy : int + Accepted for compatibility and ignored by the built-in solver. + max_calls : int + Maximum number of model evaluations per pixel. + tolerance : float + EDM tolerance (Minuit convention): the fit stops when the + estimated distance to the minimum is below 0.002 * tolerance. + compute_errors : bool + Also return Gauss-Newton parameter errors. Fixed parameters + and parameters ending on a limit report an error of 0. + )doc") .def("SetParLimits", py::overload_cast(&FM::SetParLimits), py::arg("idx"), py::arg("lo"), py::arg("hi")) @@ -83,10 +129,12 @@ template void bind_fit_model(py::module &m, const char *name) { py::array_t x, py::array_t y, py::object y_err_obj, int n_threads) -> py::object { - return fit_dispatch(self, x, y, y_err_obj, n_threads); + return fit_dispatch(self, x, y, y_err_obj, + n_threads); }, R"doc( - Fit this model to 1D or 3D data using Minuit2. + Fit this model to 1D or 3D data with the built-in + Levenberg-Marquardt solver. Parameters ---------- @@ -142,7 +190,7 @@ py::dict pack_1d_result_dict(const aare::NDArray &result, } // Helper: typed dispatch for one Model, handles 1D/3D + y_err logic -template +template py::object fit_dispatch(const aare::FitModel &model, py::array_t x, @@ -172,9 +220,9 @@ fit_dispatch(const aare::FitModel &model, new NDArray({y.shape(0), y.shape(1), npar}, 0.0); auto y_view_err = make_view_3d(y_err); - aare::fit_3d(model, x_view, y_view, y_view_err, - par_out->view(), err_out->view(), - chi2_out->view(), n_threads); + Backend::template fit_3d(model, x_view, y_view, y_view_err, + par_out->view(), err_out->view(), + chi2_out->view(), n_threads); if (model.compute_errors()) { return py::dict("par"_a = return_image_data(par_out), @@ -190,9 +238,9 @@ fit_dispatch(const aare::FitModel &model, NDView dummy_err{}; NDView dummy_err_out{}; - aare::fit_3d(model, x_view, y_view, dummy_err, - par_out->view(), dummy_err_out, - chi2_out->view(), n_threads); + Backend::template fit_3d(model, x_view, y_view, dummy_err, + par_out->view(), dummy_err_out, + chi2_out->view(), n_threads); return py::dict("par"_a = return_image_data(par_out), "chi2"_a = return_image_data(chi2_out)); @@ -214,9 +262,10 @@ fit_dispatch(const aare::FitModel &model, } auto y_view_err = make_view_1d(y_err); - result = aare::fit_pixel(model, x_view, y_view, y_view_err); + result = Backend::template fit_pixel(model, x_view, y_view, + y_view_err); } else { - result = aare::fit_pixel(model, x_view, y_view); + result = Backend::template fit_pixel(model, x_view, y_view); } return pack_1d_result_dict(result, model.compute_errors()); @@ -226,7 +275,39 @@ fit_dispatch(const aare::FitModel &model, } } -void define_fit_bindings(py::module &m) { +// Resolve the model type behind a Python object and dispatch to Backend. +template +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, py::module &experimental) { // ── Bind model classes ────────────────────────────────────────── bind_fit_model(m, "Gaussian"); bind_fit_model(m, "GaussianErfcPlateau"); @@ -245,86 +326,20 @@ 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 built-in Levenberg-Marquardt + solver. Parameters ---------- - model : Pol1, Pol2, Gaussian, RisingScurve, or FallingScurve - Configured model object. User-set limits, fixed parameters, - and start values take precedence over automatic estimates. + model : object + Configured model object: Pol1, Pol2, Gaussian, GaussianErfcPlateau, + GaussianChargeSharing, GaussianChargeSharingKb, RisingScurve or + FallingScurve. User-set limits, fixed parameters, and start values + take precedence over automatic estimates. x : array_like, shape (n_scan,) Scan points (e.g. energy or threshold values). y : array_like, shape (n_scan,) or (rows, cols, n_scan) @@ -336,16 +351,36 @@ 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 Gauss-Newton + 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); + +#ifdef AARE_MINUIT2 + experimental.def( + "fit_minuit2", + [](py::object model_obj, + py::array_t x, + py::array_t y, + py::object y_err_obj, int n_threads) -> py::object { + return dispatch_any_model(model_obj, x, y, + y_err_obj, n_threads); + }, + R"( + Validation only: fit with the Minuit2 (Migrad + Hesse) backend. + + Available when aare was built with -DAARE_MINUIT2=ON. Takes the same + arguments and returns the same dictionary as fit(). + )", + py::arg("model"), py::arg("x"), py::arg("y"), + py::arg("y_err") = py::none(), py::arg("n_threads") = 4); +#else + (void)experimental; +#endif } diff --git a/python/src/module.cpp b/python/src/module.cpp index b83deaca..a23d4c2a 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -84,7 +84,7 @@ PYBIND11_MODULE(_aare, m) { define_fast_pedestal_bindings(m, "FastPedestal_d"); define_fast_pedestal_bindings(m, "FastPedestal_f"); define_fast_pedestal_bindings(m, "FastPedestal_i16"); - define_fit_bindings(m); + define_fit_bindings(m, experimental); define_interpolation_bindings(m); define_jungfrau_data_file_io_bindings(m); diff --git a/python/tests/test_Fit.py b/python/tests/test_Fit.py index 3a477b27..551f65a5 100644 --- a/python/tests/test_Fit.py +++ b/python/tests/test_Fit.py @@ -1,8 +1,21 @@ # SPDX-License-Identifier: MPL-2.0 import numpy as np +import pytest + import aare +def gaussian(x, A, mu, sigma): + return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2) + + +@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(): x = np.linspace(-5.0, 5.0, 51) @@ -13,3 +26,159 @@ 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_fit_returns_errors_and_chi2(gaussian_data): + x, y, y_err = gaussian_data + model = aare.Gaussian(compute_errors=True) + 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(), x, y) + assert set(without_errors) == {"par", "chi2"} + + +def test_pol1_weighted_errors_match_analytic_formula(): + 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).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): + x, y, y_err = gaussian_data + model = aare.Gaussian(compute_errors=True) + 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): + x, y, _ = gaussian_data + model = aare.Gaussian() + model.SetParLimits("mu", -5.0, 0.5) + result = model.fit(x, y) + + assert result["par"][1] == pytest.approx(0.5, abs=1e-9) + 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() + 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): + x, y, _ = gaussian_data + model = aare.Gaussian() + 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): + x, y, y_err = gaussian_data + result = aare.Gaussian(max_calls=3, compute_errors=True).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_accepted_and_ignored(gaussian_data): + x, y, _ = gaussian_data + a = aare.Gaussian(strategy=0).fit(x, y) + b = aare.Gaussian(strategy=1).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): + 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) + 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): + 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) + 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_minuit2_backend_matches_builtin_solver(gaussian_data): + fit_minuit2 = getattr(aare.experimental, "fit_minuit2", None) + if fit_minuit2 is None: + pytest.skip("aare was built without AARE_MINUIT2") + x, y, y_err = gaussian_data + model = aare.Gaussian(compute_errors=True) + model.SetParLimits("mu", -5.0, 0.5) + lm = model.fit(x, y, y_err) + mn = fit_minuit2(model, x, y, y_err) + np.testing.assert_allclose(lm["par"], mn["par"], rtol=1e-3) + np.testing.assert_allclose(lm["chi2"], mn["chi2"], rtol=1e-3) diff --git a/src/Fit.cpp b/src/Fit.cpp index f11aa611..a11a142d 100644 --- a/src/Fit.cpp +++ b/src/Fit.cpp @@ -1,82 +1,41 @@ // SPDX-License-Identifier: MPL-2.0 #include "aare/Fit.hpp" -#include "Chi2.hpp" -#include "Minuit2/FunctionMinimum.h" -#include "Minuit2/MnHesse.h" -#include "Minuit2/MnMigrad.h" -#include "Minuit2/MnStrategy.h" -#include "Minuit2/MnUserParameters.h" +#include "FitHelpers.hpp" +#include "LevenbergMarquardt.hpp" #include "aare/Models.hpp" #include "aare/utils/par.hpp" #include "aare/utils/task.hpp" #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) - : impl_(std::make_unique(strategy)), max_calls_(max_calls), - tolerance_(tolerance), compute_errors_(compute_errors) { + : strategy_(strategy), max_calls_(max_calls), tolerance_(tolerance), + compute_errors_(compute_errors) { for (std::size_t i = 0; i < npar; ++i) { - const auto pi = Model::param_info[i]; - const bool has_lo = std::isfinite(pi.default_lo); - const bool has_hi = std::isfinite(pi.default_hi); - if (has_lo && has_hi) { - impl_->upar.Add(pi.name, 0.0, 1.0, pi.default_lo, pi.default_hi); - } else if (has_lo) { - impl_->upar.Add(pi.name, 0.0, 1.0, pi.default_lo, 1e6); - } else { - impl_->upar.Add(pi.name, 0.0, 1.0); - } + 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_), user_fixed_(other.user_fixed_), - user_start_(other.user_start_) {} - -template -FitModel &FitModel::operator=(const FitModel &other) { - if (this != &other) { - impl_ = std::make_unique(*other.impl_); - max_calls_ = other.max_calls_; - tolerance_ = other.tolerance_; - compute_errors_ = other.compute_errors_; - user_fixed_ = other.user_fixed_; - user_start_ = other.user_start_; - } - return *this; +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 + "'"); @@ -84,19 +43,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; } @@ -107,7 +71,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; } @@ -129,82 +94,47 @@ 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 — Levenberg-Marquardt template implementations // ============================================================================ template NDArray fit_pixel(const FitModel &model, NDView x, NDView y, NDView y_err) { - using FCN = func::Chi2Model1DGrad; constexpr std::size_t npar = Model::npar; const bool want_errors = model.compute_errors(); - const ssize_t result_size = want_errors ? (2 * npar + 1) : (npar + 1); + 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); - - if (!Model::is_valid(std::vector(start.begin(), start.end()))) { - return NDArray({result_size}, 0.0); - } - - double x_range, y_range, slope_scale; - model::compute_ranges(x, y, x_range, y_range, slope_scale); - - std::array steps{}; - Model::compute_steps(start, x_range, y_range, slope_scale, steps); - - // thread-local copy of starting parameters - auto upar_local = model.impl()->upar; - - for (std::size_t i = 0; i < npar; ++i) { - if (model.is_user_fixed(i)) - continue; - if (!model.is_user_start(i)) - upar_local.SetValue(i, start[i]); - upar_local.SetError(i, steps[i]); - } - - auto chi2 = (y_err.size() > 0) ? FCN(x, y, y_err) : FCN(x, y); - - ROOT::Minuit2::MnMigrad migrad(chi2, upar_local, model.impl()->strategy); - ROOT::Minuit2::FunctionMinimum min = - migrad(model.max_calls(), model.tolerance()); - - if (!min.IsValid()) - return NDArray({result_size}, 0.0); + const auto start = detail::start_values(model, x, y); + detail::LevenbergMarquardt solver; + std::array err{}; + const auto res = solver.fit(model, x, y, y_err, start, result.data(), + want_errors ? err.data() : nullptr); + if (!res.valid) + return result; // all zeros, like the previous implementation if (want_errors) { - ROOT::Minuit2::MnHesse hesse; - hesse(chi2, min); - - const auto &values = min.UserState().Params(); - const auto &errors = min.UserState().Errors(); - - NDArray result({result_size}); - for (std::size_t k = 0; k < npar; ++k) { - result[k] = values[k]; - result[npar + k] = errors[k]; - } - result[2 * npar] = min.Fval(); - return result; + for (std::size_t k = 0; k < npar; ++k) + result(static_cast(npar + k)) = err[k]; + result(static_cast(2 * npar)) = res.chi2; + } else { + result(static_cast(npar)) = res.chi2; } - - 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; } @@ -219,55 +149,36 @@ 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::LevenbergMarquardt solver; // buffers reused per 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); + const auto start = detail::start_values(model, x, values); + const auto res = + solver.fit(model, x, values, errors, start, par.data(), + want_par_errors ? err.data() : nullptr); 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) = res.valid ? res.chi2 : 0.0; } } }; diff --git a/src/Fit.test.cpp b/src/Fit.test.cpp index 2b0fd1bb..81f6cc8a 100644 --- a/src/Fit.test.cpp +++ b/src/Fit.test.cpp @@ -6,12 +6,26 @@ #include #include #include +#include +#include +#include +#include + +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); @@ -20,42 +34,301 @@ void fill_gaussian(aare::NDArray &x, aare::NDArray &y, } } +// Noise-free Gaussian with A = 120, mu = 0.8, sigma = 1.3 and unit errors. +Data gaussian_data() { + 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); + return d; +} + +// Noise-free S-curve on 100 points in [0, 100] with unit errors. +template Data scurve_data(const std::vector &par) { + 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); + } + return d; +} + } // 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]") { + auto d = gaussian_data(); const aare::FitModel unweighted_model; 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 = - 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)); + // Gauss-Newton errors equal the exact Hessian errors at zero residual + 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}); +TEST_CASE("Fit weighted Pol1 errors match the analytic formula", "[fit]") { + constexpr ssize_t n = 20; + NDArray x({n}); + NDArray y({n}); + NDArray y_err({n}); + double S = 0.0, Sx = 0.0, Sxx = 0.0; + for (ssize_t i = 0; i < n; ++i) { + x(i) = static_cast(i); + y(i) = 2.0 + 0.5 * x(i); + y_err(i) = 1.0 + 0.1 * static_cast(i); + const double w = 1.0 / (y_err(i) * y_err(i)); + S += w; + Sx += w * x(i); + Sxx += w * x(i) * x(i); + } + const double delta = S * Sxx - Sx * Sx; + + const aare::FitModel model(0, 100, 0.5, true); + const auto res = aare::fit_pixel(model, x.view(), y.view(), y_err.view()); + + REQUIRE(res.size() == 5); + CHECK(res(0) == Approx(2.0).epsilon(1e-6)); + CHECK(res(1) == Approx(0.5).epsilon(1e-6)); + CHECK(res(2) == Approx(std::sqrt(Sxx / delta)).epsilon(1e-6)); + CHECK(res(3) == Approx(std::sqrt(S / delta)).epsilon(1e-6)); + CHECK(res(4) == Approx(0.0).margin(1e-8)); +} + +TEST_CASE("Fit with a fixed parameter", "[fit]") { + auto d = gaussian_data(); + + SECTION("fixed at a wrong value") { + aare::FitModel model(0, 100, 0.5, 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(5e-2)); + CHECK(res(5) == 0.0); // fixed parameters report no error + CHECK(res(6) == Approx(1684.25).epsilon(1e-3)); + } + + SECTION("fixed at the true value, then released") { + aare::FitModel model(0, 100, 0.5, true); + model.FixParameter(2, 1.3); + auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(0) == Approx(120.0).epsilon(1e-4).margin(1e-6)); + CHECK(res(1) == Approx(0.8).epsilon(1e-4).margin(1e-6)); + CHECK(res(2) == 1.3); + CHECK(res(5) == 0.0); + CHECK(res(6) == Approx(0.0).margin(1e-3)); + + model.ReleaseParameter("sigma"); + REQUIRE_FALSE(model.is_user_fixed(2)); + res = aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(2) == Approx(1.3).epsilon(1e-4).margin(1e-6)); + CHECK(res(5) == Approx(0.00451).epsilon(2e-2)); + } +} + +TEST_CASE("Fit with parameter limits", "[fit]") { + auto d = gaussian_data(); + + SECTION("active upper limit on mu") { + aare::FitModel model; + 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(1e-9)); + 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") { + aare::FitModel model(0, 100, 0.5, 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(1e-9)); + CHECK(res(1) == Approx(0.8).epsilon(1e-3)); + CHECK(res(3) == Approx(0.336).epsilon(5e-2)); + 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") { + aare::FitModel model(0, 100, 0.5, 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(1e-9)); + 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(5e-2)); + CHECK(res(6) == Approx(3631.01).epsilon(1e-3)); + } + + SECTION("inactive limits do not change the result") { + const aare::FitModel free_model; + const auto reference = + aare::fit_pixel(free_model, d.x.view(), d.y.view()); + + aare::FitModel model; + 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]") { + auto d = gaussian_data(); + + SECTION("start value for mu") { + aare::FitModel model; + 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); + aare::FitModel model; + 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]") { + auto d = gaussian_data(); + + const aare::FitModel model(0, 3, 0.5, false); + 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); + + const aare::FitModel model_err(0, 3, 0.5, true); + 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.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); +} + +TEST_CASE("Fit rising and falling S-curves", "[fit]") { + 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 aare::FitModel model(0, 500, 0.5, + true); + 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); + aare::FitModel model(0, 500, 0.5, true); + 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); + aare::FitModel model(0, 500, 0.5, false); + 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(1e-9)); + CHECK(res(6) > 0.0); // chi2 of the constrained fit + } + + SECTION("falling") { + auto d = scurve_data(truth); + const aare::FitModel model(0, 500, 0.5, + false); + 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]") { + NDArray x({n_points}); + NDArray values({n_points}); fill_gaussian(x, values, 80.0, -0.6, 0.9); - aare::NDArray y({2, 2, n_points}); + 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) @@ -63,19 +336,77 @@ TEST_CASE("Minuit2 fits a Gaussian data cube in parallel", "[fit]") { } } - aare::NDArray par({2, 2, 3}); - aare::NDArray chi2({2, 2}); + NDArray par({2, 2, 3}); + 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); + 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) == 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)); + 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]") { + 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 aare::FitModel model(0, 100, 0.5, 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)); + CHECK(par(row, col, 0) == + Approx(80.0 + 10.0 * static_cast(row)) + .epsilon(1e-4) + .margin(1e-6)); + CHECK(par(row, col, 1) == + Approx(-0.6 + 0.3 * static_cast(col)) + .epsilon(1e-4) + .margin(1e-6)); + CHECK(par(row, col, 2) == + Approx(0.9 + 0.1 * static_cast(row)) + .epsilon(1e-4) + .margin(1e-6)); } } } 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..31a8c237 --- /dev/null +++ b/src/FitMinuit2.cpp @@ -0,0 +1,180 @@ +// SPDX-License-Identifier: MPL-2.0 +// Optional Minuit2 (Migrad + Hesse) fitting backend, compiled only with +// -DAARE_MINUIT2=ON. Kept for side-by-side validation of the built-in +// Levenberg-Marquardt solver in Fit.cpp. +#include "aare/FitMinuit2.hpp" +#include "Chi2.hpp" +#include "FitHelpers.hpp" +#include "Minuit2/FunctionMinimum.h" +#include "Minuit2/MnHesse.h" +#include "Minuit2/MnMigrad.h" +#include "Minuit2/MnStrategy.h" +#include "Minuit2/MnUserParameters.h" +#include "MinuitSteps.hpp" +#include "aare/Models.hpp" +#include "aare/utils/par.hpp" +#include "aare/utils/task.hpp" +#include +#include +#include + +namespace aare::minuit2 { + +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; + compute_ranges(x, y, x_range, y_range, slope_scale); + std::array steps{}; + 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)) + upar.SetLowerLimit(idx, lo); + else if (std::isfinite(hi)) + upar.SetUpperLimit(idx, hi); + if (model.is_user_fixed(idx)) + upar.Fix(idx); + } + return upar; +} + +} // 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 auto result_size = + static_cast(want_errors ? (2 * npar + 1) : (npar + 1)); + NDArray result({result_size}, 0.0); + + const auto start = detail::start_values(model, x, y); + if (!Model::is_valid(std::vector(start.begin(), start.end()))) + return result; + + const auto upar = user_parameters(model, x, y, start); + auto chi2 = (y_err.size() > 0) ? FCN(x, y, y_err) : FCN(x, y); + + ROOT::Minuit2::MnMigrad migrad(chi2, upar, + ROOT::Minuit2::MnStrategy(model.strategy())); + ROOT::Minuit2::FunctionMinimum min = + migrad(model.max_calls(), model.tolerance()); + + if (!min.IsValid()) + return result; + + if (want_errors) { + ROOT::Minuit2::MnHesse hesse; + hesse(chi2, min); + } + + const auto &values = min.UserState().Params(); + for (std::size_t k = 0; k < npar; ++k) + result(static_cast(k)) = values[k]; + + if (want_errors) { + const auto &errors = min.UserState().Errors(); + for (std::size_t k = 0; k < npar; ++k) { + const auto idx = static_cast(k); + result(static_cast(npar + k)) = + model.is_user_fixed(idx) ? 0.0 : errors[k]; + } + result(static_cast(2 * npar)) = min.Fval(); + } else { + result(static_cast(npar)) = min.Fval(); + } + return result; +} + +template +NDArray fit_pixel(const FitModel &model, NDView x, + NDView y) { + return fit_pixel(model, x, y, NDView{}); +} + +template +void fit_3d(const FitModel &model, NDView x, + NDView y, NDView y_err, + NDView par_out, NDView err_out, + NDView chi2_out, int n_threads) { + constexpr std::size_t npar = Model::npar; + detail::check_fit_3d_shapes(x, y, y_err, par_out, err_out, chi2_out); + + const bool has_errors = (y_err.size() > 0); + const bool want_par_errors = (err_out.size() > 0) && model.compute_errors(); + + auto process = [&](ssize_t first_row, ssize_t last_row) { + for (ssize_t row = first_row; row < last_row; row++) { + for (ssize_t col = 0; col < y.shape(1); col++) { + NDView values(&y(row, col, 0), {y.shape(2)}); + NDView errors = + has_errors ? NDView(&y_err(row, col, 0), + {y_err.shape(2)}) + : NDView{}; + + auto res = fit_pixel(model, x, values, errors); + + for (std::size_t k = 0; k < npar; ++k) + par_out(row, col, static_cast(k)) = + res(static_cast(k)); + + if (want_par_errors) { + for (std::size_t k = 0; k < npar; ++k) + err_out(row, col, static_cast(k)) = + res(static_cast(npar + k)); + chi2_out(row, col) = res(static_cast(2 * npar)); + } else { + chi2_out(row, col) = res(static_cast(npar)); + } + } + } + }; + + auto tasks = split_task(0, static_cast(y.shape(0)), n_threads); + RunInParallel(process, tasks); +} + +// NOLINTBEGIN +#define AARE_INSTANTIATE_MINUIT2_FIT(Model) \ + template NDArray fit_pixel( \ + const FitModel &, NDView, NDView, \ + NDView); \ + template NDArray fit_pixel( \ + const FitModel &, NDView, NDView); \ + template void fit_3d(const FitModel &, NDView, \ + NDView, NDView, \ + NDView, NDView, \ + NDView, int); + +AARE_INSTANTIATE_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::minuit2 diff --git a/src/FitMinuit2.test.cpp b/src/FitMinuit2.test.cpp new file mode 100644 index 00000000..616b0a28 --- /dev/null +++ b/src/FitMinuit2.test.cpp @@ -0,0 +1,247 @@ +// SPDX-License-Identifier: MPL-2.0 +// Cross-check of the built-in Levenberg-Marquardt solver against the optional +// Minuit2 backend. Compiled only with -DAARE_MINUIT2=ON. +#include "aare/FitMinuit2.hpp" +#include "aare/Fit.hpp" +#include "aare/FitModel.hpp" +#include "aare/Models.hpp" + +#include +#include +#include +#include +#include +#include +#include + +using aare::NDArray; +using aare::NDView; + +namespace { + +struct Data { + NDArray x; + NDArray y; + NDArray y_err; +}; + +// Noise-free Gaussian, A = 120, mu = 0.8, sigma = 1.3, unit errors. +Data gaussian_data(bool add_noise) { + constexpr ssize_t n = 61; + Data d{NDArray({n}), NDArray({n}), + NDArray({n}, 1.0)}; + for (ssize_t i = 0; i < n; ++i) { + d.x(i) = -6.0 + 0.2 * static_cast(i); + const double z = (d.x(i) - 0.8) / 1.3; + d.y(i) = 120.0 * std::exp(-0.5 * z * z); + if (add_noise) // deterministic pseudo-noise + d.y(i) += 3.0 * std::sin(7.3 * static_cast(i)); + } + return d; +} + +template +Data scurve_data(const std::vector &par, bool add_noise) { + 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; +} + +// Run the same configuration through both backends. The fits must agree on +// parameters and chi2 (relative tolerance tol), unless the built-in solver +// found a clearly lower chi2. Errors are compared for free parameters that +// do not end on a limit, with tolerance err_tol: Gauss-Newton errors match +// Hesse at a good minimum but differ when a parameter is forced to a wrong +// value, so misfit cases pass a looser err_tol. +template +void compare_backends(const std::string &name, Data d, bool weighted, + bool errors, unsigned 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; + aare::FitModel model(0, max_calls, 0.5, errors); + setup(model); + NDView ev = weighted ? d.y_err.view() : NDView{}; + const auto lm = weighted + ? aare::fit_pixel(model, d.x.view(), d.y.view(), ev) + : aare::fit_pixel(model, d.x.view(), d.y.view()); + const auto mn = + weighted ? aare::minuit2::fit_pixel(model, d.x.view(), d.y.view(), ev) + : aare::minuit2::fit_pixel(model, d.x.view(), d.y.view()); + INFO(name); + REQUIRE(lm.size() == mn.size()); + + const auto ichi = static_cast(errors ? 2 * npar : npar); + const double chi_lm = lm(ichi); + const double chi_mn = mn(ichi); + if (chi_lm > 0.0 && chi_mn > 0.0 && + chi_mn - chi_lm > tol * std::max(chi_mn, 1.0)) { + WARN(name + ": built-in solver found a lower chi2 (" + + std::to_string(chi_lm) + " vs " + std::to_string(chi_mn) + ")"); + return; + } + double scale = 0.0; + for (std::size_t k = 0; k < npar; ++k) + scale = std::max(scale, std::abs(mn(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(mn(i)), 1e-2 * scale, 1e-12}); + INFO("parameter " << model.GetParName(static_cast(k))); + CHECK(std::abs(lm(i) - mn(i)) / denom <= tol); + } + CHECK(std::abs(chi_lm - chi_mn) / std::max(std::abs(chi_mn), 1.0) <= tol); + if (!errors) + return; + for (std::size_t k = 0; k < npar; ++k) { + const auto idx = static_cast(k); + const auto i = static_cast(k); + const double v = mn(i); + const double tolb = 1e-6 * std::max(1.0, std::abs(v)); + const bool on_limit = std::abs(v - model.lower_limit(idx)) <= tolb || + std::abs(v - model.upper_limit(idx)) <= tolb; + if (model.is_user_fixed(idx) || on_limit) + continue; + const auto ie = static_cast(npar + k); + INFO("error of " << model.GetParName(idx)); + CHECK(std::abs(lm(ie) - mn(ie)) / std::max(std::abs(mn(ie)), 1e-12) <= + err_tol); + } +} + +} // namespace + +TEST_CASE("Built-in solver matches the Minuit2 backend on Gaussian fits", + "[fit][minuit2]") { + using aare::model::Gaussian; + const auto none = [](auto &) {}; + compare_backends("exact, unweighted", gaussian_data(false), false, + false, 100, 1e-4, none); + compare_backends("exact, weighted with errors", + gaussian_data(false), true, true, 100, 1e-4, + none); + compare_backends( + "sigma fixed at a wrong value", gaussian_data(false), true, true, 100, + 1e-3, [](auto &m) { m.FixParameter("sigma", 1.5); }, 0.2); + compare_backends("sigma fixed at the true value", + gaussian_data(false), true, true, 100, 1e-3, + [](auto &m) { m.FixParameter("sigma", 1.3); }); + compare_backends( + "active upper limit on mu", gaussian_data(false), false, false, 100, + 1e-3, [](auto &m) { m.SetParLimits("mu", -5.0, 0.5); }); + compare_backends( + "inactive limits", gaussian_data(false), false, false, 100, 1e-4, + [](auto &m) { m.SetParLimits("mu", -5.0, 5.0); }); + compare_backends( + "active upper limit on sigma", gaussian_data(false), true, true, 100, + 1e-3, [](auto &m) { m.SetParLimits("sigma", 0.5, 1.0); }, 0.2); + compare_backends("user start value", gaussian_data(false), false, + false, 100, 1e-4, + [](auto &m) { m.SetParameter("mu", 0.0); }); + compare_backends( + "fixed amplitude and active limit on mu", gaussian_data(false), true, + true, 100, 1e-3, + [](auto &m) { + m.FixParameter("A", 100.0); + m.SetParLimits("mu", 0.6, 0.7); + }, + 0.2); + compare_backends("max_calls exhausted", gaussian_data(false), + false, false, 3, 1e-9, none); + compare_backends("noisy, weighted with errors", + gaussian_data(true), true, true, 500, 1e-3, + none); + compare_backends( + "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("Built-in solver matches the Minuit2 backend on S-curve fits", + "[fit][minuit2]") { + using aare::model::RisingScurve; + const std::vector truth{10.0, 0.1, 50.0, 3.0, 1000.0, 2.0}; + const auto none = [](auto &) {}; + compare_backends("noisy, weighted with errors", + scurve_data(truth, true), true, + true, 500, 2e-2, none); + compare_backends("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); + }); + compare_backends( + "noisy, active limit on the threshold", + scurve_data(truth, true), false, false, 500, 1e-2, + [](auto &m) { m.SetParLimits("mu", 40.0, 48.0); }); + compare_backends( + "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("Built-in solver matches the Minuit2 backend on a data cube", + "[fit][minuit2]") { + constexpr ssize_t rows = 3; + constexpr ssize_t cols = 2; + constexpr ssize_t n = 61; + NDArray x({n}); + NDArray y({rows, cols, n}); + NDArray y_err({rows, cols, n}, 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; ++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)); + } + } + } + NDArray par_lm({rows, cols, 3}); + NDArray err_lm({rows, cols, 3}); + NDArray chi2_lm({rows, cols}); + NDArray par_mn({rows, cols, 3}); + NDArray err_mn({rows, cols, 3}); + NDArray chi2_mn({rows, cols}); + const aare::FitModel model(0, 500, 0.5, true); + + aare::fit_3d(model, x.view(), y.view(), y_err.view(), par_lm.view(), + err_lm.view(), chi2_lm.view(), 2); + aare::minuit2::fit_3d(model, x.view(), y.view(), y_err.view(), + par_mn.view(), err_mn.view(), chi2_mn.view(), 2); + + 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(par_lm(row, col, k) == + Catch::Approx(par_mn(row, col, k)).epsilon(1e-3)); + CHECK(err_lm(row, col, k) == + Catch::Approx(err_mn(row, col, k)).epsilon(2e-2)); + } + CHECK(chi2_lm(row, col) == + Catch::Approx(chi2_mn(row, col)).epsilon(1e-4)); + } + } +} diff --git a/src/LevenbergMarquardt.hpp b/src/LevenbergMarquardt.hpp new file mode 100644 index 00000000..99acdb89 --- /dev/null +++ b/src/LevenbergMarquardt.hpp @@ -0,0 +1,416 @@ +// 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. + * - 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 = static_cast(model.max_calls()); + 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