Repository: BlueBrain/libsonata (language: python)

Issue:
Retrieve spike time units in SpikePopulation
I noticed there is a "units" attribute of the `/spikes/All/timestamps` dataset in the SONATA spike report HDF5 file, but this attribute cannot be queried using `libsonata`.

I don't know how much this attribute is used, or whether values other than "ms" are allowed, but if there's a semantic for time units in the format, I think the library should support it.

A bit of context, this came up while writing some code to concatenate spike files. Time units across multiple files cannot be verified to be the same (using `libsonata`), which is unfortunate.


Discussion / hints:
Usually @jorblancoa or @alkino take care of the report reading, so perhaps they can have a look.

Fix patch:
diff --git a/include/bbp/sonata/report_reader.h b/include/bbp/sonata/report_reader.h
index 5b4cd5e1..2aa74f35 100644
--- a/include/bbp/sonata/report_reader.h
+++ b/include/bbp/sonata/report_reader.h
@@ -72,6 +72,11 @@ class SONATA_API SpikeReader
          */
         Sorting getSorting() const;
 
+        /**
+         * Return the unit of time
+         */
+        std::string getTimeUnits() const;
+
       private:
         Population(const std::string& filename, const std::string& populationName);
 
@@ -79,6 +84,7 @@ class SONATA_API SpikeReader
         Sorting sorting_ = Sorting::none;
         // Use for clamping of user values
         double tstart_, tstop_;
+        std::string time_units_;
 
         void filterNode(Spikes& spikes, const Selection& node_ids) const;
         void filterTimestamp(Spikes& spikes, double tstart, double tstop) const;
diff --git a/python/bindings.cpp b/python/bindings.cpp
index bf800e84..29c460d7 100644
--- a/python/bindings.cpp
+++ b/python/bindings.cpp
@@ -1283,7 +1283,10 @@ PYBIND11_MODULE(_libsonata, m) {
             DOC_SPIKEREADER_POP(getSorting))
         .def_property_readonly("times",
                                &SpikeReader::Population::getTimes,
-                               DOC_SPIKEREADER_POP(getTimes));
+                               DOC_SPIKEREADER_POP(getTimes))
+        .def_property_readonly("time_units",
+                               &SpikeReader::Population::getTimeUnits,
+                               DOC_REPORTREADER_POP(getTimeUnits));
     py::class_<SpikeReader>(m, "SpikeReader", DOC(bbp, sonata, SpikeReader))
         .def(py::init([](py::object h5_filepath) { return SpikeReader(py::str(h5_filepath)); }),
              "h5_filepath"_a)
diff --git a/python/generated/docstrings.h b/python/generated/docstrings.h
index db259a47..51fb0182 100644
--- a/python/generated/docstrings.h
+++ b/python/generated/docstrings.h
@@ -1324,6 +1324,8 @@ static const char *__doc_bbp_sonata_SpikeReader_Population_getSorting = R"doc(Re
 
 static const char *__doc_bbp_sonata_SpikeReader_Population_getTimes = R"doc(Return (tstart, tstop) of the population)doc";
 
+static const char *__doc_bbp_sonata_SpikeReader_Population_getTimeUnits = R"doc(Return the unit of time)doc";
+
 static const char *__doc_bbp_sonata_SpikeReader_Population_sorting = R"doc()doc";
 
 static const char *__doc_bbp_sonata_SpikeReader_Population_spike_times = R"doc()doc";
diff --git a/src/report_reader.cpp b/src/report_reader.cpp
index 1b639b59..2cd8349e 100644
--- a/src/report_reader.cpp
+++ b/src/report_reader.cpp
@@ -196,6 +196,10 @@ SpikeReader::Population::Sorting SpikeReader::Population::getSorting() const {
     return sorting_;
 }
 
