Repository: BlueBrain/libsonata (language: python)

Issue:
Add support for IO-blocks
The improvements provided in #174 and available in `libsonata v0.1.11`, introduced a specific limitation while fetching a small number of widely-disperse GIDs in very large reports. The `{min,max}` optimization, inspired by [ROMIO's Data Sieving](https://web.cels.anl.gov/~thakur/papers/romio-coll.pdf), would in such case force the implementation to retrieve large amounts of data without considering the gap between each of the position ranges of the GIDs. This means that, if we fetch only 2 GIDs that coincidentally are located in the very first entry of the dataset and in the very last one, the library would then read the complete report file just to filter in memory the two selected values.

To overcome these constraints while still trying to reduce the number of IO requests requested to the parallel file system, the current PR introduces several changes and optimizations:

- The `{min,max}` range is now replaced by a collection of ranges or IO blocks, in which each range contains a subset of GIDs. The gap between each IO block is fixed at 64MB (or 4 blocks in GPFS), allowing the code to adapt depending on the number of GIDs requested and their distribution in the file. For instance, a higher number of GIDs would generally imply a lower number of IO blocks, and vice versa.
- A flat-map inspired implementation replaces the existing `std::map` that is utilized to index the GIDs and their positions. The idea is to reduce the overhead when opening populations from very large report files with millions of GIDs. Moreover, it allows us to re-index the content based on different keys, without moving the actual data. For this purpose, the implementation uses several `std::vector` and a dedicated index that refers to the values contained in these vectors. When creating the IO blocks, it is enough to re-index the content by the position ranges and to create the groups based on this new index.
- Other minor changes are included, such as behaviour consistency between asking for all the GIDs or a subset of the GIDs, clarification comments, and others.

Using the same ~60TB report utilized in #174 and #179, the following table illustrates the performance difference between `v0.1.11` and the proposed PR by fetching a single timestep while varying the number of GIDs requested. The `Open` column reflects the cost of opening the population as part of the elapsed execution time, while the `Select` column shows the cost of the operations related to the selection of the GIDs in Python, and the `Get` column illustrates the cost of the actual retrieval of the data. The total cost is reflected in `Elapsed`. All measurements are provided in seconds:

```
v0.1.11
-------------------------------------------------------
   GIDs |   Elapsed |      Open |    Select |       Get
      1 |  5.347993 |  4.720139 |  0.626175 |  0.001667
     10 | 11.422943 |  5.391727 |  0.781538 |  5.249666
    100 | 11.814635 |  5.087994 |  0.763447 |  5.963183
   1000 | 12.090420 |  5.105253 |  0.774254 |  6.210903
  10000 | 12.084519 |  5.047278 |  0.779182 |  6.258047
 100000 | 12.789811 |  5.165905 |  0.785131 |  6.838766
1000000 | 18.066422 |  5.148455 |  0.834916 | 12.083038
4234929 | 33.317796 |  5.097557 |  0.414970 | 27.805258

PR
-------------------------------------------------------
   GIDs |   Elapsed |      Open |    Select |       Get
      1 |  1.558976 |  0.933911 |  0.623337 |  0.001715
     10 |  1.450080 |  0.733224 |  0.628178 |  0.088669
    100 |  4.454307 |  0.749991 |  0.638895 |  3.065413
   1000 |  7.777399 |  0.803403 |  0.632553 |  6.341434
  10000 |  7.836496 |  0.821947 |  0.635333 |  6.379207
 100000 |  8.150157 |  0.815197 |  0.636635 |  6.698314
1000000 | 12.450882 |  0.806348 |  0.706298 | 10.938223
4234929 | 21.047557 |  0.780823 |  0.272458 | 19.994266
```

