Multi threaded fitting and returning chi2 (#132)

Co-authored-by: Patrick <patrick.sieberer@psi.ch>
Co-authored-by: JulianHeymes <julian.heymes@psi.ch>
Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch>

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 "")
@@ -97,9 +97,6 @@ if(AARE_FETCH_LMFIT)
     
     FetchContent_MakeAvailable(lmfit)
     set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
-
-    target_include_directories (lmfit PUBLIC "${libzmq_SOURCE_DIR}/lib")
-    message(STATUS "lmfit include dir: ${lmfit_SOURCE_DIR}/lib")
 else()
     find_package(lmfit REQUIRED)
 endif()
@@ -370,7 +367,8 @@ target_link_libraries(
     ${STD_FS_LIB} # from helpers.cmake
     PRIVATE 
     aare_compiler_flags 
-    lmfit
+    "$<BUILD_INTERFACE:lmfit>"
+
 )
 
 set_target_properties(aare_core PROPERTIES

conda-recipe/meta.yaml

@@ -1,6 +1,7 @@
 package:
   name: aare
-  version: 2025.2.12 #TODO! how to not duplicate this?
+  version: 2025.2.18 #TODO! how to not duplicate this?
+
 
 
 source:

include/aare/Fit.hpp

@@ -17,6 +17,14 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
 
 } // 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;
 
 /**
@@ -33,7 +41,11 @@ NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
  * @param y y vales, 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);
+
+NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
+                            int n_threads = DEFAULT_NUM_THREADS);
+
+
 
 
 /**
@@ -45,10 +57,12 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y, int n_thre
  * @param par_err_out output error parameters
  */
 void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
-              NDView<double, 1> par_out, NDView<double, 1> par_err_out);
+              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 y y vales, layout [row, col, values]
  * @param y_err error in y, layout [row, col, values]
@@ -57,20 +71,22 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
  * @param n_threads number of threads to use
  */
 void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
-              NDView<double, 3> par_out, NDView<double, 3> par_err_out, 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, 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++
-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

include/aare/NDArray.hpp

@@ -69,6 +69,11 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         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
     NDArray(NDArray &&other) noexcept
         : 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);
 
+    //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);
 
     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*/);
 
+  
+
+
+    
+
     void sqrt() {
         for (int i = 0; i < size_; ++i) {
             data_[i] = std::sqrt(data_[i]);
@@ -318,6 +342,9 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
     return *this;
 }
 
+
+
+
 template <typename T, int64_t Ndim>
 NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
     NDArray result = *this;
@@ -418,4 +445,6 @@ NDArray<T, Ndim> load(const std::string &pathname,
     return img;
 }
 
+
+
 } // namespace aare

include/aare/NDView.hpp

@@ -1,5 +1,5 @@
 #pragma once
-
+#include "aare/defs.hpp"
 #include "aare/ArrayExpr.hpp"
 
 #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>()); }
 
+
+    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) {
         for (auto it = begin(); it != end(); ++it)
             *it = val;

include/aare/utils/par.hpp

@@ -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

pyproject.toml

@@ -4,7 +4,8 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "aare"
-version = "2025.2.12"
+version = "2025.2.18"
+
 
 [tool.scikit-build]
 cmake.verbose = true

python/CMakeLists.txt

@@ -50,10 +50,10 @@ set(PYTHON_EXAMPLES
 )
 
 
-
 # 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})
 
 

python/examples/play.py

@@ -8,38 +8,20 @@ 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
 
-base = Path('/mnt/sls_det_storage/moench_data/Julian/MOENCH05/20250113_first_xrays_redo/raw_files/')
+import aare
-cluster_file = Path('/home/l_msdetect/erik/tmp/Cu.clust')
 
-t0 = time.perf_counter()
+data = np.random.normal(10, 1, 1000)
-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
-)
 
-total_clusters = 0
+hist = bh.Histogram(bh.axis.Regular(10, 0, 20))
-with ClusterFile(cluster_file, chunk_size = 1000) as f:
+hist.fill(data)
-    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)}')
+x = hist.axes[0].centers
+y = hist.values()
+y_err = np.sqrt(y)+1
+res = aare.fit_gaus(x, y, y_err, chi2 = True)
 
