diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 45df8a0..6ec2f2d 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -282,7 +282,7 @@ template ClusterVector ClusterFile::read_clusters_with_cut(size_t n_clusters) { ClusterVector clusters; - clusters.resize(n_clusters); + clusters.reserve(n_clusters); // if there are photons left from previous frame read them first if (m_num_left) { diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 22315cc..e85a6f0 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -384,6 +384,14 @@ class ClusterVector> { void set_frame_number(uint64_t frame_number) { m_frame_number = frame_number; } + + std::vector sum() { + std::vector sums(m_data.size()); + for (size_t i = 0; i < m_data.size(); i++) { + sums[i] = m_data[i].sum(); + } + return sums; + } }; } // namespace aare \ No newline at end of file diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 09de736..9f54049 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -28,6 +28,9 @@ target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags) set( PYTHON_FILES aare/__init__.py aare/CtbRawFile.py + aare/ClusterFinder.py + aare/ClusterVector.py + aare/func.py aare/RawFile.py aare/transform.py @@ -35,6 +38,7 @@ set( PYTHON_FILES aare/utils.py ) + # Copy the python files to the build directory foreach(FILE ${PYTHON_FILES}) configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} ) diff --git a/python/aare/ClusterFinder.py b/python/aare/ClusterFinder.py new file mode 100644 index 0000000..a2042a4 --- /dev/null +++ b/python/aare/ClusterFinder.py @@ -0,0 +1,14 @@ + +from ._aare import ClusterFinder_Cluster3x3i +import numpy as np + +def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024): + """ + Factory function to create a ClusterFinder object. Provides a cleaner syntax for + the templated ClusterFinder in C++. + """ + if dtype == np.int32 and cluster_size == (3,3): + return ClusterFinder_Cluster3x3i(image_size, n_sigma = n_sigma, capacity=capacity) + else: + #TODO! add the other formats + raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.") diff --git a/python/aare/ClusterVector.py b/python/aare/ClusterVector.py new file mode 100644 index 0000000..b0dd453 --- /dev/null +++ b/python/aare/ClusterVector.py @@ -0,0 +1,11 @@ + + +from ._aare import ClusterVector_Cluster3x3i +import numpy as np + +def ClusterVector(cluster_size, dtype = np.int32): + + if dtype == np.int32 and cluster_size == (3,3): + return ClusterVector_Cluster3x3i() + else: + raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.") diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 8c51d73..b1eb604 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -11,8 +11,12 @@ from ._aare import ROI # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i +from .ClusterFinder import ClusterFinder +from .ClusterVector import ClusterVector + from ._aare import fit_gaus, fit_pol1 from ._aare import Interpolator +from ._aare import calculate_eta2 from .CtbRawFile import CtbRawFile from .RawFile import RawFile from .ScanParameters import ScanParameters diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp new file mode 100644 index 0000000..f7fa796 --- /dev/null +++ b/python/src/bind_ClusterVector.hpp @@ -0,0 +1,103 @@ +#include "aare/ClusterCollector.hpp" +#include "aare/ClusterFileSink.hpp" +#include "aare/ClusterFinder.hpp" +#include "aare/ClusterFinderMT.hpp" +#include "aare/ClusterVector.hpp" +#include "aare/NDView.hpp" +#include "aare/Pedestal.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include +#include + +namespace py = pybind11; +using pd_type = double; + +using namespace aare; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" + + +template +void define_ClusterVector(py::module &m, const std::string &typestr) { + using ClusterType = + Cluster; + auto class_name = fmt::format("ClusterVector_{}", typestr); + + py::class_, void>>( + m, class_name.c_str(), + py::buffer_protocol()) + + .def(py::init()) // TODO change!!! + + .def("push_back", + [](ClusterVector &self, const ClusterType &cluster) { + self.push_back(cluster); + }) + + .def("sum", [](ClusterVector &self) { + auto *vec = new std::vector(self.sum()); + return return_vector(vec); + }) + .def_property_readonly("size", &ClusterVector::size) + .def("item_size", &ClusterVector::item_size) + .def_property_readonly("fmt", + [typestr](ClusterVector &self) { + return fmt_format; + }) + + .def_property_readonly("cluster_size_x", + &ClusterVector::cluster_size_x) + .def_property_readonly("cluster_size_y", + &ClusterVector::cluster_size_y) + .def_property_readonly("capacity", + &ClusterVector::capacity) + .def_property("frame_number", &ClusterVector::frame_number, + &ClusterVector::set_frame_number) + .def_buffer( + [typestr](ClusterVector &self) -> py::buffer_info { + return py::buffer_info( + self.data(), /* Pointer to buffer */ + self.item_size(), /* Size of one scalar */ + fmt_format, /* Format descriptor */ + 1, /* Number of dimensions */ + {self.size()}, /* Buffer dimensions */ + {self.item_size()} /* Strides (in bytes) for each index */ + ); + }); + + // Free functions using ClusterVector + m.def("hitmap", + [](std::array image_size, ClusterVector &cv) { + + // Create a numpy array to hold the hitmap + // The shape of the array is (image_size[0], image_size[1]) + // note that the python array is passed as [row, col] which + // is the opposite of the clusters [x,y] + py::array_t hitmap(image_size); + auto r = hitmap.mutable_unchecked<2>(); + + // Initialize hitmap to 0 + for (py::ssize_t i = 0; i < r.shape(0); i++) + for (py::ssize_t j = 0; j < r.shape(1); j++) + r(i, j) = 0; + + + // Loop over the clusters and increment the hitmap + // Skip out of bound clusters + for (const auto& cluster : cv) { + auto x = cluster.x; + auto y = cluster.y; + if(x -void define_cluster_vector(py::module &m, const std::string &typestr) { - using ClusterType = - Cluster; - auto class_name = fmt::format("ClusterVector_{}", typestr); - py::class_, void>>( - m, class_name.c_str(), - py::buffer_protocol()) - - .def(py::init()) // TODO change!!! - - .def("push_back", - [](ClusterVector &self, const ClusterType &cluster) { - self.push_back(cluster); - }) - - // implement push_back - .def_property_readonly("size", &ClusterVector::size) - .def("item_size", &ClusterVector::item_size) - .def_property_readonly("fmt", - [typestr](ClusterVector &self) { - return fmt_format; - }) - - .def_property_readonly("cluster_size_x", - &ClusterVector::cluster_size_x) - .def_property_readonly("cluster_size_y", - &ClusterVector::cluster_size_y) - .def_property_readonly("capacity", - &ClusterVector::capacity) - .def_property("frame_number", &ClusterVector::frame_number, - &ClusterVector::set_frame_number) - .def_buffer( - [typestr](ClusterVector &self) -> py::buffer_info { - return py::buffer_info( - self.data(), /* Pointer to buffer */ - self.item_size(), /* Size of one scalar */ - fmt_format, /* Format descriptor */ - 1, /* Number of dimensions */ - {self.size()}, /* Buffer dimensions */ - {self.item_size()} /* Strides (in bytes) for each index */ - ); - }); -} template @@ -252,25 +206,5 @@ void define_cluster_finder_bindings(py::module &m, const std::string &typestr) { }, py::arg(), py::arg("frame_number") = 0); - m.def("hitmap", - [](std::array image_size, ClusterVector &cv) { - py::array_t hitmap(image_size); - auto r = hitmap.mutable_unchecked<2>(); - - // Initialize hitmap to 0 - for (py::ssize_t i = 0; i < r.shape(0); i++) - for (py::ssize_t j = 0; j < r.shape(1); j++) - r(i, j) = 0; - - size_t stride = cv.item_size(); - auto ptr = cv.data(); - for (size_t i = 0; i < cv.size(); i++) { - auto x = *reinterpret_cast(ptr); - auto y = *reinterpret_cast(ptr + sizeof(int16_t)); - r(y, x) += 1; - ptr += stride; - } - return hitmap; - }); } #pragma GCC diagnostic pop diff --git a/python/src/module.cpp b/python/src/module.cpp index 8d5b5ab..c2067ed 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -1,4 +1,9 @@ // Files with bindings to the different classes + +//New style file naming +#include "bind_ClusterVector.hpp" + +//TODO! migrate the other names #include "cluster.hpp" #include "cluster_file.hpp" #include "ctb_raw_file.hpp" @@ -39,12 +44,12 @@ PYBIND11_MODULE(_aare, m) { define_cluster_file_io_bindings(m, "Cluster2x2f"); define_cluster_file_io_bindings(m, "Cluster2x2d"); - define_cluster_vector(m, "Cluster3x3i"); - define_cluster_vector(m, "Cluster3x3d"); - define_cluster_vector(m, "Cluster3x3f"); - define_cluster_vector(m, "Cluster2x2i"); - define_cluster_vector(m, "Cluster2x2d"); - define_cluster_vector(m, "Cluster2x2f"); + define_ClusterVector(m, "Cluster3x3i"); + define_ClusterVector(m, "Cluster3x3d"); + define_ClusterVector(m, "Cluster3x3f"); + define_ClusterVector(m, "Cluster2x2i"); + define_ClusterVector(m, "Cluster2x2d"); + define_ClusterVector(m, "Cluster2x2f"); define_cluster_finder_bindings(m, "Cluster3x3i"); define_cluster_finder_bindings(m, "Cluster3x3d"); diff --git a/python/tests/test_Cluster.py b/python/tests/test_Cluster.py index e24bcf8..ddaa6f3 100644 --- a/python/tests/test_Cluster.py +++ b/python/tests/test_Cluster.py @@ -1,12 +1,12 @@ import pytest import numpy as np -import aare._aare as aare +from aare import _aare #import the C++ module from conftest import test_data_path def test_cluster_vector_can_be_converted_to_numpy(): - cv = aare.ClusterVector_Cluster3x3i() + cv = _aare.ClusterVector_Cluster3x3i() arr = np.array(cv, copy=False) assert arr.shape == (0,) # 4 for x, y, size, energy and 9 for the cluster data @@ -14,24 +14,23 @@ def test_cluster_vector_can_be_converted_to_numpy(): def test_ClusterVector(): """Test ClusterVector""" - clustervector = aare.ClusterVector_Cluster3x3i() + clustervector = _aare.ClusterVector_Cluster3x3i() assert clustervector.cluster_size_x == 3 assert clustervector.cluster_size_y == 3 assert clustervector.item_size() == 4+9*4 assert clustervector.frame_number == 0 - assert clustervector.capacity == 1024 assert clustervector.size == 0 - cluster = aare.Cluster3x3i(0,0,np.ones(9, dtype=np.int32)) + cluster = _aare.Cluster3x3i(0,0,np.ones(9, dtype=np.int32)) clustervector.push_back(cluster) assert clustervector.size == 1 with pytest.raises(TypeError): # Or use the appropriate exception type - clustervector.push_back(aare.Cluster2x2i(0,0,np.ones(4, dtype=np.int32))) + clustervector.push_back(_aare.Cluster2x2i(0,0,np.ones(4, dtype=np.int32))) with pytest.raises(TypeError): - clustervector.push_back(aare.Cluster3x3f(0,0,np.ones(9, dtype=np.float32))) + clustervector.push_back(_aare.Cluster3x3f(0,0,np.ones(9, dtype=np.float32))) def test_Interpolator(): """Test Interpolator""" @@ -41,13 +40,13 @@ def test_Interpolator(): ybins = np.linspace(0, 5, 30, dtype=np.float64) etacube = np.zeros(shape=[30, 30, 20], dtype=np.float64) - interpolator = aare.Interpolator(etacube, xbins, ybins, ebins) + interpolator = _aare.Interpolator(etacube, xbins, ybins, ebins) assert interpolator.get_ietax().shape == (30,30,20) assert interpolator.get_ietay().shape == (30,30,20) - clustervector = aare.ClusterVector_Cluster3x3i() + clustervector = _aare.ClusterVector_Cluster3x3i() - cluster = aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32)) + cluster = _aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32)) clustervector.push_back(cluster) interpolated_photons = interpolator.interpolate(clustervector) @@ -58,9 +57,9 @@ def test_Interpolator(): assert interpolated_photons[0]["y"] == -1 assert interpolated_photons[0]["energy"] == 4 #eta_sum = 4, dx, dy = -1,-1 m_ietax = 0, m_ietay = 0 - clustervector = aare.ClusterVector_Cluster2x2i() + clustervector = _aare.ClusterVector_Cluster2x2i() - cluster = aare.Cluster2x2i(0,0, np.ones(4, dtype=np.int32)) + cluster = _aare.Cluster2x2i(0,0, np.ones(4, dtype=np.int32)) clustervector.push_back(cluster) interpolated_photons = interpolator.interpolate(clustervector) @@ -71,28 +70,15 @@ def test_Interpolator(): assert interpolated_photons[0]["y"] == 0 assert interpolated_photons[0]["energy"] == 4 -@pytest.mark.files -def test_cluster_file(test_data_path): - """Test ClusterFile""" - cluster_file = aare.ClusterFile_Cluster3x3i(test_data_path / "clust/single_frame_97_clustrers.clust") - clustervector = cluster_file.read_clusters(10) #conversion does not work - cluster_file.close() - - assert clustervector.size == 10 - - ###reading with wrong file - with pytest.raises(TypeError): - cluster_file = aare.ClusterFile_Cluster2x2i(test_data_path / "clust/single_frame_97_clustrers.clust") - cluster_file.close() def test_calculate_eta(): """Calculate Eta""" - clusters = aare.ClusterVector_Cluster3x3i() - clusters.push_back(aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32))) - clusters.push_back(aare.Cluster3x3i(0,0, np.array([1,1,1,2,2,2,3,3,3]))) + clusters = _aare.ClusterVector_Cluster3x3i() + clusters.push_back(_aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32))) + clusters.push_back(_aare.Cluster3x3i(0,0, np.array([1,1,1,2,2,2,3,3,3]))) - eta2 = aare.calculate_eta2(clusters) + eta2 = _aare.calculate_eta2(clusters) assert eta2.shape == (2,2) assert eta2[0,0] == 0.5 @@ -103,7 +89,7 @@ def test_calculate_eta(): def test_cluster_finder(): """Test ClusterFinder""" - clusterfinder = aare.ClusterFinder_Cluster3x3i([100,100]) + clusterfinder = _aare.ClusterFinder_Cluster3x3i([100,100]) #frame = np.random.rand(100,100) frame = np.zeros(shape=[100,100]) @@ -115,18 +101,7 @@ def test_cluster_finder(): assert clusters.size == 0 -#TODO dont understand behavior -def test_cluster_collector(): - """Test ClusterCollector""" - - clusterfinder = aare.ClusterFinderMT_Cluster3x3i([100,100]) #TODO: no idea what the data is in InputQueue not zero - - clustercollector = aare.ClusterCollector_Cluster3x3i(clusterfinder) - - cluster_vectors = clustercollector.steal_clusters() - - assert len(cluster_vectors) == 1 #single thread execution - assert cluster_vectors[0].size == 0 # + diff --git a/python/tests/test_ClusterFile.py b/python/tests/test_ClusterFile.py new file mode 100644 index 0000000..4126a6c --- /dev/null +++ b/python/tests/test_ClusterFile.py @@ -0,0 +1,64 @@ + +import pytest +import numpy as np +import boost_histogram as bh +import time +from pathlib import Path +import pickle + +from aare import ClusterFile +from conftest import test_data_path + +@pytest.mark.files +def test_cluster_file(test_data_path): + """Test ClusterFile""" + f = ClusterFile(test_data_path / "clust/single_frame_97_clustrers.clust") + cv = f.read_clusters(10) #conversion does not work + + + assert cv.frame_number == 135 + assert cv.size == 10 + + #Known data + #frame_number, num_clusters [135] 97 + #[ 1 200] [0 1 2 3 4 5 6 7 8] + #[ 2 201] [ 9 10 11 12 13 14 15 16 17] + #[ 3 202] [18 19 20 21 22 23 24 25 26] + #[ 4 203] [27 28 29 30 31 32 33 34 35] + #[ 5 204] [36 37 38 39 40 41 42 43 44] + #[ 6 205] [45 46 47 48 49 50 51 52 53] + #[ 7 206] [54 55 56 57 58 59 60 61 62] + #[ 8 207] [63 64 65 66 67 68 69 70 71] + #[ 9 208] [72 73 74 75 76 77 78 79 80] + #[ 10 209] [81 82 83 84 85 86 87 88 89] + + #conversion to numpy array + arr = np.array(cv, copy = False) + + assert arr.size == 10 + for i in range(10): + assert arr[i]['x'] == i+1 + +@pytest.mark.files +def test_read_clusters_and_fill_histogram(test_data_path): + # Create the histogram + n_bins = 100 + xmin = -100 + xmax = 1e4 + hist_aare = bh.Histogram(bh.axis.Regular(n_bins, xmin, xmax)) + + fname = test_data_path / "clust/beam_En700eV_-40deg_300V_10us_d0_f0_100.clust" + + #Read clusters and fill the histogram with pixel values + with ClusterFile(fname, chunk_size = 10000) as f: + for clusters in f: + arr = np.array(clusters, copy = False) + hist_aare.fill(arr['data'].flat) + + + #Load the histogram from the pickle file + with open(fname.with_suffix('.pkl'), 'rb') as f: + hist_py = pickle.load(f) + + #Compare the two histograms + assert hist_aare == hist_py \ No newline at end of file diff --git a/python/tests/test_ClusterVector.py b/python/tests/test_ClusterVector.py new file mode 100644 index 0000000..b64aeef --- /dev/null +++ b/python/tests/test_ClusterVector.py @@ -0,0 +1,54 @@ +import pytest +import numpy as np +import boost_histogram as bh +import time +from pathlib import Path +import pickle + +from aare import ClusterFile +from aare import _aare +from conftest import test_data_path + + +def test_create_cluster_vector(): + cv = _aare.ClusterVector_Cluster3x3i() + assert cv.cluster_size_x == 3 + assert cv.cluster_size_y == 3 + assert cv.size == 0 + + +def test_push_back_on_cluster_vector(): + cv = _aare.ClusterVector_Cluster2x2i() + assert cv.cluster_size_x == 2 + assert cv.cluster_size_y == 2 + assert cv.size == 0 + + cluster = _aare.Cluster2x2i(19, 22, np.ones(4, dtype=np.int32)) + cv.push_back(cluster) + assert cv.size == 1 + + arr = np.array(cv, copy=False) + assert arr[0]['x'] == 19 + assert arr[0]['y'] == 22 + + +def test_make_a_hitmap_from_cluster_vector(): + cv = _aare.ClusterVector_Cluster3x3i() + + # Push back 4 clusters with different positions + cv.push_back(_aare.Cluster3x3i(0, 0, np.ones(9, dtype=np.int32))) + cv.push_back(_aare.Cluster3x3i(1, 1, np.ones(9, dtype=np.int32))) + cv.push_back(_aare.Cluster3x3i(1, 1, np.ones(9, dtype=np.int32))) + cv.push_back(_aare.Cluster3x3i(2, 2, np.ones(9, dtype=np.int32))) + + ref = np.zeros((5, 5), dtype=np.int32) + ref[0,0] = 1 + ref[1,1] = 2 + ref[2,2] = 1 + + + img = _aare.hitmap((5,5), cv) + # print(img) + # print(ref) + assert (img == ref).all() + \ No newline at end of file diff --git a/src/CalculateEta.test.cpp b/src/CalculateEta.test.cpp index 2bdf387..29d7ed3 100644 --- a/src/CalculateEta.test.cpp +++ b/src/CalculateEta.test.cpp @@ -37,7 +37,7 @@ auto get_test_parameters() { Eta2{3. / 5, 4. / 6, 1, 11})); } -TEST_CASE("compute_largest_2x2_subcluster", "[.eta_calculation]") { +TEST_CASE("compute_largest_2x2_subcluster", "[eta_calculation]") { auto [cluster, expected_eta] = get_test_parameters(); auto [sum, index] = std::visit( @@ -47,7 +47,7 @@ TEST_CASE("compute_largest_2x2_subcluster", "[.eta_calculation]") { CHECK(expected_eta.sum == sum); } -TEST_CASE("calculate_eta2", "[.eta_calculation]") { +TEST_CASE("calculate_eta2", "[eta_calculation]") { auto [cluster, expected_eta] = get_test_parameters(); @@ -60,3 +60,69 @@ TEST_CASE("calculate_eta2", "[.eta_calculation]") { CHECK(eta.c == expected_eta.c); CHECK(eta.sum == expected_eta.sum); } + + +//3x3 cluster layout (in case of cBottomLeft etc corner): +// 6, 7, 8 +// 3, 4, 5 +// 0, 1, 2 + + +TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in the bottom left", + "[eta_calculation]") { + + // Create a 3x3 cluster + Cluster cl; + cl.x = 0; + cl.y = 0; + cl.data[0] = 30; + cl.data[1] = 23; + cl.data[2] = 5; + cl.data[3] = 20; + cl.data[4] = 50; + cl.data[5] = 3; + cl.data[6] = 8; + cl.data[7] = 2; + cl.data[8] = 3; + + // 8, 2, 3 + // 20, 50, 3 + // 30, 23, 5 + + auto eta = calculate_eta2(cl); + CHECK(eta.c == corner::cBottomLeft); + CHECK(eta.x == 50.0 / (20 + 50)); // 4/(3+4) + CHECK(eta.y == 50.0 / (23 + 50)); // 4/(1+4) + CHECK(eta.sum == 30+23+20+50); + +} + +TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in the top left", + "[eta_calculation]") { + +// Create a 3x3 cluster +Cluster cl; +cl.x = 0; +cl.y = 0; +cl.data[0] = 8; +cl.data[1] = 12; +cl.data[2] = 5; +cl.data[3] = 77; +cl.data[4] = 80; +cl.data[5] = 3; +cl.data[6] = 82; +cl.data[7] = 91; +cl.data[8] = 3; + +// 82, 91, 3 +// 77, 80, 3 +// 8, 12, 5 + +auto eta = calculate_eta2(cl); +CHECK(eta.c == corner::cTopLeft); +CHECK(eta.x == 80. / (77 + 80)); // 4/(3+4) +CHECK(eta.y == 91.0 / (91 + 80)); // 7/(7+4) +CHECK(eta.sum == 77+80+82+91); + +} +