From 135e91acc1388e2658dcc70fba03d21e919a9163 Mon Sep 17 00:00:00 2001
From: Sam Yates <halfflat@gmail.com>
Date: Thu, 27 Aug 2020 14:38:05 +0200
Subject: [PATCH] Merge pull request #1123 from
halfflat/feature/stitch-morphology
New PR post master rollback; squashed and rebased, but reprises #1111.
* Add (forward) ordered forest implementation, tests.
* Add `segment` region expression; to ease implementation, `msegment` now also knows its own id.
* Add `stitch_builder` and `stitched_morphology`. A stitch corresponds to a labelled, linearly-interpolated segment which can be attached at any point along a parent stitch. A `stitched_morphology` takes a `stitch_builder` object and constructs a segment tree and morphology, and provides a dictionary of stitch labels to segment indices and region expressions.
* Add `import` method for `label_dict`, so that the label dictionary returned by `stitched_morphology` can be merged with an existing dictionary.
* Add section on stitch builder etc. to cable cell docs.
* Update cable cell docs to remove out of date info and to provide some context.
* Describe ordered forest datastructure and interface in a long comment at the beginning of the header file.
---
arbor/CMakeLists.txt | 1 +
arbor/include/arbor/morph/embed_pwlin.hpp | 19 +-
arbor/include/arbor/morph/label_dict.hpp | 1 +
arbor/include/arbor/morph/morphexcept.hpp | 21 +
arbor/include/arbor/morph/morphology.hpp | 2 +-
arbor/include/arbor/morph/primitives.hpp | 2 +-
arbor/include/arbor/morph/region.hpp | 3 +
arbor/include/arbor/morph/stitch.hpp | 93 +++
arbor/morph/embed_pwlin.cpp | 23 +-
arbor/morph/label_dict.cpp | 9 +
arbor/morph/locset.cpp | 2 +-
arbor/morph/morphexcept.cpp | 33 +-
arbor/morph/primitives.cpp | 6 +-
arbor/morph/region.cpp | 58 +-
arbor/morph/segment_tree.cpp | 2 +-
arbor/morph/stitch.cpp | 214 +++++++
arbor/util/ordered_forest.hpp | 668 ++++++++++++++++++++++
doc/cpp_cable_cell.rst | 204 +++++--
test/unit/CMakeLists.txt | 2 +
test/unit/test_morph_embedding.cpp | 76 +--
test/unit/test_morph_expr.cpp | 44 +-
test/unit/test_morph_stitch.cpp | 236 ++++++++
test/unit/test_ordered_forest.cpp | 513 +++++++++++++++++
23 files changed, 2074 insertions(+), 158 deletions(-)
create mode 100644 arbor/include/arbor/morph/stitch.hpp
create mode 100644 arbor/morph/stitch.cpp
create mode 100644 arbor/util/ordered_forest.hpp
create mode 100644 test/unit/test_morph_stitch.cpp
create mode 100644 test/unit/test_ordered_forest.cpp
diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 55a23170..0174af4b 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -31,6 +31,7 @@ set(arbor_sources
memory/gpu_wrappers.cpp
memory/util.cpp
morph/embed_pwlin.cpp
+ morph/stitch.cpp
morph/label_dict.cpp
morph/locset.cpp
morph/morphexcept.cpp
diff --git a/arbor/include/arbor/morph/embed_pwlin.hpp b/arbor/include/arbor/morph/embed_pwlin.hpp
index ac065005..412e8997 100644
--- a/arbor/include/arbor/morph/embed_pwlin.hpp
+++ b/arbor/include/arbor/morph/embed_pwlin.hpp
@@ -22,14 +22,17 @@ using pw_constant_fn = util::pw_elements<double>;
struct embed_pwlin {
explicit embed_pwlin(const arb::morphology& m);
- // Locations that mark the boundaries between segments.
- const std::vector<mlocation>& segment_locations() const {
- return segment_locations_;
+ // Segment queries.
+ msize_t num_segments() const {
+ return segment_cables_.size();
}
- std::pair<const mlocation*, const mlocation*> branch_segment_locations(msize_t i) const {
- const mlocation* p = segment_locations_.data();
- return std::make_pair(p+branch_segment_part_[i], p+branch_segment_part_[i+1]);
+ mcable segment(msize_t seg_id) const {
+ return segment_cables_.at(seg_id);
+ }
+
+ const mlocation_list& segment_ends() const {
+ return all_segment_ends_;
}
// Interpolated radius at location.
@@ -65,8 +68,8 @@ struct embed_pwlin {
}
private:
- std::vector<mlocation> segment_locations_;
- std::vector<msize_t> branch_segment_part_;
+ mlocation_list all_segment_ends_;
+ std::vector<mcable> segment_cables_;
std::shared_ptr<embed_pwlin_data> data_;
};
diff --git a/arbor/include/arbor/morph/label_dict.hpp b/arbor/include/arbor/morph/label_dict.hpp
index 852cef71..f15630d6 100644
--- a/arbor/include/arbor/morph/label_dict.hpp
+++ b/arbor/include/arbor/morph/label_dict.hpp
@@ -16,6 +16,7 @@ class label_dict {
reg_map regions_;
public:
+ void import(const label_dict& other);
void set(const std::string& name, locset ls);
void set(const std::string& name, region reg);
diff --git a/arbor/include/arbor/morph/morphexcept.hpp b/arbor/include/arbor/morph/morphexcept.hpp
index bd7872c8..13963c20 100644
--- a/arbor/include/arbor/morph/morphexcept.hpp
+++ b/arbor/include/arbor/morph/morphexcept.hpp
@@ -41,6 +41,27 @@ struct invalid_segment_parent: morphology_error {
msize_t tree_size;
};
+struct duplicate_stitch_id: morphology_error {
+ duplicate_stitch_id(const std::string& id);
+ std::string id;
+};
+
+struct no_such_stitch: morphology_error {
+ no_such_stitch(const std::string& id);
+ std::string id;
+};
+
+struct missing_stitch_start: morphology_error {
+ missing_stitch_start(const std::string& id);
+ std::string id;
+};
+
+struct invalid_stitch_position: morphology_error {
+ invalid_stitch_position(const std::string& id, double along);
+ std::string id;
+ double along;
+};
+
struct label_type_mismatch: morphology_error {
label_type_mismatch(const std::string& label);
std::string label;
diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index a4a6b622..18da8ee3 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -4,9 +4,9 @@
#include <ostream>
#include <vector>
-#include <arbor/util/lexcmp_def.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/morph/segment_tree.hpp>
+#include <arbor/util/lexcmp_def.hpp>
namespace arb {
diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp
index cf75c0d4..919730a1 100644
--- a/arbor/include/arbor/morph/primitives.hpp
+++ b/arbor/include/arbor/morph/primitives.hpp
@@ -46,11 +46,11 @@ enum class comp_op {
// Describe a cable segment between two adjacent samples.
struct msegment {
+ msize_t id;
mpoint prox;
mpoint dist;
int tag;
- friend bool operator==(const msegment&, const msegment&);
friend std::ostream& operator<<(std::ostream&, const msegment&);
};
diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp
index ebdedff4..f875dcdf 100644
--- a/arbor/include/arbor/morph/region.hpp
+++ b/arbor/include/arbor/morph/region.hpp
@@ -133,6 +133,9 @@ region branch(msize_t);
// Region with all segments with segment tag id.
region tagged(int id);
+// Region corresponding to a single segment.
+region segment(int id);
+
// Region up to `distance` distal from points in `start`.
region distal_interval(locset start, double distance);
diff --git a/arbor/include/arbor/morph/stitch.hpp b/arbor/include/arbor/morph/stitch.hpp
new file mode 100644
index 00000000..1d0287f9
--- /dev/null
+++ b/arbor/include/arbor/morph/stitch.hpp
@@ -0,0 +1,93 @@
+#pragma once
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/label_dict.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/util/optional.hpp>
+
+namespace arb {
+
+// Stiches represent an alternative building block for morphologies.
+//
+// A stitch describes a portion of the morphology delimited by two
+// `mpoint`s. Stitches can be attached to a parent stich at any
+// point along the parent, interpolated linearly from the end points.
+// Each stitch is associated with a unique string label, and optionally
+// an integer tag value.
+//
+// The stitch builder collects stitches and produces the corresponding
+// morphology and region/location labels.
+
+struct mstitch {
+ std::string id;
+ util::optional<mpoint> prox;
+ mpoint dist;
+ int tag;
+
+ mstitch(std::string id, mpoint prox, mpoint dist, int tag = 0):
+ id(std::move(id)), prox(std::move(prox)), dist(std::move(dist)), tag(tag)
+ {}
+
+ mstitch(std::string id, mpoint dist, int tag = 0):
+ id(std::move(id)), dist(std::move(dist)), tag(tag)
+ {}
+};
+
+struct stitch_builder_impl;
+struct stitched_morphology;
+
+struct stitch_builder {
+ stitch_builder();
+
+ stitch_builder(const stitch_builder&) = delete;
+ stitch_builder(stitch_builder&&);
+
+ stitch_builder& operator=(const stitch_builder&) = delete;
+ stitch_builder& operator=(stitch_builder&&);
+
+ // Make a new stitch in the morphology, return reference to self.
+ //
+ // If the stitch does not contained a proximal point, it will be
+ // inferred from the point where it attaches to the parent stitch.
+ // If the parent is omitted, it will be taken to be the last stitch
+ // added.
+
+ stitch_builder& add(mstitch f, const std::string& parent_id, double along = 1.);
+ stitch_builder& add(mstitch f, double along = 1.);
+
+ ~stitch_builder();
+private:
+ friend stitched_morphology;
+ std::unique_ptr<stitch_builder_impl> impl_;
+};
+
+// From stitch builder construct morphology, region expressions.
+
+struct stitched_morphology_impl;
+
+struct stitched_morphology {
+ stitched_morphology() = delete;
+ stitched_morphology(const stitch_builder&); // implicit
+ stitched_morphology(stitch_builder&&); // implicit
+
+ stitched_morphology(const stitched_morphology&) = delete;
+ stitched_morphology(stitched_morphology&&);
+
+ arb::morphology morphology() const;
+ region stitch(const std::string& id) const;
+ std::vector<msize_t> segments(const std::string& id) const;
+
+ // Create labeled regions for each stitch with label equal to the stitch id, prepended by `prefix`.
+ label_dict labels(const std::string& prefix="") const;
+
+ ~stitched_morphology();
+private:
+ std::unique_ptr<stitched_morphology_impl> impl_;
+};
+
+} // namesapce arb
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index 2e381804..ef5ab650 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/arbor/morph/embed_pwlin.cpp
@@ -201,12 +201,9 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
constexpr double pi = math::pi<double>;
msize_t n_branch = m.num_branches();
data_ = std::make_shared<embed_pwlin_data>(n_branch);
- branch_segment_part_.reserve(n_branch+1);
if (!n_branch) return;
- branch_segment_part_.push_back(0);
-
double proj_shift = m.branch_segments(0).front().prox.z;
for (msize_t bid = 0; bid<n_branch; ++bid) {
@@ -214,16 +211,15 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
auto& segments = m.branch_segments(bid);
arb_assert(segments.size());
- branch_segment_part_.push_back(branch_segment_part_.back()+segments.size()+1);
-
std::vector<double> seg_pos;
seg_pos.reserve(segments.size()+1);
seg_pos.push_back(0.);
- for (auto &seg: segments) {
+ for (const auto &seg: segments) {
double d = distance(seg.prox, seg.dist);
seg_pos.push_back(seg_pos.back()+d);
}
+
double branch_length = seg_pos.back();
double length_scale = branch_length>0? 1./branch_length: 0;
for (auto& d: seg_pos) {
@@ -232,9 +228,22 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
seg_pos.back() = 1; // Circumvent any rounding infelicities.
for (auto d: seg_pos) {
- segment_locations_.push_back({bid, d});
+ all_segment_ends_.push_back({bid, d});
+ }
+
+ // Second pass over segments to store associated cables.
+ auto pos_iter = seg_pos.begin();
+ for (const auto &seg: segments) {
+ double pos0 = *pos_iter++;
+ double pos1 = *pos_iter;
+
+ if (seg.id>=segment_cables_.size()) {
+ segment_cables_.resize(seg.id+1);
+ }
+ segment_cables_[seg.id] = mcable{bid, pos0, pos1};
}
+
double length_0 = parent==mnpos? 0: data_->length[parent].back().second[1];
data_->length[bid].push_back(0., 1, rat_element<1, 0>(length_0, length_0+branch_length));
diff --git a/arbor/morph/label_dict.cpp b/arbor/morph/label_dict.cpp
index 9b00acbc..b62fa563 100644
--- a/arbor/morph/label_dict.cpp
+++ b/arbor/morph/label_dict.cpp
@@ -25,6 +25,15 @@ void label_dict::set(const std::string& name, arb::region reg) {
regions_[name] = std::move(reg);
}
+void label_dict::import(const label_dict& other) {
+ for (const auto& entry: other.locsets()) {
+ set(entry.first, entry.second);
+ }
+ for (const auto& entry: other.regions()) {
+ set(entry.first, entry.second);
+ }
+}
+
util::optional<const region&> label_dict::region(const std::string& name) const {
auto it = regions_.find(name);
if (it==regions_.end()) return {};
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index 365c66a6..67fce43b 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -141,7 +141,7 @@ locset segment_boundaries() {
}
mlocation_list thingify_(const segments_&, const mprovider& p) {
- return p.embedding().segment_locations();
+ return p.embedding().segment_ends();
}
std::ostream& operator<<(std::ostream& o, const segments_& x) {
diff --git a/arbor/morph/morphexcept.cpp b/arbor/morph/morphexcept.cpp
index ac9222bb..bac840de 100644
--- a/arbor/morph/morphexcept.cpp
+++ b/arbor/morph/morphexcept.cpp
@@ -29,7 +29,6 @@ no_such_segment::no_such_segment(msize_t id):
sid(id)
{}
-
invalid_mcable::invalid_mcable(mcable cable):
morphology_error(pprintf("invalid mcable {}", cable)),
cable(cable)
@@ -39,6 +38,33 @@ invalid_mcable_list::invalid_mcable_list():
morphology_error("bad mcable_list")
{}
+invalid_segment_parent::invalid_segment_parent(msize_t parent, msize_t tree_size):
+ morphology_error(pprintf("invalid segment parent {} for a segment tree of size {}", msize_string(parent), tree_size)),
+ parent(parent),
+ tree_size(tree_size)
+{}
+
+duplicate_stitch_id::duplicate_stitch_id(const std::string& id):
+ morphology_error(pprintf("duplicate stitch id {}", id)),
+ id(id)
+{}
+
+no_such_stitch::no_such_stitch(const std::string& id):
+ morphology_error(pprintf("no such stitch id {}", id)),
+ id(id)
+{}
+
+missing_stitch_start::missing_stitch_start(const std::string& id):
+ morphology_error(pprintf("require proximal point for stitch id {}", id)),
+ id(id)
+{}
+
+invalid_stitch_position::invalid_stitch_position(const std::string& id, double along):
+ morphology_error(pprintf("invalid stitch position {} on stitch {}", along, id)),
+ id(id),
+ along(along)
+{}
+
label_type_mismatch::label_type_mismatch(const std::string& label):
morphology_error(pprintf("label \"{}\" is already bound to a different type of object", label)),
label(label)
@@ -59,11 +85,6 @@ circular_definition::circular_definition(const std::string& name):
name(name)
{}
-invalid_segment_parent::invalid_segment_parent(msize_t parent, msize_t tree_size):
- morphology_error(pprintf("invalid segment parent {} for a segment tree of size {}", msize_string(parent), tree_size)),
- parent(parent),
- tree_size(tree_size)
-{}
} // namespace arb
diff --git a/arbor/morph/primitives.cpp b/arbor/morph/primitives.cpp
index ced745ab..d4b780fe 100644
--- a/arbor/morph/primitives.cpp
+++ b/arbor/morph/primitives.cpp
@@ -132,16 +132,12 @@ bool test_invariants(const mcable_list& l) {
&& l.end()==std::find_if(l.begin(), l.end(), [](const mcable& c) {return !test_invariants(c);});
}
-bool operator==(const msegment& l, const msegment& r) {
- return l.prox==r.prox && l.dist==r.dist && l.tag==r.tag;
-}
-
std::ostream& operator<<(std::ostream& o, const mpoint& p) {
return o << "(point " << p.x << " " << p.y << " " << p.z << " " << p.radius << ")";
}
std::ostream& operator<<(std::ostream& o, const msegment& s) {
- return o << "(segment " << s.prox << " " << s.dist << " " << s.tag << ")";
+ return o << "(segment " << s.id << " " << s.prox << " " << s.dist << " " << s.tag << ")";
}
std::ostream& operator<<(std::ostream& o, const mlocation& l) {
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index 6658ed45..68ff60d3 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -140,42 +140,52 @@ region tagged(int id) {
}
mextent thingify_(const tagged_& reg, const mprovider& p) {
- const auto& m = p.morphology();
const auto& e = p.embedding();
- size_t nb = m.num_branches();
+ const auto& m = p.morphology();
- std::vector<mcable> L;
- L.reserve(nb);
- auto matches = [reg](auto& seg){return seg.tag==reg.tag;};
- auto not_matches = [&matches](auto& seg) {return !matches(seg);};
+ size_t nb = m.num_branches();
+ std::vector<mcable> cables;
for (msize_t i: util::make_span(nb)) {
- auto& segs = m.branch_segments(i); // Range of segments in the branch.
- auto locs = util::make_range(e.branch_segment_locations(i));
+ for (const auto& seg: m.branch_segments(i)) {
+ if (seg.tag==reg.tag) {
+ cables.push_back(e.segment(seg.id));
+ }
+ }
+ }
+ util::sort(cables);
+ return mextent(cables);
+}
- auto beg = std::cbegin(segs);
- auto end = std::cend(segs);
+std::ostream& operator<<(std::ostream& o, const tagged_& t) {
+ return o << "(tag " << t.tag << ")";
+}
- // Find the next section that matches reg.tag.
- auto start = std::find_if(beg, end, [reg](auto& seg){return seg.tag==reg.tag;});
- while (start!=end) {
- // find the next sample that does not match reg.tag
- auto last = std::find_if(start, end, not_matches);
+// Region comprising a single segment.
- auto l = locs[start-beg].pos;
- auto r = last==end? 1: locs[last-beg].pos;
- L.push_back({i, l, r});
+struct segment_: region_tag {
+ explicit segment_(int id): id(id) {}
+ int id;
+};
- // Find the next sample in the branch that matches reg.tag.
- start = std::find_if(last, end, matches);
- }
+region segment(int id) {
+ return region(segment_{id});
+}
+
+mextent thingify_(const segment_& reg, const mprovider& p) {
+ const auto& e = p.embedding();
+
+ msize_t id(reg.id);
+ if (id>=e.num_segments()) {
+ throw no_such_segment(id);
}
- return mextent(L);
+ mcable_list cables = {e.segment(id)};
+ return mextent(cables);
}
-std::ostream& operator<<(std::ostream& o, const tagged_& t) {
- return o << "(tag " << t.tag << ")";
+std::ostream& operator<<(std::ostream& o, const segment_& reg) {
+ return o << "(segment " << reg.id << ")";
}
// Region comprising whole morphology.
diff --git a/arbor/morph/segment_tree.cpp b/arbor/morph/segment_tree.cpp
index 8629cef5..d8fca989 100644
--- a/arbor/morph/segment_tree.cpp
+++ b/arbor/morph/segment_tree.cpp
@@ -24,7 +24,7 @@ msize_t segment_tree::append(msize_t p, const mpoint& prox, const mpoint& dist,
}
auto id = size();
- segments_.push_back(msegment{prox, dist, tag});
+ segments_.push_back(msegment{id, prox, dist, tag});
parents_.push_back(p);
// Set the point properties for the new point, and update those of the parent.
diff --git a/arbor/morph/stitch.cpp b/arbor/morph/stitch.cpp
new file mode 100644
index 00000000..fafd8460
--- /dev/null
+++ b/arbor/morph/stitch.cpp
@@ -0,0 +1,214 @@
+#include <memory>
+#include <unordered_map>
+#include <vector>
+
+#include <arbor/morph/stitch.hpp>
+#include <arbor/morph/morphexcept.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/util/optional.hpp>
+
+#include "util/ordered_forest.hpp"
+#include "util/maputil.hpp"
+#include "util/meta.hpp"
+#include "util/rangeutil.hpp"
+#include "util/transform.hpp"
+
+namespace arb {
+
+struct stitch_builder_impl {
+ struct stitch_segment {
+ double along_prox;
+ double along_dist;
+
+ mpoint prox;
+ mpoint dist;
+ int tag;
+ std::string stitch_id;
+ msize_t seg_id;
+ };
+
+ using forest_type = util::ordered_forest<stitch_segment>;
+
+ forest_type forest;
+ std::unordered_map<std::string, forest_type::iterator> id_to_node;
+ std::string last_id;
+
+ stitch_builder_impl() = default;
+ stitch_builder_impl(stitch_builder_impl&&) = default;
+ stitch_builder_impl(const stitch_builder_impl& other):
+ forest(other.forest),
+ last_id(other.last_id)
+ {
+ for (auto i = forest.preorder_begin(); i!=forest.preorder_end(); ++i) {
+ id_to_node.insert({i->stitch_id, i});
+ }
+ }
+
+ stitch_builder_impl& operator=(stitch_builder_impl&&) = default;
+ stitch_builder_impl& operator=(const stitch_builder_impl& other) {
+ if (this==&other) return *this;
+ return *this = stitch_builder_impl(other);
+ }
+
+ void add(mstitch f, const std::string& parent, double along) {
+ if (id_to_node.count(f.id)) throw duplicate_stitch_id(f.id);
+
+ forest_type::iterator p;
+
+ if (!(parent.empty() && forest.empty())) {
+ p = find_stitch_along(parent, along);
+ arb_assert(p);
+
+ if (along==p->along_prox) {
+ if (!f.prox) f.prox = p->prox;
+ p = p.parent();
+ }
+ else if (along<p->along_dist) {
+ // Split parent node p at along.
+ auto split = *p;
+
+ mpoint point = lerp(p->prox, p->dist, (along-p->along_prox)/(p->along_dist-p->along_prox));
+ if (!f.prox) f.prox = point;
+
+ p->dist = point;
+ p->along_dist = along;
+ split.prox = point;
+ split.along_prox = along;
+
+ auto i = forest.push_child(p, split);
+ while (i.next()) {
+ auto tmp = forest.prune_after(i);
+ forest.graft_child(i, std::move(tmp));
+ }
+ }
+ else {
+ if (!f.prox) f.prox = p->dist;
+ }
+ }
+ if (!f.prox) throw missing_stitch_start(f.id);
+
+ stitch_segment n{0., 1., f.prox.value(), f.dist, f.tag, f.id, msize_t(-1)};
+ id_to_node[f.id] = p? forest.push_child(p, n): forest.push_front(n);
+ last_id = f.id;
+ }
+
+ forest_type::iterator find_stitch_along(const std::string& id, double along) {
+ if (along<0 || along>1) throw invalid_stitch_position(id, along);
+
+ auto map_it = id_to_node.find(id);
+ if (map_it==id_to_node.end()) throw no_such_stitch(id);
+
+ auto i = map_it->second;
+ arb_assert(i->along_prox==0);
+ arb_assert(i->along_dist==1 || i.child());
+
+ while (along>i->along_dist) {
+ // Continuation is last child.
+ i = i.child();
+ arb_assert(i);
+ while (i.next()) i = i.next();
+ }
+ return i;
+ }
+};
+
+stitch_builder::stitch_builder(): impl_(new stitch_builder_impl) {}
+
+stitch_builder::stitch_builder(stitch_builder&&) = default;
+stitch_builder& stitch_builder::operator=(stitch_builder&&) = default;
+
+stitch_builder& stitch_builder::add(mstitch f, const std::string& parent_id, double along) {
+ impl_->add(std::move(f), parent_id, along);
+ return *this;
+}
+
+stitch_builder& stitch_builder::add(mstitch f, double along) {
+ return add(std::move(f), impl_->last_id, along);
+}
+
+stitch_builder::~stitch_builder() = default;
+
+
+struct stitched_morphology_impl {
+ std::unordered_multimap<std::string, msize_t> id_to_segs;
+ segment_tree stree;
+
+ stitched_morphology_impl(stitch_builder_impl bimpl) {
+ auto iter = bimpl.forest.preorder_begin();
+ auto end = bimpl.forest.preorder_end();
+
+ for (; iter!=end; ++iter) {
+ msize_t seg_parent_id = iter.parent()? iter.parent()->seg_id: mnpos;
+ iter->seg_id = stree.append(seg_parent_id, iter->prox, iter->dist, iter->tag);
+ }
+
+ for (const auto& id_node: bimpl.id_to_node) {
+ const std::string& id = id_node.first;
+ auto iter = id_node.second;
+
+ while (iter && iter->stitch_id==id) {
+ id_to_segs.insert({id, iter->seg_id});
+ iter = iter.child();
+ while (iter.next()) {
+ iter = iter.next();
+ }
+ }
+ }
+ }
+};
+
+stitched_morphology::stitched_morphology(stitch_builder&& builder):
+ impl_(new stitched_morphology_impl(std::move(*builder.impl_)))
+{}
+
+stitched_morphology::stitched_morphology(const stitch_builder& builder):
+ impl_(new stitched_morphology_impl(*builder.impl_))
+{}
+
+stitched_morphology::stitched_morphology(stitched_morphology&& other) = default;
+
+arb::morphology stitched_morphology::morphology() const {
+ return arb::morphology(impl_->stree);
+}
+
+label_dict stitched_morphology::labels(const std::string& prefix) const {
+ label_dict dict;
+
+ auto i0 = impl_->id_to_segs.begin();
+ auto end = impl_->id_to_segs.end();
+ while (i0 != end) {
+ auto i1 = i0;
+ while (i1 != end && i1->first==i0->first) ++i1;
+
+ region r = util::foldl(
+ [&](region r, const auto& elem) { return join(std::move(r), reg::segment(elem.second)); },
+ reg::nil(),
+ util::make_range(i0, i1));
+
+ dict.set(prefix+i0->first, std::move(r));
+ i0 = i1;
+ }
+
+ return dict;
+}
+
+region stitched_morphology::stitch(const std::string& id) const {
+ auto seg_ids = util::make_range(impl_->id_to_segs.equal_range(id));
+ if (seg_ids.empty()) throw no_such_stitch(id);
+
+ return util::foldl(
+ [&](region r, const auto& elem) { return join(std::move(r), reg::segment(elem.second)); },
+ reg::nil(), seg_ids);
+}
+
+std::vector<msize_t> stitched_morphology::segments(const std::string& id) const {
+ auto seg_ids = util::transform_view(util::make_range(impl_->id_to_segs.equal_range(id)), util::second);
+ if (seg_ids.empty()) throw no_such_stitch(id);
+
+ return std::vector<msize_t>(begin(seg_ids), end(seg_ids));
+}
+
+stitched_morphology::~stitched_morphology() = default;
+
+} // namespace arb
diff --git a/arbor/util/ordered_forest.hpp b/arbor/util/ordered_forest.hpp
new file mode 100644
index 00000000..567226c4
--- /dev/null
+++ b/arbor/util/ordered_forest.hpp
@@ -0,0 +1,668 @@
+#pragma once
+
+// The ordered forest data structure represents zero or more trees with an
+// arbitrary number of children per node, where the trees and children have a
+// well-defined ordering.
+//
+// This implementation presents a 'forward' interface, along the lines of
+// `std::forward_list`. Given an iterator pointing to a particular node in the
+// forest, one can directly obtain an iterator to its parent, its first child,
+// or its next sibling, but not its previous sibling.
+//
+// Common iterator methods:
+//
+// `parent()`: return an iterator to the node's parent.
+// `next()`: return an iterator to the node's next sibling.
+// `child()`: return an iterator to the node's first child.
+// `preorder_next()`: return an iterator to the next node in pre-order.
+// `postorder_next()`: return an iterator to the next node in post-order.
+//
+// Note that these methods are not covariant: they return 'base' iterators
+// that can be implicitly converted to any of the derived iterators described
+// below.
+//
+// A default-constructed iterator points to no node, and plays the role of
+// an `end` iterator.
+//
+// The derived iterators differ only in the semantics of `operator++`:
+//
+// `sibling_iterator`: proceed to next sibling.
+// `preorder_iterator`: proceed to next node in pre-order.
+// `postorder_iterator`: proceed to next node in post-order.
+//
+// Any of the iterators can be converted implicitly to any const iterator, and
+// any of the non-const iterators can be converted implicitly to any non-const
+// iterator.
+//
+// The iterators provided by `begin()` and `end()` are pre-order iterators.
+// The other ordered forest iterator methods are:
+//
+// `child_begin`: return a sibling iterator starting at the child of the given node.
+// `root_begin`: return a sibling iterator to the first tree.
+// `postorder_begin`: return a post-order iterator to the first node in post-order.
+// `preorder_begin`: return a post-order iterator to the first node in pre-order.
+//
+// Note that all ordered forest iterators are forward iterators.
+//
+// Ordered forests can be constructed from an intializer list with the aid of
+// the `ordered_forest_builder` helper class. Each item in the initalizer list
+// is either a value (corresponding to a node without children), or a pair
+// { value, IL } where value denotes the value at the node, and IL is an
+// initializer list with the same semantics, defining the node's children.
+//
+// Mutation operations are of the following types:
+//
+// * Insertions: add a new node as next sibling, first child, or first tree.
+// * Grafts: splice a forest between a node and its next sibling, or before its first child.
+// * Erasures: remove a node, replacing it with its children.
+// * Cuts: remove a sub-tree, presenting it as a new forest.
+//
+// These methods are described in more detail within the class definition below.
+//
+// The ordered forest implementation is allocator aware.
+
+#include <type_traits>
+#include <iterator>
+#include <memory>
+#include <utility>
+
+namespace arb {
+namespace util {
+
+template <typename V, typename Allocator>
+struct ordered_forest_builder;
+
+template <typename V, typename Allocator = std::allocator<V>>
+struct ordered_forest {
+private:
+ struct node {
+ V* item_ = nullptr;
+ node* parent_ = nullptr;
+ node* child_ = nullptr;
+ node* next_ = nullptr;
+ };
+
+ using node_alloc_t = typename std::allocator_traits<Allocator>::template rebind_alloc<node>;
+ using item_alloc_traits = std::allocator_traits<Allocator>;
+ using node_alloc_traits = std::allocator_traits<node_alloc_t>;
+
+public:
+ using value_type = V;
+ using allocator_type = Allocator;
+ using size_type = std::size_t;
+
+ struct iterator_base {
+ node* n_ = nullptr;
+
+ iterator_base() = default;
+ explicit iterator_base(node* n): n_(n) {}
+ explicit operator bool() const { return n_; }
+ };
+
+ // Constructor overloads for iterators permit construction of any const
+ // iterator from any other iterator, or of any mutable iterator from any
+ // other mutable iterator.
+
+ template <bool const_flag>
+ struct iterator_mc: iterator_base {
+ using pointer = std::conditional_t<const_flag, const V*, V*>;
+ using reference = std::conditional_t<const_flag, const V&, V&>;
+ using value_type = V;
+ using difference_type = std::ptrdiff_t;
+ using iterator_category = std::forward_iterator_tag;
+
+ iterator_mc() = default;
+
+ template <bool flag = const_flag, typename std::enable_if_t<flag, int> = 0>
+ iterator_mc(const iterator_mc<false>& i): iterator_base(i.n_) {}
+
+ iterator_mc parent() const { return iterator_mc{n_? n_->parent_: nullptr}; }
+ iterator_mc next() const { return iterator_mc{n_? n_->next_: nullptr}; }
+ iterator_mc child() const { return iterator_mc{n_? n_->child_: nullptr}; }
+
+ bool operator==(const iterator_base& a) const { return n_ == a.n_; }
+ bool operator!=(const iterator_base& a) const { return n_ != a.n_; }
+
+ iterator_mc preorder_next() const {
+ if (!n_) return {};
+ if (n_->child_) return iterator_mc{n_->child_};
+
+ node* x = n_;
+ while (x && !x->next_) x = x->parent_;
+ return iterator_mc{x? x->next_: nullptr};
+ }
+
+ iterator_mc postorder_next() const {
+ if (!n_) return {};
+ if (n_->next_) {
+ node* x = n_->next_;
+ while (node* c = x->child_) x = c;
+ return iterator_mc{x};
+ }
+ else return parent();
+ }
+
+ reference operator*() const { return *n_->item_; }
+ pointer operator->() const { return n_->item_; }
+
+ protected:
+ friend ordered_forest;
+
+ explicit iterator_mc(node* n): iterator_base(n) {}
+ using iterator_base::n_;
+ };
+
+ template <bool const_flag>
+ struct sibling_iterator_mc: iterator_mc<const_flag> {
+ sibling_iterator_mc() = default;
+ sibling_iterator_mc(const iterator_mc<const_flag>& i): iterator_mc<const_flag>(i) {}
+
+ sibling_iterator_mc& operator++() { return *this = this->next(); }
+ sibling_iterator_mc operator++(int) { auto p = *this; return ++*this, p; }
+ };
+
+ using sibling_iterator = sibling_iterator_mc<false>;
+ using const_sibling_iterator = sibling_iterator_mc<true>;
+
+ template <bool const_flag>
+ struct preorder_iterator_mc: iterator_mc<const_flag> {
+ preorder_iterator_mc() = default;
+ preorder_iterator_mc(const iterator_mc<const_flag>& i): iterator_mc<const_flag>(i) {}
+
+ preorder_iterator_mc& operator++() { return *this = this->preorder_next(); }
+ preorder_iterator_mc operator++(int) { auto p = *this; return ++*this, p; }
+ };
+
+ using preorder_iterator = preorder_iterator_mc<false>;
+ using const_preorder_iterator = preorder_iterator_mc<true>;
+
+ template <bool const_flag>
+ struct postorder_iterator_mc: iterator_mc<const_flag> {
+ postorder_iterator_mc() = default;
+ postorder_iterator_mc(const iterator_mc<const_flag>& i): iterator_mc<const_flag>(i) {}
+
+ postorder_iterator_mc& operator++() { return *this = this->postorder_next(); }
+ postorder_iterator_mc operator++(int) { auto p = *this; return ++*this, p; }
+ };
+
+ using postorder_iterator = postorder_iterator_mc<false>;
+ using const_postorder_iterator = postorder_iterator_mc<true>;
+
+ bool empty() const { return !first_; }
+
+ // Note: size is O(n) in the number of elements n.
+ std::size_t size() const {
+ std::size_t n = 0;
+ for (const auto& _: *this) (void)_, ++n;
+ return n;
+ }
+
+ sibling_iterator child_begin(const iterator_mc<false>& i) { return sibling_iterator{i.child()}; }
+ const_sibling_iterator child_begin(const iterator_mc<true>& i) const { return const_sibling_iterator{i.child()}; }
+
+ sibling_iterator child_end(const iterator_base&) { return {}; }
+ const_sibling_iterator child_end(const iterator_base&) const { return {}; }
+
+ sibling_iterator root_begin() { return sibling_iterator{first_else_end()}; }
+ const_sibling_iterator root_begin() const { return const_sibling_iterator{first_else_end()}; }
+
+ sibling_iterator root_end() { return sibling_iterator{}; }
+ const_sibling_iterator root_end() const { return const_sibling_iterator{}; }
+
+ postorder_iterator postorder_begin() { return postorder_iterator{first_leaf()}; }
+ const_postorder_iterator postorder_begin() const { return const_postorder_iterator{first_leaf()}; }
+
+ postorder_iterator postorder_end() { return {}; }
+ const_postorder_iterator postorder_end() const { return {}; }
+
+ preorder_iterator preorder_begin() { return preorder_iterator{first_else_end()}; }
+ const_preorder_iterator preorder_begin() const { return const_preorder_iterator{first_else_end()}; }
+
+ preorder_iterator preorder_end() { return {}; }
+ const_preorder_iterator preorder_end() const { return {}; }
+
+ // Default iteration is preorder.
+
+ using iterator = preorder_iterator;
+ using const_iterator = const_preorder_iterator;
+
+ iterator begin() { return preorder_begin(); }
+ iterator end() { return {}; }
+
+ const_iterator begin() const { return preorder_begin(); }
+ const_iterator end() const { return {}; }
+
+ const_iterator cbegin() const { return preorder_begin(); }
+ const_iterator cend() const { return {}; }
+
+ // Insertion and emplace operations:
+ //
+ // * All return an iterator to the last inserted node, or an iterator to the referenced
+ // node (or first tree, for push_front) if the collection of inserted items is empty.
+ //
+ // * The iterator argument may not be an end iterator.
+ //
+ // Insertion member functions are templated on iterator class in order to return covariantly.
+
+ // Insert/emplace item as next sibling.
+
+ template <typename Iter, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter insert_after(const Iter& i, const V& item) {
+ return assert_valid(i), splice_impl(i.n_->parent_, i.n_->next_, make_node(item));
+ }
+
+ template <typename Iter, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter insert_after(const Iter& i, V&& item) {
+ return assert_valid(i), splice_impl(i.n_->parent_, i.n_->next_, make_node(std::move(item)));
+ }
+
+ template <typename Iter, typename... Args, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter emplace_after(const Iter& i, Args&&... args) {
+ return assert_valid(i), splice_impl(i.n_->parent_, i.n_->next_, make_node(std::forward<Args>(args)...));
+ }
+
+ // Insert trees in forest as next siblings.
+
+ template <typename Iter, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter graft_after(const Iter& i, ordered_forest of) {
+ assert_valid(i);
+ if (of.empty()) return i;
+
+ // Underlying move or copy depends upon allocator equality.
+ ordered_forest f(std::move(of), get_allocator());
+ node* sp_first = f.first_;
+ f.first_ = nullptr;
+ return splice_impl(i.n_->parent_, i.n_->next_, sp_first);
+ }
+
+ // Insert item as first child.
+
+ template <typename Iter, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter push_child(const Iter& i, const V& item) {
+ return assert_valid(i), splice_impl(i.n_, i.n_->child_, make_node(item));
+ }
+
+ template <typename Iter, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter push_child(const Iter& i, V&& item) {
+ return assert_valid(i), splice_impl(i.n_, i.n_->child_, make_node(std::move(item)));
+ }
+
+ template <typename Iter, typename... Args, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter emplace_child(const Iter& i, Args&&... args) {
+ return assert_valid(i), splice_impl(i.n_, i.n_->child_, make_node(std::forward<Args>(args)...));
+ }
+
+ // Insert trees in forest as first children.
+
+ template <typename Iter, typename = std::enable_if_t<std::is_base_of<iterator_mc<false>, Iter>::value>>
+ Iter graft_child(const Iter& i, ordered_forest of) {
+ assert_valid(i);
+ if (of.empty()) return i;
+
+ // Underlying move or copy depends upon allocator equality.
+ ordered_forest f(std::move(of), get_allocator());
+ node* sp_first = f.first_;
+ f.first_ = nullptr;
+ return splice_impl(i.n_, i.n_->child_, sp_first);
+ }
+
+ // Insert item as first top-level tree.
+
+ iterator push_front(const V& item) {
+ return splice_impl(nullptr, first_, make_node(item));
+ }
+
+ iterator push_front(V&& item) {
+ return splice_impl(nullptr, first_, make_node(std::move(item)));
+ }
+
+ template <typename... Args>
+ iterator emplace_front(Args&&... args) {
+ return splice_impl(nullptr, first_, make_node(std::forward<Args>(args)...));
+ }
+
+ // Insert trees in forest as first top-level children.
+
+ iterator graft_front(ordered_forest of) {
+ if (of.empty()) return {};
+
+ // Underlying move or copy depends upon allocator equality.
+ ordered_forest f(std::move(of), get_allocator());
+ node* sp_first = f.first_;
+ f.first_ = nullptr;
+ return splice_impl(nullptr, first_, sp_first);
+ }
+
+ // Erase and cut operations:
+ //
+ // * Erase/pop operations replace a node with all of that node's children.
+ // * Prune operations remove a whole subtree, and return it as a new ordered forest.
+
+ // Erase/cut next sibling.
+
+ void erase_after(const iterator_mc<false>& i) {
+ assert_valid(i.next()), erase_impl(i.n_->next_);
+ }
+
+ ordered_forest prune_after(const iterator_mc<false>& i) {
+ return assert_valid(i.next()), prune_impl(i.n_->next_);
+ }
+
+ // Erase/cut first child.
+
+ void erase_child(const iterator_mc<false>& i) {
+ assert_valid(i.child()), erase_impl(i.n_->child_);
+ }
+
+ ordered_forest prune_child(const iterator_mc<false>& i) {
+ return assert_valid(i.child()), prune_impl(i.n_->child_);
+ }
+
+ // Erase/cut root of first tree. Precondition: forest is non-empty.
+
+ void erase_front() {
+ assert_nonempty(), erase_impl(first_);
+ }
+
+ ordered_forest prune_front() {
+ return assert_nonempty(), prune_impl(first_);
+ }
+
+ // Access by reference to root of first tree.
+
+ V& front() { return *begin(); }
+ const V& front() const { return *begin(); }
+
+ // Comparison:
+
+ bool operator==(const ordered_forest& other) const {
+ const_iterator a = begin();
+ const_iterator b = other.begin();
+
+ while (a && b) {
+ if (!a.child() != !b.child() || !a.next() != !b.next() || !(*a==*b)) return false;
+ ++a, ++b;
+ }
+ return !a && !b;
+ }
+
+ bool operator!=(const ordered_forest& other) const {
+ return !(*this==other);
+ }
+
+ // Constructors, assignment, destructors:
+
+ ordered_forest(const Allocator& alloc = Allocator()) noexcept:
+ item_alloc_(alloc),
+ node_alloc_(alloc)
+ {}
+
+ ordered_forest(ordered_forest&& other) noexcept:
+ ordered_forest(std::move(other.item_alloc_))
+ {
+ first_ = other.first_;
+ other.first_ = nullptr;
+ }
+
+ ordered_forest(ordered_forest&& other, const Allocator& alloc):
+ ordered_forest(alloc)
+ {
+ if (allocators_equal(other)) {
+ first_ = other.first_;
+ other.first_ = nullptr;
+ }
+ else {
+ copy_impl(other);
+ }
+ }
+
+ ordered_forest(std::initializer_list<ordered_forest_builder<V, Allocator>> blist, const Allocator& alloc = Allocator{}):
+ ordered_forest(alloc)
+ {
+ sibling_iterator j;
+ for (auto& b: blist) {
+ ordered_forest f(std::move(b.f_), item_alloc_);
+ j = j? graft_after(j, std::move(f)): sibling_iterator(graft_front(std::move(f)));
+ }
+ }
+
+ ordered_forest(const ordered_forest& other):
+ ordered_forest(item_alloc_traits::select_on_container_copy_construction(other.item_alloc_))
+ {
+ copy_impl(other);
+ }
+
+ ordered_forest(const ordered_forest& other, const Allocator& alloc):
+ ordered_forest(alloc)
+ {
+ copy_impl(other);
+ }
+
+ ordered_forest& operator=(const ordered_forest& other) {
+ if (this==&other) return *this;
+ delete_node(first_);
+
+ if (item_alloc_traits::propagate_on_container_copy_assignment::value) {
+ item_alloc_ = other.item_alloc_;
+ }
+ if (node_alloc_traits::propagate_on_container_copy_assignment::value) {
+ node_alloc_ = other.node_alloc_;
+ }
+
+ copy_impl(other);
+ return *this;
+ }
+
+ ordered_forest& operator=(ordered_forest&& other) {
+ if (this==&other) return *this;
+ delete_node(first_);
+
+ if (item_alloc_traits::propagate_on_container_move_assignment::value) {
+ item_alloc_ = other.item_alloc_;
+ }
+ if (node_alloc_traits::propagate_on_container_move_assignment::value) {
+ node_alloc_ = other.node_alloc_;
+ }
+
+ if (allocators_equal(other)) {
+ first_ = other.first_;
+ other.first_ = nullptr;
+ }
+ else {
+ copy_impl(other);
+ }
+ return *this;
+ }
+
+ ~ordered_forest() { delete_node(first_); }
+
+ // Swap
+
+ void swap(ordered_forest& other)
+ noexcept(
+ (item_alloc_traits::propagate_on_container_swap::value && node_alloc_traits::propagate_on_container_swap::value) ||
+ (item_alloc_traits::is_always_equal::value && node_alloc_traits::is_always_equal::value)
+ )
+ {
+ using std::swap;
+ if (item_alloc_traits::propagate_on_container_swap::value) {
+ swap(item_alloc_, other.item_alloc_);
+ }
+ if (node_alloc_traits::propagate_on_container_swap::value) {
+ swap(node_alloc_, other.node_alloc_);
+ }
+
+ swap(first_, other.first_);
+ }
+
+ friend void swap(ordered_forest& a, ordered_forest& b) {
+ a.swap(b);
+ }
+
+ // Allocator access
+
+ allocator_type get_allocator() const noexcept { return item_alloc_; }
+
+private:
+ Allocator item_alloc_;
+ node_alloc_t node_alloc_;
+ node* first_ = nullptr;
+
+ bool allocators_equal(const ordered_forest& other) const {
+ return item_alloc_==other.item_alloc_ && node_alloc_==other.node_alloc_;
+ }
+
+ iterator_mc<false> first_else_end() { return iterator_mc<false>{first_}; }
+ iterator_mc<true> first_else_end() const { return iterator_mc<true>{first_}; }
+
+ iterator_mc<false> first_leaf() { return iterator_mc<false>{first_leaf_node()}; }
+ iterator_mc<true> first_leaf() const { return iterator_mc<true>{first_leaf_node()}; }
+
+ node* first_leaf_node() const {
+ node* n = first_;
+ while (n && n->child_) n = n->child_;
+ return n;
+ }
+
+ // Throw on invalid iterator.
+ void assert_valid(iterator_base i) {
+ if (!i.n_) throw std::invalid_argument("bad iterator");
+ }
+
+ void assert_nonempty() {
+ if (!first_) throw std::invalid_argument("empty forest");
+ }
+
+ ordered_forest prune_impl(node*& next_write) {
+ node* r = next_write;
+ next_write = next_write->next_;
+ r->next_ = nullptr;
+ r->parent_ = nullptr;
+
+ ordered_forest f(get_allocator());
+ f.first_ = r;
+
+ return f;
+ }
+
+ void erase_impl(node*& next_write) {
+ ordered_forest cut = prune_impl(next_write);
+
+ node* cut_root = cut.first_;
+ if (cut_root->child_) {
+ splice_impl(next_write->parent_, next_write, cut_root->child_);
+ }
+
+ cut_root->child_ = nullptr;
+ }
+
+ iterator_mc<false> splice_impl(node* parent, node*& next_write, node* sp_first) {
+ node* sp_last = nullptr;
+
+ for (node* j = sp_first; j; j = j->next_) {
+ j->parent_ = parent;
+ sp_last = j;
+ }
+ sp_last->next_ = next_write;
+ next_write = sp_first;
+
+ return iterator_mc<false>{sp_last};
+ }
+
+ template <typename U, typename OtherAllocator>
+ void copy_impl(const ordered_forest<U, OtherAllocator>& other) {
+ auto copy_children = [&](auto& self, const auto& from, auto& to) -> void {
+ sibling_iterator j;
+ for (auto i = other.child_begin(from); i!=other.child_end(from); ++i) {
+ // TODO: explicit `this` required for g++6; remove when g++6 deprecated.
+ j = j? this->insert_after(j, *i): sibling_iterator(this->push_child(to, *i));
+ self(self, i, j);
+ }
+ };
+
+ sibling_iterator j;
+ for (auto i = other.root_begin(); i!=other.root_end(); ++i) {
+ j = j? insert_after(j, *i): sibling_iterator(push_front(*i));
+ copy_children(copy_children, i, j);
+ }
+ }
+
+ // Node creation, destruction:
+
+ template <typename... Args>
+ node* make_node(Args&&... args) {
+ node* x = node_alloc_traits::allocate(node_alloc_, 1);
+ try {
+ node_alloc_traits::construct(node_alloc_, x);
+ }
+ catch (...) {
+ node_alloc_traits::deallocate(node_alloc_, x, 1);
+ throw;
+ }
+
+ x->item_ = item_alloc_traits::allocate(item_alloc_, 1);
+ try {
+ item_alloc_traits::construct(item_alloc_, x->item_, std::forward<Args>(args)...);
+ }
+ catch (...) {
+ item_alloc_traits::deallocate(item_alloc_, x->item_, 1);
+ throw;
+ }
+ return x;
+ }
+
+ void delete_node(node* n) {
+ if (!n) return;
+
+ delete_item(n->item_);
+ delete_node(n->child_);
+ delete_node(n->next_);
+
+ node_alloc_traits::destroy(node_alloc_, n);
+ node_alloc_traits::deallocate(node_alloc_, n, 1);
+ }
+
+ void delete_item(V* item) {
+ if (!item) return;
+
+ item_alloc_traits::destroy(item_alloc_, item);
+ item_alloc_traits::deallocate(item_alloc_, item, 1);
+ }
+};
+
+// Helper class for building trees from initializer_lists. Ordered forest can be
+// constructed from an initializer list of builder objects; each builder object
+// represents a tree, constructed from a single value, or a pair: root value
+// and a sequence of builder objects represeting the children.
+//
+// While the builder takes the same allocator type as its associated forest, it
+// will build the temporary trees with a default-constructed allocator.
+
+template <typename V, typename Allocator>
+struct ordered_forest_builder {
+ template <typename X, typename std::enable_if_t<std::is_constructible<V, X&&>::value, int> = 0>
+ ordered_forest_builder(X&& x) {
+ f_.emplace_front(x);
+ }
+
+ template <typename X, typename std::enable_if_t<std::is_constructible<V, X&&>::value, int> = 0>
+ ordered_forest_builder(X&& x, std::initializer_list<ordered_forest_builder<V, Allocator>> children) {
+ using sibling_iterator = typename ordered_forest<V, Allocator>::sibling_iterator;
+
+ auto top = f_.emplace_front(x);
+ sibling_iterator j;
+
+ for (auto& g: children) {
+ ordered_forest<V, Allocator> c(std::move(g.f_));
+ j = j? f_.graft_after(j, std::move(c)): sibling_iterator(f_.graft_child(top, std::move(c)));
+ }
+ }
+
+ friend class ordered_forest<V, Allocator>;
+
+private:
+ ordered_forest<V, Allocator> f_;
+};
+
+} // namespace util
+} // namespace arb
diff --git a/doc/cpp_cable_cell.rst b/doc/cpp_cable_cell.rst
index f60e322d..8083fe7e 100644
--- a/doc/cpp_cable_cell.rst
+++ b/doc/cpp_cable_cell.rst
@@ -5,8 +5,8 @@ Cable cells
.. Warning::
The interface for building and modifying cable cell objects
- will be thoroughly revised in the near future. The documentation
- here is primarily a place holder.
+ has changed significantly; some of the documentation below is
+ out of date.
Cable cells, which use the :cpp:enum:`cell_kind` :cpp:expr:`cable`,
represent morphologically-detailed neurons as 1-d trees, with
@@ -21,38 +21,140 @@ object of type :cpp:type:`cable_cell_global_properties`.
The :cpp:type:`cable_cell` object
---------------------------------
-Cable cells are built up from a series of :cpp:type:`segment`
-objects, which themselves describe an unbranched component of the
-cell morphology. These segments are added via the methods:
+Cable cells are constructed from a :cpp:type:`morphology` and, optionally, a
+:cpp:type:`label_dict` that associates names with particular points
+(:cpp:type:`locset` objects) or subsets (:cpp:type:`region` objects) of the
+morphology.
-.. cpp:function:: soma_segment* cable_cell::add_soma(double radius)
+Morphologies are constructed from a :cpp:type:`segment_tree`, but can also
+be generated via the :cpp:type:`stitch_builder`, which offers a slightly
+higher level interface. Details are described in :ref:`morphology-construction`.
- Add the soma to the cable cell with the given radius. There
- can be only one per cell.
+Each cell has particular values for its electrical and ionic properties. These
+are determined first by the set of global defaults, then the defaults
+associated with the cell, and finally by any values specified explicitly for a
+given subsection of the morphology via the ``paint`` method
+(see :ref:`electrical-properties` and :ref:`paint-properties`).
- The soma segment has index 0, and must be added before any
- cable segments.
+Ion channels and other distributed dynamical processes are also specified
+on the cell via the ``paint`` method; while synapses, current clamps,
+gap junction locations, and the site for testing the threshold potential
+are specified via the ``place`` method. See :ref:`cable-cell-dynamics`, below.
-.. cpp:function:: cable_segment* cable_cell::add_cable(cell_lid_type index, Args&&... args)
+.. _morphology-construction:
- Add a unbranched section of the cell morphology, with its proximal
- end attached to the segment given by :cpp:expr:`index`. The
- following arguments are forwarded to the :cpp:type:`cable_segment`
- constructor.
-
-Segment indices are exactly the order in which they have been added
-to a cell, counting from zero (for the soma). Both :cpp:type:`soma_segment`
-and :cpp:type:`cable_segment` are derived from the abstract base
-class :cpp:type:`segment`.
+Constucting cell morphologies
+-----------------------------
.. todo::
- Describe cable_segment constructor arguments, unless we get to the
- replace cell building/morphology implementation first.
+ TODO: Description of segment trees is in the works.
+
+
+The stitch-builder interface
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Like the segment tree, the :cpp:type:`stich_builder` class constructs morphologies
+through attaching simple components described by a pair of :cpp:type:`mpoint` values,
+proximal and distal. These components are :cpp:type:`mstitch` objects, and
+they differ from segments in two regards:
+
+1. Stitches are identified by a unique string identifier, in addition to an optional tag value.
+
+2. Stitches can be attached to a parent stitch at either end, or anywhere in the middle.
+
+The ability to attach a stitch some way along another stitch dictates that one
+stitch may correspond to more than one morphological segment once the morphology
+is fully specified. When these attachment points are internal to a stitch, the
+corresponding geometrical point is determined by linearly interpolating between
+the proximal and distal points.
+
+The required header file is ``arbor/morph/stitch.hpp``.
+
+:cpp:type:`mstitch` has two constructors:
+
+.. code::
+
+ mstitch::mstitch(std::string id, mpoint prox, mpoint dist, int tag = 0)
+ mstitch::mstitch(std::string id, mpoint dist, int tag = 0)
+
+If the proximal point is omitted, it will be inferred from the point at which
+the stitch is attached to its parent.
+
+The :cpp:type:`stitch_builder` class collects the stitches with the ``add`` method:
+
+.. code::
+
+ stitch_builder::add(mstitch, const std::string& parent_id, double along = 1.)
+ stitch_builder::add(mstitch, double along = 1.)
+
+The first stitch will have no parent. If no parent id is specified for a subsequent
+stitch, the last stitch added will be used as parent. The ``along`` parameter
+must lie between zero and one inclusive, and determines the point of attachment
+as a relative position between the parent's proximal and distal points.
+
+A :cpp:type:`stitched_morphology` is constructed from a :cpp:type:`stitch_builder`,
+and provides both the :cpp:type:`morphology` built from the stitches, and methods
+for querying the extent of individual stitches.
+
+.. cpp:class:: stitched_morphology
+
+ .. cpp:function:: stitched_morphology(const stitch_builder&)
+ .. cpp:function:: stitched_morphology(stitch_builder&&)
+
+ Construct from a ``stitch_builder``. Note that constructing from an
+ rvalue is more efficient, as it avoids making a copy of the underlying
+ tree structure.
+
+ .. cpp:function:: arb::morphology morphology() const
+
+ Return the constructed morphology object.
+
+ .. cpp:function:: region stitch(const std::string& id) const
+
+ Return the region expression corresponding to the specified stitch.
-Each segment will inherit the electrical properties of the cell, unless
-otherwise overriden (see below).
+ .. cpp:function:: std::vector<msize_t> segments(const std::string& id) const
+ Return the collection of segments by index comprising the specified stitch.
+
+ .. cpp:function:: label_dict labels(const std::string& prefix="") const
+
+ Provide a :cpp:type:`label_dict` with a region entry for each stitch; if
+ a prefix is provided, this prefix is applied to each segment id to determine
+ the region labels.
+
+Example code, constructing a cable cell from a T-shaped morphology specified
+by two stitches:
+
+.. code::
+
+ using namespace arb;
+
+ mpoint soma0{0, 0, 0, 10};
+ mpoint soma1{20, 0, 0, 10};
+ mpoint dend_end{10, 100, 0, 1};
+
+ stitch_builder builder;
+ builder.add({"soma", soma0, soma1, 1});
+ builder.add({"dend", dend_end, 4}, "soma", 0.5);
+
+ stitched_morphology stitched(std::move(builder));
+ cable_cell cell(stitched.morphology(), stitched.labels());
+
+ cell.paint("soma", "hh");
+
+
+.. _locsets-and-regions:
+
+Identifying sites and subsets of the morphology
+-----------------------------------------------
+
+.. todo::
+
+ TODO: Region and locset documentation is under development.
+
+.. _cable-cell-dynamics:
Cell dynamics
-------------
@@ -62,8 +164,8 @@ which describe biophysical processes. These are processes
that are distributed in space, but whose behaviour is defined purely
by the state of the cell and the process at any given point.
-Cells may also have *point* mechanisms, which are added directly to the
-:cpp:type:`cable_cell` object.
+Cells may also have *point* mechanisms, describing the dynamics
+at post-synaptic sites.
A third type of mechanism, which describes ionic reversal potential
behaviour, can be specified for cells or the whole model via cell parameter
@@ -92,16 +194,19 @@ mechanism name, and mechanism parameter values then set with the
Density mechanisms are associated with a cable cell object with:
-.. cpp:function:: void segment::add_mechanism(mechanism_desc mech)
+.. cpp:function:: void cable_cell::paint(const region&, mechanism_desc)
Point mechanisms, which are associated with connection end points on a
cable cell, are attached to a cell with:
-.. cpp:function:: void cable_cell::add_synapse(mlocation loc, mechanism_desc mech)
+.. cpp:function:: void cable_cell::place(const locset&, mechanism_desc)
-where :cpp:type:`mlocation` is a simple struct holding a segment index
-and a relative position (from 0, proximal, to 1, distal) along that segment:
+.. todo::
+
+ TODO: describe other ``place``-able things: current clamps, gap junction
+ sites, threshold potential measurement point.
+.. _electrical-properties:
Electrical properities and ion values
-------------------------------------
@@ -169,39 +274,39 @@ Global properties
.. cpp:member:: const mechanism_catalogue* catalogue
- All mechanism names refer to mechanism instances in this mechanism catalogue.
- By default, this is set to point to `global_default_catalogue()`, the catalogue
- that contains all mechanisms bundled with Arbor.
+ all mechanism names refer to mechanism instances in this mechanism catalogue.
+ by default, this is set to point to `global_default_catalogue()`, the catalogue
+ that contains all mechanisms bundled with arbor.
- .. cpp:member:: double membrane_voltage_limit_mV
+ .. cpp:member:: double membrane_voltage_limit_mv
- If non-zero, check to see if the membrane voltage ever exceeds this value
- in magnitude during the course of a simulation. If so, throw an exception
+ if non-zero, check to see if the membrane voltage ever exceeds this value
+ in magnitude during the course of a simulation. if so, throw an exception
and abort the simulation.
.. cpp:member:: bool coalesce_synapses
- When synapse dynamics are sufficiently simple, the states of synapses within
- the same discretized element can be combined for better performance. This
+ when synapse dynamics are sufficiently simple, the states of synapses within
+ the same discretized element can be combined for better performance. this
is true by default.
.. cpp:member:: std::unordered_map<std::string, int> ion_species
- Every ion species used by cable cells in the simulation must have an entry in
+ every ion species used by cable cells in the simulation must have an entry in
this map, which takes an ion name to its charge, expressed as a multiple of
- the elementary charge. By default, it is set to include sodium "na" with
+ the elementary charge. by default, it is set to include sodium "na" with
charge 1, calcium "ca" with charge 2, and potassium "k" with charge 1.
.. cpp:member:: cable_cell_parameter_set default_parameters
- The default electrical and physical properties associated with each cable
- cell, unless overridden locally. In the global properties, *every
+ the default electrical and physical properties associated with each cable
+ cell, unless overridden locally. in the global properties, *every
optional field must be given a value*, and every ion must have its default
values set in :cpp:expr:`default_parameters.ion_data`.
.. cpp:function:: add_ion(const std::string& ion_name, int charge, double init_iconc, double init_econc, double init_revpot)
- Convenience function for adding a new ion to the global :cpp:expr:`ion_species`
+ convenience function for adding a new ion to the global :cpp:expr:`ion_species`
table, and setting up its default values in the `ion_data` table.
.. cpp:function:: add_ion(const std::string& ion_name, int charge, double init_iconc, double init_econc, mechanism_desc revpot_mechanism)
@@ -261,6 +366,17 @@ constants.
gprop.default_parameters.reversal_potential_method["ca"] = "nernst1998/ca";
+.. _paint-properties:
+
+Overriding properties locally
+-----------------------------
+
+.. todo::
+
+ TODO: using ``paint`` to specify electrical properties on subsections of
+ the morphology.
+
+
Cable cell probes
-----------------
@@ -291,7 +407,7 @@ The probe metadata passed to the sampler will be a const pointer to:
The type ``cable_probe_point_info`` holds metadata for a single target on a cell:
.. code::
-
+
struct cable_probe_point_info {
// Target number of point process instance on cell.
cell_lid_type target;
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index ee680570..c7dfc382 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -126,8 +126,10 @@ set(unit_sources
test_morph_expr.cpp
test_morph_place.cpp
test_morph_primitives.cpp
+ test_morph_stitch.cpp
test_multi_event_stream.cpp
test_optional.cpp
+ test_ordered_forest.cpp
test_padded.cpp
test_partition.cpp
test_partition_by_constraint.cpp
diff --git a/test/unit/test_morph_embedding.cpp b/test/unit/test_morph_embedding.cpp
index f8de6f4a..59570a6a 100644
--- a/test/unit/test_morph_embedding.cpp
+++ b/test/unit/test_morph_embedding.cpp
@@ -13,28 +13,13 @@
#include "../test/gtest.h"
#include "common.hpp"
#include "common_cells.hpp"
+#include "morph_pred.hpp"
using namespace arb;
using embedding = embed_pwlin;
-::testing::AssertionResult location_eq(const morphology& m, mlocation a, mlocation b) {
- a = canonical(m, a);
- b = canonical(m, b);
-
- if (a.branch!=b.branch) {
- return ::testing::AssertionFailure()
- << "branch ids " << a.branch << " and " << b.branch << " differ";
- }
-
- using FP = testing::internal::FloatingPoint<double>;
- if (FP(a.pos).AlmostEquals(FP(b.pos))) {
- return ::testing::AssertionSuccess();
- }
- else {
- return ::testing::AssertionFailure()
- << "location positions " << a.pos << " and " << b.pos << " differ";
- }
-}
+using testing::mlocation_eq;
+using testing::cable_eq;
TEST(embedding, segments_and_branch_length) {
using pvec = std::vector<msize_t>;
@@ -58,21 +43,20 @@ TEST(embedding, segments_and_branch_length) {
embedding em(m);
auto nloc = 5u;
- EXPECT_EQ(nloc, em.segment_locations().size());
- auto locs = arb::util::make_range(em.branch_segment_locations(0));
+ EXPECT_EQ(nloc, em.segment_ends().size());
+ const auto& locs = em.segment_ends();
EXPECT_EQ(nloc, locs.size());
- EXPECT_TRUE(location_eq(m, locs[0], (loc{0,0})));
- EXPECT_TRUE(location_eq(m, locs[1], (loc{0,0.1})));
- EXPECT_TRUE(location_eq(m, locs[2], (loc{0,0.3})));
- EXPECT_TRUE(location_eq(m, locs[3], (loc{0,0.7})));
- EXPECT_TRUE(location_eq(m, locs[4], (loc{0,1})));
+ EXPECT_TRUE(mlocation_eq(locs[0], (loc{0,0})));
+ EXPECT_TRUE(mlocation_eq(locs[1], (loc{0,0.1})));
+ EXPECT_TRUE(mlocation_eq(locs[2], (loc{0,0.3})));
+ EXPECT_TRUE(mlocation_eq(locs[3], (loc{0,0.7})));
+ EXPECT_TRUE(mlocation_eq(locs[4], (loc{0,1})));
EXPECT_EQ(10., em.branch_length(0));
}
- // Eight samples
+ // Eight samples - point indices:
//
- // sample ids:
// 0
// 1 3
// 2 4
@@ -96,24 +80,26 @@ TEST(embedding, segments_and_branch_length) {
embedding em(m);
- auto locs = arb::util::make_range(em.branch_segment_locations(0));
- EXPECT_TRUE(location_eq(m, locs[0], (loc{0,0})));
- EXPECT_TRUE(location_eq(m, locs[1], (loc{0,0.1})));
- EXPECT_TRUE(location_eq(m, locs[2], (loc{0,1})));
-
- locs = em.branch_segment_locations(1);
- EXPECT_TRUE(location_eq(m, locs[0], (loc{1,0})));
- EXPECT_TRUE(location_eq(m, locs[1], (loc{1,0.1})));
- EXPECT_TRUE(location_eq(m, locs[2], (loc{1,1})));
-
- locs = em.branch_segment_locations(2);
- EXPECT_TRUE(location_eq(m, locs[0], (loc{2,0})));
- EXPECT_TRUE(location_eq(m, locs[1], (loc{2,1})));
-
- locs = em.branch_segment_locations(3);
- EXPECT_TRUE(location_eq(m, locs[0], (loc{3,0})));
- EXPECT_TRUE(location_eq(m, locs[1], (loc{3,0.15})));
- EXPECT_TRUE(location_eq(m, locs[2], (loc{3,1})));
+ const auto& locs = em.segment_ends();
+ EXPECT_TRUE(mlocation_eq(locs[0], (loc{0,0})));
+ EXPECT_TRUE(mlocation_eq(locs[1], (loc{0,0.1})));
+ EXPECT_TRUE(mlocation_eq(locs[2], (loc{0,1})));
+ EXPECT_TRUE(mlocation_eq(locs[3], (loc{1,0})));
+ EXPECT_TRUE(mlocation_eq(locs[4], (loc{1,0.1})));
+ EXPECT_TRUE(mlocation_eq(locs[5], (loc{1,1})));
+ EXPECT_TRUE(mlocation_eq(locs[6], (loc{2,0})));
+ EXPECT_TRUE(mlocation_eq(locs[7], (loc{2,1})));
+ EXPECT_TRUE(mlocation_eq(locs[8], (loc{3,0})));
+ EXPECT_TRUE(mlocation_eq(locs[9], (loc{3,0.15})));
+ EXPECT_TRUE(mlocation_eq(locs[10], (loc{3,1})));
+
+ EXPECT_TRUE(cable_eq(mcable{0, 0. , 0.1 }, em.segment(0)));
+ EXPECT_TRUE(cable_eq(mcable{0, 0.1 , 1. }, em.segment(1)));
+ EXPECT_TRUE(cable_eq(mcable{1, 0. , 0.1 }, em.segment(2)));
+ EXPECT_TRUE(cable_eq(mcable{1, 0.1 , 1. }, em.segment(3)));
+ EXPECT_TRUE(cable_eq(mcable{2, 0. , 1. }, em.segment(4)));
+ EXPECT_TRUE(cable_eq(mcable{3, 0. , 0.15}, em.segment(5)));
+ EXPECT_TRUE(cable_eq(mcable{3, 0.15, 1. }, em.segment(6)));
EXPECT_EQ(100., em.branch_length(0));
EXPECT_EQ(100., em.branch_length(1));
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
index 0b42a541..f3e38cd1 100644
--- a/test/unit/test_morph_expr.cpp
+++ b/test/unit/test_morph_expr.cpp
@@ -353,6 +353,8 @@ TEST(region, thingify_simple_morphologies) {
auto q2 = reg::cable(0, 0.25, 0.5);
auto q3 = reg::cable(0, 0.5, 0.75);
auto q4 = reg::cable(0, 0.75, 1);
+ auto s0 = reg::segment(0);
+ auto s2 = reg::segment(2);
// Concrete cable lists
cl h1_{{0, 0.0, 0.5}};
@@ -360,8 +362,9 @@ TEST(region, thingify_simple_morphologies) {
cl t1_{{0, 0.0, 0.1}, {0, 0.3, 0.7}};
cl t2_{{0, 0.1, 0.3}, {0, 0.7, 1.0}};
cl all_{{0, 0, 1}};
- cl empty_{};
+ EXPECT_TRUE(region_eq(mp, s0, cl{{0, 0, 0.1}}));
+ EXPECT_TRUE(region_eq(mp, s2, cl{{0, 0.3, 0.7}}));
EXPECT_TRUE(region_eq(mp, join(h1, h2), all_));
EXPECT_TRUE(region_eq(mp, intersect(h1, h2), cl{{0, 0.5, 0.5}}));
EXPECT_TRUE(region_eq(mp, t1, t1_));
@@ -389,7 +392,7 @@ TEST(region, thingify_simple_morphologies) {
}
- // Test handling of spherical soma on multi-branch morphologies
+ // A simple Y-shaped morphology.
//
// sample ids:
// 0 |
@@ -417,6 +420,7 @@ TEST(region, thingify_simple_morphologies) {
auto mp = mprovider(morphology(sm));
using ls::location;
+ using reg::segment;
using reg::tagged;
using reg::distal_interval;
using reg::proximal_interval;
@@ -436,6 +440,8 @@ TEST(region, thingify_simple_morphologies) {
auto reg5_ = distal_interval(start1_, 0);
auto reg6_ = proximal_interval(start1_, 0);
+ EXPECT_TRUE(region_eq(mp, segment(3), mcable_list{{2,0,0}}));
+ EXPECT_TRUE(region_eq(mp, segment(4), mcable_list{{2,0,1}}));
EXPECT_TRUE(region_eq(mp, tagged(1), mcable_list{{0,0,1}}));
EXPECT_TRUE(region_eq(mp, tagged(2), mcable_list{{2,0,1}}));
EXPECT_TRUE(region_eq(mp, tagged(3), mcable_list{{1,0,1}}));
@@ -459,20 +465,27 @@ TEST(region, thingify_moderate_morphologies) {
// Test multi-level morphologies.
//
- // Eight samples
- //
- // sample ids:
+ // Eight points:
// 0
// 1 3
// 2 4
// 5 6
// 7
- // tags:
- // 1
- // 3 2
- // 3 2
- // 4 3
- // 3
+ //
+ // Segments [with tag]:
+ // +
+ // ./ \.
+ // 0[3] 2[2]
+ // ./ \.
+ // 1[3] 3[2]
+ // ./ \.
+ // + +
+ // ./ \.
+ // 4[4] 5[3]
+ // ./ \.
+ // + 6[3]
+ // \.
+ // +
{
pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
svec points = {
@@ -492,6 +505,7 @@ TEST(region, thingify_moderate_morphologies) {
auto mp = mprovider(morphology(sm));
using ls::location;
+ using reg::segment;
using reg::tagged;
using reg::distal_interval;
using reg::proximal_interval;
@@ -504,7 +518,6 @@ TEST(region, thingify_moderate_morphologies) {
using reg::cable;
using reg::complete;
- auto soma = tagged(1);
auto axon = tagged(2);
auto dend = tagged(3);
auto apic = tagged(4);
@@ -512,9 +525,6 @@ TEST(region, thingify_moderate_morphologies) {
auto b3 = branch(3);
auto b13 = join(b1, b3);
- cl empty_{};
- cl soma_ = empty_;
-
// Whole branches:
mcable b0_{0,0,1};
mcable b1_{1,0,1};
@@ -529,6 +539,10 @@ TEST(region, thingify_moderate_morphologies) {
EXPECT_TRUE(region_eq(mp, join(dend, apic), cl{b0_,b2_,b3_}));
EXPECT_TRUE(region_eq(mp, join(axon, join(dend, apic)), all_));
+ EXPECT_TRUE(region_eq(mp, tagged(2), join(segment(2), segment(3))));
+ EXPECT_TRUE(region_eq(mp, tagged(3), join(segment(0), segment(1), segment(5), segment(6))));
+ EXPECT_TRUE(region_eq(mp, tagged(4), segment(4)));
+
// Test that intersection correctly generates zero-length cables at
// parent-child interfaces.
EXPECT_TRUE(region_eq(mp, intersect(apic, dend), cl{}));
diff --git a/test/unit/test_morph_stitch.cpp b/test/unit/test_morph_stitch.cpp
new file mode 100644
index 00000000..1712334b
--- /dev/null
+++ b/test/unit/test_morph_stitch.cpp
@@ -0,0 +1,236 @@
+#include <unordered_map>
+#include <string>
+#include <vector>
+
+#include <arbor/morph/morphexcept.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/mprovider.hpp>
+#include <arbor/morph/place_pwlin.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/stitch.hpp>
+
+#include "../test/gtest.h"
+#include "morph_pred.hpp"
+
+using namespace arb;
+using testing::region_eq;
+
+TEST(morph, stitch_none_or_one) {
+ stitch_builder B;
+
+ stitched_morphology sm0(B);
+ EXPECT_TRUE(sm0.morphology().empty());
+
+ mpoint p1{1, 2, 3, 0.5}, p2{2, 4, 5, 1.};
+ B.add({"first", p1, p2, 3});
+
+ stitched_morphology sm1(std::move(B));
+ morphology m1 = sm1.morphology();
+
+ msegment seg0 = m1.branch_segments(0).front();
+ EXPECT_EQ(3, seg0.tag);
+ EXPECT_EQ(p1, seg0.prox);
+ EXPECT_EQ(p2, seg0.dist);
+
+ mprovider p(m1, sm1.labels("stitch:"));
+ EXPECT_TRUE(region_eq(p, "stitch:first", reg::segment(0)));
+}
+
+TEST(morph, stitch_two) {
+ {
+ // p1 ===== p2 ===== p3
+
+ mpoint p1{1, 2, 3, 0.5}, p2{2, 4, 5, 1.}, p3{3, 6, 7, 1.5};
+ stitch_builder B;
+
+ B.add({"0", p1, p2, 0})
+ .add({"1", p3, 1}, "0");
+
+ stitched_morphology sm(std::move(B));
+ morphology m = sm.morphology();
+
+ ASSERT_EQ(1u, m.num_branches());
+ ASSERT_EQ(2u, m.branch_segments(0).size());
+
+ msegment seg0 = m.branch_segments(0)[0];
+ msegment seg1 = m.branch_segments(0)[1];
+
+ EXPECT_EQ(p1, seg0.prox);
+ EXPECT_EQ(p2, seg0.dist);
+ EXPECT_EQ(p2, seg1.prox);
+ EXPECT_EQ(p3, seg1.dist);
+
+ mprovider p(m, sm.labels("stitch:"));
+ EXPECT_TRUE(region_eq(p, "stitch:0", reg::segment(0)));
+ EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(1)));
+ }
+ {
+ // p1 ===== p2
+ // \.
+ // \.
+ // p3
+
+ mpoint p1{1, 2, 3, 0.5}, p2{2, 4, 5, 1.}, p3{3, 6, 7, 1.5};
+ stitch_builder B;
+
+ B.add({"0", p1, p2, 0})
+ .add({"1", p3, 1}, "0", 0);
+
+ stitched_morphology sm(std::move(B));
+ morphology m = sm.morphology();
+
+ ASSERT_EQ(2u, m.num_branches());
+ ASSERT_EQ(1u, m.branch_segments(0).size());
+ ASSERT_EQ(1u, m.branch_segments(1).size());
+
+ EXPECT_EQ(mnpos, m.branch_parent(0));
+ EXPECT_EQ(mnpos, m.branch_parent(1));
+
+ msegment seg0 = m.branch_segments(0)[0];
+ msegment seg1 = m.branch_segments(1)[0];
+
+ EXPECT_EQ(p1, seg0.prox);
+ EXPECT_EQ(p1, seg1.prox);
+
+ mprovider p(m, sm.labels("stitch:"));
+ // Branch ordering is arbitrary, so check both possibilities:
+ if (seg0.dist == p2) {
+ EXPECT_TRUE(region_eq(p, "stitch:0", reg::segment(0)));
+ EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(1)));
+ }
+ else {
+ ASSERT_EQ(p2, seg1.dist);
+ EXPECT_TRUE(region_eq(p, "stitch:0", reg::segment(1)));
+ EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(0)));
+ }
+ }
+ {
+ // p1 ==x== p2
+ // \.
+ // \.
+ // p3
+
+ mpoint p1{1, 2, 3, 0.5}, p2{2, 4, 5, 1.}, p3{3, 6, 7, 1.5};
+ stitch_builder B;
+
+ B.add({"0", p1, p2, 0})
+ .add({"1", p3, 1}, "0", 0.5);
+
+ stitched_morphology sm(std::move(B));
+ morphology m = sm.morphology();
+
+ ASSERT_EQ(3u, m.num_branches());
+ ASSERT_EQ(1u, m.branch_segments(0).size());
+ ASSERT_EQ(1u, m.branch_segments(1).size());
+ ASSERT_EQ(1u, m.branch_segments(2).size());
+
+ msegment seg0 = m.branch_segments(0)[0];
+ msegment seg1 = m.branch_segments(1)[0];
+ msegment seg2 = m.branch_segments(2)[0];
+
+ EXPECT_EQ(p1, seg0.prox);
+ mpoint x = lerp(p1, p2, 0.5);
+ EXPECT_EQ(x, seg0.dist);
+
+ EXPECT_EQ(x, seg1.prox);
+ EXPECT_EQ(x, seg2.prox);
+
+ mprovider p(m, sm.labels("stitch:"));
+ // Branch ordering is arbitrary, so check both possibilities:
+ if (seg2.dist == p2) {
+ EXPECT_TRUE(region_eq(p, "stitch:0", join(reg::segment(0), reg::segment(2))));
+ EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(1)));
+ }
+ else {
+ ASSERT_EQ(p2, seg1.dist);
+ EXPECT_TRUE(region_eq(p, "stitch:0", join(reg::segment(0), reg::segment(1))));
+ EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(2)));
+ }
+ }
+}
+
+TEST(morph, stitch_errors) {
+ mpoint p1{1, 2, 3, 0.5}, p2{2, 4, 5, 1.}, p3{3, 6, 7, 1.5};
+ stitch_builder B;
+
+ B.add({"0", p1, p2, 0});
+ ASSERT_THROW(B.add({"0", p3, 0}, "0", 0.5), duplicate_stitch_id);
+ ASSERT_THROW(B.add({"1", p3, 0}, "x", 0.5), no_such_stitch);
+ ASSERT_THROW(B.add({"1", p3, 0}, "0", 1.5), invalid_stitch_position);
+
+ stitch_builder C;
+ ASSERT_THROW(C.add({"0", p1, 0}), missing_stitch_start);
+}
+
+TEST(morph, stitch_complex) {
+ // p[8]
+ // |
+ // p[0]---x----x----x---p[1]---x---p[2]
+ // | | | |
+ // p[3] p[4] p[5]--p[7] p[6]
+
+ mpoint p[] = {
+ {0, 1, 0, 1.},
+ {4, 1, 0, 1.},
+ {6, 1, 0, 1.},
+ {1, 0, 0, 1.},
+ {2, 0, 0, 1.},
+ {3, 0, 0, 1.},
+ {5, 0, 0, 1.},
+ {4, 0, 0, 1.},
+ {3, 2, 0, 1.}
+ };
+
+ stitch_builder B;
+
+ // (labels chosen to reflect distal point)
+ B.add({"1", p[0], p[1]})
+ .add({"2", p[2]}, "1")
+ .add({"3", p[3]}, "1", 0.25)
+ .add({"4", p[4]}, "1", 0.50)
+ .add({"5", p[5]}, "1", 0.75)
+ .add({"6", p[6]}, "2", 0.50)
+ .add({"7", p[7]}, "5")
+ .add({"8", p[8]}, "1", 0.75);
+
+ stitched_morphology sm(std::move(B));
+ morphology m = sm.morphology();
+ mprovider P(m, sm.labels());
+
+ EXPECT_EQ(10u, m.num_branches());
+
+ EXPECT_EQ(4u, sm.segments("1").size());
+ EXPECT_EQ(2u, sm.segments("2").size());
+ EXPECT_EQ(1u, sm.segments("3").size());
+ EXPECT_EQ(1u, sm.segments("4").size());
+ EXPECT_EQ(1u, sm.segments("5").size());
+ EXPECT_EQ(1u, sm.segments("6").size());
+ EXPECT_EQ(1u, sm.segments("7").size());
+ EXPECT_EQ(1u, sm.segments("8").size());
+
+ auto region_prox_eq = [&P, place = place_pwlin(m)](region r, mpoint p) {
+ mlocation_list ls = thingify(ls::most_proximal(r), P);
+ if (ls.empty()) {
+ return ::testing::AssertionFailure() << "region " << r << " is empty";
+ }
+ else if (ls.size()>1u) {
+ return ::testing::AssertionFailure() << "region " << r << " has multiple proximal points";
+ }
+
+ mpoint prox = place.at(ls.front());
+ if (prox!=p) {
+ return ::testing::AssertionFailure() << "region " << r << " proximal point " << prox << " is not equal to " << p;
+ }
+
+ return ::testing::AssertionSuccess();
+ };
+
+ EXPECT_TRUE(region_prox_eq(sm.stitch("2"), p[1]));
+ EXPECT_TRUE(region_prox_eq(sm.stitch("3"), mpoint{1, 1, 0, 1.}));
+ EXPECT_TRUE(region_prox_eq(sm.stitch("4"), mpoint{2, 1, 0, 1.}));
+ EXPECT_TRUE(region_prox_eq(sm.stitch("5"), mpoint{3, 1, 0, 1.}));
+ EXPECT_TRUE(region_prox_eq(sm.stitch("6"), mpoint{5, 1, 0, 1.}));
+ EXPECT_TRUE(region_prox_eq(sm.stitch("7"), p[5]));
+ EXPECT_TRUE(region_prox_eq(sm.stitch("8"), mpoint{3, 1, 0, 1.}));
+}
+
diff --git a/test/unit/test_ordered_forest.cpp b/test/unit/test_ordered_forest.cpp
new file mode 100644
index 00000000..1ca5aea7
--- /dev/null
+++ b/test/unit/test_ordered_forest.cpp
@@ -0,0 +1,513 @@
+#include <cstddef>
+#include <memory>
+#include <type_traits>
+#include <vector>
+
+#include "../gtest.h"
+#include "util/ordered_forest.hpp"
+
+using arb::util::ordered_forest;
+
+namespace {
+template <typename T>
+struct simple_allocator {
+ using value_type = T;
+
+ simple_allocator():
+ n_alloc_(new std::size_t()),
+ n_dealloc_(new std::size_t())
+ {}
+
+ simple_allocator(const simple_allocator&) noexcept = default;
+
+ template <typename U>
+ simple_allocator(const simple_allocator<U>& other) noexcept {
+ n_alloc_ = other.n_alloc_;
+ n_dealloc_ = other.n_dealloc_;
+ }
+
+ T* allocate(std::size_t n) {
+ auto p = new T[n];
+ ++*n_alloc_;
+ return p;
+ }
+
+ void deallocate(T* p, std::size_t) {
+ delete [] p;
+ ++*n_dealloc_;
+ }
+
+ bool operator==(const simple_allocator& other) const {
+ return other.n_alloc_ == n_alloc_ && other.n_dealloc_ == n_dealloc_;
+ }
+
+ bool operator!=(const simple_allocator& other) const {
+ return !(*this==other);
+ }
+
+ using propagate_on_container_copy_assignment = std::false_type;
+ using propagate_on_container_move_assignment = std::false_type;
+ using propagate_on_container_swap = std::true_type;
+
+ void reset_counts() {
+ *n_alloc_ = 0;
+ *n_dealloc_ = 0;
+ }
+
+ std::size_t n_alloc() const { return *n_alloc_; }
+ std::size_t n_dealloc() const { return *n_dealloc_; }
+
+ std::shared_ptr<std::size_t> n_alloc_, n_dealloc_;
+};
+
+template <typename T, typename A>
+::testing::AssertionResult forest_invariant(const ordered_forest<T, A>& f) {
+ // check parent-child confluence:
+
+ for (auto i = f.root_begin(); i!=f.root_end(); ++i) {
+ if (i.parent()) {
+ return ::testing::AssertionFailure() << "root node " << *i << " has parent " << *i.parent();
+ }
+ }
+
+ for (auto i = f.begin(); i!=f.end(); ++i) {
+ auto cb = f.child_begin(i);
+ auto ce = f.child_end(i);
+
+ for (auto j = cb; j!=ce; ++j) {
+ if (j.parent() != i) {
+ auto failure = ::testing::AssertionFailure() << "child node " << *j << " of " << *i << " has parent ";
+ return j.parent()? (failure << *j.parent()): (failure << "(null)");
+ }
+ }
+ }
+ return ::testing::AssertionSuccess();
+}
+
+} // anonymous namespace
+
+TEST(ordered_forest, empty) {
+ ordered_forest<int> f1;
+ EXPECT_EQ(0u, f1.size());
+ EXPECT_TRUE(f1.begin() == f1.end());
+
+ simple_allocator<int> alloc;
+ ASSERT_EQ(0u, alloc.n_alloc());
+ ASSERT_EQ(0u, alloc.n_dealloc());
+
+ ordered_forest<int, simple_allocator<int>> f2(alloc);
+ EXPECT_EQ(0u, alloc.n_alloc());
+ EXPECT_EQ(0u, alloc.n_dealloc());
+}
+
+TEST(ordered_forest, push) {
+ simple_allocator<int> alloc;
+
+ {
+ ordered_forest<int, simple_allocator<int>> f(alloc);
+
+ f.push_front(3);
+ auto i2 = f.push_front(2);
+ f.push_child(i2, 5);
+ f.push_child(i2, 4);
+ f.push_front(1);
+
+ ASSERT_TRUE(forest_invariant(f));
+
+ ASSERT_EQ(5u, f.size());
+ EXPECT_EQ(10u, alloc.n_alloc()); // five nodes, five items.
+
+ auto i = f.begin();
+ ASSERT_TRUE(i);
+ EXPECT_EQ(1, *i);
+ EXPECT_FALSE(i.child());
+
+ i = i.next();
+ ASSERT_TRUE(i);
+ auto j = i.child();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(4, *j);
+ j = j.next();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(5, *j);
+ EXPECT_FALSE(j.next());
+
+ i = i.next();
+ ASSERT_TRUE(i);
+ EXPECT_EQ(3, *i);
+ EXPECT_FALSE(i.child());
+ EXPECT_FALSE(i.next());
+ }
+
+ EXPECT_EQ(alloc.n_dealloc(), alloc.n_alloc());
+}
+
+TEST(ordered_forest, insert) {
+ simple_allocator<int> alloc;
+
+ {
+ ordered_forest<int, simple_allocator<int>> f(alloc);
+
+ auto r = f.push_front(1);
+ f.insert_after(r, 3);
+ r = f.insert_after(r, 2);
+ auto c = f.push_child(r, 4);
+ f.insert_after(c, 6);
+ f.insert_after(c, 5);
+
+ ASSERT_TRUE(forest_invariant(f));
+
+ ASSERT_EQ(6u, f.size());
+ EXPECT_EQ(12u, alloc.n_alloc()); // six nodes, six items.
+
+ auto i = f.begin();
+ ASSERT_TRUE(i);
+ EXPECT_EQ(1, *i);
+ EXPECT_FALSE(i.child());
+
+ i = i.next();
+ ASSERT_TRUE(i);
+ auto j = i.child();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(4, *j);
+ j = j.next();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(5, *j);
+ j = j.next();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(6, *j);
+ EXPECT_FALSE(j.next());
+
+ i = i.next();
+ ASSERT_TRUE(i);
+ EXPECT_EQ(3, *i);
+ EXPECT_FALSE(i.child());
+ EXPECT_FALSE(i.next());
+ }
+
+ EXPECT_EQ(alloc.n_dealloc(), alloc.n_alloc());
+}
+
+TEST(ordered_forest, initializer_list) {
+ ordered_forest<int> f = {1, {2, {4, 5, 6}}, 3};
+ EXPECT_EQ(6u, f.size());
+
+ ASSERT_TRUE(forest_invariant(f));
+
+ auto i = f.begin();
+ ASSERT_TRUE(i);
+ EXPECT_EQ(1, *i);
+ EXPECT_FALSE(i.child());
+
+ i = i.next();
+ ASSERT_TRUE(i);
+ auto j = i.child();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(4, *j);
+ j = j.next();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(5, *j);
+ j = j.next();
+ ASSERT_TRUE(j);
+ EXPECT_EQ(6, *j);
+ EXPECT_FALSE(j.next());
+
+ i = i.next();
+ ASSERT_TRUE(i);
+ EXPECT_EQ(3, *i);
+ EXPECT_FALSE(i.child());
+ EXPECT_FALSE(i.next());
+}
+
+TEST(ordered_forest, equality) {
+ using of = ordered_forest<int>;
+
+ EXPECT_EQ(of{}, of{});
+ EXPECT_NE(of{1}, of{});
+ EXPECT_NE(of{}, of{1});
+
+ EXPECT_EQ((of{1, 2, 3}), (of{1, 2, 3}));
+ EXPECT_NE((of{1, 2, 3}), (of{1, {2, {3}}}));
+ EXPECT_NE((of{{1, {2, 3}}}), (of{1, 2, 3}));
+
+ struct always_eq {
+ int n_;
+ always_eq(int n): n_(n) {}
+ bool operator==(const always_eq&) const { return true; }
+ bool operator!=(const always_eq&) const { return false; }
+ };
+
+ ordered_forest<always_eq> f1 = {{1, {2, 3}}, {4, {5, {6, {7}}, 8}}, 9};
+ ordered_forest<always_eq> f2 = {{3, {1, 0}}, {2, {8, {6, {4}}, 7}}, 5};
+ EXPECT_EQ(f2, f1);
+
+ ordered_forest<always_eq> f3 = {{3, {{1, {0}}}}, {2, {8, {6, {4}}, 7}}, 5};
+ EXPECT_NE(f1, f3);
+}
+
+TEST(ordered_forest, iteration) {
+ using ivector = std::vector<int>;
+
+ ordered_forest<int> f = {{1, {2, 3}}, {4, {5, {6, {7}}, 8}}, 9};
+ const ordered_forest<int>& cf = f;
+
+ ivector pre{f.preorder_begin(), f.preorder_end()};
+ EXPECT_EQ((ivector{1, 2, 3, 4, 5, 6, 7, 8, 9}), pre);
+
+ ivector cpre{cf.preorder_begin(), cf.preorder_end()};
+ EXPECT_EQ((ivector{1, 2, 3, 4, 5, 6, 7, 8, 9}), cpre);
+
+ ivector pre_default{f.begin(), f.end()};
+ EXPECT_EQ((ivector{1, 2, 3, 4, 5, 6, 7, 8, 9}), pre_default);
+
+ ivector cpre_default{cf.begin(), cf.end()};
+ EXPECT_EQ((ivector{1, 2, 3, 4, 5, 6, 7, 8, 9}), cpre_default);
+
+ ivector pre_cdefault{f.cbegin(), f.cend()};
+ EXPECT_EQ((ivector{1, 2, 3, 4, 5, 6, 7, 8, 9}), pre_cdefault);
+
+ ivector root{f.root_begin(), f.root_end()};
+ EXPECT_EQ((ivector{1, 4, 9}), root);
+
+ ivector croot{cf.root_begin(), cf.root_end()};
+ EXPECT_EQ((ivector{1, 4, 9}), croot);
+
+ ivector post{f.postorder_begin(), f.postorder_end()};
+ EXPECT_EQ((ivector{2, 3, 1, 5, 7, 6, 8, 4, 9}), post);
+
+ ivector cpost{cf.postorder_begin(), cf.postorder_end()};
+ EXPECT_EQ((ivector{2, 3, 1, 5, 7, 6, 8, 4, 9}), cpost);
+
+ auto four = std::find(f.begin(), f.end(), 4);
+ ivector child_four{f.child_begin(four), f.child_end(four)};
+ EXPECT_EQ((ivector{5, 6, 8}), child_four);
+
+ ivector cchild_four{cf.child_begin(four), cf.child_end(four)};
+ EXPECT_EQ((ivector{5, 6, 8}), cchild_four);
+
+ using preorder_iterator = ordered_forest<int>::preorder_iterator;
+ using postorder_iterator = ordered_forest<int>::postorder_iterator;
+
+ ivector pre_four{preorder_iterator(four), preorder_iterator(four.next())};
+ EXPECT_EQ((ivector{4, 5, 6, 7, 8}), pre_four);
+
+ auto seven = std::find(f.begin(), f.end(), 7);
+ auto nine = std::find(f.begin(), f.end(), 9);
+ ivector post_seven_nine{postorder_iterator(seven), postorder_iterator(nine)};
+ EXPECT_EQ((ivector{7, 6, 8, 4}), post_seven_nine);
+}
+
+TEST(ordered_forest, copy_move) {
+ simple_allocator<int> alloc;
+
+ using ivector = std::vector<int>;
+ using of = ordered_forest<int, simple_allocator<int>>;
+
+ of f1(alloc);
+ {
+ of f({{1, {2, 3}}, {4, {5, {6, {7}}, 8}}, 9}, alloc);
+ EXPECT_EQ(18u, alloc.n_alloc());
+ ASSERT_TRUE(forest_invariant(f));
+
+ f1 = f;
+ ASSERT_TRUE(forest_invariant(f1));
+ EXPECT_EQ(36u, alloc.n_alloc());
+ EXPECT_FALSE(f.empty());
+
+ of f2 = std::move(f);
+ ASSERT_TRUE(forest_invariant(f2));
+ EXPECT_EQ(36u, alloc.n_alloc());
+ EXPECT_TRUE(f.empty());
+
+ ivector elems2{f2.begin(), f2.end()};
+ EXPECT_EQ((ivector{1, 2, 3, 4, 5, 6, 7, 8, 9}), elems2);
+ }
+
+ EXPECT_EQ(36u, alloc.n_alloc());
+ EXPECT_EQ(18u, alloc.n_dealloc());
+
+ ivector elems1{f1.begin(), f1.end()};
+ EXPECT_EQ((ivector{1, 2, 3, 4, 5, 6, 7, 8, 9}), elems1);
+
+ // With a different != allocator object, require move assignment
+ // to perform copy.
+
+ simple_allocator<int> other_alloc;
+ ASSERT_NE(alloc, other_alloc);
+ ASSERT_FALSE(std::allocator_traits<simple_allocator<int>>::propagate_on_container_move_assignment::value);
+
+ of f3(other_alloc);
+
+ f3 = std::move(f1);
+ ASSERT_TRUE(forest_invariant(f3));
+ EXPECT_FALSE(f1.empty());
+
+ EXPECT_EQ(36u, alloc.n_alloc());
+ EXPECT_EQ(18u, alloc.n_dealloc());
+ EXPECT_EQ(18u, other_alloc.n_alloc());
+ EXPECT_EQ(0u, other_alloc.n_dealloc());
+}
+
+TEST(ordered_forest, erase) {
+ simple_allocator<int> alloc;
+ using of = ordered_forest<int, simple_allocator<int>>;
+
+ of f({1, 2, {3, {4, {5, {6, 7}}, 8}}, 9}, alloc);
+ ASSERT_TRUE(forest_invariant(f));
+ EXPECT_EQ(18u, alloc.n_alloc());
+
+ auto two = std::find(f.begin(), f.end(), 2);
+ f.erase_after(two);
+
+ ASSERT_TRUE(forest_invariant(f));
+ EXPECT_EQ((of{1, 2, 4, {5, {6, 7}}, 8, 9}), f);
+ EXPECT_EQ(2u, alloc.n_dealloc());
+
+ auto five = std::find(f.begin(), f.end(), 5);
+ f.erase_child(five);
+
+ ASSERT_TRUE(forest_invariant(f));
+ EXPECT_EQ((of{1, 2, 4, {5, {7}}, 8, 9}), f);
+ EXPECT_EQ(4u, alloc.n_dealloc());
+
+ auto eight = std::find(f.begin(), f.end(), 8);
+ ASSERT_THROW(f.erase_child(eight), std::invalid_argument);
+
+ auto seven = std::find(f.begin(), f.end(), 7);
+ ASSERT_THROW(f.erase_after(seven), std::invalid_argument);
+
+ f.erase_front();
+ ASSERT_TRUE(forest_invariant(f));
+ EXPECT_EQ((of{2, 4, {5, {7}}, 8, 9}), f);
+ EXPECT_EQ(6u, alloc.n_dealloc());
+
+ of empty;
+ ASSERT_THROW(empty.erase_front(), std::invalid_argument);
+}
+
+TEST(ordered_forest, prune) {
+ simple_allocator<int> alloc;
+ using of = ordered_forest<int, simple_allocator<int>>;
+
+ of f({1, 2, {3, {4, {5, {6, 7}}, 8}}, 9}, alloc);
+ ASSERT_TRUE(forest_invariant(f));
+ EXPECT_EQ(18u, alloc.n_alloc());
+ EXPECT_EQ(0u, alloc.n_dealloc());
+
+ of p1 = f.prune_front();
+ ASSERT_TRUE(forest_invariant(f));
+ ASSERT_TRUE(forest_invariant(p1));
+ EXPECT_EQ((of{2, {3, {4, {5, {6, 7}}, 8}}, 9}), f);
+ EXPECT_EQ((of{1}), p1);
+
+ of p2 = f.prune_after(std::find(f.begin(), f.end(), 4));
+ ASSERT_TRUE(forest_invariant(f));
+ ASSERT_TRUE(forest_invariant(p2));
+ EXPECT_EQ((of{2, {3, {4, 8}}, 9}), f);
+ EXPECT_EQ((of{{5, {6, 7}}}), p2);
+
+ of p3 = f.prune_child(std::find(f.begin(), f.end(), 3));
+ ASSERT_TRUE(forest_invariant(f));
+ ASSERT_TRUE(forest_invariant(p3));
+ EXPECT_EQ((of{2, {3, {8}}, 9}), f);
+ EXPECT_EQ((of{4}), p3);
+
+ EXPECT_EQ(0u, alloc.n_dealloc());
+ EXPECT_EQ(f.get_allocator(), p1.get_allocator());
+ EXPECT_EQ(f.get_allocator(), p2.get_allocator());
+ EXPECT_EQ(f.get_allocator(), p3.get_allocator());
+
+ of empty;
+ ASSERT_THROW(empty.erase_child(empty.begin()), std::invalid_argument);
+ ASSERT_THROW(empty.erase_after(empty.begin()), std::invalid_argument);
+ ASSERT_THROW(empty.erase_front(), std::invalid_argument);
+
+ of unit{1};
+ ASSERT_THROW(unit.erase_child(unit.begin()), std::invalid_argument);
+ ASSERT_THROW(unit.erase_after(unit.begin()), std::invalid_argument);
+}
+
+TEST(ordered_forest, graft) {
+ using of = ordered_forest<int, simple_allocator<int>>;
+
+ of f1{1, {2, {3, 4}}, 5};
+ auto j = f1.graft_after(f1.begin(), of{6, {7, {8}}});
+
+ ASSERT_TRUE(j);
+ ASSERT_TRUE(forest_invariant(f1));
+ EXPECT_EQ(7, *j);
+ EXPECT_EQ((of{1, 6, {7, {8}}, {2, {3, 4}}, 5}), f1);
+
+ j = f1.graft_child(std::find(f1.begin(), f1.end(), 2), of{9, 10});
+
+ ASSERT_TRUE(j);
+ ASSERT_TRUE(forest_invariant(f1));
+ EXPECT_EQ(10, *j);
+ EXPECT_EQ((of{1, 6, {7, {8}}, {2, {9, 10, 3, 4}}, 5}), f1);
+
+ j = f1.graft_front(of{{11, {12, 13}}});
+
+ ASSERT_TRUE(j);
+ ASSERT_TRUE(forest_invariant(f1));
+ EXPECT_EQ(11, *j);
+ EXPECT_EQ((of{{11, {12, 13}}, 1, 6, {7, {8}}, {2, {9, 10, 3, 4}}, 5}), f1);
+
+ simple_allocator<int> alloc1, alloc2;
+ of f2({1, 2}, alloc1);
+ of f3({3, 4}, alloc1);
+ of f4({5, 6}, alloc2);
+
+ ASSERT_TRUE(forest_invariant(f2));
+ ASSERT_TRUE(forest_invariant(f3));
+ ASSERT_TRUE(forest_invariant(f4));
+
+ EXPECT_EQ(8u, alloc1.n_alloc());
+ EXPECT_EQ(0u, alloc1.n_dealloc());
+ EXPECT_EQ(4u, alloc2.n_alloc());
+
+ f2.graft_front(std::move(f3));
+ ASSERT_TRUE(forest_invariant(f2));
+ ASSERT_TRUE(forest_invariant(f3));
+ EXPECT_EQ(8u, alloc1.n_alloc());
+ EXPECT_EQ(0u, alloc1.n_dealloc());
+
+ f2.graft_front(std::move(f4));
+ ASSERT_TRUE(forest_invariant(f2));
+ ASSERT_TRUE(forest_invariant(f4));
+ EXPECT_EQ(12u, alloc1.n_alloc());
+ EXPECT_EQ(0u, alloc1.n_dealloc());
+ EXPECT_EQ(4u, alloc2.n_dealloc());
+
+ EXPECT_EQ((of{5, 6, 3, 4, 1, 2}), f2);
+}
+
+TEST(ordered_forest, swap) {
+ simple_allocator<int> alloc1, alloc2;
+ using of = ordered_forest<int, simple_allocator<int>>;
+
+ of a({1, {2, {3, 4}}, 5}, alloc1);
+ of b({6, 7, 8, 9}, alloc2);
+
+ of a_copy(a), b_copy(b);
+
+ ASSERT_TRUE(forest_invariant(a));
+ ASSERT_TRUE(forest_invariant(b));
+ ASSERT_TRUE(forest_invariant(a_copy));
+ ASSERT_TRUE(forest_invariant(b_copy));
+
+ ASSERT_EQ(alloc1, a.get_allocator());
+ ASSERT_EQ(alloc2, b.get_allocator());
+ ASSERT_EQ(alloc1, a_copy.get_allocator());
+ ASSERT_EQ(alloc2, b_copy.get_allocator());
+ ASSERT_EQ(a_copy, a);
+ ASSERT_EQ(b_copy, b);
+
+ swap(a, b);
+
+ ASSERT_TRUE(forest_invariant(a));
+ ASSERT_TRUE(forest_invariant(b));
+ EXPECT_EQ(alloc2, a.get_allocator());
+ EXPECT_EQ(alloc1, b.get_allocator());
+ EXPECT_EQ(b_copy, a);
+ EXPECT_EQ(a_copy, b);
+}
--
GitLab