diff --git a/arbor/algorithms.hpp b/arbor/algorithms.hpp
deleted file mode 100644
index bcc3b0f1ca049ce7ff463d1c5d0d7fa73e539166..0000000000000000000000000000000000000000
--- a/arbor/algorithms.hpp
+++ /dev/null
@@ -1,341 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <iostream>
-#include <iterator>
-#include <numeric>
-#include <type_traits>
-#include <vector>
-
-#include <arbor/assert.hpp>
-
-#include <arbor/util/compat.hpp>
-#include <util/meta.hpp>
-#include <util/range.hpp>
-#include <util/rangeutil.hpp>
-
-/*
- * Some simple wrappers around stl algorithms to improve readability of code
- * that uses them.
- *
- * For example, a simple sum() wrapper for finding the sum of a container,
- * is much more readable than using std::accumulate()
- *
- */
-
-namespace arb {
-namespace algorithms {
-
-template <typename C>
-typename util::sequence_traits<C>::value_type
-mean(C const& c)
-{
-    return util::sum(c)/std::size(c);
-}
-
-// returns the prefix sum of c in the form `[0, c[0], c[0]+c[1], ..., sum(c)]`.
-// This means that the returned vector has one more element than c.
-template <typename C>
-C make_index(C const& c)
-{
-    static_assert(
-        std::is_integral<typename C::value_type>::value,
-        "make_index only applies to integral types"
-    );
-
-    C out(c.size()+1);
-    out[0] = 0;
-    std::partial_sum(c.begin(), c.end(), out.begin()+1);
-    return out;
-}
-
-/// test for membership within half-open interval
-template <typename X, typename I, typename J>
-bool in_interval(const X& x, const I& lower, const J& upper) {
-    return x>=lower && x<upper;
-}
-
-template <typename X, typename I, typename J>
-bool in_interval(const X& x, const std::pair<I, J>& bounds) {
-    return x>=bounds.first && x<bounds.second;
-}
-
-/// works like std::is_sorted(), but with stronger condition that succesive
-/// elements must be greater than those before them
-template <typename C>
-bool is_strictly_monotonic_increasing(C const& c)
-{
-    using value_type = typename C::value_type;
-    return std::is_sorted(
-        c.begin(),
-        c.end(),
-        [] (value_type const& lhs, value_type const& rhs) {
-            return lhs <= rhs;
-        }
-    );
-}
-
-template <typename C>
-bool is_strictly_monotonic_decreasing(C const& c)
-{
-    using value_type = typename C::value_type;
-    return std::is_sorted(
-        c.begin(),
-        c.end(),
-        [] (value_type const& lhs, value_type const& rhs) {
-            return lhs >= rhs;
-        }
-    );
-}
-
-// Check if c[0] == 0 and c[i] < 0 holds for i != 0
-// Also handle the valid case of c[0]==value_type(-1)
-// This means that children of a node always have larger indices than their
-// parent.
-template <
-    typename C,
-    typename = typename std::enable_if<std::is_integral<typename C::value_type>::value>
->
-bool is_minimal_degree(C const& c)
-{
-    static_assert(
-        std::is_integral<typename C::value_type>::value,
-        "is_minimal_degree only applies to integral types"
-    );
-
-    if (c.size()==0u) {
-        return true;
-    }
-
-    using value_type = typename C::value_type;
-    if (!(c[0]==value_type(0) || c[0]==value_type(-1))) {
-        return false;
-    }
-    auto i = value_type(1);
-    auto it = std::find_if(
-        c.begin()+1, c.end(), [&i](value_type v) { return v>=(i++); }
-    );
-    return it==c.end();
-}
-
-template <typename C>
-bool all_positive(const C& c) {
-    return util::all_of(c, [](auto v) { return v>decltype(v){}; });
-}
-
-template <typename C>
-bool all_negative(const C& c) {
-    return util::all_of(c, [](auto v) { return v<decltype(v){}; });
-}
-
-// returns a vector containing the number of children for each node.
-template<typename C>
-std::vector<typename C::value_type> child_count(const C& parent_index)
-{
-    using value_type = typename C::value_type;
-    static_assert(
-        std::is_integral<value_type>::value,
-        "integral type required"
-    );
-
-    std::vector<value_type> count(parent_index.size(), 0);
-    for (auto i = 0u; i < parent_index.size(); ++i) {
-        auto p = parent_index[i];
-        // -1 means no parent
-        if (p != value_type(i) && p != value_type(-1)) {
-            ++count[p];
-        }
-    }
-
-    return count;
-}
-
-template<typename C>
-bool has_contiguous_compartments(const C& parent_index)
-{
-    using value_type = typename C::value_type;
-    static_assert(
-        std::is_integral<value_type>::value,
-        "integral type required"
-    );
-
-    if (!is_minimal_degree(parent_index)) {
-        return false;
-    }
-
-    auto num_child = child_count(parent_index);
-    for (auto i=1u; i < parent_index.size(); ++i) {
-        auto p = parent_index[i];
-        if (num_child[p]==1 && p!=value_type(i-1)) {
-            return false;
-        }
-    }
-
-    return true;
-}
-
-template<typename C>
-std::vector<typename C::value_type> branches(const C& parent_index)
-{
-    static_assert(
-        std::is_integral<typename C::value_type>::value,
-        "integral type required"
-    );
-
-    arb_assert(has_contiguous_compartments(parent_index));
-
-    std::vector<typename C::value_type> branch_index;
-    if (parent_index.empty()) {
-        return branch_index;
-    }
-
-    auto num_child = child_count(parent_index);
-    branch_index.push_back(0);
-    for (std::size_t i = 1; i < parent_index.size(); ++i) {
-        auto p = parent_index[i];
-        if (num_child[p] > 1 || parent_index[p] == p) {
-            // `parent_index[p] == p` ~> parent_index[i] is the soma
-            branch_index.push_back(i);
-        }
-    }
-
-    branch_index.push_back(parent_index.size());
-    return branch_index;
-}
-
-// creates a vector that contains the branch index for each compartment.
-// e.g. {0, 1, 5, 9, 10} -> {0, 1, 1, 1, 1, 2, 2, 2, 2, 3}
-//                  indices  0  1           5           9
-template<typename C>
-std::vector<typename C::value_type> expand_branches(const C& branch_index)
-{
-    static_assert(
-        std::is_integral<typename C::value_type>::value,
-        "integral type required"
-    );
-
-    if (branch_index.empty())
-        return {};
-
-    std::vector<typename C::value_type> expanded(branch_index.back());
-    for (std::size_t i = 0; i < branch_index.size()-1; ++i) {
-        for (auto j = branch_index[i]; j < branch_index[i+1]; ++j) {
-            expanded[j] = i;
-        }
-    }
-
-    return expanded;
-}
-
-template<typename C>
-typename C::value_type find_branch(const C& branch_index,
-                                   typename C::value_type nid)
-{
-    using value_type = typename C::value_type;
-    static_assert(
-        std::is_integral<value_type>::value,
-        "integral type required"
-    );
-
-    auto it =  std::find_if(
-        branch_index.begin(), branch_index.end(),
-        [nid](const value_type& v) { return v > nid; }
-    );
-
-    return it - branch_index.begin() - 1;
-}
-
-// Find the reduced form of a tree represented as a parent index.
-// The operation of transforming a tree into a homeomorphic pre-image is called
-// 'reduction'. The 'reduced' tree is the smallest tree that maps into the tree
-// homeomorphically (i.e. preserving labels and child relationships).
-// For example, the tree represented by the following index:
-//      {0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12}
-// Has reduced form
-//      {0, 0, 0, 0, 3, 3}
-//
-// This transformation can be represented graphically:
-//
-//        0               0
-//       /|\             /|\.
-//      1 4 6      =>   1 2 3
-//     /  |  \             / \.
-//    2   5   7           4   5
-//   /         \.
-//  3           8
-//             / \.
-//            9   11
-//           /     \.
-//          10     12
-//                   \.
-//                   13
-//
-template<typename C>
-std::vector<typename C::value_type> tree_reduce(
-    const C& parent_index, const C& branch_index)
-{
-    static_assert(
-        std::is_integral<typename C::value_type>::value,
-        "integral type required"
-    );
-
-    if (parent_index.empty() && branch_index.empty()) {
-        return {};
-    }
-
-    arb_assert(parent_index.size()-branch_index.back() == 0);
-    arb_assert(has_contiguous_compartments(parent_index));
-    arb_assert(is_strictly_monotonic_increasing(branch_index));
-
-    // expand the branch index to lookup the banch id for each compartment
-    auto expanded_branch = expand_branches(branch_index);
-
-    std::vector<typename C::value_type> new_parent_index;
-    // push the first element manually as the parent of the root might be -1
-    new_parent_index.push_back(expanded_branch[0]);
-    for (std::size_t i = 1; i < branch_index.size()-1; ++i) {
-        auto p = parent_index[branch_index[i]];
-        new_parent_index.push_back(expanded_branch[p]);
-    }
-
-    return new_parent_index;
-}
-
-template <typename Seq, typename = util::enable_if_sequence_t<Seq&>>
-bool is_unique(const Seq& seq) {
-    using std::begin;
-    using std::end;
-
-    return std::adjacent_find(begin(seq), end(seq)) == end(seq);
-}
-
-
-/// Binary search, because std::binary_search doesn't return information
-/// about where a match was found.
-
-// TODO: consolidate these with rangeutil routines; make them sentinel
-
-template <typename It, typename T>
-It binary_find(It b, It e, const T& value) {
-    auto it = std::lower_bound(b, e, value);
-    return it==e ? e : (*it==value ? it : e);
-}
-
-template <typename Seq, typename T>
-auto binary_find(const Seq& seq, const T& value) {
-    using std::begin;
-    using std::end;
-
-    return binary_find(begin(seq), end(seq), value);
-}
-
-template <typename Seq, typename T>
-auto binary_find(Seq& seq, const T& value) {
-    using std::begin;
-    using std::end;
-
-    return binary_find(begin(seq), end(seq), value);
-}
-
-} // namespace algorithms
-} // namespace arb
diff --git a/arbor/backends/gpu/forest.cpp b/arbor/backends/gpu/forest.cpp
index e920147cb0fb117fe5a812a7ab5aff821fae3968..cc146af02370a6da1ca4ef6b8e35140549bfe233 100644
--- a/arbor/backends/gpu/forest.cpp
+++ b/arbor/backends/gpu/forest.cpp
@@ -1,5 +1,4 @@
 #include "backends/gpu/forest.hpp"