From the results, it can be observed that the use of multiple-blocks and the integration of a flat-map reduces the overhead of the execution time considerably, specially in small number of GIDs. After 100 GIDs, the implementation can only find a single IO block to fetch from storage and, thus, its behaviour is equivalent to the original `{min,max}` single-block implementation of `v0.1.11`. Hence, allowing us to be more flexible in lower number of GIDs, and more aggressive in large number of GIDs to prevent overheads in the parallel file system. This is clearly noticeable by observing the `Get` column of the table above (note that the performance is improved because the PR also contains the changes in #179).

The following table reflects the number of blocks created by the implementation for each of the tests, showing how we can now be more flexible in requests with lower number of GIDs and more aggressive in requests with large number of GIDs:

```
PR
-----------------
   GIDs |  Blocks
      1 |       1
      2 |       2
     10 |      10
    100 |      32
   1000 |       1
  10000 |       1
 100000 |       1
1000000 |       1
4234929 |       1
```

If we now evaluate the edge case mentioned in the beginning of this description, the entry with 2 GIDs in the following table illustrates the performance by fetching the very first and the very last GID of the same ~60TB report for a single timestep:

```
v0.1.11
-------------------------------------------------------
   GIDs |   Elapsed |      Open |    Select |       Get
      2 | 11.180127 |  4.742555 |  0.205418 |  6.232141

PR
-------------------------------------------------------
   GIDs |   Elapsed |      Open |    Select |       Get
      2 |  1.103376 |  0.893741 |  0.209070 |  0.000555
```


Discussion / hints:
Note that I'll update this PR after #182 is merged. Under my opinion, the logic plan would be to merge first #182, then try to merge #183 after it is ready, and finally refactor the code as proposed in #181. We can also combine #181 with this PR #183 and do both things together, although I'd prefer if we could keep it separated for now.

Nice speedup!

> Under my opinion, the logic plan would be to ...

That makes sense; don't worry about #181, I can just redo that when I have time.  I just merged #182, so we should be good to go.
A quick update. I have adapted the code of the PR after rebasing from `master`, including the use of `std::next` and `std::advance` with iterators instead of raw pointers. I have also replaced the ENV variable for the `block_gap_size` with a parameter, that is also exposed through Python. The parameter includes a boundary check to prevent users affecting the performance on GPFS and there are new CI tests to guarantee that we do not alter this behaviour accidentally in the future. Finally, I have also verified that the performance remains equivalent to what we had with raw pointers and the use of the ENV variable.

I will wait for further comments from your side. Also, thanks to @matz-e, @mgeplf, and @iomaganaris for the comments / insight.
> Awesome, and thanks for adding the python tests, it's appreciated.

Thank you, @mgeplf.

By the way, I forgot to add a small thing in the Python tests, can you approve once again, please?

Fix patch:
diff --git a/src/edge_index.cpp b/src/edge_index.cpp
index 58bbb5c5..c99b5b8e 100644
--- a/src/edge_index.cpp
+++ b/src/edge_index.cpp
@@ -17,6 +17,7 @@
 #include <unordered_map>
 #include <vector>
 
+#include "read_bulk.hpp"
 
 namespace bbp {
 namespace sonata {
@@ -53,7 +54,6 @@ const HighFive::Group targetIndex(const HighFive::Group& h5Root) {
     return h5Root.getGroup(TARGET_INDEX_GROUP);
 }
 
-
 Selection resolve(const HighFive::Group& indexGroup, const NodeID nodeID) {
     if (nodeID >= indexGroup.getDataSet(NODE_ID_TO_RANGES_DSET).getSpace().getDimensions()[0]) {
         // Returning empty set for out-of-range node IDs, to be aligned with SYN2 reader
@@ -90,19 +90,54 @@ Selection resolve(const HighFive::Group& indexGroup, const NodeID nodeID) {
     return Selection(std::move(ranges));
 }
 
-
 Selection resolve(const HighFive::Group& indexGroup, const std::vector<NodeID>& nodeIDs) {
+    constexpr size_t min_gap_size = SONATA_PAGESIZE / (2 * sizeof(uint64_t));
+    constexpr size_t max_aggregated_block_size = 128 * min_gap_size;
+
     if (nodeIDs.size() == 1) {
         return resolve(indexGroup, nodeIDs[0]);
     }
-    // TODO optimize: bulk read for primary index
-    // TODO optimize: range merging
-    std::set<EdgeID> merged;
-    for (NodeID nodeID : nodeIDs) {
-        const auto ids = resolve(indexGroup, nodeID).flatten();
-        merged.insert(ids.begin(), ids.end());
+
+    auto readBlock = [&indexGroup](auto& buffer, const auto& range, const std::string& dset_name) {
+        size_t i_begin = std::get<0>(range);
+        size_t i_end = std::get<1>(range);
+
+        indexGroup.getDataSet(dset_name).select({i_begin, 0}, {i_end - i_begin, 2}).read(buffer);
+    };
+
+    auto primaryRange = bulk_read::bulkRead<std::array<uint64_t, 2>>(
+        [&readBlock](auto& buffer, const auto& range) {
+            readBlock(buffer, range, NODE_ID_TO_RANGES_DSET);
+        },
+        Selection::fromValues(nodeIDs),
+        min_gap_size,
+        max_aggregated_block_size);
+
+    // Sort and eliminate empty ranges.
+    primaryRange = bulk_read::sortAndMerge(primaryRange);
+    if (primaryRange.empty()) {
+        return Selection({});
     }
-    return Selection::fromValues(merged.begin(), merged.end());
+
+    auto secondaryRange = bulk_read::bulkRead<std::array<uint64_t, 2>>(
+        [&readBlock](auto& buffer, const auto& range) {
+            readBlock(buffer, range, RANGE_TO_EDGE_ID_DSET);
+        },
+        primaryRange,
+        min_gap_size,
+        max_aggregated_block_size);
+
+    // Sort and eliminate empty ranges.
+    secondaryRange = bulk_read::sortAndMerge(secondaryRange);
+
+    // Copy `secondaryRange`, because the types don't match.
+    Selection::Ranges edgeIds;
+    edgeIds.reserve(secondaryRange.size());
+    for (const auto& range : secondaryRange) {
+        edgeIds.emplace_back(range[0], range[1]);
+    }
+
+    return Selection(std::move(edgeIds));
 }
 
 
diff --git a/src/population.cpp b/src/population.cpp
index 53d6efb1..ef93ec99 100644
--- a/src/population.cpp
+++ b/src/population.cpp
@@ -17,6 +17,7 @@
 #include <highfive/H5File.hpp>
 
 #include "population.hpp"
+#include "read_bulk.hpp"
 
 namespace bbp {
 namespace sonata {
@@ -34,25 +35,7 @@ void _checkRanges(const Ranges& ranges) {
 }
 
 Ranges _sortAndMerge(const Ranges& ranges) {
-    if (ranges.empty()) {
-        return ranges;
-    }
-    Ranges ret;
-    Ranges sorted(ranges);
-    std::sort(sorted.begin(), sorted.end());
-
-    auto it = sorted.cbegin();
-    ret.push_back(*(it++));
-
-    for (; it != sorted.cend(); ++it) {
-        auto& last = ret[ret.size() - 1].second;
-        if (last < it->first) {
-            ret.push_back(*it);
-        } else {
-            last = std::max(last, it->second);
-        }
-    }
-    return ret;
+    return bulk_read::sortAndMerge(ranges);
 }
 
 Selection intersection_(const Ranges& lhs, const Ranges& rhs) {
@@ -122,11 +105,7 @@ Selection::Values Selection::flatten() const {
 
 
 size_t Selection::flatSize() const {
-    size_t result = 0;
-    for (const auto& range : ranges_) {
-        result += (range.second - range.first);
-    }
-    return result;
+    return bulk_read::detail::flatSize(ranges_);
 }
 
 
diff --git a/src/read_bulk.hpp b/src/read_bulk.hpp
new file mode 100644
index 00000000..adbc88cc
--- /dev/null
+++ b/src/read_bulk.hpp
@@ -0,0 +1,213 @@
+#pragma once
+
+#include <bbp/sonata/population.h>
+#include <cstdint>
+
+#define SONATA_PAGESIZE (4 * 1 << 20)
+
+namespace bbp {
+namespace sonata {
+namespace bulk_read {
+
+namespace detail {
+
+/** Is the selection sorted and non-overlapping?
+ */
+template <class Range>
+bool isCanonical(const std::vector<Range>& ranges) {
+    for (size_t i = 1; i < ranges.size(); ++i) {
+        if (std::get<1>(ranges[i - 1]) > std::get<0>(ranges[i])) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+inline bool isCanonical(const Selection selection) {
+    return isCanonical(selection.ranges());
+}
+
+/** Number of elements in the selection.
+ */
+template <class Range>
+size_t flatSize(const std::vector<Range>& ranges) {
+    size_t size = 0;
+    for (const auto& range : ranges) {
+        size += std::get<1>(range) - std::get<0>(range);
+    }
+
+    return size;
+}
+
+template <class T, class Pred>
+void erase_if(std::vector<T>& v, Pred pred) {
+    auto it = std::remove_if(v.begin(), v.end(), pred);
+    v.erase(it, v.end());
+}
+}  // namespace detail
+
+/** Sort the selection and merge small gaps.
+ *
+ * The ranges of the selection are sorted and then ranges that are separated by
+ * a gap of less than `min_gap_size` are merged into a single range. A gap size
+ * of: `0` will only sort, `1` will merge ranges that touch, pass the page size
+ * to merge anything that's separated less than one page.
+ *
+ * Reading only a few elements from each page could lead to reading very large
+ * ranges into memory. To avoid this issue the one must pick a maximum
+ * aggregated block size. Consecutive blocks wont be merged if the current
+ * block exceeds the threshold `max_aggregated_block_size`.
+ *
+ * Note, that the returned selection is canonical (sorted and non-overlapping)
+ * and may have a larger `flatSize` than `ranges`. Additionally any empty
+ * ranges are removed.
+ */
+template <class Range>
+std::vector<Range> sortAndMerge(const std::vector<Range>& ranges,
+                                size_t min_gap_size = 1,
+                                size_t max_aggregated_block_size = size_t(-1)) {
+    if (ranges.empty()) {
+        return std::vector<Range>{};
+    }
+
+    std::vector<Range> ret;
+    std::vector<Range> sorted(ranges);
+    std::sort(sorted.begin(), sorted.end());
+
+    detail::erase_if(sorted,
+                     [](const auto& range) { return std::get<0>(range) >= std::get<1>(range); });
+
+    if (sorted.empty()) {
+        return std::vector<Range>{};
+    }
+
+    auto it = sorted.cbegin();
+    ret.push_back(*(it++));
+
+    for (; it != sorted.cend(); ++it) {
+        auto& current_range = ret.back();
+        auto& last = std::get<1>(current_range);
+
+        size_t current_range_size = last - std::get<0>(current_range);
+        if (last + min_gap_size <= std::get<0>(*it) ||
+            current_range_size >= max_aggregated_block_size) {
+            // Start a new range.
+            ret.push_back(*it);
+        } else {
+            // Extend the current range.
+            last = std::max(last, std::get<1>(*it));
+        }
+    }
+
+    return ret;
+}
+
+inline Selection sortAndMerge(const Selection& selection, size_t min_gap_size = 0) {
+    return Selection(sortAndMerge(selection.ranges(), min_gap_size));
+}
+
+
+/** Extract a slice of values from a block.
+ *
+ * This is one step of the merge-read-extract algorithm. It operates on one
+ * block of data, i.e. one merged range.
+ *
+ * Remember ranges refer to the indices of the array on disk. They are not the
+ * indices into either `out` or `buffer`.
+ *
+ * This copies the slice `range_small` from `buffer` to `out`. The `buffer` has
+ * size `range_large[1] - range_large[0]`. The `out` pointer must point to the
+ * first element to be written to.
+ */
+template <class T, class Range>
+void extractBlock(T* const out,
+                  T const* const buffer,
+                  const Range& range_large,
+                  const Range& range_small) {
+    size_t i0_small = std::get<0>(range_small);
+    size_t i1_small = std::get<1>(range_small);
+
+    size_t i0_large = std::get<0>(range_large);
+
+    size_t buffer_offset = i0_small - i0_large;
+    size_t count_small = i1_small - i0_small;
+
+    std::copy(buffer + buffer_offset, buffer + buffer_offset + count_small, out);
+}
+
+/** Read larger block and extract required values in memory.
+ *
+ *  This implements the read and extract parts of the merge-read-extract
+ *  algorithm.
+ *
+ *  The ranges denote the indices of the array on disk. The algorithm will read
+ *  one block, i.e. one element of `ranges` at a time. It reads a block by
+ *  calling
+ *
+ *      readBlock(buffer, range);
+ *
+ *  the function object `readBlock` must fill `buffer` (an `std::vector<T>`)
+ *  with the values for `range`. The algorithm will then extract the required
+ *  values (controlled by `subranges`).
+ *
+ *  Note that both `ranges` and `subranges` must be canonical (sorted and
+ *  non-overlapping). Additionally, any range in `subranges` must be fully
+ *  contained in exactly one range in `ranges`.
+ */
+template <class T, class F, class Range>
+std::vector<T> bulkRead(F readBlock,
+                        const std::vector<Range>& ranges,
+                        const std::vector<Range>& subranges) {
+    std::vector<T> values(detail::flatSize(subranges));
+    T* values_ptr = values.data();
+
+    std::vector<T> buffer;
+
+    size_t k_sub = 0;
+    size_t n_sub = subranges.size();
+    for (const auto& range : ranges) {
+        readBlock(buffer, range);
+
+        for (; k_sub < n_sub; ++k_sub) {
+            const auto& subrange = subranges[k_sub];
+            if (std::get<1>(subrange) > std::get<1>(range)) {
+                break;
+            }
+
+            extractBlock(values_ptr, buffer.data(), range, subrange);
+            values_ptr += std::get<1>(subrange) - std::get<0>(subrange);
+        }
+    }
+
+    return values;
+}
+
+/** Read `ranges` using merge-read-extract.
+ *
+ *  @sa `sortAndMerge` and `bulkRead`.
+ */
+template <class T, class F, class Ranges>
+std::vector<T> bulkRead(F readBlock,
+                        const std::vector<Ranges>& ranges,
+                        size_t min_gap_size,
+                        size_t max_aggregated_block_size) {
+    auto super_ranges = sortAndMerge(ranges, min_gap_size, max_aggregated_block_size);
+    return bulkRead<T>(readBlock, super_ranges, ranges);
+}
+
+/** Read `ranges` using merge-read-extract.
+ *
+ *  @sa `sortAndMerge` and `bulkRead`.
+ */
+template <class T, class F>
+std::vector<T> bulkRead(F readBlock,
+                        const Selection& selection,
+                        size_t min_gap_size,
+                        size_t max_aggregated_block_size) {
+    return bulkRead<T>(readBlock, selection.ranges(), min_gap_size, max_aggregated_block_size);
+}
+
+}  // namespace bulk_read
+}  // namespace sonata
+}  // namespace bbp


Test patch:
diff --git a/tests/test_edges.cpp b/tests/test_edges.cpp
index 61592541..48af1a66 100644
--- a/tests/test_edges.cpp
+++ b/tests/test_edges.cpp
@@ -27,6 +27,14 @@ std::ostream& operator<<(std::ostream& oss, const Selection& selection) {
 
 
 TEST_CASE("EdgePopulation", "[edges]") {
+    // E-mapping-good              Dataset {6}
+    // Data:
+    //      2, 1, 2, 0, 2, 2
+    //
+    // @library/E-mapping-good     Dataset {3}
+    // Data:
+    //      "A", "B", "C"
+
     const EdgePopulation population("./data/edges1.h5", "", "edges-AB");
 
     CHECK(population.source() == "nodes-A");
@@ -60,8 +68,21 @@ TEST_CASE("EdgePopulation", "[edges]") {
     // duplicate node IDs are ignored; order of node IDs is not relevant
     CHECK(population.connectingEdges({2, 1, 2}, {2, 1, 2}) == Selection({{0, 4}}));
 
-    CHECK(population.getAttribute<size_t>("E-mapping-good", Selection({{0, 1}, {2, 3}})) ==
-          std::vector<size_t>{2, 2});
+    // Canonical.
+    CHECK(population.getAttribute<size_t>("E-mapping-good", Selection({{0, 1}, {2, 4}})) ==
+          std::vector<size_t>{2, 2, 0});
+
+    // Non-overlapping but not sorted.
+    CHECK(population.getAttribute<size_t>("E-mapping-good", Selection({{2, 4}, {0, 1}})) ==
+          std::vector<size_t>{2, 0, 2});
+
+    // Sorted, but overlapping.
+    CHECK(population.getAttribute<size_t>("E-mapping-good", Selection({{0, 2}, {1, 4}})) ==
+          std::vector<size_t>{2, 1, 1, 2, 0});
+
+    // Overlapping and not sorted.
+    CHECK(population.getAttribute<size_t>("E-mapping-good", Selection({{1, 4}, {0, 2}})) ==
+          std::vector<size_t>{1, 2, 0, 2, 1});
 
     CHECK(population.getAttribute<std::string>("E-mapping-good", Selection({{0, 1}})) ==
           std::vector<std::string>{"C"});