mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2026-08-13 03:10:34 +02:00
Merge branch 'main' into dev/fast-pedestal
This commit is contained in:
@@ -47,7 +47,7 @@ jobs:
|
||||
source venv/bin/activate
|
||||
mkdir build
|
||||
cd build
|
||||
cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON
|
||||
cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON -DAARE_WARNINGS_AS_ERRORS=ON
|
||||
make -j 2
|
||||
make docs
|
||||
|
||||
|
||||
@@ -47,7 +47,7 @@ jobs:
|
||||
run: |
|
||||
mkdir build
|
||||
cd build
|
||||
cmake .. -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON -DAARE_TESTS=ON
|
||||
cmake .. -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON -DAARE_TESTS=ON -DAARE_WARNINGS_AS_ERRORS=ON
|
||||
make -j 4
|
||||
make docs
|
||||
|
||||
|
||||
+10
-8
@@ -65,6 +65,9 @@ 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"
|
||||
OFF)
|
||||
|
||||
# Convenience option to use system libraries only (no FetchContent)
|
||||
option(AARE_SYSTEM_LIBRARIES "Use system libraries" OFF)
|
||||
if(AARE_SYSTEM_LIBRARIES)
|
||||
@@ -152,14 +155,10 @@ endif()
|
||||
if(AARE_FETCH_MINUIT2)
|
||||
# We are building Minuit2 from sources.
|
||||
|
||||
# Patch minuit2 to avoid messing with cmake policies
|
||||
set(MINUIT2_PATCH_COMMAND git apply
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/patches/minuit2.patch)
|
||||
FetchContent_Declare(
|
||||
Minuit2
|
||||
GIT_REPOSITORY https://github.com/GooFit/Minuit2.git
|
||||
GIT_TAG master
|
||||
PATCH_COMMAND ${MINUIT2_PATCH_COMMAND}
|
||||
UPDATE_DISCONNECTED 1)
|
||||
# Disable Minuit2 extras we don't need
|
||||
set(minuit2_mpi
|
||||
@@ -172,7 +171,7 @@ if(AARE_FETCH_MINUIT2)
|
||||
OFF
|
||||
CACHE BOOL "")
|
||||
set(MINUIT2_INSTALL
|
||||
ON
|
||||
OFF
|
||||
CACHE BOOL "")
|
||||
FetchContent_MakeAvailable(Minuit2)
|
||||
|
||||
@@ -349,6 +348,10 @@ else()
|
||||
# builds
|
||||
)
|
||||
|
||||
if(AARE_WARNINGS_AS_ERRORS)
|
||||
target_compile_options(aare_compiler_flags INTERFACE -Werror)
|
||||
endif()
|
||||
|
||||
endif() # GCC/Clang specific
|
||||
|
||||
if(AARE_PYTHON_BINDINGS)
|
||||
@@ -388,7 +391,6 @@ set(PUBLICHEADERS
|
||||
include/aare/Dtype.hpp
|
||||
include/aare/File.hpp
|
||||
include/aare/Fit.hpp
|
||||
include/aare/Chi2.hpp
|
||||
include/aare/FitModel.hpp
|
||||
include/aare/Models.hpp
|
||||
include/aare/FileInterface.hpp
|
||||
@@ -452,8 +454,8 @@ target_link_libraries(
|
||||
aare_core
|
||||
PUBLIC fmt::fmt nlohmann_json::nlohmann_json ${STD_FS_LIB} # from
|
||||
# helpers.cmake
|
||||
aare::Minuit2
|
||||
PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:lmfit>)
|
||||
PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:lmfit>
|
||||
$<BUILD_INTERFACE:aare::Minuit2>)
|
||||
|
||||
target_include_directories(
|
||||
aare_core SYSTEM
|
||||
|
||||
+1
-1
@@ -1,6 +1,6 @@
|
||||
# Release notes
|
||||
|
||||
## HEAD
|
||||
## 2026.7.2
|
||||
|
||||
|
||||
### New Features:
|
||||
|
||||
@@ -1,5 +1,4 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#include "aare/Chi2.hpp"
|
||||
#include "aare/Fit.hpp"
|
||||
#include "aare/FitModel.hpp"
|
||||
#include "aare/Models.hpp"
|
||||
@@ -108,9 +107,7 @@ static void BM_FitGausMinuitGrad(benchmark::State &state) {
|
||||
|
||||
aare::NDArray<double, 1> result;
|
||||
for (auto _ : state) {
|
||||
result =
|
||||
aare::fit_pixel<aare::model::Gaussian, aare::func::Chi2Gaussian>(
|
||||
model, xv, yv);
|
||||
result = aare::fit_pixel<aare::model::Gaussian>(model, xv, yv);
|
||||
benchmark::DoNotOptimize(result.data());
|
||||
}
|
||||
|
||||
@@ -132,9 +129,7 @@ static void BM_FitGausMinuitGradHesse(benchmark::State &state) {
|
||||
|
||||
aare::NDArray<double, 1> result;
|
||||
for (auto _ : state) {
|
||||
result =
|
||||
aare::fit_pixel<aare::model::Gaussian, aare::func::Chi2Gaussian>(
|
||||
model, xv, yv, ev);
|
||||
result = aare::fit_pixel<aare::model::Gaussian>(model, xv, yv, ev);
|
||||
benchmark::DoNotOptimize(result.data());
|
||||
}
|
||||
|
||||
|
||||
+19
-190
@@ -2,21 +2,13 @@
|
||||
#pragma once
|
||||
|
||||
#include <cmath>
|
||||
#include <fmt/core.h>
|
||||
#include <vector>
|
||||
|
||||
#include "aare/Chi2.hpp"
|
||||
#include "aare/FitModel.hpp"
|
||||
#include "aare/NDArray.hpp"
|
||||
#include "aare/utils/par.hpp"
|
||||
#include "aare/utils/task.hpp"
|
||||
|
||||
#include "Minuit2/FunctionMinimum.h"
|
||||
#include "Minuit2/MnHesse.h"
|
||||
#include "Minuit2/MnMigrad.h"
|
||||
#include "Minuit2/MnPrint.h"
|
||||
#include "Minuit2/MnUserParameters.h"
|
||||
|
||||
namespace aare {
|
||||
|
||||
namespace func {
|
||||
@@ -49,7 +41,6 @@ NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
|
||||
* @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);
|
||||
|
||||
@@ -115,138 +106,43 @@ void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
|
||||
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
|
||||
int n_threads);
|
||||
|
||||
// Minuit2 fit_pixel / fit_3d object based API
|
||||
|
||||
// _____________________________________________________________________
|
||||
//
|
||||
// fit_pixel — single-pixel minimisation
|
||||
// _____________________________________________________________________
|
||||
// ---------------------------------------------------------------------------
|
||||
// Minuit2-based pixel fitting.
|
||||
// Template bodies and explicit instantiations live in src/Fit.cpp.
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
/**
|
||||
* @brief Fit a single pixel's data using Minuit2.
|
||||
*
|
||||
* The caller provides a thread-local clone of MnUserParameters so that
|
||||
* no heap allocation happens here (only SetValue/SetError stores).
|
||||
*
|
||||
* User-precedence rules:
|
||||
* - Fixed parameters: untouched (value and fixed flag preserved from clone).
|
||||
* - Fixed parameters: untouched (value and fixed flag preserved from model).
|
||||
* - User-set start: value preserved, step size auto-filled.
|
||||
* - Neither: both value and step size auto-filled from data.
|
||||
*
|
||||
* @tparam Model Model struct (Gaussian, RisingScurve, …).
|
||||
* @tparam FCN Chi2 functor type (Chi2Model1D or Chi2Model1DGrad
|
||||
* instantiation).
|
||||
*
|
||||
* @param model The FitModel configuration (read-only).
|
||||
* @param upar_local Thread-local clone of model.upar(). Modified in place.
|
||||
* @param x Scan points (shared across all pixels).
|
||||
* @param y Measured values for this pixel.
|
||||
* @param y_err Per-point uncertainties (empty view -> unweighted fit).
|
||||
* @param model The FitModel configuration (read-only).
|
||||
* @param x Scan points.
|
||||
* @param y Measured values for this pixel.
|
||||
* @param y_err Per-point uncertainties (empty view -> unweighted fit).
|
||||
*
|
||||
* @return NDArray<double,1> of size:
|
||||
* - compute_errors: [p0..pN, err0..errN, chi2] -> 2*npar + 1
|
||||
* - otherwise: [p0..pN, chi2] -> npar + 1
|
||||
*/
|
||||
template <typename Model, typename FCN>
|
||||
NDArray<double, 1> fit_pixel(const FitModel<Model> &model,
|
||||
ROOT::Minuit2::MnUserParameters &upar_local,
|
||||
NDView<double, 1> x, NDView<double, 1> y,
|
||||
NDView<double, 1> y_err) {
|
||||
|
||||
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);
|
||||
|
||||
// ──── automatic parameter estimation ─────────────
|
||||
auto start = Model::estimate_par(x, y);
|
||||
|
||||
// dead / degenerate pixel guard
|
||||
if (!Model::is_valid(std::vector<double>(start.begin(), start.end()))) {
|
||||
return NDArray<double, 1>({result_size}, 0.0);
|
||||
}
|
||||
|
||||
// ──── data-range statistics for step sizes ─────────────
|
||||
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);
|
||||
|
||||
// ── apply auto-estimates respecting user precedence ─────────────
|
||||
for (std::size_t i = 0; i < npar; ++i) {
|
||||
// fixed: do not touch at all
|
||||
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]);
|
||||
}
|
||||
|
||||
// ──── build functor ────────
|
||||
auto chi2 = (y_err.size() > 0) ? FCN(x, y, y_err) : FCN(x, y);
|
||||
|
||||
// ──── run minimizer ────────
|
||||
ROOT::Minuit2::MnMigrad migrad(chi2, upar_local, model.strategy());
|
||||
ROOT::Minuit2::FunctionMinimum min =
|
||||
migrad(model.max_calls(), model.tolerance());
|
||||
|
||||
if (!min.IsValid())
|
||||
return NDArray<double, 1>({result_size}, 0.0);
|
||||
|
||||
// ──── pack results ────────
|
||||
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;
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
// ── self-contained for 1D / standalone use ─────────
|
||||
template <typename Model, typename FCN>
|
||||
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) {
|
||||
auto upar_local = model.upar();
|
||||
return fit_pixel<Model, FCN>(model, upar_local, x, y, y_err);
|
||||
}
|
||||
NDView<double, 1> y, NDView<double, 1> y_err);
|
||||
|
||||
// Overload: uncertainties not provided
|
||||
template <typename Model, typename FCN>
|
||||
template <typename Model>
|
||||
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
|
||||
NDView<double, 1> y) {
|
||||
auto upar_local = model.upar();
|
||||
return fit_pixel<Model, FCN>(model, upar_local, x, y, NDView<double, 1>{});
|
||||
}
|
||||
NDView<double, 1> y);
|
||||
|
||||
// _____________________________________________________________________
|
||||
//
|
||||
// fit_3d — row-parallel fitting over (rows, cols) pixel grid
|
||||
// _____________________________________________________________________
|
||||
/**
|
||||
* @brief Fit all pixels in a 3D data cube (rows x cols x n_scan).
|
||||
*
|
||||
* @tparam Model Model struct.
|
||||
* @tparam FCN Chi2 functor type.
|
||||
*
|
||||
* @param model Fit configuration shared by all pixels.
|
||||
* @param x Scan points, shape `(n_scan)`.
|
||||
@@ -259,78 +155,11 @@ NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
|
||||
* @param chi2_out Output chi-squared / objective values, shape `(rows,
|
||||
* cols)`.
|
||||
* @param n_threads Number of threads used to split rows.
|
||||
*
|
||||
*/
|
||||
template <typename Model, typename FCN>
|
||||
void fit_3d(
|
||||
const FitModel<Model> &model, NDView<double, 1> x, // (n_scan)
|
||||
NDView<double, 3> y, // (rows, cols, n_scan)
|
||||
NDView<double, 3> y_err, // (rows, cols, n_scan) or empty for unweighted fit
|
||||
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;
|
||||
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);
|
||||
|
||||
// ──── checks ───────
|
||||
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 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].");
|
||||
}
|
||||
|
||||
// ──── parallel dispatch ───────
|
||||
auto process = [&](ssize_t first_row, ssize_t last_row) {
|
||||
// one clone per thread
|
||||
auto upar_local = model.upar();
|
||||
|
||||
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, FCN>(model, upar_local, x, values, errors);
|
||||
|
||||
for (std::size_t k = 0; k < npar; ++k) {
|
||||
par_out(row, col, k) = res(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);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, static_cast<int>(y.shape(0)), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
} // namespace aare
|
||||
} // namespace aare
|
||||
|
||||
+45
-80
@@ -2,16 +2,21 @@
|
||||
#pragma once
|
||||
|
||||
#include "aare/Models.hpp"
|
||||
#include <type_traits>
|
||||
|
||||
#include "Minuit2/MnStrategy.h"
|
||||
#include "Minuit2/MnUserParameters.h"
|
||||
#include <array>
|
||||
#include <memory>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
namespace aare {
|
||||
|
||||
template <typename Model> class FitModel {
|
||||
ROOT::Minuit2::MnUserParameters upar_;
|
||||
ROOT::Minuit2::MnStrategy strategy_;
|
||||
|
||||
// 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_;
|
||||
@@ -20,14 +25,7 @@ template <typename Model> class FitModel {
|
||||
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 {
|
||||
for (std::size_t i = 0; i < npar; ++i) {
|
||||
if (upar_.Name(i) == name)
|
||||
return static_cast<unsigned int>(i);
|
||||
}
|
||||
throw std::runtime_error("FitModel: unknown parameter name '" + name +
|
||||
"'");
|
||||
}
|
||||
unsigned int checked_index(const std::string &name) const;
|
||||
|
||||
public:
|
||||
static constexpr std::size_t npar = Model::npar;
|
||||
@@ -41,97 +39,64 @@ template <typename Model> class FitModel {
|
||||
* @param tolerance Minuit2 EDM tolerance.
|
||||
* @param compute_errors If true, run MnHesse after minimisation.
|
||||
*/
|
||||
FitModel(unsigned int strategy = 0, unsigned int max_calls = 100,
|
||||
double tolerance = 0.5, bool compute_errors = false)
|
||||
: 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);
|
||||
explicit FitModel(unsigned int strategy = 0, unsigned int max_calls = 100,
|
||||
double tolerance = 0.5, bool compute_errors = false);
|
||||
|
||||
// Add parameters and valid bounds
|
||||
if (has_lo && has_hi) {
|
||||
upar_.Add(pi.name, 0.0, 1.0, pi.default_lo, pi.default_hi);
|
||||
} else if (has_lo) {
|
||||
upar_.Add(pi.name, 0.0, 1.0, pi.default_lo, 1e6);
|
||||
} else {
|
||||
upar_.Add(pi.name, 0.0, 1.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
// 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.*/
|
||||
void SetParLimits(unsigned int idx, double lo, double hi) {
|
||||
upar_.SetLimits(idx, lo, hi);
|
||||
}
|
||||
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.
|
||||
*/
|
||||
void FixParameter(unsigned int idx, double val) {
|
||||
SetParameter(idx, val);
|
||||
upar_.Fix(idx);
|
||||
user_fixed_[idx] = true;
|
||||
}
|
||||
void FixParameter(unsigned int idx, double val);
|
||||
|
||||
/** @brief Release a previously fixed parameter, re-enabling auto estimates.
|
||||
*/
|
||||
void ReleaseParameter(unsigned int idx) {
|
||||
upar_.Release(idx);
|
||||
user_fixed_[idx] = false;
|
||||
}
|
||||
|
||||
void ReleaseParameter(const std::string &name) {
|
||||
ReleaseParameter(checked_index(name));
|
||||
}
|
||||
void ReleaseParameter(unsigned int idx);
|
||||
void ReleaseParameter(const std::string &name);
|
||||
|
||||
/** @brief Set an explicit starting value for parameter idx.*/
|
||||
void SetParameter(unsigned int idx, double val) {
|
||||
upar_.SetValue(idx, val);
|
||||
user_start_[idx] = true;
|
||||
}
|
||||
void SetParameter(unsigned int idx, double val);
|
||||
void SetParameter(const std::string &name, double val);
|
||||
void FixParameter(const std::string &name, double val);
|
||||
void SetParLimits(const std::string &name, double lo, double hi);
|
||||
std::string GetParName(unsigned int idx) const;
|
||||
std::vector<std::string> GetParNames() const;
|
||||
|
||||
void SetParameter(const std::string &name, double val) {
|
||||
// go through index to maintain user_start_ bookkeeping
|
||||
SetParameter(checked_index(name), val);
|
||||
}
|
||||
|
||||
void FixParameter(const std::string &name, double val) {
|
||||
// go through index to maintain user_fixed_ bookkeeping
|
||||
FixParameter(checked_index(name), val);
|
||||
}
|
||||
|
||||
void SetParLimits(const std::string &name, double lo, double hi) {
|
||||
SetParLimits(checked_index(name), lo, hi);
|
||||
}
|
||||
|
||||
std::string GetParName(unsigned int idx) const {
|
||||
return upar_.GetName(idx);
|
||||
}
|
||||
|
||||
std::vector<std::string> GetParNames() const {
|
||||
std::vector<std::string> names;
|
||||
for (std::size_t i = 0; i < npar; ++i)
|
||||
names.push_back(GetParName(i));
|
||||
return names;
|
||||
}
|
||||
static constexpr std::size_t GetNpar() noexcept { return npar; }
|
||||
|
||||
void SetMaxCalls(unsigned int n) { max_calls_ = n; }
|
||||
void SetTolerance(double t) { tolerance_ = t; }
|
||||
void SetComputeErrors(bool b) { compute_errors_ = b; }
|
||||
|
||||
// accessors
|
||||
const ROOT::Minuit2::MnUserParameters &upar() const { return upar_; }
|
||||
const ROOT::Minuit2::MnStrategy &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]; }
|
||||
|
||||
// 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(); }
|
||||
};
|
||||
|
||||
// Suppress implicit instantiation for all supported model types.
|
||||
// Definitions live in src/Fit.cpp.
|
||||
extern template class FitModel<model::Gaussian>;
|
||||
extern template class FitModel<model::GaussianErfcPlateau>;
|
||||
extern template class FitModel<model::GaussianChargeSharing>;
|
||||
extern template class FitModel<model::GaussianChargeSharingKb>;
|
||||
extern template class FitModel<model::Pol1>;
|
||||
extern template class FitModel<model::Pol2>;
|
||||
extern template class FitModel<model::RisingScurve>;
|
||||
extern template class FitModel<model::FallingScurve>;
|
||||
|
||||
} // namespace aare
|
||||
|
||||
@@ -219,6 +219,17 @@ struct Moench04 {
|
||||
constexpr static size_t superColumnWidth = 25;
|
||||
};
|
||||
|
||||
/// @brief Chip specifications for Moench05
|
||||
struct Moench05 {
|
||||
constexpr static size_t nRows = 160;
|
||||
constexpr static size_t nCols =
|
||||
150; // TODO: should probably be seperated for each adc
|
||||
|
||||
/// @brief used ADCs for moench given in relative numbers to the default
|
||||
/// enabled blocks of 4 adcs in absolute adcs 9, 13, 1 are used
|
||||
constexpr static std::array<int, 3> adcNumbers = {5, 9, 1};
|
||||
};
|
||||
|
||||
enum ReadoutMode : uint8_t {
|
||||
ANALOG_ONLY = 0,
|
||||
DIGITAL_ONLY = 1,
|
||||
|
||||
@@ -1,39 +0,0 @@
|
||||
diff --git a/CMakeLists.txt b/CMakeLists.txt
|
||||
index 4e72289..33dd58a 100644
|
||||
--- a/CMakeLists.txt
|
||||
+++ b/CMakeLists.txt
|
||||
@@ -4,7 +4,7 @@
|
||||
# For the licensing terms see $ROOTSYS/LICENSE.
|
||||
# For the list of contributors see $ROOTSYS/README/CREDITS.
|
||||
|
||||
-cmake_minimum_required(VERSION 3.1)
|
||||
+cmake_minimum_required(VERSION 3.15)
|
||||
|
||||
if(NOT CMAKE_PROJECT_NAME STREQUAL ROOT)
|
||||
project(Minuit2 LANGUAGES CXX)
|
||||
diff --git a/StandAlone.cmake b/StandAlone.cmake
|
||||
index 0a662ba..b7a05ca 100644
|
||||
--- a/StandAlone.cmake
|
||||
+++ b/StandAlone.cmake
|
||||
@@ -1,12 +1,12 @@
|
||||
-cmake_minimum_required(VERSION 3.1)
|
||||
-
|
||||
-# Tested with and supporting policies up to the following CMake version.
|
||||
-# Not using ... syntax due to parser bug in MSVC's built-in CMake server mode.
|
||||
-if(${CMAKE_VERSION} VERSION_LESS 3.12)
|
||||
- cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
|
||||
-else()
|
||||
- cmake_policy(VERSION 3.12)
|
||||
-endif()
|
||||
+cmake_minimum_required(VERSION 3.15)
|
||||
+
|
||||
+# # Tested with and supporting policies up to the following CMake version.
|
||||
+# # Not using ... syntax due to parser bug in MSVC's built-in CMake server mode.
|
||||
+# if(${CMAKE_VERSION} VERSION_LESS 3.12)
|
||||
+# cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
|
||||
+# else()
|
||||
+# cmake_policy(VERSION 3.12)
|
||||
+# endif()
|
||||
|
||||
include(FeatureSummary)
|
||||
include(CMakeDependentOption)
|
||||
@@ -17,6 +17,15 @@ class AdcSar05060708Transform64to16:
|
||||
return _aare.adc_sar_05_06_07_08decode64to16(data)
|
||||
|
||||
class Moench05Transform:
|
||||
"""
|
||||
Transforms Moench05 chip data from a buffer of bytes (uint8_t)
|
||||
to a numpy array of uint16. Assumes data taken with analog samples and assumes adc 1, 9, 13 are enabled.
|
||||
(e.g. for 10g mode adc 0,1,2,3 and 8,9,10,11 and 12,13,14,15 are enabled but only adc 1,9,13 contain relevant data)
|
||||
|
||||
.. note::
|
||||
A moench05 chip has 160 rows and 50 cols per adc and has dynamic range 16 bit. Each adc sample is encoded in 16 bits.
|
||||
The transformation thus requires 160*50*16/16 = 8000 analog samples per adc.
|
||||
"""
|
||||
#Could be moved to C++ without changing the interface
|
||||
def __init__(self):
|
||||
self.pixel_map = _aare.GenerateMoench05PixelMap()
|
||||
|
||||
@@ -22,4 +22,9 @@ void define_defs_bindings(py::module &m) {
|
||||
moench04.attr("nPixelsPerSuperColumn") = Moench04::nPixelsPerSuperColumn;
|
||||
moench04.attr("superColumnWidth") = Moench04::superColumnWidth;
|
||||
moench04.attr("adcNumbers") = Moench04::adcNumbers;
|
||||
|
||||
auto moench05 = py::class_<Moench05>(m, "Moench05");
|
||||
moench05.attr("nRows") = Moench05::nRows;
|
||||
moench05.attr("nCols") = Moench05::nCols;
|
||||
moench05.attr("adcNumbers") = Moench05::adcNumbers;
|
||||
}
|
||||
|
||||
+24
-34
@@ -5,7 +5,6 @@
|
||||
#include <pybind11/stl.h>
|
||||
#include <pybind11/stl_bind.h>
|
||||
|
||||
#include "aare/Chi2.hpp"
|
||||
#include "aare/Fit.hpp"
|
||||
#include "aare/FitModel.hpp"
|
||||
#include "aare/Models.hpp"
|
||||
@@ -13,7 +12,7 @@
|
||||
namespace py = pybind11;
|
||||
using namespace pybind11::literals;
|
||||
|
||||
template <typename Model, typename FCN>
|
||||
template <typename Model>
|
||||
py::object
|
||||
fit_dispatch(const aare::FitModel<Model> &model,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
@@ -22,7 +21,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
|
||||
|
||||
template <typename Model> void bind_fit_model(py::module &m, const char *name) {
|
||||
using FM = aare::FitModel<Model>;
|
||||
using FCN = aare::func::Chi2Model1DGrad<Model>;
|
||||
py::class_<FM>(m, name)
|
||||
.def(py::init<unsigned int, unsigned int, double, bool>(),
|
||||
py::arg("strategy") = 0, py::arg("max_calls") = 100,
|
||||
@@ -85,8 +83,7 @@ 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, FCN>(self, x, y, y_err_obj,
|
||||
n_threads);
|
||||
return fit_dispatch<Model>(self, x, y, y_err_obj, n_threads);
|
||||
},
|
||||
R"doc(
|
||||
Fit this model to 1D or 3D data using Minuit2.
|
||||
@@ -145,7 +142,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, typename FCN>
|
||||
template <typename Model>
|
||||
py::object
|
||||
fit_dispatch(const aare::FitModel<Model> &model,
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> x,
|
||||
@@ -175,9 +172,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, FCN>(model, x_view, y_view, y_view_err,
|
||||
par_out->view(), err_out->view(),
|
||||
chi2_out->view(), n_threads);
|
||||
aare::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),
|
||||
@@ -193,9 +190,9 @@ fit_dispatch(const aare::FitModel<Model> &model,
|
||||
NDView<double, 3> dummy_err{};
|
||||
NDView<double, 3> dummy_err_out{};
|
||||
|
||||
aare::fit_3d<Model, FCN>(model, x_view, y_view, dummy_err,
|
||||
par_out->view(), dummy_err_out,
|
||||
chi2_out->view(), n_threads);
|
||||
aare::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));
|
||||
@@ -217,10 +214,9 @@ fit_dispatch(const aare::FitModel<Model> &model,
|
||||
}
|
||||
|
||||
auto y_view_err = make_view_1d(y_err);
|
||||
result =
|
||||
aare::fit_pixel<Model, FCN>(model, x_view, y_view, y_view_err);
|
||||
result = aare::fit_pixel<Model>(model, x_view, y_view, y_view_err);
|
||||
} else {
|
||||
result = aare::fit_pixel<Model, FCN>(model, x_view, y_view);
|
||||
result = aare::fit_pixel<Model>(model, x_view, y_view);
|
||||
}
|
||||
|
||||
return pack_1d_result_dict<Model>(result, model.compute_errors());
|
||||
@@ -696,29 +692,26 @@ void define_fit_bindings(py::module &m) {
|
||||
py::array_t<double, py::array::c_style | py::array::forcecast> y,
|
||||
py::object y_err_obj, int n_threads) -> py::object {
|
||||
using namespace aare::model;
|
||||
using namespace aare::func;
|
||||
|
||||
// ── Polynomial of degree 1 ───────
|
||||
if (py::isinstance<aare::FitModel<Pol1>>(model_obj)) {
|
||||
const auto &mdl =
|
||||
model_obj.cast<const aare::FitModel<Pol1> &>();
|
||||
return fit_dispatch<Pol1, Chi2Pol1>(mdl, x, y, y_err_obj,
|
||||
n_threads);
|
||||
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, Chi2Pol2>(mdl, x, y, y_err_obj,
|
||||
n_threads);
|
||||
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, Chi2Gaussian>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
return fit_dispatch<Gaussian>(mdl, x, y, y_err_obj, n_threads);
|
||||
}
|
||||
|
||||
// ── GaussianErfcPlateau ───────
|
||||
@@ -727,9 +720,8 @@ void define_fit_bindings(py::module &m) {
|
||||
const auto &mdl =
|
||||
model_obj
|
||||
.cast<const aare::FitModel<GaussianErfcPlateau> &>();
|
||||
return fit_dispatch<GaussianErfcPlateau,
|
||||
Chi2GaussianErfcPlateau>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
return fit_dispatch<GaussianErfcPlateau>(mdl, x, y, y_err_obj,
|
||||
n_threads);
|
||||
}
|
||||
|
||||
// ── GaussianChargeSharing ───────
|
||||
@@ -738,9 +730,8 @@ void define_fit_bindings(py::module &m) {
|
||||
const auto &mdl =
|
||||
model_obj
|
||||
.cast<const aare::FitModel<GaussianChargeSharing> &>();
|
||||
return fit_dispatch<GaussianChargeSharing,
|
||||
Chi2GaussianChargeSharing>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
return fit_dispatch<GaussianChargeSharing>(mdl, x, y, y_err_obj,
|
||||
n_threads);
|
||||
}
|
||||
|
||||
// ── GaussianChargeSharingKb ───────
|
||||
@@ -748,8 +739,7 @@ void define_fit_bindings(py::module &m) {
|
||||
model_obj)) {
|
||||
const auto &mdl = model_obj.cast<
|
||||
const aare::FitModel<GaussianChargeSharingKb> &>();
|
||||
return fit_dispatch<GaussianChargeSharingKb,
|
||||
Chi2GaussianChargeSharingKb>(
|
||||
return fit_dispatch<GaussianChargeSharingKb>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
}
|
||||
|
||||
@@ -757,16 +747,16 @@ void define_fit_bindings(py::module &m) {
|
||||
if (py::isinstance<aare::FitModel<RisingScurve>>(model_obj)) {
|
||||
const auto &mdl =
|
||||
model_obj.cast<const aare::FitModel<RisingScurve> &>();
|
||||
return fit_dispatch<RisingScurve, Chi2RisingScurve>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
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, Chi2FallingScurve>(
|
||||
mdl, x, y, y_err_obj, n_threads);
|
||||
return fit_dispatch<FallingScurve>(mdl, x, y, y_err_obj,
|
||||
n_threads);
|
||||
}
|
||||
|
||||
throw std::runtime_error(
|
||||
|
||||
File diff suppressed because one or more lines are too long
@@ -396,7 +396,7 @@
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.12.5"
|
||||
"version": "3.11.15"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -135,20 +135,6 @@ class Chi2Model1DGrad : public ROOT::Minuit2::FCNGradientBase {
|
||||
bool weighted_;
|
||||
};
|
||||
|
||||
// ── Convenient aliases ──────────────────────────────────────────────
|
||||
|
||||
using Chi2Gaussian = Chi2Model1DGrad<aare::model::Gaussian>;
|
||||
using Chi2GaussianErfcPlateau =
|
||||
Chi2Model1DGrad<aare::model::GaussianErfcPlateau>;
|
||||
using Chi2GaussianChargeSharing =
|
||||
Chi2Model1DGrad<aare::model::GaussianChargeSharing>;
|
||||
using Chi2GaussianChargeSharingKb =
|
||||
Chi2Model1DGrad<aare::model::GaussianChargeSharingKb>;
|
||||
using Chi2RisingScurve = Chi2Model1DGrad<aare::model::RisingScurve>;
|
||||
using Chi2FallingScurve = Chi2Model1DGrad<aare::model::FallingScurve>;
|
||||
using Chi2Pol1 = Chi2Model1DGrad<aare::model::Pol1>;
|
||||
using Chi2Pol2 = Chi2Model1DGrad<aare::model::Pol2>;
|
||||
|
||||
} // namespace func
|
||||
|
||||
} // namespace aare
|
||||
+299
-1
@@ -1,12 +1,21 @@
|
||||
// SPDX-License-Identifier: MPL-2.0
|
||||
#include "aare/Fit.hpp"
|
||||
#include "aare/Chi2.hpp"
|
||||
#include "Chi2.hpp"
|
||||
#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"
|
||||
#include "aare/utils/par.hpp"
|
||||
#include "aare/utils/task.hpp"
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <lmcurve2.h>
|
||||
#include <lmfit.hpp>
|
||||
#include <memory>
|
||||
#include <stdexcept>
|
||||
#include <thread>
|
||||
#include <type_traits>
|
||||
|
||||
@@ -458,4 +467,293 @@ void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
// 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) {
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
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)
|
||||
return static_cast<unsigned int>(i);
|
||||
}
|
||||
throw std::runtime_error("FitModel: unknown parameter name '" + name + "'");
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
void FitModel<Model>::SetParLimits(unsigned int idx, double lo, double hi) {
|
||||
impl_->upar.SetLimits(idx, lo, 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);
|
||||
user_fixed_[idx] = false;
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
void FitModel<Model>::ReleaseParameter(const std::string &name) {
|
||||
ReleaseParameter(checked_index(name));
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
void FitModel<Model>::SetParameter(unsigned int idx, double val) {
|
||||
impl_->upar.SetValue(idx, val);
|
||||
user_start_[idx] = true;
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
void FitModel<Model>::SetParameter(const std::string &name, double val) {
|
||||
SetParameter(checked_index(name), val);
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
void FitModel<Model>::FixParameter(const std::string &name, double val) {
|
||||
FixParameter(checked_index(name), val);
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
void FitModel<Model>::SetParLimits(const std::string &name, double lo,
|
||||
double hi) {
|
||||
SetParLimits(checked_index(name), lo, hi);
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
std::string FitModel<Model>::GetParName(unsigned int idx) const {
|
||||
return impl_->upar.GetName(idx);
|
||||
}
|
||||
|
||||
template <typename Model>
|
||||
std::vector<std::string> FitModel<Model>::GetParNames() const {
|
||||
std::vector<std::string> names;
|
||||
for (std::size_t i = 0; i < npar; ++i)
|
||||
names.push_back(GetParName(i));
|
||||
return names;
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
// fit_pixel / fit_3d — Minuit2 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);
|
||||
|
||||
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);
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
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) {
|
||||
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].");
|
||||
|
||||
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) {
|
||||
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, k) = res(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);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
auto tasks = split_task(0, static_cast<int>(y.shape(0)), n_threads);
|
||||
RunInParallel(process, tasks);
|
||||
}
|
||||
|
||||
// ============================================================================
|
||||
// Explicit instantiations for all supported model types
|
||||
// ============================================================================
|
||||
|
||||
// NOLINTBEGIN
|
||||
#define AARE_INSTANTIATE_FIT(Model) \
|
||||
template class FitModel<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_FIT(model::Gaussian)
|
||||
AARE_INSTANTIATE_FIT(model::GaussianErfcPlateau)
|
||||
AARE_INSTANTIATE_FIT(model::GaussianChargeSharing)
|
||||
AARE_INSTANTIATE_FIT(model::GaussianChargeSharingKb)
|
||||
AARE_INSTANTIATE_FIT(model::Pol1)
|
||||
AARE_INSTANTIATE_FIT(model::Pol2)
|
||||
AARE_INSTANTIATE_FIT(model::RisingScurve)
|
||||
AARE_INSTANTIATE_FIT(model::FallingScurve)
|
||||
|
||||
#undef AARE_INSTANTIATE_FIT
|
||||
// NOLINTEND
|
||||
|
||||
} // namespace aare
|
||||
+14
-8
@@ -56,16 +56,22 @@ NDArray<ssize_t, 2> GenerateMoench04AnalogPixelMap() {
|
||||
}
|
||||
|
||||
NDArray<ssize_t, 2> GenerateMoench05PixelMap() {
|
||||
std::array<int, 3> adc_numbers = {5, 9, 1};
|
||||
NDArray<ssize_t, 2> order_map({160, 150});
|
||||
constexpr size_t num_adcs = 3; // num adcs with relevant data
|
||||
constexpr std::array<int, num_adcs> adc_numbers = Moench05::adcNumbers;
|
||||
NDArray<ssize_t, 2> order_map({Moench05::nRows, Moench05::nCols});
|
||||
constexpr size_t n_cols = Moench05::nCols / adc_numbers.size();
|
||||
int n_pixel = 0;
|
||||
for (int row = 0; row < 160; row++) {
|
||||
for (int i_col = 0; i_col < 50; i_col++) {
|
||||
n_pixel = row * 50 + i_col;
|
||||
for (int i_sc = 0; i_sc < 3; i_sc++) {
|
||||
int col = 50 * i_sc + i_col;
|
||||
|
||||
constexpr size_t num_adcs_enabled =
|
||||
12; // number of adcs enabled -> for 10g adcs are enabled in blocks of 4
|
||||
|
||||
for (size_t row = 0; row < Moench05::nRows; row++) {
|
||||
for (size_t i_col = 0; i_col < n_cols; i_col++) {
|
||||
n_pixel = row * n_cols + i_col;
|
||||
for (size_t i_sc = 0; i_sc < num_adcs; i_sc++) {
|
||||
size_t col = n_cols * i_sc + i_col;
|
||||
int adc_nr = adc_numbers[i_sc];
|
||||
int i_analog = n_pixel * 12 + adc_nr;
|
||||
int i_analog = n_pixel * num_adcs_enabled + adc_nr;
|
||||
|
||||
// analog_frame[row * 150 + col] = analog_data[i_analog] &
|
||||
// 0x3FFF;
|
||||
|
||||
@@ -58,7 +58,7 @@ PedestalTrackingPixelHistogram::PedestalTrackingPixelHistogram(
|
||||
partial_std_.reserve(n_threads_);
|
||||
for (int i = 0; i < n_threads_; ++i) {
|
||||
const auto local_rows = row_count(i);
|
||||
partial_hists_.emplace_back(local_rows, cols, n_bins, xmin, xmax);
|
||||
partial_hists_.emplace_back(local_rows, cols, n_bins, xmin_, xmax_);
|
||||
partial_pedestals_.emplace_back(static_cast<uint32_t>(local_rows),
|
||||
static_cast<uint32_t>(cols));
|
||||
partial_std_.emplace_back(NDArray<AxisType, 2>(
|
||||
|
||||
@@ -1,6 +1,7 @@
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <catch2/generators/catch_generators.hpp>
|
||||
|
||||
#include "aare/logger.hpp"
|
||||
#include <chrono>
|
||||
#include <cstdint>
|
||||
#include <random>
|
||||
@@ -25,14 +26,15 @@ namespace {
|
||||
} // namespace
|
||||
|
||||
TEST_CASE("Fill one pixel of a 5x10 histogram") {
|
||||
PixelHistogram hist(5, 10, 20, 0.0, 10.0);
|
||||
PixelHistogram hist(5, 10, 20, 0.0f, 10.0f);
|
||||
NDArray<float, 2> image(
|
||||
{5, 10}, -1.0); // Need to fill with -1 to not generate counts
|
||||
{5, 10}, -1.0f); // Need to fill with -1 to not generate counts
|
||||
|
||||
image(2, 3) = 5.7; // This should go into bin 11 (since bins are [0-0.5),
|
||||
// [0.5-1.0), ..., [9.5-10.0))
|
||||
image(2, 3) = 5.7f; // This should go into bin 11 (since bins are [0-0.5),
|
||||
// [0.5-1.0), ..., [9.5-10.0))
|
||||
|
||||
// fill_blocking(hist, image.view());
|
||||
|
||||
hist.fill_async(NDArray<float, 2>(image));
|
||||
hist.flush(); // Wait for the async fill to complete before we check the
|
||||
// results
|
||||
@@ -58,14 +60,14 @@ TEST_CASE("Fill one pixel of a 5x10 histogram") {
|
||||
}
|
||||
|
||||
TEST_CASE("Fill pixels with uneven partial histogram row slices") {
|
||||
PixelHistogram hist(5, 4, 10, 0.0, 10.0, 3);
|
||||
NDArray<float, 2> image({5, 4}, -1.0);
|
||||
PixelHistogram hist(5, 4, 10, 0.0f, 10.0f, 3);
|
||||
NDArray<float, 2> image({5, 4}, -1.0f);
|
||||
|
||||
image(0, 0) = 0.2;
|
||||
image(1, 1) = 1.2;
|
||||
image(2, 2) = 2.2;
|
||||
image(3, 3) = 3.2;
|
||||
image(4, 0) = 4.2;
|
||||
image(0, 0) = 0.2f;
|
||||
image(1, 1) = 1.2f;
|
||||
image(2, 2) = 2.2f;
|
||||
image(3, 3) = 3.2f;
|
||||
image(4, 0) = 4.2f;
|
||||
|
||||
hist.fill_async(NDArray<float, 2>(image));
|
||||
hist.flush();
|
||||
@@ -228,7 +230,7 @@ TEST_CASE("Random fills match a reference implementation") {
|
||||
}
|
||||
|
||||
TEST_CASE("fill_async with mismatched shape throws") {
|
||||
PixelHistogram hist(8, 8, 16, 0.0, 1.0, 2);
|
||||
PixelHistogram hist(8, 8, 16, 0.0f, 1.0f, 2);
|
||||
NDArray<float, 2> bad({4, 4}, 0.0f);
|
||||
CHECK_THROWS_AS(hist.fill_async(std::move(bad)), std::invalid_argument);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user