detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2025-04-20 05:40:03 +02:00
Code Issues Actions Packages Releases Activity

Added chi2 to fit results (#131)

Browse Source 
- fit_gaus and fit_pol1 now return a dict
- calculate chi2 after fit
- cleaned up code
This commit is contained in:
Erik Fröjdh 2025-02-18 21:13:27 +01:00 committed by GitHub
parent 8abfc68138
commit 6a83988485
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
10 changed files with 291 additions and 255 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -85,7 +85,7 @@ if(AARE_FETCH_LMFIT)
GIT_TAG main GIT_TAG main
PATCH_COMMAND ${lmfit_patch} PATCH_COMMAND ${lmfit_patch}
UPDATE_DISCONNECTED 1 UPDATE_DISCONNECTED 1
EXCLUDE_FROM_ALL EXCLUDE_FROM_ALL 1
) )
#Disable what we don't need from lmfit #Disable what we don't need from lmfit
set(BUILD_TESTING OFF CACHE BOOL "") set(BUILD_TESTING OFF CACHE BOOL "")

39
include/aare/Fit.hpp
View File

@ -17,6 +17,13 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
} // namespace func } // namespace func
/**
* @brief Estimate the initial parameters for a Gaussian fit
*/
std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
static constexpr int DEFAULT_NUM_THREADS = 4; static constexpr int DEFAULT_NUM_THREADS = 4;
/** /**
@ -26,14 +33,15 @@ static constexpr int DEFAULT_NUM_THREADS = 4;
*/ */
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y); 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] * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
* @param x x values * @param x x values
* @param y y vales, layout [row, col, values] * @param y y vales, layout [row, col, values]
* @param n_threads number of threads to use * @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); NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS);
/** /**
@ -45,10 +53,12 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y, int n_thre
* @param par_err_out output error parameters * @param par_err_out output error parameters
*/ */
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err, 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); 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] * @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
* [row, col, values]
* @param x x values * @param x x values
* @param y y vales, layout [row, col, values] * @param y y vales, layout [row, col, values]
* @param y_err error in y, layout [row, col, values] * @param y_err error in y, layout [row, col, values]
@ -57,20 +67,21 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
* @param n_threads number of threads to use * @param n_threads number of threads to use
*/ */
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, 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, int n_threads = DEFAULT_NUM_THREADS); 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, 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); 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, void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> y_err, NDView<double, 1> par_out, NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
NDView<double, 1> par_err_out);
//TODO! not sure we need to offer the different version in C++ // TODO! not sure we need to offer the different version in C++
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> y_err, NDView<double, 3> par_out, NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out,
NDView<double, 3> par_err_out, int n_threads = DEFAULT_NUM_THREADS); int n_threads = DEFAULT_NUM_THREADS);
} // namespace aare } // namespace aare

29
include/aare/NDArray.hpp
View File

