WIP

include/aare/ClusterFile.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "aare/defs.hpp"
+#include "aare/ClusterVector.hpp"
 #include <filesystem>
 #include <fstream>
 
@@ -38,30 +39,40 @@ struct ClusterAnalysis {
     double etay;
 };
 
-
+/*
-
+Binary cluster file. Expects data to be layed out as:
+int32_t frame_number
+uint32_t number_of_clusters
+int16_t x, int16_t y, int32_t data[9] x number_of_clusters
+int32_t frame_number
+uint32_t number_of_clusters
+....
+*/
 class ClusterFile {
     FILE *fp{};
     uint32_t m_num_left{};
     size_t m_chunk_size{};
+    const std::string m_mode;
 
   public:
-    ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000);
+    ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
+                const std::string &mode = "r");
     ~ClusterFile();
     std::vector<Cluster> read_clusters(size_t n_clusters);
     std::vector<Cluster> read_frame(int32_t &out_fnum);
+    void write_frame(int32_t frame_number, const ClusterVector<int32_t>& clusters);
     std::vector<Cluster>
     read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny);
 
     int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
-                 double *eta2x, double *eta2y, double *eta3x, double *eta3y);
+                     double *eta2x, double *eta2y, double *eta3x,
+                     double *eta3y);
     int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad,
                         double *eta2x, double *eta2y, double *eta3x,
                         double *eta3y);
 
     size_t chunk_size() const { return m_chunk_size; }
     void close();
-
 };
 
 } // namespace aare

include/aare/ClusterVector.hpp

