Merge branch 'main' into developer

CMakeLists.txt

@@ -97,7 +97,6 @@ if(AARE_FETCH_LMFIT)
     
     FetchContent_MakeAvailable(lmfit)
     set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
-
 else()
     find_package(lmfit REQUIRED)
 endif()
@@ -369,6 +368,7 @@ target_link_libraries(
     PRIVATE 
     aare_compiler_flags 
     "$<BUILD_INTERFACE:lmfit>"
+
 )
 
 set_target_properties(aare_core PROPERTIES

conda-recipe/meta.yaml

@@ -3,6 +3,7 @@ package:
   version: 2025.2.18 #TODO! how to not duplicate this?
 
 
+
 source:
   path: ..
 

include/aare/Fit.hpp

@@ -17,6 +17,7 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
 
 } // namespace func
 
+
 /**
  * @brief Estimate the initial parameters for a Gaussian fit
  */
@@ -33,17 +34,20 @@ static constexpr int DEFAULT_NUM_THREADS = 4;
  */
 NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
 
+
 /**
  * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
  * @param x x values
  * @param y y 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);
 
 
 
+
 /**
  * @brief Fit a 1D Gaussian with error estimates
  * @param x x values
@@ -84,4 +88,5 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
               NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out,
               int n_threads = DEFAULT_NUM_THREADS);
 
+
 } // namespace aare

pyproject.toml

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

python/CMakeLists.txt

@@ -49,6 +49,7 @@ 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}  )

python/examples/play.py

@@ -8,6 +8,7 @@ import numpy as np
 import boost_histogram as bh
 import time
 
+
 import aare
 
 data = np.random.normal(10, 1, 1000)
@@ -19,4 +20,31 @@ hist.fill(data)
 x = hist.axes[0].centers
 y = hist.values()
 y_err = np.sqrt(y)+1
-res = aare.fit_gaus(x, y, y_err, chi2 = True)
+res = aare.fit_gaus(x, y, y_err, chi2 = True)
+
+
+        
+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)
+

python/src/fit.hpp

@@ -9,6 +9,7 @@
 namespace py = pybind11;
 using namespace pybind11::literals;
 
+
 void define_fit_bindings(py::module &m) {
 
     // TODO! Evaluate without converting to double
@@ -54,6 +55,7 @@ void define_fit_bindings(py::module &m) {
         )",
         py::arg("x"), py::arg("par"));
 
+
     m.def(
         "fit_gaus",
         [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
@@ -76,6 +78,7 @@ void define_fit_bindings(py::module &m) {
             }
         },
         R"(
+
 Fit a 1D Gaussian to data.
 
 Parameters
@@ -95,6 +98,7 @@ n_threads : int, optional
            py::array_t<double, py::array::c_style | py::array::forcecast> y,
            py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
            int n_threads) {
+
             if (y.ndim() == 3) {
                 // Allocate memory for the output
                 // Need to have pointers to allow python to manage
@@ -142,6 +146,7 @@ n_threads : int, optional
             }
         },
         R"(
+
 Fit a 1D Gaussian to data with error estimates.
 
 Parameters
@@ -189,6 +194,7 @@ n_threads : int, optional
            int n_threads) {
             if (y.ndim() == 3) {
                 auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
+
                 auto par_err =
                     new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
 
@@ -221,6 +227,7 @@ n_threads : int, optional
                 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");
             }

src/Fit.cpp

@@ -1,13 +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 {
@@ -160,6 +160,7 @@ 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, 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) {
@@ -175,6 +176,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
  
                 fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view,
                          chi2_out(row, col));
+
             }
         }
     };
@@ -185,6 +187,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
 
 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) {
@@ -221,13 +224,16 @@ 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, 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) {
@@ -242,6 +248,7 @@ NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
     lm_status_struct status;
     lmcurve(par.size(), par.data(), x.size(), x.data(), y.data(),
             aare::func::pol1, &lm_control_double, &status);
+
     return par;
 }
 
@@ -261,6 +268,7 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
     };
 
     auto tasks = split_task(0, y.shape(0), n_threads);
+
     RunInParallel(process, tasks);
     return result;
 }