diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 2257032d3f4754a461d15d9c837a9ea95c48c05c..24e73acd1a904f874c61252fe01c46735c95871d 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -28,11 +28,13 @@ set(arbor_sources
mechcat.cpp
memory/cuda_wrappers.cpp
memory/util.cpp
- morph/em_morphology.cpp
+ morph/embed_pwlin.cpp
morph/label_dict.cpp
morph/locset.cpp
morph/mbranch.cpp
+ morph/morphexcept.cpp
morph/morphology.cpp
+ morph/mprovider.cpp
morph/primitives.cpp
morph/region.cpp
morph/sample_tree.cpp
diff --git a/arbor/algorithms.hpp b/arbor/algorithms.hpp
index ec4b2dc45df025c77ae8ab3d67a555343ea48585..722edd0724c6878c1ecf4c6c7afc97b833259d72 100644
--- a/arbor/algorithms.hpp
+++ b/arbor/algorithms.hpp
@@ -26,19 +26,11 @@
namespace arb {
namespace algorithms {
-template <typename C>
-typename util::sequence_traits<C>::value_type
-sum(C const& c)
-{
- using value_type = typename util::sequence_traits<C>::value_type;
- return std::accumulate(std::cbegin(c), std::cend(c), value_type{0});
-}
-
template <typename C>
typename util::sequence_traits<C>::value_type
mean(C const& c)
{
- return sum(c)/util::size(c);
+ return util::sum(c)/util::size(c);
}
// returns the prefix sum of c in the form `[0, c[0], c[0]+c[1], ..., sum(c)]`.
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index d90094e834713776401cc38a8cc46c0599842a23..6dea2d3991bb90191d8470dd16b64cc46dbe3d02 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -1,18 +1,23 @@
#include <sstream>
+#include <unordered_map>
#include <vector>
#include <arbor/cable_cell.hpp>
#include <arbor/morph/label_dict.hpp>
#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/mprovider.hpp>
#include <arbor/segment.hpp>
-#include "morph/em_morphology.hpp"
+#include "util/piecewise.hpp"
#include "util/rangeutil.hpp"
#include "util/span.hpp"
#include "util/strprintf.hpp"
namespace arb {
+using region_map = std::unordered_map<std::string, mcable_list>;
+using locset_map = std::unordered_map<std::string, mlocation_list>;
+
using value_type = cable_cell::value_type;
using index_type = cable_cell::index_type;
using size_type = cable_cell::size_type;
@@ -27,13 +32,10 @@ struct cable_cell_impl {
using gap_junction_instance = cable_cell::gap_junction_instance;
using detector_instance = cable_cell::detector_instance;
- using region_map = cable_cell::region_map;
- using locset_map = cable_cell::locset_map;
-
cable_cell_impl(const arb::morphology& m,
const label_dict& dictionary,
bool compartments_from_discretization):
- morph(m)
+ provider(m, dictionary)
{
using point = cable_cell::point_type;
if (!m.num_branches()) {
@@ -97,28 +99,18 @@ struct cable_cell_impl {
segments.back()->as_cable()->set_compartments(ncomp);
}
}
-
- for (auto& r: dictionary.regions()) {
- regions[r.first] = thingify(r.second, morph);
- }
-
- for (auto& l: dictionary.locsets()) {
- locations[l.first] = thingify(l.second, morph);
- }
}
cable_cell_impl(): cable_cell_impl({},{},false) {}
- cable_cell_impl(const cable_cell_impl& other) {
- parents = other.parents;
- stimuli = other.stimuli;
- synapses = other.synapses;
- gap_junction_sites = other.gap_junction_sites;
- spike_detectors = other.spike_detectors;
- regions = other.regions;
- morph = other.morph;
- locations = other.locations;
-
+ cable_cell_impl(const cable_cell_impl& other):
+ parents(other.parents),
+ stimuli(other.stimuli),
+ synapses(other.synapses),
+ gap_junction_sites(other.gap_junction_sites),
+ spike_detectors(other.spike_detectors),
+ provider(other.provider)
+ {
// unique_ptr's cannot be copy constructed, do a manual assignment
segments.reserve(other.segments.size());
for (const auto& s: other.segments) {
@@ -146,14 +138,8 @@ struct cable_cell_impl {
// the sensors
std::vector<detector_instance> spike_detectors;
- // Named regions
- region_map regions;
-
- // Named location sets
- locset_map locations;
-
- // Underlying embedded morphology
- em_morphology morph;
+ // Embedded morphology and labelled region/locset lookup.
+ mprovider provider;
template <typename Desc, typename T>
lid_range place(const mlocation_list& locs, const Desc& desc, std::vector<T>& list) {
@@ -167,13 +153,8 @@ struct cable_cell_impl {
}
template <typename Desc, typename T>
- lid_range place(const std::string& target, const Desc& desc, std::vector<T>& list) {
- const auto first = list.size();
-
- const auto it = locations.find(target);
- if (it==locations.end()) return lid_range(first, first);
-
- return place(it->second, desc, list);
+ lid_range place(const locset& locs, const Desc& desc, std::vector<T>& list) {
+ return place(thingify(locs, provider), desc, list);
}
lid_range place_gj(const mlocation_list& locs) {
@@ -184,13 +165,8 @@ struct cable_cell_impl {
return lid_range(first, gap_junction_sites.size());
}
- lid_range place_gj(const std::string& target) {
- const auto first = gap_junction_sites.size();
-
- const auto it = locations.find(target);
- if (it==locations.end()) return lid_range(first, first);
-
- return place_gj(it->second);
+ lid_range place_gj(const locset& locs) {
+ return place_gj(thingify(locs, provider));
}
void assert_valid_segment(index_type i) const {
@@ -203,16 +179,6 @@ struct cable_cell_impl {
return test_invariants(loc) && loc.branch<segments.size();
}
- template <typename F>
- void paint(const std::string& target, F&& f) {
- auto it = regions.find(target);
-
- // Nothing to do if there are no regions that match.
- if (it==regions.end()) return;
-
- paint(it->second, std::forward<F>(f));
- }
-
template <typename F>
void paint(const mcable_list& cables, F&& f) {
for (auto c: cables) {
@@ -224,6 +190,11 @@ struct cable_cell_impl {
f(segments[c.branch]);
}
}
+
+ template <typename F>
+ void paint(const region& reg, F&& f) {
+ paint(thingify(reg, provider), std::forward<F>(f));
+ }
};
using impl_ptr = std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)>;
@@ -294,33 +265,32 @@ const std::vector<cable_cell::stimulus_instance>& cable_cell::stimuli() const {
return impl_->stimuli;
}
-const em_morphology* cable_cell::morphology() const {
- return &(impl_->morph);
+const concrete_embedding& cable_cell::embedding() const {
+ return impl_->provider.embedding();
+}
+
+const arb::morphology& cable_cell::morphology() const {
+ return impl_->provider.morphology();
}
+const mprovider& cable_cell::provider() const {
+ return impl_->provider;
+}
+
+
//
// Painters.
//
// Implementation of user API for painting density channel and electrical properties on cells.
//
-void cable_cell::paint(const std::string& target, mechanism_desc desc) {
- impl_->paint(target,
- [&desc](segment_ptr& s){return s->add_mechanism(desc);});
-}
-
-void cable_cell::paint(const std::string& target, cable_cell_local_parameter_set params) {
- impl_->paint(target,
- [¶ms](segment_ptr& s){return s->parameters = params;});
-}
-
void cable_cell::paint(const region& target, mechanism_desc desc) {
- impl_->paint(thingify(target, impl_->morph),
+ impl_->paint(target,
[&desc](segment_ptr& s){return s->add_mechanism(desc);});
}
void cable_cell::paint(const region& target, cable_cell_local_parameter_set params) {
- impl_->paint(thingify(target, impl_->morph),
+ impl_->paint(target,
[¶ms](segment_ptr& s){return s->parameters = params;});
}
@@ -335,47 +305,32 @@ void cable_cell::paint(const region& target, cable_cell_local_parameter_set para
// Synapses.
//
-lid_range cable_cell::place(const std::string& target, const mechanism_desc& desc) {
- return impl_->place(target, desc, impl_->synapses);
-}
-
lid_range cable_cell::place(const locset& ls, const mechanism_desc& desc) {
- return impl_->place(thingify(ls, impl_->morph), desc, impl_->synapses);
+ return impl_->place(ls, desc, impl_->synapses);
}
//
// Stimuli.
//
-lid_range cable_cell::place(const std::string& target, const i_clamp& desc) {
- return impl_->place(target, desc, impl_->stimuli);
-}
-
lid_range cable_cell::place(const locset& ls, const i_clamp& desc) {
- return impl_->place(thingify(ls, impl_->morph), desc, impl_->stimuli);
+ return impl_->place(ls, desc, impl_->stimuli);
}
//
// Gap junctions.
//
-lid_range cable_cell::place(const std::string& target, gap_junction_site) {
- return impl_->place_gj(target);
-}
-
lid_range cable_cell::place(const locset& ls, gap_junction_site) {
- return impl_->place_gj(thingify(ls, impl_->morph));
+ return impl_->place_gj(ls);
}
//
// Spike detectors.
//
-lid_range cable_cell::place(const std::string& target, const threshold_detector& desc) {
- return impl_->place(target, desc.threshold, impl_->spike_detectors);
-}
lid_range cable_cell::place(const locset& ls, const threshold_detector& desc) {
- return impl_->place(thingify(ls, impl_->morph), desc.threshold, impl_->spike_detectors);
+ return impl_->place(ls, desc.threshold, impl_->spike_detectors);
}
//
diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp
index 9a13e1a93ac17a32c69b9387a834653698f9e904..f6843b80b462111474cfbc4cdb6f78934c1112b8 100644
--- a/arbor/cable_cell_param.cpp
+++ b/arbor/cable_cell_param.cpp
@@ -7,7 +7,6 @@
#include <arbor/cable_cell_param.hpp>
#include <arbor/morph/locset.hpp>
-#include "morph/em_morphology.hpp"
#include "util/maputil.hpp"
namespace arb {
@@ -72,8 +71,8 @@ cable_cell_local_parameter_set neuron_parameter_defaults = {
// Discretization policy implementations:
locset cv_policy_max_extent::cv_boundary_points(const cable_cell& cell) const {
- const auto& emorph = *cell.morphology();
- const unsigned nbranch = emorph.num_branches();
+ const unsigned nbranch = cell.morphology().num_branches();
+ const auto& embed = cell.embedding();
if (!nbranch || max_extent_<=0) return ls::nil();
std::vector<mlocation> points;
@@ -87,7 +86,7 @@ locset cv_policy_max_extent::cv_boundary_points(const cable_cell& cell) const {
const double oomax_extent = 1./max_extent_;
while (bidx<nbranch) {
- unsigned ncv = std::ceil(emorph.branch_length(bidx)*oomax_extent);
+ unsigned ncv = std::ceil(embed.branch_length(bidx)*oomax_extent);
double ooncv = 1./ncv;
if (flags_&cv_policy_flag::interior_forks) {
@@ -108,7 +107,7 @@ locset cv_policy_max_extent::cv_boundary_points(const cable_cell& cell) const {
}
locset cv_policy_fixed_per_branch::cv_boundary_points(const cable_cell& cell) const {
- const unsigned nbranch = cell.morphology()->num_branches();
+ const unsigned nbranch = cell.morphology().num_branches();
if (!nbranch) return ls::nil();
std::vector<mlocation> points;
diff --git a/arbor/fvm_compartment.hpp b/arbor/fvm_compartment.hpp
index 4790f595e3f15d07fa192ccd870dda2a8ead0649..9f5713fc9ffb16bb66afaec9a8291c77e8311d7b 100644
--- a/arbor/fvm_compartment.hpp
+++ b/arbor/fvm_compartment.hpp
@@ -7,7 +7,6 @@
#include <arbor/math.hpp>
#include <arbor/util/compat.hpp>
-#include "algorithms.hpp"
#include "util/iterutil.hpp"
#include "util/partition.hpp"
#include "util/rangeutil.hpp"
@@ -88,7 +87,7 @@ public:
template <typename Radii, typename Lengths>
div_compartment_by_ends(size_type n, const Radii& radii, const Lengths& lengths):
oon_(1/value_type(n)),
- length_(algorithms::sum(lengths)),
+ length_(util::sum(lengths)),
ra_(util::front(radii)),
rb_(util::back(radii))
{}
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index d6c8ff4f2129d7fdbac954f74ba832a58a8bfb90..9745d922d2f4f4b5c4c2504571ef639a6a79d156 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -1,7 +1,5 @@
#pragma once
-
-#include <unordered_map>
#include <string>
#include <vector>
@@ -11,6 +9,7 @@
#include <arbor/constants.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/morph/label_dict.hpp>
+#include <arbor/morph/mprovider.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/segment.hpp>
@@ -48,9 +47,6 @@ public:
using value_type = double;
using point_type = point<value_type>;
- using region_map = std::unordered_map<std::string, mcable_list>;
- using locset_map = std::unordered_map<std::string, mlocation_list>;
-
using gap_junction_instance = mlocation;
struct synapse_instance {
@@ -84,6 +80,11 @@ public:
const label_dict& dictionary={},
bool compartments_from_discretization=false);
+ /// Access to morphology and embedding
+ const concrete_embedding& embedding() const;
+ const arb::morphology& morphology() const;
+ const mprovider& provider() const;
+
// the number of branches in the cell
size_type num_branches() const;
@@ -120,48 +121,36 @@ public:
// LEGACY
size_type num_compartments() const;
- //
// Painters and placers.
//
// Used to describe regions and locations where density channels, stimuli,
// synapses, gap juncitons and detectors are located.
- //
// Density channels.
- void paint(const std::string&, mechanism_desc);
void paint(const region&, mechanism_desc);
// Properties.
- void paint(const std::string&, cable_cell_local_parameter_set);
void paint(const region&, cable_cell_local_parameter_set);
// Synapses.
- lid_range place(const std::string&, const mechanism_desc&);
lid_range place(const locset&, const mechanism_desc&);
// Stimuli.
- lid_range place(const std::string&, const i_clamp&);
lid_range place(const locset&, const i_clamp&);
// Gap junctions.
- lid_range place(const std::string&, gap_junction_site);
lid_range place(const locset&, gap_junction_site);
- // spike detectors
- lid_range place(const std::string&, const threshold_detector&);
+ // Spike detectors.
lid_range place(const locset&, const threshold_detector&);
- //
- // access to placed items
- //
+ // Access to placed items.
const std::vector<synapse_instance>& synapses() const;
const std::vector<gap_junction_instance>& gap_junction_sites() const;
const std::vector<detector_instance>& detectors() const;
const std::vector<stimulus_instance>& stimuli() const;
- const em_morphology* morphology() const;
-
// Checks that two cells have the same
// - number and type of segments
// - volume and area properties of each segment
@@ -184,7 +173,6 @@ public:
const cable_cell_parameter_set& global_defaults) const;
private:
-
std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)> impl_;
};
diff --git a/arbor/include/arbor/morph/embed_pwlin.hpp b/arbor/include/arbor/morph/embed_pwlin.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..82d83a570bc6f95d163acd8d838869f1a47a9591
--- /dev/null
+++ b/arbor/include/arbor/morph/embed_pwlin.hpp
@@ -0,0 +1,56 @@
+#pragma once
+
+// Embedding of cell morphology as 1-d tree with piecewise linear radius.
+
+#include <vector>
+
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+
+namespace arb {
+
+struct embed_pwlin_data;
+
+namespace util {
+template <typename X> struct pw_elements;
+}
+
+// Piecewise-constant functions are represented as scalar (double)
+// values defined over contiguous intervals.
+using pw_constant_fn = util::pw_elements<double>;
+
+struct embed_pwlin {
+ explicit embed_pwlin(const arb::morphology& m);
+
+ mlocation sample_location(msize_t sid) const {
+ return sample_locations_.at(sid);
+ }
+
+ // Interpolated radius at location.
+ double radius(mlocation) const;
+
+ // Computed length of mcable.
+ double integrate_length(mcable c) const;
+ double integrate_length(msize_t bid, const pw_constant_fn&) const;
+
+ // Membrane surface area of given mcable.
+ double integrate_area(mcable c) const;
+ double integrate_area(msize_t bid, const pw_constant_fn&) const;
+
+ // Integrated inverse cross-sectional area of given mcable.
+ double integrate_ixa(mcable c) const;
+ double integrate_ixa(msize_t bid, const pw_constant_fn&) const;
+
+ // Length of whole branch.
+ double branch_length(msize_t bid) const {
+ return integrate_length(mcable{bid, 0, 1});
+ }
+
+private:
+ std::vector<mlocation> sample_locations_;
+ std::shared_ptr<embed_pwlin_data> data_;
+};
+
+} // namespace arb
+
+
diff --git a/arbor/include/arbor/morph/error.hpp b/arbor/include/arbor/morph/error.hpp
deleted file mode 100644
index 14c40be496cc6be0b8873072a3e18433876c3a4d..0000000000000000000000000000000000000000
--- a/arbor/include/arbor/morph/error.hpp
+++ /dev/null
@@ -1,15 +0,0 @@
-#pragma once
-
-#include <stdexcept>
-
-#include <arbor/arbexcept.hpp>
-
-namespace arb {
-
-struct morphology_error: public arbor_exception {
- morphology_error(const char* what): arbor_exception(what) {}
- morphology_error(const std::string& what): arbor_exception(what) {}
-};
-
-} // namespace arb
-
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
index a9a64cbb663122577d07fe19ff076eb115efac26..2e39b39f03096c2450f51b24fa364e3a07df5d66 100644
--- a/arbor/include/arbor/morph/locset.hpp
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -13,9 +13,7 @@
namespace arb {
-// Forward declare the backend em_morphology type, required for defining the
-// interface for concretising locsets.
-class em_morphology;
+struct mprovider;
class locset;
@@ -46,6 +44,10 @@ public:
// Construct an explicit location set with a single location.
locset(mlocation other);
+ // Implicitly convert string to named locset expression.
+ locset(std::string label);
+ locset(const char* label);
+
template <typename Impl,
typename X=std::enable_if_t<!std::is_same<std::decay_t<Impl>, locset>::value>>
locset& operator=(Impl&& other) {
@@ -59,7 +61,7 @@ public:
return *this;
}
- friend mlocation_list thingify(const locset& p, const em_morphology& m) {
+ friend mlocation_list thingify(const locset& p, const mprovider& m) {
return p.impl_->thingify(m);
}
@@ -88,7 +90,7 @@ private:
virtual ~interface() {}
virtual std::unique_ptr<interface> clone() = 0;
virtual std::ostream& print(std::ostream&) = 0;
- virtual mlocation_list thingify(const em_morphology&) = 0;
+ virtual mlocation_list thingify(const mprovider&) = 0;
};
std::unique_ptr<interface> impl_;
@@ -102,7 +104,7 @@ private:
return std::unique_ptr<interface>(new wrap<Impl>(wrapped));
}
- virtual mlocation_list thingify(const em_morphology& m) override {
+ virtual mlocation_list thingify(const mprovider& m) override {
return thingify_(wrapped, m);
}
@@ -128,6 +130,9 @@ locset terminal();
// The root node of a morphology.
locset root();
+// Named locset.
+locset named(std::string);
+
// The null (empty) set.
locset nil();
diff --git a/arbor/include/arbor/morph/morphexcept.hpp b/arbor/include/arbor/morph/morphexcept.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c97947daeb420b9cbd01b76221778160cda8e47
--- /dev/null
+++ b/arbor/include/arbor/morph/morphexcept.hpp
@@ -0,0 +1,56 @@
+#pragma once
+
+#include <string>
+
+#include <arbor/arbexcept.hpp>
+#include <arbor/morph/primitives.hpp>
+
+namespace arb {
+
+struct morphology_error: public arbor_exception {
+ morphology_error(const std::string& what): arbor_exception(what) {}
+};
+
+struct invalid_mlocation: morphology_error {
+ invalid_mlocation(mlocation loc);
+ mlocation loc;
+};
+
+struct no_such_branch: morphology_error {
+ no_such_branch(msize_t bid);
+ msize_t bid;
+};
+
+struct invalid_mcable: morphology_error {
+ invalid_mcable(mcable cable);
+ mcable cable;
+};
+
+struct invalid_sample_parent: morphology_error {
+ invalid_sample_parent(msize_t parent, msize_t tree_size);
+ msize_t parent;
+ msize_t tree_size;
+};
+
+struct label_type_mismatch: morphology_error {
+ label_type_mismatch(const std::string& label);
+ std::string label;
+};
+
+struct incomplete_branch: morphology_error {
+ incomplete_branch(msize_t bid);
+ msize_t bid;
+};
+
+struct unbound_name: morphology_error {
+ unbound_name(const std::string& name);
+ std::string name;
+};
+
+struct circular_definition: morphology_error {
+ circular_definition(const std::string& name);
+ std::string name;
+};
+
+} // namespace arb
+
diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index 4e90d30a483cd03ebf53da76450622f02403cb64..bd36d65fb78e7ba3ee718b30b8e1a44b68294840 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -41,6 +41,9 @@ public:
// The child branches of branch b.
const std::vector<msize_t>& branch_children(msize_t b) const;
+ // Branches with no children.
+ const std::vector<msize_t>& terminal_branches() const;
+
// Range of indexes into the sample points in branch b.
mindex_range branch_indexes(msize_t b) const;
@@ -56,4 +59,9 @@ public:
friend std::ostream& operator<<(std::ostream&, const morphology&);
};
+// Morphology utility functions.
+mlocation_list minset(const morphology&, const mlocation_list&);
+
+mlocation canonical(const morphology&, mlocation);
+
} // namespace arb
diff --git a/arbor/include/arbor/morph/mprovider.hpp b/arbor/include/arbor/morph/mprovider.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..727550ccc957cce8700ef872e85bcbff2a8fb5be
--- /dev/null
+++ b/arbor/include/arbor/morph/mprovider.hpp
@@ -0,0 +1,47 @@
+#pragma once
+
+#include <string>
+#include <unordered_map>
+
+#include <arbor/morph/embed_pwlin.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/label_dict.hpp>
+#include <arbor/util/either.hpp>
+
+namespace arb {
+
+using concrete_embedding = embed_pwlin;
+
+struct mprovider {
+ mprovider(arb::morphology m, const label_dict& dict): mprovider(m, &dict) {}
+ explicit mprovider(arb::morphology m): mprovider(m, nullptr) {}
+
+ // Throw exception on missing or recursive definition.
+ const mcable_list& region(const std::string& name) const;
+ const mlocation_list& locset(const std::string& name) const;
+
+ // Read-only access to morphology and constructed embedding.
+ const auto& morphology() const { return morphology_; }
+ const auto& embedding() const { return embedding_; }
+
+private:
+ mprovider(arb::morphology m, const label_dict* ldptr):
+ morphology_(m), embedding_(m), label_dict_ptr(ldptr) { init(); }
+
+ arb::morphology morphology_;
+ concrete_embedding embedding_;
+
+ struct circular_def {};
+
+ // Maps are mutated only during initialization phase of mprovider.
+ mutable std::unordered_map<std::string, util::either<mcable_list, circular_def>> regions_;
+ mutable std::unordered_map<std::string, util::either<mlocation_list, circular_def>> locsets_;
+
+ // Non-null only during initialization phase.
+ mutable const label_dict* label_dict_ptr;
+
+ // Perform greedy initialization of concrete region, locset maps.
+ void init();
+};
+
+} // namespace arb
diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp
index 52f9ea295b4d754682ca18b0d6aa6311e0ee4297..4b9e223b993db2659e94c00d78796c9b5f53f76c 100644
--- a/arbor/include/arbor/morph/primitives.hpp
+++ b/arbor/include/arbor/morph/primitives.hpp
@@ -69,6 +69,11 @@ std::ostream& operator<<(std::ostream& o, const mlocation_list& l);
// and that the locations in the vector are ordered.
bool test_invariants(const mlocation_list&);
+// Multiset operations on location lists.
+mlocation_list sum(const mlocation_list&, const mlocation_list&);
+mlocation_list join(const mlocation_list&, const mlocation_list&);
+mlocation_list intersection(const mlocation_list&, const mlocation_list&);
+
// Describe an unbranched cable in the morphology.
struct mcable {
// The id of the branch on which the cable lies.
diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp
index 1cb884b025bde88fdb3440e4f0bd14a1e9dcb8a9..6f0fd63e2358633f34a464b9eff0af8799234035 100644
--- a/arbor/include/arbor/morph/region.hpp
+++ b/arbor/include/arbor/morph/region.hpp
@@ -13,9 +13,7 @@
namespace arb {
-// Forward declare the backend em_morphology type, required for defining the
-// interface for concretising locsets.
-class em_morphology;
+struct mprovider;
class region {
public:
@@ -54,7 +52,11 @@ public:
return *this;
}
- friend mcable_list thingify(const region& r, const em_morphology& m) {
+ // Implicitly convert string to named region expression.
+ region(std::string label);
+ region(const char* label);
+
+ friend mcable_list thingify(const region& r, const mprovider& m) {
return r.impl_->thingify(m);
}
@@ -83,7 +85,7 @@ private:
virtual ~interface() {}
virtual std::unique_ptr<interface> clone() = 0;
virtual std::ostream& print(std::ostream&) = 0;
- virtual mcable_list thingify(const em_morphology&) = 0;
+ virtual mcable_list thingify(const mprovider&) = 0;
};
std::unique_ptr<interface> impl_;
@@ -97,7 +99,7 @@ private:
return std::unique_ptr<interface>(new wrap<Impl>(wrapped));
}
- virtual mcable_list thingify(const em_morphology& m) override {
+ virtual mcable_list thingify(const mprovider& m) override {
return thingify_(wrapped, m);
}
@@ -128,6 +130,9 @@ region tagged(int id);
// Region with all segments in a cell.
region all();
+// Region associated with a name.
+region named(std::string);
+
} // namespace reg
} // namespace arb
diff --git a/arbor/util/either.hpp b/arbor/include/arbor/util/either.hpp
similarity index 99%
rename from arbor/util/either.hpp
rename to arbor/include/arbor/util/either.hpp
index faea135a63a4e6c04e49071872f6d5e41c00df3c..08cc66122f8e20c7f604b28fec20572cf5990bdd 100644
--- a/arbor/util/either.hpp
+++ b/arbor/include/arbor/util/either.hpp
@@ -13,8 +13,6 @@
#include <arbor/util/uninitialized.hpp>
-#include "util/meta.hpp"
-
namespace arb {
namespace util {
@@ -97,7 +95,7 @@ namespace detail {
};
} // namespace detail
-constexpr std::size_t variant_npos = static_cast<std::size_t>(-1); // emulating C++17 variant type
+constexpr std::size_t variant_npos = static_cast<std::size_t>(-1); // emulating C++17 variant type
template <typename A, typename B>
class either: public detail::either_data<A, B> {
diff --git a/arbor/mechcat.cpp b/arbor/mechcat.cpp
index 1c2a5d0bc213671a6ee656d648d48533f4ec656f..8d2a9c19ab2cd5a2220d02ea289b55f3ac9d4296 100644
--- a/arbor/mechcat.cpp
+++ b/arbor/mechcat.cpp
@@ -6,8 +6,8 @@
#include <arbor/arbexcept.hpp>
#include <arbor/mechcat.hpp>
+#include <arbor/util/either.hpp>
-#include "util/either.hpp"
#include "util/maputil.hpp"
/* Notes on implementation:
diff --git a/arbor/morph/em_morphology.cpp b/arbor/morph/em_morphology.cpp
deleted file mode 100644
index 8e49c3a9d42401d95038d5fbbeea28edeb23237a..0000000000000000000000000000000000000000
--- a/arbor/morph/em_morphology.cpp
+++ /dev/null
@@ -1,194 +0,0 @@
-#include <mutex>
-#include <stack>
-#include <vector>
-
-#include <arbor/morph/error.hpp>
-#include <arbor/morph/morphology.hpp>
-#include <arbor/morph/primitives.hpp>
-
-#include "morph/em_morphology.hpp"
-#include "util/rangeutil.hpp"
-#include "util/span.hpp"
-#include "util/strprintf.hpp"
-
-namespace arb {
-
-em_morphology::em_morphology(const morphology& m):
- morph_(m)
-{
- using util::count_along;
- using util::make_span;
-
- auto& parents = morph_.sample_parents();
- auto& samples = morph_.samples();
-
- const auto ns = morph_.num_samples();
- const auto nb = morph_.num_branches();
-
- if (!ns) return;
-
- // Cache distance of each sample from the root.
- dist2root_.resize(ns);
- dist2root_[0] = 0.;
- for (auto i: make_span(1, ns)) {
- const auto p = parents[i];
- dist2root_[i] = dist2root_[p] + distance(samples[p], samples[i]);
- }
- // Cache the legth of each branch.
- branch_lengths_.reserve(nb);
- for (auto i: make_span(nb)) {
- auto idx = util::make_range(morph_.branch_indexes(i));
- branch_lengths_.push_back(dist2root_[idx.back()]- dist2root_[idx.front()]);
- }
- if (morph_.spherical_root()) {
- branch_lengths_[0] = samples[0].loc.radius*2;
- }
-
- // Cache the sample locations.
- // Iterate backwards over branches distal to root, so that the parent branch at
- // fork points will label its distal sample.
- sample_locs_.resize(ns);
- for (int b=nb-1; b>=0; --b) {
- auto idx = util::make_range(morph_.branch_indexes(b));
- double len = branch_lengths_[b];
- // Handle 0 length branch.
- len = len==0.? 1.: len;
- double start = dist2root_[idx.front()];
- for (auto i: idx) {
- sample_locs_[i] = {msize_t(b), (dist2root_[i]-start)/len};
- }
- // Ensure that all non-spherical branches have their last sample
- if (idx.size()>1u) {
- sample_locs_[idx.back()] = mlocation{msize_t(b), 1};
- }
- }
- sample_locs_[0] = mlocation{0, 0.};
-
- // Cache the location of terminal and fork points.
- auto& props = morph_.sample_props();
- for (auto i: count_along(props)) {
- auto p = props[i];
- if (is_terminal(p)) {
- terminals_.push_back(sample2loc(i));
- }
- if (is_fork(p)) {
- forks_.push_back(sample2loc(i));
- }
- }
-}
-
-em_morphology::em_morphology():
- em_morphology(morphology())
-{}
-
-const morphology& em_morphology::morph() const {
- return morph_;
-}
-
-mlocation em_morphology::root() const {
- return {0,0};
-}
-
-mlocation em_morphology::sample2loc(msize_t sid) const {
- if (sid>=morph_.num_samples()) {
- throw morphology_error(util::pprintf("Sample {} does not exist in morpology", sid));
- }
- return sample_locs_[sid];
-}
-
-mlocation_list em_morphology::terminals() const {
- return terminals_;
-}
-
-mlocation_list em_morphology::cover(mlocation loc, bool include_loc) const {
- mlocation_list L{};
- if (include_loc) L.push_back(loc);
-
- // If the location is not at the end of a branch, it is its own cover.
- if (loc.pos>0. && loc.pos<1.) return L;
-
- // First location is {0,0} on a spherical root: nothing more to do.
- if (loc==mlocation{0,0} && morph_.spherical_root()) {
- return L;
- }
-
- if (loc.pos==1) {
- // The location is at the end of a branch:
- // add the location at the start of each child branch.
- for (auto b: morph_.branch_children(loc.branch)) {
- L.push_back({b, 0});
- }
- }
- else if (loc.pos==0) {
- // The location is at the start of a branch:
- // add the location at the end of the parent branch, and locations
- // at the start of each child of the parent branch.
- auto p = morph_.branch_parent(loc.branch);
- if (p!=mnpos) L.push_back({p, 1});
- for (auto b: morph_.branch_children(p)) {
- if (b!=loc.branch) L.push_back({b, 0});
- }
- }
-
- util::sort(L);
-
- return L;
-}
-
-void em_morphology::assert_valid_location(mlocation loc) const {
- if (!test_invariants(loc)) {
- throw morphology_error(util::pprintf("Invalid location {}", loc));
- }
- if (loc.branch>=morph_.num_branches()) {
- throw morphology_error(util::pprintf("Location {} does not exist in morpology", loc));
- }
-}
-
-mlocation em_morphology::canonicalize(mlocation loc) const {
- assert_valid_location(loc);
-
- // Test if location is at the start of a branch.
- if (loc.pos==0.) {
- auto p = morph_.branch_parent(loc.branch);
- return p==mnpos? root(): mlocation{p, 1};
- }
- return loc;
-}
-
-mlocation_list em_morphology::minset(const mlocation_list& in) const {
- mlocation_list L;
-
- std::stack<msize_t> stack;
-
- // All root branches must be searched.
- for (auto c: morph_.branch_children(mnpos)) {
- stack.push(c);
- }
-
- // Depth-first traversal of the branch tree.
- while (!stack.empty()) {
- auto branch = stack.top();
- stack.pop();
-
- // Search for a location on branch.
- auto it = std::lower_bound(in.begin(), in.end(), mlocation{branch, 0});
-
- // If found, insert to the minset and skip the rest of this sub-tree.
- if (it!=in.end() && it->branch==branch) {
- L.push_back(*it);
- continue;
- }
-
- // No location on this branch, so continue searching in this sub-tree.
- for (auto c: morph_.branch_children(branch)) {
- stack.push(c);
- }
- }
-
- util::sort(L);
-
- return L;
-}
-
-} // namespace arb
-
diff --git a/arbor/morph/em_morphology.hpp b/arbor/morph/em_morphology.hpp
deleted file mode 100644
index 1feb88cda374c17985a269856e2857bce91591d4..0000000000000000000000000000000000000000
--- a/arbor/morph/em_morphology.hpp
+++ /dev/null
@@ -1,57 +0,0 @@
-#pragma once
-
-#include <vector>
-
-#include <arbor/morph/morphology.hpp>
-#include <arbor/morph/primitives.hpp>
-
-namespace arb {
-
-class em_morphology {
- morphology morph_;
-
- mlocation_list terminals_;
- mlocation_list forks_;
- mlocation_list sample_locs_;
-
- // distance from sample to root
- std::vector<double> dist2root_;
- std::vector<double> branch_lengths_;
-
-public:
- em_morphology();
- em_morphology(const morphology& m);
-
- const morphology& morph() const;
-
- // Convenience methods for morphology access
- // that are forwarded directly to the morphology object:
-
- auto empty() const { return morph_.empty(); }
- auto spherical_root() const { return morph_.spherical_root(); }
- auto num_branches() const { return morph_.num_branches(); }
- auto num_samples() const { return morph_.num_samples(); }
- auto branch_parent(msize_t b) const { return morph_.branch_parent(b); }
- auto branch_children(msize_t b) const { return morph_.branch_children(b); }
-
- // Access to computed and cached data:
-
- mlocation_list terminals() const;
- mlocation root() const;
-
- mlocation sample2loc(msize_t sid) const;
-
- void assert_valid_location(mlocation) const;
- mlocation canonicalize(mlocation) const;
-
- // Find all locations on the morphology that share the same canonoical
- // representation of loc.
- // If include_loc is false, the input location is excluded from the result.
- mlocation_list cover(mlocation, bool include_loc=true) const;
-
- mlocation_list minset(const mlocation_list&) const;
-
- double branch_length(msize_t bid) const { return branch_lengths_.at(bid); }
-};
-
-} // namespace arb
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e04436b4598224ab3a975ac60f5c527e40fd11b2
--- /dev/null
+++ b/arbor/morph/embed_pwlin.cpp
@@ -0,0 +1,204 @@
+#include <cstddef>
+#include <utility>
+#include <vector>
+
+#include <arbor/morph/embed_pwlin.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+
+#include "util/piecewise.hpp"
+#include "util/range.hpp"
+#include "util/rangeutil.hpp"
+#include "util/ratelem.hpp"
+#include "util/span.hpp"
+
+namespace arb {
+
+using util::rat_element;
+
+template <unsigned p, unsigned q>
+using pw_ratpoly = util::pw_elements<rat_element<p, q>>;
+
+template <unsigned p, unsigned q>
+using branch_pw_ratpoly = std::vector<pw_ratpoly<p, q>>;
+
+template <unsigned p, unsigned q>
+double interpolate(const branch_pw_ratpoly<p, q>& f, unsigned bid, double pos) {
+ const auto& pw = f.at(bid);
+ unsigned index = pw.index_of(pos);
+
+ const auto& element = pw.element(index);
+ std::pair<double, double> bounds = pw.interval(index);
+
+ if (bounds.first==bounds.second) return element[0];
+ else {
+ double x = (pos-bounds.first)/(bounds.second-bounds.first);
+ return element(x);
+ }
+}
+
+// Length, area, and ixa are polynomial or rational polynomial functions of branch position,
+// continuos and monotonically increasing with respect to distance from root.
+//
+// Integration wrt a piecewise constant function is performed by taking the difference between
+// interpolated values at the end points of each constant interval.
+
+template <unsigned p, unsigned q>
+double integrate(const branch_pw_ratpoly<p, q>& f, unsigned bid, const pw_constant_fn& g) {
+ double accum = 0;
+ for (msize_t i = 0; i<g.size(); ++i) {
+ std::pair<double, double> interval = g.interval(i);
+ accum += g.element(i)*(interpolate(f, bid, interval.second)-interpolate(f, bid, interval.first));
+ }
+ return accum;
+}
+
+struct embed_pwlin_data {
+ branch_pw_ratpoly<1, 0> length;
+ branch_pw_ratpoly<1, 0> radius;
+ branch_pw_ratpoly<2, 0> area;
+ branch_pw_ratpoly<1, 1> ixa;
+
+ explicit embed_pwlin_data(msize_t n_branch):
+ length(n_branch),
+ radius(n_branch),
+ area(n_branch),
+ ixa(n_branch)
+ {}
+};
+
+double embed_pwlin::radius(mlocation loc) const {
+ return interpolate(data_->radius, loc.branch, loc.pos);
+}
+
+double embed_pwlin::integrate_length(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(data_->length, bid, g);
+}
+
+double embed_pwlin::integrate_area(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(data_->area, bid, g);
+}
+
+double embed_pwlin::integrate_ixa(msize_t bid, const pw_constant_fn& g) const {
+ return integrate(data_->ixa, bid, g);
+}
+
+// Cable versions of integration methods:
+
+double embed_pwlin::integrate_length(mcable c) const {
+ return integrate_length(c.branch, pw_constant_fn{{c.prox_pos, c.dist_pos}, {1.}});
+}
+
+double embed_pwlin::integrate_area(mcable c) const {
+ return integrate_area(c.branch, pw_constant_fn{{c.prox_pos, c.dist_pos}, {1.}});
+}
+
+double embed_pwlin::integrate_ixa(mcable c) const {
+ return integrate_ixa(c.branch, pw_constant_fn{{c.prox_pos, c.dist_pos}, {1.}});
+}
+
+// Initialization, creation of geometric data.
+
+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);
+
+ if (!n_branch) return;
+
+ const auto& samples = m.samples();
+ sample_locations_.resize(m.num_samples());
+
+ for (msize_t bid = 0; bid<n_branch; ++bid) {
+ unsigned parent = m.branch_parent(bid);
+ auto sample_indices = util::make_range(m.branch_indexes(bid));
+
+ if (bid==0 && m.spherical_root()) {
+ arb_assert(sample_indices.size()==1);
+
+ // Treat spherical root as area-equivalent cylinder.
+ double r = samples[0].loc.radius;
+
+ data_->length[bid].push_back(0., 1., rat_element<1, 0>(0, r*2));
+ data_->radius[bid].push_back(0., 1., rat_element<1, 0>(r, r));
+
+ double cyl_area = 4*pi*r*r;
+ data_->area[bid].push_back(0., 1., rat_element<2, 0>(0., cyl_area*0.5, cyl_area));
+
+ double cyl_ixa = 2.0/(pi*r);
+ data_->ixa[bid].push_back(0., 1., rat_element<1, 1>(0., cyl_ixa*0.5, cyl_ixa));
+
+ sample_locations_[0] = mlocation{0, 0.5};
+ }
+ else {
+ arb_assert(sample_indices.size()>1);
+
+ std::vector<double> sample_distance;
+ sample_distance.reserve(samples.size());
+ sample_distance.push_back(0.);
+
+ for (auto i: util::count_along(sample_indices)) {
+ if (!i) continue;
+
+ double d = distance(samples[sample_indices[i-1]], samples[sample_indices[i]]);
+ sample_distance.push_back(sample_distance.back()+d);
+ }
+
+ double branch_length = sample_distance.back();
+ double length_scale = branch_length>0? 1./branch_length: 0;
+
+ for (auto i: util::count_along(sample_indices)) {
+ sample_locations_[sample_indices[i]] = mlocation{bid, length_scale*sample_distance[i]};
+ }
+ sample_locations_[sample_indices.back()].pos = 1.; // Circumvent any rounding infelicities.
+
+ 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));
+
+ double area_0 = parent=mnpos? 0: data_->area[parent].back().second[1];
+ double ixa_0 = parent=mnpos? 0: data_->ixa[parent].back().second[1];
+
+ if (length_scale==0) {
+ // Zero-length branch? Weird, but make best show of it.
+ double r = samples[sample_indices[0]].loc.radius;
+ data_->radius[bid].push_back(0., 1., rat_element<1, 0>(r, r));
+ data_->area[bid].push_back(0., 1., rat_element<2, 0>(area_0, area_0, area_0));
+ data_->ixa[bid].push_back(0., 1., rat_element<1, 1>(ixa_0, ixa_0, ixa_0));
+ }
+ else {
+ for (auto i: util::count_along(sample_indices)) {
+ if (!i) continue;
+
+ double p0 = i>1? sample_locations_[sample_indices[i-1]].pos: 0;
+ double p1 = sample_locations_[sample_indices[i]].pos;
+ if (p0==p1) continue;
+
+ double r0 = samples[sample_indices[i-1]].loc.radius;
+ double r1 = samples[sample_indices[i]].loc.radius;
+ data_->radius[bid].push_back(p0, p1, rat_element<1, 0>(r0, r1));
+
+ double dx = (p1-p0)*branch_length;
+ double dr = r1-r0;
+ double c = pi*std::sqrt(dr*dr+dx*dx);
+ double area_half = area_0 + (0.75*r0+0.25*r1)*c;
+ double area_1 = area_0 + (r0+r1)*c;
+ data_->area[bid].push_back(p0, p1, rat_element<2, 0>(area_0, area_half, area_1));
+ area_0 = area_1;
+
+ double ixa_half = ixa_0 + dx/(pi*r0*(r0+r1));
+ double ixa_1 = ixa_0 + dx/(pi*r0*r1);
+ data_->ixa[bid].push_back(p0, p1, rat_element<1, 1>(ixa_0, ixa_half, ixa_1));
+ ixa_0 = ixa_1;
+ }
+ }
+
+ arb_assert((data_->radius[bid].size()>0));
+ arb_assert((data_->radius[bid].bounds()==std::pair<double, double>(0., 1.)));
+ arb_assert((data_->area[bid].bounds()==std::pair<double, double>(0., 1.)));
+ arb_assert((data_->ixa[bid].bounds()==std::pair<double, double>(0., 1.)));
+ arb_assert(sample_locations_.size()==m.samples().size());
+ }
+ }
+};
+
+} // namespace arb
diff --git a/arbor/morph/label_dict.cpp b/arbor/morph/label_dict.cpp
index e72c40067ebf50e5618fe65d25ca03e1b7f61a03..9b00acbc65ecf27be82984551f8de1f5c43e55b2 100644
--- a/arbor/morph/label_dict.cpp
+++ b/arbor/morph/label_dict.cpp
@@ -1,18 +1,10 @@
-#include <algorithm>
-#include <set>
-#include <fstream>
#include <unordered_map>
-#include <arbor/morph/error.hpp>
-#include <arbor/morph/locset.hpp>
+#include <arbor/morph/morphexcept.hpp>
#include <arbor/morph/label_dict.hpp>
+#include <arbor/morph/locset.hpp>
#include <arbor/morph/region.hpp>
-#include "util/filter.hpp"
-#include "util/span.hpp"
-#include "util/strprintf.hpp"
-#include "io/sepval.hpp"
-
namespace arb {
size_t label_dict::size() const {
@@ -21,16 +13,14 @@ size_t label_dict::size() const {
void label_dict::set(const std::string& name, arb::locset ls) {
if (regions_.count(name)) {
- throw morphology_error(util::pprintf(
- "Attempt to add a locset \"{}\" to a label dictionary that already contains a region with the same name.", name));
+ throw label_type_mismatch(name);
}
locsets_[name] = std::move(ls);
}
void label_dict::set(const std::string& name, arb::region reg) {
if (locsets_.count(name)) {
- throw morphology_error(util::pprintf(
- "Attempt to add a region \"{}\" to a label dictionary that already contains a locset with the same name.", name));
+ throw label_type_mismatch(name);
}
regions_[name] = std::move(reg);
}
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index c21bb3d69ded5d2b146f525d1217ebcb99c3a50f..dd81e56e4ef6236a69a23fa9dd0cfae1d68c51c0 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -2,109 +2,34 @@
#include <iostream>
#include <numeric>
-#include <arbor/morph/error.hpp>
#include <arbor/morph/locset.hpp>
-#include <arbor/morph/region.hpp>
+#include <arbor/morph/morphexcept.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/mprovider.hpp>
+#include <arbor/morph/primitives.hpp>
#include "util/rangeutil.hpp"
#include "util/strprintf.hpp"
-#include "morph/em_morphology.hpp"
-
namespace arb {
namespace ls {
-//
-// Functions for taking the sum, union and intersection of location_lists (multisets).
-//
-
-using it_t = mlocation_list::iterator;
-using const_it_t = mlocation_list::const_iterator;
-
-// Advance an iterator to the first value that is not equal to its current
-// value, or end, whichever comes first.
-template <typename T>
-T next_unique(T& it, T end) {
- const auto& x = *it;
- ++it;
- while (it!=end && *it==x) ++it;
- return it;
-};
-
-// Return the number of times that the value at it is repeated. Advances the
-// iterator to the first value not equal to its current value, or end,
-// whichever comse first.
-template <typename T>
-int multiplicity(T& it, T end) {
- const auto b = it;
- return std::distance(b, next_unique(it, end));
-};
-
-mlocation_list sum(const mlocation_list& lhs, const mlocation_list& rhs) {
- mlocation_list v;
- v.resize(lhs.size() + rhs.size());
- std::merge(lhs.begin(), lhs.end(), rhs.begin(), rhs.end(), v.begin());
- return v;
-}
-
-mlocation_list join(const mlocation_list& lhs, const mlocation_list& rhs) {
- mlocation_list L;
- L.reserve(lhs.size()+rhs.size());
-
- auto l = lhs.begin();
- auto lend = lhs.end();
- auto r = rhs.begin();
- auto rend = rhs.end();
-
- auto at_end = [&]() { return l==lend || r==rend; };
- while (!at_end()) {
- auto x = (*l<*r) ? *l: *r;
- auto count = (*l<*r)? multiplicity(l, lend):
- (*r<*l)? multiplicity(r, rend):
- std::max(multiplicity(l, lend), multiplicity(r, rend));
- L.insert(L.end(), count, x);
+// Throw on invalid mlocation.
+void assert_valid(mlocation x) {
+ if (!test_invariants(x)) {
+ throw invalid_mlocation(x);
}
- L.insert(L.end(), l, lend);
- L.insert(L.end(), r, rend);
-
- return L;
-}
-
-mlocation_list intersection(const mlocation_list& lhs, const mlocation_list& rhs) {
- mlocation_list L;
- L.reserve(lhs.size()+rhs.size());
-
- auto l = lhs.begin();
- auto lend = lhs.end();
- auto r = rhs.begin();
- auto rend = rhs.end();
-
- auto at_end = [&]() { return l==lend || r==rend; };
- while (!at_end()) {
- if (*l==*r) {
- auto x = *l;
- auto count = std::min(multiplicity(l, lend), multiplicity(r, rend));
- L.insert(L.end(), count, x);
- }
- else if (*l<*r) {
- next_unique(l, lend);
- }
- else {
- next_unique(r, rend);
- }
- }
-
- return L;
}
-// Null set
+// Empty locset.
+
struct nil_ {};
locset nil() {
return locset{nil_{}};
}
-mlocation_list thingify_(const nil_& x, const em_morphology& m) {
+mlocation_list thingify_(const nil_& x, const mprovider&) {
return {};
}
@@ -112,20 +37,22 @@ std::ostream& operator<<(std::ostream& o, const nil_& x) {
return o << "nil";
}
-// An explicit location
+// An explicit location.
+
struct location_ {
mlocation loc;
};
locset location(mlocation loc) {
- if (!test_invariants(loc)) {
- throw morphology_error(util::pprintf("invalid location {}", loc));
- }
+ assert_valid(loc);
return locset{location_{loc}};
}
-mlocation_list thingify_(const location_& x, const em_morphology& m) {
- m.assert_valid_location(x.loc);
+mlocation_list thingify_(const location_& x, const mprovider& p) {
+ assert_valid(x.loc);
+ if (x.loc.branch>=p.morphology().num_branches()) {
+ throw no_such_branch(x.loc.branch);
+ }
return {x.loc};
}
@@ -134,7 +61,8 @@ std::ostream& operator<<(std::ostream& o, const location_& x) {
}
-// Location corresponding to a sample id
+// Location corresponding to a sample id.
+
struct sample_ {
msize_t index;
};
@@ -143,83 +71,111 @@ locset sample(msize_t index) {
return locset{sample_{index}};
}
-mlocation_list thingify_(const sample_& x, const em_morphology& m) {
- return {m.sample2loc(x.index)};
+mlocation_list thingify_(const sample_& x, const mprovider& p) {
+ return {canonical(p.morphology(), p.embedding().sample_location(x.index))};
}
std::ostream& operator<<(std::ostream& o, const sample_& x) {
return o << "(sample " << x.index << ")";
}
-// set of terminal nodes on a morphology
+// Set of terminal points (most distal points).
+
struct terminal_ {};
locset terminal() {
return locset{terminal_{}};
}
-mlocation_list thingify_(const terminal_&, const em_morphology& m) {
- return m.terminals();
+mlocation_list thingify_(const terminal_&, const mprovider& p) {
+ mlocation_list locs;
+ util::assign(locs, util::transform_view(p.morphology().terminal_branches(),
+ [](msize_t bid) { return mlocation{bid, 1.}; }));
+
+ return locs;
}
std::ostream& operator<<(std::ostream& o, const terminal_& x) {
return o << "terminal";
}
-// the root node of a morphology
+// Root location (most proximal point).
+
struct root_ {};
locset root() {
return locset{root_{}};
}
-mlocation_list thingify_(const root_&, const em_morphology& m) {
- return {m.root()};
+mlocation_list thingify_(const root_&, const mprovider& p) {
+ return {mlocation{0, 0.}};
}
std::ostream& operator<<(std::ostream& o, const root_& x) {
return o << "root";
}
-// intersection of two point sets
+// Named locset.
+
+struct named_ {
+ std::string name;
+};
+
+locset named(std::string name) {
+ return locset(named_{std::move(name)});
+}
+
+mlocation_list thingify_(const named_& n, const mprovider& p) {
+ return p.locset(n.name);
+}
+
+std::ostream& operator<<(std::ostream& o, const named_& x) {
+ return o << "(named \"" << x.name << "\")";
+}
+
+
+// Intersection of two point sets.
+
struct land {
locset lhs;
locset rhs;
land(locset lhs, locset rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {}
};
-mlocation_list thingify_(const land& P, const em_morphology& m) {
- return intersection(thingify(P.lhs, m), thingify(P.rhs, m));
+mlocation_list thingify_(const land& P, const mprovider& p) {
+ return intersection(thingify(P.lhs, p), thingify(P.rhs, p));
}
std::ostream& operator<<(std::ostream& o, const land& x) {
return o << "(intersect " << x.lhs << " " << x.rhs << ")";
}
-// union of two point sets
+// Union of two point sets.
+
struct lor {
locset lhs;
locset rhs;
lor(locset lhs, locset rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {}
};
-mlocation_list thingify_(const lor& P, const em_morphology& m) {
- return join(thingify(P.lhs, m), thingify(P.rhs, m));
+mlocation_list thingify_(const lor& P, const mprovider& p) {
+ return join(thingify(P.lhs, p), thingify(P.rhs, p));
}
std::ostream& operator<<(std::ostream& o, const lor& x) {
return o << "(join " << x.lhs << " " << x.rhs << ")";
}
-// sum of two point sets
+// Sum of two point sets.
+
struct lsum {
locset lhs;
locset rhs;
lsum(locset lhs, locset rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {}
};
-mlocation_list thingify_(const lsum& P, const em_morphology& m) {
- return sum(thingify(P.lhs, m), thingify(P.rhs, m));
+mlocation_list thingify_(const lsum& P, const mprovider& p) {
+ return sum(thingify(P.lhs, p), thingify(P.rhs, p));
}
std::ostream& operator<<(std::ostream& o, const lsum& x) {
@@ -248,8 +204,17 @@ locset::locset() {
*this = ls::nil();
}
-locset::locset(mlocation other) {
- *this = ls::location(other);
+locset::locset(mlocation loc) {
+ *this = ls::location(loc);
+}
+
+locset::locset(std::string name) {
+ *this = ls::named(std::move(name));
+}
+
+locset::locset(const char* name) {
+ *this = ls::named(name);
}
+
} // namespace arb
diff --git a/arbor/morph/morphexcept.cpp b/arbor/morph/morphexcept.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..21463f20c185bacee26c211109f16397946dc00e
--- /dev/null
+++ b/arbor/morph/morphexcept.cpp
@@ -0,0 +1,55 @@
+#include <string>
+#include <sstream>
+
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/morphexcept.hpp>
+
+#include "util/strprintf.hpp"
+
+namespace arb {
+
+using arb::util::pprintf;
+
+invalid_mlocation::invalid_mlocation(mlocation loc):
+ morphology_error(pprintf("invalid mlocation {}", loc)),
+ loc(loc)
+{}
+
+no_such_branch::no_such_branch(msize_t bid):
+ morphology_error(pprintf("no such branch id {}", bid)),
+ bid(bid)
+{}
+
+invalid_mcable::invalid_mcable(mcable cable):
+ morphology_error(pprintf("invalid mcable {}", cable)),
+ cable(cable)
+{}
+
+invalid_sample_parent::invalid_sample_parent(msize_t parent, msize_t tree_size):
+ morphology_error(pprintf("invalid sample parent {} for a sample tree of size {}", parent, tree_size)),
+ parent(parent),
+ tree_size(tree_size)
+{}
+
+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)
+{}
+
+incomplete_branch::incomplete_branch(msize_t bid):
+ morphology_error(pprintf("insufficent samples to define branch id {}", bid)),
+ bid(bid)
+{}
+
+unbound_name::unbound_name(const std::string& name):
+ morphology_error(pprintf("no definition for '{}'", name)),
+ name(name)
+{}
+
+circular_definition::circular_definition(const std::string& name):
+ morphology_error(pprintf("definition of '{}' requires a definition for '{}'", name, name)),
+ name(name)
+{}
+
+} // namespace arb
+
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
index 34253041750bba913d66d0acf25bd4adc4401b94..3b59fef868ee5c0999ba03186427ab52dda21d59 100644
--- a/arbor/morph/morphology.cpp
+++ b/arbor/morph/morphology.cpp
@@ -1,42 +1,31 @@
-#include <cmath>
#include <iostream>
-#include <unordered_map>
+#include <stack>
#include <utility>
-#include <arbor/math.hpp>
-#include <arbor/morph/error.hpp>
+#include <arbor/morph/morphexcept.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/sample_tree.hpp>
#include <arbor/morph/primitives.hpp>
-#include "algorithms.hpp"
-#include "io/sepval.hpp"
#include "morph/mbranch.hpp"
+#include "util/rangeutil.hpp"
#include "util/span.hpp"
-#include "util/strprintf.hpp"
-namespace arb {
+using arb::util::make_span;
+namespace arb {
namespace impl{
std::vector<mbranch> branches_from_parent_index(const std::vector<msize_t>& parents,
const std::vector<point_prop>& props,
bool spherical_root)
{
- using util::make_span;
-
- const char* errstr_single_sample_root =
- "A morphology with only one sample must have a spherical root";
- const char* errstr_incomplete_cable =
- "A branch must contain at least two samples";
-
- if (parents.empty()) return {};
-
auto nsamp = parents.size();
+ if (!nsamp) return {};
// Enforce that a morphology with one sample has a spherical root.
if (!spherical_root && nsamp==1u) {
- throw morphology_error(errstr_single_sample_root);
+ throw incomplete_branch(0);
}
std::vector<int> bids(nsamp);
@@ -70,7 +59,7 @@ std::vector<mbranch> branches_from_parent_index(const std::vector<msize_t>& pare
if (spherical_root) {
for (auto i: make_span(1, nbranches)) { // skip the root.
if (branches[i].size()<2u) {
- throw morphology_error(errstr_incomplete_cable);
+ throw incomplete_branch(i);
}
}
}
@@ -109,6 +98,7 @@ struct morphology_impl {
std::vector<impl::mbranch> branches_;
std::vector<msize_t> branch_parents_;
std::vector<msize_t> root_children_;
+ std::vector<msize_t> terminal_branches_;
std::vector<std::vector<msize_t>> branch_children_;
morphology_impl(sample_tree m, bool use_spherical_root);
@@ -134,11 +124,7 @@ morphology_impl::morphology_impl(sample_tree m):
}
void morphology_impl::init() {
- using util::make_span;
- using util::count_along;
-
auto nsamp = samples_.size();
-
if (!nsamp) return;
// Generate branches.
@@ -158,6 +144,15 @@ void morphology_impl::init() {
root_children_.push_back(i);
}
}
+
+ // Collect terminal branches.
+ terminal_branches_.reserve(nbranch);
+ for (auto i: make_span(nbranch)) {
+ if (branch_children_[i].empty()) {
+ terminal_branches_.push_back(i);
+ }
+ }
+ terminal_branches_.shrink_to_fit();
}
std::ostream& operator<<(std::ostream& o, const morphology_impl& m) {
@@ -205,6 +200,10 @@ const std::vector<msize_t>& morphology::branch_children(msize_t b) const {
return b==mnpos? impl_->root_children_: impl_->branch_children_[b];
}
+const std::vector<msize_t>& morphology::terminal_branches() const {
+ return impl_->terminal_branches_;
+}
+
// Whether the root of the morphology is spherical.
bool morphology::spherical_root() const {
return impl_->spherical_root_;
@@ -236,5 +235,49 @@ std::ostream& operator<<(std::ostream& o, const morphology& m) {
return o << *m.impl_;
}
+// Utilities.
+
+mlocation_list minset(const morphology& m, const mlocation_list& in) {
+ mlocation_list L;
+
+ std::stack<msize_t> stack;
+
+ // All root branches must be searched.
+ for (auto c: m.branch_children(mnpos)) {
+ stack.push(c);
+ }
+
+ // Depth-first traversal of the branch tree.
+ while (!stack.empty()) {
+ auto branch = stack.top();
+ stack.pop();
+
+ // Search for a location on branch.
+ auto it = std::lower_bound(in.begin(), in.end(), mlocation{branch, 0});
+
+ // If found, insert to the minset and skip the rest of this sub-tree.
+ if (it!=in.end() && it->branch==branch) {
+ L.push_back(*it);
+ continue;
+ }
+
+ // No location on this branch, so continue searching in this sub-tree.
+ for (auto c: m.branch_children(branch)) {
+ stack.push(c);
+ }
+ }
+
+ util::sort(L);
+ return L;
+}
+
+mlocation canonical(const morphology& m, mlocation loc) {
+ if (loc.pos==0) {
+ msize_t parent = m.branch_parent(loc.branch);
+ return parent==mnpos? mlocation{0, 0.}: mlocation{parent, 1.};
+ }
+ return loc;
+}
+
} // namespace arb
diff --git a/arbor/morph/mprovider.cpp b/arbor/morph/mprovider.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8b8cb216c91e519e9007b05d62728412581c71c9
--- /dev/null
+++ b/arbor/morph/mprovider.cpp
@@ -0,0 +1,75 @@
+#include <string>
+#include <utility>
+
+#include <arbor/morph/label_dict.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/morphexcept.hpp>
+#include <arbor/morph/mprovider.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/region.hpp>
+
+namespace arb {
+
+void mprovider::init() {
+ // Evaluate each named region or locset in provided dictionary
+ // to populate concrete regions_, locsets_ maps.
+
+ if (!label_dict_ptr) return;
+
+ for (const auto& pair: label_dict_ptr->regions()) {
+ (void)region(pair.first);
+ }
+
+ for (const auto& pair: label_dict_ptr->locsets()) {
+ (void)locset(pair.first);
+ }
+
+ label_dict_ptr = nullptr;
+}
+
+// Evaluation of a named region or locset requires the recursive evaluation of
+// any component regions or locsets in its definition.
+//
+// During the initialization phase, 'named' expressions will be looked up in the
+// provided label_dict, and the maps updated accordingly. Post-initialization,
+// label_dict_ptr will be null, and concrete regions/locsets will only be retrieved
+// from the maps established during initialization.
+
+template <typename RegOrLocMap, typename LabelDictMap, typename Err>
+static const auto& try_lookup(const mprovider& provider, const std::string& name, RegOrLocMap& map, const LabelDictMap* dict_ptr, Err errval) {
+ auto it = map.find(name);
+ if (it==map.end()) {
+ if (dict_ptr) {
+ map.emplace(name, errval);
+
+ auto it = dict_ptr->find(name);
+ if (it==dict_ptr->end()) {
+ throw unbound_name(name);
+ }
+
+ return (map[name] = thingify(it->second, provider)).first();
+ }
+ else {
+ throw unbound_name(name);
+ }
+ }
+ else if (!it->second) {
+ throw circular_definition(name);
+ }
+ else {
+ return it->second.first();
+ }
+}
+
+const mcable_list& mprovider::region(const std::string& name) const {
+ const auto* regions_ptr = label_dict_ptr? &(label_dict_ptr->regions()): nullptr;
+ return try_lookup(*this, name, regions_, regions_ptr, circular_def{});
+}
+
+const mlocation_list& mprovider::locset(const std::string& name) const {
+ const auto* locsets_ptr = label_dict_ptr? &(label_dict_ptr->locsets()): nullptr;
+ return try_lookup(*this, name, locsets_, locsets_ptr, circular_def{});
+}
+
+
+} // namespace arb
diff --git a/arbor/morph/primitives.cpp b/arbor/morph/primitives.cpp
index 486959bb60cd2cef279dd54ec5d6806c433a1d51..f98ddfcf5d0423b4f736f8fb64fce390d01372b9 100644
--- a/arbor/morph/primitives.cpp
+++ b/arbor/morph/primitives.cpp
@@ -1,3 +1,4 @@
+#include <algorithm>
#include <ostream>
#include <arbor/math.hpp>
@@ -9,6 +10,30 @@
namespace arb {
+namespace {
+
+// Advance an iterator to the first value that is not equal to its current
+// value, or end, whichever comes first.
+template <typename T>
+T next_unique(T& it, T end) {
+ const auto& x = *it;
+ ++it;
+ while (it!=end && *it==x) ++it;
+ return it;
+};
+
+// Return the number of times that the value at it is repeated. Advances the
+// iterator to the first value not equal to its current value, or end,
+// whichever comse first.
+template <typename T>
+int multiplicity(T& it, T end) {
+ const auto b = it;
+ return std::distance(b, next_unique(it, end));
+};
+
+} // anonymous namespace
+
+
// interpolate between two points.
mpoint lerp(const mpoint& a, const mpoint& b, double u) {
return { math::lerp(a.x, b.x, u),
@@ -43,6 +68,63 @@ bool test_invariants(const mlocation& l) {
return (0.<=l.pos && l.pos<=1.) && l.branch!=mnpos;
}
+mlocation_list sum(const mlocation_list& lhs, const mlocation_list& rhs) {
+ mlocation_list v;
+ v.resize(lhs.size() + rhs.size());
+ std::merge(lhs.begin(), lhs.end(), rhs.begin(), rhs.end(), v.begin());
+ return v;
+}
+
+mlocation_list join(const mlocation_list& lhs, const mlocation_list& rhs) {
+ mlocation_list L;
+ L.reserve(lhs.size()+rhs.size());
+
+ auto l = lhs.begin();
+ auto lend = lhs.end();
+ auto r = rhs.begin();
+ auto rend = rhs.end();
+
+ auto at_end = [&]() { return l==lend || r==rend; };
+ while (!at_end()) {
+ auto x = (*l<*r) ? *l: *r;
+ auto count = (*l<*r)? multiplicity(l, lend):
+ (*r<*l)? multiplicity(r, rend):
+ std::max(multiplicity(l, lend), multiplicity(r, rend));
+ L.insert(L.end(), count, x);
+ }
+ L.insert(L.end(), l, lend);
+ L.insert(L.end(), r, rend);
+
+ return L;
+}
+
+mlocation_list intersection(const mlocation_list& lhs, const mlocation_list& rhs) {
+ mlocation_list L;
+ L.reserve(lhs.size()+rhs.size());
+
+ auto l = lhs.begin();
+ auto lend = lhs.end();
+ auto r = rhs.begin();
+ auto rend = rhs.end();
+
+ auto at_end = [&]() { return l==lend || r==rend; };
+ while (!at_end()) {
+ if (*l==*r) {
+ auto x = *l;
+ auto count = std::min(multiplicity(l, lend), multiplicity(r, rend));
+ L.insert(L.end(), count, x);
+ }
+ else if (*l<*r) {
+ next_unique(l, lend);
+ }
+ else {
+ next_unique(r, rend);
+ }
+ }
+
+ return L;
+}
+
bool test_invariants(const mcable& c) {
return (0.<=c.prox_pos && c.prox_pos<=c.dist_pos && c.dist_pos<=1.) && c.branch!=mnpos;
}
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index b456c42e19da81d06fa8d615be573d7e65dede88..f0536d6110cfb6fefc38ef6305dd2e9e7245880b 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -2,12 +2,12 @@
#include <string>
#include <vector>
-#include <arbor/morph/error.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/morphexcept.hpp>
+#include <arbor/morph/mprovider.hpp>
#include <arbor/morph/region.hpp>
-#include "morph/em_morphology.hpp"
#include "util/span.hpp"
#include "util/strprintf.hpp"
#include "util/range.hpp"
@@ -16,6 +16,15 @@
namespace arb {
namespace reg {
+// Head and tail of an mcable as mlocations.
+inline mlocation head(mcable c) {
+ return mlocation{c.branch, c.prox_pos};
+}
+
+inline mlocation tail(mcable c) {
+ return mlocation{c.branch, c.dist_pos};
+}
+
// Returns true iff cable sections a and b:
// 1. are on the same branch
// 2. overlap, i.e. their union is not empty.
@@ -57,36 +66,61 @@ mcable_list merge(const mcable_list& v) {
return L;
}
+// List of colocated mlocations, excluding the parameter.
+mlocation_list colocated(mlocation loc, const morphology& m) {
+ mlocation_list L{};
+
+ // Note: if the location is not at the end of a branch, there are no
+ // other colocated points.
+
+ if (loc.pos==0) {
+ // Include head of each branch with same parent,
+ // and end of parent branch if not mnpos.
+
+ auto p = m.branch_parent(loc.branch);
+ if (p!=mnpos) L.push_back({p, 1});
+
+ for (auto b: m.branch_children(p)) {
+ if (b!=loc.branch) L.push_back({b, 0});
+ }
+ }
+ else if (loc.pos==1) {
+ // Include head of each child branch.
+
+ for (auto b: m.branch_children(loc.branch)) {
+ L.push_back({b, 0});
+ }
+ }
+
+ return L;
+}
+
// Insert a zero-length region at the start of each child branch for every cable
// that includes the end of a branch.
// Insert a zero-length region at the end of the parent branch for each cable
// that includes the start of the branch.
-mcable_list cover(mcable_list cables, const em_morphology& m) {
- mcable_list L;
+mcable_list cover(mcable_list cables, const morphology& m) {
+ mcable_list L = cables;
+
for (auto& c: cables) {
- if (c.prox_pos==0) {
- for (auto& x: m.cover(mlocation{c.branch, 0}, false)) {
- L.push_back({x.branch, x.pos, x.pos});
- }
+ for (auto& x: colocated(head(c), m)) {
+ L.push_back({x.branch, x.pos, x.pos});
}
- if (c.dist_pos==1) {
- for (auto& x: m.cover(mlocation{c.branch, 1}, false)) {
- L.push_back({x.branch, x.pos, x.pos});
- }
+ for (auto& x: colocated(tail(c), m)) {
+ L.push_back({x.branch, x.pos, x.pos});
}
}
- L.insert(L.end(), cables.begin(), cables.end());
- util::sort(L);
+ util::sort(L);
return L;
}
-mcable_list remove_cover(mcable_list cables, const em_morphology& m) {
+mcable_list remove_cover(mcable_list cables, const morphology& m) {
// Find all zero-length cables at the end of cables, and convert to
// their canonical representation.
for (auto& c: cables) {
if (c.dist_pos==0 || c.prox_pos==1) {
- auto cloc = m.canonicalize({c.branch, c.prox_pos});
+ auto cloc = canonical(m, head(c));
c = {cloc.branch, cloc.pos, cloc.pos};
}
}
@@ -98,16 +132,16 @@ mcable_list remove_cover(mcable_list cables, const em_morphology& m) {
return merge(cables);
}
-//
-// Null/empty region
-//
+
+// Empty region.
+
struct nil_ {};
region nil() {
return region{nil_{}};
}
-mcable_list thingify_(const nil_& x, const em_morphology& m) {
+mcable_list thingify_(const nil_& x, const mprovider&) {
return {};
}
@@ -115,9 +149,8 @@ std::ostream& operator<<(std::ostream& o, const nil_& x) {
return o << "nil";
}
-//
-// Explicit cable section
-//
+
+// Explicit cable section.
struct cable_ {
mcable cable;
@@ -125,7 +158,7 @@ struct cable_ {
region cable(mcable c) {
if (!test_invariants(c)) {
- throw morphology_error(util::pprintf("Invalid cable section {}", c));
+ throw invalid_mcable(c);
}
return region(cable_{c});
}
@@ -134,13 +167,10 @@ region branch(msize_t bid) {
return cable({bid, 0, 1});
}
-mcable_list thingify_(const cable_& reg, const em_morphology& em) {
- auto& m = em.morph();
-
- if (reg.cable.branch>=m.num_branches()) {
- throw morphology_error(util::pprintf("Branch {} does not exist in morpology", reg.cable.branch));
+mcable_list thingify_(const cable_& reg, const mprovider& p) {
+ if (reg.cable.branch>=p.morphology().num_branches()) {
+ throw no_such_branch(reg.cable.branch);
}
-
return {reg.cable};
}
@@ -148,9 +178,9 @@ std::ostream& operator<<(std::ostream& o, const cable_& c) {
return o << c.cable;
}
-//
-// region with all segments with the same numeric tag
-//
+
+// Region with all segments with the same numeric tag.
+
struct tagged_ {
int tag;
};
@@ -159,8 +189,9 @@ region tagged(int id) {
return region(tagged_{id});
}
-mcable_list thingify_(const tagged_& reg, const em_morphology& em) {
- auto& m = em.morph();
+mcable_list thingify_(const tagged_& reg, const mprovider& p) {
+ const auto& m = p.morphology();
+ const auto& e = p.embedding();
size_t nb = m.num_branches();
std::vector<mcable> L;
@@ -193,8 +224,8 @@ mcable_list thingify_(const tagged_& reg, const em_morphology& em) {
auto first = start-1;
auto last = std::find_if(start, end, not_matches);
- auto l = first==beg? 0.: em.sample2loc(*first).pos;
- auto r = last==end? 1.: em.sample2loc(*(last-1)).pos;
+ auto l = first==beg? 0.: e.sample_location(*first).pos;
+ auto r = last==end? 1.: e.sample_location(*(last-1)).pos;
L.push_back({i, l, r});
// Find the next sample in the branch that matches reg.tag.
@@ -211,17 +242,17 @@ std::ostream& operator<<(std::ostream& o, const tagged_& t) {
return o << "(tag " << t.tag << ")";
}
-//
-// region with all segments in a cell
-//
+
+// Region comprising whole morphology.
+
struct all_ {};
region all() {
return region(all_{});
}
-mcable_list thingify_(const all_&, const em_morphology& m) {
- auto nb = m.morph().num_branches();
+mcable_list thingify_(const all_&, const mprovider& p) {
+ auto nb = p.morphology().num_branches();
mcable_list branches;
branches.reserve(nb);
for (auto i: util::make_span(nb)) {
@@ -234,21 +265,42 @@ std::ostream& operator<<(std::ostream& o, const all_& t) {
return o << "all";
}
-//
-// intersection of two regions.
-//
+
+// Named region.
+
+struct named_ {
+ std::string name;
+};
+
+region named(std::string name) {
+ return region(named_{std::move(name)});
+}
+
+mcable_list thingify_(const named_& n, const mprovider& p) {
+ return p.region(n.name);
+}
+
+std::ostream& operator<<(std::ostream& o, const named_& x) {
+ return o << "(named \"" << x.name << "\")";
+}
+
+
+// Intersection of two regions.
+
struct reg_and {
region lhs;
region rhs;
reg_and(region lhs, region rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {}
};
-mcable_list thingify_(const reg_and& P, const em_morphology& m) {
+mcable_list thingify_(const reg_and& P, const mprovider& p) {
+ auto& m = p.morphology();
+
using cable_it = std::vector<mcable>::const_iterator;
using cable_it_pair = std::pair<cable_it, cable_it>;
- auto lhs = cover(thingify(P.lhs, m), m);
- auto rhs = cover(thingify(P.rhs, m), m);
+ auto lhs = cover(thingify(P.lhs, p), m);
+ auto rhs = cover(thingify(P.rhs, p), m);
// Perform intersection
cable_it_pair it{lhs.begin(), rhs.begin()};
@@ -279,18 +331,18 @@ std::ostream& operator<<(std::ostream& o, const reg_and& x) {
return o << "(intersect " << x.lhs << " " << x.rhs << ")";
}
-//
-// union of two point sets
-//
+
+// Union of two regions.
+
struct reg_or {
region lhs;
region rhs;
reg_or(region lhs, region rhs): lhs(std::move(lhs)), rhs(std::move(rhs)) {}
};
-mcable_list thingify_(const reg_or& P, const em_morphology& m) {
- auto lhs = thingify(P.lhs, m);
- auto rhs = thingify(P.rhs, m);
+mcable_list thingify_(const reg_or& P, const mprovider& p) {
+ auto lhs = thingify(P.lhs, p);
+ auto rhs = thingify(P.rhs, p);
mcable_list L;
L.resize(lhs.size() + rhs.size());
@@ -321,4 +373,12 @@ region::region() {
*this = reg::nil();
}
+region::region(std::string label) {
+ *this = reg::named(std::move(label));
+}
+
+region::region(const char* label) {
+ *this = reg::named(label);
+}
+
} // namespace arb
diff --git a/arbor/morph/sample_tree.cpp b/arbor/morph/sample_tree.cpp
index d8f82c10b33652407a32f7cde86c6331653d2417..15ec4d8cf64cd84fbd44701a7d22409e240e7454 100644
--- a/arbor/morph/sample_tree.cpp
+++ b/arbor/morph/sample_tree.cpp
@@ -1,12 +1,10 @@
-#include <cmath>
+#include <iostream>
+#include <stdexcept>
#include <vector>
-#include <arbor/math.hpp>
-
-#include <arbor/morph/error.hpp>
+#include <arbor/morph/morphexcept.hpp>
#include <arbor/morph/sample_tree.hpp>
-#include "algorithms.hpp"
#include "io/sepval.hpp"
#include "util/span.hpp"
@@ -14,8 +12,7 @@ namespace arb {
sample_tree::sample_tree(std::vector<msample> samples, std::vector<msize_t> parents) {
if (samples.size()!=parents.size()) {
- throw std::runtime_error(
- "The same number of samples and parent indices used to create a sample morphology");
+ throw std::invalid_argument("sample and parent vectors differ in size");
}
reserve(samples.size());
for (auto i: util::make_span(samples.size())) {
@@ -30,14 +27,11 @@ void sample_tree::reserve(msize_t n) {
}
msize_t sample_tree::append(msize_t p, const msample& s) {
- if (empty()) {
- if (p!=mnpos) throw morphology_error("Parent id of root sample must be mnpos");
- }
- else if (p>=size()) {
- throw morphology_error("Parent id of a sample must be less than the sample id");
+ if ((empty() && p!=mnpos) || p>=size()) {
+ if (p!=mnpos) throw invalid_sample_parent(p, size());
}
- auto id = size();
+ auto id = size();
samples_.push_back(s);
parents_.push_back(p);
diff --git a/arbor/profile/meter_manager.cpp b/arbor/profile/meter_manager.cpp
index d3e4097eff06cfc05273a9dbe8e3d66adcb38442..6b17ac44bc647c4877953e393452fa5579074faa 100644
--- a/arbor/profile/meter_manager.cpp
+++ b/arbor/profile/meter_manager.cpp
@@ -171,7 +171,7 @@ std::ostream& operator<<(std::ostream& o, const meter_report& report) {
// Energy measurements are per "per node", so only normalise
// by the number of ranks per node. TODO, this is an
// approximation: better reduce a subset of measurements.
- double energy = algorithms::sum(m.measurements[cp_index])/doms_per_host*1e-3;
+ double energy = util::sum(m.measurements[cp_index])/doms_per_host*1e-3;
sums[m_index] += energy;
o << strprintf("%16.3f", energy);
}
diff --git a/arbor/util/partition.hpp b/arbor/util/partition.hpp
index 0f360fd2bd78bf2cadda452f80b53ea5268a5c1b..e4190b262ebb677336071878d048d4ba91b79b06 100644
--- a/arbor/util/partition.hpp
+++ b/arbor/util/partition.hpp
@@ -5,7 +5,8 @@
#include <stdexcept>
#include <type_traits>
-#include "util/either.hpp"
+#include <arbor/util/either.hpp>
+
#include "util/meta.hpp"
#include "util/partition_iterator.hpp"
#include "util/range.hpp"
diff --git a/arbor/util/piecewise.hpp b/arbor/util/piecewise.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7a11ab029374f2321e6a5e6ee92762011d85877
--- /dev/null
+++ b/arbor/util/piecewise.hpp
@@ -0,0 +1,481 @@
+#pragma once
+
+// Create/manipulate 1-d piece-wise defined objects.
+//
+// Using vectors everywhere here for ease; consider making
+// something more container/sequence-generic later.
+
+#include <initializer_list>
+#include <iterator>
+#include <type_traits>
+#include <vector>
+
+#include "util/meta.hpp"
+#include "util/partition.hpp"
+
+namespace arb {
+namespace util {
+
+using pw_size_type = unsigned;
+constexpr pw_size_type pw_npos = -1;
+
+// Generic random access const iterator for a collection
+// providing operator[].
+
+template <typename T>
+struct indexed_const_iterator {
+ using size_type = decltype(util::size(std::declval<T>()));
+ using difference_type = std::make_signed_t<size_type>;
+
+ using value_type = decltype(std::declval<T>()[0]);
+ struct pointer {
+ value_type v;
+ const value_type* operator->() const { return &v; }
+ };
+
+ using reference = value_type;
+ using iterator_category = std::random_access_iterator_tag;
+
+ const T* ptr_ = nullptr;
+ size_type i_ = 0;
+
+ bool operator==(indexed_const_iterator x) const { return ptr_ == x.ptr_ && i_ == x.i_; }
+ bool operator!=(indexed_const_iterator x) const { return !(*this==x); }
+ bool operator<=(indexed_const_iterator x) const { return i_<=x.i_; }
+ bool operator<(indexed_const_iterator x) const { return i_<x.i_; }
+ bool operator>=(indexed_const_iterator x) const { return i_>=x.i_; }
+ bool operator>(indexed_const_iterator x) const { return i_>x.i_; }
+
+ difference_type operator-(indexed_const_iterator x) const { return i_-x.i_; }
+ indexed_const_iterator& operator++() { return ++i_, *this; }
+ indexed_const_iterator& operator--() { return --i_, *this; }
+ indexed_const_iterator operator++(int) { auto x = *this; return ++i_, x; }
+ indexed_const_iterator operator--(int) { auto x = *this; return --i_, x; }
+
+ indexed_const_iterator operator+(difference_type n) { return indexed_const_iterator{ptr_, i_+n}; }
+ indexed_const_iterator operator-(difference_type n) { return indexed_const_iterator{ptr_, i_-n}; }
+ indexed_const_iterator& operator+=(difference_type n) { return i_+=n, *this; }
+ indexed_const_iterator& operator-=(difference_type n) { return i_-=n, *this; }
+
+ friend indexed_const_iterator operator+(difference_type n, indexed_const_iterator x) {
+ indexed_const_iterator r(std::move(x));
+ return r+=n;
+ }
+
+ reference operator*() const { return (*ptr_)[i_]; }
+ pointer operator->() const { return pointer{(*ptr_)[i_]}; }
+};
+
+
+template <typename X = void>
+struct pw_elements {
+ using size_type = pw_size_type;
+ static constexpr size_type npos = pw_npos;
+
+ using value_type = std::pair<std::pair<double, double>, X>;
+
+ // Consistency requirements:
+ // 1. empty() || element.size()+1 = vertex.size()
+ // 2. vertex[i]<=vertex[j] for i<=j.
+
+ std::vector<double> vertex_;
+ std::vector<X> element_;
+
+ // Ctors and assignment:
+
+ pw_elements() = default;
+
+ template <typename VSeq, typename ESeq>
+ pw_elements(const VSeq& vs, const ESeq& es) {
+ assign(vs, es);
+ }
+
+ pw_elements(std::initializer_list<double> vs, std::initializer_list<X> es) {
+ assign(vs, es);
+ }
+
+ pw_elements(pw_elements&&) = default;
+ pw_elements(const pw_elements&) = default;
+
+ template <typename Y>
+ explicit pw_elements(const pw_elements<Y>& from):
+ vertex_(from.vertex_), element_(from.element_.begin(), from.element_.end())
+ {}
+
+ pw_elements& operator=(pw_elements&&) = default;
+ pw_elements& operator=(const pw_elements&) = default;
+
+ // Access:
+
+ auto intervals() const { return util::partition_view(vertex_); }
+ auto interval(size_type i) const { return intervals()[i]; }
+
+ auto bounds() const { return intervals().bounds(); }
+
+ size_type size() const { return element_.size(); }
+ bool empty() const { return size()==0; }
+
+ bool operator==(const pw_elements& x) const {
+ return vertex_==x.vertex_ && element_==x.element_;
+ }
+
+ bool operator!=(const pw_elements& x) const { return !(*this==x); }
+
+ const auto& elements() const { return element_; }
+ const auto& vertices() const { return vertex_; }
+
+ X& element(size_type i) & { return element_[i]; }
+ const X& element(size_type i) const & { return element_[i]; }
+ value_type operator[](size_type i) const { return value_type{interval(i), element(i)}; }
+
+ using const_iterator = indexed_const_iterator<pw_elements<X>>;
+ using iterator = const_iterator;
+
+ const_iterator cbegin() const { return const_iterator{this, 0}; }
+ const_iterator begin() const { return cbegin(); }
+ const_iterator cend() const { return const_iterator{this, size()}; }
+ const_iterator end() const { return cend(); }
+ value_type front() const { return (*this)[0]; }
+ value_type back() const { return (*this)[size()-1]; }
+
+ size_type index_of(double x) const {
+ if (empty()) return npos;
+
+ auto partn = intervals();
+ if (x == partn.bounds().second) return size()-1;
+ else return partn.index(x);
+ }
+
+ // mutating operations:
+
+ void reserve(size_type n) {
+ vertex_.reserve(n+1);
+ element_.reserve(n);
+ }
+
+ void clear() {
+ vertex_.clear();
+ element_.clear();
+ }
+
+ template <typename U>
+ void push_back(double left, double right, U&& elem) {
+ if (!empty() && left!=vertex_.back()) {
+ throw std::runtime_error("noncontiguous element");
+ }
+
+ if (right<left) {
+ throw std::runtime_error("inverted element");
+ }
+
+ // Extend element_ first in case a conversion/copy/move throws.
+
+ element_.push_back(std::forward<U>(elem));
+ if (vertex_.empty()) vertex_.push_back(left);
+ vertex_.push_back(right);
+ }
+
+ template <typename U>
+ void push_back(double right, U&& elem) {
+ if (empty()) {
+ throw std::runtime_error("require initial left vertex for element");
+ }
+
+ push_back(vertex_.back(), right, elem);
+ }
+
+ void assign(std::initializer_list<double> vs, std::initializer_list<X> es) {
+ using util::make_range;
+ assign(make_range(vs.begin(), vs.end()), make_range(es.begin(), es.end()));
+ }
+
+ template <typename Seq1, typename Seq2>
+ void assign(const Seq1& vertices, const Seq2& elements) {
+ using std::begin;
+ using std::end;
+
+ auto vi = begin(vertices);
+ auto ve = end(vertices);
+
+ auto ei = begin(elements);
+ auto ee = end(elements);
+
+ if (ei==ee) { // empty case
+ if (vi!=ve) {
+ throw std::runtime_error("vertex list too long");
+ }
+ clear();
+ return;
+ }
+
+ double left = *vi++;
+ if (vi==ve) {
+ throw std::runtime_error("vertex list too short");
+ }
+
+ clear();
+
+ double right = *vi++;
+ push_back(left, right, *ei++);
+
+ while (ei!=ee) {
+ if (vi==ve) {
+ throw std::runtime_error("vertex list too short");
+ }
+ double right = *vi++;
+ push_back(right, *ei++);
+ }
+
+ if (vi!=ve) {
+ throw std::runtime_error("vertex list too long");
+ }
+ }
+};
+
+// With X = void, present the element intervals only,
+// keeping othewise the same interface.
+
+template <> struct pw_elements<void> {
+ using size_type = pw_size_type;
+ static constexpr size_type npos = pw_npos;
+
+ std::vector<double> vertex_;
+
+ using value_type = std::pair<double, double>;
+
+ // ctors and assignment:
+
+ template <typename VSeq>
+ explicit pw_elements(const VSeq& vs) { assign(vs); }
+
+ pw_elements(std::initializer_list<double> vs) { assign(vs); }
+
+ pw_elements() = default;
+ pw_elements(pw_elements&&) = default;
+ pw_elements(const pw_elements&) = default;
+
+ template <typename Y>
+ explicit pw_elements(const pw_elements<Y>& from):
+ vertex_(from.vertex_) {}
+
+ pw_elements& operator=(pw_elements&&) = default;
+ pw_elements& operator=(const pw_elements&) = default;
+
+ // access:
+
+ auto intervals() const { return util::partition_view(vertex_); }
+ auto interval(size_type i) const { return intervals()[i]; }
+ value_type operator[](size_type i) const { return interval(i); }
+
+ auto bounds() const { return intervals().bounds(); }
+
+ size_type size() const { return vertex_.empty()? 0: vertex_.size()-1; }
+ bool empty() const { return vertex_.empty(); }
+
+ bool operator==(const pw_elements& x) const { return vertex_==x.vertex_; }
+ bool operator!=(const pw_elements& x) const { return !(*this==x); }
+
+ const auto& vertices() const { return vertex_; }
+
+ size_type index_of(double x) const {
+ if (empty()) return npos;
+
+ auto partn = intervals();
+ if (x == partn.bounds().second) return size()-1;
+ else return partn.index(x);
+ }
+
+ using const_iterator = indexed_const_iterator<pw_elements<void>>;
+ using iterator = const_iterator;
+
+ const_iterator cbegin() const { return const_iterator{this, 0}; }
+ const_iterator begin() const { return cbegin(); }
+ const_iterator cend() const { return const_iterator{this, size()}; }
+ const_iterator end() const { return cend(); }
+ value_type front() const { return (*this)[0]; }
+ value_type back() const { return (*this)[size()-1]; }
+
+ // mutating operations:
+
+ void reserve(size_type n) { vertex_.reserve(n+1); }
+ void clear() { vertex_.clear(); }
+
+ void push_back(double left, double right) {
+ if (!empty() && left!=vertex_.back()) {
+ throw std::runtime_error("noncontiguous element");
+ }
+
+ if (right<left) {
+ throw std::runtime_error("inverted element");
+ }
+
+ if (vertex_.empty()) vertex_.push_back(left);
+ vertex_.push_back(right);
+ }
+
+ void push_back(double right) {
+ if (empty()) {
+ throw std::runtime_error("require initial left vertex for element");
+ }
+ vertex_.push_back(right);
+ }
+
+ void assign(std::initializer_list<double> vs) {
+ assign(util::make_range(vs.begin(), vs.end()));
+ }
+
+ template <typename Seq1>
+ void assign(const Seq1& vertices) {
+ using std::begin;
+ using std::end;
+
+ auto vi = begin(vertices);
+ auto ve = end(vertices);
+
+ if (vi==ve) {
+ clear();
+ return;
+ }
+
+ double left = *vi++;
+ if (vi==ve) {
+ throw std::runtime_error("vertex list too short");
+ }
+
+ clear();
+
+ double right = *vi++;
+ push_back(left, right);
+
+ while (vi!=ve) {
+ double right = *vi++;
+ push_back(right);
+ }
+ }
+};
+
+template <typename X>
+using pw_element = typename pw_elements<X>::value_type;
+
+namespace impl {
+ template <typename A, typename B>
+ struct piecewise_pairify {
+ std::pair<A, B> operator()(
+ double left, double right,
+ const pw_element<A> a_elem,
+ const pw_element<B> b_elem) const
+ {
+ return {a_elem.second, b_elem.second};
+ }
+ };
+
+ template <typename X>
+ struct piecewise_pairify<X, void> {
+ X operator()(
+ double left, double right,
+ const pw_element<X> a_elem,
+ const pw_element<void> b_elem) const
+ {
+ return a_elem.second;
+ }
+ };
+
+ template <typename X>
+ struct piecewise_pairify<void, X> {
+ X operator()(
+ double left, double right,
+ const pw_element<void> a_elem,
+ const pw_element<X> b_elem) const
+ {
+ return b_elem.second;
+ }
+ };
+}
+
+// TODO: Consider making a lazy `zip_view` version of zip.
+
+// Combine functional takes four arguments:
+// double left, double right, pw_elements<A>::value_type, pw_elements<B>::value_type b>
+//
+// Default combine functional returns std::pair<A, B>, unless one of A and B is void.
+
+template <typename A, typename B, typename Combine = impl::piecewise_pairify<A, B>>
+auto zip(const pw_elements<A>& a, const pw_elements<B>& b, Combine combine = {})
+{
+ using Out = decltype(combine(0., 0., a.front(), b.front()));
+ pw_elements<Out> z;
+ if (a.empty() || b.empty()) return z;
+
+ double lmax = std::max(a.bounds().first, b.bounds().first);
+ double rmin = std::min(a.bounds().second, b.bounds().second);
+ if (rmin<lmax) return z;
+
+ double left = lmax;
+ pw_size_type ai = a.intervals().index(left);
+ pw_size_type bi = b.intervals().index(left);
+
+ if (rmin==left) {
+ z.push_back(left, left, combine(left, left, a[ai], b[bi]));
+ return z;
+ }
+
+ double a_right = a.interval(ai).second;
+ double b_right = b.interval(bi).second;
+
+ while (left<rmin) {
+ double right = std::min(a_right, b_right);
+ right = std::min(right, rmin);
+
+ z.push_back(left, right, combine(left, right, a[ai], b[bi]));
+ if (a_right<=right) {
+ a_right = a.interval(++ai).second;
+ }
+ if (b_right<=right) {
+ b_right = b.interval(++bi).second;
+ }
+
+ left = right;
+ }
+ return z;
+}
+
+inline pw_elements<void> zip(const pw_elements<void>& a, const pw_elements<void>& b) {
+ pw_elements<void> z;
+ if (a.empty() || b.empty()) return z;
+
+ double lmax = std::max(a.bounds().first, b.bounds().first);
+ double rmin = std::min(a.bounds().second, b.bounds().second);
+ if (rmin<lmax) return z;
+
+ double left = lmax;
+ pw_size_type ai = a.intervals().index(left);
+ pw_size_type bi = b.intervals().index(left);
+
+ if (rmin==left) {
+ z.push_back(left, left);
+ return z;
+ }
+
+ double a_right = a.interval(ai).second;
+ double b_right = b.interval(bi).second;
+
+ while (left<rmin) {
+ double right = std::min(a_right, b_right);
+ right = std::min(right, rmin);
+
+ z.push_back(left, right);
+ if (a_right<=right) {
+ a_right = a.interval(++ai).second;
+ }
+ if (b_right<=right) {
+ b_right = b.interval(++bi).second;
+ }
+
+ left = right;
+ }
+ return z;
+}
+
+
+} // namespace util
+} // namespace arb
diff --git a/arbor/util/range.hpp b/arbor/util/range.hpp
index 9f175abf98e06e971e8280e28755bf6d8ffe41a7..7eaddb5369e0c14424ff629620d4babf2c114096 100644
--- a/arbor/util/range.hpp
+++ b/arbor/util/range.hpp
@@ -28,7 +28,6 @@
#include <arbor/assert.hpp>
-#include "util/either.hpp"
#include "util/counter.hpp"
#include "util/iterutil.hpp"
#include "util/meta.hpp"
diff --git a/arbor/util/rangeutil.hpp b/arbor/util/rangeutil.hpp
index 25e74a0ca83384cef258c46eef172cdc25a3d240..224d228753858837b083736f65a4a0bb43a8ce5e 100644
--- a/arbor/util/rangeutil.hpp
+++ b/arbor/util/rangeutil.hpp
@@ -193,6 +193,17 @@ bool any_of(const Seq& seq, const Predicate& pred) {
return std::any_of(canon.begin(), canon.end(), pred);
}
+// Accumulate over range
+
+template <
+ typename Seq,
+ typename Value = typename util::sequence_traits<Seq>::value_type
+>
+Value sum(const Seq& seq, Value base = Value{}) {
+ auto canon = canonical_view(seq);
+ return std::accumulate(canon.begin(), canon.end(), base);
+}
+
// Accumulate by projection `proj`
template <
diff --git a/arbor/util/ratelem.hpp b/arbor/util/ratelem.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7a01cada42df3e244f7a1ca726cb8a87e52b9796
--- /dev/null
+++ b/arbor/util/ratelem.hpp
@@ -0,0 +1,173 @@
+#pragma once
+
+// An element representing a segment of a rational polynomial function
+// of order p, q, as determined by its values on n = p+q+1 nodes at
+// [0, 1/n, ..., 1].
+//
+// Rational polynomial interpolation scheme from:
+// F. M. Larkin (1967). Some techniques for rational interpolation.
+// The Computer Journal 10(2), pp. 178–187.
+//
+// TODO: Consider implementing a more generally robust scheme, e.g.
+// S. L. Loi and A. W. McInnes (1983). An algorithm for generalized
+// rational interpolation. BIT Numerical Mathematics 23(1),
+// pp. 105–117. doi:10.1007/BF01937330
+//
+// Current implementation should be sufficient for interpolation of
+// monotonic ratpoly functions of low order, but we can revisit
+// this is if the limitations of the Larkin method bite us.
+
+#include <algorithm>
+#include <array>
+#include <functional>
+#include <type_traits>
+#include <utility>
+
+#include <arbor/math.hpp>
+
+namespace arb {
+namespace util {
+
+namespace impl {
+
+template <unsigned n, unsigned sz>
+struct array_init_n {
+ template <typename A, typename X, typename... Tail>
+ static void set(A& array, X value, Tail... tail) {
+ array[sz-n] = std::move(value);
+ array_init_n<n-1, sz>::set(array, std::forward<Tail>(tail)...);
+ }
+};
+
+template <unsigned sz>
+struct array_init_n<0, sz> {
+ template <typename A>
+ static void set(A& array) {}
+};
+
+// TODO: C++17. The following is much nicer with if constexpr...
+
+template <unsigned a, unsigned c, unsigned k, bool upper>
+struct rat_eval {
+ static double eval(const std::array<double, 1+a+c>& g, double x) {
+ std::array<double, a+c> h;
+ interpolate(h, g, x);
+ return rat_eval<a-1, c, k+1, upper>::eval(h, g, x);
+ }
+
+ static double eval(const std::array<double, 1+a+c>& g, const std::array<double, 2+a+c>&, double x) {
+ std::array<double, a+c> h;
+ interpolate(h, g, x);
+ return rat_eval<a-1, c, k+1, upper>::eval(h, g, x);
+ }
+
+ static void interpolate(std::array<double, a+c>& h, const std::array<double, a+c+1>& g, double x) {
+ constexpr double ook = 1./k;
+ for (unsigned i = 0; i<a+c; ++i) {
+ if (upper) {
+ h[i] = ook*((x - i)*g[i+1] + (i+k - x)*g[i]);
+ }
+ else {
+ // Using h[i] = k/(g[i+1]/(x - i) + g[i]/(i+k - x)) is more robust to
+ // singularities, but the expense should not be necessary if we stick
+ // to strictly monotonic elements for rational polynomials.
+ h[i] = k*g[i]*g[i+1]/(g[i]*(x - i) + g[i+1]*(i+k - x));
+ }
+ }
+ }
+};
+
+template <unsigned k, bool upper>
+struct rat_eval<0, 0, k, upper> {
+ static double eval(const std::array<double, 1>& g, double) {
+ return g[0];
+ }
+
+ static double eval(const std::array<double, 1>& g, const std::array<double, 2>&, double) {
+ return g[0];
+ }
+};
+
+template <unsigned c, unsigned k, bool upper>
+struct rat_eval<0, c, k, upper> {
+ static double eval(const std::array<double, 1+c>& g, const std::array<double, 2+c>& p, double x) {
+ // 'rhombus' interpolation
+ std::array<double, c> h;
+ for (unsigned i = 0; i<c; ++i) {
+ h[i] = p[i+1] + k/((x - i)/(g[i+1]-p[i+1]) + (i+k - x)/(g[i]-p[i+1]));
+ }
+
+ return rat_eval<0, c-1, k+1, upper>::eval(h, g, x);
+ }
+};
+
+} // namespace impl
+
+template <unsigned p, unsigned q>
+struct rat_element {
+ // Construct from function evaluated on nodes.
+ template <
+ typename F,
+ typename _ = std::enable_if_t<std::is_convertible<decltype(std::declval<F>()(0.0)), double>::value>
+ >
+ explicit rat_element(F&& fn, _* = nullptr) {
+ if (size()>1) {
+ for (unsigned i = 0; i<size(); ++i) data_[i] = fn(i/(size()-1.0));
+ }
+ else {
+ data_[0] = fn(0);
+ }
+ }
+
+ // Construct from node values.
+ template <typename... Tail>
+ rat_element(double y0, Tail... tail) {
+ impl::array_init_n<p+q+1, p+q+1>::set(data_, y0, tail...);
+ }
+
+ // Construct from node values in array or std::array.
+ template <typename Y>
+ rat_element(const std::array<Y, 1+p+q>& a) { unchecked_range_init(a); }
+
+ template <typename Y>
+ rat_element(Y (&a)[1+p+q]) { unchecked_range_init(a); }
+
+ // Number of nodes.
+ static constexpr unsigned size() { return 1+p+q; }
+
+ // Rational interpolation at x.
+ double operator()(double x) const {
+ // upper j => interpolate polynomials first;
+ // !upper => interpolate reciprocal polynomials first;
+ // a is the number of (reciprocal) polynomial interpolation steps;
+ // c is the number of 'rhombus' interpolation steps.
+
+ constexpr bool upper = p>=q;
+ constexpr unsigned a = upper? p-q+(q>0): q-p+(p>0);
+ constexpr unsigned c = p+q-a;
+
+ return impl::rat_eval<a, c, 1, upper>::eval(data_, x*(p+q));
+ }
+
+ // Node values.
+ double operator[](unsigned i) const {
+ return data_.at(i);
+ }
+
+ double& operator[](unsigned i) {
+ return data_.at(i);
+ }
+
+private:
+ template <typename C>
+ void unchecked_range_init(const C& a) {
+ using std::begin;
+ using std::end;
+ std::copy(begin(a), end(a), data_.begin());
+ }
+
+ std::array<double, 1+p+q> data_;
+};
+
+} // namespace util
+} // namespace arb
diff --git a/arbor/util/sentinel.hpp b/arbor/util/sentinel.hpp
index 09b985130a584398f4fa285996b368dddaca3258..7a47227e80931b672573efaeb71812b66f15016b 100644
--- a/arbor/util/sentinel.hpp
+++ b/arbor/util/sentinel.hpp
@@ -2,7 +2,8 @@
#include <type_traits>
-#include <util/meta.hpp>
+#include <arbor/util/either.hpp>
+#include "util/meta.hpp"
/*
* Use a proxy iterator to present a range as having the same begin and
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 3e557a08a44667baa0f9861184e02ad6e2bd1da6..8deea24dcce1ab01983d8a176b770427ad338faa 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -89,7 +89,6 @@ set(unit_sources
test_dry_run_context.cpp
test_double_buffer.cpp
test_either.cpp
- test_em_morphology.cpp
test_event_binner.cpp
test_event_delivery.cpp
test_event_generators.cpp
@@ -113,15 +112,20 @@ set(unit_sources
test_mechinfo.cpp
test_merge_events.cpp
test_morphology.cpp
+ test_morph_embedding.cpp
+ test_morph_expr.cpp
+ test_morph_primitives.cpp
test_multi_event_stream.cpp
test_optional.cpp
test_padded.cpp
test_partition.cpp
test_partition_by_constraint.cpp
test_path.cpp
+ test_piecewise.cpp
test_point.cpp
test_probe.cpp
test_range.cpp
+ test_ratelem.cpp
test_segment.cpp
test_schedule.cpp
test_spike_source.cpp
diff --git a/test/unit/test_algorithms.cpp b/test/unit/test_algorithms.cpp
index e77e3f9978007d18d3189d0b2f31fd227f18fdba..ab975e0cef8ab41a292856e02875dc6bd332837b 100644
--- a/test/unit/test_algorithms.cpp
+++ b/test/unit/test_algorithms.cpp
@@ -9,28 +9,13 @@
#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, sum)
-{
- // sum of 10 times 2 is 20
- std::vector<int> v1(10, 2);
- EXPECT_EQ(10*2, arb::algorithms::sum(v1));
-
- // make an array 1:20 and sum it up using formula for arithmetic sequence
- auto n = 20;
- std::vector<int> v2(n);
- // can't use iota because the Intel compiler optimizes it out, despite
- // the result being required in EXPECT_EQ
- // std::iota(v2.begin(), v2.end(), 1);
- for(auto i=0; i<n; ++i) { v2[i] = i+1; }
- EXPECT_EQ((n+1)*n/2, arb::algorithms::sum(v2));
-}
-
TEST(algorithms, make_index)
{
{
@@ -39,7 +24,7 @@ TEST(algorithms, make_index)
EXPECT_EQ(index.size(), 11u);
EXPECT_EQ(index.front(), 0);
- EXPECT_EQ(index.back(), arb::algorithms::sum(v));
+ EXPECT_EQ(index.back(), arb::util::sum(v));
}
{
@@ -49,7 +34,7 @@ TEST(algorithms, make_index)
EXPECT_EQ(index.size(), 11u);
EXPECT_EQ(index.front(), 0);
- EXPECT_EQ(index.back(), arb::algorithms::sum(v));
+ EXPECT_EQ(index.back(), arb::util::sum(v));
}
}
diff --git a/test/unit/test_compartments.cpp b/test/unit/test_compartments.cpp
index a86ef3e65f0850894e1af0ed7134abcbc29bb4a0..deb5996ee18d806e63a4d26fa9d6552af3fbc136 100644
--- a/test/unit/test_compartments.cpp
+++ b/test/unit/test_compartments.cpp
@@ -5,8 +5,8 @@
#include <arbor/math.hpp>
-#include "algorithms.hpp"
#include "fvm_compartment.hpp"
+#include "util/rangeutil.hpp"
#include "util/span.hpp"
#include "util/transform.hpp"
@@ -15,7 +15,7 @@ using namespace arb::math;
using arb::util::make_span;
using arb::util::transform_view;
-using arb::algorithms::sum;
+using arb::util::sum;
// Divided compartments
// (FVM-friendly compartment data)
diff --git a/test/unit/test_cv_policy.cpp b/test/unit/test_cv_policy.cpp
index 273a413ef8423372ff8a63a9987bc616de4a67db..aa11e16782001a8dc4bfcbc496f26a1ab366a2f1 100644
--- a/test/unit/test_cv_policy.cpp
+++ b/test/unit/test_cv_policy.cpp
@@ -8,8 +8,8 @@
#include <arbor/cable_cell_param.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/locset.hpp>
+#include <arbor/morph/mprovider.hpp>
-#include "morph/em_morphology.hpp"
#include "util/filter.hpp"
#include "util/rangeutil.hpp"
#include "util/span.hpp"
@@ -72,7 +72,7 @@ TEST(cv_policy, explicit_policy) {
cable_cell cell(m);
locset result = pol.cv_boundary_points(cell);
- EXPECT_EQ(thingify(lset, *cell.morphology()), thingify(result, *cell.morphology()));
+ EXPECT_EQ(thingify(lset, cell.provider()), thingify(result, cell.provider()));
}
}
@@ -90,10 +90,10 @@ TEST(cv_policy, empty_morphology) {
};
cable_cell cell(m_empty);
- auto empty_loclist = thingify(ls::nil(), *cell.morphology());
+ auto empty_loclist = thingify(ls::nil(), cell.provider());
for (auto& pol: policies) {
- EXPECT_EQ(empty_loclist, thingify(pol.cv_boundary_points(cell), *cell.morphology()));
+ EXPECT_EQ(empty_loclist, thingify(pol.cv_boundary_points(cell), cell.provider()));
}
}
@@ -115,8 +115,8 @@ TEST(cv_policy, single_root_cv) {
cable_cell cell(morph);
for (auto& polpair: policy_pairs) {
- mlocation_list p1 = thingify(polpair.first.cv_boundary_points(cell), *cell.morphology());
- mlocation_list p2 = thingify(polpair.second.cv_boundary_points(cell), *cell.morphology());
+ mlocation_list p1 = thingify(polpair.first.cv_boundary_points(cell), cell.provider());
+ mlocation_list p2 = thingify(polpair.second.cv_boundary_points(cell), cell.provider());
auto p1_no_b0 = util::filter(p1, [](mlocation l) { return l.branch>0; });
mlocation_list expected = {{0,0}, {0,1}};
@@ -139,14 +139,14 @@ TEST(cv_policy, fixed_per_branch) {
cv_policy pol = cv_policy_fixed_per_branch(4);
locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(L{0, 0}, L{0, 0.25}, L{0, 0.5}, L{0, 0.75}, L{0, 1});
- EXPECT_EQ(thingify(expected, *cell.morphology()), thingify(points, *cell.morphology()));
+ EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
}
{
// interior fork points
cv_policy pol = cv_policy_fixed_per_branch(4, interior_forks);
locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(L{0, 0.125}, L{0, 0.375}, L{0, 0.625}, L{0, 0.875});
- EXPECT_EQ(thingify(expected, *cell.morphology()), thingify(points, *cell.morphology()));
+ EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
}
}
@@ -163,7 +163,7 @@ TEST(cv_policy, fixed_per_branch) {
L{0, 0}, L{0, 0.5}, L{0,1}, L{1, 0}, L{1, 0.5}, L{1,1}, L{2, 0}, L{2, 0.5}, L{2,1},
L{3, 0}, L{3, 0.5}, L{3,1}, L{4, 0}, L{4, 0.5}, L{4,1}, L{5, 0}, L{5, 0.5}, L{5,1}
);
- EXPECT_EQ(thingify(expected, *cell.morphology()), thingify(points, *cell.morphology()));
+ EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
}
{
// interior fork points
@@ -173,7 +173,7 @@ TEST(cv_policy, fixed_per_branch) {
L{0, 0.25}, L{0, 0.75}, L{1, 0.25}, L{1, 0.75}, L{2, 0.25}, L{2, 0.75},
L{3, 0.25}, L{3, 0.75}, L{4, 0.25}, L{4, 0.75}, L{5, 0.25}, L{5, 0.75}
);
- EXPECT_EQ(thingify(expected, *cell.morphology()), thingify(points, *cell.morphology()));
+ EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
}
}
}
@@ -185,7 +185,7 @@ TEST(cv_policy, max_extent) {
// root branch only
for (auto& morph: {m_sph_b1, m_reg_b1}) {
cable_cell cell(morph);
- ASSERT_EQ(1.0, cell.morphology()->branch_length(0));
+ ASSERT_EQ(1.0, cell.embedding().branch_length(0));
{
// extent of 0.25 should give exact fp calculation, giving
@@ -193,31 +193,31 @@ TEST(cv_policy, max_extent) {
cv_policy pol = cv_policy_max_extent(0.25);
locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(L{0, 0}, L{0, 0.25}, L{0, 0.5}, L{0, 0.75}, L{0, 1});
- EXPECT_EQ(thingify(expected, *cell.morphology()), thingify(points, *cell.morphology()));
+ EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
}
{
cv_policy pol = cv_policy_max_extent(0.25, interior_forks);
locset points = pol.cv_boundary_points(cell);
locset expected = as_locset(L{0, 0.125}, L{0, 0.375}, L{0, 0.625}, L{0, 0.875});
- EXPECT_EQ(thingify(expected, *cell.morphology()), thingify(points, *cell.morphology()));
+ EXPECT_EQ(thingify(expected, cell.provider()), thingify(points, cell.provider()));
}
}
// cell with varying branch lengths; extent not exact fraction.
{
cable_cell cell(m_mlt_b6);
- ASSERT_EQ(1.0, cell.morphology()->branch_length(0));
- ASSERT_EQ(1.0, cell.morphology()->branch_length(1));
- ASSERT_EQ(2.0, cell.morphology()->branch_length(2));
- ASSERT_EQ(4.0, cell.morphology()->branch_length(3));
- ASSERT_EQ(1.0, cell.morphology()->branch_length(4));
- ASSERT_EQ(2.0, cell.morphology()->branch_length(5));
+ ASSERT_EQ(1.0, cell.embedding().branch_length(0));
+ ASSERT_EQ(1.0, cell.embedding().branch_length(1));
+ ASSERT_EQ(2.0, cell.embedding().branch_length(2));
+ ASSERT_EQ(4.0, cell.embedding().branch_length(3));
+ ASSERT_EQ(1.0, cell.embedding().branch_length(4));
+ ASSERT_EQ(2.0, cell.embedding().branch_length(5));
{
// max extent of 0.6 should give two CVs on branches of length 1,
// four CVs on branches of length 2, and seven CVs on the branch of length 4.
cv_policy pol = cv_policy_max_extent(0.6);
- mlocation_list points = thingify(pol.cv_boundary_points(cell), *cell.morphology());
+ mlocation_list points = thingify(pol.cv_boundary_points(cell), cell.provider());
mlocation_list points_b012 = util::assign_from(util::filter(points, [](mlocation l) { return l.branch<3; }));
mlocation_list expected_b012 = {
diff --git a/test/unit/test_either.cpp b/test/unit/test_either.cpp
index 79344e211edccf0fdec3d9585245f6d7a55e10c3..4d18d47c60f1b126c5a09cca7d52a723cf4ed47b 100644
--- a/test/unit/test_either.cpp
+++ b/test/unit/test_either.cpp
@@ -2,8 +2,9 @@
#include <array>
#include <algorithm>
+#include <arbor/util/either.hpp>
+
#include "../gtest.h"
-#include "util/either.hpp"
// TODO: coverage!
diff --git a/test/unit/test_em_morphology.cpp b/test/unit/test_em_morphology.cpp
deleted file mode 100644
index 36cbf05f5b1e0a686445c70db71cd47d3097f70a..0000000000000000000000000000000000000000
--- a/test/unit/test_em_morphology.cpp
+++ /dev/null
@@ -1,694 +0,0 @@
-/*
- * Unit tests for em_morphology, region, locset, label_dict.
- */
-
-#include <cmath>
-#include <string>
-#include <vector>
-
-#include "../test/gtest.h"
-
-#include <arbor/math.hpp>
-#include <arbor/morph/error.hpp>
-#include <arbor/morph/morphology.hpp>
-#include <arbor/morph/sample_tree.hpp>
-#include <arbor/morph/locset.hpp>
-#include <arbor/morph/region.hpp>
-
-#include "io/sepval.hpp"
-#include "morph/em_morphology.hpp"
-#include "util/span.hpp"
-
-template <typename T>
-std::ostream& operator<<(std::ostream& o, const std::vector<T>& v) {
- return o << "[" << arb::io::csv(v) << "]";
-}
-
-// Test the morphology meta-data that is cached on construction of
-// em_morpholgy, e.g. interpolation information and terminal nodes.
-TEST(em_morphology, cache) {
- using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
- using loc = arb::mlocation;
- constexpr auto npos = arb::mnpos;
-
- // A single unbranched cable with 5 sample points.
- // The cable has length 10 μm, with samples located at
- // 0 μm, 1 μm, 3 μm, 7 μm and 10 μm.
- {
- pvec parents = {npos, 0, 1, 2, 3};
- svec samples = {
- {{ 0, 0, 0, 2}, 1},
- {{ 1, 0, 0, 2}, 1},
- {{ 3, 0, 0, 2}, 1},
- {{ 7, 0, 0, 2}, 1},
- {{ 10, 0, 0, 2}, 1},
- };
- arb::sample_tree sm(samples, parents);
- arb::em_morphology em(arb::morphology(sm, false));
- EXPECT_EQ(em.sample2loc(0), (loc{0,0}));
- EXPECT_EQ(em.sample2loc(1), (loc{0,0.1}));
- EXPECT_EQ(em.sample2loc(2), (loc{0,0.3}));
- EXPECT_EQ(em.sample2loc(3), (loc{0,0.7}));
- EXPECT_EQ(em.sample2loc(4), (loc{0,1}));
-
- EXPECT_EQ(em.root(), (loc{0,0}));
- EXPECT_EQ(em.terminals(), (arb::mlocation_list{{0,1}}));
-
- EXPECT_EQ(10., em.branch_length(0));
- }
-
- // Eight samples
- //
- // sample ids:
- // 0
- // 1 3
- // 2 4
- // 5 6
- // 7
- { // Spherical root.
- pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
-
- svec samples = {
- {{ 0, 0, 0, 10}, 1},
- {{ 10, 0, 0, 2}, 3},
- {{100, 0, 0, 2}, 3},
- {{ 0, 10, 0, 2}, 3},
- {{ 0,100, 0, 2}, 3},
- {{100,100, 0, 2}, 3},
- {{ 0,200, 0, 2}, 3},
- {{ 0,300, 0, 2}, 3},
- };
- arb::sample_tree sm(samples, parents);
- arb::em_morphology em(arb::morphology(sm, true));
-
- EXPECT_EQ(em.sample2loc(0), (loc{0,0}));
- EXPECT_EQ(em.sample2loc(1), (loc{1,0}));
- EXPECT_EQ(em.sample2loc(2), (loc{1,1}));
- EXPECT_EQ(em.sample2loc(3), (loc{2,0}));
- EXPECT_EQ(em.sample2loc(4), (loc{2,1}));
- EXPECT_EQ(em.sample2loc(5), (loc{3,1}));
- EXPECT_EQ(em.sample2loc(6), (loc{4,0.5}));
- EXPECT_EQ(em.sample2loc(7), (loc{4,1}));
-
- EXPECT_EQ(em.root(), (loc{0,0}));
- EXPECT_EQ(em.terminals(), (arb::mlocation_list{{1,1}, {3,1}, {4,1}}));
-
- ASSERT_EQ(5u, em.morph().num_branches());
- EXPECT_EQ(20., em.branch_length(0));
- EXPECT_EQ(90., em.branch_length(1));
- EXPECT_EQ(90., em.branch_length(2));
- EXPECT_EQ(100., em.branch_length(3));
- EXPECT_EQ(200., em.branch_length(4));
- }
- { // No Spherical root
- pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
-
- svec samples = {
- {{ 0, 0, 0, 2}, 1},
- {{ 10, 0, 0, 2}, 3},
- {{100, 0, 0, 2}, 3},
- {{ 0, 10, 0, 2}, 3},
- {{ 0,100, 0, 2}, 3},
- {{100,100, 0, 2}, 3},
- {{ 0,130, 0, 2}, 3},
- {{ 0,300, 0, 2}, 3},
- };
- arb::sample_tree sm(samples, parents);
- arb::em_morphology em(arb::morphology(sm, false));
-
- EXPECT_EQ(em.sample2loc(0), (loc{0,0}));
- EXPECT_EQ(em.sample2loc(1), (loc{0,0.1}));
- EXPECT_EQ(em.sample2loc(2), (loc{0,1}));
- EXPECT_EQ(em.sample2loc(3), (loc{1,0.1}));
- EXPECT_EQ(em.sample2loc(4), (loc{1,1}));
- EXPECT_EQ(em.sample2loc(5), (loc{2,1}));
- EXPECT_EQ(em.sample2loc(6), (loc{3,0.15}));
- EXPECT_EQ(em.sample2loc(7), (loc{3,1}));
-
- EXPECT_EQ(em.root(), (loc{0,0}));
- EXPECT_EQ(em.terminals(), (arb::mlocation_list{{0,1}, {2,1}, {3,1}}));
-
- ASSERT_EQ(4u, em.morph().num_branches());
- EXPECT_EQ(100., em.branch_length(0));
- EXPECT_EQ(100., em.branch_length(1));
- EXPECT_EQ(100., em.branch_length(2));
- EXPECT_EQ(200., em.branch_length(3));
- }
-}
-
-TEST(em_morphology, cover) {
- using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
- using ll = arb::mlocation_list;
- constexpr auto npos = arb::mnpos;
- using arb::util::make_span;
-
- // Eight samples
- // 0
- // 1 3
- // 2 4
- // 5 6
- // 7
- pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
-
- svec samples = {
- {{ 0, 0, 0, 10}, 1},
- {{ 10, 0, 0, 2}, 3},
- {{100, 0, 0, 2}, 3},
- {{ 0, 10, 0, 2}, 3},
- {{ 0,100, 0, 2}, 3},
- {{100,100, 0, 2}, 3},
- {{ 0,200, 0, 2}, 3},
- {{ 0,300, 0, 2}, 3},
- };
- arb::sample_tree sm(samples, parents);
-
- // non-spherical root.
- {
- arb::em_morphology em(arb::morphology(sm, false));
-
- for (arb::msize_t i: make_span(em.morph().num_branches())) {
- for (auto j: make_span(1,9)) {
- arb::mlocation l{i, j/10.};
- EXPECT_EQ(em.cover(l), (ll{l}));
- }
- }
-
- EXPECT_EQ(em.cover({0,0}), (ll{{0,0}, {1,0}}));
- EXPECT_EQ(em.cover({0,1}), (ll{{0,1}}));
- EXPECT_EQ(em.cover({1,1}), (ll{{1,1}, {2,0}, {3,0}}));
- EXPECT_EQ(em.cover({2,1}), (ll{{2,1}}));
- EXPECT_EQ(em.cover({3,1}), (ll{{3,1}}));
-
- EXPECT_EQ(em.cover({1,0}), (ll{{0,0}, {1,0}}));
- EXPECT_EQ(em.cover({2,0}), (ll{{1,1}, {2,0}, {3,0}}));
- EXPECT_EQ(em.cover({3,0}), (ll{{1,1}, {2,0}, {3,0}}));
- }
- // spherical root.
- {
- arb::em_morphology em(arb::morphology(sm, true));
-
- EXPECT_EQ(em.cover({0,0}), (ll{{0,0}}));
- EXPECT_EQ(em.cover({0,1}), (ll{{0,1}, {1,0}, {2,0}}));
- EXPECT_EQ(em.cover({1,1}), (ll{{1,1}}));
- EXPECT_EQ(em.cover({2,1}), (ll{{2,1}, {3,0}, {4,0}}));
- EXPECT_EQ(em.cover({3,1}), (ll{{3,1}}));
- EXPECT_EQ(em.cover({4,1}), (ll{{4,1}}));
- EXPECT_EQ(em.cover({1,0}), (ll{{0,1}, {1,0}, {2,0}}));
- EXPECT_EQ(em.cover({2,0}), (ll{{0,1}, {1,0}, {2,0}}));
- EXPECT_EQ(em.cover({3,0}), (ll{{2,1}, {3,0}, {4,0}}));
- EXPECT_EQ(em.cover({4,0}), (ll{{2,1}, {3,0}, {4,0}}));
- }
-}
-
-TEST(em_morphology, minset) {
- using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
- using ll = arb::mlocation_list;
- constexpr auto npos = arb::mnpos;
-
- // Eight samples
- // no-sphere sphere
- // sample branch branch
- // 0 0 0
- // 1 3 0 1 1 2
- // 2 4 0 1 1 2
- // 5 6 2 3 3 4
- // 7 3 4
- pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
- svec samples = {
- {{ 0, 0, 0, 2}, 3},
- {{ 10, 0, 0, 2}, 3},
- {{100, 0, 0, 2}, 3},
- {{ 0, 10, 0, 2}, 3},
- {{ 0,100, 0, 2}, 3},
- {{100,100, 0, 2}, 3},
- {{ 0,200, 0, 2}, 3},
- {{ 0,300, 0, 2}, 3},
- };
- arb::sample_tree sm(samples, parents);
- {
- arb::em_morphology em(arb::morphology(sm, false));
-
- EXPECT_EQ((ll{}), em.minset(ll{}));
- EXPECT_EQ((ll{{2,0.1}}), em.minset(ll{{2,0.1}}));
- EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), em.minset(ll{{0,0.5}, {1,0.5}}));
- EXPECT_EQ((ll{{0,0.5}}), em.minset(ll{{0,0.5}}));
- EXPECT_EQ((ll{{0,0}, {1,0}}), em.minset(ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
- EXPECT_EQ((ll{{0,0}, {1,0.5}}), em.minset(ll{{0,0}, {0,0.5}, {1,0.5}, {2,0.5}}));
- EXPECT_EQ((ll{{0,0}, {2,0.5}}), em.minset(ll{{0,0}, {0,0.5}, {2,0.5}}));
- EXPECT_EQ((ll{{0,0}, {2,0.5}, {3,0}}), em.minset(ll{{0,0}, {0,0.5}, {2,0.5}, {3,0}, {3,1}}));
- EXPECT_EQ((ll{{0,0}, {1,0}}), em.minset(ll{{0,0}, {0,0.5}, {1,0}, {2,0.5}, {3,0}, {3,1}}));
- }
- {
- arb::em_morphology em(arb::morphology(sm, true));
-
- EXPECT_EQ((ll{}), em.minset(ll{}));
- EXPECT_EQ((ll{{2,0.1}}), em.minset(ll{{2,0.1}}));
- EXPECT_EQ((ll{{0,0.5}}), em.minset(ll{{0,0.5}, {1,0.5}}));
- EXPECT_EQ((ll{{0,0.5}}), em.minset(ll{{0,0.5}}));
- EXPECT_EQ((ll{{0,0}}), em.minset(ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
- EXPECT_EQ((ll{{1,0.5}, {3,0.1}, {4,0.5}}), em.minset(ll{{1,0.5}, {1,1}, {3,0.1}, {4,0.5}, {4,0.7}}));
- }
-}
-
-// Test expressions that describe location sets (locset).
-TEST(locset, expressions) {
- auto root = arb::ls::root();
- auto term = arb::ls::terminal();
- auto samp = arb::ls::sample(42);
- auto loc = arb::ls::location({2, 0.5});
-
- auto to_string = [](auto&& x) {
- std::stringstream s;
- s << x;
- return s.str();
- };
-
- EXPECT_EQ(to_string(root), "root");
- EXPECT_EQ(to_string(term), "terminal");
- EXPECT_EQ(to_string(sum(root, term)), "(sum root terminal)");
- EXPECT_EQ(to_string(sum(root, term, samp)),
- "(sum (sum root terminal) (sample 42))");
- EXPECT_EQ(to_string(sum(root, term, samp, loc)),
- "(sum (sum (sum root terminal) (sample 42)) (location 2 0.5))");
- EXPECT_EQ(to_string(samp), "(sample 42)");
- EXPECT_EQ(to_string(loc), "(location 2 0.5)");
-
- // Location positions have to be in the range [0,1].
- // Assert that an exception is thrown if and invalide location is requested.
- EXPECT_THROW(arb::ls::location({2, 1.5}), arb::morphology_error);
- EXPECT_THROW(arb::ls::location({arb::mnpos, 0.5}), arb::morphology_error);
-}
-
-// Test expressions that describe regions.
-TEST(region, expressions) {
- using arb::reg::cable;
- auto to_string = [](auto&& x) {
- std::stringstream s;
- s << x;
- return s.str();
- };
-
- auto c1 = arb::reg::cable({1,0,1});
- auto c2 = arb::reg::cable({4,0.1,0.5});
- auto c3 = join(cable({4,0.1,0.5}), cable({3,0,1}));
- auto b1 = arb::reg::branch(1);
- auto t1 = arb::reg::tagged(1);
- auto t2 = arb::reg::tagged(2);
- auto t3 = arb::reg::tagged(3);
- auto all = arb::reg::all();
-
- EXPECT_EQ(to_string(c1), "(cable 1 0 1)");
- EXPECT_EQ(to_string(c2), "(cable 4 0.1 0.5)");
- EXPECT_EQ(to_string(c3), "(join (cable 4 0.1 0.5) (cable 3 0 1))");
- EXPECT_EQ(to_string(b1), "(cable 1 0 1)");
- EXPECT_EQ(to_string(t1), "(tag 1)");
- EXPECT_EQ(to_string(t2), "(tag 2)");
- EXPECT_EQ(to_string(intersect(c1, t2)), "(intersect (cable 1 0 1) (tag 2))");
- EXPECT_EQ(to_string(join(c1, t2)), "(join (cable 1 0 1) (tag 2))");
- EXPECT_EQ(to_string(join(t1, t2, t3)), "(join (join (tag 1) (tag 2)) (tag 3))");
- EXPECT_EQ(to_string(intersect(t1, t2, t3)), "(intersect (intersect (tag 1) (tag 2)) (tag 3))");
- EXPECT_EQ(to_string(intersect(join(c1, t2), c2)), "(intersect (join (cable 1 0 1) (tag 2)) (cable 4 0.1 0.5))");
- EXPECT_EQ(to_string(all), "all");
-
- EXPECT_THROW(arb::reg::cable({1, 0, 1.1}), arb::morphology_error);
- EXPECT_THROW(arb::reg::branch(-1), arb::morphology_error);
-}
-
-TEST(locset, thingify) {
- using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
- using ll = arb::mlocation_list;
- const auto npos = arb::mnpos;
-
- auto root = arb::ls::root();
- auto term = arb::ls::terminal();
- auto samp = arb::ls::sample(4);
- auto midb2 = arb::ls::location({2, 0.5});
- auto midb1 = arb::ls::location({1, 0.5});
- auto begb0 = arb::ls::location({0, 0});
- auto begb1 = arb::ls::location({1, 0});
- auto begb2 = arb::ls::location({2, 0});
- auto begb3 = arb::ls::location({3, 0});
- auto begb4 = arb::ls::location({4, 0});
-
- // Eight samples
- //
- // 0
- // 1 3
- // 2 4
- // 5 6
- // 7
- pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
- svec samples = {
- {{ 0, 0, 0, 2}, 3},
- {{ 10, 0, 0, 2}, 3},
- {{100, 0, 0, 2}, 3},
- {{ 0, 10, 0, 2}, 3},
- {{ 0,100, 0, 2}, 3},
- {{100,100, 0, 2}, 3},
- {{ 0,200, 0, 2}, 3},
- {{ 0,300, 0, 2}, 3},
- };
- arb::sample_tree sm(samples, parents);
-
- {
- arb::em_morphology em(arb::morphology(sm, true));
-
- EXPECT_EQ(thingify(root, em), (ll{{0,0}}));
- EXPECT_EQ(thingify(term, em), (ll{{1,1},{3,1},{4,1}}));
- EXPECT_EQ(thingify(samp, em), (ll{{2,1}}));
- EXPECT_EQ(thingify(midb2, em), (ll{{2,0.5}}));
- EXPECT_EQ(thingify(midb1, em), (ll{{1,0.5}}));
- EXPECT_EQ(thingify(begb0, em), (ll{{0,0}}));
- EXPECT_EQ(thingify(begb1, em), (ll{{1,0}}));
- EXPECT_EQ(thingify(begb2, em), (ll{{2,0}}));
- EXPECT_EQ(thingify(begb3, em), (ll{{3,0}}));
- EXPECT_EQ(thingify(begb4, em), (ll{{4,0}}));
- }
- {
- arb::em_morphology em(arb::morphology(sm, false));
-
- EXPECT_EQ(thingify(root, em), (ll{{0,0}}));
- EXPECT_EQ(thingify(term, em), (ll{{0,1},{2,1},{3,1}}));
- EXPECT_EQ(thingify(samp, em), (ll{{1,1}}));
- EXPECT_EQ(thingify(midb2, em), (ll{{2,0.5}}));
- EXPECT_EQ(thingify(midb1, em), (ll{{1,0.5}}));
- EXPECT_EQ(thingify(begb0, em), (ll{{0,0}}));
- EXPECT_EQ(thingify(begb1, em), (ll{{1,0}}));
- EXPECT_EQ(thingify(begb2, em), (ll{{2,0}}));
- EXPECT_EQ(thingify(begb3, em), (ll{{3,0}}));
- // In the absence of a spherical root, there is no branch 4.
- EXPECT_THROW(thingify(begb4, em), arb::morphology_error);
- }
-}
-
-// Forward declare implementation of join, union & intersect on location lists
-// for testing.
-namespace arb { namespace ls {
- mlocation_list intersection(const mlocation_list& lhs, const mlocation_list& rhs);
- mlocation_list join(const mlocation_list& lhs, const mlocation_list& rhs);
- mlocation_list sum(const mlocation_list& lhs, const mlocation_list& rhs);
-}}
-
-arb::mlocation_list ml(std::vector<int> bids) {
- arb::mlocation_list L;
- for (arb::msize_t b: bids) L.push_back({b, 0});
- return L;
-}
-
-
-// Test out multiset joins, intersections and sums
-TEST(locset, join_intersect_sum) {
- using ll = arb::mlocation_list;
- using arb::ls::intersection;
- using arb::ls::sum;
- using arb::ls::join;
-
- {
- ll lhs{};
- ll rhs{};
- EXPECT_EQ(sum(lhs, rhs), ll{});
- EXPECT_EQ(join(lhs, rhs), ll{});
- EXPECT_EQ(intersection(lhs, rhs), ll{});
- }
- {
- ll lhs{};
- ll rhs = ml({0,1});
- EXPECT_EQ(sum(lhs, rhs), rhs);
- EXPECT_EQ(join(lhs, rhs), rhs);
- EXPECT_EQ(intersection(lhs, rhs), ll{});
- }
- {
- ll lhs = ml({1});
- ll rhs = ml({1});
- EXPECT_EQ(sum(lhs, rhs), ml({1,1}));
- EXPECT_EQ(join(lhs, rhs), ml({1}));
- EXPECT_EQ(intersection(lhs, rhs), ml({1}));
- }
- {
- ll lhs = ml({1,1});
- ll rhs = ml({1});
- EXPECT_EQ(sum(lhs, rhs), ml({1,1,1}));
- EXPECT_EQ(join(lhs, rhs), ml({1,1}));
- EXPECT_EQ(intersection(lhs, rhs), ml({1}));
- }
- {
- ll lhs = ml({0,3});
- ll rhs = ml({1,2});
- EXPECT_EQ(sum(lhs, rhs), ml({0,1,2,3}));
- EXPECT_EQ(join(lhs, rhs), ml({0,1,2,3}));
- EXPECT_EQ(intersection(lhs, rhs), ll{});
- }
- {
- ll lhs = ml({0,1,3});
- ll rhs = ml({0,1,3});
- EXPECT_EQ(sum(lhs, rhs), ml({0,0,1,1,3,3}));
- EXPECT_EQ(join(lhs, rhs), lhs);
- EXPECT_EQ(intersection(lhs, rhs), lhs);
- }
- {
- ll lhs = ml({0,1,3});
- ll rhs = ml({1,2});
- EXPECT_EQ(sum(lhs, rhs), ml({0,1,1,2,3}));
- EXPECT_EQ(join(lhs, rhs), ml({0,1,2,3}));
- EXPECT_EQ(intersection(lhs, rhs), ml({1}));
- }
- {
- ll lhs = ml({0,1,1,3});
- ll rhs = ml({1,2});
- EXPECT_EQ(sum(lhs, rhs), ml({0,1,1,1,2,3}));
- EXPECT_EQ(join(lhs, rhs), ml({0,1,1,2,3}));
- EXPECT_EQ(intersection(lhs, rhs), ml({1}));
- }
- {
- ll lhs = ml({0,1,1,3,5,5,12});
- ll rhs = ml({1,2,2,5,5,5});
- EXPECT_EQ(sum(lhs, rhs), ml({0,1,1,1,2,2,3,5,5,5,5,5,12}));
- EXPECT_EQ(join(lhs, rhs), ml({0,1,1,2,2,3,5,5,5,12}));
- EXPECT_EQ(intersection(lhs, rhs), ml({1,5,5}));
- }
-}
-
-TEST(region, thingify) {
- using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
- //using cab = arb::mcable;
- using cl = arb::mcable_list;
- constexpr auto npos = arb::mnpos;
-
- // A single unbranched cable with 5 sample points.
- // The cable has length 10 μm, with samples located at
- // 0 μm, 1 μm, 3 μm, 7 μm and 10 μm.
- {
- pvec parents = {npos, 0, 1, 2, 3};
- svec samples = {
- {{ 0, 0, 0, 2}, 1},
- {{ 1, 0, 0, 2}, 1},
- {{ 3, 0, 0, 2}, 2},
- {{ 7, 0, 0, 2}, 1},
- {{ 10, 0, 0, 2}, 2},
- };
- arb::sample_tree sm(samples, parents);
- arb::em_morphology em(arb::morphology(sm, false));
-
- auto h1 = arb::reg::cable({0, 0, 0.5});
- auto h2 = arb::reg::cable({0, 0.5, 1});
- auto t1 = arb::reg::tagged(1);
- auto t2 = arb::reg::tagged(2);
- auto all = arb::reg::all();
-
- // Concrete cable lists
- cl h1_{{0, 0.0, 0.5}};
- cl h2_{{0, 0.5, 1.0}};
- 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_EQ(thingify(h1, em), h1_);
- EXPECT_EQ(thingify(h2, em), h2_);
- EXPECT_EQ(thingify(t1, em), t1_);
- EXPECT_EQ(thingify(t2, em), t2_);
- EXPECT_EQ(thingify(join(h1, h2), em), all_);
-
- EXPECT_EQ(thingify(intersect(h1, h2), em), (cl{{0, 0.5, 0.5}}));
-
- EXPECT_EQ(thingify(intersect(h1, h1), em), h1_);
- EXPECT_EQ(thingify(intersect(t1, t1), em), t1_);
- EXPECT_EQ(thingify(join(t1, t2), em), all_);
- EXPECT_EQ(thingify(intersect(all, t1), em), t1_);
- EXPECT_EQ(thingify(intersect(all, t2), em), t2_);
- EXPECT_EQ(thingify(join(all, t1), em), all_);
- EXPECT_EQ(thingify(join(all, t2), em), all_);
- EXPECT_EQ(thingify(join(h1, t1), em), (cl{{0, 0, 0.7}}));
- EXPECT_EQ(thingify(join(h1, t2), em), (cl{{0, 0, 0.5}, {0, 0.7, 1}}));
- EXPECT_EQ(thingify(intersect(h2, t1), em), (cl{{0, 0.5, 0.7}}));
- }
-
-
- // Test handling of spherical soma on multi-branch morphologies
- //
- // sample ids:
- // 0 |
- // 1 3 |
- // 2 4 |
- // tags:
- // 1 |
- // 3 2 |
- // 3 2 |
- {
- pvec parents = {npos, 0, 1, 0, 3};
- svec samples = {
- {{ 0, 0, 0, 2}, 1},
- {{ 10, 0, 0, 2}, 3},
- {{100, 0, 0, 2}, 3},
- {{ 0, 10, 0, 2}, 2},
- {{ 0,100, 0, 2}, 2},
- };
-
- // with a spherical root
- arb::sample_tree sm(samples, parents);
- arb::em_morphology em(arb::morphology(sm, true));
-
- using arb::reg::tagged;
- using arb::reg::branch;
- using arb::reg::all;
-
- EXPECT_EQ(thingify(tagged(1), em), (arb::mcable_list{{0,0,1}}));
- EXPECT_EQ(thingify(tagged(2), em), (arb::mcable_list{{2,0,1}}));
- EXPECT_EQ(thingify(tagged(3), em), (arb::mcable_list{{1,0,1}}));
- EXPECT_EQ(thingify(join(tagged(1), tagged(2), tagged(3)), em), (arb::mcable_list{{0,0,1}, {1,0,1}, {2,0,1}}));
- EXPECT_EQ(thingify(join(tagged(1), tagged(2), tagged(3)), em), thingify(all(), em));
- }
-
- // Test multi-level morphologies.
- //
- // Eight samples
- //
- // sample ids:
- // 0
- // 1 3
- // 2 4
- // 5 6
- // 7
- // tags:
- // 1
- // 3 2
- // 3 2
- // 4 3
- // 3
- {
- pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
- svec samples = {
- {{ 0, 0, 0, 2}, 1},
- {{ 10, 0, 0, 2}, 3},
- {{100, 0, 0, 2}, 3},
- {{ 0, 10, 0, 2}, 2},
- {{ 0,100, 0, 2}, 2},
- {{100,100, 0, 2}, 4},
- {{ 0,200, 0, 2}, 3},
- {{ 0,300, 0, 2}, 3},
- };
- arb::sample_tree sm(samples, parents);
-
- // Without spherical root
- arb::em_morphology em(arb::morphology(sm, false));
-
- using arb::reg::tagged;
- using arb::reg::branch;
- using arb::reg::all;
- using arb::reg::cable;
- using arb::mcable;
- auto soma = tagged(1);
- auto axon = tagged(2);
- auto dend = tagged(3);
- auto apic = tagged(4);
- auto b1 = branch(1);
- 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};
- mcable b2_{2,0,1};
- mcable b3_{3,0,1};
- cl all_ = {b0_,b1_,b2_,b3_};
-
- mcable end1_{1,1,1};
- mcable root_{0,0,0};
-
- EXPECT_EQ(thingify(all(), em), all_);
- EXPECT_EQ(thingify(soma, em), empty_);
- EXPECT_EQ(thingify(axon, em), (cl{b1_}));
- EXPECT_EQ(thingify(dend, em), (cl{b0_,b3_}));
- EXPECT_EQ(thingify(apic, em), (cl{b2_}));
- EXPECT_EQ(thingify(join(dend, apic), em), (cl{b0_,b2_,b3_}));
- EXPECT_EQ(thingify(join(axon, join(dend, apic)), em), all_);
-
- // Test that intersection correctly generates zero-length cables at
- // parent-child interfaces.
- EXPECT_EQ(thingify(intersect(apic, dend), em), (cl{end1_}));
- EXPECT_EQ(thingify(intersect(apic, axon), em), (cl{end1_}));
- EXPECT_EQ(thingify(intersect(axon, dend), em), (cl{root_, end1_}));
-
- // Test some more interesting intersections and unions.
-
- // 123456789 123456789
- // |---------|---------| lhs
- // | ----- | --- | rhs
- // | xxxxx | xxx | rand
- // |xxxxxxxxx|xxxxxxxxx| ror
- auto lhs = b13;
- auto rhs = join(cable({1,.2,.7}), cable({3,.3,.6}));
- auto rand = cl{ {1,.2,.7}, {3,.3,.6}};
- auto ror = cl{ {1,.0,1.}, {3,.0,1.}};
- EXPECT_EQ(thingify(intersect(lhs, rhs), em), rand);
- EXPECT_EQ(thingify(join(lhs, rhs), em), ror);
-
- // Assert communtativity
- std::swap(lhs, rhs);
- EXPECT_EQ(thingify(intersect(lhs, rhs), em), rand);
- EXPECT_EQ(thingify(join(lhs, rhs), em), ror);
-
- // 123456789 123456789
- // | ----- | ---- | lhs
- // | ----- | --- | rhs
- // | xxxx | xx | rand
- // | xxxxxx | xxxxx | ror
- lhs = join(cable({1,.3,.8}), cable({3,.1,.5}));
- rhs = join(cable({1,.2,.7}), cable({3,.3,.6}));
- rand = cl{ {1,.3,.7}, {3,.3,.5}};
- ror = cl{ {1,.2,.8}, {3,.1,.6}};
- EXPECT_EQ(thingify(intersect(lhs, rhs), em), rand);
- EXPECT_EQ(thingify(join(lhs, rhs), em), ror);
-
- // Assert communtativity
- std::swap(lhs, rhs);
- EXPECT_EQ(thingify(intersect(lhs, rhs), em), rand);
- EXPECT_EQ(thingify(join(lhs, rhs), em), ror);
-
- // 123456789 123456789
- // | -- - | --- --- | lhs
- // | ----- | --- | rhs
- // | x x | x x | rand
- // | xxxxxx | xxxxxxx | ror
- lhs = join(cable({1,.1,.3}), cable({1,.4,.5}), cable({3,.1,.4}), cable({3,.5,.9}));
- rhs = join(cable({1,.2,.7}), cable({3,.3,.6}));
- rand = cl{ {1,.2,.3}, {1,.4,.5}, {3,.3,.4}, {3,.5,.6}};
- ror = cl{ {1,.1,.7}, {3,.1,.9}};
- EXPECT_EQ(thingify(intersect(lhs, rhs), em), rand);
- EXPECT_EQ(thingify(join(lhs, rhs), em), ror);
-
- // Assert communtativity
- std::swap(lhs, rhs);
- EXPECT_EQ(thingify(intersect(lhs, rhs), em), rand);
- EXPECT_EQ(thingify(join(lhs, rhs), em), ror);
- }
-}
-
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index 2808d3ea8fb379791b938989fb8e411d9a67bd8f..e1a0a5f82ea84b4bba3156d15955b8765d53bc62 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -11,7 +11,6 @@
#include "util/rangeutil.hpp"
#include "util/span.hpp"
#include "io/sepval.hpp"
-#include "morph/em_morphology.hpp"
#include "common.hpp"
#include "unit_test_catalogue.hpp"
diff --git a/test/unit/test_morph_embedding.cpp b/test/unit/test_morph_embedding.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c72b62fec69a9819ae8e83fa7082c3f740ebd61
--- /dev/null
+++ b/test/unit/test_morph_embedding.cpp
@@ -0,0 +1,241 @@
+#include <cmath>
+#include <vector>
+
+#include <arbor/math.hpp>
+#include <arbor/morph/embed_pwlin.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/sample_tree.hpp>
+
+#include "util/piecewise.hpp"
+
+#include "../test/gtest.h"
+#include "common.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";
+ }
+}
+
+TEST(embedding, samples_and_branch_length) {
+ using pvec = std::vector<msize_t>;
+ using svec = std::vector<msample>;
+ using loc = mlocation;
+
+ // A single unbranched cable with 5 sample points.
+ // The cable has length 10 μm, with samples located at
+ // 0 μm, 1 μm, 3 μm, 7 μm and 10 μm.
+ {
+ pvec parents = {mnpos, 0, 1, 2, 3};
+ svec samples = {
+ {{ 0, 0, 0, 2}, 1},
+ {{ 1, 0, 0, 2}, 1},
+ {{ 3, 0, 0, 2}, 1},
+ {{ 7, 0, 0, 2}, 1},
+ {{ 10, 0, 0, 2}, 1},
+ };
+ sample_tree sm(samples, parents);
+ morphology m(sm, false);
+
+ embedding em(m);
+
+ EXPECT_TRUE(location_eq(m, em.sample_location(0), (loc{0,0})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(1), (loc{0,0.1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(2), (loc{0,0.3})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(3), (loc{0,0.7})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(4), (loc{0,1})));
+
+ EXPECT_EQ(10., em.branch_length(0));
+ }
+
+ // Eight samples
+ //
+ // sample ids:
+ // 0
+ // 1 3
+ // 2 4
+ // 5 6
+ // 7
+ { // Spherical root.
+ pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
+
+ svec samples = {
+ {{ 0, 0, 0, 10}, 1},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 3},
+ {{ 0,100, 0, 2}, 3},
+ {{100,100, 0, 2}, 3},
+ {{ 0,200, 0, 2}, 3},
+ {{ 0,300, 0, 2}, 3},
+ };
+ sample_tree sm(samples, parents);
+ morphology m(sm, true);
+ ASSERT_EQ(5u, m.num_branches());
+
+ embedding em(m);
+
+ EXPECT_TRUE(location_eq(m, em.sample_location(0), (loc{0,0.5})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(1), (loc{1,0})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(2), (loc{1,1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(3), (loc{2,0})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(4), (loc{2,1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(5), (loc{3,1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(6), (loc{4,0.5})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(7), (loc{4,1})));
+
+ EXPECT_EQ(20., em.branch_length(0));
+ EXPECT_EQ(90., em.branch_length(1));
+ EXPECT_EQ(90., em.branch_length(2));
+ EXPECT_EQ(100., em.branch_length(3));
+ EXPECT_EQ(200., em.branch_length(4));
+ }
+ { // No Spherical root
+ pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
+
+ svec samples = {
+ {{ 0, 0, 0, 2}, 1},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 3},
+ {{ 0,100, 0, 2}, 3},
+ {{100,100, 0, 2}, 3},
+ {{ 0,130, 0, 2}, 3},
+ {{ 0,300, 0, 2}, 3},
+ };
+ sample_tree sm(samples, parents);
+ morphology m(sm, false);
+ ASSERT_EQ(4u, m.num_branches());
+
+ embedding em(m);
+
+ EXPECT_TRUE(location_eq(m, em.sample_location(0), (loc{0,0})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(1), (loc{0,0.1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(2), (loc{0,1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(3), (loc{1,0.1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(4), (loc{1,1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(5), (loc{2,1})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(6), (loc{3,0.15})));
+ EXPECT_TRUE(location_eq(m, em.sample_location(7), (loc{3,1})));
+
+ EXPECT_EQ(100., em.branch_length(0));
+ EXPECT_EQ(100., em.branch_length(1));
+ EXPECT_EQ(100., em.branch_length(2));
+ EXPECT_EQ(200., em.branch_length(3));
+ }
+}
+
+// TODO: integrator tests
+
+
+TEST(embedding, partial_branch_length) {
+ using pvec = std::vector<msize_t>;
+ using svec = std::vector<msample>;
+ using util::pw_elements;
+
+ pvec parents = {mnpos, 0, 1, 2, 2};
+ svec samples = {
+ {{ 0, 0, 0, 10}, 1},
+ {{ 10, 0, 0, 20}, 1},
+ {{ 30, 0, 0, 10}, 1},
+ {{ 30, 10, 0, 5}, 2},
+ {{ 30, 0, 50, 5}, 2}
+ };
+
+ morphology m(sample_tree(samples, parents), false);
+ embedding em(m);
+
+ EXPECT_DOUBLE_EQ(30., em.branch_length(0));
+ EXPECT_DOUBLE_EQ(30., em.integrate_length(mcable{0, 0., 1.}));
+ EXPECT_DOUBLE_EQ(15., em.integrate_length(mcable{0, 0.25, 0.75}));
+
+ EXPECT_DOUBLE_EQ(10., em.branch_length(1));
+ EXPECT_DOUBLE_EQ(10., em.integrate_length(mcable{1, 0., 1.}));
+ EXPECT_DOUBLE_EQ(7.5, em.integrate_length(mcable{1, 0.25, 1.0}));
+
+ // expect 2*0.25+3*0.5 = 2.0 times corresponding cable length.
+ pw_elements<double> pw({0.25, 0.5, 1.}, {2., 3.});
+ EXPECT_DOUBLE_EQ(20., em.integrate_length(1, pw));
+}
+
+TEST(embedding, partial_area) {
+ using pvec = std::vector<msize_t>;
+ using svec = std::vector<msample>;
+ using util::pw_elements;
+ using testing::near_relative;
+
+ pvec parents = {mnpos, 0, 1, 2, 2};
+ svec samples = {
+ {{ 0, 0, 0, 10}, 1},
+ {{ 10, 0, 0, 20}, 1},
+ {{ 30, 0, 0, 10}, 1},
+ {{ 30, 10, 0, 5}, 2},
+ {{ 30, 0, 50, 5}, 2}
+ };
+
+ morphology m(sample_tree(samples, parents), false);
+ embedding em(m);
+
+ // Cable 1: single truncated cone, length L = 10,
+ // radius râ‚€ = 10 (pos 0) to râ‚ = 5 (pos 1).
+ //
+ // Expect cable area = 2πLr√(1 + m²)
+ // where m = δr/L and r = (râ‚€+râ‚)/2 = râ‚€ + δr/2.
+
+ constexpr double pi = math::pi<double>;
+ double cable1_area = 2*pi*10*7.5*std::sqrt(1.25);
+
+ constexpr double reltol = 1e-10;
+
+ EXPECT_TRUE(near_relative(cable1_area, em.integrate_area(mcable{1, 0., 1.}), reltol));
+
+ // Weighted area within cable 0:
+ // a) proportional segment [0.1, 0.3]:
+ // truncated cone length 6,
+ // râ‚€ = 13; râ‚ = 19, slope = 1
+ //
+ // b) proportional segment [0.3, 0.9]:
+ // truncated cone length 1,
+ // râ‚€ = 19, râ‚ = 20, slope = 1
+ // truncated cone length 17
+ // râ‚€ = 20, râ‚ = 11.5, slope = -0.5
+
+ EXPECT_TRUE(near_relative(13., em.radius(mlocation{0, 0.1}), reltol));
+ EXPECT_TRUE(near_relative(19., em.radius(mlocation{0, 0.3}), reltol));
+ EXPECT_TRUE(near_relative(11.5, em.radius(mlocation{0, 0.9}), reltol));
+
+ pw_elements<double> pw({0.1, 0.3, 0.9}, {5., 7.});
+ double sub_area1 = pi*6*(13+19)*std::sqrt(2);
+ double sub_area2 = pi*1*(19+20)*std::sqrt(2);
+ double sub_area3 = pi*17*(20+11.5)*std::sqrt(1.25);
+
+ EXPECT_TRUE(near_relative(sub_area1, em.integrate_area(mcable{0, 0.1, 0.3}), reltol));
+ EXPECT_TRUE(near_relative(sub_area2, em.integrate_area(mcable{0, 0.3, 1/3.}), reltol));
+ EXPECT_TRUE(near_relative(sub_area3, em.integrate_area(mcable{0, 1/3., 0.9}), reltol));
+
+ double expected_pw_area = 5.*sub_area1+7.*(sub_area2+sub_area3);
+ EXPECT_TRUE(near_relative(expected_pw_area, em.integrate_area(0, pw), reltol));
+
+ // Integrated inverse cross-sectional area in cable 1 from 0.1 to 0.4:
+ // radius râ‚€ = 9.5, râ‚ = 8, length = 3.
+
+ double expected_ixa = 3/(9.5*8)/pi;
+ EXPECT_TRUE(near_relative(expected_ixa, em.integrate_ixa(mcable{1, 0.1, 0.4}), reltol));
+}
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6304baa47e765974f9f727cc2765865f85626c7f
--- /dev/null
+++ b/test/unit/test_morph_expr.cpp
@@ -0,0 +1,475 @@
+#include "../test/gtest.h"
+
+#include <vector>
+
+#include <arbor/morph/embed_pwlin.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/morphexcept.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/mprovider.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/morph/sample_tree.hpp>
+
+#include "util/span.hpp"
+#include "util/strprintf.hpp"
+
+using namespace arb;
+using embedding = embed_pwlin;
+
+::testing::AssertionResult cable_eq(mcable a, mcable b) {
+ if (a.branch!=b.branch) {
+ return ::testing::AssertionFailure()
+ << "cables " << a << " and " << b << " differ";
+ }
+
+ using FP = testing::internal::FloatingPoint<double>;
+ if (FP(a.prox_pos).AlmostEquals(FP(b.prox_pos)) && FP(a.dist_pos).AlmostEquals(FP(b.dist_pos))) {
+ return ::testing::AssertionSuccess();
+ }
+ else {
+ return ::testing::AssertionFailure()
+ << "cables " << a << " and " << b << " differ";
+ }
+}
+
+::testing::AssertionResult cablelist_eq(const mcable_list& as, const mcable_list& bs) {
+ if (as.size()!=bs.size()) {
+ return ::testing::AssertionFailure()
+ << "cablelists " << as << " and " << bs << " differ";
+ }
+
+ for (auto i: util::count_along(as)) {
+ auto result = cable_eq(as[i], bs[i]);
+ if (!result) return ::testing::AssertionFailure()
+ << "cablelists " << as << " and " << bs << " differ";
+ }
+ return ::testing::AssertionSuccess();
+}
+
+TEST(region, expr_repn) {
+ using util::to_string;
+
+ auto c1 = reg::cable({1, 0, 1});
+ auto c2 = reg::cable({4, 0.125, 0.5});
+ auto c3 = join(reg::cable({4, 0.125, 0.5}), reg::cable({3, 0, 1}));
+ auto b1 = reg::branch(1);
+ auto t1 = reg::tagged(1);
+ auto t2 = reg::tagged(2);
+ auto t3 = reg::tagged(3);
+ auto all = reg::all();
+
+ EXPECT_EQ(to_string(c1), "(cable 1 0 1)");
+ EXPECT_EQ(to_string(c2), "(cable 4 0.125 0.5)");
+ EXPECT_EQ(to_string(c3), "(join (cable 4 0.125 0.5) (cable 3 0 1))");
+ EXPECT_EQ(to_string(b1), "(cable 1 0 1)");
+ EXPECT_EQ(to_string(t1), "(tag 1)");
+ EXPECT_EQ(to_string(t2), "(tag 2)");
+ EXPECT_EQ(to_string(intersect(c1, t2)), "(intersect (cable 1 0 1) (tag 2))");
+ EXPECT_EQ(to_string(join(c1, t2)), "(join (cable 1 0 1) (tag 2))");
+ EXPECT_EQ(to_string(join(t1, t2, t3)), "(join (join (tag 1) (tag 2)) (tag 3))");
+ EXPECT_EQ(to_string(intersect(t1, t2, t3)), "(intersect (intersect (tag 1) (tag 2)) (tag 3))");
+ EXPECT_EQ(to_string(intersect(join(c1, t2), c2)), "(intersect (join (cable 1 0 1) (tag 2)) (cable 4 0.125 0.5))");
+ EXPECT_EQ(to_string(all), "all");
+}
+
+TEST(region, invalid_mcable) {
+ EXPECT_NO_THROW(reg::cable({123, 0.5, 0.8}));
+ EXPECT_THROW(reg::cable({1, 0, 1.1}), invalid_mcable);
+ EXPECT_THROW(reg::branch(-1), invalid_mcable);
+}
+
+TEST(locset, expr_repn) {
+ using util::to_string;
+
+ auto root = ls::root();
+ auto term = ls::terminal();
+ auto samp = ls::sample(42);
+ auto loc = ls::location({2, 0.5});
+
+ EXPECT_EQ(to_string(root), "root");
+ EXPECT_EQ(to_string(term), "terminal");
+ EXPECT_EQ(to_string(sum(root, term)), "(sum root terminal)");
+ EXPECT_EQ(to_string(sum(root, term, samp)), "(sum (sum root terminal) (sample 42))");
+ EXPECT_EQ(to_string(sum(root, term, samp, loc)), "(sum (sum (sum root terminal) (sample 42)) (location 2 0.5))");
+ EXPECT_EQ(to_string(samp), "(sample 42)");
+ EXPECT_EQ(to_string(loc), "(location 2 0.5)");
+}
+
+TEST(region, invalid_mlocation) {
+ // Location positions have to be in the range [0,1].
+ EXPECT_NO_THROW(ls::location({123, 0.0}));
+ EXPECT_NO_THROW(ls::location({123, 0.02}));
+ EXPECT_NO_THROW(ls::location({123, 1.0}));
+
+ EXPECT_THROW(ls::location({0, 1.5}), invalid_mlocation);
+ EXPECT_THROW(ls::location({unsigned(-1), 0.}), invalid_mlocation);
+}
+
+// Name evaluation (thingify) tests:
+
+TEST(locset, thingify_named) {
+ using pvec = std::vector<msize_t>;
+ using svec = std::vector<msample>;
+
+ locset banana = ls::root();
+ locset cake = ls::terminal();
+
+ sample_tree sm(svec{ {{0,0,0,1},1}, {{10,0,0,1},1} }, pvec{mnpos, 0});
+ {
+ label_dict dict;
+ dict.set("banana", banana);
+ dict.set("cake", cake);
+
+ mprovider mp(morphology(sm, false), dict);
+ EXPECT_EQ(thingify(locset("cake"), mp), thingify(cake, mp));
+ EXPECT_EQ(thingify(locset("banana"), mp), thingify(banana, mp));
+
+ EXPECT_THROW(thingify(locset("durian"), mp), unbound_name);
+ }
+ {
+ label_dict dict;
+ dict.set("banana", banana);
+ dict.set("cake", cake);
+ dict.set("topping", locset("fruit"));
+ dict.set("fruit", locset("strawberry"));
+
+ EXPECT_THROW(mprovider(morphology(sm, false), dict), unbound_name);
+ }
+ {
+ label_dict dict;
+ dict.set("banana", banana);
+ dict.set("cake", cake);
+ dict.set("topping", locset("fruit"));
+ dict.set("fruit", sum(locset("banana"), locset("topping")));
+
+ EXPECT_THROW(mprovider(morphology(sm, false), dict), circular_definition);
+ }
+}
+
+TEST(region, thingify_named) {
+ using pvec = std::vector<msize_t>;
+ using svec = std::vector<msample>;
+
+ region banana = reg::branch(0);
+ region cake = reg::cable(mcable{0, 0.2, 0.3});
+
+ // copy-paste ftw
+
+ sample_tree sm(svec{ {{0,0,0,1},1}, {{10,0,0,1},1} }, pvec{mnpos, 0});
+ {
+ label_dict dict;
+ dict.set("banana", banana);
+ dict.set("cake", cake);
+
+ mprovider mp(morphology(sm, false), dict);
+ EXPECT_EQ(thingify(region("cake"), mp), thingify(cake, mp));
+ EXPECT_EQ(thingify(region("banana"), mp), thingify(banana, mp));
+
+ EXPECT_THROW(thingify(region("durian"), mp), unbound_name);
+ }
+ {
+ label_dict dict;
+ dict.set("banana", banana);
+ dict.set("cake", cake);
+ dict.set("topping", region("fruit"));
+ dict.set("fruit", region("strawberry"));
+
+ EXPECT_THROW(mprovider(morphology(sm, false), dict), unbound_name);
+ }
+ {
+ label_dict dict;
+ dict.set("banana", banana);
+ dict.set("cake", cake);
+ dict.set("topping", region("fruit"));
+ dict.set("fruit", join(region("cake"), region("topping")));
+
+ EXPECT_THROW(mprovider(morphology(sm, false), dict), circular_definition);
+ }
+}
+
+// Embedded evaluation (thingify) tests:
+
+TEST(locset, thingify) {
+ using pvec = std::vector<msize_t>;
+ using svec = std::vector<msample>;
+ using ll = mlocation_list;
+ auto root = ls::root();
+ auto term = ls::terminal();
+ auto samp = ls::sample(4);
+ auto midb2 = ls::location({2, 0.5});
+ auto midb1 = ls::location({1, 0.5});
+ auto begb0 = ls::location({0, 0});
+ auto begb1 = ls::location({1, 0});
+ auto begb2 = ls::location({2, 0});
+ auto begb3 = ls::location({3, 0});
+ auto begb4 = ls::location({4, 0});
+
+ // Eight samples
+ //
+ // 0
+ // 1 3
+ // 2 4
+ // 5 6
+ // 7
+ pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
+ svec samples = {
+ {{ 0, 0, 0, 2}, 3},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 3},
+ {{ 0,100, 0, 2}, 3},
+ {{100,100, 0, 2}, 3},
+ {{ 0,200, 0, 2}, 3},
+ {{ 0,300, 0, 2}, 3},
+ };
+ sample_tree sm(samples, parents);
+
+ {
+ mprovider mp(morphology(sm, true));
+
+ EXPECT_EQ(thingify(root, mp), (ll{{0,0}}));
+ EXPECT_EQ(thingify(term, mp), (ll{{1,1},{3,1},{4,1}}));
+ EXPECT_EQ(thingify(samp, mp), (ll{{2,1}}));
+ EXPECT_EQ(thingify(midb2, mp), (ll{{2,0.5}}));
+ EXPECT_EQ(thingify(midb1, mp), (ll{{1,0.5}}));
+ EXPECT_EQ(thingify(begb0, mp), (ll{{0,0}}));
+ EXPECT_EQ(thingify(begb1, mp), (ll{{1,0}}));
+ EXPECT_EQ(thingify(begb2, mp), (ll{{2,0}}));
+ EXPECT_EQ(thingify(begb3, mp), (ll{{3,0}}));
+ EXPECT_EQ(thingify(begb4, mp), (ll{{4,0}}));
+ }
+ {
+ mprovider mp(morphology(sm, false));
+
+ EXPECT_EQ(thingify(root, mp), (ll{{0,0}}));
+ EXPECT_EQ(thingify(term, mp), (ll{{0,1},{2,1},{3,1}}));
+ EXPECT_EQ(thingify(samp, mp), (ll{{1,1}}));
+ EXPECT_EQ(thingify(midb2, mp), (ll{{2,0.5}}));
+ EXPECT_EQ(thingify(midb1, mp), (ll{{1,0.5}}));
+ EXPECT_EQ(thingify(begb0, mp), (ll{{0,0}}));
+ EXPECT_EQ(thingify(begb1, mp), (ll{{1,0}}));
+ EXPECT_EQ(thingify(begb2, mp), (ll{{2,0}}));
+ EXPECT_EQ(thingify(begb3, mp), (ll{{3,0}}));
+
+ // In the absence of a spherical root, there is no branch 4.
+ EXPECT_THROW(thingify(begb4, mp), no_such_branch);
+ }
+}
+
+TEST(region, thingify) {
+ using pvec = std::vector<msize_t>;
+ using svec = std::vector<msample>;
+ using cl = mcable_list;
+
+ // A single unbranched cable with 5 sample points.
+ // The cable has length 10 μm, with samples located at
+ // 0 μm, 1 μm, 3 μm, 7 μm and 10 μm.
+ {
+ pvec parents = {mnpos, 0, 1, 2, 3};
+ svec samples = {
+ {{ 0, 0, 0, 2}, 1},
+ {{ 1, 0, 0, 2}, 1},
+ {{ 3, 0, 0, 2}, 2},
+ {{ 7, 0, 0, 2}, 1},
+ {{ 10, 0, 0, 2}, 2},
+ };
+ sample_tree sm(samples, parents);
+ mprovider mp(morphology(sm, false));
+
+ auto h1 = reg::cable({0, 0, 0.5});
+ auto h2 = reg::cable({0, 0.5, 1});
+ auto t1 = reg::tagged(1);
+ auto t2 = reg::tagged(2);
+ auto all = reg::all();
+
+ // Concrete cable lists
+ cl h1_{{0, 0.0, 0.5}};
+ cl h2_{{0, 0.5, 1.0}};
+ 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_EQ(thingify(h1, mp), h1_);
+ EXPECT_EQ(thingify(h2, mp), h2_);
+ EXPECT_EQ(thingify(join(h1, h2), mp), all_);
+ EXPECT_EQ(thingify(intersect(h1, h2), mp), (cl{{0, 0.5, 0.5}}));
+
+ EXPECT_TRUE(cablelist_eq(thingify(t1, mp), t1_));
+ EXPECT_TRUE(cablelist_eq(thingify(t2, mp), t2_));
+ EXPECT_TRUE(cablelist_eq(thingify(intersect(h1, h1), mp), h1_));
+ EXPECT_TRUE(cablelist_eq(thingify(intersect(t1, t1), mp), t1_));
+ EXPECT_TRUE(cablelist_eq(thingify(join(t1, t2), mp), all_));
+ EXPECT_TRUE(cablelist_eq(thingify(intersect(all, t1), mp), t1_));
+ EXPECT_TRUE(cablelist_eq(thingify(intersect(all, t2), mp), t2_));
+ EXPECT_TRUE(cablelist_eq(thingify(join(all, t1), mp), all_));
+ EXPECT_TRUE(cablelist_eq(thingify(join(all, t2), mp), all_));
+ EXPECT_TRUE(cablelist_eq(thingify(join(h1, t1), mp), (cl{{0, 0, 0.7}})));
+ EXPECT_TRUE(cablelist_eq(thingify(join(h1, t2), mp), (cl{{0, 0, 0.5}, {0, 0.7, 1}})));
+ EXPECT_TRUE(cablelist_eq(thingify(intersect(h2, t1), mp), (cl{{0, 0.5, 0.7}})));
+ }
+
+
+ // Test handling of spherical soma on multi-branch morphologies
+ //
+ // sample ids:
+ // 0 |
+ // 1 3 |
+ // 2 4 |
+ // tags:
+ // 1 |
+ // 3 2 |
+ // 3 2 |
+ {
+ pvec parents = {mnpos, 0, 1, 0, 3};
+ svec samples = {
+ {{ 0, 0, 0, 2}, 1},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 2},
+ {{ 0,100, 0, 2}, 2},
+ };
+
+ // with a spherical root
+ sample_tree sm(samples, parents);
+ mprovider mp(morphology(sm, true));
+
+ using reg::tagged;
+ using reg::branch;
+ using reg::all;
+
+ EXPECT_EQ(thingify(tagged(1), mp), (mcable_list{{0,0,1}}));
+ EXPECT_EQ(thingify(tagged(2), mp), (mcable_list{{2,0,1}}));
+ EXPECT_EQ(thingify(tagged(3), mp), (mcable_list{{1,0,1}}));
+ EXPECT_EQ(thingify(join(tagged(1), tagged(2), tagged(3)), mp), (mcable_list{{0,0,1}, {1,0,1}, {2,0,1}}));
+ EXPECT_EQ(thingify(join(tagged(1), tagged(2), tagged(3)), mp), thingify(all(), mp));
+ }
+
+ // Test multi-level morphologies.
+ //
+ // Eight samples
+ //
+ // sample ids:
+ // 0
+ // 1 3
+ // 2 4
+ // 5 6
+ // 7
+ // tags:
+ // 1
+ // 3 2
+ // 3 2
+ // 4 3
+ // 3
+ {
+ pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
+ svec samples = {
+ {{ 0, 0, 0, 2}, 1},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 2},
+ {{ 0,100, 0, 2}, 2},
+ {{100,100, 0, 2}, 4},
+ {{ 0,200, 0, 2}, 3},
+ {{ 0,300, 0, 2}, 3},
+ };
+ sample_tree sm(samples, parents);
+
+ // Without spherical root
+ mprovider mp(morphology(sm, false));
+
+ using reg::tagged;
+ using reg::branch;
+ using reg::all;
+ using reg::cable;
+
+ auto soma = tagged(1);
+ auto axon = tagged(2);
+ auto dend = tagged(3);
+ auto apic = tagged(4);
+ auto b1 = branch(1);
+ 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};
+ mcable b2_{2,0,1};
+ mcable b3_{3,0,1};
+ cl all_ = {b0_,b1_,b2_,b3_};
+
+ mcable end1_{1,1,1};
+ mcable root_{0,0,0};
+
+ EXPECT_EQ(thingify(all(), mp), all_);
+ EXPECT_EQ(thingify(soma, mp), empty_);
+ EXPECT_EQ(thingify(axon, mp), (cl{b1_}));
+ EXPECT_EQ(thingify(dend, mp), (cl{b0_,b3_}));
+ EXPECT_EQ(thingify(apic, mp), (cl{b2_}));
+ EXPECT_EQ(thingify(join(dend, apic), mp), (cl{b0_,b2_,b3_}));
+ EXPECT_EQ(thingify(join(axon, join(dend, apic)), mp), all_);
+
+ // Test that intersection correctly generates zero-length cables at
+ // parent-child interfaces.
+ EXPECT_EQ(thingify(intersect(apic, dend), mp), (cl{end1_}));
+ EXPECT_EQ(thingify(intersect(apic, axon), mp), (cl{end1_}));
+ EXPECT_EQ(thingify(intersect(axon, dend), mp), (cl{root_, end1_}));
+
+ // Test some more interesting intersections and unions.
+
+ // 123456789 123456789
+ // |---------|---------| lhs
+ // | ----- | --- | rhs
+ // | xxxxx | xxx | rand
+ // |xxxxxxxxx|xxxxxxxxx| ror
+ auto lhs = b13;
+ auto rhs = join(cable({1,.2,.7}), cable({3,.3,.6}));
+ auto rand = cl{ {1,.2,.7}, {3,.3,.6}};
+ auto ror = cl{ {1,.0,1.}, {3,.0,1.}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // Assert communtativity
+ std::swap(lhs, rhs);
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // 123456789 123456789
+ // | ----- | ---- | lhs
+ // | ----- | --- | rhs
+ // | xxxx | xx | rand
+ // | xxxxxx | xxxxx | ror
+ lhs = join(cable({1,.3,.8}), cable({3,.1,.5}));
+ rhs = join(cable({1,.2,.7}), cable({3,.3,.6}));
+ rand = cl{ {1,.3,.7}, {3,.3,.5}};
+ ror = cl{ {1,.2,.8}, {3,.1,.6}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // Assert communtativity
+ std::swap(lhs, rhs);
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // 123456789 123456789
+ // | -- - | --- --- | lhs
+ // | ----- | --- | rhs
+ // | x x | x x | rand
+ // | xxxxxx | xxxxxxx | ror
+ lhs = join(cable({1,.1,.3}), cable({1,.4,.5}), cable({3,.1,.4}), cable({3,.5,.9}));
+ rhs = join(cable({1,.2,.7}), cable({3,.3,.6}));
+ rand = cl{ {1,.2,.3}, {1,.4,.5}, {3,.3,.4}, {3,.5,.6}};
+ ror = cl{ {1,.1,.7}, {3,.1,.9}};
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+
+ // Assert communtativity
+ std::swap(lhs, rhs);
+ EXPECT_EQ(thingify(intersect(lhs, rhs), mp), rand);
+ EXPECT_EQ(thingify(join(lhs, rhs), mp), ror);
+ }
+}
diff --git a/test/unit/test_morph_primitives.cpp b/test/unit/test_morph_primitives.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c627e589b9fe31bb10c7f25f2b6eea7012ef3a6b
--- /dev/null
+++ b/test/unit/test_morph_primitives.cpp
@@ -0,0 +1,122 @@
+#include <iostream>
+#include <vector>
+
+#include <arbor/morph/primitives.hpp>
+
+#include "../gtest.h"
+
+using namespace arb;
+
+::testing::AssertionResult mpoint_eq(const mpoint& a, const mpoint& b) {
+ if (a.x==b.x && a.y==b.y && a.z==b.z && a.radius==b.radius) {
+ return ::testing::AssertionSuccess();
+ }
+ else {
+ return ::testing::AssertionFailure()
+ << "mpoints " << a << " and " << b << " differ.";
+ }
+}
+
+TEST(morph_primitives, lerp) {
+ mpoint a{1., 2., 3., 4.};
+ mpoint b{2., 4., 7., 12.};
+
+ EXPECT_TRUE(mpoint_eq(a, lerp(a, b, 0)));
+ EXPECT_TRUE(mpoint_eq(b, lerp(a, b, 1)));
+
+ EXPECT_TRUE(mpoint_eq(mpoint{1.5, 3., 5., 8.}, lerp(a, b, 0.5)));
+ EXPECT_TRUE(mpoint_eq(mpoint{1.25, 2.5, 4., 6.}, lerp(a, b, 0.25)));
+}
+
+TEST(morph_primitives, is_collocated) {
+ EXPECT_TRUE(is_collocated(mpoint{1., 2., 3., 4.}, mpoint{1., 2., 3., 4.}));
+ EXPECT_TRUE(is_collocated(mpoint{1., 2., 3., 4.}, mpoint{1., 2., 3., 4.5}));
+ EXPECT_FALSE(is_collocated(mpoint{1., 2., 2.5, 4.}, mpoint{1., 2., 3., 4}));
+ EXPECT_FALSE(is_collocated(mpoint{1., 2.5, 3., 4.}, mpoint{1., 2., 3., 4}));
+ EXPECT_FALSE(is_collocated(mpoint{2.5, 2, 3., 4.}, mpoint{1., 2., 3., 4}));
+
+ EXPECT_TRUE(is_collocated(msample{{1., 2., 3., 4.},3}, msample{{1., 2., 3., 4.},3}));
+ EXPECT_TRUE(is_collocated(msample{{1., 2., 3., 4.},3}, msample{{1., 2., 3., 4.5},1}));
+ EXPECT_FALSE(is_collocated(msample{{1., 2., 2.5, 4.},3}, msample{{1., 2., 3., 4},3}));
+ EXPECT_FALSE(is_collocated(msample{{1., 2.5, 3., 4.},3}, msample{{1., 2., 3., 4},3}));
+ EXPECT_FALSE(is_collocated(msample{{2.5, 2, 3., 4.},3}, msample{{1., 2., 3., 4},3}));
+}
+
+TEST(morph_primitives, distance) {
+ // 2² + 3² + 6² = 7²
+ EXPECT_DOUBLE_EQ(7., distance(mpoint{1.5, 2.5, 3.5, 4.5}, mpoint{3.5, 5.5, 9.5, 0.31}));
+}
+
+TEST(morph_primitives, join_intersect_sum) {
+ auto ml = [](const std::vector<int>& bids) {
+ mlocation_list L;
+ for (msize_t b: bids) L.push_back({b, 0});
+ return L;
+ };
+
+ using ll = mlocation_list;
+
+ {
+ ll lhs{};
+ ll rhs{};
+ EXPECT_EQ(sum(lhs, rhs), ll{});
+ EXPECT_EQ(join(lhs, rhs), ll{});
+ EXPECT_EQ(intersection(lhs, rhs), ll{});
+ }
+ {
+ ll lhs{};
+ ll rhs = ml({0,1});
+ EXPECT_EQ(sum(lhs, rhs), rhs);
+ EXPECT_EQ(join(lhs, rhs), rhs);
+ EXPECT_EQ(intersection(lhs, rhs), ll{});
+ }
+ {
+ ll lhs = ml({1});
+ ll rhs = ml({1});
+ EXPECT_EQ(sum(lhs, rhs), ml({1,1}));
+ EXPECT_EQ(join(lhs, rhs), ml({1}));
+ EXPECT_EQ(intersection(lhs, rhs), ml({1}));
+ }
+ {
+ ll lhs = ml({1,1});
+ ll rhs = ml({1});
+ EXPECT_EQ(sum(lhs, rhs), ml({1,1,1}));
+ EXPECT_EQ(join(lhs, rhs), ml({1,1}));
+ EXPECT_EQ(intersection(lhs, rhs), ml({1}));
+ }
+ {
+ ll lhs = ml({0,3});
+ ll rhs = ml({1,2});
+ EXPECT_EQ(sum(lhs, rhs), ml({0,1,2,3}));
+ EXPECT_EQ(join(lhs, rhs), ml({0,1,2,3}));
+ EXPECT_EQ(intersection(lhs, rhs), ll{});
+ }
+ {
+ ll lhs = ml({0,1,3});
+ ll rhs = ml({0,1,3});
+ EXPECT_EQ(sum(lhs, rhs), ml({0,0,1,1,3,3}));
+ EXPECT_EQ(join(lhs, rhs), lhs);
+ EXPECT_EQ(intersection(lhs, rhs), lhs);
+ }
+ {
+ ll lhs = ml({0,1,3});
+ ll rhs = ml({1,2});
+ EXPECT_EQ(sum(lhs, rhs), ml({0,1,1,2,3}));
+ EXPECT_EQ(join(lhs, rhs), ml({0,1,2,3}));
+ EXPECT_EQ(intersection(lhs, rhs), ml({1}));
+ }
+ {
+ ll lhs = ml({0,1,1,3});
+ ll rhs = ml({1,2});
+ EXPECT_EQ(sum(lhs, rhs), ml({0,1,1,1,2,3}));
+ EXPECT_EQ(join(lhs, rhs), ml({0,1,1,2,3}));
+ EXPECT_EQ(intersection(lhs, rhs), ml({1}));
+ }
+ {
+ ll lhs = ml({0,1,1,3,5,5,12});
+ ll rhs = ml({1,2,2,5,5,5});
+ EXPECT_EQ(sum(lhs, rhs), ml({0,1,1,1,2,2,3,5,5,5,5,5,12}));
+ EXPECT_EQ(join(lhs, rhs), ml({0,1,1,2,2,3,5,5,5,12}));
+ EXPECT_EQ(intersection(lhs, rhs), ml({1,5,5}));
+ }
+}
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
index eb34f28939ec7574e02bbd9fa51d3cb4f6cf3c49..d4d51b7e572ada70f2dbf3d73a54d84ed2491c2c 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -1,7 +1,3 @@
-/*
- * Unit tests for sample_tree and morphology.
- */
-
#include <fstream>
#include <cmath>
#include <string>
@@ -9,24 +5,13 @@
#include "../test/gtest.h"
-#include <arbor/math.hpp>
-#include <arbor/morph/error.hpp>
+#include <arbor/morph/morphexcept.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/sample_tree.hpp>
#include <arbor/cable_cell.hpp>
-#include "algorithms.hpp"
-#include "io/sepval.hpp"
#include "morph/mbranch.hpp"
-#include "morph/em_morphology.hpp"
-#include "util/range.hpp"
#include "util/span.hpp"
-#include "util/strprintf.hpp"
-
-template <typename T>
-std::ostream& operator<<(std::ostream& o, const std::vector<T>& v) {
- return o << "[" << arb::io::csv(v) << "]";
-}
// Test basic functions on properties of mpoint
TEST(morphology, mpoint) {
@@ -175,7 +160,7 @@ TEST(morphology, branches_from_parent_index) {
// A cable morphology can't be constructed from a single sample.
EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), false),
- arb::morphology_error);
+ arb::incomplete_branch);
}
{
pvec parents = {npos, 0};
@@ -186,7 +171,7 @@ TEST(morphology, branches_from_parent_index) {
// A morphology can't be constructed with a spherical soma from two samples.
EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true),
- arb::morphology_error);
+ arb::incomplete_branch);
}
{
@@ -211,7 +196,7 @@ TEST(morphology, branches_from_parent_index) {
auto tree = make_tree(parents);
// A spherical root is not valid: each cable branch would have only one sample.
- EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true), arb::morphology_error);
+ EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true), arb::incomplete_branch);
// Two cables, with two samples each, with the first sample in each being the root
auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
@@ -248,7 +233,7 @@ TEST(morphology, branches_from_parent_index) {
EXPECT_EQ(mb({0,3},npos), bc[1]);
// A spherical root is not valid: the second cable branch would have only one sample.
- EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true), arb::morphology_error);
+ EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true), arb::incomplete_branch);
}
{
@@ -362,7 +347,7 @@ TEST(morphology, construction) {
{{0.0, 5.0, 0.0, 1.0}, 2}};
arb::sample_tree sm(s, p);
- EXPECT_THROW((arb::morphology(sm)), arb::morphology_error);
+ EXPECT_THROW((arb::morphology(sm)), arb::incomplete_branch);
}
}
{
@@ -392,6 +377,17 @@ TEST(morphology, branches) {
using pvec = std::vector<arb::msize_t>;
using svec = std::vector<arb::msample>;
auto npos = arb::mnpos;
+
+ auto check_terminal_branches = [](const arb::morphology& m) {
+ pvec expected;
+ arb::msize_t n = m.num_branches();
+
+ for (arb::msize_t i = 0; i<n; ++i) {
+ if (m.branch_children(i).empty()) expected.push_back(i);
+ }
+ EXPECT_EQ(expected, m.terminal_branches());
+ };
+
{
// 0
pvec parents = {npos};
@@ -404,6 +400,8 @@ TEST(morphology, branches) {
EXPECT_EQ(1u, m.num_branches());
EXPECT_EQ(npos, m.branch_parent(0));
EXPECT_EQ(pvec{}, m.branch_children(0));
+
+ check_terminal_branches(m);
}
{
// 0 - 1
@@ -418,6 +416,8 @@ TEST(morphology, branches) {
EXPECT_EQ(1u, m.num_branches());
EXPECT_EQ(npos, m.branch_parent(0));
EXPECT_EQ(pvec{}, m.branch_children(0));
+
+ check_terminal_branches(m);
}
{
// 0 - 1 - 2
@@ -435,6 +435,8 @@ TEST(morphology, branches) {
EXPECT_EQ(1u, m.num_branches());
EXPECT_EQ(npos, m.branch_parent(0));
EXPECT_EQ(pvec{}, m.branch_children(0));
+
+ check_terminal_branches(m);
}
{
// First sample has unique tag -> spherical soma attached to a single-segment cable.
@@ -451,6 +453,8 @@ TEST(morphology, branches) {
EXPECT_EQ(0u, m.branch_parent(1));
EXPECT_EQ(pvec{1}, m.branch_children(0));
EXPECT_EQ(pvec{}, m.branch_children(1));
+
+ check_terminal_branches(m);
}
}
{
@@ -470,6 +474,8 @@ TEST(morphology, branches) {
EXPECT_EQ(npos, m.branch_parent(1));
EXPECT_EQ(pvec{}, m.branch_children(0));
EXPECT_EQ(pvec{}, m.branch_children(1));
+
+ check_terminal_branches(m);
}
{
// 0 - 1 - 2 - 3
@@ -521,6 +527,8 @@ TEST(morphology, branches) {
EXPECT_EQ((pvec{3,4}), m.branch_children(2));
EXPECT_EQ((pvec{}), m.branch_children(3));
EXPECT_EQ((pvec{}), m.branch_children(4));
+
+ check_terminal_branches(m);
}
{
svec samples = {
@@ -545,6 +553,8 @@ TEST(morphology, branches) {
EXPECT_EQ((pvec{2,3}), m.branch_children(1));
EXPECT_EQ((pvec{}), m.branch_children(2));
EXPECT_EQ((pvec{}), m.branch_children(3));
+
+ check_terminal_branches(m);
}
}
}
@@ -574,3 +584,54 @@ TEST(morphology, swc) {
EXPECT_EQ(31u, c.num_branches());
}
+TEST(morphology, minset) {
+ using pvec = std::vector<arb::msize_t>;
+ using svec = std::vector<arb::msample>;
+ using ll = arb::mlocation_list;
+ constexpr auto npos = arb::mnpos;
+
+ // Eight samples
+ // no-sphere sphere
+ // sample branch branch
+ // 0 0 0
+ // 1 3 0 1 1 2
+ // 2 4 0 1 1 2
+ // 5 6 2 3 3 4
+ // 7 3 4
+ pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
+ svec samples = {
+ {{ 0, 0, 0, 2}, 3},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 3},
+ {{ 0,100, 0, 2}, 3},
+ {{100,100, 0, 2}, 3},
+ {{ 0,200, 0, 2}, 3},
+ {{ 0,300, 0, 2}, 3},
+ };
+ arb::sample_tree sm(samples, parents);
+ {
+ arb::morphology m(sm, false);
+
+ EXPECT_EQ((ll{}), minset(m, ll{}));
+ EXPECT_EQ((ll{{2,0.1}}), minset(m, ll{{2,0.1}}));
+ EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), minset(m, ll{{0,0.5}, {1,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}}), minset(m, ll{{0,0.5}}));
+ EXPECT_EQ((ll{{0,0}, {1,0}}), minset(m, ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
+ EXPECT_EQ((ll{{0,0}, {1,0.5}}), minset(m, ll{{0,0}, {0,0.5}, {1,0.5}, {2,0.5}}));
+ EXPECT_EQ((ll{{0,0}, {2,0.5}}), minset(m, ll{{0,0}, {0,0.5}, {2,0.5}}));
+ EXPECT_EQ((ll{{0,0}, {2,0.5}, {3,0}}), minset(m, ll{{0,0}, {0,0.5}, {2,0.5}, {3,0}, {3,1}}));
+ EXPECT_EQ((ll{{0,0}, {1,0}}), minset(m, ll{{0,0}, {0,0.5}, {1,0}, {2,0.5}, {3,0}, {3,1}}));
+ }
+ {
+ arb::morphology m(sm, true);
+
+ EXPECT_EQ((ll{}), minset(m, ll{}));
+ EXPECT_EQ((ll{{2,0.1}}), minset(m, ll{{2,0.1}}));
+ EXPECT_EQ((ll{{0,0.5}}), minset(m, ll{{0,0.5}, {1,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}}), minset(m, ll{{0,0.5}}));
+ EXPECT_EQ((ll{{0,0}}), minset(m, ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
+ EXPECT_EQ((ll{{1,0.5}, {3,0.1}, {4,0.5}}), minset(m, ll{{1,0.5}, {1,1}, {3,0.1}, {4,0.5}, {4,0.7}}));
+ }
+}
+
diff --git a/test/unit/test_piecewise.cpp b/test/unit/test_piecewise.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..03047be6c591964c5bdc28c9d0ac9c799de24db9
--- /dev/null
+++ b/test/unit/test_piecewise.cpp
@@ -0,0 +1,332 @@
+#include "../gtest.h"
+
+#include <list>
+
+#include "util/piecewise.hpp"
+#include "util/rangeutil.hpp"
+
+using namespace arb;
+using util::pw_elements;
+using util::pw_element;
+using util::pw_npos;
+
+TEST(piecewise, eq) {
+ pw_elements<int> p1((double[3]){1., 1.5, 2.}, (int[2]){3, 4});
+ pw_elements<int> p2(p1);
+ EXPECT_EQ(p2, p1);
+
+ pw_elements<int> p3((double[3]){1., 1.7, 2.}, (int[2]){3, 4});
+ EXPECT_NE(p3, p1);
+
+ pw_elements<int> p4((double[3]){1., 1.5, 2.}, (int[2]){3, 5});
+ EXPECT_NE(p4, p1);
+
+ pw_elements<int> p5(p1);
+ p5.push_back(2., 7);
+ EXPECT_NE(p5, p1);
+
+ pw_elements<> v1((double[3]){1., 1.5, 2.});
+ pw_elements<> v2(v1);
+ EXPECT_EQ(v2, v1);
+
+ pw_elements<> v3((double[3]){1., 1.7, 2.});
+ EXPECT_NE(v3, v1);
+
+ pw_elements<> v5(v1);
+ v5.push_back(2.);
+ EXPECT_NE(v5, v1);
+}
+
+TEST(piecewise, generic_ctors) {
+ pw_elements<int> p1({1., 1.5, 2.}, {3, 4});
+ pw_elements<int> p2(std::list<float>{1.f, 1.5f, 2.f}, std::list<long>{3, 4});
+
+ EXPECT_EQ(p1, p2);
+
+ pw_elements<> v1(p1);
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.}), v1.vertices());
+
+ pw_elements<> v3(std::list<float>{1.f, 1.5f, 2.f});
+ EXPECT_EQ(v1, v3);
+
+ EXPECT_THROW(pw_elements<>({2.}), std::runtime_error);
+ EXPECT_THROW(pw_elements<int>({2.}, {1}), std::runtime_error);
+ EXPECT_THROW(pw_elements<int>({1., 2.}, {}), std::runtime_error);
+ EXPECT_THROW(pw_elements<int>({1., 2.}, {1, 3}), std::runtime_error);
+}
+
+TEST(piecewise, size) {
+ pw_elements<int> p1;
+ EXPECT_EQ(0u, p1.size());
+ EXPECT_TRUE(p1.empty());
+
+ pw_elements<int> p2({1., 2.}, {1});
+ EXPECT_EQ(1u, p2.size());
+ EXPECT_FALSE(p2.empty());
+
+ pw_elements<> v1;
+ EXPECT_EQ(0u, v1.size());
+ EXPECT_TRUE(v1.empty());
+
+ pw_elements<> v2({1., 2.});
+ EXPECT_EQ(1u, v2.size());
+ EXPECT_FALSE(v2.empty());
+}
+
+TEST(piecewise, assign) {
+ pw_elements<int> p;
+
+ double v[5] = {1., 1.5, 2., 2.5, 3.};
+ int x[4] = {10, 8, 9, 4};
+ p.assign(v, x);
+
+ ASSERT_EQ(4u, p.size());
+
+ EXPECT_EQ(10, p[0].second);
+ EXPECT_EQ( 8, p[1].second);
+ EXPECT_EQ( 9, p[2].second);
+ EXPECT_EQ( 4, p[3].second);
+
+ using dp = std::pair<double, double>;
+ EXPECT_EQ(dp(1.0, 1.5), p.interval(0));
+ EXPECT_EQ(dp(1.5, 2.0), p.interval(1));
+ EXPECT_EQ(dp(2.0, 2.5), p.interval(2));
+ EXPECT_EQ(dp(2.5, 3.0), p.interval(3));
+
+ pw_elements<int> q1(p);
+ pw_elements<int> q2;
+ q2 = p;
+ pw_elements<int> q3(p);
+ q3.assign(p.vertices(), p.elements());
+
+ EXPECT_EQ((std::vector<double>{1.0, 1.5, 2.0, 2.5, 3.0}), p.vertices());
+ EXPECT_EQ((std::vector<int>{10, 8, 9, 4}), p.elements());
+
+ EXPECT_EQ(q1.vertices(), p.vertices());
+ EXPECT_EQ(q2.vertices(), p.vertices());
+ EXPECT_EQ(q3.vertices(), p.vertices());
+
+ EXPECT_EQ(q1.elements(), p.elements());
+ EXPECT_EQ(q2.elements(), p.elements());
+ EXPECT_EQ(q3.elements(), p.elements());
+
+ q3.assign({}, {});
+ EXPECT_TRUE(q3.empty());
+
+ q3.assign({1., 2.}, {3});
+ pw_elements<int> q4(q3);
+ ASSERT_EQ(q3, q4);
+
+ // bad assign should throw but preserve value
+ ASSERT_THROW(q4.assign({3.}, {}), std::runtime_error);
+ EXPECT_EQ(q3, q4);
+}
+
+TEST(piecewise, assign_void) {
+ pw_elements<> p;
+
+ double v[5] = {1., 1.5, 2., 2.5, 3.};
+ p.assign(v);
+
+ ASSERT_EQ(4u, p.size());
+
+ using dp = std::pair<double, double>;
+ EXPECT_EQ(dp(1.0, 1.5), p.interval(0));
+ EXPECT_EQ(dp(1.5, 2.0), p.interval(1));
+ EXPECT_EQ(dp(2.0, 2.5), p.interval(2));
+ EXPECT_EQ(dp(2.5, 3.0), p.interval(3));
+
+ pw_elements<> q1(p);
+ pw_elements<> q2;
+ q2 = p;
+ pw_elements<> q3(p);
+ q3.assign(p.vertices());
+
+ EXPECT_EQ((std::vector<double>{1.0, 1.5, 2.0, 2.5, 3.0}), p.vertices());
+
+ EXPECT_EQ(q1.vertices(), p.vertices());
+ EXPECT_EQ(q2.vertices(), p.vertices());
+ EXPECT_EQ(q3.vertices(), p.vertices());
+
+ q3.assign({});
+ EXPECT_TRUE(q3.empty());
+
+ q3.assign({1., 2.});
+ pw_elements<> q4(q3);
+ ASSERT_EQ(q3, q4);
+
+ // bad assign should throw but preserve value
+ ASSERT_THROW(q4.assign({3.}), std::runtime_error);
+ EXPECT_EQ(q3, q4);
+}
+
+TEST(piecewise, access) {
+ pw_elements<int> p;
+
+ double v[5] = {1., 1.5, 2., 2.5, 3.};
+ int x[4] = {10, 8, 9, 4};
+ p.assign(v, x);
+
+ for (unsigned i = 0; i<4; ++i) {
+ EXPECT_EQ(v[i], p[i].first.first);
+ EXPECT_EQ(v[i+1], p[i].first.second);
+
+ EXPECT_EQ(v[i], p.interval(i).first);
+ EXPECT_EQ(v[i+1], p.interval(i).second);
+
+ EXPECT_EQ(x[i], p[i].second);
+ EXPECT_EQ(x[i], p.element(i));
+ }
+
+ EXPECT_EQ(p[0], p.front());
+ EXPECT_EQ(p[3], p.back());
+
+ unsigned j = 0;
+ for (auto entry: p) {
+ EXPECT_EQ(p[j++], entry);
+ }
+}
+
+TEST(piecewise, bounds) {
+ pw_elements<int> p{{1., 1.5, 2., 2.5, 3.}, {10, 8, 9, 4}};
+
+ EXPECT_EQ(1., p.bounds().first);
+ EXPECT_EQ(3., p.bounds().second);
+
+ pw_elements<> v{{1., 1.5, 2., 2.5, 3.}};
+
+ EXPECT_EQ(1., v.bounds().first);
+ EXPECT_EQ(3., v.bounds().second);
+}
+
+TEST(piecewise, index_of) {
+ pw_elements<int> p{{1., 1.5, 2., 2.5, 3.}, {10, 8, 9, 4}};
+
+ EXPECT_EQ(pw_npos, p.index_of(0.3));
+ EXPECT_EQ(0u, p.index_of(1.));
+ EXPECT_EQ(0u, p.index_of(1.1));
+ EXPECT_EQ(1u, p.index_of(1.5));
+ EXPECT_EQ(1u, p.index_of(1.6));
+ EXPECT_EQ(2u, p.index_of(2));
+ EXPECT_EQ(3u, p.index_of(2.9));
+ EXPECT_EQ(3u, p.index_of(3));
+ EXPECT_EQ(pw_npos, p.index_of(3.1));
+
+ pw_elements<> v(p);
+
+ EXPECT_EQ(pw_npos, v.index_of(0.3));
+ EXPECT_EQ(0u, v.index_of(1.));
+ EXPECT_EQ(0u, v.index_of(1.1));
+ EXPECT_EQ(1u, v.index_of(1.5));
+ EXPECT_EQ(1u, v.index_of(1.6));
+ EXPECT_EQ(2u, v.index_of(2));
+ EXPECT_EQ(3u, v.index_of(2.9));
+ EXPECT_EQ(3u, v.index_of(3));
+ EXPECT_EQ(pw_npos, v.index_of(3.1));
+
+ pw_elements<int> p0;
+ pw_elements<> v0;
+
+ EXPECT_EQ(pw_npos, p0.index_of(0.));
+ EXPECT_EQ(pw_npos, v0.index_of(0.));
+}
+
+TEST(piecewise, push) {
+ pw_elements<int> q;
+ using dp = std::pair<double, double>;
+
+ // Need left hand side!
+ EXPECT_THROW(q.push_back(3.1, 4), std::runtime_error);
+
+ q.clear();
+ q.push_back(1.1, 3.1, 4);
+ q.push_back(3.1, 4.3, 5);
+ EXPECT_EQ(dp(1.1, 3.1), q.interval(0));
+ EXPECT_EQ(dp(3.1, 4.3), q.interval(1));
+ EXPECT_EQ(4, q[0].second);
+ EXPECT_EQ(5, q[1].second);
+
+ q.push_back(7.2, 6);
+ EXPECT_EQ(dp(4.3, 7.2), q.interval(2));
+ EXPECT_EQ(6, q[2].second);
+
+ // Supplied left side doesn't match current right.
+ EXPECT_THROW(q.push_back(7.4, 9.1, 7), std::runtime_error);
+}
+
+TEST(piecewise, push_void) {
+ pw_elements<> p;
+ using dp = std::pair<double, double>;
+
+ // Need left hand side!
+ EXPECT_THROW(p.push_back(3.1), std::runtime_error);
+
+ p.clear();
+ p.push_back(0.1, 0.2);
+ p.push_back(0.2, 0.4);
+ p.push_back(0.5);
+
+ EXPECT_EQ(3u, p.size());
+ EXPECT_EQ((std::vector<double>{0.1, 0.2, 0.4, 0.5}), p.vertices());
+ EXPECT_EQ(dp(0.2,0.4), p.interval(1));
+
+ // Supplied left side doesn't match current right.
+ EXPECT_THROW(p.push_back(0.7, 0.9), std::runtime_error);
+}
+
+TEST(piecewise, zip) {
+ pw_elements<int> p03;
+ p03.assign((double [3]){0., 1.5, 3.}, (int [2]){10, 11});
+
+ pw_elements<int> p14;
+ p14.assign((double [5]){1, 2.25, 3., 3.5, 4.}, (int [4]){3, 4, 5, 6});
+
+ using ii = std::pair<int, int>;
+ pw_elements<ii> p03_14 = zip(p03, p14);
+ EXPECT_EQ(1., p03_14.bounds().first);
+ EXPECT_EQ(3., p03_14.bounds().second);
+
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), p03_14.vertices());
+ EXPECT_EQ((std::vector<ii>{ii(10, 3), ii(11, 3), ii(11, 4)}), p03_14.elements());
+
+ pw_elements<ii> p14_03 = zip(p14, p03);
+ EXPECT_EQ(p03_14.vertices(), p14_03.vertices());
+
+ std::vector<ii> flipped = util::assign_from(util::transform_view(p14_03.elements(),
+ [](ii p) { return ii{p.second, p.first}; }));
+ EXPECT_EQ(p03_14.elements(), flipped);
+
+ pw_elements<> v03;
+ v03.assign((double [3]){0., 1.5, 3.});
+
+ EXPECT_EQ((std::vector<int>{3, 3, 4}), zip(v03, p14).elements());
+ EXPECT_EQ((std::vector<int>{3, 3, 4}), zip(p14, v03).elements());
+
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), zip(v03, p14).vertices());
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), zip(p14, v03).vertices());
+
+ auto project = [](double l, double r, pw_element<void>, const pw_element<int>& b) -> double {
+ double b_width = b.first.second-b.first.first;
+ return b.second*(r-l)/b_width;
+ };
+
+ pw_elements<void> vxx; // elements cover bounds of p14
+ vxx.assign((double [6]){0.2, 1.7, 1.95, 2.325, 2.45, 4.9});
+
+ pw_elements<double> pxx = zip(vxx, p14, project);
+ double p14_sum = util::sum(util::transform_view(p14, [](auto v) { return v.second; }));
+ double pxx_sum = util::sum(util::transform_view(pxx, [](auto v) { return v.second; }));
+ EXPECT_DOUBLE_EQ(p14_sum, pxx_sum);
+
+}
+
+TEST(piecewise, zip_void) {
+ pw_elements<> p03;
+ p03.assign((double [3]){0., 1.5, 3.});
+
+ pw_elements<> p14;
+ p14.assign((double [5]){1, 2.25, 3., 3.5, 4.});
+
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), zip(p03, p14).vertices());
+ EXPECT_EQ((std::vector<double>{1., 1.5, 2.25, 3.}), zip(p14, p03).vertices());
+}
diff --git a/test/unit/test_range.cpp b/test/unit/test_range.cpp
index b15c86abb54523ed7c38a262c822d1d2752c3fad..39a80f45ba027f057fffb3c09dc4a1d06797bd31 100644
--- a/test/unit/test_range.cpp
+++ b/test/unit/test_range.cpp
@@ -478,6 +478,26 @@ TEST(range, sort) {
EXPECT_EQ(X, (std::vector<foo>{{0, 5}, {1, 4}, {2, 3}, {3, 2}, {4, 1}, {5, 0}}));
}
+TEST(range, sum) {
+ std::string words[] = { "fish", "cakes", "!" };
+
+ auto result = util::sum(words);
+ EXPECT_EQ("fishcakes!", result);
+
+ result = util::sum(words, "tasty"s);
+ EXPECT_EQ("tastyfishcakes!", result);
+
+ struct uwrap {
+ unsigned value = 0;
+ uwrap(unsigned v): value(v) {}
+
+ uwrap operator+(const std::string& s) { return value+s.size(); }
+ };
+
+ auto count = util::sum(words, uwrap{3});
+ EXPECT_EQ(3u+4u+5u+1u, count.value);
+}
+
TEST(range, sum_by) {
std::string words[] = { "fish", "cakes", "!" };
auto prepend_ = [](const std::string& x) { return "_"+x; };
diff --git a/test/unit/test_ratelem.cpp b/test/unit/test_ratelem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f393ae75456522ea83a7320c4d643f650baef7a6
--- /dev/null
+++ b/test/unit/test_ratelem.cpp
@@ -0,0 +1,115 @@
+#include "../gtest.h"
+
+#include <list>
+
+#include "util/ratelem.hpp"
+#include "util/rangeutil.hpp"
+
+using namespace arb;
+using util::rat_element;
+
+TEST(ratelem, direct_ctor) {
+ rat_element<0, 0> x00(3.5);
+ EXPECT_EQ(1u, x00.size());
+ EXPECT_EQ(3.5, x00[0]);
+
+ rat_element<1, 3> x13(1.1, 2.2, 3.3, 4.4, 5.5);
+ EXPECT_EQ(5u, x13.size());
+ EXPECT_EQ(1.1, x13[0]);
+ EXPECT_EQ(2.2, x13[1]);
+ EXPECT_EQ(3.3, x13[2]);
+ EXPECT_EQ(4.4, x13[3]);
+ EXPECT_EQ(5.5, x13[4]);
+
+ std::array<float, 4> x21_arr{1.25f, 1.5f, 0.5f, 2.25f};
+ rat_element<2, 1> x21(x21_arr);
+ EXPECT_EQ(4u, x21.size());
+ EXPECT_EQ(1.25, x21[0]);
+ EXPECT_EQ(1.5, x21[1]);
+ EXPECT_EQ(0.5, x21[2]);
+ EXPECT_EQ(2.25, x21[3]);
+
+ int x20_arr[3] = {3, 2, 4};
+ rat_element<2, 0> x20(x20_arr);
+ EXPECT_EQ(3u, x20.size());
+ EXPECT_EQ(3., x20[0]);
+ EXPECT_EQ(2., x20[1]);
+ EXPECT_EQ(4., x20[2]);
+}
+
+TEST(ratelem, fn_ctor) {
+ auto f = [](double x) { return 1+x*x; };
+
+ rat_element<0, 0> x00(f);
+ EXPECT_EQ(1., x00[0]);
+
+ rat_element<1, 3> x12(f);
+ EXPECT_EQ(f(0.00), x12[0]);
+ EXPECT_EQ(f(0.25), x12[1]);
+ EXPECT_EQ(f(0.50), x12[2]);
+ EXPECT_EQ(f(0.75), x12[3]);
+ EXPECT_EQ(f(1.00), x12[4]);
+}
+
+
+TEST(ratelem, constants) {
+ // (Only expect this to work for polynomial interpolators).
+ auto k = [](double c) { return [c](double x) { return c; }; };
+
+ rat_element<0, 0> k00(k(2.));
+ rat_element<1, 0> k10(k(3.));
+ rat_element<2, 0> k20(k(4.));
+ rat_element<3, 0> k30(k(5.));
+
+ double xs[] = {0., 0.1, 0.3, 0.5, 0.7, 0.9, 1.};
+ for (auto x: xs) {
+ EXPECT_DOUBLE_EQ(2., k00(x));
+ EXPECT_DOUBLE_EQ(3., k10(x));
+ EXPECT_DOUBLE_EQ(4., k20(x));
+ EXPECT_DOUBLE_EQ(5., k30(x));
+ }
+}
+
+template <unsigned p_, unsigned q_>
+struct wrap_pq {
+ static constexpr unsigned p = p_;
+ static constexpr unsigned q = q_;
+};
+
+template <typename PQ>
+struct ratelem_pq: public testing::Test {};
+
+using pq_types = ::testing::Types<
+ wrap_pq<0, 0>, wrap_pq<1, 0>, wrap_pq<2, 0>, wrap_pq<3, 0>,
+ wrap_pq<0, 1>, wrap_pq<1, 1>, wrap_pq<2, 1>, wrap_pq<3, 1>,
+ wrap_pq<0, 2>, wrap_pq<1, 2>, wrap_pq<2, 2>, wrap_pq<3, 2>
+>;
+
+// Compute 1+x+x^2+...+x^n
+template <unsigned n>
+constexpr double upoly(double x) { return x*upoly<n-1>(x)+1.; }
+
+template <>
+constexpr double upoly<0>(double x) { return 1.; }
+
+TYPED_TEST_CASE(ratelem_pq, pq_types);
+
+TYPED_TEST(ratelem_pq, interpolate_monotonic) {
+ constexpr unsigned p = TypeParam::p;
+ constexpr unsigned q = TypeParam::q;
+
+ // Pick f to have irreducible order p, q.
+ auto f = [](double x) { return upoly<p>(x)/upoly<q>(2*x); };
+ rat_element<p, q> fpq(f);
+
+ // Check values both on and off nodes.
+ for (unsigned i = 0; i<=1+p+q; ++i) {
+ double x = (double)i/(1+p+q);
+ EXPECT_DOUBLE_EQ(f(x), fpq(x));
+ }
+
+ for (unsigned i = 1; i<p+q; ++i) {
+ double x = (double)i/(p+q);
+ EXPECT_DOUBLE_EQ(f(x), fpq(x));
+ }
+}