Merge pull request #45 from Fluorescence-Tools/development

V0.25

.github/workflows/conda-build.yml

@@ -24,9 +24,7 @@ jobs:
       - name: Set up Conda
         uses: conda-incubator/setup-miniconda@v3
         with:
-          auto-update-conda: true
-          python-version: ${{ matrix.python-version }}
-          mamba-version: "*"
+          miniconda-version: "latest"
           channels: conda-forge,defaults
           channel-priority: true
 
@@ -46,8 +44,8 @@ jobs:
       - name: Build Conda Package
         shell: bash -el {0}
         run: |
-          mamba install conda-build boa
+          conda install conda-build boa mamba
           cd conda-recipe
-          curl -sLO  https://raw.githubusercontent.com/conda-forge/conda-forge-pinning-feedstock/084b098a28a7a66e9a0967d156bc55b9ebfc6726/recipe/conda_build_config.yaml
-          # curl -sLO https://github.com/conda-forge/conda-forge-pinning-feedstock/raw/main/recipe/conda_build_config.yaml
+          # curl -sLO  https://raw.githubusercontent.com/conda-forge/conda-forge-pinning-feedstock/084b098a28a7a66e9a0967d156bc55b9ebfc6726/recipe/conda_build_config.yaml
+          curl -sLO https://github.com/conda-forge/conda-forge-pinning-feedstock/raw/main/recipe/conda_build_config.yaml
           conda mambabuild .

.github/workflows/release.yml

@@ -24,9 +24,7 @@ jobs:
       - name: Set up Conda
         uses: conda-incubator/setup-miniconda@v3
         with:
-          auto-update-conda: true
-          python-version: ${{ matrix.python-version }}
-          mamba-version: "*"
+          miniconda-version: "latest"
           channels: conda-forge,defaults
           channel-priority: true
 
@@ -46,7 +44,7 @@ jobs:
       - name: Build and Upload Conda Package
         shell: bash -el {0}
         run: |
-          mamba install conda-build boa anaconda-client
+          conda install conda-build mamba boa anaconda-client
           cd conda-recipe
           # Use older conda_build_config -> boost-cpp 1.78
           # curl -sLO  https://raw.githubusercontent.com/conda-forge/conda-forge-pinning-feedstock/084b098a28a7a66e9a0967d156bc55b9ebfc6726/recipe/conda_build_config.yaml

conda-recipe/meta.yaml

@@ -24,7 +24,7 @@ requirements:
     - ninja
     - make               # [linux]
     - doxygen
-    - swig
+    - swig 4.2.*
     - python
     - numpy
   host:

doc/whats_new/changes.rst

@@ -2,14 +2,17 @@
 
 .. _changes_0_23:
 
-Version 0.24
+Version 0.25
 ============
+* Add simpler option for burst search
+* Added get_supported_container_names
 
+Version 0.24
+-----------
 * Add support for sm files
 * Add support for CZ CF3 FCS files
 * Improved type tttr file type inference to mitigate crashes
 
-
 Version 0.23
 ------------
 * Consider background with certain fraction in mean lifetime (0.23.6)

ext/python/TTTR.py

@@ -10,19 +10,22 @@ def event_types(self):
 
 def __getattr__(self, item):
     """
-    If an attribute `attribute` is accesses that does not exist
-    the corresponding getter by calling 'get_attribute' is called
-
-    :param self:
-    :param item:
-    :return:
+    If an attribute `attribute` is accessed that does not exist,
+    the corresponding getter method ('get_attribute') is called.
+    Works for both instance and static methods.
     """
     item = "get_" + str(item)
+    # Check if the static method or instance method exists in the class
     if hasattr(self.__class__, item):
-        call = getattr(self, item)
-        return call()
+        call = getattr(self.__class__, item)
+        if isinstance(call, staticmethod):
+            # If it's a static method, call it directly from the class
+            return call.__get__(None, self.__class__)()
+        else:
+            # Otherwise, assume it's an instance method
+            return call(self)
     else:
