From 6a8398848519c307200ad01bb7c4e96b805196f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Tue, 18 Feb 2025 21:13:27 +0100 Subject: [PATCH] Added chi2 to fit results (#131) - fit_gaus and fit_pol1 now return a dict - calculate chi2 after fit - cleaned up code --- CMakeLists.txt | 2 +- include/aare/Fit.hpp | 39 +++-- include/aare/NDArray.hpp | 29 ++++ include/aare/NDView.hpp | 11 +- include/aare/utils/par.hpp | 18 +++ python/CMakeLists.txt | 3 +- python/examples/play.py | 62 ++------ python/src/fit.hpp | 70 +++++---- src/Fit.cpp | 284 ++++++++++++++++--------------------- src/NDArray.test.cpp | 28 ++++ 10 files changed, 291 insertions(+), 255 deletions(-) create mode 100644 include/aare/utils/par.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 624259a..a5a576a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -85,7 +85,7 @@ if(AARE_FETCH_LMFIT) GIT_TAG main PATCH_COMMAND ${lmfit_patch} UPDATE_DISCONNECTED 1 - EXCLUDE_FROM_ALL + EXCLUDE_FROM_ALL 1 ) #Disable what we don't need from lmfit set(BUILD_TESTING OFF CACHE BOOL "") diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index 20ef4ef..6400edd 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -17,6 +17,13 @@ NDArray pol1(NDView x, NDView par); } // namespace func +/** + * @brief Estimate the initial parameters for a Gaussian fit + */ +std::array gaus_init_par(const NDView x, const NDView y); + +std::array pol1_init_par(const NDView x, const NDView y); + static constexpr int DEFAULT_NUM_THREADS = 4; /** @@ -26,14 +33,15 @@ static constexpr int DEFAULT_NUM_THREADS = 4; */ NDArray fit_gaus(NDView x, NDView y); - /** * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values] * @param x x values * @param y y vales, layout [row, col, values] * @param n_threads number of threads to use */ -NDArray fit_gaus(NDView x, NDView y, int n_threads = DEFAULT_NUM_THREADS); +NDArray fit_gaus(NDView x, NDView y, + int n_threads = DEFAULT_NUM_THREADS); + /** @@ -45,10 +53,12 @@ NDArray fit_gaus(NDView x, NDView y, int n_thre * @param par_err_out output error parameters */ void fit_gaus(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out); + NDView par_out, NDView 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 y y vales, layout [row, col, values] * @param y_err error in y, layout [row, col, values] @@ -57,20 +67,21 @@ void fit_gaus(NDView x, NDView y, NDView y_err, * @param n_threads number of threads to use */ void fit_gaus(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, int n_threads = DEFAULT_NUM_THREADS); - + NDView par_out, NDView par_err_out, NDView chi2_out, + int n_threads = DEFAULT_NUM_THREADS + ); NDArray fit_pol1(NDView x, NDView y); -NDArray fit_pol1(NDView x, NDView y, int n_threads = DEFAULT_NUM_THREADS); +NDArray fit_pol1(NDView x, NDView y, + int n_threads = DEFAULT_NUM_THREADS); -void fit_pol1(NDView x, NDView y, - NDView y_err, NDView par_out, - NDView par_err_out); +void fit_pol1(NDView x, NDView y, NDView y_err, + NDView par_out, NDView par_err_out, double& chi2); -//TODO! not sure we need to offer the different version in C++ -void fit_pol1(NDView x, NDView y, - NDView y_err, NDView par_out, - NDView par_err_out, int n_threads = DEFAULT_NUM_THREADS); +// TODO! not sure we need to offer the different version in C++ +void fit_pol1(NDView x, NDView y, NDView y_err, + NDView par_out, NDView par_err_out,NDView chi2_out, + int n_threads = DEFAULT_NUM_THREADS); } // namespace aare \ No newline at end of file diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 15beb02..cfa5b5c 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -69,6 +69,11 @@ class NDArray : public ArrayExpr, Ndim> { std::copy(v.begin(), v.end(), begin()); } + template + NDArray(const std::array& arr) : NDArray({Size}) { + std::copy(arr.begin(), arr.end(), begin()); + } + // Move constructor NDArray(NDArray &&other) noexcept : shape_(other.shape_), strides_(c_strides(shape_)), @@ -105,6 +110,20 @@ class NDArray : public ArrayExpr, Ndim> { NDArray &operator-=(const NDArray &other); NDArray &operator*=(const NDArray &other); + //Write directly to the data array, or create a new one + template + NDArray& operator=(const std::array &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); template NDArray &operator/=(const NDArray &other) { @@ -135,6 +154,11 @@ class NDArray : public ArrayExpr, Ndim> { NDArray &operator&=(const T & /*mask*/); + + + + + void sqrt() { for (int i = 0; i < size_; ++i) { data_[i] = std::sqrt(data_[i]); @@ -318,6 +342,9 @@ NDArray &NDArray::operator+=(const T &value) { return *this; } + + + template NDArray NDArray::operator+(const T &value) { NDArray result = *this; @@ -418,4 +445,6 @@ NDArray load(const std::string &pathname, return img; } + + } // namespace aare \ No newline at end of file diff --git a/include/aare/NDView.hpp b/include/aare/NDView.hpp index e3a6d30..f53f758 100644 --- a/include/aare/NDView.hpp +++ b/include/aare/NDView.hpp @@ -1,5 +1,5 @@ #pragma once - +#include "aare/defs.hpp" #include "aare/ArrayExpr.hpp" #include @@ -99,6 +99,15 @@ template class NDView : public ArrayExpr()); } + + template + NDView& operator=(const std::array &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) { for (auto it = begin(); it != end(); ++it) *it = val; diff --git a/include/aare/utils/par.hpp b/include/aare/utils/par.hpp new file mode 100644 index 0000000..efb1c77 --- /dev/null +++ b/include/aare/utils/par.hpp @@ -0,0 +1,18 @@ +#include +#include +#include + +namespace aare { + + template + void RunInParallel(F func, const std::vector>& tasks) { + // auto tasks = split_task(0, y.shape(0), n_threads); + std::vector threads; + for (auto &task : tasks) { + threads.push_back(std::thread(func, task.first, task.second)); + } + for (auto &thread : threads) { + thread.join(); + } + } +} // namespace aare \ No newline at end of file diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 2aaa222..4cbd69d 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -49,11 +49,10 @@ set(PYTHON_EXAMPLES examples/fits.py ) - - # Copy the python examples to the build directory foreach(FILE ${PYTHON_EXAMPLES}) configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} ) + message(STATUS "Copying ${FILE} to ${CMAKE_BINARY_DIR}/${FILE}") endforeach(FILE ${PYTHON_EXAMPLES}) diff --git a/python/examples/play.py b/python/examples/play.py index f1a869b..82aaa7d 100644 --- a/python/examples/play.py +++ b/python/examples/play.py @@ -8,61 +8,15 @@ import numpy as np import boost_histogram as bh import time -<<<<<<< HEAD -from aare import File, ClusterFinder, VarClusterFinder, ClusterFile, CtbRawFile -from aare import gaus, fit_gaus +import aare -base = Path('/mnt/sls_det_storage/moench_data/Julian/MOENCH05/20250113_first_xrays_redo/raw_files/') -cluster_file = Path('/home/l_msdetect/erik/tmp/Cu.clust') +data = np.random.normal(10, 1, 1000) -t0 = time.perf_counter() -offset= -0.5 -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 -) +hist = bh.Histogram(bh.axis.Regular(10, 0, 20)) +hist.fill(data) -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 -print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s') - -histogram_data = hist3d.counts() -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) +x = hist.axes[0].centers +y = hist.values() +y_err = np.sqrt(y)+1 +res = aare.fit_gaus(x, y, y_err, chi2 = True) \ No newline at end of file diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 60cdecc..2506e9b 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -7,6 +7,7 @@ #include "aare/Fit.hpp" namespace py = pybind11; +using namespace pybind11::literals; 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. 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")); + )", + py::arg("x"), py::arg("par")); m.def( "pol1", @@ -49,7 +51,8 @@ void define_fit_bindings(py::module &m) { 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")); + )", + py::arg("x"), py::arg("par")); m.def( "fit_gaus", @@ -72,7 +75,7 @@ void define_fit_bindings(py::module &m) { throw std::runtime_error("Data must be 1D or 3D"); } }, -R"( + R"( Fit a 1D Gaussian to data. Parameters @@ -90,8 +93,8 @@ n_threads : int, optional "fit_gaus", [](py::array_t x, py::array_t y, - py::array_t - y_err, int n_threads) { + py::array_t y_err, + int n_threads) { if (y.ndim() == 3) { // Allocate memory for the output // Need to have pointers to allow python to manage @@ -99,15 +102,20 @@ n_threads : int, optional auto par = new NDArray({y.shape(0), y.shape(1), 3}); auto par_err = new NDArray({y.shape(0), y.shape(1), 3}); + auto chi2 = new NDArray({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(), n_threads); - // return return_image_data(par); - return py::make_tuple(return_image_data(par), - 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) - 3); } else if (y.ndim() == 1) { // Allocate memory for the output // 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 x_view = make_view_1d(x); + + double chi2 = 0; aare::fit_gaus(x_view, y_view, y_view_err, par->view(), - par_err->view()); - return py::make_tuple(return_image_data(par), - return_image_data(par_err)); + 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"( + R"( Fit a 1D Gaussian to data with error estimates. Parameters @@ -172,11 +185,10 @@ n_threads : int, optional "fit_pol1", [](py::array_t x, py::array_t y, - py::array_t - y_err, int n_threads) { + py::array_t y_err, + int n_threads) { if (y.ndim() == 3) { - auto par = - new NDArray({y.shape(0), y.shape(1), 2}); + auto par = new NDArray({y.shape(0), y.shape(1), 2}); auto par_err = new NDArray({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 x_view = make_view_1d(x); - aare::fit_pol1(x_view, y_view,y_view_err, par->view(), - par_err->view(), n_threads); - return py::make_tuple(return_image_data(par), - return_image_data(par_err)); + auto chi2 = new NDArray({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({2}); @@ -197,15 +214,18 @@ n_threads : int, optional 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()); - return py::make_tuple(return_image_data(par), - return_image_data(par_err)); + 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"( + R"( Fit a 1D polynomial to data with error estimates. Parameters diff --git a/src/Fit.cpp b/src/Fit.cpp index 08ecaec..c8ce178 100644 --- a/src/Fit.cpp +++ b/src/Fit.cpp @@ -1,10 +1,12 @@ #include "aare/Fit.hpp" #include "aare/utils/task.hpp" +#include "aare/utils/par.hpp" #include #include #include +#include namespace aare { @@ -35,33 +37,11 @@ NDArray pol1(NDView x, NDView par) { } // namespace func NDArray fit_gaus(NDView x, NDView y) { - NDArray result({3}, 0); - lm_control_struct control = lm_control_double; + NDArray result = gaus_init_par(x, y); + lm_status_struct status; - // Estimate the initial parameters for the fit - std::vector start_par{0, 0, 0}; - 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]; + lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(), + aare::func::gaus, &lm_control_double, &status); return result; } @@ -81,65 +61,17 @@ NDArray fit_gaus(NDView x, NDView y, } } }; - auto tasks = split_task(0, y.shape(0), n_threads); - std::vector 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; } -void fit_gaus(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_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 y_view(&y(row, col, 0), {y.shape(2)}); - NDView y_err_view(&y_err(row, col, 0), - {y_err.shape(2)}); - NDView par_out_view(&par_out(row, col, 0), - {par_out.shape(2)}); - NDView 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 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 x, NDView y, NDView y_err, - NDView par_out, NDView 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 start_par{0, 0, 0}; - std::vector start_par_err{0, 0, 0}; - std::vector start_cov{0, 0, 0, 0, 0, 0, 0, 0, 0}; - +std::array gaus_init_par(const NDView x, const NDView y) { + std::array start_par{0, 0, 0}; 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 @@ -152,66 +84,82 @@ void fit_gaus(NDView x, NDView y, NDView y_err, [e, delta](double val) { return val > *e / 2; }) * delta / 2.35; - lmfit::result_t res(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]; + return start_par; } -void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out) { + +std::array pol1_init_par(const NDView x, const NDView y){ + // Estimate the initial parameters for the fit + std::array 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 x, NDView y, NDView y_err, + NDView par_out, NDView 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) { + 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 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 - std::vector start_par{0, 0}; - std::vector start_par_err{0, 0}; - std::vector start_cov{0, 0, 0, 0}; + // /* 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; - 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 + lm_status_struct status; + par_out = gaus_init_par(x, y); + std::array cov{0, 0, 0, 0, 0, 0, 0 , 0 , 0}; - lmfit::result_t res(start_par); - lmfit::result_t res_err(start_par_err); - lmfit::result_t cov(start_cov); + // 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(res.par.size(), res.par.data(), res_err.par.data(), cov.par.data(), - x.size(), x.data(), y.data(), y_err.data(), aare::func::pol1, - &control, &res.status); + 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); - par_out(0) = res.par[0]; - par_out(1) = res.par[1]; - par_err_out(0) = res_err.par[0]; - par_err_out(1) = res_err.par[1]; + // Calculate chi2 + chi2 = 0; + for (size_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_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, +void fit_gaus(NDView x, NDView y, NDView y_err, + NDView par_out, NDView par_err_out, NDView chi2_out, int n_threads) { auto process = [&](ssize_t first_row, ssize_t last_row) { @@ -224,21 +172,64 @@ void fit_pol1(NDView x, NDView y, NDView y_err, {par_out.shape(2)}); NDView 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); + + 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); - std::vector threads; - for (auto &task : tasks) { - threads.push_back(std::thread(process, task.first, task.second)); + RunInParallel(process, tasks); +} + +void fit_pol1(NDView x, NDView y, NDView y_err, + NDView par_out, NDView 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 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 x, NDView y, NDView y_err, + NDView par_out, NDView par_err_out, NDView 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 y_view(&y(row, col, 0), {y.shape(2)}); + NDView y_err_view(&y_err(row, col, 0), + {y_err.shape(2)}); + NDView par_out_view(&par_out(row, col, 0), + {par_out.shape(2)}); + NDView 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 fit_pol1(NDView x, NDView y) { // // Check that we have the correct sizes // if (y.size() != x.size() || y.size() != y_err.size() || @@ -246,28 +237,11 @@ NDArray fit_pol1(NDView x, NDView y) { // throw std::runtime_error("Data, x, data_err must have the same size " // "and par_out, par_err_out must have size 2"); // } - NDArray par({2}, 0); + NDArray par = pol1_init_par(x, y); - lm_control_struct control = lm_control_double; - - // Estimate the initial parameters for the fit - std::vector 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]; + lm_status_struct status; + lmcurve(par.size(), par.data(), x.size(), x.data(), y.data(), + aare::func::pol1, &lm_control_double, &status); return par; } @@ -287,13 +261,7 @@ NDArray fit_pol1(NDView x, NDView y, }; auto tasks = split_task(0, y.shape(0), n_threads); - std::vector threads; - for (auto &task : tasks) { - threads.push_back(std::thread(process, task.first, task.second)); - } - for (auto &thread : threads) { - thread.join(); - } + RunInParallel(process, tasks); return result; } diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 54099fd..942481c 100644 --- a/src/NDArray.test.cpp +++ b/src/NDArray.test.cpp @@ -379,4 +379,32 @@ TEST_CASE("Elementwise operations on images") { REQUIRE(A(i) == a_val); } } +} + +TEST_CASE("Assign an std::array to a 1D NDArray") { + NDArray a{{5}, 0}; + std::array 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 a{{3}, 0}; + std::array 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 b{1, 2, 3, 4, 5}; + NDArray a(b); + for (uint32_t i = 0; i < a.size(); ++i) { + REQUIRE(a(i) == b[i]); + } } \ No newline at end of file