diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index e7d9740df549268667da5a7c4b21c5fcf6b9c12a..c8e7c339d6ad689183977fa2e6660fdcce7afd82 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -28,9 +28,14 @@ set(arbor_sources
mechcat.cpp
memory/cuda_wrappers.cpp
memory/util.cpp
+ morph/em_morphology.cpp
+ morph/label_dict.cpp
+ morph/locset.cpp
+ morph/mbranch.cpp
morph/morphology.cpp
- morph/sample_tree.cpp
morph/primitives.cpp
+ morph/region.cpp
+ morph/sample_tree.cpp
merge_events.cpp
simulation.cpp
partition_load_balance.cpp
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index 50ca10c05c968e6200e51d7823b402c618797608..90dcd74438148e6adf4d9cd262af4ea004303a21 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -1,9 +1,15 @@
+#include <sstream>
+#include <vector>
+
#include <arbor/cable_cell.hpp>
+#include <arbor/morph/label_dict.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/segment.hpp>
+#include "morph/em_morphology.hpp"
#include "util/rangeutil.hpp"
#include "util/span.hpp"
+#include "util/strprintf.hpp"
namespace arb {
@@ -72,6 +78,17 @@ bool cable_cell::has_soma() const {
return !segment(0)->is_placeholder();
}
+void cable_cell::paint(const std::string& target, const std::string& description) {
+ auto it = regions_.find(target);
+
+ // Nothing to do if there are no regions that match.
+ if (it==regions_.end()) return;
+
+ for (auto i: it->second) {
+ segment(i)->add_mechanism(description);
+ }
+}
+
soma_segment* cable_cell::soma() {
return has_soma()? segment(0)->as_soma(): nullptr;
}
@@ -106,12 +123,12 @@ size_type cable_cell::num_compartments() const {
[](const segment_ptr& s) { return s->num_compartments(); });
}
-void cable_cell::add_stimulus(segment_location loc, i_clamp stim) {
- (void)segment(loc.segment); // assert loc.segment in range
+void cable_cell::add_stimulus(mlocation loc, i_clamp stim) {
+ (void)segment(loc.branch); // assert loc.segment in range
stimuli_.push_back({loc, std::move(stim)});
}
-void cable_cell::add_detector(segment_location loc, double threshold) {
+void cable_cell::add_detector(mlocation loc, double threshold) {
spike_detectors_.push_back({loc, threshold});
}
@@ -164,24 +181,27 @@ value_type cable_cell::segment_mean_attenuation(
return 2*std::sqrt(math::pi<double>*tau_per_um*frequency)*length_factor; // [1/µm]
}
-cable_cell make_cable_cell(const morphology& morph, bool compartments_from_discretization) {
+cable_cell make_cable_cell(const morphology& m,
+ const label_dict& dictionary,
+ bool compartments_from_discretization)
+{
using point3d = cable_cell::point_type;
cable_cell newcell;
- if (!morph.num_branches()) {
+ if (!m.num_branches()) {
return newcell;
}
// Add the soma.
- auto loc = morph.samples()[0].loc; // location of soma.
+ auto loc = m.samples()[0].loc; // location of soma.
// If there is no spherical root/soma use a zero-radius soma.
- double srad = morph.spherical_root()? loc.radius: 0.;
+ double srad = m.spherical_root()? loc.radius: 0.;
newcell.add_soma(srad, point3d(loc.x, loc.y, loc.z));
- auto& samples = morph.samples();
- for (auto i: util::make_span(1, morph.num_branches())) {
- auto index = util::make_range(morph.branch_sample_span(i));
+ auto& samples = m.samples();
+ for (auto i: util::make_span(1, m.num_branches())) {
+ auto index = util::make_range(m.branch_indexes(i));
// find kind for the branch. Use the tag of the last sample in the branch.
int tag = samples[index.back()].tag;
@@ -207,9 +227,9 @@ cable_cell make_cable_cell(const morphology& morph, bool compartments_from_discr
}
// Find the id of this branch's parent.
- auto pid = morph.branch_parent(i);
+ auto pid = m.branch_parent(i);
// Adjust pid if a zero-radius soma was used.
- if (!morph.spherical_root()) {
+ if (!m.spherical_root()) {
pid = pid==mnpos? 0: pid+1;
}
auto cable = newcell.add_cable(pid, kind, radii, points);
@@ -218,6 +238,26 @@ cable_cell make_cable_cell(const morphology& morph, bool compartments_from_discr
}
}
+ // Construct concrete regions.
+ // Ignores the pointsets, for now.
+ auto em = em_morphology(m); // for converting labels to "segments"
+
+ std::unordered_map<std::string, std::vector<msize_t>> regions;
+ for (auto r: dictionary.regions()) {
+ mcable_list L = thingify(r.second, em);
+ std::vector<msize_t> bids;
+ for (auto c: L) {
+ if (c.prox_pos!=0 || c.dist_pos!=1) {
+ throw cable_cell_error(util::pprintf(
+ "cable_cell does not support regions with partial branches: \"{}\": {}",
+ r.first, c));
+ }
+ bids.push_back(c.branch);
+ }
+ regions[r.first] = bids;
+ }
+ newcell.set_regions(std::move(regions));
+
return newcell;
}
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index f035b70b88f8b2290a281a8af1f5a3ce4d08d8c8..aa0e9b3af846fcc9a2b0011432d7d3ab68c7eeb9 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -612,15 +612,15 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
for (const auto& cellsyn: cell.synapses()) {
const mechanism_desc& desc = cellsyn.mechanism;
- size_type cv = D.segment_location_cv(cell_idx, cellsyn.location);
+ size_type cv = D.branch_location_cv(cell_idx, cellsyn.location);
const auto& name = desc.name();
point_mech_data& entry = point_mech_table[name];
update_paramset_and_validate(desc, entry.info, entry.paramset);
entry.points.push_back({cv, target_id++, &desc});
- const segment_ptr& seg = cell.segments()[cellsyn.location.segment];
- size_type segment_idx = D.cell_segment_bounds[cell_idx]+cellsyn.location.segment;
+ const segment_ptr& seg = cell.segments()[cellsyn.location.branch];
+ size_type segment_idx = D.cell_segment_bounds[cell_idx]+cellsyn.location.branch;
for (const auto& ion_entry: entry.info->ions) {
const std::string& ion_name = ion_entry.first;
@@ -635,7 +635,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
}
for (const auto& stimulus: cell.stimuli()) {
- size_type cv = D.segment_location_cv(cell_idx, stimulus.location);
+ size_type cv = D.branch_location_cv(cell_idx, stimulus.location);
stimuli.push_back({cv, stimulus.clamp});
}
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index c0172e06e21ba10a15cba1c117a6bcc1d6298b70..879be39b806e9e2cf455655eb988d1087a0786c6 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -85,12 +85,12 @@ struct fvm_discretization {
return util::partition_view(cell_cv_bounds);
}
- size_type segment_location_cv(size_type cell_index, segment_location segloc) const {
+ size_type branch_location_cv(size_type cell_index, mlocation loc) const {
auto cell_segs = cell_segment_part()[cell_index];
- size_type seg = segloc.segment+cell_segs.first;
+ size_type seg = loc.branch+cell_segs.first;
arb_assert(seg<cell_segs.second);
- return segments[seg].cv_by_position(segloc.position);
+ return segments[seg].cv_by_position(loc.pos);
}
};
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index a466077b19611dcc2132e511a8c49c1ab7f30290..b016793e68aa2c559161d1170a5fd7f2c468e0df 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -507,7 +507,7 @@ void fvm_lowered_cell_impl<B>::initialize(
cell_gid_type gid = gids[cell_idx];
for (auto detector: cells[cell_idx].detectors()) {
- detector_cv.push_back(D.segment_location_cv(cell_idx, detector.location));
+ detector_cv.push_back(D.branch_location_cv(cell_idx, detector.location));
detector_threshold.push_back(detector.threshold);
}
@@ -515,7 +515,7 @@ void fvm_lowered_cell_impl<B>::initialize(
probe_info pi = rec.get_probe({gid, j});
auto where = any_cast<cell_probe_address>(pi.address);
- auto cv = D.segment_location_cv(cell_idx, where.location);
+ auto cv = D.branch_location_cv(cell_idx, where.location);
probe_handle handle;
switch (where.kind) {
@@ -555,7 +555,7 @@ std::vector<fvm_gap_junction> fvm_lowered_cell_impl<B>::fvm_gap_junctions(
auto cell_gj = cells[cell_idx].gap_junction_sites();
for (auto gj : cell_gj) {
- auto cv = D.segment_location_cv(cell_idx, gj);
+ auto cv = D.branch_location_cv(cell_idx, gj);
gid_to_cvs[gids[cell_idx]].push_back(cv);
}
}
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index 09baa5de0507dc6a8e28f51e3f59112892b45ba0..2303aaff7f63e70277b005b0c38e7fd1f9c307ea 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -9,28 +9,13 @@
#include <arbor/common_types.hpp>
#include <arbor/constants.hpp>
#include <arbor/mechcat.hpp>
+#include <arbor/morph/label_dict.hpp>
#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
#include <arbor/segment.hpp>
namespace arb {
-// Location specification for point processes.
-
-struct segment_location {
- segment_location(cell_lid_type s, double l):
- segment(s), position(l)
- {
- arb_assert(position>=0. && position<=1.);
- }
-
- bool operator==(segment_location other) const {
- return segment==other.segment && position==other.position;
- }
-
- cell_lid_type segment;
- double position;
-};
-
// Current clamp description for stimulus specification.
struct i_clamp {
@@ -52,7 +37,7 @@ struct cell_probe_address {
membrane_voltage, membrane_current
};
- segment_location location;
+ mlocation location;
probe_kind kind;
};
@@ -64,20 +49,22 @@ public:
using value_type = double;
using point_type = point<value_type>;
- using gap_junction_instance = segment_location;
+ using gap_junction_instance = mlocation;
+
+ using region_map = std::unordered_map<std::string, std::vector<msize_t>>;
struct synapse_instance {
- segment_location location;
+ mlocation location;
mechanism_desc mechanism;
};
struct stimulus_instance {
- segment_location location;
+ mlocation location;
i_clamp clamp;
};
struct detector_instance {
- segment_location location;
+ mlocation location;
double threshold;
};
@@ -93,7 +80,8 @@ public:
stimuli_(other.stimuli_),
synapses_(other.synapses_),
gap_junction_sites_(other.gap_junction_sites_),
- spike_detectors_(other.spike_detectors_)
+ spike_detectors_(other.spike_detectors_),
+ regions_(other.regions_)
{
// unique_ptr's cannot be copy constructed, do a manual assignment
segments_.reserve(other.segments_.size());
@@ -158,7 +146,7 @@ public:
//////////////////
// stimuli
//////////////////
- void add_stimulus(segment_location loc, i_clamp stim);
+ void add_stimulus(mlocation loc, i_clamp stim);
std::vector<stimulus_instance>&
stimuli() {
@@ -170,10 +158,15 @@ public:
return stimuli_;
}
+ //////////////////
+ // painters
+ //////////////////
+ void paint(const std::string& target, const std::string& description);
+
//////////////////
// synapses
//////////////////
- void add_synapse(segment_location loc, mechanism_desc p)
+ void add_synapse(mlocation loc, mechanism_desc p)
{
synapses_.push_back(synapse_instance{loc, std::move(p)});
}
@@ -184,7 +177,7 @@ public:
//////////////////
// gap-junction
//////////////////
- void add_gap_junction(segment_location location)
+ void add_gap_junction(mlocation location)
{
gap_junction_sites_.push_back(location);
}
@@ -195,7 +188,7 @@ public:
//////////////////
// spike detectors
//////////////////
- void add_detector(segment_location loc, double threshold);
+ void add_detector(mlocation loc, double threshold);
std::vector<detector_instance>&
detectors() {
@@ -207,6 +200,10 @@ public:
return spike_detectors_;
}
+ void set_regions(region_map r) {
+ regions_ = std::move(r);
+ }
+
// Checks that two cells have the same
// - number and type of segments
// - volume and area properties of each segment
@@ -253,6 +250,9 @@ private:
// the sensors
std::vector<detector_instance> spike_detectors_;
+
+ // Named regions, oh my.
+ region_map regions_;
};
// create a cable by forwarding cable construction parameters provided by the user
@@ -273,6 +273,8 @@ cable_segment* cable_cell::add_cable(cable_cell::index_type parent, Args&&... ar
// If compartments_from_discretization is true, set number of compartments
// in each segment to be the number of piecewise linear sections in the
// corresponding section of the morphology.
-cable_cell make_cable_cell(const morphology& morph, bool compartments_from_discretization);
+cable_cell make_cable_cell(const morphology& morph,
+ const label_dict& labels={},
+ bool compartments_from_discretization=false);
} // namespace arb
diff --git a/arbor/include/arbor/morph/label_dict.hpp b/arbor/include/arbor/morph/label_dict.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..852cef7121fb4d1a3a6b69b3bb1a71d7e1287a50
--- /dev/null
+++ b/arbor/include/arbor/morph/label_dict.hpp
@@ -0,0 +1,32 @@
+#pragma once
+
+#include <memory>
+#include <unordered_map>
+
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/util/optional.hpp>
+
+namespace arb {
+
+class label_dict {
+ using ps_map = std::unordered_map<std::string, arb::locset>;
+ using reg_map = std::unordered_map<std::string, arb::region>;
+ ps_map locsets_;
+ reg_map regions_;
+
+public:
+
+ void set(const std::string& name, locset ls);
+ void set(const std::string& name, region reg);
+
+ util::optional<const arb::region&> region(const std::string& name) const;
+ util::optional<const arb::locset&> locset(const std::string& name) const;
+
+ const ps_map& locsets() const;
+ const reg_map& regions() const;
+
+ std::size_t size() const;
+};
+
+} //namespace arb
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..00c7e97386025aa9dbbeea782bb0ed178e01733a
--- /dev/null
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -0,0 +1,122 @@
+#pragma once
+
+#include <cstdlib>
+#include <memory>
+#include <set>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/morphology.hpp>
+
+namespace arb {
+
+// Forward declare the backend em_morphology type, required for defining the
+// interface for concretising locsets.
+class em_morphology;
+
+class locset {
+public:
+ locset() = delete;
+
+ template <typename Impl,
+ typename X=std::enable_if_t<!std::is_same<std::decay_t<Impl>, locset>::value>>
+ explicit locset(Impl&& impl):
+ impl_(new wrap<Impl>(std::forward<Impl>(impl))) {}
+
+ template <typename Impl>
+ explicit locset(const Impl& impl):
+ impl_(new wrap<Impl>(impl)) {}
+
+ locset(locset&& other) = default;
+
+ locset(const locset& other):
+ impl_(other.impl_->clone()) {}
+
+ locset& operator=(const locset& other) {
+ impl_ = other.impl_->clone();
+ return *this;
+ }
+
+ template <typename Impl,
+ typename X=std::enable_if_t<!std::is_same<std::decay_t<Impl>, locset>::value>>
+ locset& operator=(Impl&& other) {
+ impl_ = new wrap<Impl>(std::forward<Impl>(other));
+ return *this;
+ }
+
+ template <typename Impl>
+ locset& operator=(const Impl& other) {
+ impl_ = new wrap<Impl>(other);
+ return *this;
+ }
+
+ friend mlocation_list thingify(const locset& p, const em_morphology& m) {
+ return p.impl_->thingify(m);
+ }
+
+ friend std::ostream& operator<<(std::ostream& o, const locset& p) {
+ return p.impl_->print(o);
+ }
+
+ // The sum of two location sets.
+ friend locset sum(locset, locset);
+
+ template <typename ...Args>
+ friend locset sum(locset l, locset r, Args... args) {
+ return sum(sum(std::move(l), std::move(r)), std::move(args)...);
+ }
+
+private:
+ struct interface {
+ virtual ~interface() {}
+ virtual std::unique_ptr<interface> clone() = 0;
+ virtual std::ostream& print(std::ostream&) = 0;
+ virtual mlocation_list thingify(const em_morphology&) = 0;
+ };
+
+ std::unique_ptr<interface> impl_;
+
+ template <typename Impl>
+ struct wrap: interface {
+ explicit wrap(const Impl& impl): wrapped(impl) {}
+ explicit wrap(Impl&& impl): wrapped(std::move(impl)) {}
+
+ virtual std::unique_ptr<interface> clone() override {
+ return std::unique_ptr<interface>(new wrap<Impl>(wrapped));
+ }
+
+ virtual mlocation_list thingify(const em_morphology& m) override {
+ return thingify_(wrapped, m);
+ }
+
+ virtual std::ostream& print(std::ostream& o) override {
+ return o << wrapped;
+ }
+
+ Impl wrapped;
+ };
+};
+
+namespace ls {
+
+// Location of a sample.
+locset location(mlocation);
+
+// Location of a sample.
+locset sample(msize_t);
+
+// Set of terminal nodes on a morphology.
+locset terminal();
+
+// The root node of a morphology.
+locset root();
+
+// The null (empty) set.
+locset nil();
+
+} // namespace ls
+
+} // namespace arb
diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index 4ad07955f68af5ea976bdb8d89ecdf74fbffa280..6eed5ef9abc081b61961dc7fc1f849c3e650f30b 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <memory>
#include <ostream>
#include <vector>
@@ -9,38 +10,13 @@
namespace arb {
-// An unbranched cable segment that has root, terminal or fork point at each end.
-struct mbranch {
- std::vector<msize_t> index; // sample index
- msize_t parent_id = mnpos; // branch index
+class morphology_impl;
- mbranch() = default;
- mbranch(std::vector<msize_t> idx, msize_t parent):
- index(std::move(idx)), parent_id(parent) {}
-
- bool is_sphere() const { return size()==1u; }
- msize_t size() const { return index.size(); }
- bool has_parent() const { return parent_id!=mnpos;}
-
- friend bool operator==(const mbranch& l, const mbranch& r);
- friend std::ostream& operator<<(std::ostream& o, const mbranch& b);
-};
+using mindex_range = std::pair<const msize_t*, const msize_t*>;
class morphology {
- // The sample tree of sample points and their parent-child relationships.
- sample_tree samples_;
-
- // Indicates whether the soma is a sphere.
- bool spherical_root_ = false;
-
- // Branch state.
- std::vector<mbranch> branches_;
- std::vector<msize_t> branch_parents_;
- std::vector<std::vector<msize_t>> branch_children_;
-
- using index_range = std::pair<const msize_t*, const msize_t*>;
-
- void init();
+ // Hold an immutable copy of the morphology implementation.
+ std::shared_ptr<const morphology_impl> impl_;
public:
morphology(sample_tree m, bool use_spherical_root);
@@ -52,8 +28,8 @@ public:
// The number of branches in the morphology.
msize_t num_branches() const;
- // The parent sample of sample i.
- const std::vector<msize_t>& sample_parents() const;
+ // The number of samples in the morphology.
+ msize_t num_samples() const;
// The parent branch of branch b.
msize_t branch_parent(msize_t b) const;
@@ -62,11 +38,17 @@ public:
const std::vector<msize_t>& branch_children(msize_t b) const;
// Range of indexes into the sample points in branch b.
- index_range branch_sample_span(msize_t b) const;
+ mindex_range branch_indexes(msize_t b) const;
- // Range of the samples in branch b.
+ // All of the samples in the morphology.
const std::vector<msample>& samples() const;
+ // The parent sample of sample i.
+ const std::vector<msize_t>& sample_parents() const;
+
+ // Point properties of samples in the morphology.
+ const std::vector<point_prop>& sample_props() const;
+
friend std::ostream& operator<<(std::ostream&, const morphology&);
};
diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp
index 42baa9b7f46c6d24d28df259fd9478cfb8daf60b..ae58c511f7eb2cb0afd076b01a3c5f5c7e1eb92a 100644
--- a/arbor/include/arbor/morph/primitives.hpp
+++ b/arbor/include/arbor/morph/primitives.hpp
@@ -48,6 +48,53 @@ struct msample {
bool is_collocated(const msample& a, const msample& b);
double distance(const msample& a, const msample& b);
+// Describe a specific location on a morpholology.
+struct mlocation {
+ // The id of the branch.
+ msize_t branch;
+ // The relative position on the branch ∈ [0,1].
+ double pos;
+
+ // branch ≠npos and 0 ≤ pos ≤ 1
+ friend bool test_invariants(const mlocation&);
+ friend std::ostream& operator<<(std::ostream&, const mlocation&);
+};
+
+ARB_DEFINE_LEXICOGRAPHIC_ORDERING(mlocation, (a.branch,a.pos), (b.branch,b.pos));
+
+using mlocation_list = std::vector<mlocation>;
+std::ostream& operator<<(std::ostream& o, const mlocation_list& l);
+
+// Tests whether each location in the list satisfies the invariants for a location,
+// and that the locations in the vector are ordered.
+bool test_invariants(const mlocation_list&);
+
+// Describe an unbranched cable in the morphology.
+struct mcable {
+ // The id of the branch on which the cable lies.
+ msize_t branch;
+
+ // Relative location of the end points on the branch.
+ // 0 ≤ prox_pos ≤ dist_pos ≤ 1
+ double prox_pos; // ∈ [0,1]
+ double dist_pos; // ∈ [0,1]
+
+ friend mlocation prox_loc(const mcable&);
+ friend mlocation dist_loc(const mcable&);
+
+ // branch ≠npos, and 0 ≤ prox_pos ≤ dist_pos ≤ 1
+ friend bool test_invariants(const mcable&);
+ friend std::ostream& operator<<(std::ostream&, const mcable&);
+};
+
+ARB_DEFINE_LEXICOGRAPHIC_ORDERING(mcable, (a.branch,a.prox_pos,a.dist_pos), (b.branch,b.prox_pos,b.dist_pos));
+
+using mcable_list = std::vector<mcable>;
+std::ostream& operator<<(std::ostream& o, const mcable_list& c);
+// Tests whether each cable in the list satisfies the invariants for a cable,
+// and that the cables in the vector are ordered.
+bool test_invariants(const mcable_list&);
+
using point_prop = std::uint8_t;
enum point_prop_mask: point_prop {
point_prop_mask_none = 0,
diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d3d7f735d553202da1070300be948a57619bed47
--- /dev/null
+++ b/arbor/include/arbor/morph/region.hpp
@@ -0,0 +1,129 @@
+#pragma once
+
+#include <cstdlib>
+#include <memory>
+#include <set>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/morphology.hpp>
+
+namespace arb {
+
+// Forward declare the backend em_morphology type, required for defining the
+// interface for concretising locsets.
+class em_morphology;
+
+class region {
+public:
+ region() = delete;
+
+ template <typename Impl,
+ typename X=std::enable_if_t<!std::is_same<std::decay_t<Impl>, region>::value>>
+ explicit region(Impl&& impl):
+ impl_(new wrap<Impl>(std::forward<Impl>(impl))) {}
+
+ template <typename Impl>
+ explicit region(const Impl& impl):
+ impl_(new wrap<Impl>(impl)) {}
+
+ region(region&& other) = default;
+
+ region(const region& other):
+ impl_(other.impl_->clone()) {}
+
+ region& operator=(const region& other) {
+ impl_ = other.impl_->clone();
+ return *this;
+ }
+
+ template <typename Impl,
+ typename X=std::enable_if_t<!std::is_same<std::decay_t<Impl>, region>::value>>
+ region& operator=(Impl&& other) {
+ impl_ = new wrap<Impl>(std::forward<Impl>(other));
+ return *this;
+ }
+
+ template <typename Impl>
+ region& operator=(const Impl& other) {
+ impl_ = new wrap<Impl>(other);
+ return *this;
+ }
+
+ friend mcable_list thingify(const region& r, const em_morphology& m) {
+ return r.impl_->thingify(m);
+ }
+
+ friend std::ostream& operator<<(std::ostream& o, const region& p) {
+ return p.impl_->print(o);
+ }
+
+ // The union of regions.
+ friend region join(region, region);
+
+ template <typename ...Args>
+ friend region join(region l, region r, Args... args) {
+ return join(join(std::move(l), std::move(r)), std::move(args)...);
+ }
+
+ // The intersection of regions.
+ friend region intersect(region, region);
+
+ template <typename ...Args>
+ friend region intersect(region l, region r, Args... args) {
+ return intersect(intersect(std::move(l), std::move(r)), std::move(args)...);
+ }
+
+private:
+ struct interface {
+ virtual ~interface() {}
+ virtual std::unique_ptr<interface> clone() = 0;
+ virtual std::ostream& print(std::ostream&) = 0;
+ virtual mcable_list thingify(const em_morphology&) = 0;
+ };
+
+ std::unique_ptr<interface> impl_;
+
+ template <typename Impl>
+ struct wrap: interface {
+ explicit wrap(const Impl& impl): wrapped(impl) {}
+ explicit wrap(Impl&& impl): wrapped(std::move(impl)) {}
+
+ virtual std::unique_ptr<interface> clone() override {
+ return std::unique_ptr<interface>(new wrap<Impl>(wrapped));
+ }
+
+ virtual mcable_list thingify(const em_morphology& m) override {
+ return thingify_(wrapped, m);
+ }
+
+ virtual std::ostream& print(std::ostream& o) override {
+ return o << wrapped;
+ }
+
+ Impl wrapped;
+ };
+};
+
+namespace reg {
+
+// An explicit cable section.
+region cable(mcable);
+
+region interval(mlocation, mlocation);
+
+// An explicit branch.
+region branch(msize_t);
+
+// Region with all segments with segment tag id.
+region tagged(int id);
+
+// Region with all segments in a cell.
+region all();
+
+} // namespace reg
+
+} // namespace arb
diff --git a/arbor/morph/em_morphology.cpp b/arbor/morph/em_morphology.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e482c7e062e68eefc4545e0129b7441e1aa7042
--- /dev/null
+++ b/arbor/morph/em_morphology.cpp
@@ -0,0 +1,144 @@
+#include <mutex>
+#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();
+
+ // 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()]);
+ }
+
+ // 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};
+ }
+ // For 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));
+ }
+ }
+}
+
+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;
+}
+
+mlocation em_morphology::canonicalize(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));
+ }
+
+ // 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;
+}
+
+} // namespace arb
diff --git a/arbor/morph/em_morphology.hpp b/arbor/morph/em_morphology.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5b2f3dc4a5fab0ed403bb998f7d07352372907a7
--- /dev/null
+++ b/arbor/morph/em_morphology.hpp
@@ -0,0 +1,39 @@
+#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(const morphology& m);
+
+ const morphology& morph() const;
+
+ mlocation_list terminals() const;
+ mlocation root() const;
+
+ mlocation sample2loc(msize_t sid) 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;
+};
+
+} // namespace arb
diff --git a/arbor/morph/label_dict.cpp b/arbor/morph/label_dict.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6bf10c8a546f83bd0fed9866eb00a0fff0d85e2
--- /dev/null
+++ b/arbor/morph/label_dict.cpp
@@ -0,0 +1,76 @@
+#include <algorithm>
+#include <set>
+#include <fstream>
+#include <unordered_map>
+
+#include <arbor/morph/error.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/label_dict.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 {
+ return locsets_.size() + regions_.size();
+}
+
+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));
+ }
+ // First remove an entry with the same name if it exists.
+ // Has to be this way, because insert_or_assign() is C++17, and we
+ // can't use operator[] because locset is not default constructable.
+ auto it = locsets_.find(name);
+ if (it!=locsets_.end()) {
+ it->second = std::move(ls);
+ }
+ else {
+ locsets_.emplace(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));
+ }
+ // First remove an entry with the same name if it exists.
+ // Has to be this way, because insert_or_assign() is C++17, and we
+ // can't use operator[] because region is not default constructable.
+ auto it = regions_.find(name);
+ if (it!=regions_.end()) {
+ it->second = std::move(reg);
+ }
+ else {
+ regions_.emplace(name, std::move(reg));
+ }
+}
+
+util::optional<const region&> label_dict::region(const std::string& name) const {
+ auto it = regions_.find(name);
+ if (it==regions_.end()) return {};
+ return it->second;
+}
+
+util::optional<const locset&> label_dict::locset(const std::string& name) const {
+ auto it = locsets_.find(name);
+ if (it==locsets_.end()) return {};
+ return it->second;
+}
+
+const std::unordered_map<std::string, locset>& label_dict::locsets() const {
+ return locsets_;
+}
+
+const std::unordered_map<std::string, region>& label_dict::regions() const {
+ return regions_;
+}
+
+} // namespace arb
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f1215c4eaf2fa301988eff86736e373b7d8cda09
--- /dev/null
+++ b/arbor/morph/locset.cpp
@@ -0,0 +1,247 @@
+#include <algorithm>
+#include <iostream>
+#include <numeric>
+
+#include <arbor/morph/error.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/region.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);
+ }
+ 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
+struct nil_ {};
+
+locset nil() {
+ return locset{nil_{}};
+}
+
+mlocation_list thingify_(const nil_& x, const em_morphology& m) {
+ return {};
+}
+
+std::ostream& operator<<(std::ostream& o, const nil_& x) {
+ return o << "nil";
+}
+
+// An explicit location
+struct location_ {
+ mlocation loc;
+};
+
+locset location(mlocation loc) {
+ if (!test_invariants(loc)) {
+ throw morphology_error(util::pprintf("invalid location {}", loc));
+ }
+ return locset{location_{loc}};
+}
+
+mlocation_list thingify_(const location_& x, const em_morphology& m) {
+ // canonicalize will throw if the location is not present.
+ return {m.canonicalize(x.loc)};
+}
+
+std::ostream& operator<<(std::ostream& o, const location_& x) {
+ return o << "(location " << x.loc.branch << " " << x.loc.pos << ")";
+}
+
+
+// Location corresponding to a sample id
+struct sample_ {
+ msize_t index;
+};
+
+locset sample(msize_t index) {
+ return locset{sample_{index}};
+}
+
+mlocation_list thingify_(const sample_& x, const em_morphology& m) {
+ return {m.sample2loc(x.index)};
+}
+
+std::ostream& operator<<(std::ostream& o, const sample_& x) {
+ return o << "(sample " << x.index << ")";
+}
+
+// set of terminal nodes on a morphology
+struct terminal_ {};
+
+locset terminal() {
+ return locset{terminal_{}};
+}
+
+mlocation_list thingify_(const terminal_&, const em_morphology& m) {
+ return m.terminals();
+}
+
+std::ostream& operator<<(std::ostream& o, const terminal_& x) {
+ return o << "terminal";
+}
+
+// the root node of a morphology
+struct root_ {};
+
+locset root() {
+ return locset{root_{}};
+}
+
+mlocation_list thingify_(const root_&, const em_morphology& m) {
+ return {m.root()};
+}
+
+std::ostream& operator<<(std::ostream& o, const root_& x) {
+ return o << "root";
+}
+
+// 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));
+}
+
+std::ostream& operator<<(std::ostream& o, const land& x) {
+ return o << "(intersect " << x.lhs << " " << x.rhs << ")";
+}
+
+// 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));
+}
+
+std::ostream& operator<<(std::ostream& o, const lor& x) {
+ return o << "(join " << x.lhs << " " << x.rhs << ")";
+}
+
+// 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));
+}
+
+std::ostream& operator<<(std::ostream& o, const lsum& x) {
+ return o << "(sum " << x.lhs << " " << x.rhs << ")";
+}
+
+} // namespace ls
+
+// The intersect and join operations in the arb:: namespace with locset so that
+// ADL allows for construction of expressions with locsets without having
+// to namespace qualify the intersect/join.
+
+locset intersect(locset lhs, locset rhs) {
+ return locset(ls::land(std::move(lhs), std::move(rhs)));
+}
+
+locset join(locset lhs, locset rhs) {
+ return locset(ls::lor(std::move(lhs), std::move(rhs)));
+}
+
+locset sum(locset lhs, locset rhs) {
+ return locset(ls::lsum(std::move(lhs), std::move(rhs)));
+}
+
+} // namespace arb
diff --git a/arbor/morph/mbranch.cpp b/arbor/morph/mbranch.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b52b1e1372912300fd3da78e528b1d4548e3dc54
--- /dev/null
+++ b/arbor/morph/mbranch.cpp
@@ -0,0 +1,29 @@
+#include <ostream>
+
+#include <arbor/morph/primitives.hpp>
+
+#include "io/sepval.hpp"
+#include "morph/mbranch.hpp"
+
+namespace arb {
+namespace impl{
+
+//
+// mbranch implementation
+//
+
+bool operator==(const mbranch& l, const mbranch& r) {
+ return l.parent_id==r.parent_id && l.index==r.index;
+}
+
+std::ostream& operator<<(std::ostream& o, const mbranch& b) {
+ o <<"mbranch([" << io::csv(b.index) << "], ";
+ if (b.parent_id==mnpos) o << "none)";
+ else o << b.parent_id << ")";
+ return o;
+}
+
+
+} // namespace impl
+} // namespace arb
+
diff --git a/arbor/morph/mbranch.hpp b/arbor/morph/mbranch.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d680d1d23af1ac9fd7de79f8ba302c868643042
--- /dev/null
+++ b/arbor/morph/mbranch.hpp
@@ -0,0 +1,38 @@
+#pragma once
+
+/*
+ * Definitions and prototypes used in the morphology implementation.
+ * This header is used to simplify unit testing of the morphology implementation.
+ */
+
+#include <cmath>
+#include <vector>
+
+#include <arbor/morph/primitives.hpp>
+
+namespace arb {
+namespace impl{
+
+// An unbranched cable segment that has root, terminal or fork point at each end.
+struct mbranch {
+ std::vector<msize_t> index; // sample index
+ msize_t parent_id = mnpos; // branch index
+
+ mbranch() = default;
+ mbranch(std::vector<msize_t> idx, msize_t parent):
+ index(std::move(idx)), parent_id(parent) {}
+
+ bool is_sphere() const { return size()==1u; }
+ msize_t size() const { return index.size(); }
+ bool has_parent() const { return parent_id!=mnpos;}
+
+ friend bool operator==(const mbranch& l, const mbranch& r);
+ friend std::ostream& operator<<(std::ostream& o, const mbranch& b);
+};
+
+std::vector<mbranch> branches_from_parent_index(const std::vector<msize_t>& parents,
+ const std::vector<point_prop>& props,
+ bool spherical_root);
+
+} // namespace impl
+} // namespace arb
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
index 23d38c1d2cb7448d863828f0ba83ee6124b42c33..35a73def3deffcf56ec3c0b45146cabde9f6c6d4 100644
--- a/arbor/morph/morphology.cpp
+++ b/arbor/morph/morphology.cpp
@@ -11,6 +11,7 @@
#include "algorithms.hpp"
#include "io/sepval.hpp"
+#include "morph/mbranch.hpp"
#include "util/span.hpp"
#include "util/strprintf.hpp"
@@ -93,32 +94,46 @@ bool root_sample_has_same_tag_as_child(const sample_tree& st) {
} // namespace impl
-bool operator==(const mbranch& l, const mbranch& r) {
- return l.parent_id==r.parent_id && l.index==r.index;
-}
+//
+// morphology_impl definition and implementation
+//
-std::ostream& operator<<(std::ostream& o, const mbranch& b) {
- o <<"mbranch([" << io::csv(b.index) << "], ";
- if (b.parent_id==mnpos) o << "none)";
- else o << b.parent_id << ")";
- return o;
-}
+struct morphology_impl {
+ // The sample tree of sample points and their parent-child relationships.
+ sample_tree samples_;
-morphology::morphology(sample_tree m, bool use_spherical_root):
+ // Indicates whether the soma is a sphere.
+ bool spherical_root_ = false;
+
+ // Branch state.
+ std::vector<impl::mbranch> branches_;
+ std::vector<msize_t> branch_parents_;
+ std::vector<msize_t> root_children_;
+ std::vector<std::vector<msize_t>> branch_children_;
+
+ morphology_impl(sample_tree m, bool use_spherical_root);
+ morphology_impl(sample_tree m);
+
+ void init();
+
+ friend std::ostream& operator<<(std::ostream&, const morphology_impl&);
+};
+
+morphology_impl::morphology_impl(sample_tree m, bool use_spherical_root):
samples_(std::move(m)),
spherical_root_(use_spherical_root)
{
init();
}
-morphology::morphology(sample_tree m):
+morphology_impl::morphology_impl(sample_tree m):
samples_(std::move(m)),
spherical_root_(impl::root_sample_has_same_tag_as_child(samples_))
{
init();
}
-void morphology::init() {
+void morphology_impl::init() {
using util::make_span;
using util::count_along;
@@ -139,50 +154,78 @@ void morphology::init() {
if (id!=mnpos) {
branch_children_[id].push_back(i);
}
+ else {
+ root_children_.push_back(i);
+ }
}
}
+std::ostream& operator<<(std::ostream& o, const morphology_impl& m) {
+ o << "morphology: "
+ << m.samples_.size() << " samples, "
+ << m.branches_.size() << " branches.";
+ for (auto i: util::make_span(m.branches_.size()))
+ o << "\n branch " << i << ": " << m.branches_[i];
+
+ return o;
+}
+
+//
+// morphology implementation
+//
+
+morphology::morphology(sample_tree m, bool use_spherical_root):
+ impl_(std::make_shared<const morphology_impl>(std::move(m), use_spherical_root))
+{}
+
+morphology::morphology(sample_tree m):
+ impl_(std::make_shared<const morphology_impl>(std::move(m)))
+{}
+
// The parent branch of branch b.
msize_t morphology::branch_parent(msize_t b) const {
- return branch_parents_[b];
+ return impl_->branch_parents_[b];
}
// The parent sample of sample i.
const std::vector<msize_t>& morphology::sample_parents() const {
- return samples_.parents();
+ return impl_->samples_.parents();
}
// The child branches of branch b.
const std::vector<msize_t>& morphology::branch_children(msize_t b) const {
- return branch_children_[b];
+ return b==mnpos? impl_->root_children_: impl_->branch_children_[b];
}
// Whether the root of the morphology is spherical.
bool morphology::spherical_root() const {
- return spherical_root_;
+ return impl_->spherical_root_;
}
-morphology::index_range morphology::branch_sample_span(msize_t b) const {
- const auto& idx = branches_[b].index;
+mindex_range morphology::branch_indexes(msize_t b) const {
+ const auto& idx = impl_->branches_[b].index;
return std::make_pair(idx.data(), idx.data()+idx.size());
}
const std::vector<msample>& morphology::samples() const {
- return samples_.samples();
+ return impl_->samples_.samples();
+}
+
+// Point properties of samples in the morphology.
+const std::vector<point_prop>& morphology::sample_props() const {
+ return impl_->samples_.properties();
+}
+
+msize_t morphology::num_samples() const {
+ return impl_->samples_.size();
}
msize_t morphology::num_branches() const {
- return branches_.size();
+ return impl_->branches_.size();
}
std::ostream& operator<<(std::ostream& o, const morphology& m) {
- o << "morphology: "
- << m.samples_.size() << " samples, "
- << m.num_branches() << " branches.";
- for (auto i: util::make_span(m.num_branches()))
- o << "\n branch " << i << ": " << m.branches_[i];
-
- return o;
+ return o << *m.impl_;
}
} // namespace arb
diff --git a/arbor/morph/primitives.cpp b/arbor/morph/primitives.cpp
index 6decddbe1c93fa5c2ee35da12d59a9001fff77e8..486959bb60cd2cef279dd54ec5d6806c433a1d51 100644
--- a/arbor/morph/primitives.cpp
+++ b/arbor/morph/primitives.cpp
@@ -28,7 +28,7 @@ double distance(const mpoint& a, const mpoint& b) {
double dy = a.y - b.y;
double dz = a.z - b.z;
- return std::sqrt(dx*dx + dy*dy * dz*dz);
+ return std::sqrt(dx*dx + dy*dy + dz*dz);
}
bool is_collocated(const msample& a, const msample& b) {
@@ -39,12 +39,49 @@ double distance(const msample& a, const msample& b) {
return distance(a.loc, b.loc);
}
+bool test_invariants(const mlocation& l) {
+ return (0.<=l.pos && l.pos<=1.) && l.branch!=mnpos;
+}
+
+bool test_invariants(const mcable& c) {
+ return (0.<=c.prox_pos && c.prox_pos<=c.dist_pos && c.dist_pos<=1.) && c.branch!=mnpos;
+}
+
+mlocation prox_loc(const mcable& c) {
+ return {c.branch, c.prox_pos};
+}
+
+mlocation dist_loc(const mcable& c) {
+ return {c.branch, c.dist_pos};
+}
+
+bool test_invariants(const mcable_list& l) {
+ return std::is_sorted(l.begin(), l.end())
+ && l.end()==std::find_if(l.begin(), l.end(), [](const mcable& c) {return !test_invariants(c);});
+}
+
std::ostream& operator<<(std::ostream& o, const mpoint& p) {
- return o << "mpoint(" << p.x << "," << p.y << "," << p.z << "," << p.radius << ")";
+ return o << "(point " << p.x << " " << p.y << " " << p.z << " " << p.radius << ")";
}
std::ostream& operator<<(std::ostream& o, const msample& s) {
- return o << "msample(" << s.loc << ", " << s.tag << ")";
+ return o << "(sample " << s.loc << " " << s.tag << ")";
+}
+
+std::ostream& operator<<(std::ostream& o, const mlocation& l) {
+ return o << "(location " << l.branch << " " << l.pos << ")";
+}
+
+std::ostream& operator<<(std::ostream& o, const mlocation_list& l) {
+ return o << "(list " << io::sepval(l, ' ') << ")";
+}
+
+std::ostream& operator<<(std::ostream& o, const mcable& c) {
+ return o << "(cable " << c.branch << " " << c.prox_pos << " " << c.dist_pos << ")";
+}
+
+std::ostream& operator<<(std::ostream& o, const mcable_list& c) {
+ return o << "(list " << io::sepval(c, ' ') << ")";
}
} // namespace arb
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..90d999b7fe0a5fcbf4a78aff428987e23c187d61
--- /dev/null
+++ b/arbor/morph/region.cpp
@@ -0,0 +1,303 @@
+#include <set>
+#include <string>
+#include <vector>
+
+#include <arbor/morph/error.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/region.hpp>
+
+#include "morph/em_morphology.hpp"
+#include "util/span.hpp"
+#include "util/strprintf.hpp"
+#include "util/range.hpp"
+#include "util/rangeutil.hpp"
+
+namespace arb {
+namespace reg {
+
+// Returns true iff cable sections a and b:
+// 1. are on the same branch
+// 2. overlap, i.e. their union is not empty.
+bool is_disjoint(const mcable& a, const mcable& b) {
+ if (a.branch!=b.branch) return true;
+ return a<b? a.dist_pos<b.prox_pos: b.dist_pos<a.prox_pos;
+}
+
+// Take the union of two cable sections that are not disjoint.
+mcable make_union(const mcable& a, const mcable& b) {
+ assert(!is_disjoint(a,b));
+ return mcable{a.branch, std::min(a.prox_pos, b.prox_pos), std::max(a.dist_pos, b.dist_pos)};
+}
+
+// Take the intersection of two cable sections that are not disjoint.
+mcable make_intersection(const mcable& a, const mcable& b) {
+ assert(!is_disjoint(a,b));
+
+ return mcable{a.branch, std::max(a.prox_pos, b.prox_pos), std::min(a.dist_pos, b.dist_pos)};
+}
+
+mcable_list merge(const mcable_list& v) {
+ if (v.size()<2) return v;
+ std::vector<mcable> L;
+ L.reserve(v.size());
+ auto c = v.front();
+ auto it = v.begin()+1;
+ while (it!=v.end()) {
+ if (!is_disjoint(c, *it)) {
+ c = make_union(c, *it);
+ }
+ else {
+ L.push_back(c);
+ c = *it;
+ }
+ ++it;
+ }
+ L.push_back(c);
+ 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;
+ 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});
+ }
+ }
+ if (c.dist_pos==1) {
+ for (auto& x: m.cover(mlocation{c.branch, 1}, false)) {
+ L.push_back({x.branch, x.pos, x.pos});
+ }
+ }
+ }
+ L.insert(L.end(), cables.begin(), cables.end());
+ util::sort(L);
+
+ return L;
+}
+
+mcable_list remove_cover(mcable_list cables, const em_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});
+ c = {cloc.branch, cloc.pos, cloc.pos};
+ }
+ }
+ // Some new zero-length cables may be out of order, so sort
+ // the cables.
+ util::sort(cables);
+
+ // Remove multiple copies of zero-length cables if present.
+ return merge(cables);
+}
+
+//
+// Explicit cable section
+//
+
+struct cable_ {
+ mcable cable;
+};
+
+region cable(mcable c) {
+ if (!test_invariants(c)) {
+ throw morphology_error(util::pprintf("Invalid cable section {}", c));
+ }
+ return region(cable_{c});
+}
+
+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));
+ }
+
+ return {reg.cable};
+}
+
+std::ostream& operator<<(std::ostream& o, const cable_& c) {
+ return o << c.cable;
+}
+
+//
+// region with all segments with the same numeric tag
+//
+struct tagged_ {
+ int tag;
+};
+
+region tagged(int id) {
+ return region(tagged_{id});
+}
+
+mcable_list thingify_(const tagged_& reg, const em_morphology& em) {
+ auto& m = em.morph();
+ size_t nb = m.num_branches();
+
+ std::vector<mcable> L;
+ L.reserve(nb);
+ auto& samples = m.samples();
+ auto matches = [reg, &samples](msize_t i) {return samples[i].tag==reg.tag;};
+ auto not_matches = [&matches](msize_t i) {return !matches(i);};
+
+ for (msize_t i: util::make_span(nb)) {
+ auto ids = util::make_range(m.branch_indexes(i)); // range of sample ids in branch.
+ size_t ns = util::size(ids); // number of samples in branch.
+
+ if (ns==1) {
+ // The branch is a spherical soma
+ if (samples[0].tag==reg.tag) {
+ L.push_back({0,0,1});
+ }
+ continue;
+ }
+
+ // The branch has at least 2 samples.
+ // Start at begin+1 because a segment gets its tag from its distal sample.
+ auto beg = std::cbegin(ids);
+ auto end = std::cend(ids);
+
+ // Find the next sample that matches reg.tag.
+ auto start = std::find_if(beg+1, end, matches);
+ while (start!=end) {
+ // find the next sample that does not match reg.tag
+ 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;
+ L.push_back({i, l, r});
+
+ // Find the next sample in the branch that matches reg.tag.
+ start = std::find_if(last, end, matches);
+ }
+ }
+ if (L.size()<L.capacity()/4) {
+ L.shrink_to_fit();
+ }
+ return L;
+}
+
+std::ostream& operator<<(std::ostream& o, const tagged_& t) {
+ return o << "(tag " << t.tag << ")";
+}
+
+//
+// region with all segments in a cell
+//
+struct all_ {};
+
+region all() {
+ return region(all_{});
+}
+
+mcable_list thingify_(const all_&, const em_morphology& m) {
+ auto nb = m.morph().num_branches();
+ mcable_list branches;
+ branches.reserve(nb);
+ for (auto i: util::make_span(nb)) {
+ branches.push_back({i,0,1});
+ }
+ return branches;
+}
+
+std::ostream& operator<<(std::ostream& o, const all_& t) {
+ return o << "all";
+}
+
+//
+// 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) {
+ 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);
+
+ // Perform intersection
+ cable_it_pair it{lhs.begin(), rhs.begin()};
+ cable_it_pair end{lhs.end(), rhs.end()};
+ std::vector<mcable> L;
+
+ bool at_end = it.first==end.first || it.second==end.second;
+ while (!at_end) {
+ bool first_less = *(it.first) < *(it.second);
+ auto& lhs = first_less? it.first: it.second;
+ auto& rhs = first_less? it.second: it.first;
+ if (!is_disjoint(*lhs, *rhs)) {
+ L.push_back(make_intersection(*lhs, *rhs));
+ }
+ if (dist_loc(*lhs) < dist_loc(*rhs)) {
+ ++lhs;
+ }
+ else {
+ ++rhs;
+ }
+ at_end = it.first==end.first || it.second==end.second;
+ }
+
+ return remove_cover(L, m);
+}
+
+std::ostream& operator<<(std::ostream& o, const reg_and& x) {
+ return o << "(intersect " << x.lhs << " " << x.rhs << ")";
+}
+
+//
+// union of two point sets
+//
+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 L;
+ L.resize(lhs.size() + rhs.size());
+
+ std::merge(lhs.begin(), lhs.end(), rhs.begin(), rhs.end(), L.begin());
+
+ return merge(L);
+}
+
+std::ostream& operator<<(std::ostream& o, const reg_or& x) {
+ return o << "(join " << x.lhs << " " << x.rhs << ")";
+}
+
+} // namespace reg
+
+// The intersect and join operations in the arb:: namespace with region so that
+// ADL allows for construction of expressions with regions without having
+// to namespace qualify the intersect/join.
+
+region intersect(region l, region r) {
+ return region{reg::reg_and(std::move(l), std::move(r))};
+}
+
+region join(region l, region r) {
+ return region{reg::reg_or(std::move(l), std::move(r))};
+}
+
+} // namespace arb
diff --git a/arbor/morph/sample_tree.cpp b/arbor/morph/sample_tree.cpp
index 7a687e52ca9eb966470fbdd23f2c9d14ff010c93..d8f82c10b33652407a32f7cde86c6331653d2417 100644
--- a/arbor/morph/sample_tree.cpp
+++ b/arbor/morph/sample_tree.cpp
@@ -70,7 +70,6 @@ msize_t sample_tree::append(msize_t p, const msample& s) {
}
msize_t sample_tree::append(const msample& s) {
- std::cout << "decision : " << empty() << " " << (empty()? mnpos: size()-1) << "\n";
return append(empty()? mnpos: size()-1, s);
}
diff --git a/doc/cpp_cable_cell.rst b/doc/cpp_cable_cell.rst
index 4cf1dce867ae1eeae7a954ce6afde27f96db56a6..56c42a4223ba3d40295a11c200930170260c0553 100644
--- a/doc/cpp_cable_cell.rst
+++ b/doc/cpp_cable_cell.rst
@@ -97,13 +97,11 @@ Density mechanisms are associated with a cable cell object with:
Point mechanisms, which are associated with connection end points on a
cable cell, are attached to a cell with:
-.. cpp:function:: void cable_cell::add_synapse(segment_location loc, mechanism_desc mech)
+.. cpp:function:: void cable_cell::add_synapse(mlocation loc, mechanism_desc mech)
-where :cpp:type:`segment_location` is a simple struct holding a segment index
+where :cpp:type:`mlocation` is a simple struct holding a segment index
and a relative position (from 0, proximal, to 1, distal) along that segment:
-.. cpp:function:: segment_location::segment_location(cell_lid_type index, double position)
-
Electrical properities and ion values
-------------------------------------
diff --git a/example/dryrun/dryrun.cpp b/example/dryrun/dryrun.cpp
index 42b4ca7e4ef77673aaa57105fea3d5848dfc715f..e54bc1a70faeb29bc30a3e001e2dbeb1f18678d3 100644
--- a/example/dryrun/dryrun.cpp
+++ b/example/dryrun/dryrun.cpp
@@ -12,8 +12,9 @@
#include <arbor/assert_macro.hpp>
#include <arbor/common_types.hpp>
#include <arbor/context.hpp>
-#include <arbor/load_balance.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/load_balance.hpp>
+#include <arbor/morph/primitives.hpp>
#include <arbor/profile/meter_manager.hpp>
#include <arbor/profile/profiler.hpp>
#include <arbor/simple_sampler.hpp>
@@ -140,7 +141,7 @@ public:
// Get the appropriate kind for measuring voltage.
cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma.
- arb::segment_location loc(0, 0.0);
+ arb::mlocation loc{0, 0.0};
return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
}
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index f3e1e26ee9c867189abdd823747ddf41d287d33e..8212beec6b8d181de0cf96d089b6dd59c96d007a 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -103,7 +103,7 @@ public:
// Get the appropriate kind for measuring voltage.
cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma.
- arb::segment_location loc(0, 1);
+ arb::mlocation loc{0, 1.};
return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
}
diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp
index b7acd24a06cea07a0bdd2a8e418c2f03f8217b00..9c18c5deba2430698db2ee81c934ce1de2335b62 100644
--- a/example/generators/generators.cpp
+++ b/example/generators/generators.cpp
@@ -123,7 +123,7 @@ public:
// Get the appropriate kind for measuring voltage
cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma
- arb::segment_location loc(0, 0.0);
+ arb::mlocation loc{0, 0.0};
return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
}
diff --git a/example/generators/readme.md b/example/generators/readme.md
index 0f2e71f220d83f1d4ffbfbe6996b230b794d9b6d..d8858740da6fa318d8174c899a39f1fad4d84191 100644
--- a/example/generators/readme.md
+++ b/example/generators/readme.md
@@ -156,9 +156,9 @@ In our case, the cell has one probe, which refers to the voltage at the soma.
// Get the appropriate kind for measuring voltage
cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
- // The location at which we measure: position 0 on segment 0.
- // The cell has only one segment, segment 0, which is the soma.
- arb::segment_location loc(0, 0.0);
+ // The location at which we measure: position 0 on branch 0.
+ // The cell has only one branch, branch 0, which is the soma.
+ arb::mlocation loc(0, 0.0);
// Put this together into a `probe_info`
return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
index 6f4260320fc31fa9ed385c0b9fc436e39a77f28d..4cc95891bbe7a41e0c24933a031c1368b7173376 100644
--- a/example/ring/ring.cpp
+++ b/example/ring/ring.cpp
@@ -134,7 +134,7 @@ public:
// Get the appropriate kind for measuring voltage.
cell_probe_address::probe_kind kind = cell_probe_address::membrane_voltage;
// Measure at the soma.
- arb::segment_location loc(0, 0.0);
+ arb::mlocation loc{0, 0.0};
return arb::probe_info{id, kind, cell_probe_address{loc, kind}};
}
diff --git a/example/single/single.cpp b/example/single/single.cpp
index df614950c3cf6b2060e48593dc71c92a283c8aca..6673ef09460ca8352d774391c5b0d42b0531e58c 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -36,7 +36,7 @@ struct single_recipe: public arb::recipe {
arb::cell_size_type num_targets(arb::cell_gid_type) const override { return 1; }
arb::probe_info get_probe(arb::cell_member_type probe_id) const override {
- arb::segment_location mid_soma = {0, 0.5};
+ arb::mlocation mid_soma = {0, 0.5};
arb::cell_probe_address probe = {mid_soma, arb::cell_probe_address::membrane_voltage};
// Probe info consists of: the probe id, a tag value to distinguish this probe
@@ -54,27 +54,29 @@ struct single_recipe: public arb::recipe {
}
arb::util::unique_any get_cell_description(arb::cell_gid_type) const override {
- arb::cable_cell c = make_cable_cell(morpho, false);
+ arb::label_dict dict;
+ using arb::reg::tagged;
+ dict.set("soma", tagged(1));
+ dict.set("dend", join(tagged(3), tagged(4), tagged(42)));
+ arb::cable_cell c = make_cable_cell(morpho, dict, false);
// Add HH mechanism to soma, passive channels to dendrites.
+ c.paint("soma", "hh");
+ c.paint("dend", "pas");
+
// Discretize dendrites according to the NEURON d-lambda rule.
+ for (std::size_t i=1; i<c.num_segments(); ++i) {
+ arb::cable_segment* branch = c.cable(i);
- c.soma()->add_mechanism("hh");
- auto& segments = c.segments();
- for (std::size_t i = 1; i<segments.size(); ++i) {
double dx = c.segment_length_constant(100., i, gprop.default_parameters)*0.3;
-
- arb::cable_segment* seg = c.cable(i);
- seg->add_mechanism("pas");
-
- unsigned n = std::ceil(seg->length()/dx);
- seg->set_compartments(n);
+ unsigned n = std::ceil(branch->length()/dx);
+ branch->set_compartments(n);
}
// Add synapse to last branch.
- arb::cell_lid_type last_segment = morpho.num_branches()-1;
- arb::segment_location end_last_segment = { last_segment, 1. };
+ arb::cell_lid_type last_segment = c.num_segments()-1;
+ arb::mlocation end_last_segment = { last_segment, 1. };
c.add_synapse(end_last_segment, "exp2syn");
return c;
@@ -151,7 +153,6 @@ arb::morphology default_morphology() {
arb::sample_tree samples;
auto p = samples.append(arb::msample{{ 0.0, 0.0, 0.0, 6.3}, 1});
- std::cout << "p is " << p << "\n";
p = samples.append(p, arb::msample{{ 6.3, 0.0, 0.0, 0.5}, 3});
p = samples.append(p, arb::msample{{206.3, 0.0, 0.0, 0.2}, 3});
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 4bebb1919fb99f20bbb3944a5765071d66ba458a..468264f5c12f9c73b7743bb700d2381859dddc66 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -72,13 +72,15 @@ set(unit_sources
test_algorithms.cpp
test_any.cpp
test_backend.cpp
- test_double_buffer.cpp
- test_dry_run_context.cpp
+ test_cable_cell.cpp
test_compartments.cpp
test_counter.cpp
test_cycle.cpp
test_domain_decomposition.cpp
+ 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
@@ -87,23 +89,23 @@ set(unit_sources
test_fvm_layout.cpp
test_fvm_lowered.cpp
test_glob_basic.cpp
- test_mc_cell_group.cpp
test_kinetic_linear.cpp
test_lexcmp.cpp
test_lif_cell_group.cpp
+ test_local_context.cpp
test_maputil.cpp
test_mask_stream.cpp
test_math.cpp
test_matrix.cpp
- test_cable_cell.cpp
+ test_mc_cell_group.cpp
test_mechanisms.cpp
- test_mech_temperature.cpp
test_mechcat.cpp
+ test_mechinfo.cpp
+ test_mech_temperature.cpp
test_merge_events.cpp
+ test_morphology.cpp
test_multi_event_stream.cpp
test_optional.cpp
- test_mechinfo.cpp
- test_morphology.cpp
test_padded.cpp
test_partition.cpp
test_partition_by_constraint.cpp
@@ -114,7 +116,6 @@ set(unit_sources
test_segment.cpp
test_schedule.cpp
test_spike_source.cpp
- test_local_context.cpp
test_scope_exit.cpp
test_simd.cpp
test_span.cpp
diff --git a/test/unit/test_em_morphology.cpp b/test/unit/test_em_morphology.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c724f5d1f9b772df2893d1c1a7e71a13ceccf094
--- /dev/null
+++ b/test/unit/test_em_morphology.cpp
@@ -0,0 +1,627 @@
+/*
+ * 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}}));
+ }
+
+ // 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}}));
+ }
+ { // 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}}));
+ }
+}
+
+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 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{{0,1}}));
+ EXPECT_EQ(thingify(begb2, em), (ll{{0,1}}));
+ EXPECT_EQ(thingify(begb3, em), (ll{{2,1}}));
+ EXPECT_EQ(thingify(begb4, em), (ll{{2,1}}));
+ }
+ {
+ 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{{0,0}}));
+ EXPECT_EQ(thingify(begb2, em), (ll{{1,1}}));
+ EXPECT_EQ(thingify(begb3, em), (ll{{1,1}}));
+ // 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_morphology.cpp b/test/unit/test_morphology.cpp
index f0af3a4a80e78121496013f73f6009e892f0e90b..fddbd74434a83d13bf1981d9a664170c36b33f03 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -3,9 +3,7 @@
*/
#include <fstream>
-#include <cassert>
#include <cmath>
-#include <stdexcept>
#include <string>
#include <vector>
@@ -19,21 +17,38 @@
#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"
-// Forward declare non-public functions that are used internally to build
-// morphologies so that they can be tested.
-namespace arb { namespace impl {
- std::vector<mbranch> branches_from_parent_index(const std::vector<arb::msize_t>&, const std::vector<point_prop>&, bool);
-}}
-
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) {
+ using mp = arb::mpoint;
+
+ EXPECT_EQ(arb::distance(mp{0,0,0},mp{0 , 0 , 0 }), 0.);
+ EXPECT_EQ(arb::distance(mp{0,0,0},mp{3.14, 0 , 0 }), 3.14);
+ EXPECT_EQ(arb::distance(mp{0,0,0},mp{0 , 3.14, 0 }), 3.14);
+ EXPECT_EQ(arb::distance(mp{0,0,0},mp{0 , 0 , 3.14}), 3.14);
+ EXPECT_EQ(arb::distance(mp{0,0,0},mp{1 , 2 , 3 }), std::sqrt(14.));
+
+ EXPECT_TRUE(arb::is_collocated(mp{0,0,0}, mp{0,0,0}));
+ EXPECT_TRUE(arb::is_collocated(mp{3,0,0}, mp{3,0,0}));
+ EXPECT_TRUE(arb::is_collocated(mp{0,3,0}, mp{0,3,0}));
+ EXPECT_TRUE(arb::is_collocated(mp{0,0,3}, mp{0,0,3}));
+ EXPECT_TRUE(arb::is_collocated(mp{2,0,3}, mp{2,0,3}));
+
+ EXPECT_FALSE(arb::is_collocated(mp{1,0,3}, mp{2,0,3}));
+ EXPECT_FALSE(arb::is_collocated(mp{2,1,3}, mp{2,0,3}));
+ EXPECT_FALSE(arb::is_collocated(mp{2,0,1}, mp{2,0,3}));
+}
+
TEST(morphology, point_props) {
arb::point_prop p = arb::point_prop_mask_none;
@@ -142,7 +157,7 @@ TEST(sample_tree, properties) {
TEST(morphology, branches_from_parent_index) {
const auto npos = arb::mnpos;
using pvec = std::vector<arb::msize_t>;
- using mb = arb::mbranch;
+ using mb = arb::impl::mbranch;
// make a sample tree from a parent vector with non-collocated points.
auto make_tree = [] (const pvec& parents) {
@@ -555,7 +570,7 @@ TEST(morphology, swc) {
EXPECT_EQ(31u, m.num_branches());
// Confirm that converting to a cable_cell generates the same number of branches.
- auto c = arb::make_cable_cell(m, false);
+ auto c = arb::make_cable_cell(m, {}, false);
EXPECT_EQ(31u, c.num_segments());
}
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 989a6cc822cfd3617959e3adc7821c838c168f7b..5afba4f3e89a1f040eaf2e6764a4d58b49192464 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -28,9 +28,9 @@ TEST(probe, fvm_lowered_cell) {
cable1d_recipe rec(bs);
- segment_location loc0{0, 0};
- segment_location loc1{1, 1};
- segment_location loc2{1, 0.3};
+ mlocation loc0{0, 0};
+ mlocation loc1{1, 1};
+ mlocation loc2{1, 0.3};
rec.add_probe(0, 10, cell_probe_address{loc0, cell_probe_address::membrane_voltage});
rec.add_probe(0, 20, cell_probe_address{loc1, cell_probe_address::membrane_voltage});
diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp
index 56700e8352fb9311bb6f8a740ba03ab678cda8c3..b7983f1ebadbb90924da1a461c3c1b52b2a68ee6 100644
--- a/test/unit/test_synapses.cpp
+++ b/test/unit/test_synapses.cpp
@@ -45,16 +45,16 @@ TEST(synapses, add_to_cell) {
EXPECT_EQ(3u, cell.synapses().size());
const auto& syns = cell.synapses();
- EXPECT_EQ(syns[0].location.segment, 0u);
- EXPECT_EQ(syns[0].location.position, 0.1);
+ EXPECT_EQ(syns[0].location.branch, 0u);
+ EXPECT_EQ(syns[0].location.pos, 0.1);
EXPECT_EQ(syns[0].mechanism.name(), "expsyn");
- EXPECT_EQ(syns[1].location.segment, 1u);
- EXPECT_EQ(syns[1].location.position, 0.2);
+ EXPECT_EQ(syns[1].location.branch, 1u);
+ EXPECT_EQ(syns[1].location.pos, 0.2);
EXPECT_EQ(syns[1].mechanism.name(), "exp2syn");
- EXPECT_EQ(syns[2].location.segment, 0u);
- EXPECT_EQ(syns[2].location.position, 0.3);
+ EXPECT_EQ(syns[2].location.branch, 0u);
+ EXPECT_EQ(syns[2].location.pos, 0.3);
EXPECT_EQ(syns[2].mechanism.name(), "expsyn");
}
diff --git a/test/unit/test_threading_exceptions.cpp b/test/unit/test_threading_exceptions.cpp
index c83896406c05b6acf76a07684b7e01ce11a70143..681a0ce4310783e3abd25bbcd6c4b66dfc9fe1ea 100644
--- a/test/unit/test_threading_exceptions.cpp
+++ b/test/unit/test_threading_exceptions.cpp
@@ -3,11 +3,7 @@
#include <csignal>
#include <arbor/domain_decomposition.hpp>
-#include <arbor/load_balance.hpp>
-#include <arbor/recipe.hpp>
-#include <arbor/simulation.hpp>
#include <arbor/context.hpp>
-#include <arbor/cable_cell.hpp>
#include <arbor/version.hpp>
#include "threading/threading.hpp"
diff --git a/test/validation/validate_ball_and_stick.cpp b/test/validation/validate_ball_and_stick.cpp
index a7843adcb7983ae82790557fb30f070be2a64495..4f60de9e945b0b568d5d4fbb1cef29dbe32c563b 100644
--- a/test/validation/validate_ball_and_stick.cpp
+++ b/test/validation/validate_ball_and_stick.cpp
@@ -27,7 +27,7 @@ using namespace arb;
struct probe_point {
const char* label;
- segment_location where;
+ mlocation where;
};
template <typename ProbePointSeq>