mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2026-09-04 21:50:54 +02:00
Removing deprecated use of lmfit (#344)
Removed the lmfit based fitting and the old API for calling functions closes #296
This commit is contained in:
Erik Fröjdh
+4
-60
@@ -62,7 +62,6 @@ 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_LMFIT "Use FetchContent to download lmfit" ON)
|
||||
option(AARE_FETCH_MINUIT2 "Use FetchContent to download Minuit2" ON)
|
||||
|
||||
option(AARE_WARNINGS_AS_ERRORS "Treat warnings as errors during compilation"
|
||||
@@ -87,8 +86,8 @@ if(AARE_SYSTEM_LIBRARIES)
|
||||
set(AARE_FETCH_ZMQ
|
||||
OFF
|
||||
CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
|
||||
# Still fetch lmfit and Minuit2 when setting AARE_SYSTEM_LIBRARIES since these
|
||||
# are not available on conda-forge
|
||||
# Still fetch Minuit2 when setting AARE_SYSTEM_LIBRARIES since it is not
|
||||
# available on conda-forge
|
||||
endif()
|
||||
|
||||
if(AARE_BENCHMARKS)
|
||||
@@ -97,61 +96,6 @@ endif()
|
||||
|
||||
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
|
||||
|
||||
if(AARE_FETCH_LMFIT)
|
||||
# TODO! Should we fetch lmfit from the web or inlcude a tar.gz in the repo?
|
||||
set(LMFIT_PATCH_COMMAND git apply
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch)
|
||||
|
||||
# For cmake < 3.28 we can't supply EXCLUDE_FROM_ALL to FetchContent_Declare so
|
||||
# we need this workaround
|
||||
if(${CMAKE_VERSION} VERSION_LESS "3.28")
|
||||
FetchContent_Declare(
|
||||
lmfit
|
||||
GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
|
||||
GIT_TAG main
|
||||
PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
|
||||
UPDATE_DISCONNECTED 1)
|
||||
else()
|
||||
FetchContent_Declare(
|
||||
lmfit
|
||||
GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
|
||||
GIT_TAG main
|
||||
PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
|
||||
UPDATE_DISCONNECTED 1
|
||||
EXCLUDE_FROM_ALL 1)
|
||||
endif()
|
||||
|
||||
# Disable what we don't need from lmfit
|
||||
set(BUILD_TESTING
|
||||
OFF
|
||||
CACHE BOOL "")
|
||||
set(LMFIT_CPPTEST
|
||||
OFF
|
||||
CACHE BOOL "")
|
||||
set(LIB_MAN
|
||||
OFF
|
||||
CACHE BOOL "")
|
||||
set(LMFIT_CPPTEST
|
||||
OFF
|
||||
CACHE BOOL "")
|
||||
set(BUILD_SHARED_LIBS
|
||||
OFF
|
||||
CACHE BOOL "")
|
||||
|
||||
if(${CMAKE_VERSION} VERSION_LESS "3.28")
|
||||
if(NOT lmfit_POPULATED)
|
||||
FetchContent_Populate(lmfit)
|
||||
add_subdirectory(${lmfit_SOURCE_DIR} ${lmfit_BINARY_DIR} EXCLUDE_FROM_ALL)
|
||||
endif()
|
||||
else()
|
||||
FetchContent_MakeAvailable(lmfit)
|
||||
endif()
|
||||
|
||||
set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
|
||||
else()
|
||||
find_package(lmfit REQUIRED)
|
||||
endif()
|
||||
|
||||
if(AARE_FETCH_MINUIT2)
|
||||
# We are building Minuit2 from sources.
|
||||
|
||||
@@ -457,8 +401,7 @@ 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:lmfit>
|
||||
$<BUILD_INTERFACE:aare::Minuit2>)
|
||||
PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:aare::Minuit2>)
|
||||
|
||||
target_include_directories(
|
||||
aare_core SYSTEM
|
||||
@@ -490,6 +433,7 @@ if(AARE_TESTS)
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolation.test.cpp
|
||||
|
||||
+12
-1
@@ -1,5 +1,17 @@
|
||||
# Release notes
|
||||
|
||||
## Next
|
||||
|
||||
### API Changes:
|
||||
|
||||
- Removed the lmfit dependency and the legacy ``fit_gaus``, ``fit_pol1``,
|
||||
``fit_scurve``, and ``fit_scurve2`` APIs. Use ``Gaussian``, ``Pol1``,
|
||||
``RisingScurve``, or ``FallingScurve`` and call ``model.fit(...)`` (or
|
||||
``fit(model, ...)``) instead.
|
||||
- Removed the legacy ``gaus``, ``pol1``, ``scurve``, and ``scurve2`` function
|
||||
evaluators. Model objects are callable and provide the replacement, for
|
||||
example ``Gaussian()(x, par)``.
|
||||
|
||||
## 2026.7.2
|
||||
|
||||
|
||||
@@ -151,4 +163,3 @@ dhanya.thattil@psi.ch
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -75,24 +75,7 @@ static void report_accuracy(benchmark::State &state, const TestCase &tc,
|
||||
// Benchmarks
|
||||
// ----------
|
||||
|
||||
// 1. lmcurve
|
||||
static void BM_FitGausLm(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();
|
||||
|
||||
aare::NDArray<double, 1> result;
|
||||
for (auto _ : state) {
|
||||
result = aare::fit_gaus(xv, yv);
|
||||
benchmark::DoNotOptimize(result.data());
|
||||
}
|
||||
|
||||
report_accuracy(state, tc, result);
|
||||
state.SetLabel(tc.name);
|
||||
}
|
||||
|
||||
// 2. Minuit2, analytic gradient (no Hesse)
|
||||
// Minuit2, analytic gradient (no Hesse)
|
||||
static void BM_FitGausMinuitGrad(benchmark::State &state) {
|
||||
const auto &tc = get_test_cases()[state.range(0)];
|
||||
auto data = generate_gaussian_data(tc);
|
||||
@@ -115,7 +98,7 @@ static void BM_FitGausMinuitGrad(benchmark::State &state) {
|
||||
state.SetLabel(tc.name);
|
||||
}
|
||||
|
||||
// 3. Minuit2, analytic gradient + Hesse
|
||||
// Minuit2, analytic gradient + Hesse
|
||||
static void BM_FitGausMinuitGradHesse(benchmark::State &state) {
|
||||
const auto &tc = get_test_cases()[state.range(0)];
|
||||
auto data = generate_gaussian_data(tc);
|
||||
@@ -145,8 +128,6 @@ static void BM_FitGausMinuitGradHesse(benchmark::State &state) {
|
||||
state.SetLabel(tc.name);
|
||||
}
|
||||
|
||||
BENCHMARK(BM_FitGausLm)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond);
|
||||
|
||||
BENCHMARK(BM_FitGausMinuitGrad)
|
||||
->DenseRange(0, 5)
|
||||
->Unit(benchmark::kMicrosecond);
|
||||
@@ -155,4 +136,4 @@ BENCHMARK(BM_FitGausMinuitGradHesse)
|
||||
->DenseRange(0, 5)
|
||||
->Unit(benchmark::kMicrosecond);
|
||||
|
||||
BENCHMARK_MAIN();
|
||||
BENCHMARK_MAIN();
|
||||
|
||||
@@ -14,7 +14,7 @@ Requirements
|
||||
To simplify deployment we build and statically link a few libraries.
|
||||
|
||||
- fmt
|
||||
- lmfit - https://jugit.fz-juelich.de/mlz/lmfit
|
||||
- Minuit2
|
||||
- nlohmann_json
|
||||
- pybind11
|
||||
- ZeroMQ
|
||||
@@ -23,4 +23,4 @@ To simplify deployment we build and statically link a few libraries.
|
||||
|
||||
- Sphinx
|
||||
- Breathe
|
||||
- Doxygen
|
||||
- Doxygen
|
||||
|
||||
+20
-9
@@ -1,19 +1,30 @@
|
||||
|
||||
Fit
|
||||
========
|
||||
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::
|
||||
|
||||
**Functions**
|
||||
model = Gaussian(compute_errors=True)
|
||||
result = model.fit(x, y, y_err)
|
||||
|
||||
.. autofunction:: gaus
|
||||
The model object is also callable, which evaluates it at the supplied points::
|
||||
|
||||
.. autofunction:: pol1
|
||||
fitted_y = model(x, result["par"])
|
||||
|
||||
The available models are ``Gaussian``, ``GaussianErfcPlateau``,
|
||||
``GaussianChargeSharing``, ``GaussianChargeSharingKb``, ``Pol1``, ``Pol2``,
|
||||
``RisingScurve``, and ``FallingScurve``. The module-level :func:`fit` function
|
||||
accepts the same model objects when a functional interface is preferred.
|
||||
|
||||
**Fitting**
|
||||
For three-dimensional data, pass an array with shape
|
||||
``(rows, columns, scan_points)`` and select the worker count with
|
||||
``n_threads``::
|
||||
|
||||
.. autofunction:: fit_gaus
|
||||
result = model.fit(x, image_data, image_errors, n_threads=8)
|
||||
|
||||
.. autofunction:: fit_pol1
|
||||
The result dictionary contains ``par`` and ``chi2``. It also contains
|
||||
``par_err`` when ``compute_errors`` is enabled.
|
||||
|
||||
.. autofunction:: fit
|
||||
|
||||
@@ -1,111 +1,11 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#pragma once
|
||||
|
||||
#include <cmath>
|
||||
#include <vector>
|
||||
|
||||
#include "aare/FitModel.hpp"
|
||||
#include "aare/NDArray.hpp"
|
||||
#include "aare/utils/par.hpp"
|
||||
#include "aare/utils/task.hpp"
|
||||
|
||||
namespace aare {
|
||||
|
||||
namespace func {
|
||||
double gaus(const double x, const double *par);
|
||||
NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par);
|
||||
|
||||
double pol1(const double x, const double *par);
|
||||
NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
|
||||
|
||||
double scurve(const double x, const double *par);
|
||||
NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par);
|
||||
|
||||
double scurve2(const double x, const double *par);
|
||||
NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
|
||||
|
||||
} // namespace func
|
||||
|
||||
static constexpr int DEFAULT_NUM_THREADS = 4;
|
||||
|
||||
/**
|
||||
* @brief Fit a 1D Gaussian to data.
|
||||
* @param data data to fit
|
||||
* @param x x values
|
||||
*/
|
||||
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
|
||||
|
||||
/**
|
||||
* @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
|
||||
* @param x x values
|
||||
* @param y y values, layout [row, col, values]
|
||||
* @param n_threads number of threads to use
|
||||
*/
|
||||
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads = DEFAULT_NUM_THREADS);
|
||||
|
||||
/**
|
||||
* @brief Fit a 1D Gaussian with error estimates
|
||||
* @param x x values
|
||||
* @param y y values, layout [row, col, values]
|
||||
* @param y_err error in y, layout [row, col, values]
|
||||
* @param par_out output parameters
|
||||
* @param par_err_out output error parameters
|
||||
*/
|
||||
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
|
||||
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
|
||||
double &chi2);
|
||||
|
||||
/**
|
||||
* @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
|
||||
* [row, col, values]
|
||||
* @param x x values
|
||||
* @param y y values, layout [row, col, values]
|
||||
* @param y_err error in y, layout [row, col, values]
|
||||
* @param par_out output parameters, layout [row, col, values]
|
||||
* @param par_err_out output parameter errors, layout [row, col, values]
|
||||
* @param n_threads number of threads to use
|
||||
*/
|
||||
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
|
||||
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
|
||||
NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
|
||||
|
||||
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y);
|
||||
|
||||
NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads = DEFAULT_NUM_THREADS);
|
||||
|
||||
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
|
||||
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
|
||||
double &chi2);
|
||||
|
||||
// TODO! not sure we need to offer the different version in C++
|
||||
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
|
||||
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
|
||||
NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
|
||||
|
||||
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y);
|
||||
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads);
|
||||
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
|
||||
NDView<double, 1> y_err, NDView<double, 1> par_out,
|
||||
NDView<double, 1> par_err_out, double &chi2);
|
||||
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
|
||||
NDView<double, 3> y_err, NDView<double, 3> par_out,
|
||||
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
|
||||
int n_threads);
|
||||
|
||||
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y);
|
||||
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads);
|
||||
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
|
||||
NDView<double, 1> y_err, NDView<double, 1> par_out,
|
||||
NDView<double, 1> par_err_out, double &chi2);
|
||||
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
|
||||
NDView<double, 3> y_err, NDView<double, 3> par_out,
|
||||
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
|
||||
int n_threads);
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Minuit2-based pixel fitting.
|
||||
// Template bodies and explicit instantiations live in src/Fit.cpp.
|
||||
|
||||
@@ -1,13 +0,0 @@
|
||||
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
|
||||
index 4efb7ed..6533660 100644
|
||||
--- a/lib/CMakeLists.txt
|
||||
+++ b/lib/CMakeLists.txt
|
||||
@@ -11,7 +11,7 @@ target_compile_definitions(${lib} PRIVATE "LMFIT_EXPORT") # for Windows DLL expo
|
||||
|
||||
target_include_directories(${lib}
|
||||
PUBLIC
|
||||
- $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/>
|
||||
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/>
|
||||
$<INSTALL_INTERFACE:include/>
|
||||
)
|
||||
|
||||
@@ -41,7 +41,6 @@ set(PYTHON_FILES
|
||||
aare/Cluster.py
|
||||
aare/calibration.py
|
||||
aare/experimental.py
|
||||
aare/func.py
|
||||
aare/RawFile.py
|
||||
aare/transform.py
|
||||
aare/ScanParameters.py
|
||||
|
||||
@@ -26,7 +26,6 @@ from .Cluster import Cluster
|
||||
|
||||
from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb
|
||||
from ._aare import fit
|
||||
from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
|
||||
from ._aare import Interpolator
|
||||
from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2
|
||||
from ._aare import reduce_to_2x2, reduce_to_3x3
|
||||
@@ -42,9 +41,6 @@ from .ScanParameters import ScanParameters
|
||||
from .utils import random_pixels, random_pixel, flat_list, add_colorbar, Timer
|
||||
|
||||
|
||||
#make functions available in the top level API
|
||||
from .func import *
|
||||
|
||||
from .calibration import *
|
||||
from ._aare import apply_calibration, count_switching_pixels
|
||||
from ._aare import calculate_pedestal, calculate_pedestal_float, calculate_pedestal_g0, calculate_pedestal_g0_float
|
||||
|
||||
@@ -1,2 +0,0 @@
|
||||
# SPDX-License-Identifier: MPL-2.0
|
||||
from ._aare import gaus, pol1, scurve, scurve2
|
||||
+42
-70
@@ -1,106 +1,78 @@
|
||||
# SPDX-License-Identifier: MPL-2.0
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import sys
|
||||
sys.path.insert(0, '/home/kferjaoui/sw/aare/build')
|
||||
from aare import fit_gaus, fit_pol1
|
||||
from aare import Gaussian, fit
|
||||
from aare import pol1
|
||||
|
||||
textpm = f"±" #
|
||||
textmu = f"μ" #
|
||||
textsigma = f"σ" #
|
||||
from aare import Gaussian, Pol1
|
||||
|
||||
textpm = "±"
|
||||
textmu = "μ"
|
||||
textsigma = "σ"
|
||||
|
||||
|
||||
# ================================= Gauss fit =================================
|
||||
# Parameters
|
||||
mu = np.random.uniform(1, 100) # Mean of Gaussian
|
||||
sigma = np.random.uniform(4, 20) # Standard deviation
|
||||
num_points = 10000 # Number of points for smooth distribution
|
||||
# noise_sigma = 10
|
||||
|
||||
# Generate Gaussian distribution
|
||||
data = np.random.normal(mu, sigma, num_points)
|
||||
mu = np.random.uniform(1, 100)
|
||||
sigma = np.random.uniform(4, 20)
|
||||
data = np.random.normal(mu, sigma, 10000)
|
||||
counts, edges = np.histogram(data, bins=100)
|
||||
|
||||
x = 0.5 * (edges[:-1] + edges[1:]) # proper bin centers
|
||||
x = 0.5 * (edges[:-1] + edges[1:])
|
||||
y = counts.astype(np.float64)
|
||||
|
||||
# Poisson noise
|
||||
yerr = np.sqrt(np.maximum(y, 1))
|
||||
# yerr = np.abs(np.random.normal(0, noise_sigma, len(x)))
|
||||
|
||||
# Create subplot
|
||||
fig0, ax0 = plt.subplots(1, 1, num=0, figsize=(12, 8))
|
||||
|
||||
# Add the errors as error bars in the step plot
|
||||
ax0.errorbar(x, y, yerr=yerr, fmt=". ", capsize=5)
|
||||
ax0.grid()
|
||||
|
||||
# Fit with lmfit
|
||||
result_lm = fit_gaus(x, y, yerr)
|
||||
par_lm = result_lm["par"]
|
||||
err_lm = result_lm["par_err"]
|
||||
chi2_lm = result_lm["chi2"]
|
||||
print("[lmfit] fit_gaus: ", par_lm, err_lm, chi2_lm)
|
||||
gaussian = Gaussian(compute_errors=True)
|
||||
result = gaussian.fit(x, y, yerr)
|
||||
par = result["par"]
|
||||
err = result["par_err"]
|
||||
chi2 = result["chi2"]
|
||||
print(f"Gaussian.fit: par={par}, err={err}, chi2={chi2}")
|
||||
|
||||
# Fit with Minuit2 + analytic gradient + Hesse errors
|
||||
gaussian = Gaussian()
|
||||
gaussian.compute_errors = True
|
||||
result_m2 = gaussian.fit(x, y, yerr)
|
||||
par_m2 = result_m2['par']
|
||||
err_m2 = result_m2['par_err']
|
||||
chi2_m2 = result_m2['chi2']
|
||||
print(f"[minuit2] gaussian.fit: par={par_m2}, err={err_m2}, chi2={chi2_m2}")
|
||||
|
||||
x = np.linspace(x[0], x[-1], 1000)
|
||||
ax0.plot(x, gaussian(x, par_lm), marker="", label="fit_gaus")
|
||||
ax0.plot(x, gaussian(x, par_m2), marker="", linestyle=":", label="fit_gaus_minuit_grad")
|
||||
x_plot = np.linspace(x[0], x[-1], 1000)
|
||||
ax0.plot(x_plot, gaussian(x_plot, par), marker="", label="Gaussian.fit")
|
||||
ax0.legend()
|
||||
ax0.set(xlabel="x", ylabel="Counts",
|
||||
ax0.set(
|
||||
xlabel="x",
|
||||
ylabel="Counts",
|
||||
title=(
|
||||
f"fit_gaus: A={par_lm[0]:0.2f}{textpm}{err_lm[0]:0.2f} "
|
||||
f"{textmu}={par_lm[1]:0.2f}{textpm}{err_lm[1]:0.2f} "
|
||||
f"{textsigma}={par_lm[2]:0.2f}{textpm}{err_lm[2]:0.2f}\n"
|
||||
f"minuit_grad: A={par_m2[0]:0.2f}{textpm}{err_m2[0]:0.2f} "
|
||||
f"{textmu}={par_m2[1]:0.2f}{textpm}{err_m2[1]:0.2f} "
|
||||
f"{textsigma}={par_m2[2]:0.2f}{textpm}{err_m2[2]:0.2f}\n"
|
||||
f"A={par[0]:0.2f}{textpm}{err[0]:0.2f} "
|
||||
f"{textmu}={par[1]:0.2f}{textpm}{err[1]:0.2f} "
|
||||
f"{textsigma}={par[2]:0.2f}{textpm}{err[2]:0.2f}\n"
|
||||
f"(truth: {textmu}={mu:0.2f}, {textsigma}={sigma:0.2f})"
|
||||
),
|
||||
)
|
||||
fig0.tight_layout()
|
||||
|
||||
|
||||
|
||||
# ================================= pol1 fit =================================
|
||||
# Parameters
|
||||
# ================================= Pol1 fit =================================
|
||||
n_points = 40
|
||||
|
||||
# Generate random slope and intercept (origin)
|
||||
slope = np.random.uniform(-10, 10) # Random slope between 0.5 and 2.0
|
||||
intercept = np.random.uniform(-10, 10) # Random intercept between -10 and 10
|
||||
|
||||
# Generate random x values
|
||||
slope = np.random.uniform(-10, 10)
|
||||
intercept = np.random.uniform(-10, 10)
|
||||
x_values = np.random.uniform(-10, 10, n_points)
|
||||
|
||||
# Calculate y values based on the linear function y = mx + b + error
|
||||
errors = np.abs(np.random.normal(0, np.random.uniform(1, 5), n_points))
|
||||
var_points = np.random.normal(0, np.random.uniform(0.1, 2), n_points)
|
||||
y_values = slope * x_values + intercept + var_points
|
||||
y_values = slope * x_values + intercept + np.random.normal(0, 1, n_points)
|
||||
|
||||
fig1, ax1 = plt.subplots(1, 1, num=1, figsize=(12, 8))
|
||||
ax1.errorbar(x_values, y_values, yerr=errors, fmt=". ", capsize=5)
|
||||
result_pol = fit_pol1(x_values, y_values, errors)
|
||||
par = result_pol["par"]
|
||||
err = result_pol["par_err"]
|
||||
|
||||
x = np.linspace(np.min(x_values), np.max(x_values), 1000)
|
||||
ax1.plot(x, pol1(x, par), marker="")
|
||||
ax1.set(xlabel="x", ylabel="y", title=f"a = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
|
||||
f"b = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
|
||||
f"(init: {slope:0.2f}, {intercept:0.2f})")
|
||||
pol1 = Pol1(compute_errors=True)
|
||||
result = pol1.fit(x_values, y_values, errors)
|
||||
par = result["par"]
|
||||
err = result["par_err"]
|
||||
|
||||
x_plot = np.linspace(np.min(x_values), np.max(x_values), 1000)
|
||||
ax1.plot(x_plot, pol1(x_plot, par), marker="")
|
||||
ax1.set(
|
||||
xlabel="x",
|
||||
ylabel="y",
|
||||
title=(
|
||||
f"intercept = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
|
||||
f"slope = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
|
||||
f"(truth: {intercept:0.2f}, {slope:0.2f})"
|
||||
),
|
||||
)
|
||||
fig1.tight_layout()
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
+1
-447
@@ -227,452 +227,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
|
||||
}
|
||||
|
||||
void define_fit_bindings(py::module &m) {
|
||||
|
||||
// TODO! Evaluate without converting to double
|
||||
m.def(
|
||||
"gaus",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
|
||||
auto x_view = make_view_1d(x);
|
||||
auto par_view = make_view_1d(par);
|
||||
auto y = new NDArray<double, 1>{aare::func::gaus(x_view, par_view)};
|
||||
return return_image_data(y);
|
||||
},
|
||||
R"(
|
||||
Evaluate a 1D Gaussian function for all points in x using parameters par.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The points at which to evaluate the Gaussian function.
|
||||
par : array_like
|
||||
The parameters of the Gaussian function. The first element is the amplitude, the second element is the mean, and the third element is the standard deviation.
|
||||
)",
|
||||
py::arg("x"), py::arg("par"));
|
||||
|
||||
m.def(
|
||||
"pol1",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
|
||||
auto x_view = make_view_1d(x);
|
||||
auto par_view = make_view_1d(par);
|
||||
auto y = new NDArray<double, 1>{aare::func::pol1(x_view, par_view)};
|
||||
return return_image_data(y);
|
||||
},
|
||||
R"(
|
||||
Evaluate a 1D polynomial function for all points in x using parameters par. (p0+p1*x)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The points at which to evaluate the polynomial function.
|
||||
par : array_like
|
||||
The parameters of the polynomial function. The first element is the intercept, and the second element is the slope.
|
||||
)",
|
||||
py::arg("x"), py::arg("par"));
|
||||
|
||||
m.def(
|
||||
"scurve",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
|
||||
auto x_view = make_view_1d(x);
|
||||
auto par_view = make_view_1d(par);
|
||||
auto y =
|
||||
new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
|
||||
return return_image_data(y);
|
||||
},
|
||||
R"(
|
||||
Evaluate a 1D scurve function for all points in x using parameters par.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The points at which to evaluate the scurve function.
|
||||
par : array_like
|
||||
The parameters of the scurve function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.
|
||||
)",
|
||||
py::arg("x"), py::arg("par"));
|
||||
|
||||
m.def(
|
||||
"scurve2",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
|
||||
auto x_view = make_view_1d(x);
|
||||
auto par_view = make_view_1d(par);
|
||||
auto y =
|
||||
new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
|
||||
return return_image_data(y);
|
||||
},
|
||||
R"(
|
||||
Evaluate a 1D scurve2 function for all points in x using parameters par.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The points at which to evaluate the scurve function.
|
||||
par : array_like
|
||||
The parameters of the scurve2 function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.
|
||||
)",
|
||||
py::arg("x"), py::arg("par"));
|
||||
|
||||
m.def(
|
||||
"fit_gaus",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
auto par = new NDArray<double, 3>{};
|
||||
auto y_view = make_view_3d(y);
|
||||
auto x_view = make_view_1d(x);
|
||||
*par = aare::fit_gaus(x_view, y_view, n_threads);
|
||||
return return_image_data(par);
|
||||
} else if (y.ndim() == 1) {
|
||||
auto par = new NDArray<double, 1>{};
|
||||
auto y_view = make_view_1d(y);
|
||||
auto x_view = make_view_1d(x);
|
||||
*par = aare::fit_gaus(x_view, y_view);
|
||||
return return_image_data(par);
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
R"(
|
||||
Fit a 1D Gaussian to data.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The x values.
|
||||
y : array_like
|
||||
The y values.
|
||||
n_threads : int, optional
|
||||
The number of threads to use. Default is 4.
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
|
||||
|
||||
m.def(
|
||||
"fit_gaus",
|
||||
[](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::array_t<double, py::array::c_style | py::array::forcecast> y_err,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
// Allocate memory for the output
|
||||
// Need to have pointers to allow python to manage
|
||||
// the memory
|
||||
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
|
||||
auto par_err =
|
||||
new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
|
||||
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
|
||||
|
||||
// Make views of the numpy arrays
|
||||
auto y_view = make_view_3d(y);
|
||||
auto y_view_err = make_view_3d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2->view(), n_threads);
|
||||
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = return_image_data(chi2),
|
||||
"Ndf"_a = y.shape(2) - 3);
|
||||
} else if (y.ndim() == 1) {
|
||||
// Allocate memory for the output
|
||||
// Need to have pointers to allow python to manage
|
||||
// the memory
|
||||
auto par = new NDArray<double, 1>({3});
|
||||
auto par_err = new NDArray<double, 1>({3});
|
||||
|
||||
// Decode the numpy arrays
|
||||
auto y_view = make_view_1d(y);
|
||||
auto y_view_err = make_view_1d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
double chi2 = 0;
|
||||
aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2);
|
||||
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = chi2, "Ndf"_a = y.size() - 3);
|
||||
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
R"(
|
||||
Fit a 1D Gaussian to data with error estimates.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The x values.
|
||||
y : array_like
|
||||
The y values.
|
||||
y_err : array_like
|
||||
The error in the y values.
|
||||
n_threads : int, optional
|
||||
The number of threads to use. Default is 4.
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
|
||||
|
||||
m.def(
|
||||
"fit_pol1",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
auto par = new NDArray<double, 3>{};
|
||||
|
||||
auto x_view = make_view_1d(x);
|
||||
auto y_view = make_view_3d(y);
|
||||
*par = aare::fit_pol1(x_view, y_view, n_threads);
|
||||
return return_image_data(par);
|
||||
} else if (y.ndim() == 1) {
|
||||
auto par = new NDArray<double, 1>{};
|
||||
auto x_view = make_view_1d(x);
|
||||
auto y_view = make_view_1d(y);
|
||||
*par = aare::fit_pol1(x_view, y_view);
|
||||
return return_image_data(par);
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
|
||||
|
||||
m.def(
|
||||
"fit_pol1",
|
||||
[](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::array_t<double, py::array::c_style | py::array::forcecast> y_err,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
|
||||
|
||||
auto par_err =
|
||||
new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
|
||||
|
||||
auto y_view = make_view_3d(y);
|
||||
auto y_view_err = make_view_3d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
|
||||
|
||||
aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2->view(), n_threads);
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = return_image_data(chi2),
|
||||
"Ndf"_a = y.shape(2) - 2);
|
||||
|
||||
} else if (y.ndim() == 1) {
|
||||
auto par = new NDArray<double, 1>({2});
|
||||
auto par_err = new NDArray<double, 1>({2});
|
||||
|
||||
auto y_view = make_view_1d(y);
|
||||
auto y_view_err = make_view_1d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
double chi2 = 0;
|
||||
|
||||
aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2);
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
|
||||
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
R"(
|
||||
Fit a 1D polynomial to data with error estimates.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The x values.
|
||||
y : array_like
|
||||
The y values.
|
||||
y_err : array_like
|
||||
The error in the y values.
|
||||
n_threads : int, optional
|
||||
The number of threads to use. Default is 4.
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
|
||||
|
||||
//=========
|
||||
m.def(
|
||||
"fit_scurve",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
auto par = new NDArray<double, 3>{};
|
||||
|
||||
auto x_view = make_view_1d(x);
|
||||
auto y_view = make_view_3d(y);
|
||||
*par = aare::fit_scurve(x_view, y_view, n_threads);
|
||||
return return_image_data(par);
|
||||
} else if (y.ndim() == 1) {
|
||||
auto par = new NDArray<double, 1>{};
|
||||
auto x_view = make_view_1d(x);
|
||||
auto y_view = make_view_1d(y);
|
||||
*par = aare::fit_scurve(x_view, y_view);
|
||||
return return_image_data(par);
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
|
||||
|
||||
m.def(
|
||||
"fit_scurve",
|
||||
[](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::array_t<double, py::array::c_style | py::array::forcecast> y_err,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
|
||||
|
||||
auto par_err =
|
||||
new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
|
||||
|
||||
auto y_view = make_view_3d(y);
|
||||
auto y_view_err = make_view_3d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
|
||||
|
||||
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2->view(), n_threads);
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = return_image_data(chi2),
|
||||
"Ndf"_a = y.shape(2) - 2);
|
||||
|
||||
} else if (y.ndim() == 1) {
|
||||
auto par = new NDArray<double, 1>({2});
|
||||
auto par_err = new NDArray<double, 1>({2});
|
||||
|
||||
auto y_view = make_view_1d(y);
|
||||
auto y_view_err = make_view_1d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
double chi2 = 0;
|
||||
|
||||
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2);
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
|
||||
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
R"(
|
||||
Fit a 1D polynomial to data with error estimates.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The x values.
|
||||
y : array_like
|
||||
The y values.
|
||||
y_err : array_like
|
||||
The error in the y values.
|
||||
n_threads : int, optional
|
||||
The number of threads to use. Default is 4.
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
|
||||
|
||||
m.def(
|
||||
"fit_scurve2",
|
||||
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
auto par = new NDArray<double, 3>{};
|
||||
|
||||
auto x_view = make_view_1d(x);
|
||||
auto y_view = make_view_3d(y);
|
||||
*par = aare::fit_scurve2(x_view, y_view, n_threads);
|
||||
return return_image_data(par);
|
||||
} else if (y.ndim() == 1) {
|
||||
auto par = new NDArray<double, 1>{};
|
||||
auto x_view = make_view_1d(x);
|
||||
auto y_view = make_view_1d(y);
|
||||
*par = aare::fit_scurve2(x_view, y_view);
|
||||
return return_image_data(par);
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
|
||||
|
||||
m.def(
|
||||
"fit_scurve2",
|
||||
[](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::array_t<double, py::array::c_style | py::array::forcecast> y_err,
|
||||
int n_threads) {
|
||||
if (y.ndim() == 3) {
|
||||
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
|
||||
|
||||
auto par_err =
|
||||
new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
|
||||
|
||||
auto y_view = make_view_3d(y);
|
||||
auto y_view_err = make_view_3d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
|
||||
|
||||
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2->view(), n_threads);
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = return_image_data(chi2),
|
||||
"Ndf"_a = y.shape(2) - 2);
|
||||
|
||||
} else if (y.ndim() == 1) {
|
||||
auto par = new NDArray<double, 1>({6});
|
||||
auto par_err = new NDArray<double, 1>({6});
|
||||
|
||||
auto y_view = make_view_1d(y);
|
||||
auto y_view_err = make_view_1d(y_err);
|
||||
auto x_view = make_view_1d(x);
|
||||
|
||||
double chi2 = 0;
|
||||
|
||||
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
|
||||
par_err->view(), chi2);
|
||||
return py::dict("par"_a = return_image_data(par),
|
||||
"par_err"_a = return_image_data(par_err),
|
||||
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
|
||||
|
||||
} else {
|
||||
throw std::runtime_error("Data must be 1D or 3D");
|
||||
}
|
||||
},
|
||||
R"(
|
||||
Fit a 1D polynomial to data with error estimates.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
The x values.
|
||||
y : array_like
|
||||
The y values.
|
||||
y_err : array_like
|
||||
The error in the y values.
|
||||
n_threads : int, optional
|
||||
The number of threads to use. Default is 4.
|
||||
)",
|
||||
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
|
||||
|
||||
// ── Bind model classes ──────────────────────────────────────────
|
||||
bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
|
||||
bind_fit_model<aare::model::GaussianErfcPlateau>(m, "GaussianErfcPlateau");
|
||||
@@ -794,4 +348,4 @@ void define_fit_bindings(py::module &m) {
|
||||
)",
|
||||
py::arg("model"), py::arg("x"), py::arg("y"),
|
||||
py::arg("y_err") = py::none(), py::arg("n_threads") = 4);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -0,0 +1,15 @@
|
||||
# SPDX-License-Identifier: MPL-2.0
|
||||
import numpy as np
|
||||
import aare
|
||||
|
||||
|
||||
|
||||
def test_gaussian_model_evaluates_and_fits_data():
|
||||
x = np.linspace(-5.0, 5.0, 51)
|
||||
expected = np.array([20.0, 0.5, 1.2])
|
||||
model = aare.Gaussian()
|
||||
|
||||
y = model(x, expected)
|
||||
result = model.fit(x, y)
|
||||
|
||||
np.testing.assert_allclose(result["par"], expected, atol=2e-3)
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+1
-452
@@ -4,7 +4,6 @@
|
||||
#include "Minuit2/FunctionMinimum.h"
|
||||
#include "Minuit2/MnHesse.h"
|
||||
#include "Minuit2/MnMigrad.h"
|
||||
#include "Minuit2/MnPrint.h"
|
||||
#include "Minuit2/MnStrategy.h"
|
||||
#include "Minuit2/MnUserParameters.h"
|
||||
#include "aare/Models.hpp"
|
||||
@@ -12,461 +11,11 @@
|
||||
#include "aare/utils/task.hpp"
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <lmcurve2.h>
|
||||
#include <lmfit.hpp>
|
||||
#include <memory>
|
||||
#include <stdexcept>
|
||||
#include <thread>
|
||||
#include <type_traits>
|
||||
|
||||
namespace aare {
|
||||
|
||||
namespace func {
|
||||
|
||||
double gaus(const double x, const double *par) {
|
||||
return par[0] * exp(-pow(x - par[1], 2) / (2 * pow(par[2], 2)));
|
||||
}
|
||||
|
||||
NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par) {
|
||||
NDArray<double, 1> y({x.shape(0)}, 0);
|
||||
for (ssize_t i = 0; i < x.size(); i++) {
|
||||
y(i) = gaus(x(i), par.data());
|
||||
}
|
||||
return y;
|
||||
}
|
||||
|
||||
double pol1(const double x, const double *par) { return par[0] * x + par[1]; }
|
||||
|
||||
NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
|
||||
NDArray<double, 1> y({x.shape()}, 0);
|
||||
for (ssize_t i = 0; i < x.size(); i++) {
|
||||
y(i) = pol1(x(i), par.data());
|
||||
}
|
||||
return y;
|
||||
}
|
||||
|
||||
double scurve(const double x, const double *par) {
|
||||
return (par[0] + par[1] * x) +
|
||||
0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) *
|
||||
(par[4] + par[5] * (x - par[2]));
|
||||
}
|
||||
|
||||
NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
|
||||
NDArray<double, 1> y({x.shape()}, 0);
|
||||
for (ssize_t i = 0; i < x.size(); i++) {
|
||||
y(i) = scurve(x(i), par.data());
|
||||
}
|
||||
return y;
|
||||
}
|
||||
|
||||
double scurve2(const double x, const double *par) {
|
||||
return (par[0] + par[1] * x) +
|
||||
0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) *
|
||||
(par[4] + par[5] * (x - par[2]));
|
||||
}
|
||||
|
||||
NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
|
||||
NDArray<double, 1> y({x.shape()}, 0);
|
||||
for (ssize_t i = 0; i < x.size(); i++) {
|
||||
y(i) = scurve2(x(i), par.data());
|
||||
}
|
||||
return y;
|
||||
}
|
||||
|
||||
} // namespace func
|
||||
|
||||
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
|
||||
NDArray<double, 1> result = model::Gaussian::estimate_par(x, y);
|
||||
lm_status_struct status;
|
||||
|
||||
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
|
||||
aare::func::gaus, &lm_control_double, &status);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads) {
|
||||
NDArray<double, 3> result({y.shape(0), y.shape(1), 3}, 0);
|
||||
|
||||
auto process = [&x, &y, &result](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)});
|
||||
auto res = fit_gaus(x, values);
|
||||
result(row, col, 0) = res(0);
|
||||
result(row, col, 1) = res(1);
|
||||
result(row, col, 2) = res(2);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
return result;
|
||||
}
|
||||
|
||||
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
|
||||
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
|
||||
double &chi2) {
|
||||
|
||||
// Check that we have the correct sizes
|
||||
if (y.size() != x.size() || y.size() != y_err.size() ||
|
||||
par_out.size() != 3 || par_err_out.size() != 3) {
|
||||
throw std::runtime_error("Data, x, data_err must have the same size "
|
||||
"and par_out, par_err_out must have size 3");
|
||||
}
|
||||
|
||||
// /* Collection of output parameters for status info. */
|
||||
// typedef struct {
|
||||
// double fnorm; /* norm of the residue vector fvec. */
|
||||
// int nfev; /* actual number of iterations. */
|
||||
// int outcome; /* Status indicator. Nonnegative values are used as
|
||||
// index
|
||||
// for the message text lm_infmsg, set in lmmin.c. */
|
||||
// int userbreak; /* Set when function evaluation requests termination.
|
||||
// */
|
||||
// } lm_status_struct;
|
||||
|
||||
lm_status_struct status;
|
||||
par_out = model::Gaussian::estimate_par(x, y);
|
||||
std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0, 0, 0};
|
||||
|
||||
// void lmcurve2( const int n_par, double *par, double *parerr, double
|
||||
// *covar, const int m_dat, const double *t, const double *y, const double
|
||||
// *dy, double (*f)( const double ti, const double *par ), const
|
||||
// lm_control_struct *control, lm_status_struct *status); n_par - Number of
|
||||
// free variables. Length of parameter vector par. par - Parameter vector.
|
||||
// On input, it must contain a reasonable guess. On output, it contains the
|
||||
// solution found to minimize ||r||. parerr - Parameter uncertainties
|
||||
// vector. Array of length n_par or NULL. On output, unless it or covar is
|
||||
// NULL, it contains the weighted parameter uncertainties for the found
|
||||
// parameters. covar - Covariance matrix. Array of length n_par * n_par or
|
||||
// NULL. On output, unless it is NULL, it contains the covariance matrix.
|
||||
// m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy
|
||||
// n_par <= m_dat. t - Array of length m_dat. Contains the abcissae (time,
|
||||
// or "x") for which function f will be evaluated. y - Array of length
|
||||
// m_dat. Contains the ordinate values that shall be fitted. dy - Array of
|
||||
// length m_dat. Contains the standard deviations of the values y. f - A
|
||||
// user-supplied parametric function f(ti;par). control - Parameter
|
||||
// collection for tuning the fit procedure. In most cases, the default
|
||||
// &lm_control_double is adequate. If f is only computed with
|
||||
// single-precision accuracy, &lm_control_float should be used. Parameters
|
||||
// are explained in lmmin2(3). status - A record used to return information
|
||||
// about the minimization process: For details, see lmmin2(3).
|
||||
|
||||
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
|
||||
x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
|
||||
&lm_control_double, &status);
|
||||
|
||||
// Calculate chi2
|
||||
chi2 = 0;
|
||||
for (ssize_t i = 0; i < y.size(); i++) {
|
||||
chi2 +=
|
||||
std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
|
||||
}
|
||||
}
|
||||
|
||||
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
|
||||
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
|
||||
NDView<double, 2> chi2_out,
|
||||
|
||||
int n_threads) {
|
||||
|
||||
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> y_view(&y(row, col, 0), {y.shape(2)});
|
||||
NDView<double, 1> y_err_view(&y_err(row, col, 0),
|
||||
{y_err.shape(2)});
|
||||
NDView<double, 1> par_out_view(&par_out(row, col, 0),
|
||||
{par_out.shape(2)});
|
||||
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
|
||||
{par_err_out.shape(2)});
|
||||
|
||||
fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view,
|
||||
chi2_out(row, col));
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
|
||||
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
|
||||
double &chi2) {
|
||||
|
||||
// Check that we have the correct sizes
|
||||
if (y.size() != x.size() || y.size() != y_err.size() ||
|
||||
par_out.size() != 2 || par_err_out.size() != 2) {
|
||||
throw std::runtime_error("Data, x, data_err must have the same size "
|
||||
"and par_out, par_err_out must have size 2");
|
||||
}
|
||||
|
||||
lm_status_struct status;
|
||||
par_out = model::Pol1::estimate_par(x, y);
|
||||
std::array<double, 4> cov{0, 0, 0, 0};
|
||||
|
||||
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
|
||||
x.size(), x.data(), y.data(), y_err.data(), aare::func::pol1,
|
||||
&lm_control_double, &status);
|
||||
|
||||
// Calculate chi2
|
||||
chi2 = 0;
|
||||
for (ssize_t i = 0; i < y.size(); i++) {
|
||||
chi2 +=
|
||||
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
|
||||
}
|
||||
}
|
||||
|
||||
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
|
||||
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
|
||||
NDView<double, 2> chi2_out, int n_threads) {
|
||||
|
||||
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> y_view(&y(row, col, 0), {y.shape(2)});
|
||||
NDView<double, 1> y_err_view(&y_err(row, col, 0),
|
||||
{y_err.shape(2)});
|
||||
NDView<double, 1> par_out_view(&par_out(row, col, 0),
|
||||
{par_out.shape(2)});
|
||||
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
|
||||
{par_err_out.shape(2)});
|
||||
|
||||
fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view,
|
||||
chi2_out(row, col));
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
|
||||
// // Check that we have the correct sizes
|
||||
// if (y.size() != x.size() || y.size() != y_err.size() ||
|
||||
// par_out.size() != 2 || par_err_out.size() != 2) {
|
||||
// throw std::runtime_error("Data, x, data_err must have the same size "
|
||||
// "and par_out, par_err_out must have size 2");
|
||||
// }
|
||||
NDArray<double, 1> par = model::Pol1::estimate_par(x, y);
|
||||
|
||||
lm_status_struct status;
|
||||
lmcurve(par.size(), par.data(), x.size(), x.data(), y.data(),
|
||||
aare::func::pol1, &lm_control_double, &status);
|
||||
|
||||
return par;
|
||||
}
|
||||
|
||||
NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads) {
|
||||
NDArray<double, 3> result({y.shape(0), y.shape(1), 2}, 0);
|
||||
|
||||
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)});
|
||||
auto res = fit_pol1(x, values);
|
||||
result(row, col, 0) = res(0);
|
||||
result(row, col, 1) = res(1);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
|
||||
RunInParallel(process, tasks);
|
||||
return result;
|
||||
}
|
||||
|
||||
// ~~ S-CURVES ~~
|
||||
|
||||
// - No error
|
||||
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
|
||||
NDArray<double, 1> result = model::RisingScurve::estimate_par(x, y);
|
||||
lm_status_struct status;
|
||||
|
||||
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
|
||||
aare::func::scurve, &lm_control_double, &status);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads) {
|
||||
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
|
||||
|
||||
auto process = [&x, &y, &result](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)});
|
||||
auto res = fit_scurve(x, values);
|
||||
result(row, col, 0) = res(0);
|
||||
result(row, col, 1) = res(1);
|
||||
result(row, col, 2) = res(2);
|
||||
result(row, col, 3) = res(3);
|
||||
result(row, col, 4) = res(4);
|
||||
result(row, col, 5) = res(5);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
return result;
|
||||
}
|
||||
|
||||
// - Error
|
||||
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
|
||||
NDView<double, 1> y_err, NDView<double, 1> par_out,
|
||||
NDView<double, 1> par_err_out, double &chi2) {
|
||||
|
||||
// Check that we have the correct sizes
|
||||
if (y.size() != x.size() || y.size() != y_err.size() ||
|
||||
par_out.size() != 6 || par_err_out.size() != 6) {
|
||||
throw std::runtime_error("Data, x, data_err must have the same size "
|
||||
"and par_out, par_err_out must have size 6");
|
||||
}
|
||||
|
||||
lm_status_struct status;
|
||||
par_out = model::RisingScurve::estimate_par(x, y);
|
||||
std::array<double, 36> cov = {0}; // size 6x6
|
||||
// std::array<double, 4> cov{0, 0, 0, 0};
|
||||
|
||||
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
|
||||
x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve,
|
||||
&lm_control_double, &status);
|
||||
|
||||
// Calculate chi2
|
||||
chi2 = 0;
|
||||
for (ssize_t i = 0; i < y.size(); i++) {
|
||||
chi2 +=
|
||||
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
|
||||
}
|
||||
}
|
||||
|
||||
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
|
||||
NDView<double, 3> y_err, NDView<double, 3> par_out,
|
||||
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
|
||||
int n_threads) {
|
||||
|
||||
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> y_view(&y(row, col, 0), {y.shape(2)});
|
||||
NDView<double, 1> y_err_view(&y_err(row, col, 0),
|
||||
{y_err.shape(2)});
|
||||
NDView<double, 1> par_out_view(&par_out(row, col, 0),
|
||||
{par_out.shape(2)});
|
||||
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
|
||||
{par_err_out.shape(2)});
|
||||
|
||||
fit_scurve(x, y_view, y_err_view, par_out_view,
|
||||
par_err_out_view, chi2_out(row, col));
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
// SCURVE2 ---
|
||||
|
||||
// - No error
|
||||
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
|
||||
NDArray<double, 1> result = model::FallingScurve::estimate_par(x, y);
|
||||
lm_status_struct status;
|
||||
|
||||
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
|
||||
aare::func::scurve2, &lm_control_double, &status);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
|
||||
int n_threads) {
|
||||
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
|
||||
|
||||
auto process = [&x, &y, &result](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)});
|
||||
auto res = fit_scurve2(x, values);
|
||||
result(row, col, 0) = res(0);
|
||||
result(row, col, 1) = res(1);
|
||||
result(row, col, 2) = res(2);
|
||||
result(row, col, 3) = res(3);
|
||||
result(row, col, 4) = res(4);
|
||||
result(row, col, 5) = res(5);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
return result;
|
||||
}
|
||||
|
||||
// - Error
|
||||
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
|
||||
NDView<double, 1> y_err, NDView<double, 1> par_out,
|
||||
NDView<double, 1> par_err_out, double &chi2) {
|
||||
|
||||
// Check that we have the correct sizes
|
||||
if (y.size() != x.size() || y.size() != y_err.size() ||
|
||||
par_out.size() != 6 || par_err_out.size() != 6) {
|
||||
throw std::runtime_error("Data, x, data_err must have the same size "
|
||||
"and par_out, par_err_out must have size 6");
|
||||
}
|
||||
|
||||
lm_status_struct status;
|
||||
par_out = model::FallingScurve::estimate_par(x, y);
|
||||
std::array<double, 36> cov = {0}; // size 6x6
|
||||
// std::array<double, 4> cov{0, 0, 0, 0};
|
||||
|
||||
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
|
||||
x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve2,
|
||||
&lm_control_double, &status);
|
||||
|
||||
// Calculate chi2
|
||||
chi2 = 0;
|
||||
for (ssize_t i = 0; i < y.size(); i++) {
|
||||
chi2 +=
|
||||
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
|
||||
}
|
||||
}
|
||||
|
||||
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
|
||||
NDView<double, 3> y_err, NDView<double, 3> par_out,
|
||||
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
|
||||
int n_threads) {
|
||||
|
||||
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> y_view(&y(row, col, 0), {y.shape(2)});
|
||||
NDView<double, 1> y_err_view(&y_err(row, col, 0),
|
||||
{y_err.shape(2)});
|
||||
NDView<double, 1> par_out_view(&par_out(row, col, 0),
|
||||
{par_out.shape(2)});
|
||||
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
|
||||
{par_err_out.shape(2)});
|
||||
|
||||
fit_scurve2(x, y_view, y_err_view, par_out_view,
|
||||
par_err_out_view, chi2_out(row, col));
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, y.shape(0), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
// FitModel<Model> — method definitions
|
||||
// (constructor, destructor, copy, and all methods that touch Minuit2 state)
|
||||
@@ -756,4 +305,4 @@ AARE_INSTANTIATE_FIT(model::FallingScurve)
|
||||
#undef AARE_INSTANTIATE_FIT
|
||||
// NOLINTEND
|
||||
|
||||
} // namespace aare
|
||||
} // namespace aare
|
||||
|
||||
@@ -0,0 +1,81 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#include "aare/Fit.hpp"
|
||||
#include "aare/FitModel.hpp"
|
||||
#include "aare/Models.hpp"
|
||||
|
||||
#include <catch2/catch_approx.hpp>
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <cmath>
|
||||
|
||||
namespace {
|
||||
|
||||
constexpr ssize_t n_points = 61;
|
||||
|
||||
void fill_gaussian(aare::NDArray<double, 1> &x, aare::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);
|
||||
const double z = (x(i) - mean) / sigma;
|
||||
y(i) = amplitude * std::exp(-0.5 * z * z);
|
||||
}
|
||||
}
|
||||
|
||||
} // 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);
|
||||
|
||||
const aare::FitModel<aare::model::Gaussian> unweighted_model;
|
||||
const auto unweighted =
|
||||
aare::fit_pixel(unweighted_model, x.view(), y.view());
|
||||
|
||||
REQUIRE(unweighted.size() == 4);
|
||||
CHECK(unweighted(0) == Catch::Approx(120.0).epsilon(1e-5));
|
||||
CHECK(unweighted(1) == Catch::Approx(0.8).epsilon(1e-5));
|
||||
CHECK(unweighted(2) == Catch::Approx(1.3).epsilon(1e-5));
|
||||
CHECK(unweighted(3) == Catch::Approx(0.0).margin(1e-4));
|
||||
|
||||
const aare::FitModel<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());
|
||||
|
||||
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));
|
||||
}
|
||||
|
||||
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});
|
||||
fill_gaussian(x, values, 80.0, -0.6, 0.9);
|
||||
|
||||
aare::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)
|
||||
y(row, col, i) = values(i);
|
||||
}
|
||||
}
|
||||
|
||||
aare::NDArray<double, 3> par({2, 2, 3});
|
||||
aare::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);
|
||||
|
||||
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));
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user