-for i in range(10):
-    header, img = f.read_frame()
-    print(header['frameNumber'], img.shape)
->>>>>>> developer
 
         
 t_elapsed = time.perf_counter()-t0

python/src/fit.hpp

@@ -7,6 +7,8 @@
 #include "aare/Fit.hpp"
 
 namespace py = pybind11;
+using namespace pybind11::literals;
+
 
 void define_fit_bindings(py::module &m) {
 
@@ -29,7 +31,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 +52,9 @@ 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",
@@ -73,6 +78,7 @@ void define_fit_bindings(py::module &m) {
             }
         },
         R"(
+
 Fit a 1D Gaussian to data.
 
 Parameters
@@ -90,8 +96,9 @@ n_threads : int, optional
         "fit_gaus",
         [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
            py::array_t<double, py::array::c_style | py::array::forcecast> y,
-           py::array_t<double, py::array::c_style | py::array::forcecast>
+           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) {
                 // Allocate memory for the output
                 // Need to have pointers to allow python to manage
@@ -99,15 +106,20 @@ n_threads : int, optional
                 auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
                 auto par_err =
                     new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
 
+                // Make views of the numpy arrays
                 auto y_view = make_view_3d(y);
                 auto y_view_err = make_view_3d(y_err);
                 auto x_view = make_view_1d(x);
+
                 aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
-                               par_err->view(), 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) {
                 // Allocate memory for the output
                 // Need to have pointers to allow python to manage
@@ -120,15 +132,21 @@ 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());
+                               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 {
                 throw std::runtime_error("Data must be 1D or 3D");
             }
         },
         R"(
+
 Fit a 1D Gaussian to data with error estimates.
 
 Parameters
@@ -172,11 +190,11 @@ n_threads : int, optional
         "fit_pol1",
         [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
            py::array_t<double, py::array::c_style | py::array::forcecast> y,
-           py::array_t<double, py::array::c_style | py::array::forcecast>
+           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) {
-                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 =
                     new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
 
@@ -184,10 +202,15 @@ n_threads : int, optional
                 auto y_view_err = make_view_3d(y_err);
                 auto x_view = make_view_1d(x);
 
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
+
                 aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
-                               par_err->view(), n_threads);
+                               par_err->view(), chi2->view(), n_threads);
-                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) - 2);
+
 
             } else if (y.ndim() == 1) {
                 auto par = new NDArray<double, 1>({2});
@@ -197,10 +220,14 @@ 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());
+                               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 {
                 throw std::runtime_error("Data must be 1D or 3D");
             }

src/Fit.cpp

@@ -1,11 +1,13 @@
 #include "aare/Fit.hpp"
 #include "aare/utils/task.hpp"
-
+#include "aare/utils/par.hpp"
 #include <lmcurve2.h>
 #include <lmfit.hpp>
-
 #include <thread>
 
+#include <array>
+
+
 namespace aare {
 
 namespace func {
@@ -35,33 +37,11 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
 } // namespace func
 
 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;
 }
@@ -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;
 }
 
-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 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,38 +84,14 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
                       [e, delta](double val) { return val > *e / 2; }) *
         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) {
-    // Check that we have the correct sizes
-    if (y.size() != x.size() || y.size() != y_err.size() ||
-        par_out.size() != 2 || par_err_out.size() != 2) {
-        throw std::runtime_error("Data, x, data_err must have the same size "
-                                 "and par_out, par_err_out must have size 2");
-    }
-
-    lm_control_struct control = lm_control_double;
 
+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::vector<double> start_par{0, 0};
+        std::array<double, 2> start_par{0, 0};
-    std::vector<double> start_par_err{0, 0};
+
-    std::vector<double> start_cov{0, 0, 0, 0};
     
         auto y2 = std::max_element(y.begin(), y.end());
         auto x2 = x[std::distance(y.begin(), y2)];
@@ -195,23 +103,64 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
         start_par[1] =
             *y1 - ((*y2 - *y1) / (x2 - x1)) *
                       x1; // For the mean we use the x value of the maximum value
-
+        return start_par;
-    lmfit::result_t res(start_par);
-    lmfit::result_t res_err(start_par_err);
-    lmfit::result_t cov(start_cov);
-
-    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);
-
-    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];
 }
 