-        raise AttributeError
+        raise AttributeError(f"'{self.__class__.__name__}' object has no attribute '{item}'")
 
 def __len__(self):
     return self.get_n_valid_events()

ext/python/misc_types.i

@@ -33,8 +33,9 @@ import_array();
 %template(VectorUint128) std::vector<unsigned long long>;
 %template(VectorUint32_3D) std::vector<std::vector<std::vector<unsigned int>>>;
 %template(VectorDouble_2D) std::vector<std::vector<double>>;
-%template(MapStringString) std::map<std::string, std::string>;
 %template(MapShortVectorDouble) std::map<short, std::vector<double>>;
+%template(MapStringString) std::map<std::string, std::string>;
+%template(VectorString) std::vector<std::string>;
 
 %typemap(out) std::vector< double,std::allocator< double > > * {
 $result = swig::from(static_cast<std::vector< double,std::allocator< double > >>(*($1)));

include/DecayFit.h

@@ -107,7 +107,7 @@ struct DecayFitIntegrateSignals{
     double Fp(){
         double g = 1.0;
         if(corrections != nullptr){
-            g = corrections->g;
+            g = corrections->gamma;
         }
         if(g == 1.0){
             return (Sp - Bp);
@@ -117,15 +117,23 @@ struct DecayFitIntegrateSignals{
     }
 
     double Fs(){
-        double g = 1.0;
+        double g = 1.0, r;
         if(corrections != nullptr){
-            g = corrections->g;
+            g = corrections->gamma;
         }
         if(g == 1.0){
-            return (Ss - Bs);
+            r = (Ss - Bs);
         } else{
-            return (Ss - g * Bs) / (1. - g);
+            r = (Ss - g * Bs) / (1. - g);
         }
+#ifdef VERBOSE_TTTRLIB
+        std::cout << "Fs()" << std::endl;
+        std::cout << "g:" << g << std::endl;
+        std::cout << "Ss:" << Ss << std::endl;
+        std::cout << "Bs:" << Bs << std::endl;
+        std::cout << "Fs:" << r << std::endl;
+#endif
+        return r;
     }
 
     double r(){
@@ -137,6 +145,15 @@ struct DecayFitIntegrateSignals{
             l1 = corrections->l1;
             l2 = corrections->l2;
         }
+
+#ifdef VERBOSE_TTTRLIB
+        std::cout << "fp:" << fp << std::endl;
+        std::cout << "fs:" << fs << std::endl;
+        std::cout << "g:" << g << std::endl;
+        std::cout << "l1:" << l1 << std::endl;
+        std::cout << "l2:" << l2 << std::endl;
+#endif
+
         double nom = (fp - g * fs);
         double denom = (fp * (1. - 3. * l2) + (2. - 3. * l1) * g * fs);
         return nom / denom;

include/TTTR.h

@@ -211,7 +211,21 @@ class TTTR : public std::enable_shared_from_this<TTTR>{
     TTTRHeader *header = nullptr;
 
     /// map to translates string container types to int container types
-    boost::bimap<std::string, int> container_names = {};
+    static boost::bimap<std::string, int> container_names;
+
+    // Static function to initialize the container_names
+    static boost::bimap<std::string, int> initialize_container_names() {
+        boost::bimap<std::string, int> m;
+        m.insert({std::string("PTU"), PQ_PTU_CONTAINER});
+        m.insert({std::string("HT3"), PQ_HT3_CONTAINER});
+        m.insert({std::string("SPC-130"), BH_SPC130_CONTAINER});
+        m.insert({std::string("SPC-600_256"), BH_SPC600_256_CONTAINER});
+        m.insert({std::string("SPC-600_4096"), BH_SPC600_4096_CONTAINER});
+        m.insert({std::string("PHOTON-HDF5"), PHOTON_HDF_CONTAINER});
+        m.insert({std::string("CZ-RAW"), CZ_CONFOCOR3_CONTAINER});
+        m.insert({std::string("SM"), SM_CONTAINER});
+        return m;
+    }
 
     typedef bool (*processRecord_t)(
             uint32_t&,  // input
@@ -536,16 +550,36 @@ class TTTR : public std::enable_shared_from_this<TTTR>{
       return get_n_valid_events();
      }
 
-     /*!
-      * \brief Retrieves the used routing channel numbers from the TTTR data.
-      *
-      * This function populates the provided output array with the routing channel
-      * numbers that are used in the TTTR file. The number of elements in the output
-      * array is stored in the n_output parameter.
-      *
-      * @param output Pointer to the output array to be populated.
-      * @param n_output Pointer to the number of elements in the output array.
-      */
+    /**
+    * Sliding window burst search.
+    *
+    * Finds bursts in the macro time array. A burst starts when the photon rate
+    * is above a minimum threshold, and ends when the rate falls below the same
+    * threshold. The rate-threshold is defined by the ratio `m`/`T` (`m` photons
+    * in a time interval `T`). A burst is discarded if it has less than `L`
+    * photons.
+    *
+    * Arguments:
+    *     L (int): minimum number of photons in a burst. Bursts with size
+    *         (or counts) < L are discarded.
+    *     m (int): number of consecutive photons used to compute the rate.
+    *     T (double): max time separation of `m` photons to be inside a burst (in seconds).
+    *
+    * Returns:
+    *     vector<int64_t>: A vector of interleaved start and stop indices.
+    */
+    std::vector<long long> burst_search(int L, int m, double T);
+
+    /*!
+     * \brief Retrieves the used routing channel numbers from the TTTR data.
+     *
+     * This function populates the provided output array with the routing channel
+     * numbers that are used in the TTTR file. The number of elements in the output
+     * array is stored in the n_output parameter.
+     *
+     * @param output Pointer to the output array to be populated.
+     * @param n_output Pointer to the number of elements in the output array.
+     */
      void get_used_routing_channels(signed char **output, int *n_output);
 
      /*!
@@ -640,6 +674,26 @@ class TTTR : public std::enable_shared_from_this<TTTR>{
         return tttr_container_type_str;
     }
 
+    /**
+     * @brief Retrieves a list of supported container names.
+     *
+     * This function returns a list of supported container names, e.g. 'PTU', 'SPC-130', etc..
+     *
+     * @return std::vector<std::string>
+     * A list of supported container names.
+     *
+     */
+     static std::vector<std::string> get_supported_container_names(){
+
+        // Extract container names from the right side of the bimap
+        std::vector<std::string> supported_container_names;
+        for (const auto& element : container_names.left) {
+            supported_container_names.push_back(element.first);
+        }
+
+        return supported_container_names;
+     }
+
      /*!
       * \brief Creates a new TTTR object by selecting specific events based on the provided indices.
       *

pyproject.toml

@@ -1,5 +1,5 @@
 [project]
-version = "0.24.4"
+version = "0.25.0"
 name = "tttrlib"
 requires-python = ">=3.8"
 description = "Read, write & process time-tagged time-resolved (TTTR) data."

src/DecayFit.cpp

@@ -27,7 +27,7 @@ void DecayFitIntegrateSignals::compute_signal_and_background(MParam *p) {
     std::cout << "-- Nchannels_exp:" << Nchannels_exp << std::endl;
     std::cout << "-- Bp, Bs: " << Bp << ", " << Bs << std::endl;
     std::cout << "-- Sp, Ss: " << Sp << ", " << Ss << std::endl;
-    std::cout << "-- Bexpected, Bs: " << Bp << ", " << Bs << std::endl;
+    std::cout << "-- Bexpected: " << Bexpected << std::endl;
 #endif
 }
 

src/TTTR.cpp

@@ -1,5 +1,8 @@
 #include "TTTR.h"
 
+// Static member definition outside the class
+boost::bimap<std::string, int> TTTR::container_names = TTTR::initialize_container_names();
+
 
 TTTR::TTTR() :
         // private
@@ -20,14 +23,6 @@ TTTR::TTTR() :
         n_records_read(0),
         n_valid_events(0),
         processRecord(nullptr){
-    container_names.insert({std::string("PTU"), PQ_PTU_CONTAINER});
-    container_names.insert({std::string("HT3"), PQ_HT3_CONTAINER});
-    container_names.insert({std::string("SPC-130"), BH_SPC130_CONTAINER});
-    container_names.insert({std::string("SPC-600_256"), BH_SPC600_256_CONTAINER});
-    container_names.insert({std::string("SPC-600_4096"), BH_SPC600_4096_CONTAINER});
-    container_names.insert({std::string("PHOTON-HDF5"), PHOTON_HDF_CONTAINER});
-    container_names.insert({std::string("CZ-RAW"), CZ_CONFOCOR3_CONTAINER});
-    container_names.insert({std::string("SM"), SM_CONTAINER});
     header = new TTTRHeader(tttr_container_type);
     allocate_memory_for_records(0);
 }
@@ -122,8 +117,9 @@ TTTR::TTTR(const TTTR &p2){
 
 TTTR::TTTR(const char *fn, int container_type, bool read_input) : TTTR(){
     if(container_type >= 0){
-        filename.assign(fn);
+        tttr_container_type_str = container_names.right.at(container_type);
         tttr_container_type = container_type;
+        filename.assign(fn);
         if(read_input){
             if(read_file())
                 find_used_routing_channels();
@@ -133,17 +129,6 @@ TTTR::TTTR(const char *fn, int container_type, bool read_input) : TTTR(){
     }
 }
 
-TTTR::TTTR(const char *fn, int container_type) : TTTR(fn, container_type, true) {
-    try {
-        tttr_container_type_str.assign(
-                container_names.right.at(container_type)
-        );
-    }
-    catch(...) {
-        std::cerr << "TTTR::TTTR(const char *fn, int container_type): Container type " << container_type << " not supported." << std::endl;
-    }
-}
-
 TTTR::TTTR(const char *fn, const char *container_type) : TTTR() {
     try {
         tttr_container_type_str.assign(container_type);
@@ -157,6 +142,17 @@ TTTR::TTTR(const char *fn, const char *container_type) : TTTR() {
     }
 }
 
+TTTR::TTTR(const char *fn, int container_type) : TTTR(fn, container_type, true) {
+    try {
+        tttr_container_type_str.assign(
+                container_names.right.at(container_type)
+        );
+    }
+    catch(...) {
+        std::cerr << "TTTR::TTTR(const char *fn, int container_type): Container type " << container_type << " not supported." << std::endl;
+    }
+}
+
 TTTR::TTTR(const char* filename) : TTTR(filename, inferTTTRFileType(filename), true) {}
 
 
@@ -752,6 +748,42 @@ void selection_by_count_rate(
     }
 }
 
+
+std::vector<long long> TTTR::burst_search(int L, int m, double T) {
+
+    int64_t i, i_start, i_stop;
+    uint8_t in_burst = 0;
+
+    std::vector<long long> bursts;  // Now storing interleaved start and stop indices
+    long long Ti = T / header->get_macro_time_resolution();
+
+    for (i = 0; i <= size() - m; ++i) {
+        if (macro_times[i + m - 1] - macro_times[i] <= Ti) {
+            if (!in_burst) {
+                in_burst = 1;
+                i_start = i;
+            }
+        } else if (in_burst) {
+            in_burst = 0;
+            i_stop = i + m - 2;
+            if (i_stop - i_start + 1 >= L) {
+                bursts.push_back(i_start);
+                bursts.push_back(i_stop);
+            }
+        }
+    }
+
+    if (in_burst) {
+        i_stop = i + m - 1;
+        if (i_stop - i_start + 1 >= L) {
+            bursts.push_back(i_start);
+            bursts.push_back(i_stop);
+        }
+    }
+
+    return bursts;
+}
+
 unsigned int TTTR::get_number_of_micro_time_channels(){
     return header->get_effective_number_of_micro_time_channels();
 }