@ -69,6 +69,11 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::copy(v.begin(), v.end(), begin()); std::copy(v.begin(), v.end(), begin());
} }
template<size_t Size>
NDArray(const std::array<T, Size>& arr) : NDArray<T,1>({Size}) {
std::copy(arr.begin(), arr.end(), begin());
}
// Move constructor // Move constructor
NDArray(NDArray &&other) noexcept NDArray(NDArray &&other) noexcept
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)), : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
@ -105,6 +110,20 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
NDArray &operator-=(const NDArray &other); NDArray &operator-=(const NDArray &other);
NDArray &operator*=(const NDArray &other); NDArray &operator*=(const NDArray &other);
//Write directly to the data array, or create a new one
template<size_t Size>
NDArray<T,1>& operator=(const std::array<T,Size> &other){
if(Size != size_){
delete[] data_;
size_ = Size;
data_ = new T[size_];
}
for (size_t i = 0; i < Size; ++i) {
data_[i] = other[i];
}
return *this;
}
// NDArray& operator/=(const NDArray& other); // NDArray& operator/=(const NDArray& other);
template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) { template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
@ -135,6 +154,11 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
NDArray &operator&=(const T & /*mask*/); NDArray &operator&=(const T & /*mask*/);
void sqrt() { void sqrt() {
for (int i = 0; i < size_; ++i) { for (int i = 0; i < size_; ++i) {
data_[i] = std::sqrt(data_[i]); data_[i] = std::sqrt(data_[i]);
@ -318,6 +342,9 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
return *this; return *this;
} }
template <typename T, int64_t Ndim> template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) { NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
NDArray result = *this; NDArray result = *this;
@ -418,4 +445,6 @@ NDArray<T, Ndim> load(const std::string &pathname,
return img; return img;
} }
} // namespace aare } // namespace aare

11
include/aare/NDView.hpp
View File

@ -1,5 +1,5 @@
#pragma once #pragma once
#include "aare/defs.hpp"
#include "aare/ArrayExpr.hpp" #include "aare/ArrayExpr.hpp"
#include <algorithm> #include <algorithm>
@ -99,6 +99,15 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides<T>()); } NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides<T>()); }
template<size_t Size>
NDView& operator=(const std::array<T, Size> &arr) {
if(size() != arr.size())
throw std::runtime_error(LOCATION + "Array and NDView size mismatch");
std::copy(arr.begin(), arr.end(), begin());
return *this;
}
NDView &operator=(const T val) { NDView &operator=(const T val) {
for (auto it = begin(); it != end(); ++it) for (auto it = begin(); it != end(); ++it)
*it = val; *it = val;

18
include/aare/utils/par.hpp Normal file
View File

@ -0,0 +1,18 @@
#include <thread>
#include <vector>
#include <utility>
namespace aare {
template<typename F>
void RunInParallel(F func, const std::vector<std::pair<int, int>>& tasks) {
// auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(func, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
}
} // namespace aare

3
python/CMakeLists.txt
View File

@ -49,11 +49,10 @@ set(PYTHON_EXAMPLES
examples/fits.py examples/fits.py
) )
# Copy the python examples to the build directory # Copy the python examples to the build directory
foreach(FILE ${PYTHON_EXAMPLES}) foreach(FILE ${PYTHON_EXAMPLES})
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} ) configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} )
message(STATUS "Copying ${FILE} to ${CMAKE_BINARY_DIR}/${FILE}")
endforeach(FILE ${PYTHON_EXAMPLES}) endforeach(FILE ${PYTHON_EXAMPLES})

62
python/examples/play.py
View File

@ -8,61 +8,15 @@ import numpy as np
import boost_histogram as bh import boost_histogram as bh
import time import time
<<<<<<< HEAD import aare
from aare import File, ClusterFinder, VarClusterFinder, ClusterFile, CtbRawFile
from aare import gaus, fit_gaus
base = Path('/mnt/sls_det_storage/moench_data/Julian/MOENCH05/20250113_first_xrays_redo/raw_files/') data = np.random.normal(10, 1, 1000)
cluster_file = Path('/home/l_msdetect/erik/tmp/Cu.clust')
t0 = time.perf_counter() hist = bh.Histogram(bh.axis.Regular(10, 0, 20))
offset= -0.5 hist.fill(data)
hist3d = bh.Histogram(
bh.axis.Regular(160, 0+offset, 160+offset), #x
bh.axis.Regular(150, 0+offset, 150+offset), #y
bh.axis.Regular(200, 0, 6000), #ADU
)
total_clusters = 0
with ClusterFile(cluster_file, chunk_size = 1000) as f:
for i, clusters in enumerate(f):
arr = np.array(clusters)
total_clusters += clusters.size
hist3d.fill(arr['y'],arr['x'], clusters.sum_2x2()) #python talks [row, col] cluster finder [x,y]
=======
from aare import RawFile
f = RawFile('/mnt/sls_det_storage/jungfrau_data1/vadym_tests/jf12_M431/laser_scan/laserScan_pedestal_G0_master_0.json')
print(f'{f.frame_number(1)}')
for i in range(10):
header, img = f.read_frame()
print(header['frameNumber'], img.shape)
>>>>>>> developer
t_elapsed = time.perf_counter()-t0 x = hist.axes[0].centers
print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s') y = hist.values()
y_err = np.sqrt(y)+1
histogram_data = hist3d.counts() res = aare.fit_gaus(x, y, y_err, chi2 = True)
x = hist3d.axes[2].edges[:-1]
y = histogram_data[100,100,:]
xx = np.linspace(x[0], x[-1])
# fig, ax = plt.subplots()
# ax.step(x, y, where = 'post')
y_err = np.sqrt(y)
y_err = np.zeros(y.size)
y_err += 1
# par = fit_gaus2(y,x, y_err)
# ax.plot(xx, gaus(xx,par))
# print(par)
res = fit_gaus(y,x)
res2 = fit_gaus(y,x, y_err)
print(res)
print(res2)

70
python/src/fit.hpp
View File