-#include "tree.hpp"
 #include "util/span.hpp"
 
 namespace arb {
@@ -31,7 +30,7 @@ forest::forest(const std::vector<size_type>& p, const std::vector<size_type>& ce
         }
 
         // find the index of the first node for each branch
-        auto branch_starts = algorithms::branches(fine_tree.parents());
+        auto branch_starts = branches(fine_tree.parents());
 
         // compute branch length and apply permutation
         std::vector<unsigned> branch_lengths(branch_starts.size() - 1);
@@ -43,7 +42,7 @@ forest::forest(const std::vector<size_type>& p, const std::vector<size_type>& ce
         // we need to convert to cell_lid_type, required to construct a tree.
         std::vector<cell_lid_type> branch_p =
             util::assign_from(
-                algorithms::tree_reduce(fine_tree.parents(), branch_starts));
+                tree_reduce(fine_tree.parents(), branch_starts));
         // build tree structure that describes the branch topology
         auto cell_tree = tree(branch_p);
 
diff --git a/arbor/backends/gpu/forest.hpp b/arbor/backends/gpu/forest.hpp
index 439008585dd0fa787736b2dcc19645abbb16ad9c..1bc13bc858151a904b7b2c65ca5b1b55ff0343fe 100644
--- a/arbor/backends/gpu/forest.hpp
+++ b/arbor/backends/gpu/forest.hpp
@@ -135,6 +135,153 @@ private:
     unsigned only_on_level;
 };
 
+// Works like std::is_sorted(), but with stronger condition that successive
+// elements must be greater than those before them
+template <typename C>
+bool is_strictly_monotonic_increasing(C const& c)
+{
+    using value_type = typename C::value_type;
+    return std::is_sorted(
+        c.begin(),
+        c.end(),
+        [] (value_type const& lhs, value_type const& rhs) {
+          return lhs <= rhs;
+        }
+    );
+}
+
+template<typename C>
+bool has_contiguous_compartments(const C& parent_index)
+{
+    using value_type = typename C::value_type;
+    static_assert(
+        std::is_integral<value_type>::value,
+        "integral type required"
+    );
+
+    if (!is_minimal_degree(parent_index)) {
+        return false;
+    }
+
+    auto num_child = child_count(parent_index);
+    for (auto i=1u; i < parent_index.size(); ++i) {
+        auto p = parent_index[i];
+        if (num_child[p]==1 && p!=value_type(i-1)) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+template<typename C>
+std::vector<typename C::value_type> branches(const C& parent_index)
+{
+    static_assert(
+        std::is_integral<typename C::value_type>::value,
+        "integral type required"
+    );
+
+    arb_assert(has_contiguous_compartments(parent_index));
+
+    std::vector<typename C::value_type> branch_index;
+    if (parent_index.empty()) {
+        return branch_index;
+    }
+
+    auto num_child = child_count(parent_index);
+    branch_index.push_back(0);
+    for (std::size_t i = 1; i < parent_index.size(); ++i) {
+        auto p = parent_index[i];
+        if (num_child[p] > 1 || parent_index[p] == p) {
+            // `parent_index[p] == p` ~> parent_index[i] is the soma
+            branch_index.push_back(i);
+        }
+    }
+
+    branch_index.push_back(parent_index.size());
+    return branch_index;
+}
+
+// creates a vector that contains the branch index for each compartment.
+// e.g. {0, 1, 5, 9, 10} -> {0, 1, 1, 1, 1, 2, 2, 2, 2, 3}
+//                  indices  0  1           5           9
+template<typename C>
+std::vector<typename C::value_type> expand_branches(const C& branch_index)
+{
+    static_assert(
+        std::is_integral<typename C::value_type>::value,
+        "integral type required"
+    );
+
+    if (branch_index.empty())
+        return {};
+
+    std::vector<typename C::value_type> expanded(branch_index.back());
+    for (std::size_t i = 0; i < branch_index.size()-1; ++i) {
+        for (auto j = branch_index[i]; j < branch_index[i+1]; ++j) {
+            expanded[j] = i;
+        }
+    }
+
+    return expanded;
+}
+
+// Find the reduced form of a tree represented as a parent index.
+// The operation of transforming a tree into a homeomorphic pre-image is called
+// 'reduction'. The 'reduced' tree is the smallest tree that maps into the tree
+// homeomorphically (i.e. preserving labels and child relationships).
+// For example, the tree represented by the following index:
+//      {0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12}
+// Has reduced form
+//      {0, 0, 0, 0, 3, 3}
+//
+// This transformation can be represented graphically:
+//
+//        0               0
+//       /|\             /|\.
+//      1 4 6      =>   1 2 3
+//     /  |  \             / \.
+//    2   5   7           4   5
+//   /         \.
+//  3           8
+//             / \.
+//            9   11
+//           /     \.
+//          10     12
+//                   \.
+//                   13
+//
+template<typename C>
+std::vector<typename C::value_type> tree_reduce(
+    const C& parent_index, const C& branch_index)
+{
+    static_assert(
+        std::is_integral<typename C::value_type>::value,
+        "integral type required"
+    );
+
+    if (parent_index.empty() && branch_index.empty()) {
+        return {};
+    }
+
+    arb_assert(parent_index.size()-branch_index.back() == 0);
+    arb_assert(has_contiguous_compartments(parent_index));
+    arb_assert(is_strictly_monotonic_increasing(branch_index));
+
+    // expand the branch index to lookup the banch id for each compartment
+    auto expanded_branch = expand_branches(branch_index);
+
+    std::vector<typename C::value_type> new_parent_index;
+    // push the first element manually as the parent of the root might be -1
+    new_parent_index.push_back(expanded_branch[0]);
+    for (std::size_t i = 1; i < branch_index.size()-1; ++i) {
+        auto p = parent_index[branch_index[i]];
+        new_parent_index.push_back(expanded_branch[p]);
+    }
+
+    return new_parent_index;
+}
 
 } // namespace gpu
 } // namespace arb
diff --git a/arbor/backends/gpu/matrix_state_fine.hpp b/arbor/backends/gpu/matrix_state_fine.hpp
index 4d0d46c50dbbdfd8dd7feb141a9462d51ef97274..d460597cfbaa86052d2fa87f2939bebb19467692 100644
--- a/arbor/backends/gpu/matrix_state_fine.hpp
+++ b/arbor/backends/gpu/matrix_state_fine.hpp
@@ -7,7 +7,6 @@
 
 #include <arbor/common_types.hpp>
 
-#include "algorithms.hpp"
 #include "memory/memory.hpp"
 #include "util/partition.hpp"
 #include "util/rangeutil.hpp"
diff --git a/arbor/backends/multicore/multi_event_stream.hpp b/arbor/backends/multicore/multi_event_stream.hpp
index 8b4c283e5a8f4b33ff5a857e608f52ed90e0cbba..efc4d482e196c74b6a06a1dd6c13765a2c5a1afa 100644
--- a/arbor/backends/multicore/multi_event_stream.hpp
+++ b/arbor/backends/multicore/multi_event_stream.hpp
@@ -13,7 +13,6 @@
 
 #include "backends/event.hpp"
 #include "backends/multi_event_stream_state.hpp"
-#include "algorithms.hpp"
 #include "util/range.hpp"
 #include "util/rangeutil.hpp"
 #include "util/strprintf.hpp"
diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp
index 638ad093f4f6dc1b25597e6333ee2c1e9560dca6..032f09c43416a3f15d8f170175009fba4f4ccf34 100644
--- a/arbor/communication/communicator.cpp
+++ b/arbor/communication/communicator.cpp
@@ -8,7 +8,6 @@
 #include <arbor/spike.hpp>
 #include <include/arbor/arbexcept.hpp>
 
