diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index 06825572..960e66f5 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -10,13 +10,57 @@ namespace aare { namespace func { double gauss(const double x, const double *par); +NDArray gauss(NDView x, NDView par); + +double affine(const double x, const double *par); +NDArray affine(NDView x, NDView par); + } // namespace func +/** + * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values] + * @param data data to fit + * @param x x values + */ NDArray fit_gaus(NDView data, NDView x); + +/** + * @brief Fit a 1D Gaussian to data. + * @param data data to fit + * @param x x values + */ NDArray fit_gaus(NDView data, NDView x); -NDArray fit_gaus2(NDView data, NDView x, NDView data_err); -NDArray fit_gaus2(NDView data, NDView x, NDView data_err); + +/** + * @brief Fit a 1D gaus function to each pixel. Data layout [row, col, values] + * @param data data to fit + * @param x x values + * @param data_err error in data + */ +NDArray fit_gaus(NDView data, NDView x, NDView data_err); + +/** + * @brief Fit a 1D Gaussian to data. + * @param data data to fit + * @param x x values + * @param data_err error in data + */ +NDArray fit_gaus(NDView data, NDView x, NDView data_err); + +/** + * @brief Fit an affine function to each pixel. Data layout [row, col, values] + * @param data data to fit + * @param x x values + * @param data_err error in data + */ NDArray fit_affine(NDView data, NDView x, NDView data_err); + +/** + * @brief Fit an affine function to data. + * @param data data to fit + * @param x x values + * @param data_err error in data + */ NDArray fit_affine(NDView data, NDView x, NDView data_err); } // namespace aare \ No newline at end of file diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 89ad5e70..2aaa222a 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -28,6 +28,7 @@ target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags) set( PYTHON_FILES aare/__init__.py aare/CtbRawFile.py + aare/func.py aare/RawFile.py aare/transform.py aare/ScanParameters.py @@ -43,10 +44,17 @@ set_target_properties(_aare PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/aare ) +set(PYTHON_EXAMPLES + examples/play.py + examples/fits.py +) + -# Copy the examples/scripts to the build directory -configure_file(examples/play.py ${CMAKE_BINARY_DIR}/play.py) +# Copy the python examples to the build directory +foreach(FILE ${PYTHON_EXAMPLES}) + configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} ) +endforeach(FILE ${PYTHON_EXAMPLES}) if(AARE_INSTALL_PYTHONEXT) diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 2d03f07c..56a244b9 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -10,10 +10,14 @@ from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i -from ._aare import fit_gaus, fit_gaus2, fit_affine +from ._aare import fit_gaus, fit_affine from .CtbRawFile import CtbRawFile from .RawFile import RawFile from .ScanParameters import ScanParameters -from .utils import random_pixels, random_pixel \ No newline at end of file +from .utils import random_pixels, random_pixel, flat_list + + +#make functions available in the top level API +from .func import * diff --git a/python/aare/func.py b/python/aare/func.py new file mode 100644 index 00000000..d4a9c54a --- /dev/null +++ b/python/aare/func.py @@ -0,0 +1 @@ +from ._aare import gaus, affine \ No newline at end of file diff --git a/python/aare/utils.py b/python/aare/utils.py index d53f844c..de90e51c 100644 --- a/python/aare/utils.py +++ b/python/aare/utils.py @@ -20,4 +20,7 @@ def random_pixel(xmin=0, xmax=512, ymin=0, ymax=1024): Returns: tuple: (row, col) """ - return random_pixels(1, xmin, xmax, ymin, ymax)[0] \ No newline at end of file + return random_pixels(1, xmin, xmax, ymin, ymax)[0] + +def flat_list(xss): + return [x for xs in xss for x in xs] \ No newline at end of file diff --git a/python/examples/play.py b/python/examples/play.py index 121f835c..50d9add0 100644 --- a/python/examples/play.py +++ b/python/examples/play.py @@ -9,137 +9,48 @@ import time 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_xrays_scan_energy/raw_files/') +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') -pedestal_file = base/'moench05_noise_Cu_10_us_master_0.json' -data_file = base/'moench05_xray_Cu_10_us_master_0.json' -from aare._aare import fit_gaus -from aare.transform import moench05 - -offset = -0.5 +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(100, 3000, 4000), #ADU + 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] -rows = np.zeros(160*150) -cols = np.zeros(160*150) - -for row in range(160): - for col in range(150): - rows[row*150+col] = row - cols[row*150+col] = col - -#TODO how can we speed up the filling of the 3d histogram? -with CtbRawFile(pedestal_file, transform=moench05) as f: - for i in range(1000): - print(f'{i}', end='\r') - frame_number, frame = f.read_frame() - hist3d.fill(rows, cols, frame.flat) - -data = hist3d.values() -x = hist3d.axes[2].centers -a = fit_gaus(data, x) - - - - -# from aare._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink - - -# cf = ClusterFinderMT((400,400), (3,3), n_threads = 3) -# # collector = ClusterCollector(cf) -# out_file = ClusterFileSink(cf, "test.clust") - -# for i in range(1000): -# img = f.read_frame() -# cf.push_pedestal_frame(img) -# print('Pedestal done') -# cf.sync() - -# for i in range(100): -# img = f.read_frame() -# cf.find_clusters(img) - - -# # time.sleep(1) -# cf.stop() -# time.sleep(1) -# print('Second run') -# cf.start() -# for i in range(100): -# img = f.read_frame() -# cf.find_clusters(img) - -# cf.stop() -# print('Third run') -# cf.start() -# for i in range(129): -# img = f.read_frame() -# cf.find_clusters(img) - -# cf.stop() -# out_file.stop() -# print('Done') - - -# cfile = ClusterFile("test.clust") -# i = 0 -# while True: -# try: -# cv = cfile.read_frame() -# i+=1 -# except RuntimeError: -# break -# print(f'Read {i} frames') - - - - -# # cf = ClusterFinder((400,400), (3,3)) -# # for i in range(1000): -# # cf.push_pedestal_frame(f.read_frame()) - -# # fig, ax = plt.subplots() -# # im = ax.imshow(cf.pedestal()) -# # cf.pedestal() -# # cf.noise() - - - -# # N = 500 -# # t0 = time.perf_counter() -# # hist1 = bh.Histogram(bh.axis.Regular(40, -2, 4000)) -# # f.seek(0) - -# # t0 = time.perf_counter() -# # data = f.read_n(N) -# # t_elapsed = time.perf_counter()-t0 - - -# # n_bytes = data.itemsize*data.size - -# # print(f'Reading {N} frames took {t_elapsed:.3f}s {N/t_elapsed:.0f} FPS, {n_bytes/1024**2:.4f} GB/s') - - -# # for frame in data: -# # a = cf.find_clusters(frame) + +t_elapsed = time.perf_counter()-t0 +print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s') -# # clusters = cf.steal_clusters() +histogram_data = hist3d.counts() +x = hist3d.axes[2].edges[:-1] -# # t_elapsed = time.perf_counter()-t0 -# # print(f'Clustering {N} frames took {t_elapsed:.2f}s {N/t_elapsed:.0f} FPS') +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 -# # t0 = time.perf_counter() -# # total_clusters = clusters.size +# par = fit_gaus2(y,x, y_err) +# ax.plot(xx, gaus(xx,par)) +# print(par) -# # hist1.fill(clusters.sum()) +res = fit_gaus(y,x) +res2 = fit_gaus(y,x, y_err) +print(res) +print(res2) -# # t_elapsed = time.perf_counter()-t0 -# # print(f'Filling histogram with the sum of {total_clusters} clusters took: {t_elapsed:.3f}s, {total_clusters/t_elapsed:.3g} clust/s') -# # print(f'Average number of clusters per frame {total_clusters/N:.3f}') \ No newline at end of file diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 630e09b5..4f03ad66 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -9,81 +9,100 @@ namespace py = pybind11; void define_fit_bindings(py::module &m) { + + //TODO! Evaluate without converting to double + m.def( + "gaus", + [](py::array_t x, + py::array_t par) { + auto x_view = make_view_1d(x); + auto par_view = make_view_1d(par); + auto y = new NDArray{}; + *y = aare::func::gauss(x_view, par_view); + return return_image_data(y); + }); + + m.def("affine", + [](py::array_t x, + py::array_t par) { + auto x_view = make_view_1d(x); + auto par_view = make_view_1d(par); + auto y = new NDArray{}; + *y = aare::func::affine(x_view, par_view); + return return_image_data(y); + }); + m.def( "fit_gaus", [](py::array_t data, py::array_t x) { - - if(data.ndim() == 3){ + if (data.ndim() == 3) { auto par = new NDArray{}; auto data_view = make_view_3d(data); auto x_view = make_view_1d(x); *par = aare::fit_gaus(data_view, x_view); return return_image_data(par); - }else if(data.ndim() == 1){ + } else if (data.ndim() == 1) { auto par = new NDArray{}; auto data_view = make_view_1d(data); auto x_view = make_view_1d(x); *par = aare::fit_gaus(data_view, x_view); return return_image_data(par); - }else{ + } else { throw std::runtime_error("Data must be 1D or 3D"); } - }, py::arg("data"), py::arg("x")); - m.def( - "fit_gaus2", + "fit_gaus", [](py::array_t data, - py::array_t x, - py::array_t data_err) { - - if(data.ndim() == 3){ + py::array_t x, + py::array_t + data_err) { + if (data.ndim() == 3) { auto par = new NDArray{}; auto data_view = make_view_3d(data); auto data_view_err = make_view_3d(data_err); auto x_view = make_view_1d(x); - *par = aare::fit_gaus2(data_view, x_view, data_view_err); + *par = aare::fit_gaus(data_view, x_view, data_view_err); return return_image_data(par); - }else if(data.ndim() == 1){ + } else if (data.ndim() == 1) { auto par = new NDArray{}; auto data_view = make_view_1d(data); auto data_view_err = make_view_1d(data_err); auto x_view = make_view_1d(x); - *par = aare::fit_gaus2(data_view, x_view, data_view_err); + *par = aare::fit_gaus(data_view, x_view, data_view_err); return return_image_data(par); - }else{ + } else { throw std::runtime_error("Data must be 1D or 3D"); } - }, py::arg("data"), py::arg("x"), py::arg("data_err")); - m.def( + + m.def( "fit_affine", [](py::array_t data, - py::array_t x, - py::array_t data_err) { - - if(data.ndim() == 3){ + py::array_t x, + py::array_t + data_err) { + if (data.ndim() == 3) { auto par = new NDArray{}; auto data_view = make_view_3d(data); auto data_view_err = make_view_3d(data_err); auto x_view = make_view_1d(x); *par = aare::fit_affine(data_view, x_view, data_view_err); return return_image_data(par); - }else if(data.ndim() == 1){ + } else if (data.ndim() == 1) { auto par = new NDArray{}; auto data_view = make_view_1d(data); auto data_view_err = make_view_1d(data_err); auto x_view = make_view_1d(x); *par = aare::fit_affine(data_view, x_view, data_view_err); return return_image_data(par); - }else{ + } else { throw std::runtime_error("Data must be 1D or 3D"); } - }, py::arg("data"), py::arg("x"), py::arg("data_err")); } \ No newline at end of file diff --git a/src/Fit.cpp b/src/Fit.cpp index f39ef711..f4e14005 100644 --- a/src/Fit.cpp +++ b/src/Fit.cpp @@ -10,72 +10,48 @@ double gauss(const double x, const double *par) { return par[0] * exp(-pow(x - par[1], 2) / (2 * pow(par[2], 2))); } +NDArray gauss(NDView x, NDView par) { + NDArray y({x.size()}, 0); + for (size_t i = 0; i < x.size(); i++) { + y(i) = gauss(x(i), par.data()); + } + return y; +} + double affine(const double x, const double *par) { return par[0] * x + par[1]; } +NDArray affine(NDView x, NDView par) { + NDArray y({x.size()}, 0); + for (size_t i = 0; i < x.size(); i++) { + y(i) = affine(x(i), par.data()); + } + return y; +} + } // namespace func NDArray fit_gaus(NDView data, NDView x) { - // std::vector start_par{300, 3300, 50}; NDArray result({data.shape(0), data.shape(1), 3}, 0); - lm_control_struct control = lm_control_double; - - // fmt::print("data.shape(0): {}\n", data.shape(0)); - // fmt::print("data.shape(1): {}\n", data.shape(1)); - // fmt::print("data.shape(2): {}\n", data.shape(2)); - - auto ptr = data.data(); - for (ssize_t row = 0; row < data.shape(0); row++) { for (ssize_t col = 0; col < data.shape(1); col++) { - - //Estimate the initial parameters for the fit - std::vector start_par{0, 0, 0}; - auto e = std::max_element(ptr, ptr + data.shape(2)); - auto idx = std::distance(ptr, 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(ptr, ptr + data.shape(2), - [e, delta](double val) { return val > *e / 2; }) * delta / 2.35; - - // if (row == 0 && col == 0) { - // fmt::print("start_par: {} {} {}\n", start_par[0], start_par[1], start_par[2]); - // for (size_t i = 0; i < data.shape(2); i++) { - // fmt::print("data({},{},{}): {} {}\n", row, col, i, x[i], data(row, col, i)); - // } - // } - - lmfit::result_t res(start_par); - lmcurve(res.par.size(), res.par.data(), x.size(), x.data(), - &data(row, col,0), aare::func::gauss, &control, &res.status); - // fmt::print("par: {} {} {}\n", res.par[0], res.par[1], res.par[2]); - - result(row, col, 0) = res.par[0]; - result(row, col, 1) = res.par[1]; - result(row, col, 2) = res.par[2]; - - ptr += data.shape(2); + NDView y(&data(row, col, 0), {data.shape(2)}); + auto res = fit_gaus(y, x); + result(row, col, 0) = res(0); + result(row, col, 1) = res(1); + result(row, col, 2) = res(2); } - } return result; } NDArray fit_gaus(NDView data, NDView x) { - // std::vector start_par{300, 3300, 50}; - NDArray result({3}, 0); lm_control_struct control = lm_control_double; - //Estimate the initial parameters for the fit std::vector start_par{0, 0, 0}; auto e = std::max_element(data.begin(), data.end()); @@ -90,7 +66,7 @@ NDArray fit_gaus(NDView data, NDView x) { start_par[2] = std::count_if(data.begin(), data.end(), [e, delta](double val) { return val > *e / 2; }) * delta / 2.35; - + // fmt::print("start_par: {} {} {}\n", start_par[0], start_par[1], start_par[2]); lmfit::result_t res(start_par); lmcurve(res.par.size(), res.par.data(), x.size(), x.data(), data.data(), aare::func::gauss, &control, &res.status); @@ -101,79 +77,34 @@ NDArray fit_gaus(NDView data, NDView x) { return result; } -// ============================================================================================ -// With errors -NDArray fit_gaus2(NDView data, NDView x, NDView data_err) { - // std::vector start_par{300, 3300, 50}; - NDArray result({data.shape(0), data.shape(1), 6}, 0); - // NDArray result_error({data.shape(0), data.shape(1), 3}, 0); - lm_control_struct control = lm_control_double; - - // fmt::print("data.shape(0): {}\n", data.shape(0)); - // fmt::print("data.shape(1): {}\n", data.shape(1)); - // fmt::print("data.shape(2): {}\n", data.shape(2)); - - auto ptr = data.data(); +NDArray fit_gaus(NDView data, NDView x, NDView data_err) { + NDArray result({data.shape(0), data.shape(1), 6}, 0); for (ssize_t row = 0; row < data.shape(0); row++) { for (ssize_t col = 0; col < data.shape(1); col++) { - - //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}; - - auto e = std::max_element(ptr, ptr + data.shape(2)); - auto idx = std::distance(ptr, 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(ptr, ptr + data.shape(2), - [e, delta](double val) { return val > *e / 2; }) * delta / 2.35; - - // if (row == 0 && col == 0) { - // fmt::print("start_par: {} {} {}\n", start_par[0], start_par[1], start_par[2]); - // for (size_t i = 0; i < data.shape(2); i++) { - // fmt::print("data({},{},{}): {} {}\n", row, col, i, x[i], data(row, col, i)); - // } - // } - - 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(), - &data(row, col,0), &data_err(row, col,0), aare::func::gauss, &control, &res.status); - // fmt::print("par: {} {} {}\n", res.par[0], res.par[1], res.par[2]); - - result(row, col, 0) = res.par[0]; - result(row, col, 1) = res.par[1]; - result(row, col, 2) = res.par[2]; - result(row, col, 3) = res_err.par[0]; - result(row, col, 4) = res_err.par[1]; - result(row, col, 5) = res_err.par[2]; - ptr += data.shape(2); + NDView y(&data(row, col, 0), {data.shape(2)}); + NDView y_err(&data_err(row, col, 0), {data.shape(2)}); + auto res = fit_gaus(y, x, y_err); + result(row, col, 0) = res[0]; + result(row, col, 1) = res[1]; + result(row, col, 2) = res[2]; + result(row, col, 3) = res[3]; + result(row, col, 4) = res[4]; + result(row, col, 5) = res[5]; } } return result; } -// With errors -NDArray fit_gaus2(NDView data, NDView x, NDView data_err) { - // std::vector start_par{300, 3300, 50}; + +NDArray fit_gaus(NDView data, NDView x, NDView data_err) { NDArray result({6}, 0); - // NDArray result_err({3}, 0); 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}; @@ -181,7 +112,6 @@ NDArray fit_gaus2(NDView data, NDView x, NDView auto e = std::max_element(data.begin(), data.end()); auto idx = std::distance(data.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 @@ -199,6 +129,7 @@ NDArray fit_gaus2(NDView data, NDView x, NDView lmcurve2(res.par.size(), res.par.data(), res_err.par.data(), cov.par.data(), x.size(), x.data(), data.data(), data_err.data(), aare::func::gauss, &control, &res.status); + result(0) = res.par[0]; result(1) = res.par[1]; result(2) = res.par[2]; @@ -208,56 +139,19 @@ NDArray fit_gaus2(NDView data, NDView x, NDView return result; } -// ---- AFFINE FUNCTION ---- -NDArray fit_affine(NDView data, NDView x, NDView data_err) { - // std::vector start_par{300, 3300, 50}; +NDArray fit_affine(NDView data, NDView x, NDView data_err) { NDArray result({data.shape(0), data.shape(1), 4}, 0); - // NDArray result_error({data.shape(0), data.shape(1), 3}, 0); - - lm_control_struct control = lm_control_double; - - // fmt::print("data.shape(0): {}\n", data.shape(0)); - // fmt::print("data.shape(1): {}\n", data.shape(1)); - // fmt::print("data.shape(2): {}\n", data.shape(2)); - - auto ptr = data.data(); - for (ssize_t row = 0; row < data.shape(0); row++) { for (ssize_t col = 0; col < data.shape(1); col++) { - - //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}; - - - - auto y2 = std::max_element(ptr, ptr + data.shape(2)); - auto x2 = x[std::distance(ptr, y2)]; - - auto y1 = std::min_element(ptr, ptr + data.shape(2)); - auto x1 = x[std::distance(ptr, 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 - - - 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(), - &data(row, col,0), &data_err(row, col,0), aare::func::affine, &control, &res.status); - // fmt::print("par: {} {} {}\n", res.par[0], res.par[1], res.par[2]); - - result(row, col, 0) = res.par[0]; - result(row, col, 1) = res.par[1]; - result(row, col, 2) = res_err.par[0]; - result(row, col, 3) = res_err.par[1]; - ptr += data.shape(2); + NDView y(&data(row, col, 0), {data.shape(2)}); + NDView y_err(&data_err(row, col, 0), {data.shape(2)}); + auto res = fit_affine(y, x, y_err); + result(row, col, 0) = res[0]; + result(row, col, 1) = res[1]; + result(row, col, 2) = res[2]; + result(row, col, 3) = res[3]; } - } return result; }