@@ -6,14 +6,25 @@
 #include <fmt/core.h>
 
 template <typename T> class ClusterVector {
-    int32_t m_cluster_size_x;
+    using value_type = T;
-    int32_t m_cluster_size_y;
+    using coord_t = int16_t;
+    coord_t m_cluster_size_x;
+    coord_t m_cluster_size_y;
     std::byte *m_data{};
     size_t m_size{0};
     size_t m_capacity;
+    /*
+    Format string used in the python bindings to create a numpy
+    array from the buffer
+    = - native byte order
+    h - short
+    d - double
+    i - int
+    */
+    constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:" ;
 
   public:
-    ClusterVector(int32_t cluster_size_x, int32_t cluster_size_y,
+    ClusterVector(coord_t cluster_size_x, coord_t cluster_size_y,
                   size_t capacity = 1024)
         : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
           m_capacity(capacity) {
@@ -25,15 +36,8 @@ template <typename T> class ClusterVector {
         delete[] m_data;
     }
 
-    // ClusterVector(const ClusterVector & other){
-    //     m_cluster_size_x = other.m_cluster_size_x;
-    //     m_cluster_size_y = other.m_cluster_size_y;
-    //     m_data = new std::byte[other.m_capacity];
-    //     std::copy(other.m_data, other.m_data + other.m_capacity, m_data);
-    //     m_size = other.m_size;
-    //     m_capacity = other.m_capacity;
-    // }
    
+    //Move constructor
     ClusterVector(ClusterVector &&other) noexcept
         : m_cluster_size_x(other.m_cluster_size_x),
           m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data),
@@ -66,15 +70,15 @@ template <typename T> class ClusterVector {
     }
 
     // data better hold data of the right size!
-    void push_back(int32_t x, int32_t y, const std::byte *data) {
+    void push_back(coord_t x, coord_t y, const std::byte *data) {
         if (m_size == m_capacity) {
             allocate_buffer(m_capacity * 2);
         }
         std::byte *ptr = element_ptr(m_size);
-        *reinterpret_cast<int32_t *>(ptr) = x;
+        *reinterpret_cast<coord_t *>(ptr) = x;
-        ptr += sizeof(int32_t);
+        ptr += sizeof(coord_t);
-        *reinterpret_cast<int32_t *>(ptr) = y;
+        *reinterpret_cast<coord_t *>(ptr) = y;
-        ptr += sizeof(int32_t);
+        ptr += sizeof(coord_t);
 
         std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
                   ptr);
@@ -85,7 +89,7 @@ template <typename T> class ClusterVector {
         std::vector<T> sums(m_size);
         const size_t stride = element_offset();
         const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
-        std::byte *ptr = m_data + 2 * sizeof(int32_t); // skip x and y
+        std::byte *ptr = m_data + 2 * sizeof(coord_t); // skip x and y
 
         for (size_t i = 0; i < m_size; i++) {
             sums[i] =
@@ -109,6 +113,12 @@ template <typename T> class ClusterVector {
     int16_t cluster_size_y() const { return m_cluster_size_y; }
 
     std::byte *data() { return m_data; }
+    const std::byte *data() const { return m_data; }
+
+    const std::string_view fmt_base() const {
+        //TODO! how do we match on coord_t?
+        return m_fmt_base;
+    }
 
   private:
     void allocate_buffer(size_t new_capacity) {

include/aare/File.hpp

@@ -44,6 +44,7 @@ class File {
     void read_into(std::byte *image_buf);
     void read_into(std::byte *image_buf, size_t n_frames);
     
+    size_t frame_number();                    //!< get the frame number at the current position
     size_t frame_number(size_t frame_index);  //!< get the frame number at the given frame index
     size_t bytes_per_frame() const;           
     size_t pixels_per_frame() const;    

python/src/cluster.hpp

@@ -21,25 +21,25 @@ void define_cluster_vector(py::module &m, const std::string &typestr) {
         .def("element_offset",
              py::overload_cast<>(&ClusterVector<T>::element_offset, py::const_))
         .def_property_readonly("fmt",
-             [typestr](ClusterVector<T> &v) {
+             [typestr](ClusterVector<T> &self) {
                  return fmt::format(
-                     "i:x:\ni:y:\n({},{}){}:data:", v.cluster_size_x(),
+                     self.fmt_base(), self.cluster_size_x(),
-                     v.cluster_size_y(), typestr);
+                     self.cluster_size_y(), typestr);
              })
         .def("sum", [](ClusterVector<T> &self) {
             auto *vec = new std::vector<T>(self.sum());
             return return_vector(vec);
         })
-        .def_buffer([typestr](ClusterVector<T> &v) -> py::buffer_info {
+        .def_buffer([typestr](ClusterVector<T> &self) -> py::buffer_info {
             return py::buffer_info(
-                v.data(),           /* Pointer to buffer */
+                self.data(),           /* Pointer to buffer */
-                v.element_offset(), /* Size of one scalar */
+                self.element_offset(), /* Size of one scalar */
-                fmt::format("i:x:\ni:y:\n{}{}:data:", v.cluster_size_x()*
+                fmt::format(self.fmt_base(), self.cluster_size_x(),
-                            v.cluster_size_y(),
+                            self.cluster_size_y(),
                             typestr), /* Format descriptor */
                 1,                    /* Number of dimensions */
-                {v.size()},           /* Buffer dimensions */
+                {self.size()},           /* Buffer dimensions */
-                {v.element_offset()}  /* Strides (in bytes) for each index */
+                {self.element_offset()}  /* Strides (in bytes) for each index */
             );
         });
 }

python/src/cluster_file.hpp

@@ -1,7 +1,6 @@
 #include "aare/ClusterFile.hpp"
 #include "aare/defs.hpp"
 
-
 #include <cstdint>
 #include <filesystem>
 #include <pybind11/iostream.h>
@@ -18,33 +17,47 @@ void define_cluster_file_io_bindings(py::module &m) {
     PYBIND11_NUMPY_DTYPE(Cluster, x, y, data);
 
     py::class_<ClusterFile>(m, "ClusterFile")
-        .def(py::init<const std::filesystem::path &, size_t>(), py::arg(), py::arg("chunk_size") = 1000)
+        .def(py::init<const std::filesystem::path &, size_t,
+                      const std::string &>(),
+             py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r")
         .def("read_clusters",
              [](ClusterFile &self, size_t n_clusters) {
-                 auto* vec = new std::vector<Cluster>(self.read_clusters(n_clusters));
+                 auto *vec =
+                     new std::vector<Cluster>(self.read_clusters(n_clusters));
                  return return_vector(vec);
              })
         .def("read_frame",
              [](ClusterFile &self) {
                  int32_t frame_number;
-                 auto* vec = new std::vector<Cluster>(self.read_frame(frame_number));
+                 auto *vec =
+                     new std::vector<Cluster>(self.read_frame(frame_number));
                  return py::make_tuple(frame_number, return_vector(vec));
              })
+        .def("write_frame", &ClusterFile::write_frame)
         .def("read_cluster_with_cut",
-             [](ClusterFile &self, size_t n_clusters, py::array_t<double> noise_map, int nx, int ny) {
+             [](ClusterFile &self, size_t n_clusters,
+                py::array_t<double> noise_map, int nx, int ny) {
                  auto view = make_view_2d(noise_map);
-                 auto* vec = new std::vector<Cluster>(self.read_cluster_with_cut(n_clusters, view.data(), nx, ny));
+                 auto *vec =
+                     new std::vector<Cluster>(self.read_cluster_with_cut(
+                         n_clusters, view.data(), nx, ny));
                  return return_vector(vec);
              })
         .def("__enter__", [](ClusterFile &self) { return &self; })
-        .def("__exit__", [](ClusterFile &self) {  self.close();})
+        .def("__exit__",
+             [](ClusterFile &self,
+                const std::optional<pybind11::type> &exc_type,
+                const std::optional<pybind11::object> &exc_value,
+                const std::optional<pybind11::object> &traceback) {
+                 self.close();
+             })
         .def("__iter__", [](ClusterFile &self) { return &self; })
         .def("__next__", [](ClusterFile &self) {
-            auto vec =  new std::vector<Cluster>(self.read_clusters(self.chunk_size()));
+            auto vec =
+                new std::vector<Cluster>(self.read_clusters(self.chunk_size()));
             if (vec->size() == 0) {
                 throw py::stop_iteration();
             }
             return return_vector(vec);
         });
-
 }

python/src/file.hpp

@@ -51,7 +51,8 @@ void define_file_io_bindings(py::module &m) {
         .def(py::init<const std::filesystem::path &, const std::string &,
                       const FileConfig &>())
 
-        .def("frame_number", &File::frame_number)
+        .def("frame_number", py::overload_cast<>(&File::frame_number))
+        .def("frame_number", py::overload_cast<size_t>(&File::frame_number))
         .def_property_readonly("bytes_per_frame", &File::bytes_per_frame)
         .def_property_readonly("pixels_per_frame", &File::pixels_per_frame)
         .def("seek", &File::seek)

src/ClusterFile.cpp

@@ -2,17 +2,28 @@
 
 namespace aare {
 
-ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size): m_chunk_size(chunk_size) {
+ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
+                         const std::string &mode)
+    : m_chunk_size(chunk_size), m_mode(mode) {
+
+    if (mode == "r") {
         fp = fopen(fname.c_str(), "rb");
         if (!fp) {
-        throw std::runtime_error("Could not open file: " + fname.string());
+            throw std::runtime_error("Could not open file for reading: " + fname.string());
         }
+    } else if (mode == "w") {
+        fp = fopen(fname.c_str(), "wb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for writing: " + fname.string());
+        }
+    } else {
+        throw std::runtime_error("Unsupported mode: " + mode);
     }
 
-ClusterFile::~ClusterFile() {
-    close();
 }
 
+ClusterFile::~ClusterFile() { close(); }
+
 void ClusterFile::close() {
     if (fp) {
         fclose(fp);
@@ -20,7 +31,25 @@ void ClusterFile::close(){
     }
 }
 
+void ClusterFile::write_frame(int32_t frame_number, const ClusterVector<int32_t>& clusters){
+    if (m_mode != "w") {
+        throw std::runtime_error("File not opened for writing");
+    }
+    if(!(clusters.cluster_size_x()==3) && !(clusters.cluster_size_y()==3)){
+        throw std::runtime_error("Only 3x3 clusters are supported");
+    }
+    fwrite(&frame_number, sizeof(frame_number), 1, fp);
+    uint32_t n_clusters = clusters.size();
+    fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
+    fwrite(clusters.data(), clusters.element_offset(), clusters.size(), fp);
+    // write clusters
+    // fwrite(clusters.data(), sizeof(Cluster), clusters.size(), fp);
+}
+
 std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
     std::vector<Cluster> clusters(n_clusters);
 
     int32_t iframe = 0; // frame number needs to be 4 bytes!
@@ -38,7 +67,8 @@ std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
         } else {
             nn = nph;
         }
-        nph_read += fread(reinterpret_cast<void *>(buf + nph_read), sizeof(Cluster), nn, fp);
+        nph_read += fread(reinterpret_cast<void *>(buf + nph_read),
+                          sizeof(Cluster), nn, fp);
         m_num_left = nph - nn; // write back the number of photons left
     }
 
@@ -52,8 +82,8 @@ std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
                 else
                     nn = nph;
 
-                nph_read +=
+                nph_read += fread(reinterpret_cast<void *>(buf + nph_read),
-                    fread(reinterpret_cast<void *>(buf + nph_read), sizeof(Cluster), nn, fp);
+                                  sizeof(Cluster), nn, fp);
                 m_num_left = nph - nn;
             }
             if (nph_read >= n_clusters)
@@ -68,8 +98,12 @@ std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
 }
 
 std::vector<Cluster> ClusterFile::read_frame(int32_t &out_fnum) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
     if (m_num_left) {
-        throw std::runtime_error("There are still photons left in the last frame");
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
     }
 
     if (fread(&out_fnum, sizeof(out_fnum), 1, fp) != 1) {
@@ -82,17 +116,19 @@ std::vector<Cluster> ClusterFile::read_frame(int32_t &out_fnum) {
     }
     std::vector<Cluster> clusters(n_clusters);
 
-    if (fread(clusters.data(), sizeof(Cluster), n_clusters, fp) != static_cast<size_t>(n_clusters)) {
+    if (fread(clusters.data(), sizeof(Cluster), n_clusters, fp) !=
+        static_cast<size_t>(n_clusters)) {
         throw std::runtime_error("Could not read clusters");
     }
     return clusters;
-
 }
 
 std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
                                                         double *noise_map,
                                                         int nx, int ny) {
-
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
     std::vector<Cluster> clusters(n_clusters);
     // size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
     //                               uint32_t *n_left, double *noise_map, int
@@ -124,7 +160,8 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
         }
         for (size_t iph = 0; iph < nn; iph++) {
             // read photons 1 by 1
-            size_t n_read = fread(reinterpret_cast<void *>(ptr), sizeof(Cluster), 1, fp);
+            size_t n_read =
+                fread(reinterpret_cast<void *>(ptr), sizeof(Cluster), 1, fp);
             if (n_read != 1) {
                 clusters.resize(nph_read);
                 return clusters;
@@ -159,7 +196,8 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
         }
     }
     if (nph_read < n_clusters) {
-    //         // keep on reading frames and photons until reaching n_clusters
+        //         // keep on reading frames and photons until reaching
+        //         n_clusters
         while (fread(&iframe, sizeof(iframe), 1, fp)) {
             //             // printf("%d\n",nph_read);
 
@@ -168,8 +206,8 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
                 m_num_left = nph;
                 for (size_t iph = 0; iph < nph; iph++) {
                     //                     // read photons 1 by 1
-                        size_t n_read =
+                    size_t n_read = fread(reinterpret_cast<void *>(ptr),
-                            fread(reinterpret_cast<void *>(ptr), sizeof(Cluster), 1, fp);
+                                          sizeof(Cluster), 1, fp);
                     if (n_read != 1) {
                         clusters.resize(nph_read);
                         return clusters;
@@ -180,8 +218,7 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
                         if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 &&
                             ptr->y < ny) {
                             tot1 = ptr->data[4];
-                                analyze_cluster(*ptr, &t2max, &tot3, NULL,
+                            analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL,
-                                NULL,
                                             NULL, NULL, NULL);
                             // noise = noise_map[ptr->y * nx + ptr->x];
                             noise = noise_map[ptr->y + ny * ptr->x];
@@ -191,8 +228,8 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
                             } else
                                 good = 0;
                         } else {
-                                printf("Bad pixel number %d %d\n", ptr->x,
+                            printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
-                                ptr->y); good = 0;
+                            good = 0;
                         }
                     }
                     if (good) {
@@ -211,18 +248,18 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
     // printf("%d\n",nph_read);
     clusters.resize(nph_read);
     return clusters;
-
 }
 
-int ClusterFile::analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad,
+int ClusterFile::analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3,
-                    double *eta2x, double *eta2y, double *eta3x,
+                                 char *quad, double *eta2x, double *eta2y,
-                    double *eta3y) {
+                                 double *eta3x, double *eta3y) {
 
     return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
 }
 
-int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
+int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3,
-                 double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
+                              char *quad, double *eta2x, double *eta2y,
+                              double *eta3x, double *eta3y) {
 
     int ok = 1;
 
@@ -263,7 +300,8 @@ int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *qua
                 c = i;
             }
         }
-	//printf("*** %d %d %d %d -- %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
+        // printf("*** %d %d %d %d --
+        // %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
         if (quad)
             *quad = c;
         if (t2)
@@ -318,6 +356,4 @@ int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *qua
     return ok;
 }
 
-
-
 } // namespace aare

src/File.cpp

@@ -58,6 +58,8 @@ void File::read_into(std::byte *image_buf) { file_impl->read_into(image_buf); }
 void File::read_into(std::byte *image_buf, size_t n_frames) {
     file_impl->read_into(image_buf, n_frames);
 }
+
+size_t File::frame_number() { return file_impl->frame_number(tell()); }
 size_t File::frame_number(size_t frame_index) {
     return file_impl->frame_number(frame_index);
 }