@ -7,6 +7,7 @@
#include "aare/Fit.hpp" #include "aare/Fit.hpp"
namespace py = pybind11; namespace py = pybind11;
using namespace pybind11::literals;
void define_fit_bindings(py::module &m) { void define_fit_bindings(py::module &m) {
@ -29,7 +30,8 @@ void define_fit_bindings(py::module &m) {
The points at which to evaluate the Gaussian function. The points at which to evaluate the Gaussian function.
par : array_like 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. 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")); )",
py::arg("x"), py::arg("par"));
m.def( m.def(
"pol1", "pol1",
@ -49,7 +51,8 @@ void define_fit_bindings(py::module &m) {
The points at which to evaluate the polynomial function. The points at which to evaluate the polynomial function.
par : array_like par : array_like
The parameters of the polynomial function. The first element is the intercept, and the second element is the slope. 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")); )",
py::arg("x"), py::arg("par"));
m.def( m.def(
"fit_gaus", "fit_gaus",
@ -72,7 +75,7 @@ void define_fit_bindings(py::module &m) {
throw std::runtime_error("Data must be 1D or 3D"); throw std::runtime_error("Data must be 1D or 3D");
} }
}, },
R"( R"(
Fit a 1D Gaussian to data. Fit a 1D Gaussian to data.
Parameters Parameters
@ -90,8 +93,8 @@ n_threads : int, optional
"fit_gaus", "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> 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,
py::array_t<double, py::array::c_style | py::array::forcecast> py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
y_err, int n_threads) { int n_threads) {
if (y.ndim() == 3) { if (y.ndim() == 3) {
// Allocate memory for the output // Allocate memory for the output
// Need to have pointers to allow python to manage // Need to have pointers to allow python to manage
@ -99,15 +102,20 @@ n_threads : int, optional
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 3}); auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
auto par_err = auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 3}); 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 = make_view_3d(y);
auto y_view_err = make_view_3d(y_err); auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x); auto x_view = make_view_1d(x);
aare::fit_gaus(x_view, y_view, y_view_err, par->view(), aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view(), n_threads); par_err->view(), chi2->view(), n_threads);
// return return_image_data(par);
return py::make_tuple(return_image_data(par), return py::dict("par"_a = return_image_data(par),
return_image_data(par_err)); "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) { } else if (y.ndim() == 1) {
// Allocate memory for the output // Allocate memory for the output
// Need to have pointers to allow python to manage // Need to have pointers to allow python to manage
@ -120,15 +128,20 @@ n_threads : int, optional
auto y_view_err = make_view_1d(y_err); auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x); auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_gaus(x_view, y_view, y_view_err, par->view(), aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view()); par_err->view(), chi2);
return py::make_tuple(return_image_data(par),
return_image_data(par_err)); 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 { } else {
throw std::runtime_error("Data must be 1D or 3D"); throw std::runtime_error("Data must be 1D or 3D");
} }
}, },
R"( R"(
Fit a 1D Gaussian to data with error estimates. Fit a 1D Gaussian to data with error estimates.
Parameters Parameters
@ -172,11 +185,10 @@ n_threads : int, optional
"fit_pol1", "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> 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,
py::array_t<double, py::array::c_style | py::array::forcecast> py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
y_err, int n_threads) { int n_threads) {
if (y.ndim() == 3) { if (y.ndim() == 3) {
auto par = auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
auto par_err = auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 2}); new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
@ -184,10 +196,15 @@ n_threads : int, optional
auto y_view_err = make_view_3d(y_err); auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x); auto x_view = make_view_1d(x);
aare::fit_pol1(x_view, y_view,y_view_err, par->view(), auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
par_err->view(), n_threads);
return py::make_tuple(return_image_data(par), aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
return_image_data(par_err)); 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) { } else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2}); auto par = new NDArray<double, 1>({2});
@ -197,15 +214,18 @@ n_threads : int, optional
auto y_view_err = make_view_1d(y_err); auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x); auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_pol1(x_view, y_view, y_view_err, par->view(), aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
par_err->view()); par_err->view(), chi2);
return py::make_tuple(return_image_data(par), return py::dict("par"_a = return_image_data(par),
return_image_data(par_err)); "par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
} else { } else {
throw std::runtime_error("Data must be 1D or 3D"); throw std::runtime_error("Data must be 1D or 3D");
} }
}, },
R"( R"(
Fit a 1D polynomial to data with error estimates. Fit a 1D polynomial to data with error estimates.
Parameters Parameters

284
src/Fit.cpp
View File

@ -1,10 +1,12 @@
#include "aare/Fit.hpp" #include "aare/Fit.hpp"
#include "aare/utils/task.hpp" #include "aare/utils/task.hpp"
#include "aare/utils/par.hpp"
#include <lmcurve2.h> #include <lmcurve2.h>
#include <lmfit.hpp> #include <lmfit.hpp>
#include <thread> #include <thread>
#include <array>
namespace aare { namespace aare {
@ -35,33 +37,11 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
} // namespace func } // namespace func
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) { NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
NDArray<double, 1> result({3}, 0); NDArray<double, 1> result = gaus_init_par(x, y);
lm_control_struct control = lm_control_double; lm_status_struct status;
// Estimate the initial parameters for the fit lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
std::vector<double> start_par{0, 0, 0}; aare::func::gaus, &lm_control_double, &status);
auto e = std::max_element(y.begin(), y.end());
auto idx = std::distance(y.begin(), e);
start_par[0] = *e; // For amplitude we use the maximum value
start_par[1] =
x[idx]; // For the mean we use the x value of the maximum value
// For sigma we estimate the fwhm and divide by 2.35
// assuming equally spaced x values
auto delta = x[1] - x[0];
start_par[2] =
std::count_if(y.begin(), y.end(),
[e, delta](double val) { return val > *e / 2; }) *
delta / 2.35;
lmfit::result_t res(start_par);
lmcurve(res.par.size(), res.par.data(), x.size(), x.data(), y.data(),
aare::func::gaus, &control, &res.status);
result(0) = res.par[0];
result(1) = res.par[1];
result(2) = res.par[2];
return result; return result;
} }
@ -81,65 +61,17 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
} }
} }
}; };
auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(process, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
return result; return result;
} }
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err, std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y) {
NDView<double, 3> par_out, NDView<double, 3> par_err_out, std::array<double, 3> start_par{0, 0, 0};
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);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(process, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
}
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) {
// 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");
}
lm_control_struct control = lm_control_double;
// Estimate the initial parameters for the fit
std::vector<double> start_par{0, 0, 0};
std::vector<double> start_par_err{0, 0, 0};
std::vector<double> start_cov{0, 0, 0, 0, 0, 0, 0, 0, 0};
auto e = std::max_element(y.begin(), y.end()); auto e = std::max_element(y.begin(), y.end());
auto idx = std::distance(y.begin(), e); auto idx = std::distance(y.begin(), e);
start_par[0] = *e; // For amplitude we use the maximum value start_par[0] = *e; // For amplitude we use the maximum value
start_par[1] = start_par[1] =
x[idx]; // For the mean we use the x value of the maximum value x[idx]; // For the mean we use the x value of the maximum value
@ -152,66 +84,82 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
[e, delta](double val) { return val > *e / 2; }) * [e, delta](double val) { return val > *e / 2; }) *
delta / 2.35; delta / 2.35;
lmfit::result_t res(start_par); return start_par;
lmfit::result_t res_err(start_par_err);
lmfit::result_t cov(start_cov);
// TODO can we make lmcurve write the result directly where is should be?
lmcurve2(res.par.size(), res.par.data(), res_err.par.data(), cov.par.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
&control, &res.status);
par_out(0) = res.par[0];
par_out(1) = res.par[1];
par_out(2) = res.par[2];
par_err_out(0) = res_err.par[0];
par_err_out(1) = res_err.par[1];
par_err_out(2) = res_err.par[2];
} }
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) { std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
// Estimate the initial parameters for the fit
std::array<double, 2> start_par{0, 0};
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
auto x1 = x[std::distance(y.begin(), y1)];
start_par[0] =
(*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value
start_par[1] =
*y1 - ((*y2 - *y1) / (x2 - x1)) *
x1; // For the mean we use the x value of the maximum value
return start_par;
}
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 // Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() || if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 2 || par_err_out.size() != 2) { par_out.size() != 3 || par_err_out.size() != 3) {
throw std::runtime_error("Data, x, data_err must have the same size " throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 2"); "and par_out, par_err_out must have size 3");
} }
lm_control_struct control = lm_control_double;
// Estimate the initial parameters for the fit // /* Collection of output parameters for status info. */
std::vector<double> start_par{0, 0}; // typedef struct {
std::vector<double> start_par_err{0, 0}; // double fnorm; /* norm of the residue vector fvec. */
std::vector<double> start_cov{0, 0, 0, 0}; // 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;
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
auto x1 = x[std::distance(y.begin(), y1)];
start_par[0] = lm_status_struct status;
(*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value par_out = gaus_init_par(x, y);
start_par[1] = std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0 , 0 , 0};
*y1 - ((*y2 - *y1) / (x2 - x1)) *
x1; // For the mean we use the x value of the maximum value
lmfit::result_t res(start_par); // 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);
lmfit::result_t res_err(start_par_err); // n_par - Number of free variables. Length of parameter vector par.
lmfit::result_t cov(start_cov); // 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(res.par.size(), res.par.data(), res_err.par.data(), cov.par.data(), 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, x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
&control, &res.status); &lm_control_double, &status);
par_out(0) = res.par[0]; // Calculate chi2
par_out(1) = res.par[1]; chi2 = 0;
par_err_out(0) = res_err.par[0]; for (size_t i = 0; i < y.size(); i++) {
par_err_out(1) = res_err.par[1]; chi2 += std::pow((y(i) - func::gaus(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, 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, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) { int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) { auto process = [&](ssize_t first_row, ssize_t last_row) {
@ -224,21 +172,64 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
{par_out.shape(2)}); {par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0), NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)}); {par_err_out.shape(2)});
fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view);
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); auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads; RunInParallel(process, tasks);
for (auto &task : tasks) { }
threads.push_back(std::thread(process, task.first, task.second));
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");
} }
for (auto &thread : threads) {
thread.join(); lm_status_struct status;
par_out = pol1_init_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 (size_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) { NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
// // Check that we have the correct sizes // // Check that we have the correct sizes
// if (y.size() != x.size() || y.size() != y_err.size() || // if (y.size() != x.size() || y.size() != y_err.size() ||
@ -246,28 +237,11 @@ NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
// throw std::runtime_error("Data, x, data_err must have the same size " // throw std::runtime_error("Data, x, data_err must have the same size "
// "and par_out, par_err_out must have size 2"); // "and par_out, par_err_out must have size 2");
// } // }
NDArray<double, 1> par({2}, 0); NDArray<double, 1> par = pol1_init_par(x, y);
lm_control_struct control = lm_control_double; lm_status_struct status;
lmcurve(par.size(), par.data(), x.size(), x.data(), y.data(),
// Estimate the initial parameters for the fit aare::func::pol1, &lm_control_double, &status);
std::vector<double> start_par{0, 0};
auto y2 = std::max_element(y.begin(), y.end());
auto x2 = x[std::distance(y.begin(), y2)];
auto y1 = std::min_element(y.begin(), y.end());
auto x1 = x[std::distance(y.begin(), y1)];
start_par[0] = (*y2 - *y1) / (x2 - x1);
start_par[1] = *y1 - ((*y2 - *y1) / (x2 - x1)) * x1;
lmfit::result_t res(start_par);
lmcurve(res.par.size(), res.par.data(), x.size(), x.data(), y.data(),
aare::func::pol1, &control, &res.status);
par(0) = res.par[0];
par(1) = res.par[1];
return par; return par;
} }
@ -287,13 +261,7 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
}; };
auto tasks = split_task(0, y.shape(0), n_threads); auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads; RunInParallel(process, tasks);
for (auto &task : tasks) {
threads.push_back(std::thread(process, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
return result; return result;
} }

28
src/NDArray.test.cpp
View File

@ -380,3 +380,31 @@ TEST_CASE("Elementwise operations on images") {
} }
} }
} }
TEST_CASE("Assign an std::array to a 1D NDArray") {
NDArray<int, 1> a{{5}, 0};
std::array<int, 5> b{1, 2, 3, 4, 5};
a = b;
for (uint32_t i = 0; i < a.size(); ++i) {
REQUIRE(a(i) == b[i]);
}
}
TEST_CASE("Assign an std::array to a 1D NDArray of a different size") {
NDArray<int, 1> a{{3}, 0};
std::array<int, 5> b{1, 2, 3, 4, 5};
a = b;
REQUIRE(a.size() == 5);
for (uint32_t i = 0; i < a.size(); ++i) {
REQUIRE(a(i) == b[i]);
}
}
TEST_CASE("Construct an NDArray from an std::array") {
std::array<int, 5> b{1, 2, 3, 4, 5};
NDArray<int, 1> a(b);
for (uint32_t i = 0; i < a.size(); ++i) {
REQUIRE(a(i) == b[i]);
}
}