-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, 1> y, NDView<double, 1> y_err,
-              NDView<double, 3> par_out, NDView<double, 3> par_err_out,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out,
+              double &chi2) {
+
+    // Check that we have the correct sizes
+    if (y.size() != x.size() || y.size() != y_err.size() ||
+        par_out.size() != 3 || par_err_out.size() != 3) {
+        throw std::runtime_error("Data, x, data_err must have the same size "
+                                 "and par_out, par_err_out must have size 3");
+    }
+
+
+    // /* Collection of output parameters for status info. */
+    // typedef struct {
+    //     double fnorm;  /* norm of the residue vector fvec. */
+    //     int nfev;      /* actual number of iterations. */
+    //     int outcome;   /* Status indicator. Nonnegative values are used as
+    //     index
+    //                     for the message text lm_infmsg, set in lmmin.c. */
+    //     int userbreak; /* Set when function evaluation requests termination.
+    //     */
+    // } lm_status_struct;
+
+
+    lm_status_struct status;
+    par_out = gaus_init_par(x, y);
+    std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0 , 0 , 0};
+
+    // void lmcurve2( const int n_par, double *par, double *parerr, double *covar, const int m_dat, const double *t, const double *y, const double *dy, double (*f)( const double ti, const double *par ), const lm_control_struct *control, lm_status_struct *status);
+    // n_par - Number of free variables. Length of parameter vector par.
+    // par - Parameter vector. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||r||.
+    // parerr - Parameter uncertainties vector. Array of length n_par or NULL. On output, unless it or covar is NULL, it contains the weighted parameter uncertainties for the found parameters.
+    // covar - Covariance matrix. Array of length n_par * n_par or NULL. On output, unless it is NULL, it contains the covariance matrix.
+    // m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy n_par <= m_dat.
+    // t - Array of length m_dat. Contains the abcissae (time, or "x") for which function f will be evaluated.
+    // y - Array of length m_dat. Contains the ordinate values that shall be fitted.
+    // dy - Array of length m_dat. Contains the standard deviations of the values y.
+    // f - A user-supplied parametric function f(ti;par).
+    // control - Parameter collection for tuning the fit procedure. In most cases, the default &lm_control_double is adequate. If f is only computed with single-precision accuracy, &lm_control_float should be used. Parameters are explained in lmmin2(3).
+    // status - A record used to return information about the minimization process: For details, see lmmin2(3).
+
+    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
+             x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
+             &lm_control_double, &status);
+
+    // Calculate chi2
+    chi2 = 0;
+    for (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_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out,  NDView<double, 2> chi2_out,
+
               int n_threads) {
 
     auto process = [&](ssize_t first_row, ssize_t last_row) {
@@ -224,19 +173,67 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
                                                {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);
+ 
+                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<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();
+void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
+
+    // Check that we have the correct sizes
+    if (y.size() != x.size() || y.size() != y_err.size() ||
+        par_out.size() != 2 || par_err_out.size() != 2) {
+        throw std::runtime_error("Data, x, data_err must have the same size "
+                                 "and par_out, par_err_out must have size 2");
     }
+
+    lm_status_struct status;
+    par_out = 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) {
@@ -246,28 +243,12 @@ 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 "
     //                      "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(),
+            aare::func::pol1, &lm_control_double, &status);
 
-    // Estimate the initial parameters for the fit
-    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;
 }
 
@@ -287,13 +268,8 @@ NDArray<double, 3> fit_pol1(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) {
+    RunInParallel(process, tasks);
-        threads.push_back(std::thread(process, task.first, task.second));
-    }
-    for (auto &thread : threads) {
-        thread.join();
-    }
     return result;
 }
 

src/NDArray.test.cpp

@@ -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]);
+    }
+}