+std::string SpikeReader::Population::getTimeUnits() const {
+    return time_units_;
+}
+
 SpikeReader::Population::Population(const std::string& filename,
                                     const std::string& populationName) {
     HighFive::File file(filename, HighFive::File::ReadOnly);
@@ -206,6 +210,7 @@ SpikeReader::Population::Population(const std::string& filename,
 
     pop.getDataSet("node_ids").read(node_ids);
     pop.getDataSet("timestamps").read(timestamps);
+    pop.getDataSet("timestamps").getAttribute("units").read(time_units_);
 
     if (node_ids.size() != timestamps.size()) {
         throw SonataError(


Test patch:
diff --git a/python/tests/test_reports.py b/python/tests/test_reports.py
index 2624a0dc..d67bec65 100644
--- a/python/tests/test_reports.py
+++ b/python/tests/test_reports.py
@@ -44,6 +44,7 @@ def test_get_spikes_from_population(self):
         self.assertEqual(self.test_obj['All'].sorting, "by_time")
         self.assertEqual(self.test_obj['spikes1'].sorting, "by_id")
         self.assertEqual(self.test_obj['spikes2'].sorting, "none")
+        self.assertEqual(self.test_obj['spikes2'].time_units, 'ms')
         self.assertEqual(self.test_obj['empty'].get(), [])
 
         self.assertEqual(len(self.test_obj['All'].get(node_ids=[])), 0)
diff --git a/tests/data/generate.py b/tests/data/generate.py
index e03e798c..0ea7c228 100755
--- a/tests/data/generate.py
+++ b/tests/data/generate.py
@@ -204,13 +204,15 @@ def write_spikes(filepath):
         gpop_all = h5f.create_group('/spikes/' + population_names[0])
         gpop_all.attrs.create('sorting', data=2, dtype=sorting_type)
         timestamps, node_ids = zip(*sorted(zip(timestamps_base, node_ids_base)))
-        set = gpop_all.create_dataset('timestamps', data=timestamps, dtype=np.double)
+        dtimestamps = gpop_all.create_dataset('timestamps', data=timestamps, dtype=np.double)
+        dtimestamps.attrs.create('units', data="ms", dtype=string_dtype)
         gpop_all.create_dataset('node_ids', data=node_ids, dtype=np.uint64)
 
         gpop_spikes1 = h5f.create_group('/spikes/' + population_names[1])
         gpop_spikes1.attrs.create('sorting', data=1, dtype=sorting_type)
         node_ids, timestamps = zip(*sorted(zip(node_ids_base, timestamps_base)))
-        gpop_spikes1.create_dataset('timestamps', data=timestamps, dtype=np.double)
+        dtimestamps = gpop_spikes1.create_dataset('timestamps', data=timestamps, dtype=np.double)
+        dtimestamps.attrs.create('units', data="ms", dtype=string_dtype)
         gpop_spikes1.create_dataset('node_ids', data=node_ids, dtype=np.uint64)
 
         gpop_spikes2 = h5f.create_group('/spikes/' + population_names[2])
@@ -222,6 +224,7 @@ def write_spikes(filepath):
         gpop_empty = h5f.create_group('/spikes/' + population_names[3])
         gpop_empty.attrs.create('sorting', data=1, dtype=sorting_type)
         dtimestamps = gpop_empty.create_dataset('timestamps', data=[], dtype=np.double)
+        dtimestamps.attrs.create('units', data="ms", dtype=string_dtype)
         gpop_empty.create_dataset('node_ids', data=[], dtype=np.uint64)
 
 if __name__ == '__main__':
diff --git a/tests/data/spikes.h5 b/tests/data/spikes.h5
index b80cac94..6ede470d 100644
Binary files a/tests/data/spikes.h5 and b/tests/data/spikes.h5 differ
diff --git a/tests/test_report_reader.cpp b/tests/test_report_reader.cpp
index b269d63e..5acdf849 100644
--- a/tests/test_report_reader.cpp
+++ b/tests/test_report_reader.cpp
@@ -35,6 +35,7 @@ TEST_CASE("SpikeReader", "[base]") {
     REQUIRE(reader.openPopulation("empty").get() == std::vector<std::pair<uint64_t, double>>{});
 
     REQUIRE(reader.openPopulation("All").getTimes() == std::make_tuple(0.1, 1.3));
+    REQUIRE(reader.openPopulation("All").getTimeUnits() == "ms");
 }
 
 TEST_CASE("SomaReportReader limits", "[base]") {