detectors/aare
25
0 Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-11-06 09:43:26 +01:00
Code Issues Actions Packages Releases Activity

first attempt on custom minimizer

Browse Source
This commit is contained in:
Erik Frojdh committed 2026-09-24 17:36:30 +02:00
1 parent 4677e1003e
commit f001188dcf
28 files changed
+2078 -533
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

AGENTS.md
+16 -5
View File
@@ -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<Model>`, 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
CMakeLists.txt
+92 -47
View File
@@ -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 $<BUILD_INTERFACE:aare::Minuit2>)
target_include_directories(
aare_core SYSTEM
PRIVATE $<TARGET_PROPERTY:aare::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>)
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 $<BUILD_INTERFACE:aare::Minuit2>)
target_include_directories(
aare_core SYSTEM
PRIVATE $<TARGET_PROPERTY:aare::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>)
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"
README.md
+5 -4
View File
@@ -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
RELEASE.md
+26 -7
View File
@@ -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
THIRD-PARTY-NOTICES.txt
+5 -34
View File
@@ -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}
benchmarks/CMakeLists.txt
-3
View File
@@ -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 $<TARGET_PROPERTY:Minuit2::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>)
set_target_properties(fit_benchmark PROPERTIES RUNTIME_OUTPUT_DIRECTORY
${CMAKE_BINARY_DIR})
benchmarks/fit_benchmark.cpp
+43 -14
View File
@@ -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 <benchmark/benchmark.h>
#include <cmath>
@@ -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 <typename... Args>
static aare::NDArray<double, 1> fit(Args &&...args) {
return aare::fit_pixel<aare::model::Gaussian>(
std::forward<Args>(args)...);
}
};
#ifdef AARE_MINUIT2
struct Minuit2 {
template <typename... Args>
static aare::NDArray<double, 1> fit(Args &&...args) {
return aare::minuit2::fit_pixel<aare::model::Gaussian>(
std::forward<Args>(args)...);
}
};
#endif
// ----------
// Benchmarks
// ----------
// Minuit2, analytic gradient (no Hesse)
static void BM_FitGausMinuitGrad(benchmark::State &state) {
// Gaussian, unweighted, no errors
template <typename Backend>
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<double, 1> result;
for (auto _ : state) {
result = aare::fit_pixel<aare::model::Gaussian>(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 <typename Backend>
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<aare::model::Gaussian>(
0, 500, 0.5, true); // compute_errors = true -> Runs Hesse and provides
// errors on fitted params
const auto model = aare::FitModel<aare::model::Gaussian>(0, 500, 0.5, true);
aare::NDArray<double, 1> result;
for (auto _ : state) {
result = aare::fit_pixel<aare::model::Gaussian>(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();
conda-recipe/meta.yaml
+1 -2
View File
@@ -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
docs/src/Installation.rst
+9
View File
@@ -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.
docs/src/Requirements.rst
+3 -1
View File
@@ -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
docs/src/pyFit.rst
+29 -2
View File
@@ -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
include/aare/Fit.hpp
+3 -2
View File
@@ -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).
include/aare/FitMinuit2.hpp
+33
View File
@@ -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 <typename Model>
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y, NDView<double, 1> y_err);
template <typename Model>
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y);
template <typename Model>
void fit_3d(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> err_out,
NDView<double, 2> chi2_out, int n_threads);
} // namespace aare::minuit2
include/aare/FitModel.hpp
+46 -34
View File
@@ -3,59 +3,54 @@
#include "aare/Models.hpp"
#include <array>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
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 <typename Model> 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<FitModelImpl> impl_;
unsigned int max_calls_;
double tolerance_;
bool compute_errors_;
std::array<bool, Model::npar> user_fixed_{};
std::array<bool, Model::npar> 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 <typename Model> 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<double, npar> lower_{};
std::array<double, npar> upper_{};
std::array<double, npar> value_{};
std::array<bool, npar> user_fixed_{};
std::array<bool, npar> user_start_{};
};
// Suppress implicit instantiation for all supported model types.
include/aare/Models.hpp
+1 -103
View File
@@ -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<double>::infinity();
/**
* @brief Compute data-range statistics used by step-size estimators.
*
* Model-independent, called once per pixel.
*/
inline void compute_ranges(NDView<double, 1> x, NDView<double, 1> 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<double, npar> &start,
[[maybe_unused]] double x_range, double y_range,
double slope_scale,
std::array<double, npar> &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<double, npar> &start,
double x_range, double y_range,
double slope_scale,
std::array<double, npar> &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<double, npar> &start,
double x_range, double y_range,
double /*slope_scale*/,
std::array<double, npar> &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<double, npar> &start,
double x_range, double y_range,
double /*slope_scale*/,
std::array<double, npar> &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<double, npar> &start,
double x_range, double y_range,
double slope_scale,
std::array<double, npar> &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<double, npar> &start,
double x_range, double y_range,
double slope_scale,
std::array<double, npar> &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<double, npar> &start,
double x_range, double y_range,
double slope_scale,
std::array<double, npar> &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<double, npar> &start,
double x_range, double y_range,
double slope_scale,
std::array<double, npar> &steps) {
RisingScurve::compute_steps(start, x_range, y_range, slope_scale,
steps);
}
};
} // namespace aare::model
pyproject.toml
+1 -1
View File
@@ -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"]
python/CMakeLists.txt
-4
View File
@@ -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 $<TARGET_PROPERTY:aare::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>)
# List of python files to be copied to the build directory
set(PYTHON_FILES
aare/__init__.py
python/aare/experimental.py
+7
View File
@@ -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
python/src/fit.hpp
+133 -98
View File
@@ -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 <typename Model>
// 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 <typename Model, typename... Args>
static aare::NDArray<double, 1> fit_pixel(Args &&...args) {
return aare::fit_pixel<Model>(std::forward<Args>(args)...);
}
template <typename Model, typename... Args>
static void fit_3d(Args &&...args) {
aare::fit_3d<Model>(std::forward<Args>(args)...);
}
};
#ifdef AARE_MINUIT2
struct Minuit2Backend {
template <typename Model, typename... Args>
static aare::NDArray<double, 1> fit_pixel(Args &&...args) {
return aare::minuit2::fit_pixel<Model>(std::forward<Args>(args)...);
}
template <typename Model, typename... Args>
static void fit_3d(Args &&...args) {
aare::minuit2::fit_3d<Model>(std::forward<Args>(args)...);
}
};
#endif
template <typename Model, typename Backend>
py::object
fit_dispatch(const aare::FitModel<Model> &model,
py::array_t<double, py::array::c_style | py::array::forcecast> x,
@@ -24,7 +54,23 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
py::class_<FM>(m, name)
.def(py::init<unsigned int, unsigned int, double, bool>(),
py::arg("strategy") = 0, py::arg("max_calls") = 100,
py::arg("tolerance") = 0.5, py::arg("compute_errors") = false)
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<unsigned int, double, double>(&FM::SetParLimits),
py::arg("idx"), py::arg("lo"), py::arg("hi"))
@@ -83,10 +129,12 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::object y_err_obj, int n_threads) -> py::object {
return fit_dispatch<Model>(self, x, y, y_err_obj, n_threads);
return fit_dispatch<Model, LmBackend>(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<double, 1> &result,
}
// Helper: typed dispatch for one Model, handles 1D/3D + y_err logic
template <typename Model>
template <typename Model, typename Backend>
py::object
fit_dispatch(const aare::FitModel<Model> &model,
py::array_t<double, py::array::c_style | py::array::forcecast> x,
@@ -172,9 +220,9 @@ fit_dispatch(const aare::FitModel<Model> &model,
new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
auto y_view_err = make_view_3d(y_err);
aare::fit_3d<Model>(model, x_view, y_view, y_view_err,
par_out->view(), err_out->view(),
chi2_out->view(), n_threads);
Backend::template fit_3d<Model>(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> &model,
NDView<double, 3> dummy_err{};
NDView<double, 3> dummy_err_out{};
aare::fit_3d<Model>(model, x_view, y_view, dummy_err,
par_out->view(), dummy_err_out,
chi2_out->view(), n_threads);
Backend::template fit_3d<Model>(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> &model,
}
auto y_view_err = make_view_1d(y_err);
result = aare::fit_pixel<Model>(model, x_view, y_view, y_view_err);
result = Backend::template fit_pixel<Model>(model, x_view, y_view,
y_view_err);
} else {
result = aare::fit_pixel<Model>(model, x_view, y_view);
result = Backend::template fit_pixel<Model>(model, x_view, y_view);
}
return pack_1d_result_dict<Model>(result, model.compute_errors());
@@ -226,7 +275,39 @@ fit_dispatch(const aare::FitModel<Model> &model,
}
}
void define_fit_bindings(py::module &m) {
// Resolve the model type behind a Python object and dispatch to Backend.
template <typename Backend>
py::object dispatch_any_model(
py::object model_obj,
py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::object y_err_obj, int n_threads) {
using namespace aare::model;
#define AARE_DISPATCH_MODEL(Model) \
if (py::isinstance<aare::FitModel<Model>>(model_obj)) { \
const auto &mdl = model_obj.cast<const aare::FitModel<Model> &>(); \
return fit_dispatch<Model, Backend>(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<aare::model::Gaussian>(m, "Gaussian");
bind_fit_model<aare::model::GaussianErfcPlateau>(m, "GaussianErfcPlateau");
@@ -245,86 +326,20 @@ void define_fit_bindings(py::module &m) {
py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::object y_err_obj, int n_threads) -> py::object {
using namespace aare::model;
// ── Polynomial of degree 1 ───────
if (py::isinstance<aare::FitModel<Pol1>>(model_obj)) {
const auto &mdl =
model_obj.cast<const aare::FitModel<Pol1> &>();
return fit_dispatch<Pol1>(mdl, x, y, y_err_obj, n_threads);
}
// ── Polynomial of degree 2 ───────
if (py::isinstance<aare::FitModel<Pol2>>(model_obj)) {
const auto &mdl =
model_obj.cast<const aare::FitModel<Pol2> &>();
return fit_dispatch<Pol2>(mdl, x, y, y_err_obj, n_threads);
}
// ── Gaussian ───────
if (py::isinstance<aare::FitModel<Gaussian>>(model_obj)) {
const auto &mdl =
model_obj.cast<const aare::FitModel<Gaussian> &>();
return fit_dispatch<Gaussian>(mdl, x, y, y_err_obj, n_threads);
}
// ── GaussianErfcPlateau ───────
if (py::isinstance<aare::FitModel<GaussianErfcPlateau>>(
model_obj)) {
const auto &mdl =
model_obj
.cast<const aare::FitModel<GaussianErfcPlateau> &>();
return fit_dispatch<GaussianErfcPlateau>(mdl, x, y, y_err_obj,
n_threads);
}
// ── GaussianChargeSharing ───────
if (py::isinstance<aare::FitModel<GaussianChargeSharing>>(
model_obj)) {
const auto &mdl =
model_obj
.cast<const aare::FitModel<GaussianChargeSharing> &>();
return fit_dispatch<GaussianChargeSharing>(mdl, x, y, y_err_obj,
n_threads);
}
// ── GaussianChargeSharingKb ───────
if (py::isinstance<aare::FitModel<GaussianChargeSharingKb>>(
model_obj)) {
const auto &mdl = model_obj.cast<
const aare::FitModel<GaussianChargeSharingKb> &>();
return fit_dispatch<GaussianChargeSharingKb>(
mdl, x, y, y_err_obj, n_threads);
}
// ── Rising Scurve ───────
if (py::isinstance<aare::FitModel<RisingScurve>>(model_obj)) {
const auto &mdl =
model_obj.cast<const aare::FitModel<RisingScurve> &>();
return fit_dispatch<RisingScurve>(mdl, x, y, y_err_obj,
n_threads);
}
// ── Falling Scurve ───────
if (py::isinstance<aare::FitModel<FallingScurve>>(model_obj)) {
const auto &mdl =
model_obj.cast<const aare::FitModel<FallingScurve> &>();
return fit_dispatch<FallingScurve>(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<LmBackend>(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<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::object y_err_obj, int n_threads) -> py::object {
return dispatch_any_model<Minuit2Backend>(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
}
python/src/module.cpp
+1 -1
View File
@@ -84,7 +84,7 @@ PYBIND11_MODULE(_aare, m) {
define_fast_pedestal_bindings<double>(m, "FastPedestal_d");
define_fast_pedestal_bindings<float>(m, "FastPedestal_f");
define_fast_pedestal_bindings<int16_t>(m, "FastPedestal_i16");
define_fit_bindings(m);
define_fit_bindings(m, experimental);
define_interpolation_bindings(m);
define_jungfrau_data_file_io_bindings(m);
python/tests/test_Fit.py
+169
View File
@@ -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)
src/Fit.cpp
+54 -143
View File
@@ -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 <array>
#include <cmath>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <string>
namespace aare {
// ============================================================================
// FitModel<Model> — method definitions
// (constructor, destructor, copy, and all methods that touch Minuit2 state)
// ============================================================================
template <typename Model> struct FitModel<Model>::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 <typename Model>
FitModel<Model>::FitModel(unsigned int strategy, unsigned int max_calls,
double tolerance, bool compute_errors)
: impl_(std::make_unique<FitModelImpl>(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 <typename Model> FitModel<Model>::~FitModel() = default;
template <typename Model>
FitModel<Model>::FitModel(const FitModel &other)
: impl_(std::make_unique<FitModelImpl>(*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 <typename Model>
FitModel<Model> &FitModel<Model>::operator=(const FitModel &other) {
if (this != &other) {
impl_ = std::make_unique<FitModelImpl>(*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<Model>::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 <typename Model>
unsigned int FitModel<Model>::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<unsigned int>(i);
}
throw std::runtime_error("FitModel: unknown parameter name '" + name + "'");
@@ -84,19 +43,24 @@ unsigned int FitModel<Model>::checked_index(const std::string &name) const {
template <typename Model>
void FitModel<Model>::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 <typename Model>
void FitModel<Model>::FixParameter(unsigned int idx, double val) {
SetParameter(idx, val);
impl_->upar.Fix(idx);
user_fixed_[idx] = true;
}
template <typename Model>
void FitModel<Model>::ReleaseParameter(unsigned int idx) {
impl_->upar.Release(idx);
check_index(idx);
user_fixed_[idx] = false;
}
@@ -107,7 +71,8 @@ void FitModel<Model>::ReleaseParameter(const std::string &name) {
template <typename Model>
void FitModel<Model>::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<Model>::SetParLimits(const std::string &name, double lo,
template <typename Model>
std::string FitModel<Model>::GetParName(unsigned int idx) const {
return impl_->upar.GetName(idx);
check_index(idx);
return Model::param_info[idx].name;
}
template <typename Model>
std::vector<std::string> FitModel<Model>::GetParNames() const {
std::vector<std::string> 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 <typename Model>
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y, NDView<double, 1> y_err) {
using FCN = func::Chi2Model1DGrad<Model>;
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<ssize_t>(want_errors ? (2 * npar + 1) : (npar + 1));
NDArray<double, 1> result({result_size}, 0.0);
auto start = Model::estimate_par(x, y);
if (!Model::is_valid(std::vector<double>(start.begin(), start.end()))) {
return NDArray<double, 1>({result_size}, 0.0);
}
double x_range, y_range, slope_scale;
model::compute_ranges(x, y, x_range, y_range, slope_scale);
std::array<double, npar> 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<double, 1>({result_size}, 0.0);
const auto start = detail::start_values(model, x, y);
detail::LevenbergMarquardt<Model> solver;
std::array<double, npar> 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<double, 1> 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<ssize_t>(npar + k)) = err[k];
result(static_cast<ssize_t>(2 * npar)) = res.chi2;
} else {
result(static_cast<ssize_t>(npar)) = res.chi2;
}
const auto &values = min.UserState().Params();
NDArray<double, 1> 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> &model, NDView<double, 1> x,
NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> err_out,
NDView<double, 2> 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<Model>(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<Model> solver; // buffers reused per task
std::array<double, npar> par{};
std::array<double, npar> err{};
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> errors =
has_errors ? NDView<double, 1>(&y_err(row, col, 0),
{y_err.shape(2)})
: NDView<double, 1>{};
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<ssize_t>(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<ssize_t>(k)) = err[k];
}
chi2_out(row, col) = res.valid ? res.chi2 : 0.0;
}
}
};
src/Fit.test.cpp
+359 -28
View File
@@ -6,12 +6,26 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
using aare::NDArray;
using aare::NDView;
using Catch::Approx;
namespace {
constexpr ssize_t n_points = 61;
void fill_gaussian(aare::NDArray<double, 1> &x, aare::NDArray<double, 1> &y,
struct Data {
NDArray<double, 1> x;
NDArray<double, 1> y;
NDArray<double, 1> y_err;
};
void fill_gaussian(NDArray<double, 1> &x, NDArray<double, 1> &y,
double amplitude, double mean, double sigma) {
for (ssize_t i = 0; i < n_points; ++i) {
x(i) = -6.0 + 0.2 * static_cast<double>(i);
@@ -20,42 +34,301 @@ void fill_gaussian(aare::NDArray<double, 1> &x, aare::NDArray<double, 1> &y,
}
}
// Noise-free Gaussian with A = 120, mu = 0.8, sigma = 1.3 and unit errors.
Data gaussian_data() {
Data d{NDArray<double, 1>({n_points}), NDArray<double, 1>({n_points}),
NDArray<double, 1>({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 <typename Model> Data scurve_data(const std::vector<double> &par) {
constexpr ssize_t n = 100;
Data d{NDArray<double, 1>({n}), NDArray<double, 1>({n}),
NDArray<double, 1>({n}, 1.0)};
for (ssize_t i = 0; i < n; ++i) {
d.x(i) = 100.0 * static_cast<double>(i) / static_cast<double>(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<double, 1> x({n_points});
aare::NDArray<double, 1> y({n_points});
aare::NDArray<double, 1> 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<aare::model::Gaussian> 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<aare::model::Gaussian> 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<double, 1> x({n_points});
aare::NDArray<double, 1> values({n_points});
TEST_CASE("Fit weighted Pol1 errors match the analytic formula", "[fit]") {
constexpr ssize_t n = 20;
NDArray<double, 1> x({n});
NDArray<double, 1> y({n});
NDArray<double, 1> 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<double>(i);
y(i) = 2.0 + 0.5 * x(i);
y_err(i) = 1.0 + 0.1 * static_cast<double>(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<aare::model::Pol1> 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<aare::model::Gaussian> 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<aare::model::Gaussian> 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<aare::model::Gaussian> 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<aare::model::Gaussian> 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<aare::model::Gaussian> 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<aare::model::Gaussian> free_model;
const auto reference =
aare::fit_pixel(free_model, d.x.view(), d.y.view());
aare::FitModel<aare::model::Gaussian> model;
model.SetParLimits("mu", -5.0, 5.0);
model.SetParLimits("sigma", 0.5,
std::numeric_limits<double>::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<aare::model::Gaussian> 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<double, 1>({n_points}), NDArray<double, 1>({n_points}),
NDArray<double, 1>({n_points}, 1.0)};
fill_gaussian(n.x, n.y, -120.0, 0.8, 1.3);
aare::FitModel<aare::model::Gaussian> 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<aare::model::Gaussian> 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<aare::model::Gaussian> 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<aare::model::Gaussian> model;
CHECK(model.GetParNames() == std::vector<std::string>{"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<double> truth{10.0, 0.1, 50.0, 3.0, 1000.0, 2.0};
SECTION("rising, free fit with errors") {
auto d = scurve_data<aare::model::RisingScurve>(truth);
const aare::FitModel<aare::model::RisingScurve> 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<ssize_t>(k)) ==
Approx(truth[k]).epsilon(1e-3));
CHECK(res(static_cast<ssize_t>(6 + k)) > 0.0);
}
CHECK(res(12) == Approx(0.0).margin(1e-3));
}
SECTION("rising, fixed slope and limited width") {
auto d = scurve_data<aare::model::RisingScurve>(truth);
aare::FitModel<aare::model::RisingScurve> 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<aare::model::RisingScurve>(truth);
aare::FitModel<aare::model::RisingScurve> 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<aare::model::FallingScurve>(truth);
const aare::FitModel<aare::model::FallingScurve> 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<ssize_t>(k)) ==
Approx(truth[k]).epsilon(1e-3));
}
}
TEST_CASE("Fit a Gaussian data cube in parallel", "[fit]") {
NDArray<double, 1> x({n_points});
NDArray<double, 1> values({n_points});
fill_gaussian(x, values, 80.0, -0.6, 0.9);
aare::NDArray<double, 3> y({2, 2, n_points});
NDArray<double, 3> 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<double, 3> par({2, 2, 3});
aare::NDArray<double, 2> chi2({2, 2});
NDArray<double, 3> par({2, 2, 3});
NDArray<double, 2> chi2({2, 2});
const aare::FitModel<aare::model::Gaussian> model;
aare::fit_3d(model, x.view(), y.view(), aare::NDView<double, 3>{},
par.view(), aare::NDView<double, 3>{}, chi2.view(), 2);
aare::fit_3d(model, x.view(), y.view(), NDView<double, 3>{}, par.view(),
NDView<double, 3>{}, 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<double, 1>(x.data(), {10}),
y.view(), NDView<double, 3>{}, par.view(),
NDView<double, 3>{}, 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<double, 1> x({n_points});
NDArray<double, 3> y({rows, cols, n_points});
NDArray<double, 3> y_err({rows, cols, n_points}, 1.0);
NDArray<double, 1> 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<double>(row),
-0.6 + 0.3 * static_cast<double>(col),
0.9 + 0.1 * static_cast<double>(row));
for (ssize_t i = 0; i < n_points; ++i)
y(row, col, i) = tmp(i);
}
}
NDArray<double, 3> par({rows, cols, 3});
NDArray<double, 3> err({rows, cols, 3});
NDArray<double, 2> chi2({rows, cols});
const aare::FitModel<aare::model::Gaussian> 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<double, 1> yv(&y(row, col, 0), {n_points});
NDView<double, 1> 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<double>(row))
.epsilon(1e-4)
.margin(1e-6));
CHECK(par(row, col, 1) ==
Approx(-0.6 + 0.3 * static_cast<double>(col))
.epsilon(1e-4)
.margin(1e-6));
CHECK(par(row, col, 2) ==
Approx(0.9 + 0.1 * static_cast<double>(row))
.epsilon(1e-4)
.margin(1e-6));
}
}
}
src/FitHelpers.hpp
+73
View File
@@ -0,0 +1,73 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/FitModel.hpp"
#include "aare/NDView.hpp"
#include <algorithm>
#include <array>
#include <stdexcept>
// 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 <typename Model>
std::array<double, Model::npar> start_values(const FitModel<Model> &model,
NDView<double, 1> x,
NDView<double, 1> y) {
auto start = Model::estimate_par(x, y);
for (std::size_t k = 0; k < Model::npar; ++k) {
const auto idx = static_cast<unsigned int>(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 <typename Model>
void check_fit_3d_shapes(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> err_out,
NDView<double, 2> chi2_out) {
constexpr auto npar = static_cast<ssize_t>(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
src/FitMinuit2.cpp
+180
View File
@@ -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 <array>
#include <cmath>
#include <vector>
namespace aare::minuit2 {
namespace {
// Translate the FitModel state into Minuit2 user parameters for one pixel.
template <typename Model>
ROOT::Minuit2::MnUserParameters
user_parameters(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y,
const std::array<double, Model::npar> &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<double, npar> steps{};
Steps<Model>::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<unsigned int>(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 <typename Model>
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y, NDView<double, 1> y_err) {
using FCN = func::Chi2Model1DGrad<Model>;
constexpr std::size_t npar = Model::npar;
const bool want_errors = model.compute_errors();
const auto result_size =
static_cast<ssize_t>(want_errors ? (2 * npar + 1) : (npar + 1));
NDArray<double, 1> result({result_size}, 0.0);
const auto start = detail::start_values(model, x, y);
if (!Model::is_valid(std::vector<double>(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<ssize_t>(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<unsigned int>(k);
result(static_cast<ssize_t>(npar + k)) =
model.is_user_fixed(idx) ? 0.0 : errors[k];
}
result(static_cast<ssize_t>(2 * npar)) = min.Fval();
} else {
result(static_cast<ssize_t>(npar)) = min.Fval();
}
return result;
}
template <typename Model>
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y) {
return fit_pixel(model, x, y, NDView<double, 1>{});
}
template <typename Model>
void fit_3d(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> err_out,
NDView<double, 2> chi2_out, int n_threads) {
constexpr std::size_t npar = Model::npar;
detail::check_fit_3d_shapes<Model>(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<double, 1> values(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> errors =
has_errors ? NDView<double, 1>(&y_err(row, col, 0),
{y_err.shape(2)})
: NDView<double, 1>{};
auto res = fit_pixel(model, x, values, errors);
for (std::size_t k = 0; k < npar; ++k)
par_out(row, col, static_cast<ssize_t>(k)) =
res(static_cast<ssize_t>(k));
if (want_par_errors) {
for (std::size_t k = 0; k < npar; ++k)
err_out(row, col, static_cast<ssize_t>(k)) =
res(static_cast<ssize_t>(npar + k));
chi2_out(row, col) = res(static_cast<ssize_t>(2 * npar));
} else {
chi2_out(row, col) = res(static_cast<ssize_t>(npar));
}
}
}
};
auto tasks = split_task(0, static_cast<int>(y.shape(0)), n_threads);
RunInParallel(process, tasks);
}
// NOLINTBEGIN
#define AARE_INSTANTIATE_MINUIT2_FIT(Model) \
template NDArray<double, 1> fit_pixel<Model>( \
const FitModel<Model> &, NDView<double, 1>, NDView<double, 1>, \
NDView<double, 1>); \
template NDArray<double, 1> fit_pixel<Model>( \
const FitModel<Model> &, NDView<double, 1>, NDView<double, 1>); \
template void fit_3d<Model>(const FitModel<Model> &, NDView<double, 1>, \
NDView<double, 3>, NDView<double, 3>, \
NDView<double, 3>, NDView<double, 3>, \
NDView<double, 2>, 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
src/FitMinuit2.test.cpp
+247
View File
@@ -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 <algorithm>
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <cmath>
#include <functional>
#include <string>
#include <vector>
using aare::NDArray;
using aare::NDView;
namespace {
struct Data {
NDArray<double, 1> x;
NDArray<double, 1> y;
NDArray<double, 1> 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<double, 1>({n}), NDArray<double, 1>({n}),
NDArray<double, 1>({n}, 1.0)};
for (ssize_t i = 0; i < n; ++i) {
d.x(i) = -6.0 + 0.2 * static_cast<double>(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<double>(i));
}
return d;
}
template <typename Model>
Data scurve_data(const std::vector<double> &par, bool add_noise) {
constexpr ssize_t n = 100;
Data d{NDArray<double, 1>({n}), NDArray<double, 1>({n}),
NDArray<double, 1>({n}, 1.0)};
for (ssize_t i = 0; i < n; ++i) {
d.x(i) = 100.0 * static_cast<double>(i) / static_cast<double>(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<double>(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 <typename Model, typename Setup>
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> model(0, max_calls, 0.5, errors);
setup(model);
NDView<double, 1> ev = weighted ? d.y_err.view() : NDView<double, 1>{};
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<ssize_t>(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<ssize_t>(k))));
for (std::size_t k = 0; k < npar; ++k) {
const auto i = static_cast<ssize_t>(k);
const double denom = std::max({std::abs(mn(i)), 1e-2 * scale, 1e-12});
INFO("parameter " << model.GetParName(static_cast<unsigned int>(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<unsigned int>(k);
const auto i = static_cast<ssize_t>(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<ssize_t>(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<Gaussian>("exact, unweighted", gaussian_data(false), false,
false, 100, 1e-4, none);
compare_backends<Gaussian>("exact, weighted with errors",
gaussian_data(false), true, true, 100, 1e-4,
none);
compare_backends<Gaussian>(
"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<Gaussian>("sigma fixed at the true value",
gaussian_data(false), true, true, 100, 1e-3,
[](auto &m) { m.FixParameter("sigma", 1.3); });
compare_backends<Gaussian>(
"active upper limit on mu", gaussian_data(false), false, false, 100,
1e-3, [](auto &m) { m.SetParLimits("mu", -5.0, 0.5); });
compare_backends<Gaussian>(
"inactive limits", gaussian_data(false), false, false, 100, 1e-4,
[](auto &m) { m.SetParLimits("mu", -5.0, 5.0); });
compare_backends<Gaussian>(
"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<Gaussian>("user start value", gaussian_data(false), false,
false, 100, 1e-4,
[](auto &m) { m.SetParameter("mu", 0.0); });
compare_backends<Gaussian>(
"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<Gaussian>("max_calls exhausted", gaussian_data(false),
false, false, 3, 1e-9, none);
compare_backends<Gaussian>("noisy, weighted with errors",
gaussian_data(true), true, true, 500, 1e-3,
none);
compare_backends<Gaussian>(
"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<double> truth{10.0, 0.1, 50.0, 3.0, 1000.0, 2.0};
const auto none = [](auto &) {};
compare_backends<RisingScurve>("noisy, weighted with errors",
scurve_data<RisingScurve>(truth, true), true,
true, 500, 2e-2, none);
compare_backends<RisingScurve>("noisy, fixed slope and limited width",
scurve_data<RisingScurve>(truth, true), true,
true, 500, 1e-2, [](auto &m) {
m.FixParameter("p1", 0.1);
m.SetParLimits("sigma", 1.0, 4.0);
});
compare_backends<RisingScurve>(
"noisy, active limit on the threshold",
scurve_data<RisingScurve>(truth, true), false, false, 500, 1e-2,
[](auto &m) { m.SetParLimits("mu", 40.0, 48.0); });
compare_backends<RisingScurve>(
"noisy, two fixed, limited amplitude, start value",
scurve_data<RisingScurve>(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<double, 1> x({n});
NDArray<double, 3> y({rows, cols, n});
NDArray<double, 3> 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<double>(col);
const double sigma = 0.9 + 0.1 * static_cast<double>(row);
for (ssize_t i = 0; i < n; ++i) {
x(i) = -6.0 + 0.2 * static_cast<double>(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<double>(i));
}
}
}
NDArray<double, 3> par_lm({rows, cols, 3});
NDArray<double, 3> err_lm({rows, cols, 3});
NDArray<double, 2> chi2_lm({rows, cols});
NDArray<double, 3> par_mn({rows, cols, 3});
NDArray<double, 3> err_mn({rows, cols, 3});
NDArray<double, 2> chi2_mn({rows, cols});
const aare::FitModel<aare::model::Gaussian> 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));
}
}
}
src/LevenbergMarquardt.hpp
+416
View File
@@ -0,0 +1,416 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/FitModel.hpp"
#include "aare/NDView.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
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 <typename Model> class LevenbergMarquardt {
public:
static constexpr int npar = static_cast<int>(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> &model, NDView<double, 1> x,
NDView<double, 1> y, NDView<double, 1> s,
const std::array<double, Model::npar> &start, double *par_out,
double *err_out) {
x_ = x;
y_ = y;
s_ = s;
weighted_ = s.size() > 0;
const auto n = static_cast<std::size_t>(x.size());
nfree_ = 0;
for (int k = 0; k < npar; ++k) {
const auto idx = static_cast<unsigned int>(k);
lower_[k] = model.lower_limit(idx);
upper_[k] = model.upper_limit(idx);
p_[k] = start[static_cast<std::size_t>(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<std::size_t>(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<int>(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<double, Model::npar> &q,
std::vector<double> &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<double, Model::npar> 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<std::size_t>(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<std::size_t>(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<double>::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<double, 1> x_, y_, s_;
bool weighted_ = false;
int nfree_ = 0;
int free_[npar] = {};
bool active_[npar] = {};
double lower_[npar] = {};
double upper_[npar] = {};
std::array<double, Model::npar> p_{};
std::array<double, Model::npar> q_{};
std::array<double, Model::npar> q_new_{};
std::array<double, Model::npar> 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<double> r_;
std::vector<double> r_new_;
std::vector<double> J_;
std::vector<double> pvec_;
};
} // namespace aare::detail
src/MinuitSteps.hpp
+126
View File
@@ -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 <algorithm>
#include <array>
#include <cmath>
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<double, 1> x, NDView<double, 1> 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 <typename Model> struct Steps;
template <> struct Steps<model::Pol1> {
static void compute(const std::array<double, 2> &start,
[[maybe_unused]] double x_range, double y_range,
double slope_scale, std::array<double, 2> &steps) {
steps[0] = std::max(0.1 * std::abs(start[0]), 0.1 * y_range);
steps[1] = 0.1 * slope_scale;
}
};
template <> struct Steps<model::Pol2> {
static void compute(const std::array<double, 3> &start, double x_range,
double y_range, double slope_scale,
std::array<double, 3> &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<model::Gaussian> {
static void compute(const std::array<double, 3> &start, double x_range,
double y_range, double /*slope_scale*/,
std::array<double, 3> &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<model::GaussianErfcPlateau> {
static void compute(const std::array<double, 4> &start, double x_range,
double y_range, double /*slope_scale*/,
std::array<double, 4> &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<model::GaussianChargeSharing> {
static void compute(const std::array<double, 6> &start, double x_range,
double y_range, double slope_scale,
std::array<double, 6> &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<model::GaussianChargeSharingKb> {
static void compute(const std::array<double, 8> &start, double x_range,
double y_range, double slope_scale,
std::array<double, 8> &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<model::RisingScurve> {
static void compute(const std::array<double, 6> &start, double x_range,
double y_range, double slope_scale,
std::array<double, 6> &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<model::FallingScurve> {
static void compute(const std::array<double, 6> &start, double x_range,
double y_range, double slope_scale,
std::array<double, 6> &steps) {
Steps<model::RisingScurve>::compute(start, x_range, y_range,
slope_scale, steps);
}
};
} // namespace aare::minuit2