-#include "algorithms.hpp"
 #include "communication/gathered_vector.hpp"
 #include "connection.hpp"
 #include "distributed_context.hpp"
@@ -104,7 +103,7 @@ communicator::communicator(const recipe& rec,
     // The loop above gave the information required to construct in place
     // the connections as partitioned by the domain of their source gid.
     connections_.resize(n_cons);
-    connection_part_ = algorithms::make_index(src_counts);
+    util::make_partition(connection_part_, src_counts);
     auto offsets = connection_part_;
     std::size_t pos = 0;
     for (const auto& cell: gid_infos) {
diff --git a/arbor/communication/mpi.hpp b/arbor/communication/mpi.hpp
index 1f23b24afa1dc7d8625c340400a21bd8d7e65bfe..f07412af6c1a548fb7b7ae7427963f1f913041b5 100644
--- a/arbor/communication/mpi.hpp
+++ b/arbor/communication/mpi.hpp
@@ -10,10 +10,9 @@
 #include <arbor/assert.hpp>
 #include <arbor/communication/mpi_error.hpp>
 
-#include "algorithms.hpp"
 #include "communication/gathered_vector.hpp"
 #include "profile/profiler_macro.hpp"
-
+#include "util/partition.hpp"
 
 namespace arb {
 namespace mpi {
@@ -100,8 +99,9 @@ std::vector<T> gather_all(T value, MPI_Comm comm) {
 inline std::vector<std::string> gather(std::string str, int root, MPI_Comm comm) {
     using traits = mpi_traits<char>;
 
-    auto counts = gather_all(int(str.size()), comm);
-    auto displs = algorithms::make_index(counts);
+    std::vector<int> counts, displs;
+    counts = gather_all(int(str.size()), comm);
+    util::make_partition(displs, counts);
 
     std::vector<char> buffer(displs.back());
 
@@ -127,11 +127,12 @@ template <typename T>
 std::vector<T> gather_all(const std::vector<T>& values, MPI_Comm comm) {
 
     using traits = mpi_traits<T>;
-    auto counts = gather_all(int(values.size()), comm);
+    std::vector<int> counts, displs;
+    counts = gather_all(int(values.size()), comm);
     for (auto& c : counts) {
         c *= traits::count();
     }
-    auto displs = algorithms::make_index(counts);
+    util::make_partition(displs, counts);
 
     std::vector<T> buffer(displs.back()/traits::count());
     MPI_OR_THROW(MPI_Allgatherv,
@@ -154,11 +155,12 @@ gathered_vector<T> gather_all_with_partition(const std::vector<T>& values, MPI_C
     // We have to use int for the count and displs vectors instead
     // of count_type because these are used as arguments to MPI_Allgatherv
     // which expects int arguments.
-    auto counts = gather_all(int(values.size()), comm);
+    std::vector<int> counts, displs;
+    counts = gather_all(int(values.size()), comm);
     for (auto& c : counts) {
         c *= traits::count();
     }
-    auto displs = algorithms::make_index(counts);
+    util::make_partition(displs, counts);
 
     std::vector<T> buffer(displs.back()/traits::count());
 
diff --git a/arbor/profile/meter_manager.cpp b/arbor/profile/meter_manager.cpp
index 6b17ac44bc647c4877953e393452fa5579074faa..dc2e12804cd8fa44078dfdec04422b052164e278 100644
--- a/arbor/profile/meter_manager.cpp
+++ b/arbor/profile/meter_manager.cpp
@@ -6,7 +6,6 @@
 #include "memory_meter.hpp"
 #include "power_meter.hpp"
 
-#include "algorithms.hpp"
 #include "execution_context.hpp"
 #include "util/hostname.hpp"
 #include "util/strprintf.hpp"
@@ -18,6 +17,11 @@ namespace profile {
 using timer_type = timer<>;
 using util::strprintf;
 
+template <typename C>
+double mean(const C& c) {
+    return util::sum(c)/(double)std::size(c);
+}
+
 measurement::measurement(std::string n, std::string u,
                          const std::vector<double>& readings,
                          const context& ctx):
@@ -66,7 +70,6 @@ void meter_manager::start(const context& ctx) {
     start_time_ = timer_type::tic();
 };
 
-
 void meter_manager::checkpoint(std::string name, const context& ctx) {
     arb_assert(started_);
 
@@ -155,13 +158,13 @@ std::ostream& operator<<(std::ostream& o, const meter_report& report) {
         for (const auto& m: report.meters) {
             if (m.name=="time") {
                 // Calculate the average time per rank in s.
-                double time = algorithms::mean(m.measurements[cp_index]);
+                double time = mean(m.measurements[cp_index]);
                 sums[m_index] += time;
                 o << strprintf("%16.3f", time);
             }
             else if (m.name.find("memory")!=std::string::npos) {
                 // Calculate the average memory per rank in MB.
-                double mem = algorithms::mean(m.measurements[cp_index])*1e-6;
+                double mem = mean(m.measurements[cp_index])*1e-6;
                 sums[m_index] += mem;
                 o << strprintf("%16.3f", mem);
             }
@@ -176,7 +179,7 @@ std::ostream& operator<<(std::ostream& o, const meter_report& report) {
                 o << strprintf("%16.3f", energy);
             }
             else {
-                double value = algorithms::mean(m.measurements[cp_index]);
+                double value = mean(m.measurements[cp_index]);
                 sums[m_index] += value;
                 o << strprintf("%16.3f", value);
             }
diff --git a/arbor/tree.cpp b/arbor/tree.cpp
index a4bc28d62a2999211dad0423e0cf984716644572..a46718699176774fa3ee5af1a921928d8e72208b 100644
--- a/arbor/tree.cpp
+++ b/arbor/tree.cpp
@@ -6,16 +6,16 @@
 
 #include <arbor/common_types.hpp>
 
-#include "algorithms.hpp"
 #include "memory/memory.hpp"
 #include "tree.hpp"
+#include "util/partition.hpp"
 #include "util/span.hpp"
 
 namespace arb {
 
 tree::tree(std::vector<tree::int_type> parent_index) {
     // validate the input
-    if(!algorithms::is_minimal_degree(parent_index)) {
+    if(!is_minimal_degree(parent_index)) {
         throw std::domain_error(
             "parent index used to build a tree did not satisfy minimal degree ordering"
         );
@@ -29,7 +29,7 @@ tree::tree(std::vector<tree::int_type> parent_index) {
     parents_[0] = no_parent;
 
     // compute offsets into children_ array
-    memory::copy(algorithms::make_index(algorithms::child_count(parents_)), child_index_);
+    arb::util::make_partition(child_index_, child_count(parents_));
 
     std::vector<int_type> pos(parents_.size(), 0);
     for (auto i = 1u; i < parents_.size(); ++i) {
@@ -181,7 +181,7 @@ tree::iarray tree::select_new_root(int_type root) {
     iarray branch_ix (num_nodes, 0);
     // we cannot use the existing `children_` array as we only updated the
     // parent structure yet
-    auto new_num_children = algorithms::child_count(parents_);
+    auto new_num_children = child_count(parents_);
     for (auto n: make_span(num_nodes)) {
         branch_ix[n] = n;
         auto prev = n;
@@ -241,7 +241,7 @@ tree::iarray tree::select_new_root(int_type root) {
 
     // recompute the children array
     memory::copy(new_parents, parents_);
-    memory::copy(algorithms::make_index(algorithms::child_count(parents_)), child_index_);
+    arb::util::make_partition(child_index_, child_count(parents_));
 
     std::vector<int_type> pos(parents_.size(), 0);
     for (auto i = 1u; i < parents_.size(); ++i) {
diff --git a/arbor/tree.hpp b/arbor/tree.hpp
index da55f2a94b9f2dd4d6e1a384d6d74c7fcd30a489..23f44fbdab084166b5dd3a7ef7158a09ff1cddd3 100644
--- a/arbor/tree.hpp
+++ b/arbor/tree.hpp
@@ -8,8 +8,8 @@
 
 #include <arbor/common_types.hpp>
 
-#include "algorithms.hpp"
 #include "memory/memory.hpp"
+#include "util/rangeutil.hpp"
 #include "util/span.hpp"
 
 namespace arb {
@@ -115,49 +115,55 @@ private:
 // The root has depth 0, it's children have depth 1, and so on.
 tree::iarray depth_from_root(const tree& t);
 
-template <typename C>
-std::vector<tree::int_type> make_parent_index(tree const& t, C const& counts)
+// Check if c[0] == 0 and c[i] < 0 holds for i != 0
+// Also handle the valid case of c[0]==value_type(-1)
+// This means that children of a node always have larger indices than their
+// parent.
+template <
+    typename C,
+    typename = typename std::enable_if<std::is_integral<typename C::value_type>::value>
+>
+bool is_minimal_degree(C const& c)
 {
-    using util::make_span;
-    using int_type = tree::int_type;
-    constexpr auto no_parent = tree::no_parent;
-
-    if (!algorithms::all_positive(counts) || counts.size() != t.num_segments()) {
-        throw std::domain_error(
-            "make_parent_index requires one non-zero count per segment"
-        );
-    }
-    auto index = algorithms::make_index(counts);
-    auto num_compartments = index.back();
-    std::vector<int_type> parent_index(num_compartments);
-    int_type pos = 0;
-    for (int_type i : make_span(0, t.num_segments())) {
-        // get the parent of this segment
-        // taking care for the case where the root node has -1 as its parent
-        auto parent = t.parent(i);
-        parent = parent!=no_parent ? parent : 0;
-
-        // the index of the first compartment in the segment
-        // is calculated differently for the root (i.e when i==parent)
-        if (i!=parent) {
-            parent_index[pos++] = index[parent+1]-1;
-        }
-        else {
-            parent_index[pos++] = parent;
-        }
-        // number the remaining compartments in the segment consecutively
-        while (pos<index[i+1]) {
-            parent_index[pos] = pos-1;
-            pos++;
-        }
+    static_assert(
+        std::is_integral<typename C::value_type>::value,
+        "is_minimal_degree only applies to integral types"
+    );
+
+    if (c.size()==0u) {
+        return true;
     }
 
-    // if one of these assertions is tripped, we have to improve
-    // the input validation above
-    assert(pos==num_compartments);
-    assert(algorithms::is_minimal_degree(parent_index));
+    using value_type = typename C::value_type;
+    if (!(c[0]==value_type(0) || c[0]==value_type(-1))) {
+        return false;
+    }
+    auto i = value_type(1);
+    auto it = std::find_if(
+        c.begin()+1, c.end(), [&i](value_type v) { return v>=(i++); }
+    );
+    return it==c.end();
+}
 
-    return parent_index;
+// returns a vector containing the number of children for each node.
+template<typename C>
+std::vector<typename C::value_type> child_count(const C& parent_index)
+{
+    using value_type = typename C::value_type;
+    static_assert(
+        std::is_integral<value_type>::value,
+        "integral type required"
+    );
+
+    std::vector<value_type> count(parent_index.size(), 0);
+    for (auto i = 0u; i < parent_index.size(); ++i) {
+        auto p = parent_index[i];
+        // -1 means no parent
+        if (p != value_type(i) && p != value_type(-1)) {
+            ++count[p];
+        }
+    }
+    return count;
 }
 
 } // namespace arb
diff --git a/arbor/util/maputil.hpp b/arbor/util/maputil.hpp
index 693d91454b55104d3b242049a44bc24cbee88f37..f6b4a939d7feedfa62f935e36b67f6fd66f667ab 100644
--- a/arbor/util/maputil.hpp
+++ b/arbor/util/maputil.hpp
@@ -143,7 +143,6 @@ auto ptr_by_key(C&& c, const Key& k) {
 
 // Find the index into an ordered sequence of a value by binary search;
 // returns optional<size_type> for the size_type associated with the sequence.
-// (Note: this is pretty much all we use algorthim::binary_find for.) 
 
 template <typename C, typename Key>
 std::optional<typename sequence_traits<C>::difference_type> binary_search_index(const C& c, const Key& key) {
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 81de2e85a8c13e78144ee8f68feac185ce965f62..dde25eb7793662268e62666d627523ba0c7e86bf 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -85,7 +85,6 @@ endforeach()
 
 set(unit_sources
     ../common_cells.cpp
-    test_algorithms.cpp
     test_any_cast.cpp
     test_any_ptr.cpp
     test_any_visitor.cpp
@@ -105,8 +104,10 @@ set(unit_sources
     test_event_queue.cpp
     test_expected.cpp
     test_filter.cpp
+    test_forest.cpp
     test_fvm_layout.cpp
     test_fvm_lowered.cpp
+    test_index.cpp
     test_kinetic_linear.cpp
     test_lexcmp.cpp
     test_lif_cell_group.cpp
diff --git a/test/unit/test_algorithms.cpp b/test/unit/test_algorithms.cpp
deleted file mode 100644
index 303996dfb1f8e3f291d70061b1211e64cbdc714f..0000000000000000000000000000000000000000
--- a/test/unit/test_algorithms.cpp
+++ /dev/null
@@ -1,771 +0,0 @@
-#include <forward_list>
-#include <iterator>
-#include <random>
-#include <string>
-#include <vector>
-
-#include "../gtest.h"
-
-#include "algorithms.hpp"
-#include "util/index_into.hpp"
-#include "util/meta.hpp"
-#include "util/rangeutil.hpp"
-
-// (Pending abstraction of threading interface)
-#include <arbor/version.hpp>
-#include "threading/threading.hpp"
-#include "common.hpp"
-
-TEST(algorithms, make_index)
-{
-    {
-        std::vector<int> v(10, 1);
-        auto index = arb::algorithms::make_index(v);
-
-        EXPECT_EQ(index.size(), 11u);
-        EXPECT_EQ(index.front(), 0);
-        EXPECT_EQ(index.back(), arb::util::sum(v));
-    }
-
-    {
-        std::vector<int> v(10);
-        std::iota(v.begin(), v.end(), 1);
-        auto index = arb::algorithms::make_index(v);
-
-        EXPECT_EQ(index.size(), 11u);
-        EXPECT_EQ(index.front(), 0);
-        EXPECT_EQ(index.back(), arb::util::sum(v));
-    }
-}
-
-TEST(algorithms, minimal_degree)
-{
-    {
-        std::vector<int> v = {0};
-        EXPECT_TRUE(arb::algorithms::is_minimal_degree(v));
-    }
-
-    {
-        std::vector<int> v = {0, 0, 1, 2, 3, 4};
-        EXPECT_TRUE(arb::algorithms::is_minimal_degree(v));
-    }
-
-    {
-        std::vector<int> v = {0, 0, 1, 2, 0, 4};
-        EXPECT_TRUE(arb::algorithms::is_minimal_degree(v));
-    }
-
-    {
-        std::vector<int> v = {0, 0, 1, 2, 0, 4, 5, 4};
-        EXPECT_TRUE(arb::algorithms::is_minimal_degree(v));
-    }
-
-    {
-        std::vector<int> v = {1};
-        EXPECT_FALSE(arb::algorithms::is_minimal_degree(v));
-    }
-
-    {
-        std::vector<int> v = {0, 2};
-        EXPECT_FALSE(arb::algorithms::is_minimal_degree(v));
-    }
-
-    {
-        std::vector<int> v = {0, 1, 2};
-        EXPECT_FALSE(arb::algorithms::is_minimal_degree(v));
-    }
-}
-
-TEST(algorithms, is_strictly_monotonic_increasing)
-{
-    EXPECT_TRUE(
-        arb::algorithms::is_strictly_monotonic_increasing(
-            std::vector<int>{0}
-        )
-    );
-    EXPECT_TRUE(
-        arb::algorithms::is_strictly_monotonic_increasing(
-            std::vector<int>{0, 1, 2, 3}
-        )
-    );
-    EXPECT_TRUE(
-        arb::algorithms::is_strictly_monotonic_increasing(
-            std::vector<int>{8, 20, 42, 89}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_strictly_monotonic_increasing(
-            std::vector<int>{0, 0}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_strictly_monotonic_increasing(
-            std::vector<int>{8, 20, 20, 89}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_strictly_monotonic_increasing(
-            std::vector<int>{3, 2, 1, 0}
-        )
-    );
-}
-
-TEST(algorithms, is_strictly_monotonic_decreasing)
-{
-    EXPECT_TRUE(
-        arb::algorithms::is_strictly_monotonic_decreasing(
-            std::vector<int>{0}
-        )
-    );
-    EXPECT_TRUE(
-        arb::algorithms::is_strictly_monotonic_decreasing(
-            std::vector<int>{3, 2, 1, 0}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_strictly_monotonic_decreasing(
-            std::vector<int>{0, 1, 2, 3}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_strictly_monotonic_decreasing(
-            std::vector<int>{8, 20, 42, 89}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_strictly_monotonic_decreasing(
-            std::vector<int>{0, 0}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_strictly_monotonic_decreasing(
-            std::vector<int>{8, 20, 20, 89}
-        )
-    );
-}
-
-TEST(algorithms, all_positive) {
-    using arb::algorithms::all_positive;
-
-    EXPECT_TRUE(all_positive(std::vector<int>{}));
-    EXPECT_TRUE(all_positive(std::vector<int>{3, 2, 1}));
-    EXPECT_FALSE(all_positive(std::vector<int>{3, 2, 1, 0}));
-    EXPECT_FALSE(all_positive(std::vector<int>{-1}));
-
-    EXPECT_TRUE(all_positive((double []){1., 2.}));
-    EXPECT_FALSE(all_positive((double []){1., 0.}));
-    EXPECT_FALSE(all_positive((double []){NAN}));
-
-    EXPECT_TRUE(all_positive((std::string []){"a", "b"}));
-    EXPECT_FALSE(all_positive((std::string []){"a", "", "b"}));
-}
-
-TEST(algorithms, all_negative) {
-    using arb::algorithms::all_negative;
-
-    EXPECT_TRUE(all_negative(std::vector<int>{}));
-    EXPECT_TRUE(all_negative(std::vector<int>{-3, -2, -1}));
-    EXPECT_FALSE(all_negative(std::vector<int>{-3, -2, -1, 0}));
-    EXPECT_FALSE(all_negative(std::vector<int>{1}));
-
-    double negzero = std::copysign(0., -1.);
-
-    EXPECT_TRUE(all_negative((double []){-1., -2.}));
-    EXPECT_FALSE(all_negative((double []){-1., 0.}));
-    EXPECT_FALSE(all_negative((double []){-1., negzero}));
-    EXPECT_FALSE(all_negative((double []){NAN}));
-
-    EXPECT_FALSE(all_negative((std::string []){"", "b"}));
-    EXPECT_FALSE(all_negative((std::string []){""}));
-}
-
-TEST(algorithms, has_contiguous_compartments)
-{
-    //
-    //       0
-    //       |
-    //       1
-    //       |
-    //       2
-    //      /|\.
-    //     3 7 4
-    //    /     \.
-    //   5       6
-    //
-    EXPECT_FALSE(
-        arb::algorithms::has_contiguous_compartments(
-            std::vector<int>{0, 0, 1, 2, 2, 3, 4, 2}
-        )
-    );
-
-    //
-    //       0
-    //       |
-    //       1
-    //       |
-    //       2
-    //      /|\.
-    //     3 6 5
-    //    /     \.
-    //   4       7
-    //
-    EXPECT_FALSE(
-        arb::algorithms::has_contiguous_compartments(
-            std::vector<int>{0, 0, 1, 2, 3, 2, 2, 5}
-        )
-    );
-
-    //
-    //       0
-    //       |
-    //       1
-    //       |
-    //       2
-    //      /|\.
-    //     3 7 5
-    //    /     \.
-    //   4       6
-    //
-    EXPECT_TRUE(
-        arb::algorithms::has_contiguous_compartments(
-            std::vector<int>{0, 0, 1, 2, 3, 2, 5, 2}
-        )
-    );
-
-    //
-    //         0
-    //         |
-    //         1
-    //        / \.
-    //       2   7
-    //      / \.
-    //     3   5
-    //    /     \.
-    //   4       6
-    //
-    EXPECT_TRUE(
-        arb::algorithms::has_contiguous_compartments(
-            std::vector<int>{0, 0, 1, 2, 3, 2, 5, 1}
-        )
-    );
-
-    //
-    //     0
-    //    / \.
-    //   1   2
-    //  / \.
-    // 3   4
-    //
-    EXPECT_TRUE(
-        arb::algorithms::has_contiguous_compartments(
-            std::vector<int>{0, 0, 0, 1, 1}
-        )
-    );
-
-    // Soma-only list
-    EXPECT_TRUE(
-        arb::algorithms::has_contiguous_compartments(
-            std::vector<int>{0}
-        )
-    );
-
-    // Empty list
-    EXPECT_TRUE(
-        arb::algorithms::has_contiguous_compartments(
-            std::vector<int>{}
-        )
-    );
-}
-
-TEST(algorithms, is_unique)
-{
-    EXPECT_TRUE(
-        arb::algorithms::is_unique(
-            std::vector<int>{}
-        )
-    );
-    EXPECT_TRUE(
-        arb::algorithms::is_unique(
-            std::vector<int>{0}
-        )
-    );
-    EXPECT_TRUE(
-        arb::algorithms::is_unique(
-            std::vector<int>{0,1,100}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_unique(
-            std::vector<int>{0,0}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_unique(
-            std::vector<int>{0,1,2,2,3,4}
-        )
-    );
-    EXPECT_FALSE(
-        arb::algorithms::is_unique(
-            std::vector<int>{0,1,2,3,4,4}
-        )
-    );
-}
-
-TEST(algorithms, child_count)
-{
-    {
-        //
-        //        0
-        //       /|\.
-        //      1 4 6
-        //     /  |  \.
-        //    2   5   7
-        //   /         \.
-        //  3           8
-        //             / \.
-        //            9   11
-        //           /     \.
-        //          10      12
-        //                   \.
-        //                    13
-        //
-        std::vector<int> parent_index =
-            { 0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12 };
-        std::vector<int> expected_child_count =
-            { 3, 1, 1, 0, 1, 0, 1, 1, 2, 1, 0, 1, 1, 0 };
-
-        // auto count = arb::algorithms::child_count(parent_index);
-        EXPECT_EQ(expected_child_count,
-                  arb::algorithms::child_count(parent_index));
-    }
-
-}
-
-TEST(algorithms, branches)
-{
-    using namespace arb;
-
-    {
-        //
-        //        0                        0
-        //       /|\.                     /|\.
-        //      1 4 6                    1 2 3
-        //     /  |  \.           =>        / \.
-        //    2   5   7                    4   5
-        //   /         \.
-        //  3           8
-        //             / \.
-        //            9   11
-        //           /     \.
-        //          10      12
-        //                   \.
-        //                    13
-        //
-        std::vector<int> parent_index =
-            { 0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12 };
-        std::vector<int> expected_branches =
-            { 0, 1, 4, 6, 9, 11, 14 };
-        std::vector<int> expected_parent_index =
-            { 0, 0, 0, 0, 3, 3 };
-
-        auto actual_branches = algorithms::branches(parent_index);
-        EXPECT_EQ(expected_branches, actual_branches);
-
-        auto actual_parent_index =
-            algorithms::tree_reduce(parent_index, actual_branches);
-        EXPECT_EQ(expected_parent_index, actual_parent_index);
-
-        // Check find_branch
-        EXPECT_EQ(0, algorithms::find_branch(actual_branches,  0));
-        EXPECT_EQ(1, algorithms::find_branch(actual_branches,  1));
-        EXPECT_EQ(1, algorithms::find_branch(actual_branches,  2));
-        EXPECT_EQ(1, algorithms::find_branch(actual_branches,  3));
-        EXPECT_EQ(2, algorithms::find_branch(actual_branches,  4));
-        EXPECT_EQ(2, algorithms::find_branch(actual_branches,  5));
-        EXPECT_EQ(3, algorithms::find_branch(actual_branches,  6));
-        EXPECT_EQ(3, algorithms::find_branch(actual_branches,  7));
-        EXPECT_EQ(3, algorithms::find_branch(actual_branches,  8));
-        EXPECT_EQ(4, algorithms::find_branch(actual_branches,  9));
-        EXPECT_EQ(4, algorithms::find_branch(actual_branches, 10));
-        EXPECT_EQ(5, algorithms::find_branch(actual_branches, 11));
-        EXPECT_EQ(5, algorithms::find_branch(actual_branches, 12));
-        EXPECT_EQ(5, algorithms::find_branch(actual_branches, 13));
-        EXPECT_EQ(6, algorithms::find_branch(actual_branches, 55));
-
-        // Check expand_branches
-        auto expanded = algorithms::expand_branches(actual_branches);
-        EXPECT_EQ(parent_index.size(), expanded.size());
-        for (std::size_t i = 0; i < parent_index.size(); ++i) {
-            EXPECT_EQ(algorithms::find_branch(actual_branches, i),
-                      expanded[i]);
-        }
-    }
-
-    {
-        //
-        //    0      0
-        //    |      |
-        //    1  =>  1
-        //    |
-        //    2
-        //    |
-        //    3
-        //
-        std::vector<int> parent_index          = { 0, 0, 1, 2 };
-        std::vector<int> expected_branches     = { 0, 1, 4 };
-        std::vector<int> expected_parent_index = { 0, 0 };
-
-        auto actual_branches = algorithms::branches(parent_index);
-        EXPECT_EQ(expected_branches, actual_branches);
-
-        auto actual_parent_index =
-            algorithms::tree_reduce(parent_index, actual_branches);
-        EXPECT_EQ(expected_parent_index, actual_parent_index);
-    }
-
-    {
-        //
-        //    0           0
-        //    |           |
-        //    1     =>    1
-        //    |          / \.
-        //    2         2   3
-        //   / \.
-        //  3   4
-        //       \.
-        //        5
-        //
-        std::vector<int> parent_index          = { 0, 0, 1, 2, 2, 4 };
-        std::vector<int> expected_branches     = { 0, 1, 3, 4, 6 };
-        std::vector<int> expected_parent_index = { 0, 0, 1, 1 };
-
-        auto actual_branches = algorithms::branches(parent_index);
-        EXPECT_EQ(expected_branches, actual_branches);
-
-        auto actual_parent_index =
-            algorithms::tree_reduce(parent_index, actual_branches);
-        EXPECT_EQ(expected_parent_index, actual_parent_index);
-    }
-
-    {
-        std::vector<int> parent_index          = { 0 };
-        std::vector<int> expected_branches     = { 0, 1 };
-        std::vector<int> expected_parent_index = { 0 };
-
-        auto actual_branches = algorithms::branches(parent_index);
-        EXPECT_EQ(expected_branches, actual_branches);
-
-        auto actual_parent_index =
-            algorithms::tree_reduce(parent_index, actual_branches);
-        EXPECT_EQ(expected_parent_index, actual_parent_index);
-    }
-}
-
-struct test_index_into {
-    template <typename R1, typename R2, typename R3>
-    bool operator() (const R1& sub, const R2& super, const R3& index) {
-        using value_type = typename R1::value_type;
-
-        if(sub.size()!=index.size()) return false;
-        auto index_it = index.begin();
-        for(auto i=0u; i<sub.size(); ++i, ++index_it) {
-            auto idx = *index_it;
-            if(idx>=value_type(super.size())) return false;
-            if(super[idx]!=sub[i]) return false;
-        }
-
-        return true;
-    }
-};
-
-template <typename Sub, typename Sup>
-::testing::AssertionResult validate_index_into(const Sub& sub, const Sup& sup) {
-    using namespace arb;
-
-    auto indices = util::index_into(sub, sup);
-    auto n_indices = std::size(indices);
-    auto n_sub  = std::size(sub);
-    if (std::size(indices)!=std::size(sub)) {
-        return ::testing::AssertionFailure()
-             << "index_into size " << n_indices << " does not equal sub-sequence size " << n_sub;
-    }
-
-    using std::begin;
-    using std::end;
-
-    auto sub_i = begin(sub);
-    auto sup_i = begin(sup);
-    auto sup_end = end(sup);
-    std::ptrdiff_t sup_idx = 0;
-
-    for (auto i: indices) {
-        if (sup_idx>i) {
-            return ::testing::AssertionFailure() << "indices in index_into sequence not monotonic";
-        }
-
-        while (sup_idx<i && sup_i!=sup_end) ++sup_idx, ++sup_i;
-
-        if (sup_i==sup_end) {
-            return ::testing::AssertionFailure() << "index " << i << "in index_into sequence is past the end";
-        }
-
-        if (!(*sub_i==*sup_i)) {
-            return ::testing::AssertionFailure()
-                << "value mismatch: sub-sequence element " << *sub_i
-                << " not equal to super-sequence element " << *sup_i << " at index " << i;
-        }
-
-        ++sub_i;
-    }
-
-    return ::testing::AssertionSuccess();
-}
-
-template <typename I>
-arb::util::range<std::reverse_iterator<I>> reverse_range(arb::util::range<I> r) {
-    using reviter = std::reverse_iterator<I>;
-    return arb::util::make_range(reviter(r.end()), reviter(r.begin()));
-}
-
-TEST(algorithms, index_into)
-{
-    using ivector = std::vector<std::ptrdiff_t>;
-    using std::size;
-    using arb::util::index_into;
-    using arb::util::assign_from;
-    using arb::util::make_range;
-    using arb::util::all_of;
-
-    std::vector<std::pair<std::vector<int>, std::vector<int>>> vector_tests = {
-        // Empty sequences:
-        {{}, {}},
-        {{100}, {}},
-        // Strictly monotonic sequences:
-        {{0, 7}, {0, 7}},
-        {{0, 1, 3, 4, 6, 7, 10, 11}, {0, 4, 6, 7, 11}},
-        {{0, 1, 3, 4, 6, 7, 10, 11}, {0}},
-        {{0, 1, 3, 4, 6, 7, 10, 11}, {11}},
-        {{0, 1, 3, 4, 6, 7, 10, 11}, {4}},
-        {{0, 1, 3, 4, 6, 7, 10, 11}, {0, 11}},
-        {{0, 1, 3, 4, 6, 7, 10, 11}, {4, 11}},
-        {{0, 1, 3, 4, 6, 7, 10, 11}, {0, 1, 3, 4, 6, 7, 10, 11}},
-        // Sequences with duplicates:
-        {{8, 8, 10, 10, 12, 12, 12, 13}, {8, 10, 13}},
-        {{8, 8, 10, 10, 12, 12, 12, 13}, {10, 10, 13}},
-        // Unordered sequences:
-        {{10, 3, 7, -8, 11, -8, 1, 2}, {3, -8, -8, 1}}
-    };
-
-    for (auto& testcase: vector_tests) {
-        EXPECT_TRUE(validate_index_into(testcase.second, testcase.first));
-    }
-
-    // Test across array types.
-
-    int subarr1[] = {2, 3, 9, 5};
-    int suparr1[] = {10, 2, 9, 3, 8, 5, 9, 3, 5, 10};
-    EXPECT_TRUE(validate_index_into(subarr1, suparr1));
-
-    // Test bidirectionality.
-
-    auto arr_indices = index_into(subarr1, suparr1);
-    ivector arridx = assign_from(arr_indices);
-
-    ivector revidx;
-    for (auto i = arr_indices.end(); i!=arr_indices.begin(); ) {
-        revidx.push_back(*--i);
-    }
-
-    std::vector<std::ptrdiff_t> expected(arridx);
-    std::reverse(expected.begin(), expected.end());
-    EXPECT_EQ(expected, revidx);
-
-    int subarr2[] = {8, 8, 8, 8, 8};
-    int suparr2[] = {8};
-
-    auto z_indices = index_into(subarr2, suparr2);
-    EXPECT_TRUE(all_of(z_indices, [](std::ptrdiff_t n) { return n==0; }));
-    EXPECT_EQ(0, z_indices.back());
-
-    // Test: strictly forward sequences; heterogenous sequences; sentinel-terminated ranges.
-
-    std::forward_list<double> sup_flist = {10.0, 2.1, 8.0, 3.8, 4.0, 4.0, 7.0, 1.0};
-    std::forward_list<int> sub_flist = {8, 4, 4, 1};
-
-    auto flist_indices = index_into(sub_flist, sup_flist);
-    ivector idx_flist = assign_from(flist_indices);
-    EXPECT_EQ((ivector{2, 4, 4, 7}), idx_flist);
-
-    const char* hello_world = "hello world";
-    const char* lol = "lol";
-    auto sup_cstr = make_range(hello_world, testing::null_terminated);
-    auto sub_cstr = make_range(lol, testing::null_terminated);
-    auto cstr_indices = index_into(sub_cstr, sup_cstr);
-    ivector idx_cstr = assign_from(cstr_indices);
-    EXPECT_EQ((ivector{2, 4, 9}), idx_cstr);
-}
-
-TEST(algorithms, binary_find)
-{
-    using arb::algorithms::binary_find;
-
-    // empty containers
-    {
-        std::vector<int> v;
-        EXPECT_TRUE(binary_find(v, 100) == std::end(v));
-    }
-
-    // value not present and greater than all entries
-    {
-        int a[] = {1, 10, 15};
-        EXPECT_TRUE(binary_find(a, 100) == std::end(a));
-
-        std::vector<int> v{1, 10, 15};
-        EXPECT_TRUE(binary_find(v, 100) == std::end(v));
-    }
-
-    // value not present and less than all entries
-    {
-        int a[] = {1, 10, 15};
-        EXPECT_TRUE(binary_find(a, -1) == std::end(a));
-
-        std::vector<int> v{1, 10, 15};
-        EXPECT_TRUE(binary_find(v, -1) == std::end(v));
-    }
-
-    // value not present and inside lower-upper bounds
-    {
-        int a[] = {1, 10, 15};
-        EXPECT_TRUE(binary_find(a, 4) == std::end(a));
-
-        std::vector<int> v{1, 10, 15};
-        EXPECT_TRUE(binary_find(v, 4) == std::end(v));
-    }
-
-    // value is first in range
-    {
-        int a[] = {1, 10, 15};
-        auto ita = binary_find(a, 1);
-        auto found = ita!=std::end(a);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(a), ita), 0);
-        if (found) {
-            EXPECT_EQ(*ita, 1);
-        }
-
-        std::vector<int> v{1, 10, 15};
-        auto itv = binary_find(v, 1);
-        found = itv!=std::end(v);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(v), itv), 0);
-        if (found) {
-            EXPECT_EQ(*itv, 1);
-        }
-    }
-
-    // value is last in range
-    {
-        int a[] = {1, 10, 15};
-        auto ita = binary_find(a, 15);
-        auto found = ita!=std::end(a);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(a), ita), 2);
-        if (found) {
-            EXPECT_EQ(*ita, 15);
-        }
-
-        std::vector<int> v{1, 10, 15};
-        auto itv = binary_find(v, 15);
-        found = itv!=std::end(v);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(v), itv), 2);
-        if (found) {
-            EXPECT_EQ(*itv, 15);
-        }
-    }
-
-    // value is last present and neither first nor last in range
-    {
-        int a[] = {1, 10, 15};
-        auto ita = binary_find(a, 10);
-        auto found = ita!=std::end(a);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(a), ita), 1);
-        if (found) {
-            EXPECT_EQ(*ita, 10);
-        }
-
-        std::vector<int> v{1, 10, 15};
-        auto itv = binary_find(v, 10);
-        found = itv!=std::end(v);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(v), itv), 1);
-        if (found) {
-            EXPECT_EQ(*itv, 10);
-        }
-    }
-
-    // value is last present and neither first nor last in range and range has even size
-    {
-        int a[] = {1, 10, 15, 27};
-        auto ita = binary_find(a, 10);
-        auto found = ita!=std::end(a);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(a), ita), 1);
-        if (found) {
-            EXPECT_EQ(*ita, 10);
-        }
-
-        std::vector<int> v{1, 10, 15, 27};
-        auto itv = binary_find(v, 10);
-        found = itv!=std::end(v);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::begin(v), itv), 1);
-        if (found) {
-            EXPECT_EQ(*itv, 10);
-        }
-    }
-
-    // test for const types
-    // i.e. iterators returned from passing in a const reference to a container
-    // can be compared to a const iterator from the container
-    {
-        std::vector<int> v{1, 10, 15};
-        auto const& vr = v;
-        auto itv = binary_find(vr, 10);
-        auto found = itv!=std::end(vr);
-        EXPECT_TRUE(found);
-        EXPECT_EQ(std::distance(std::cbegin(v), itv), 1);
-        if (found) {
-            EXPECT_EQ(*itv, 10);
-        }
-    }
-}
-
-struct int_string {
-    int value;
-
-    friend bool operator<(const int_string& lhs, const std::string& rhs) {
-        return lhs.value<std::stoi(rhs);
-    }
-    friend bool operator<(const std::string& lhs, const int_string& rhs) {
-        return std::stoi(lhs)<rhs.value;
-    }
-    friend bool operator==(const int_string& lhs, const std::string& rhs) {
-        return lhs.value==std::stoi(rhs);
-    }
-    friend bool operator==(const std::string& lhs, const int_string& rhs) {
-        return std::stoi(lhs)==rhs.value;
-    }
-};
-
-TEST(algorithms, binary_find_convert)
-{
-    using arb::algorithms::binary_find;
-
-    std::vector<std::string> values = {"0", "10", "20", "30"};
-    auto it = arb::algorithms::binary_find(values, int_string{20});
-
-    EXPECT_TRUE(it!=values.end());
-    EXPECT_TRUE(std::distance(values.begin(), it)==2u);
-}
diff --git a/test/unit/test_forest.cpp b/test/unit/test_forest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..18c82128632ad8b321dcd730e84385b3768470e4
--- /dev/null
+++ b/test/unit/test_forest.cpp
@@ -0,0 +1,242 @@
+#include <vector>
+
+#include "../gtest.h"
+
+#include <backends/gpu/forest.hpp>
+
+using namespace arb::gpu;
+TEST(forest, is_strictly_monotonic_increasing)
+{
+    EXPECT_TRUE(
+        is_strictly_monotonic_increasing(
+            std::vector<int>{0}
+        )
+    );
+    EXPECT_TRUE(
+        is_strictly_monotonic_increasing(
+            std::vector<int>{0, 1, 2, 3}
+        )
+    );
+    EXPECT_TRUE(
+        is_strictly_monotonic_increasing(
+            std::vector<int>{8, 20, 42, 89}
+        )
+    );
+    EXPECT_FALSE(
+        is_strictly_monotonic_increasing(
+            std::vector<int>{0, 0}
+        )
+    );
+    EXPECT_FALSE(
+        is_strictly_monotonic_increasing(
+            std::vector<int>{8, 20, 20, 89}
+        )
+    );
+    EXPECT_FALSE(
+        is_strictly_monotonic_increasing(
+            std::vector<int>{3, 2, 1, 0}
+        )
+    );
+}
+
+TEST(forest, has_contiguous_compartments)
+{
+    //
+    //       0
+    //       |
+    //       1
+    //       |
+    //       2
+    //      /|\.
+    //     3 7 4
+    //    /     \.
+    //   5       6
+    //
+    EXPECT_FALSE(
+        has_contiguous_compartments(
+            std::vector<int>{0, 0, 1, 2, 2, 3, 4, 2}
+        )
+    );
+
+    //
+    //       0
+    //       |
+    //       1
+    //       |
+    //       2
+    //      /|\.
+    //     3 6 5
+    //    /     \.
+    //   4       7
+    //
+    EXPECT_FALSE(
+        has_contiguous_compartments(
+            std::vector<int>{0, 0, 1, 2, 3, 2, 2, 5}
+        )
+    );
+
+    //
+    //       0
+    //       |
+    //       1
+    //       |
+    //       2
+    //      /|\.
+    //     3 7 5
+    //    /     \.
+    //   4       6
+    //
+    EXPECT_TRUE(
+        has_contiguous_compartments(
+            std::vector<int>{0, 0, 1, 2, 3, 2, 5, 2}
+        )
+    );
+
+    //
+    //         0
+    //         |
+    //         1
+    //        / \.
+    //       2   7
+    //      / \.
+    //     3   5
+    //    /     \.
+    //   4       6
+    //
+    EXPECT_TRUE(
+        has_contiguous_compartments(
+            std::vector<int>{0, 0, 1, 2, 3, 2, 5, 1}
+        )
+    );
+
+    //
+    //     0
+    //    / \.
+    //   1   2
+    //  / \.
+    // 3   4
+    //
+    EXPECT_TRUE(
+        has_contiguous_compartments(
+            std::vector<int>{0, 0, 0, 1, 1}
+        )
+    );
+
+    // Soma-only list
+    EXPECT_TRUE(
+        has_contiguous_compartments(
+            std::vector<int>{0}
+        )
+    );
+
+    // Empty list
+    EXPECT_TRUE(
+        has_contiguous_compartments(
+            std::vector<int>{}
+        )
+    );
+}
+
+TEST(forest, branches) {
+    using namespace arb;
+    {
+        //
+        //        0                        0
+        //       /|\.                     /|\.
+        //      1 4 6                    1 2 3
+        //     /  |  \.           =>        / \.
+        //    2   5   7                    4   5
+        //   /         \.
+        //  3           8
+        //             / \.
+        //            9   11
+        //           /     \.
+        //          10      12
+        //                   \.
+        //                    13
+        //
+        std::vector<int> parent_index =
+            { 0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12 };
+        std::vector<int> expected_branches =
+            { 0, 1, 4, 6, 9, 11, 14 };
+        std::vector<int> expected_parent_index =
+            { 0, 0, 0, 0, 3, 3 };
+
+        auto actual_branches = branches(parent_index);
+        EXPECT_EQ(expected_branches, actual_branches);
+
+        auto actual_parent_index =
+            tree_reduce(parent_index, actual_branches);
+        EXPECT_EQ(expected_parent_index, actual_parent_index);
+
+        // Check expand_branches
+        auto expanded = gpu::expand_branches(actual_branches);
+        EXPECT_EQ(parent_index.size(), expanded.size());
+        for (std::size_t i = 0; i < parent_index.size(); ++i) {
+            auto it =  std::find_if(
+                actual_branches.begin(), actual_branches.end(),
+                [i](const int& v) { return v > (int)i; }
+            );
+            auto branch = it - actual_branches.begin() - 1;
+            EXPECT_EQ(branch,expanded[i]);
+        }
+    }
+    {
+        //
+        //    0      0
+        //    |      |
+        //    1  =>  1
+        //    |
+        //    2
+        //    |
+        //    3
+        //
+        std::vector<int> parent_index          = { 0, 0, 1, 2 };
+        std::vector<int> expected_branches     = { 0, 1, 4 };
+        std::vector<int> expected_parent_index = { 0, 0 };
+
+        auto actual_branches = gpu::branches(parent_index);
+        EXPECT_EQ(expected_branches, actual_branches);
+
+        auto actual_parent_index =
+            gpu::tree_reduce(parent_index, actual_branches);
+        EXPECT_EQ(expected_parent_index, actual_parent_index);
+    }
+
+    {
+        //
+        //    0           0
+        //    |           |
+        //    1     =>    1
+        //    |          / \.
+        //    2         2   3
+        //   / \.
+        //  3   4
+        //       \.
+        //        5
+        //
+        std::vector<int> parent_index          = { 0, 0, 1, 2, 2, 4 };
+        std::vector<int> expected_branches     = { 0, 1, 3, 4, 6 };
+        std::vector<int> expected_parent_index = { 0, 0, 1, 1 };
+
+        auto actual_branches = gpu::branches(parent_index);
+        EXPECT_EQ(expected_branches, actual_branches);
+
+        auto actual_parent_index =
+            gpu::tree_reduce(parent_index, actual_branches);
+        EXPECT_EQ(expected_parent_index, actual_parent_index);
+    }
+
+    {
+        std::vector<int> parent_index          = { 0 };
+        std::vector<int> expected_branches     = { 0, 1 };
+        std::vector<int> expected_parent_index = { 0 };
+
+        auto actual_branches = gpu::branches(parent_index);
+        EXPECT_EQ(expected_branches, actual_branches);
+
+        auto actual_parent_index =
+            gpu::tree_reduce(parent_index, actual_branches);
+        EXPECT_EQ(expected_parent_index, actual_parent_index);
+    }
+}
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 82b6bba4576243cef1ed183dbcd17f7670aa965c..19650d69ebfd4adeca375bb3d7bf4320302bb0e7 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -18,7 +18,6 @@
 
 #include <arborenv/concurrency.hpp>
 
-#include "algorithms.hpp"
 #include "backends/multicore/fvm.hpp"
 #include "backends/multicore/mechanism.hpp"
 #include "execution_context.hpp"
diff --git a/test/unit/test_index.cpp b/test/unit/test_index.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..36a08bbb64087923f6a3fbfb2cf2b2a13b8c1b39
--- /dev/null
+++ b/test/unit/test_index.cpp
@@ -0,0 +1,131 @@
+#include <forward_list>
+#include <random>
+#include <vector>
+
+#include "../gtest.h"
+
+#include "util/index_into.hpp"
+#include "util/rangeutil.hpp"
+
+#include "common.hpp"
+
+template <typename Sub, typename Sup>
+::testing::AssertionResult validate_index_into(const Sub& sub, const Sup& sup) {
+    using namespace arb;
+
+    auto indices = util::index_into(sub, sup);
+    auto n_indices = std::size(indices);
+    auto n_sub  = std::size(sub);
+    if (std::size(indices)!=std::size(sub)) {
+        return ::testing::AssertionFailure()
+             << "index_into size " << n_indices << " does not equal sub-sequence size " << n_sub;
+    }
+
+    using std::begin;
+    using std::end;
+
+    auto sub_i = begin(sub);
+    auto sup_i = begin(sup);
+    auto sup_end = end(sup);
+    std::ptrdiff_t sup_idx = 0;
+
+    for (auto i: indices) {
+        if (sup_idx>i) {
+            return ::testing::AssertionFailure() << "indices in index_into sequence not monotonic";
+        }
+
+        while (sup_idx<i && sup_i!=sup_end) ++sup_idx, ++sup_i;
+
+        if (sup_i==sup_end) {
+            return ::testing::AssertionFailure() << "index " << i << "in index_into sequence is past the end";
+        }
+
+        if (!(*sub_i==*sup_i)) {
+            return ::testing::AssertionFailure()
+                << "value mismatch: sub-sequence element " << *sub_i
+                << " not equal to super-sequence element " << *sup_i << " at index " << i;
+        }
+
+        ++sub_i;
+    }
+
+    return ::testing::AssertionSuccess();
+}
+
+TEST(util, index_into)
+{
+    using ivector = std::vector<std::ptrdiff_t>;
+    using std::size;
+    using arb::util::index_into;
+    using arb::util::assign_from;
+    using arb::util::make_range;
+    using arb::util::all_of;
+
+    std::vector<std::pair<std::vector<int>, std::vector<int>>> vector_tests = {
+        // Empty sequences:
+        {{}, {}},
+        {{100}, {}},
+        // Strictly monotonic sequences:
+        {{0, 7}, {0, 7}},
+        {{0, 1, 3, 4, 6, 7, 10, 11}, {0, 4, 6, 7, 11}},
+        {{0, 1, 3, 4, 6, 7, 10, 11}, {0}},
+        {{0, 1, 3, 4, 6, 7, 10, 11}, {11}},
+        {{0, 1, 3, 4, 6, 7, 10, 11}, {4}},
+        {{0, 1, 3, 4, 6, 7, 10, 11}, {0, 11}},
+        {{0, 1, 3, 4, 6, 7, 10, 11}, {4, 11}},
+        {{0, 1, 3, 4, 6, 7, 10, 11}, {0, 1, 3, 4, 6, 7, 10, 11}},
+        // Sequences with duplicates:
+        {{8, 8, 10, 10, 12, 12, 12, 13}, {8, 10, 13}},
+        {{8, 8, 10, 10, 12, 12, 12, 13}, {10, 10, 13}},
+        // Unordered sequences:
+        {{10, 3, 7, -8, 11, -8, 1, 2}, {3, -8, -8, 1}}
+    };
+
+    for (auto& testcase: vector_tests) {
+        EXPECT_TRUE(validate_index_into(testcase.second, testcase.first));
+    }
+
+    // Test across array types.
+
+    int subarr1[] = {2, 3, 9, 5};
+    int suparr1[] = {10, 2, 9, 3, 8, 5, 9, 3, 5, 10};
+    EXPECT_TRUE(validate_index_into(subarr1, suparr1));
+
+    // Test bidirectionality.
+
+    auto arr_indices = index_into(subarr1, suparr1);
+    ivector arridx = assign_from(arr_indices);
+
+    ivector revidx;
+    for (auto i = arr_indices.end(); i!=arr_indices.begin(); ) {
+        revidx.push_back(*--i);
+    }
+
+    std::vector<std::ptrdiff_t> expected(arridx);
+    std::reverse(expected.begin(), expected.end());
+    EXPECT_EQ(expected, revidx);
+
+    int subarr2[] = {8, 8, 8, 8, 8};
+    int suparr2[] = {8};
+
+    auto z_indices = index_into(subarr2, suparr2);
+    EXPECT_TRUE(all_of(z_indices, [](std::ptrdiff_t n) { return n==0; }));
+    EXPECT_EQ(0, z_indices.back());
+
+    // Test: strictly forward sequences; heterogenous sequences; sentinel-terminated ranges.
+
+    std::forward_list<double> sup_flist = {10.0, 2.1, 8.0, 3.8, 4.0, 4.0, 7.0, 1.0};
+    std::forward_list<int> sub_flist = {8, 4, 4, 1};
+
+    auto flist_indices = index_into(sub_flist, sup_flist);
+    ivector idx_flist = assign_from(flist_indices);
+    EXPECT_EQ((ivector{2, 4, 4, 7}), idx_flist);
+
+    const char* hello_world = "hello world";
+    const char* lol = "lol";
+    auto sup_cstr = make_range(hello_world, testing::null_terminated);
+    auto sub_cstr = make_range(lol, testing::null_terminated);
+    auto cstr_indices = index_into(sub_cstr, sup_cstr);
+    ivector idx_cstr = assign_from(cstr_indices);
+    EXPECT_EQ((ivector{2, 4, 9}), idx_cstr);
+}
\ No newline at end of file
diff --git a/test/unit/test_matrix_cpuvsgpu.cpp b/test/unit/test_matrix_cpuvsgpu.cpp
index ed751b21a2d957cc4b189b6e5788dba005cab147..9212f3b0c95a968d4eac054065d15c64a3ae47e3 100644
--- a/test/unit/test_matrix_cpuvsgpu.cpp
+++ b/test/unit/test_matrix_cpuvsgpu.cpp
@@ -4,7 +4,6 @@
 
 #include <arbor/math.hpp>
 
-#include "algorithms.hpp"
 #include "matrix.hpp"
 #include "memory/memory.hpp"
 #include "util/span.hpp"
diff --git a/test/unit/test_matrix_gpu.cpp b/test/unit/test_matrix_gpu.cpp
index 66499801e207ac36cd5cce4467734736f0b87d93..a2f8e2fd46c46103669de088b202865a230ca8a3 100644
--- a/test/unit/test_matrix_gpu.cpp
+++ b/test/unit/test_matrix_gpu.cpp
@@ -12,7 +12,6 @@
 
 #include <arbor/math.hpp>
 
-#include "algorithms.hpp"
 #include "matrix.hpp"
 #include "memory/memory.hpp"
 #include "util/span.hpp"
diff --git a/test/unit/test_tree.cpp b/test/unit/test_tree.cpp
index 094e77f1e18e5a98f2aa095797c3b73edd67be8a..0dc6b91e3df6d85474a6ca6c3490a067ba276633 100644
--- a/test/unit/test_tree.cpp
+++ b/test/unit/test_tree.cpp
@@ -1,6 +1,3 @@
-#include <fstream>
-#include <iostream>
-#include <numeric>
 #include <vector>
 
 #include "../gtest.h"
@@ -11,6 +8,72 @@ using namespace arb;
 using int_type = tree::int_type;
 using iarray = tree::iarray;
 
+TEST(tree, minimal_degree)
+{
+    {
+        std::vector<int> v = {0};
+        EXPECT_TRUE(is_minimal_degree(v));
+    }
+
+    {
+        std::vector<int> v = {0, 0, 1, 2, 3, 4};
+        EXPECT_TRUE(is_minimal_degree(v));
+    }
+
+    {
+        std::vector<int> v = {0, 0, 1, 2, 0, 4};
+        EXPECT_TRUE(is_minimal_degree(v));
+    }
+
+    {
+        std::vector<int> v = {0, 0, 1, 2, 0, 4, 5, 4};
+        EXPECT_TRUE(is_minimal_degree(v));
+    }
+
+    {
+        std::vector<int> v = {1};
+        EXPECT_FALSE(is_minimal_degree(v));
+    }
+
+    {
+        std::vector<int> v = {0, 2};
+        EXPECT_FALSE(is_minimal_degree(v));
+    }
+
+    {
+        std::vector<int> v = {0, 1, 2};
+        EXPECT_FALSE(is_minimal_degree(v));
+    }
+}
+
+TEST(tree, child_count)
+{
+    {
+        //
+        //        0
+        //       /|\.
+        //      1 4 6
+        //     /  |  \.
+        //    2   5   7
+        //   /         \.
+        //  3           8
+        //             / \.
+        //            9   11
+        //           /     \.
+        //          10      12
+        //                   \.
+        //                    13
+        //
+        std::vector<int> parent_index =
+            { 0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12 };
+        std::vector<int> expected_child_count =
+            { 3, 1, 1, 0, 1, 0, 1, 1, 2, 1, 0, 1, 1, 0 };
+
+        EXPECT_EQ(expected_child_count,child_count(parent_index));
+    }
+
+}
+
 TEST(tree, from_segment_index) {
     auto no_parent = tree::no_parent;
 
@@ -256,4 +319,3 @@ TEST(tree, depth_from_root) {
         EXPECT_EQ(expected, depth_from_root(tree(parent_index)));
     }
 }
-