Repository: BlueBrain/libsonata (language: python)

Issue:
Improve mem. ops. + NodeID-ElementID mapping bindings
The PR improves the copy of the returned `ids` when fetching the data from the reports, as well as the Python bindings related to the NodeID-ElementID mapping (https://github.com/BlueBrain/libsonata/pull/168):

- In the former case, the `std::vector` is [swapped in a constant-time operation instead of a full copy](https://www.cplusplus.com/reference/vector/vector/swap). This is particularly useful when a large number of GIDs is requested. 
- In the latter case, the binding of the NodeID-ElementID mapping from is improved, returning a `numpy.array` instead of a list of tuples. Hence, providing directly the result generated from the C++ code without any additional copies.

Using the same ~60TB report from https://github.com/BlueBrain/libsonata/pull/174 that contains 4.2M GIDs (1.5B compartments), the following table illustrates the performance obtained by fetching a single timestep and varying the number of GIDs requested:

```
 # GIDS |            v0.1.11 |                 PR
 ------- -------------------- -------------------
      1 |  0.010879278182983 |  0.001698255538940
     10 |  5.343438386917114 |  5.512796163558960
    100 |  5.943392753601074 |  5.604890823364258
   1000 |  6.406604766845703 |  6.062043905258179
  10000 |  6.451685428619385 |  6.049137830734253
 100000 |  7.101217269897461 |  6.462664365768433
1000000 | 12.762017250061035 |  9.955524444580078
4234929 | 27.868659019470215 | 18.158316850662230
```

On the other hand, when using in Python the function `get_node_id_element_id_mapping()` from https://github.com/BlueBrain/libsonata/pull/168, the following table shows the performance obtained varying the number of GIDs requested:

```
 # GIDS |            v0.1.11 |                 PR
 ------- -------------------- -------------------
      1 |    0.0158734321594 |  0.001111268997192
     10 |    2.6325545310974 |  2.509858608245849
    100 |    2.9563860893249 |  2.782548427581787
   1000 |    3.2792141437530 |  2.988783836364746
  10000 |    5.4664847850799 |  3.832525491714477
 100000 |   24.0672097206115 |  3.506523370742798
1000000 |  225.1766004562378 |  6.213826656341553
4234929 | 1017.6884262561798 | 11.914409875869751
```

Note that there has been other attempts to improve the copies of the voltages and even the tuples for the NodeID-ElementID mapping container, but the benefits turned out to be negligible compared to the source code complexity that was required.


Improve mem. ops. + NodeID-ElementID mapping bindings
The PR improves the copy of the returned `ids` when fetching the data from the reports, as well as the Python bindings related to the NodeID-ElementID mapping (https://github.com/BlueBrain/libsonata/pull/168):

- In the former case, the `std::vector` is [swapped in a constant-time operation instead of a full copy](https://www.cplusplus.com/reference/vector/vector/swap). This is particularly useful when a large number of GIDs is requested. 
- In the latter case, the binding of the NodeID-ElementID mapping from is improved, returning a `numpy.array` instead of a list of tuples. Hence, providing directly the result generated from the C++ code without any additional copies.

Using the same ~60TB report from https://github.com/BlueBrain/libsonata/pull/174 that contains 4.2M GIDs (1.5B compartments), the following table illustrates the performance obtained by fetching a single timestep and varying the number of GIDs requested:

```
 # GIDS |            v0.1.11 |                 PR
 ------- -------------------- -------------------
      1 |  0.010879278182983 |  0.001698255538940
     10 |  5.343438386917114 |  5.512796163558960
    100 |  5.943392753601074 |  5.604890823364258
   1000 |  6.406604766845703 |  6.062043905258179
  10000 |  6.451685428619385 |  6.049137830734253
 100000 |  7.101217269897461 |  6.462664365768433
1000000 | 12.762017250061035 |  9.955524444580078
4234929 | 27.868659019470215 | 18.158316850662230
```

On the other hand, when using in Python the function `get_node_id_element_id_mapping()` from https://github.com/BlueBrain/libsonata/pull/168, the following table shows the performance obtained varying the number of GIDs requested:

```
 # GIDS |            v0.1.11 |                 PR
 ------- -------------------- -------------------
      1 |    0.0158734321594 |  0.001111268997192
     10 |    2.6325545310974 |  2.509858608245849
    100 |    2.9563860893249 |  2.782548427581787
   1000 |    3.2792141437530 |  2.988783836364746
  10000 |    5.4664847850799 |  3.832525491714477
 100000 |   24.0672097206115 |  3.506523370742798
1000000 |  225.1766004562378 |  6.213826656341553
4234929 | 1017.6884262561798 | 11.914409875869751
```

Note that there has been other attempts to improve the copies of the voltages and even the tuples for the NodeID-ElementID mapping container, but the benefits turned out to be negligible compared to the source code complexity that was required.


Fix/assertion error
Fixes an assertion error thrown when `_GLIBCXX_ASSERTIONS` is defined, and enables the option by default


Discussion / hints:
Thank you very much, @matz-e and @mgeplf!
Thanks for the change
Thank you very much, @matz-e and @mgeplf!
Thanks for the change
I'm not sure what the corresponding flags for clang are, but I see that the CI builds completed fine so perhaps these extra definitions are simply ignored?

Fix patch:
diff --git a/.gitignore b/.gitignore
index 9e1bbb6d..b79074d4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,10 +1,12 @@
 .cache
 .clangd
 build
+spack-build*
 .eggs/
 .tox/
 libsonata.cpython-*.so
 python/__pycache__/
+python/libsonata.egg-info/
 *.pyc
 venv/
 cpp/
diff --git a/include/bbp/sonata/report_reader.h b/include/bbp/sonata/report_reader.h
index 2baa051f..80ae7c1b 100644
--- a/include/bbp/sonata/report_reader.h
+++ b/include/bbp/sonata/report_reader.h
@@ -28,16 +28,6 @@ struct SONATA_API DataFrame {
 using Spike = std::pair<NodeID, double>;
 using Spikes = std::vector<Spike>;
 
-/**
-  const SpikeReader file(filename);
-  auto pops = file.getPopulationNames();
-  for (const auto& data: file[pops.openPopulation(0).get(Selection{12UL, 34UL, 58UL})) {
-      NodeID node_id;
-      double timestamp;
-      std::tie(node_id, timestamp) = data;
-      std::cout << "[" << timestamp << "] " << node_id << std::endl;
-  }
-*/
 class SONATA_API SpikeReader
 {
   public:
@@ -135,14 +125,16 @@ class SONATA_API ReportReader
          * - For Soma report reader, the return value will be the Node ID to which the report
          *   value belongs to.
          * - For Element/full compartment readers, the return value will be an array with 2
-         *   elements, the first element is the Node ID and the second element is the 
+         *   elements, the first element is the Node ID and the second element is the
          *   compartment ID of the given Node.
          *
          * \param node_ids limit the report to the given selection. If nullptr, all nodes in the
          * report are used
+         * \param block_gap_limit gap limit between each IO block while fetching data from storage
          */
         typename DataFrame<KeyType>::DataType getNodeIdElementIdMapping(
-            const nonstd::optional<Selection>& node_ids = nonstd::nullopt) const;
+            const nonstd::optional<Selection>& node_ids = nonstd::nullopt,
+            const nonstd::optional<size_t>& block_gap_limit = nonstd::nullopt) const;
 
         /**
          * \param node_ids limit the report to the given selection.
@@ -150,17 +142,22 @@ class SONATA_API ReportReader
          * indicates no limit. \param tstop return voltages occurring on or before tstop.
          * tstop=nonstd::nullopt indicates no limit. \param tstride indicates every how many
          * timesteps we read data. tstride=nonstd::nullopt indicates that all timesteps are read.
+         * \param block_gap_limit gap limit between each IO block while fetching data from storage.
          */
-        DataFrame<KeyType> get(const nonstd::optional<Selection>& node_ids = nonstd::nullopt,
-                               const nonstd::optional<double>& tstart = nonstd::nullopt,
-                               const nonstd::optional<double>& tstop = nonstd::nullopt,
-                               const nonstd::optional<size_t>& tstride = nonstd::nullopt) const;
+        DataFrame<KeyType> get(
+            const nonstd::optional<Selection>& node_ids = nonstd::nullopt,
+            const nonstd::optional<double>& tstart = nonstd::nullopt,
+            const nonstd::optional<double>& tstop = nonstd::nullopt,
+            const nonstd::optional<size_t>& tstride = nonstd::nullopt,
+            const nonstd::optional<size_t>& block_gap_limit = nonstd::nullopt) const;
 
       private:
         struct NodeIdElementLayout {
             typename DataFrame<KeyType>::DataType ids;
             Selection::Ranges node_ranges;
-            Selection::Range min_max_range;
+            std::vector<uint64_t> node_offsets;
+            std::vector<uint64_t> node_index;
+            Selection::Ranges min_max_blocks;
         };
 
         Population(const H5::File& file, const std::string& populationName);
@@ -173,18 +170,22 @@ class SONATA_API ReportReader
          *
          * \param node_ids limit the report to the given selection. If nullptr, all nodes in the
          * report are used
+         * \param block_gap_limit gap limit between each IO block while fetching data from storage
          */
         NodeIdElementLayout getNodeIdElementLayout(
-            const nonstd::optional<Selection>& node_ids = nonstd::nullopt) const;
+            const nonstd::optional<Selection>& node_ids = nonstd::nullopt,
+            const nonstd::optional<size_t>& block_gap_limit = nonstd::nullopt) const;
 
-        std::map<NodeID, Selection::Range> node_ranges_;
         H5::Group pop_group_;
-        std::vector<NodeID> nodes_ids_;
+        std::vector<NodeID> node_ids_;
+        std::vector<Selection::Range> node_ranges_;
+        std::vector<uint64_t> node_offsets_;
+        std::vector<uint64_t> node_index_;
         double tstart_, tstop_, tstep_;
         std::vector<std::pair<size_t, double>> times_index_;
         std::string time_units_;
         std::string data_units_;
-        bool nodes_ids_sorted_ = false;
+        bool is_node_ids_sorted_;
 
         friend ReportReader;
     };
diff --git a/python/bindings.cpp b/python/bindings.cpp
index d8cab8e8..ed07225f 100644
--- a/python/bindings.cpp
+++ b/python/bindings.cpp
@@ -358,18 +358,21 @@ void bindReportReader(py::module& m, const std::string& prefix) {
              "node_ids"_a = nonstd::nullopt,
              "tstart"_a = nonstd::nullopt,
              "tstop"_a = nonstd::nullopt,
-             "tstride"_a = nonstd::nullopt)
+             "tstride"_a = nonstd::nullopt,
+             "block_gap_limit"_a = nonstd::nullopt)
         .def("get_node_ids",
              &ReportType::Population::getNodeIds,
              "Return the list of nodes ids for this population")
         .def(
             "get_node_id_element_id_mapping",
             [](const typename ReportType::Population& population,
-               const nonstd::optional<Selection>& selection) {
-                return asArray(population.getNodeIdElementIdMapping(selection));
+               const nonstd::optional<Selection>& selection,
+               const nonstd::optional<size_t>& block_gap_limit) {
+                return asArray(population.getNodeIdElementIdMapping(selection, block_gap_limit));
             },
             DOC_REPORTREADER_POP(getNodeIdElementIdMapping),
-            "selection"_a = nonstd::nullopt)
+            "selection"_a = nonstd::nullopt,
+            "block_gap_limit"_a = nonstd::nullopt)
         .def_property_readonly("sorted",
                                &ReportType::Population::getSorted,
                                DOC_REPORTREADER_POP(getSorted))
diff --git a/src/report_reader.cpp b/src/report_reader.cpp
index 34098ea1..4f27f388 100644
--- a/src/report_reader.cpp
+++ b/src/report_reader.cpp
@@ -2,7 +2,7 @@
 #include <fmt/format.h>
 
 #include <algorithm>  // std::copy, std::find, std::lower_bound, std::upper_bound
-#include <iterator>   // std::next
+#include <iterator>   // std::advance, std::next
 
 constexpr double EPSILON = 1e-6;
 
@@ -81,17 +81,13 @@ void filterTimestampSorted(Spikes& spikes, double tstart, double tstop) {
     spikes.erase(spikes.begin(), begin);
 }
 
-template <typename T>
-T make_key(NodeID node_id, ElementID element_id);
-
-template <>
-NodeID make_key(NodeID node_id, ElementID /* element_id */) {
-    return node_id;
+inline void emplace_ids(NodeID& key, NodeID node_id, ElementID /* element_id */) {
+    key = node_id;
 }
 
-template <>
-CompartmentID make_key(NodeID node_id, ElementID element_id) {
-    return {node_id, element_id};
+inline void emplace_ids(CompartmentID& key, NodeID node_id, ElementID element_id) {
+    key[0] = node_id;
+    key[1] = element_id;
 }
 
 }  // anonymous namespace
@@ -225,36 +221,55 @@ auto ReportReader<T>::openPopulation(const std::string& populationName) const ->
 
 template <typename T>
 ReportReader<T>::Population::Population(const H5::File& file, const std::string& populationName)
-    : pop_group_(file.getGroup(std::string("/report/") + populationName)) {
-    {
-        const auto mapping_group = pop_group_.getGroup("mapping");
-        mapping_group.getDataSet("node_ids").read(nodes_ids_);
+    : pop_group_(file.getGroup(std::string("/report/") + populationName))
+    , is_node_ids_sorted_(false) {
+    const auto mapping_group = pop_group_.getGroup("mapping");
+    mapping_group.getDataSet("node_ids").read(node_ids_);
 
-        std::vector<uint64_t> index_pointers;
-        mapping_group.getDataSet("index_pointers").read(index_pointers);
+    std::vector<uint64_t> index_pointers;
+    mapping_group.getDataSet("index_pointers").read(index_pointers);
 
-        for (size_t i = 0; i < nodes_ids_.size(); ++i) {
-            node_ranges_.emplace(nodes_ids_[i],
-                                 std::make_pair(index_pointers[i], index_pointers[i + 1]));
-        }
+    if (index_pointers.size() != (node_ids_.size() + 1)) {
+        throw SonataError("'index_pointers' dataset size must be 'node_ids' size plus one");
+    }
 
-        {  // Get times
-            std::vector<double> times;
-            mapping_group.getDataSet("time").read(times);
-            tstart_ = times[0];
-            tstop_ = times[1];
-            tstep_ = times[2];
-            mapping_group.getDataSet("time").getAttribute("units").read(time_units_);
-            size_t i = 0;
-            for (double t = tstart_; t < tstop_ - EPSILON; t += tstep_, ++i) {
-                times_index_.emplace_back(i, t);
-            }
-        }
+    // Expand the pointers into tuples that define the range of each GID
+    size_t element_ids_count = 0;
+    for (size_t i = 0; i < node_ids_.size(); ++i) {
+        node_ranges_.emplace_back(index_pointers[i], index_pointers[i + 1]);  // Range of GID
+        node_offsets_.emplace_back(element_ids_count);                        // Offset in output
+        node_index_.emplace_back(i);                                          // Index of previous
 
+        element_ids_count += (index_pointers[i + 1] - index_pointers[i]);
+    }
+    node_offsets_.emplace_back(element_ids_count);
+
+    {  // Sort the index according to the GIDs, if not sorted in file
         if (mapping_group.getDataSet("node_ids").hasAttribute("sorted")) {
-            uint8_t sorted;
+            uint8_t sorted = 0;
             mapping_group.getDataSet("node_ids").getAttribute("sorted").read(sorted);
-            nodes_ids_sorted_ = sorted != 0;
+            is_node_ids_sorted_ = (sorted != 0);
+        }
+
+        if (!is_node_ids_sorted_) {
+            // Note: The idea is to sort the positions to access the values, allowing us to
+            //       maintain all vectors intact, while still being able to index the data
+            std::sort(node_index_.begin(), node_index_.end(), [&](const size_t i, const size_t j) {
+                return node_ids_[i] < node_ids_[j];
+            });
+        }
+    }
+
+    {  // Get times
+        std::vector<double> times;
+        mapping_group.getDataSet("time").read(times);
+        tstart_ = times[0];
+        tstop_ = times[1];
+        tstep_ = times[2];
+        mapping_group.getDataSet("time").getAttribute("units").read(time_units_);
+        size_t i = 0;
+        for (double t = tstart_; t < tstop_ - EPSILON; t += tstep_, ++i) {
+            times_index_.emplace_back(i, t);
         }
     }
 
@@ -278,79 +293,112 @@ std::string ReportReader<T>::Population::getDataUnits() const {
 
 template <typename T>
 bool ReportReader<T>::Population::getSorted() const {
-    return nodes_ids_sorted_;
+    return is_node_ids_sorted_;
 }
 
 template <typename T>
 std::vector<NodeID> ReportReader<T>::Population::getNodeIds() const {
-    return nodes_ids_;
+    return node_ids_;
 }
 
 template <typename T>
 typename ReportReader<T>::Population::NodeIdElementLayout
 ReportReader<T>::Population::getNodeIdElementLayout(
-    const nonstd::optional<Selection>& node_ids) const {
-    std::vector<NodeID> concrete_node_ids;
+    const nonstd::optional<Selection>& node_ids,
+    const nonstd::optional<size_t>& _block_gap_limit) const {
     NodeIdElementLayout result;
-    result.min_max_range = {std::numeric_limits<uint64_t>::max(),
-                            std::numeric_limits<uint64_t>::min()};
+    std::vector<NodeID> concrete_node_ids;
     size_t element_ids_count = 0;
 
-    // Helper function to update layout, alongside the min / max values and count
-    const auto update_layout = [&](const NodeID& node_id, const Selection::Range& range) {
-        concrete_node_ids.emplace_back(node_id);
-        result.node_ranges.emplace_back(range);
+    // Set the gap between IO blocks while fetching data (Default: 64MB / 4 x GPFS blocks)
+    const size_t block_gap_limit = _block_gap_limit.value_or(16777216);
 
-        result.min_max_range.first = std::min(result.min_max_range.first, range.first);
-        result.min_max_range.second = std::max(result.min_max_range.second, range.second);
-        element_ids_count += (range.second - range.first);
-    };
+    if (block_gap_limit < 4194304) {
+        throw SonataError("block_gap_limit must be at least 4194304 (16MB / 1 x GPFS block)");
+    }
 
     // Take all nodes if no selection is provided
     if (!node_ids) {
-        concrete_node_ids.reserve(node_ranges_.size());
-        result.node_ranges.reserve(node_ranges_.size());
-
-        for (const auto& node_range : node_ranges_) {
-            update_layout(node_range.first, node_range.second);
-        }
+        concrete_node_ids = node_ids_;
+        result.node_ranges = node_ranges_;
+        result.node_offsets = node_offsets_;
+        result.node_index = node_index_;
+        element_ids_count = node_offsets_.back();
     } else if (!node_ids->empty()) {
         const auto selected_node_ids = node_ids->flatten();
 
-        // Reserve space for all the requested IDs and shrink afterwards
-        concrete_node_ids.reserve(selected_node_ids.size());
-        result.node_ranges.reserve(selected_node_ids.size());
-
         for (const auto node_id : selected_node_ids) {
-            const auto it = node_ranges_.find(node_id);
-            if (it != node_ranges_.end()) {
-                update_layout(it->first, it->second);
+            const auto it = std::lower_bound(node_index_.begin(),
+                                             node_index_.end(),
+                                             node_id,
+                                             [&](const size_t i, const NodeID node_id) {
+                                                 return node_ids_[i] < node_id;
+                                             });
+
+            if (it != node_index_.end() && node_ids_[*it] == node_id) {
+                const auto& range = node_ranges_[*it];
+
+                concrete_node_ids.emplace_back(node_id);
+                result.node_ranges.emplace_back(range);
+                result.node_offsets.emplace_back(element_ids_count);
+                result.node_index.emplace_back(result.node_index.size());
+
+                element_ids_count += (range.second - range.first);
             }
         }
-
-        concrete_node_ids.shrink_to_fit();
-        result.node_ranges.shrink_to_fit();
     } else {
         // node_ids Selection exists, but is empty
     }
 
     // Extract the ElementIDs from the GIDs
-    if (!result.node_ranges.empty()) {
-        const auto min = result.min_max_range.first;
-        const auto max = result.min_max_range.second;
+    if (!concrete_node_ids.empty()) {
+        // Sort the index by the selected ranges
+        std::sort(result.node_index.begin(),
+                  result.node_index.end(),
+                  [&](const size_t i, const size_t j) {
+                      return result.node_ranges[i].first < result.node_ranges[j].first;
+                  });
+
+        // Generate the {min,max} IO blocks for the requests
+        size_t offset = 0;
+        for (size_t i = 0; (i + 1) < result.node_index.size(); i++) {
+            const auto index = result.node_index[i];
+            const auto index_next = result.node_index[i + 1];
+            const auto max = result.node_ranges[index].second;
+            const auto min_next = result.node_ranges[index_next].first;
+
+            if ((min_next - max) > block_gap_limit) {
+                result.min_max_blocks.emplace_back(std::make_pair(offset, (i + 1)));
+                offset = (i + 1);
+            }
+        }
+        result.min_max_blocks.emplace_back(std::make_pair(offset, result.node_index.size()));
 
+        // Fill the GID-ElementID mapping in blocks to reduce the file system overhead
         std::vector<ElementID> element_ids;
         auto dataset_elem_ids = pop_group_.getGroup("mapping").getDataSet("element_ids");
-        dataset_elem_ids.select({min}, {max - min}).read(element_ids);
 
-        result.ids.reserve(element_ids_count);
+        result.ids.resize(element_ids_count);
+
+        for (const auto& min_max_block : result.min_max_blocks) {
+            const auto first_index = result.node_index[min_max_block.first];
+            const auto last_index = result.node_index[min_max_block.second - 1];
+            const auto min = result.node_ranges[first_index].first;
+            const auto max = result.node_ranges[last_index].second;
+
+            dataset_elem_ids.select({min}, {max - min}).read(element_ids);
 
-        for (size_t i = 0; i < concrete_node_ids.size(); ++i) {
-            const auto node_id = concrete_node_ids[i];
-            const auto& range = result.node_ranges[i];
+            // Copy the values for each of the GIDs assigned into this block
+            for (size_t i = min_max_block.first; i < min_max_block.second; ++i) {
+                const auto index = result.node_index[i];
+                const auto node_id = concrete_node_ids[index];
+                const auto range = Selection::Range(result.node_ranges[index].first - min,
+                                                    result.node_ranges[index].second - min);
 
-            for (auto i = (range.first - min); i < (range.second - min); i++) {
-                result.ids.emplace_back(make_key<T>(node_id, element_ids[i]));
+                auto offset = result.node_offsets[index];
+                for (auto i = range.first; i < range.second; i++, offset++) {
+                    emplace_ids(result.ids[offset], node_id, element_ids[i]);
+                }
             }
         }
 
@@ -359,7 +407,9 @@ ReportReader<T>::Population::getNodeIdElementLayout(
         //            another large range, the next IOps take an extra few seconds.
         //            We observed that fooling HDF5 hides the issue, but we should
         //            verify this behaviour once new releases of HDF5 are available.
-        dataset_elem_ids.select({min}, {1}).read(element_ids);
+        const auto min_max_block = result.min_max_blocks.back();
+        const auto index = result.node_index[min_max_block.first];
+        dataset_elem_ids.select({result.node_ranges[index].first}, {1}).read(element_ids);
     }
 
     return result;
@@ -402,15 +452,18 @@ std::pair<size_t, size_t> ReportReader<T>::Population::getIndex(
 
 template <typename T>
 typename DataFrame<T>::DataType ReportReader<T>::Population::getNodeIdElementIdMapping(
-    const nonstd::optional<Selection>& node_ids) const {
-    return getNodeIdElementLayout(node_ids).ids;
+    const nonstd::optional<Selection>& node_ids,
+    const nonstd::optional<size_t>& block_gap_limit) const {
+    return getNodeIdElementLayout(node_ids, block_gap_limit).ids;
 }
 
 template <typename T>
-DataFrame<T> ReportReader<T>::Population::get(const nonstd::optional<Selection>& node_ids,
-                                              const nonstd::optional<double>& tstart,
-                                              const nonstd::optional<double>& tstop,
-                                              const nonstd::optional<size_t>& tstride) const {
+DataFrame<T> ReportReader<T>::Population::get(
+    const nonstd::optional<Selection>& node_ids,
+    const nonstd::optional<double>& tstart,
+    const nonstd::optional<double>& tstop,
+    const nonstd::optional<size_t>& tstride,
+    const nonstd::optional<size_t>& block_gap_limit) const {
     size_t index_start = 0;
     size_t index_stop = 0;
     std::tie(index_start, index_stop) = getIndex(tstart, tstop);
@@ -422,12 +475,12 @@ DataFrame<T> ReportReader<T>::Population::get(const nonstd::optional<Selection>&
         throw SonataError("tstart should be <= to tstop");
     }
 
-    // min and max offsets of the node_ids requested are calculated
-    // to reduce the amount of IO that is brought to memory
-    auto node_id_element_layout = getNodeIdElementLayout(node_ids);
+    // Retrieve the GID-ElementID layout, alongside the {min,max} blocks
+    auto node_id_element_layout = getNodeIdElementLayout(node_ids, block_gap_limit);
     const auto& node_ranges = node_id_element_layout.node_ranges;
-    const auto min = node_id_element_layout.min_max_range.first;
-    const auto max = node_id_element_layout.min_max_range.second;
+    const auto& node_offsets = node_id_element_layout.node_offsets;
+    const auto& node_index = node_id_element_layout.node_index;
+    const auto& min_max_blocks = node_id_element_layout.min_max_blocks;
 
     if (node_id_element_layout.ids.empty()) {  // At the end no data available (wrong node_ids?)
         return DataFrame<T>{{}, {}, {}};
@@ -443,9 +496,9 @@ DataFrame<T> ReportReader<T>::Population::get(const nonstd::optional<Selection>&
     data_frame.ids.swap(node_id_element_layout.ids);
 
     // Fill .data member
-    size_t n_time_entries = ((index_stop - index_start) / stride) + 1;
-    size_t n_ids = data_frame.ids.size();
-    data_frame.data.resize(n_time_entries * n_ids);
+    const size_t n_time_entries = ((index_stop - index_start) / stride) + 1;
+    const size_t element_ids_count = data_frame.ids.size();
+    data_frame.data.resize(n_time_entries * element_ids_count);
 
     auto dataset = pop_group_.getDataSet("data");
     auto dataset_type = dataset.getDataType();
@@ -455,30 +508,39 @@ DataFrame<T> ReportReader<T>::Population::get(const nonstd::optional<Selection>&
                         dataset_type.string()));
     }
 
-    std::vector<float> buffer(max - min);
+    std::vector<float> buffer;
+    auto data_start = data_frame.data.begin();
     for (size_t timer_index = index_start; timer_index <= index_stop; timer_index += stride) {
-        // Note: The code assumes that the file is chunked by rows and not by columns
-        // (i.e., if the chunking changes in the future, the reading method must also be adapted)
-        dataset.select({timer_index, min}, {1, max - min}).read(buffer.data());
-
-        size_t offset = 0;
-        size_t data_offset = (timer_index - index_start) / stride * n_ids;
-        auto data_start = std::next(data_frame.data.begin(), data_offset);
-        for (const auto& range : node_ranges) {
-            uint64_t elements_per_gid = range.second - range.first;
-            uint64_t gid_start = range.first - min;
-
-            // Soma report
-            if (elements_per_gid == 1) {
-                data_start[offset] = buffer[gid_start];
-            } else {  // Elements report
-                uint64_t gid_end = range.second - min;
-                std::copy(std::next(buffer.begin(), gid_start),
-                          std::next(buffer.begin(), gid_end),
-                          std::next(data_start, offset));
+        // Access the data in blocks to reduce the file system overhead
+        for (const auto& min_max_block : min_max_blocks) {
+            const auto first_index = node_index[min_max_block.first];
+            const auto last_index = node_index[min_max_block.second - 1];
+            const auto min = node_ranges[first_index].first;
+            const auto max = node_ranges[last_index].second;
+
+            dataset.select({timer_index, min}, {1, max - min}).read(buffer);
+
+            // Copy the values for each of the GIDs assigned into this block
+            const auto buffer_start = buffer.begin();
+            for (size_t i = min_max_block.first; i < min_max_block.second; ++i) {
+                const auto index = node_index[i];
+                const auto range = Selection::Range(node_ranges[index].first - min,
+                                                    node_ranges[index].second - min);
+                const auto elements_per_gid = (range.second - range.first);
+                const auto offset = node_offsets[index];
+
+                // Soma report
+                if (elements_per_gid == 1) {
+                    data_start[offset] = buffer_start[range.first];
+                } else {  // Elements report
+                    std::copy(std::next(buffer_start, range.first),
+                              std::next(buffer_start, range.second),
+                              std::next(data_start, offset));
+                }
             }
-            offset += elements_per_gid;
         }
+
+        std::advance(data_start, element_ids_count);
     }
 
     return data_frame;


Test patch:
diff --git a/python/tests/test.py b/python/tests/test.py
index 4f38796f..d5662a5a 100644
--- a/python/tests/test.py
+++ b/python/tests/test.py
@@ -327,9 +327,17 @@ def test_get_reports_from_population(self):
 
         sel_all = self.test_obj['All'].get()
         keys_all = sel_all.ids
-        keys_all_ref = np.asarray([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20])
+        keys_all_ref = np.asarray([10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9])
         self.assertTrue((keys_all == keys_all_ref).all())
 
+        data_begin = sel_all.data[0][:10]
+        data_begin_ref = np.asarray([10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0], dtype=np.float32)
+        self.assertTrue((data_begin == data_begin_ref).all())
+
+        data_end = sel_all.data[-1][-10:]
+        data_end_ref = np.asarray([20.9, 1.9, 2.9, 3.9, 4.9, 5.9, 6.9, 7.9, 8.9, 9.9], dtype=np.float32)
+        self.assertTrue((data_end == data_end_ref).all())
+
         sel_empty = self.test_obj['All'].get(node_ids=[])
         np.testing.assert_allclose(sel_empty.data, np.empty(shape=(0, 0)))
 
@@ -338,6 +346,14 @@ def test_get_node_id_element_id_mapping(self):
         ids_mapping_ref = np.asarray([3, 4])
         self.assertTrue((ids_mapping == ids_mapping_ref).all())
 
+    def test_block_gap_limit(self):
+        ids_mapping = self.test_obj['All'].get_node_id_element_id_mapping([[3, 5]], block_gap_limit=4194304) # >= 1 x GPFS block
+        ids_mapping_ref = np.asarray([3, 4])
+        self.assertTrue((ids_mapping == ids_mapping_ref).all())
+
+        with self.assertRaises(SonataError):
+            self.test_obj['All'].get_node_id_element_id_mapping([[3, 5]], block_gap_limit=4194303) # < 1 x GPFS block
+
 
 class TestElementReportPopulation(unittest.TestCase):
     def setUp(self):
@@ -395,6 +411,14 @@ def test_get_node_id_element_id_mapping(self):
         ids_mapping_ref = np.asarray([[3, 5], [3, 5], [3, 6], [3, 6], [3, 7], [4, 7], [4, 8], [4, 8], [4, 9], [4, 9]])
         self.assertTrue((ids_mapping == ids_mapping_ref).all())
 
+    def test_block_gap_limit(self):
+        ids_mapping = self.test_obj['All'].get_node_id_element_id_mapping([[3, 5]], block_gap_limit=4194304) # >= 1 x GPFS block
+        ids_mapping_ref = np.asarray([[3, 5], [3, 5], [3, 6], [3, 6], [3, 7], [4, 7], [4, 8], [4, 8], [4, 9], [4, 9]])
+        self.assertTrue((ids_mapping == ids_mapping_ref).all())
+
+        with self.assertRaises(SonataError):
+            self.test_obj['All'].get_node_id_element_id_mapping([[3, 5]], block_gap_limit=4194303) # < 1 x GPFS block
+
 
 class TestNodePopulationFailure(unittest.TestCase):
     def test_CorrectStructure(self):
diff --git a/tests/test_report_reader.cpp b/tests/test_report_reader.cpp
index c3f3af9a..7ce102af 100644
--- a/tests/test_report_reader.cpp
+++ b/tests/test_report_reader.cpp
@@ -89,18 +89,23 @@ TEST_CASE("SomaReportReader", "[base]") {
     REQUIRE(data.data == std::vector<float>{3.2f, 4.2f, 3.3f, 4.3f, 3.4f, 4.4f, 3.5f, 4.5f});
 
     auto data_all = pop.get();
-    REQUIRE(data_all.ids == DataFrame<NodeID>::DataType{{1,  2,  3,  4,  5,  6,  7,  8,  9,  10,
-                                                         11, 12, 13, 14, 15, 16, 17, 18, 19, 20}});
+    REQUIRE(data_all.ids == DataFrame<NodeID>::DataType{{10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
+                                                         20, 1,  2,  3,  4,  5,  6,  7,  8,  9}});
     REQUIRE(std::vector<float>(data_all.data.begin(), data_all.data.begin() + 10)  == 
-            std::vector<float>{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f,});
+            std::vector<float>{10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f});
     REQUIRE(std::vector<float>(data_all.data.end() - 10, data_all.data.end()) == 
-            std::vector<float>{11.9f, 12.9f, 13.9f, 14.9f, 15.9f, 16.9f, 17.9f, 18.9f, 19.9f, 20.9f});
+            std::vector<float>{20.9f, 1.9f, 2.9f, 3.9f, 4.9f, 5.9f, 6.9f, 7.9f, 8.9f, 9.9f});
 
     auto data_empty = pop.get(Selection({}));
     REQUIRE(data_empty.data == std::vector<float>{});
 
     auto ids = pop.getNodeIdElementIdMapping(Selection({{3, 5}}));
     REQUIRE(ids == std::vector<NodeID>{3, 4});
+
+    ids = pop.getNodeIdElementIdMapping(Selection({{3, 5}}), 4194304); // >= 1 x GPFS block
+    REQUIRE(ids == std::vector<NodeID>{3, 4});
+
+    REQUIRE_THROWS(pop.getNodeIdElementIdMapping(Selection({{3, 5}}), 4194303)); // < 1 x GPFS block
 }
 
 TEST_CASE("ElementReportReader limits", "[base]") {
@@ -166,4 +171,9 @@ TEST_CASE("ElementReportReader", "[base]") {
 
     auto ids = pop.getNodeIdElementIdMapping(Selection({{3, 5}}));
     REQUIRE(ids == std::vector<CompartmentID>{{3, 5}, {3, 5}, {3, 6}, {3, 6}, {3, 7}, {4, 7}, {4, 8}, {4, 8}, {4, 9}, {4, 9}});
+
+    ids = pop.getNodeIdElementIdMapping(Selection({{3, 5}}), 4194304); // >= 1 x GPFS block
+    REQUIRE(ids == std::vector<CompartmentID>{{3, 5}, {3, 5}, {3, 6}, {3, 6}, {3, 7}, {4, 7}, {4, 8}, {4, 8}, {4, 9}, {4, 9}});
+
+    REQUIRE_THROWS(pop.getNodeIdElementIdMapping(Selection({{3, 5}}), 4194303)); // < 1 x GPFS block
 }