diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index a6a497a9d21c0bb9ab0ea718f4138e43e2115c11..55a23170326052d30ee745672826b8bea291393a 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -33,14 +33,13 @@ set(arbor_sources
morph/embed_pwlin.cpp
morph/label_dict.cpp
morph/locset.cpp
- morph/mbranch.cpp
morph/morphexcept.cpp
morph/morphology.cpp
morph/mprovider.cpp
morph/place_pwlin.cpp
morph/primitives.cpp
morph/region.cpp
- morph/sample_tree.cpp
+ morph/segment_tree.cpp
merge_events.cpp
simulation.cpp
partition_load_balance.cpp
diff --git a/arbor/cv_policy.cpp b/arbor/cv_policy.cpp
index fa39de7d481c78f0eac865f6b763498636eb1b34..468d6a7067cee907583d578535b9ca8330053522 100644
--- a/arbor/cv_policy.cpp
+++ b/arbor/cv_policy.cpp
@@ -139,20 +139,13 @@ locset cv_policy_fixed_per_branch::cv_boundary_points(const cable_cell& cell) co
return unique_sum(locset(std::move(points)), ls::cboundary(domain_));
}
-locset cv_policy_every_sample::cv_boundary_points(const cable_cell& cell) const {
+locset cv_policy_every_segment::cv_boundary_points(const cable_cell& cell) const {
const unsigned nbranch = cell.morphology().num_branches();
if (!nbranch) return ls::nil();
- // Always include branch proximal points, so that forks are trivial.
-
- return
- unique_sum(ls::cboundary(domain_),
- ls::restrict(
- util::foldl(
- [](locset l, msize_t sidx) { return sum(std::move(l), ls::sample(sidx)); },
- ls::on_branches(0),
- util::make_span(cell.morphology().num_samples())),
- domain_));
+ return unique_sum(
+ ls::cboundary(domain_),
+ ls::restrict(ls::segment_boundaries(), domain_));
}
} // namespace arb
diff --git a/arbor/include/arbor/cv_policy.hpp b/arbor/include/arbor/cv_policy.hpp
index f59fe418fc6a027daca96dbc3e975a2c4d66c45a..b8869f02e36c1cf22684a5996aa929c800de85e6 100644
--- a/arbor/include/arbor/cv_policy.hpp
+++ b/arbor/include/arbor/cv_policy.hpp
@@ -183,12 +183,12 @@ private:
cv_policy_flag::value flags_;
};
-struct cv_policy_every_sample: cv_policy_base {
- explicit cv_policy_every_sample(region domain = reg::all()):
+struct cv_policy_every_segment: cv_policy_base {
+ explicit cv_policy_every_segment(region domain = reg::all()):
domain_(std::move(domain)) {}
cv_policy_base_ptr clone() const override {
- return cv_policy_base_ptr(new cv_policy_every_sample(*this));
+ return cv_policy_base_ptr(new cv_policy_every_segment(*this));
}
locset cv_boundary_points(const cable_cell&) const override;
@@ -196,7 +196,6 @@ struct cv_policy_every_sample: cv_policy_base {
private:
region domain_;
- cv_policy_flag::value flags_;
};
inline cv_policy default_cv_policy() {
diff --git a/arbor/include/arbor/math.hpp b/arbor/include/arbor/math.hpp
index 08a5ff4e73a200c72c0608eea1c2ccaf99ae4481..cd228184f392b1e510e5ac075719c098c7f028c3 100644
--- a/arbor/include/arbor/math.hpp
+++ b/arbor/include/arbor/math.hpp
@@ -45,18 +45,6 @@ T constexpr volume_frustrum(T L, T r1, T r2) {
return pi<T>/T(3) * (square(r1+r2) - r1*r2) * L;
}
-// Volume of a sphere radius r.
-template <typename T>
-T constexpr volume_sphere(T r) {
- return T(4)/T(3) * pi<T> * cube(r);
-}
-
-// Surface area of a sphere radius r.
-template <typename T>
-T constexpr area_sphere(T r) {
- return T(4) * pi<T> * square(r);
-}
-
// Linear interpolation by u in interval [a,b]: (1-u)*a + u*b.
template <typename T, typename U>
T constexpr lerp(T a, T b, U u) {
diff --git a/arbor/include/arbor/morph/embed_pwlin.hpp b/arbor/include/arbor/morph/embed_pwlin.hpp
index 0a352a5708029cfe2d5ee7e40ef81f3caa8d3aa7..ac065005e1a3a7f7390d7284c6c8a23fc1870472 100644
--- a/arbor/include/arbor/morph/embed_pwlin.hpp
+++ b/arbor/include/arbor/morph/embed_pwlin.hpp
@@ -22,8 +22,14 @@ 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);
+ // Locations that mark the boundaries between segments.
+ const std::vector<mlocation>& segment_locations() const {
+ return segment_locations_;
+ }
+
+ std::pair<const mlocation*, const mlocation*> branch_segment_locations(msize_t i) const {
+ const mlocation* p = segment_locations_.data();
+ return std::make_pair(p+branch_segment_part_[i], p+branch_segment_part_[i+1]);
}
// Interpolated radius at location.
@@ -58,9 +64,9 @@ struct embed_pwlin {
return integrate_length(mcable{bid, 0, 1});
}
-
private:
- std::vector<mlocation> sample_locations_;
+ std::vector<mlocation> segment_locations_;
+ std::vector<msize_t> branch_segment_part_;
std::shared_ptr<embed_pwlin_data> data_;
};
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
index d6e3772b5970a83443530d70f2ba4c943eedd326..5d7d35e09512b0b13589f83f111e827ddc22dae1 100644
--- a/arbor/include/arbor/morph/locset.hpp
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -121,9 +121,6 @@ namespace ls {
// Explicit location on morphology.
locset location(msize_t branch, double pos);
-// Location of a sample.
-locset sample(msize_t);
-
// Set of terminal nodes on a morphology.
locset terminal();
@@ -152,6 +149,9 @@ locset cboundary(region reg);
// Returns all locations in a locset that are also in the region.
locset restrict(locset ls, region reg);
+// Returns locations that mark the segments.
+locset segment_boundaries();
+
// A range `left` to `right` of randomly selected locations with a
// uniform distribution from region `reg` generated using `seed`
locset uniform(region reg, unsigned left, unsigned right, uint64_t seed);
diff --git a/arbor/include/arbor/morph/morphexcept.hpp b/arbor/include/arbor/morph/morphexcept.hpp
index 9620f6dc02fb1a71973ec4173037a7a5cc0a32d4..bd7872c8aa6ad5febd66e00ddedc471c6d0a5934 100644
--- a/arbor/include/arbor/morph/morphexcept.hpp
+++ b/arbor/include/arbor/morph/morphexcept.hpp
@@ -21,6 +21,11 @@ struct no_such_branch: morphology_error {
msize_t bid;
};
+struct no_such_segment: arbor_exception {
+ explicit no_such_segment(msize_t sid);
+ msize_t sid;
+};
+
struct invalid_mcable: morphology_error {
invalid_mcable(mcable cable);
mcable cable;
@@ -30,8 +35,8 @@ struct invalid_mcable_list: morphology_error {
invalid_mcable_list();
};
-struct invalid_sample_parent: morphology_error {
- invalid_sample_parent(msize_t parent, msize_t tree_size);
+struct invalid_segment_parent: morphology_error {
+ invalid_segment_parent(msize_t parent, msize_t tree_size);
msize_t parent;
msize_t tree_size;
};
diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
index 208c9926f6941ddd7d5cc3e1355921b2132b3d34..a4a6b622a48577d2138f9a02d65766a10603645d 100644
--- a/arbor/include/arbor/morph/morphology.hpp
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -6,35 +6,26 @@
#include <arbor/util/lexcmp_def.hpp>
#include <arbor/morph/primitives.hpp>
-#include <arbor/morph/sample_tree.hpp>
+#include <arbor/morph/segment_tree.hpp>
namespace arb {
struct morphology_impl;
-using mindex_range = std::pair<const msize_t*, const msize_t*>;
-
class morphology {
// Hold an immutable copy of the morphology implementation.
std::shared_ptr<const morphology_impl> impl_;
public:
- morphology(sample_tree m, bool use_spherical_root);
- morphology(sample_tree m);
+ morphology(segment_tree m);
morphology();
// Empty/default-constructed morphology?
bool empty() const;
- // Whether the root of the morphology is spherical.
- bool spherical_root() const;
-
// The number of branches in the morphology.
msize_t num_branches() const;
- // The number of samples in the morphology.
- msize_t num_samples() const;
-
// The parent branch of branch b.
// Return mnpos if branch has no parent.
msize_t branch_parent(msize_t b) const;
@@ -46,17 +37,8 @@ public:
// 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;
-
- // 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;
+ // Range of segments in a branch.
+ const std::vector<msegment>& branch_segments(msize_t b) 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 f325cc20e6baba6641bbd3f8415f79012ada1fbb..cf75c0d4c5078397f1b39bfeff066f5bd181ea18 100644
--- a/arbor/include/arbor/morph/primitives.hpp
+++ b/arbor/include/arbor/morph/primitives.hpp
@@ -22,6 +22,7 @@ struct mpoint {
double x, y, z; // [µm]
double radius; // [μm]
+ friend bool operator==(const mpoint& l, const mpoint& r);
friend std::ostream& operator<<(std::ostream&, const mpoint&);
friend bool operator==(const mpoint& a, const mpoint& b) {
return a.x==b.x && a.y==b.y && a.z==b.z && a.radius==b.radius;
@@ -43,25 +44,16 @@ enum class comp_op {
ge
};
-// A morphology sample consists of a location and an integer tag.
-// When loaded from an SWC file, the tag will correspond to the SWC label,
-// which are standardised as follows:
-// 1 - soma
-// 2 - axon
-// 3 - (basal) dendrite
-// 4 - apical dendrite
-// http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
-// However, any positive integer tag can be provided and labelled dynamically.
-
-struct msample {
- mpoint loc;
+// Describe a cable segment between two adjacent samples.
+struct msegment {
+ mpoint prox;
+ mpoint dist;
int tag;
- friend std::ostream& operator<<(std::ostream&, const msample&);
+ friend bool operator==(const msegment&, const msegment&);
+ friend std::ostream& operator<<(std::ostream&, const msegment&);
};
-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 {
@@ -125,31 +117,6 @@ std::ostream& operator<<(std::ostream& o, const mcable_list& c);
// 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,
- point_prop_mask_root = 1,
- point_prop_mask_fork = 2,
- point_prop_mask_terminal = 4,
- point_prop_mask_collocated = 8
-};
-
-#define ARB_PROP(prop) \
-constexpr bool is_##prop(point_prop p) {\
- return p&point_prop_mask_##prop;\
-} \
-inline void set_##prop(point_prop& p) {\
- p |= point_prop_mask_##prop;\
-} \
-inline void unset_##prop(point_prop& p) {\
- p &= ~point_prop_mask_##prop;\
-}
-
-ARB_PROP(root)
-ARB_PROP(fork)
-ARB_PROP(terminal)
-ARB_PROP(collocated)
-
} // namespace arb
ARB_DEFINE_HASH(arb::mcable, a.branch, a.prox_pos, a.dist_pos);
diff --git a/arbor/include/arbor/morph/sample_tree.hpp b/arbor/include/arbor/morph/sample_tree.hpp
deleted file mode 100644
index 2ae690ebdc13cc9e95e440c763515de9e696c3a1..0000000000000000000000000000000000000000
--- a/arbor/include/arbor/morph/sample_tree.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-#pragma once
-
-#include <cassert>
-#include <functional>
-#include <vector>
-#include <string>
-
-#include <arbor/swcio.hpp>
-#include <arbor/morph/primitives.hpp>
-
-namespace arb {
-
-/// Morphology composed of samples.
-class sample_tree {
- std::vector<msample> samples_;
- std::vector<msize_t> parents_;
- std::vector<point_prop> props_;
-
-public:
- sample_tree() = default;
- sample_tree(std::vector<msample>, std::vector<msize_t>);
-
- // Reserve space for n samples.
- void reserve(msize_t n);
-
- // The append functions return a handle to the last sample appended by the call.
-
- // Append a single sample.
- msize_t append(msize_t p, const msample& s); // to sample p.
- msize_t append(const msample& s); // to the last sample in the tree.
-
- // Append a sequence of samples.
- msize_t append(msize_t p, const std::vector<msample>& slist); // to sample p.
- msize_t append(const std::vector<msample>& slist); // to the last sample in the tree.
-
- // The number of samples in the tree.
- msize_t size() const;
- bool empty() const;
-
- // The samples in the tree.
- const std::vector<msample>& samples() const;
-
- // The parent index of the samples.
- const std::vector<msize_t>& parents() const;
-
- // The properties of the samples.
- const std::vector<point_prop>& properties() const;
-
- friend std::ostream& operator<<(std::ostream&, const sample_tree&);
-};
-
-/// Build a sample tree from a sequence of swc records.
-sample_tree swc_as_sample_tree(const std::vector<swc_record>& swc_records);
-
-} // namesapce arb
-
diff --git a/arbor/include/arbor/morph/segment_tree.hpp b/arbor/include/arbor/morph/segment_tree.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f3686965d9f992abad2aa63b20091c45f67e0eab
--- /dev/null
+++ b/arbor/include/arbor/morph/segment_tree.hpp
@@ -0,0 +1,62 @@
+#pragma once
+
+#include <cassert>
+#include <functional>
+#include <vector>
+#include <string>
+
+#include <arbor/swcio.hpp>
+#include <arbor/morph/primitives.hpp>
+
+namespace arb {
+
+/// Morphology composed of segments.
+class segment_tree {
+ struct child_prop {
+ int count;
+ bool is_fork() const { return count>1; }
+ bool is_terminal() const { return count==0; }
+ int increment() { return ++count;}
+ };
+
+ std::vector<msegment> segments_;
+ std::vector<msize_t> parents_;
+ std::vector<child_prop> seg_children_;
+
+public:
+ segment_tree() = default;
+
+ // Reserve space for n segments.
+ void reserve(msize_t n);
+
+ // The append functions return a handle to the last segment appended by the call.
+
+ // Append a single segment.
+ msize_t append(msize_t p, const mpoint& prox, const mpoint& dist, int tag);
+ msize_t append(msize_t p, const mpoint& dist, int tag);
+
+ // The number of segments in the tree.
+ msize_t size() const;
+ bool empty() const;
+
+ // The segments in the tree.
+ const std::vector<msegment>& segments() const;
+
+ // The parent index of the segments.
+ const std::vector<msize_t>& parents() const;
+
+ // Interfaces for querying the properties of the segments by index.
+
+ bool is_fork(msize_t i) const;
+ bool is_terminal(msize_t i) const;
+ bool is_root(msize_t i) const;
+
+ friend std::ostream& operator<<(std::ostream&, const segment_tree&);
+};
+
+/// Build a sample tree from a sequence of swc records.
+segment_tree swc_as_segment_tree(const std::vector<swc_record>& swc_records);
+
+} // namesapce arb
+
+
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index 64971a03e62d5f7b6780ca4cc7e4eeade5b3aae3..2e381804a8ae99f010e684d52c057ce450b183ac 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/arbor/morph/embed_pwlin.cpp
@@ -201,108 +201,91 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
constexpr double pi = math::pi<double>;
msize_t n_branch = m.num_branches();
data_ = std::make_shared<embed_pwlin_data>(n_branch);
+ branch_segment_part_.reserve(n_branch+1);
if (!n_branch) return;
- const auto& samples = m.samples();
- sample_locations_.resize(m.num_samples());
+ branch_segment_part_.push_back(0);
- double proj_shift = samples.front().loc.z;
+ double proj_shift = m.branch_segments(0).front().prox.z;
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);
+ auto& segments = m.branch_segments(bid);
+ arb_assert(segments.size());
- // Treat spherical root as area-equivalent cylinder.
- double r = samples[0].loc.radius;
+ branch_segment_part_.push_back(branch_segment_part_.back()+segments.size()+1);
- data_->directed_projection[bid].push_back(0., 1., rat_element<1, 0>(-r, r));
- 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));
+ std::vector<double> seg_pos;
+ seg_pos.reserve(segments.size()+1);
+ seg_pos.push_back(0.);
- 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};
+ for (auto &seg: segments) {
+ double d = distance(seg.prox, seg.dist);
+ seg_pos.push_back(seg_pos.back()+d);
}
- else {
- arb_assert(sample_indices.size()>1);
+ double branch_length = seg_pos.back();
+ double length_scale = branch_length>0? 1./branch_length: 0;
+ for (auto& d: seg_pos) {
+ d *= length_scale;
+ }
+ seg_pos.back() = 1; // Circumvent any rounding infelicities.
- std::vector<double> sample_distance;
- sample_distance.reserve(samples.size());
- sample_distance.push_back(0.);
+ for (auto d: seg_pos) {
+ segment_locations_.push_back({bid, d});
+ }
- for (auto i: util::count_along(sample_indices)) {
- if (!i) continue;
+ 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 d = distance(samples[sample_indices[i-1]], samples[sample_indices[i]]);
- sample_distance.push_back(sample_distance.back()+d);
- }
+ double area_0 = parent==mnpos? 0: data_->area[parent].back().second[2];
+ double ixa_0 = parent==mnpos? 0: data_->ixa[parent].back().second[2];
- 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[2];
- double ixa_0 = parent==mnpos? 0: data_->ixa[parent].back().second[2];
-
- if (length_scale==0) {
- // Zero-length branch? Weird, but make best show of it.
- double r = samples[sample_indices[0]].loc.radius;
- double z = samples[sample_indices[0]].loc.z;
- data_->radius[bid].push_back(0., 1., rat_element<1, 0>(r, r));
- data_->directed_projection[bid].push_back(0., 1., rat_element<1, 0>(z-proj_shift, z-proj_shift));
- 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 z0 = samples[sample_indices[i-1]].loc.z - proj_shift;
- double z1 = samples[sample_indices[i]].loc.z - proj_shift;
- data_->directed_projection[bid].push_back(p0, p1, rat_element<1, 0>(z0, z1));
-
- 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;
- }
+ if (length_scale==0) {
+ // Zero-length branch? Weird, but make best show of it.
+ auto& s = segments.front().prox;
+ double r = s.radius;
+ double z = s.z;
+ data_->radius[bid].push_back(0., 1., rat_element<1, 0>(r, r));
+ data_->directed_projection[bid].push_back(0., 1., rat_element<1, 0>(z-proj_shift, z-proj_shift));
+ 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(segments)) {
+ auto prox = segments[i].prox;
+ auto dist = segments[i].dist;
+
+ double p0 = seg_pos[i];
+ double p1 = seg_pos[i+1];
+ if (p0==p1) continue;
+
+ double z0 = prox.z - proj_shift;
+ double z1 = dist.z - proj_shift;
+ data_->directed_projection[bid].push_back(p0, p1, rat_element<1, 0>(z0, z1));
+
+ double r0 = prox.radius;
+ double r1 = dist.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());
}
}
};
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index 393f9b025a35939fcd6ae1430b641e806eb5acfb..dca7c313a770e23c39b0e86dd6ae56ee57892f84 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -96,25 +96,6 @@ std::ostream& operator<<(std::ostream& o, const location_list_& x) {
return o << ')';
}
-// Location corresponding to a sample id.
-
-struct sample_: locset_tag {
- explicit sample_(msize_t index): index(index) {}
- msize_t index;
-};
-
-locset sample(msize_t index) {
- return locset{sample_{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 points (most distal points).
struct terminal_: locset_tag {};
@@ -151,6 +132,23 @@ std::ostream& operator<<(std::ostream& o, const root_& x) {
return o << "(root)";
}
+// Locations that mark interface between segments.
+
+struct segments_: locset_tag {};
+
+locset segment_boundaries() {
+ return locset{segments_{}};
+}
+
+mlocation_list thingify_(const segments_&, const mprovider& p) {
+ return p.embedding().segment_locations();
+}
+
+std::ostream& operator<<(std::ostream& o, const segments_& x) {
+ return o << "(segment_boundaries)";
+}
+
+
// Proportional location on every branch.
struct on_branches_ { double pos; };
diff --git a/arbor/morph/mbranch.cpp b/arbor/morph/mbranch.cpp
deleted file mode 100644
index b52b1e1372912300fd3da78e528b1d4548e3dc54..0000000000000000000000000000000000000000
--- a/arbor/morph/mbranch.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#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
deleted file mode 100644
index 5d680d1d23af1ac9fd7de79f8ba302c868643042..0000000000000000000000000000000000000000
--- a/arbor/morph/mbranch.hpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#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/morphexcept.cpp b/arbor/morph/morphexcept.cpp
index 3364127573b08269f343f197d3dc08fad05c8f0c..ac9222bbc8898afaec79c1d04683808df294e500 100644
--- a/arbor/morph/morphexcept.cpp
+++ b/arbor/morph/morphexcept.cpp
@@ -10,16 +10,26 @@ namespace arb {
using arb::util::pprintf;
+static std::string msize_string(msize_t x) {
+ return x==mnpos? "mnpos": pprintf("{}", x);
+}
+
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)),
+ morphology_error(pprintf("no such branch id {}", msize_string(bid))),
bid(bid)
{}
+no_such_segment::no_such_segment(msize_t id):
+ arbor_exception(pprintf("segment {} out of bounds", id)),
+ sid(id)
+{}
+
+
invalid_mcable::invalid_mcable(mcable cable):
morphology_error(pprintf("invalid mcable {}", cable)),
cable(cable)
@@ -29,19 +39,13 @@ invalid_mcable_list::invalid_mcable_list():
morphology_error("bad mcable_list")
{}
-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)),
+ morphology_error(pprintf("insufficent samples to define branch id {}", msize_string(bid))),
bid(bid)
{}
@@ -55,5 +59,11 @@ circular_definition::circular_definition(const std::string& name):
name(name)
{}
+invalid_segment_parent::invalid_segment_parent(msize_t parent, msize_t tree_size):
+ morphology_error(pprintf("invalid segment parent {} for a segment tree of size {}", msize_string(parent), tree_size)),
+ parent(parent),
+ tree_size(tree_size)
+{}
+
} // namespace arb
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
index 623eb7d9bc0d51c79fde9724796606f075c08dfc..ee861e749461f1bb49961e5216a62d1ee5410ae9 100644
--- a/arbor/morph/morphology.cpp
+++ b/arbor/morph/morphology.cpp
@@ -5,84 +5,67 @@
#include <arbor/morph/morphexcept.hpp>
#include <arbor/morph/morphology.hpp>
-#include <arbor/morph/sample_tree.hpp>
+#include <arbor/morph/segment_tree.hpp>
#include <arbor/morph/primitives.hpp>
-#include "morph/mbranch.hpp"
+#include "io/sepval.hpp"
#include "util/mergeview.hpp"
#include "util/rangeutil.hpp"
#include "util/span.hpp"
+#include "util/transform.hpp"
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)
-{
- auto nsamp = parents.size();
- if (!nsamp) return {};
-
- // Enforce that a morphology with one sample has a spherical root.
- if (!spherical_root && nsamp==1u) {
- throw incomplete_branch(0);
- }
-
- std::vector<int> bids(nsamp);
- int nbranches = spherical_root? 1: 0;
- for (auto i: make_span(1, nsamp)) {
- auto p = parents[i];
- bool first = is_root(props[p]) || is_fork(props[p]);
- bids[i] = first? nbranches++: bids[p];
- }
-
- std::vector<mbranch> branches(nbranches);
- for (auto i: make_span(nsamp)) {
- auto p = parents[i];
- auto& branch = branches[bids[i]];
- if (!branch.size()) {
- bool null_root = spherical_root? p==mnpos: p==mnpos||p==0;
- branch.parent_id = null_root? mnpos: bids[p];
-
- // Add the distal sample from the parent branch if the parent is not
- // a spherical root or mnpos.
- if (p!=mnpos && !(p==0 && spherical_root)) {
- branch.index.push_back(p);
- }
+auto branches_from_segment_tree(const segment_tree& tree) {
+ using rtype = std::pair<std::vector<msize_t>,
+ std::vector<std::vector<msegment>>>;
+ rtype branches;
+ auto& branch_parents = branches.first;
+ auto& branch_segs = branches.second;
+ auto& seg_parents = tree.parents();
+ auto& segs = tree.segments();
+
+ auto nsegs = seg_parents.size();
+ if (!nsegs) return branches;
+
+ // Determine which branch each segment belongs to while counting the number
+ // of branches in the morphology.
+ std::vector<msize_t> bids(nsegs);
+ int nbranches = 1;
+ bids[0] = 0;
+ for (auto i: make_span(1, nsegs)) {
+ if (tree.is_root(i)) {
+ bids[i] = nbranches++;
+ }
+ else {
+ auto p = seg_parents[i];
+ bids[i] = tree.is_fork(p)? nbranches++: bids[p];
}
-
- branch.index.push_back(i);
}
- // Enforce that all cable branches that are potentially connected to a spherical
- // root contain at least two samples.
- if (spherical_root) {
- for (auto i: make_span(1, nbranches)) { // skip the root.
- if (branches[i].size()<2u) {
- throw incomplete_branch(i);
- }
+ // A working vector used to track whether the first segment in a branch has been visited.
+ std::vector<char> visited(nbranches);
+
+ branch_segs.resize(nbranches);
+ branch_parents.resize(nbranches);
+ // Construct all of the cable segments for all of the branches.
+ for (auto i: make_span(nsegs)) {
+ auto p = seg_parents[i];
+ auto b = bids[i];
+ // If this is the first sample in the branch, set the branch's parent branch.
+ if (!visited[b]) {
+ branch_parents[b] = p==mnpos? mnpos: bids[p];
+ visited[b] = 1;
}
+ branch_segs[b].push_back(segs[i]);
}
return branches;
}
-// Returns false if one of the root's children has the same tag as the root.
-bool root_sample_tag_differs_from_children(const sample_tree& st) {
- if (st.empty()) return false;
- auto& P = st.parents();
- auto& S = st.samples();
- auto root_tag = S.front().tag;
- for (auto i: util::make_span(1, st.size())) {
- if (0u==P[i] && S[i].tag==root_tag) {
- return false;
- }
- }
- return true;
-}
-
} // namespace impl
//
@@ -90,55 +73,35 @@ bool root_sample_tag_differs_from_children(const sample_tree& st) {
//
struct morphology_impl {
- // 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<impl::mbranch> branches_;
+ std::vector<std::vector<msegment>> 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);
- morphology_impl(sample_tree m);
+ morphology_impl(const segment_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_impl::morphology_impl(sample_tree m):
- samples_(std::move(m)),
- spherical_root_(impl::root_sample_tag_differs_from_children(samples_))
-{
- init();
-}
-
-void morphology_impl::init() {
- auto nsamp = samples_.size();
+morphology_impl::morphology_impl(const segment_tree& tree) {
+ auto nsamp = tree.size();
if (!nsamp) return;
// Generate branches.
- branches_ = impl::branches_from_parent_index(samples_.parents(), samples_.properties(), spherical_root_);
+ auto B = impl::branches_from_segment_tree(tree);
+ branches_ = std::move(B.second);
+ branch_parents_ = std::move(B.first);
auto nbranch = branches_.size();
// Generate branch tree.
branch_children_.resize(nbranch);
branch_parents_.reserve(nbranch);
for (auto i: make_span(nbranch)) {
- auto id = branches_[i].parent_id;
- branch_parents_.push_back(id);
+ auto id = branch_parents_[i];
if (id!=mnpos) {
branch_children_[id].push_back(i);
}
@@ -158,29 +121,29 @@ void morphology_impl::init() {
}
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;
+ if (m.branches_.empty()) {
+ return o << "(morphology ())";
+ }
+ bool first = true;
+ o << "(morphology\n (";
+ for (auto i: util::make_span(m.branches_.size())) {
+ if (!first) o << "\n ";
+ o << "(" << m.branch_parents_[i] << " (" << io::sepval(m.branches_[i], " ") << "))";
+ first = false;
+ }
+ 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):
+morphology::morphology(segment_tree m):
impl_(std::make_shared<const morphology_impl>(std::move(m)))
{}
morphology::morphology():
- morphology(sample_tree())
+ morphology(segment_tree())
{}
bool morphology::empty() const {
@@ -192,11 +155,6 @@ msize_t morphology::branch_parent(msize_t b) const {
return impl_->branch_parents_[b];
}
-// The parent sample of sample i.
-const std::vector<msize_t>& morphology::sample_parents() const {
- return impl_->samples_.parents();
-}
-
// The child branches of branch b.
const std::vector<msize_t>& morphology::branch_children(msize_t b) const {
return b==mnpos? impl_->root_children_: impl_->branch_children_[b];
@@ -206,27 +164,8 @@ 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_;
-}
-
-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 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();
+const std::vector<msegment>& morphology::branch_segments(msize_t b) const {
+ return impl_->branches_[b];
}
msize_t morphology::num_branches() const {
diff --git a/arbor/morph/place_pwlin.cpp b/arbor/morph/place_pwlin.cpp
index 4f7a78ef55ea7f2daf752afd26f6142c5f4475ec..9fe9afdb1643e88886fc1e8820fe27cec00dcb4e 100644
--- a/arbor/morph/place_pwlin.cpp
+++ b/arbor/morph/place_pwlin.cpp
@@ -8,6 +8,7 @@
#include "morph/pwlin_common.hpp"
#include "util/piecewise.hpp"
+#include "util/rangeutil.hpp"
#include "util/ratelem.hpp"
#include "util/span.hpp"
@@ -36,106 +37,48 @@ place_pwlin::place_pwlin(const arb::morphology& m, const isometry& iso) {
if (!n_branch) return;
- const auto& samples = m.samples();
-
std::vector<double> sample_pos_on_branch;
+ std::vector<double> seg_pos;
for (msize_t bid = 0; bid<n_branch; ++bid) {
- auto sample_indices = util::make_range(m.branch_indexes(bid));
- if (bid==0 && m.spherical_root()) {
- arb_assert(sample_indices.size()==1);
- arb_assert(sample_indices.front()==0);
-
- // Use the next distinct sample, if it exists, to determine the
- // proximal-distal axis for a spherical root.
-
- mpoint c = samples[0].loc;
- mpoint d = c;
- for (msize_t i = 1; i<samples.size(); ++i) {
- const auto p = samples[i].loc;
- if (p.x!=c.x || p.y!=c.y || p.z!=c.z) {
- d = p;
- break;
- }
- }
-
- mpoint u0 = c, u1 = c;
+ auto& segments = m.branch_segments(bid);
+ arb_assert(!segments.empty());
- if (c.x!=d.x || c.y!=d.y || c.z!=d.z) {
- double dx = d.x-c.x;
- double dy = d.y-c.y;
- double dz = d.z-c.z;
- double scale = c.radius/std::sqrt(dx*dx+dy*dy+dz*dz);
-
- u0.x -= dx*scale;
- u0.y -= dy*scale;
- u0.z -= dz*scale;
+ seg_pos.reserve(segments.size()+1);
+ seg_pos = {0};
+ for (auto& seg: segments) {
+ seg_pos.push_back(seg_pos.back()+distance(seg.prox, seg.dist));
+ }
- u1.x += dx*scale;
- u1.y += dy*scale;
- u1.z += dz*scale;
- }
- else {
- // No luck, so pick distal as negative z-axis.
- u0.z += c.radius;
- u1.z -= c.radius;
- }
+ double branch_length = seg_pos.back();
+ double length_scale = branch_length>0? 1./branch_length: 0;
+ for (auto& x: seg_pos) {
+ x *= length_scale;
+ }
- u0 = iso.apply(u0);
- u1 = iso.apply(u1);
+ if (length_scale==0) {
+ // Zero-length branch case?
+ mpoint p = iso.apply(segments[0].prox);
- data_->x[bid].push_back(0., 1., rat_element<1, 0>(u0.x, u1.x));
- data_->y[bid].push_back(0., 1., rat_element<1, 0>(u0.y, u1.y));
- data_->z[bid].push_back(0., 1., rat_element<1, 0>(u0.z, u1.z));
- data_->r[bid].push_back(0., 1., rat_element<1, 0>(u0.radius, u1.radius));
+ data_->x[bid].push_back(0., 1., rat_element<1, 0>(p.x, p.x));
+ data_->y[bid].push_back(0., 1., rat_element<1, 0>(p.y, p.y));
+ data_->z[bid].push_back(0., 1., rat_element<1, 0>(p.z, p.z));
+ data_->r[bid].push_back(0., 1., rat_element<1, 0>(p.radius, p.radius));
}
else {
- arb_assert(sample_indices.size()>1);
-
- sample_pos_on_branch.reserve(samples.size());
- sample_pos_on_branch = {0};
-
- for (auto i: util::count_along(sample_indices)) {
- if (!i) continue;
-
- sample_pos_on_branch.push_back(
- sample_pos_on_branch.back()+
- distance(samples[sample_indices[i-1]], samples[sample_indices[i]]));
- }
-
- double branch_length = sample_pos_on_branch.back();
- double length_scale = branch_length>0? 1./branch_length: 0;
-
- for (auto& x: sample_pos_on_branch) {
- x *= length_scale;
- }
-
- if (length_scale==0) {
- // Zero-length branch case?
-
- mpoint p = iso.apply(samples[sample_indices[0]].loc);
-
- data_->x[bid].push_back(0., 1., rat_element<1, 0>(p.x, p.x));
- data_->y[bid].push_back(0., 1., rat_element<1, 0>(p.y, p.y));
- data_->z[bid].push_back(0., 1., rat_element<1, 0>(p.z, p.z));
- data_->r[bid].push_back(0., 1., rat_element<1, 0>(p.radius, p.radius));
- }
- else {
- for (auto i: util::count_along(sample_indices)) {
- if (!i) continue;
-
- double p0 = i>1? sample_pos_on_branch[i-1]: 0;
- double p1 = sample_pos_on_branch[i];
- if (p0==p1) continue;
-
- mpoint u0 = iso.apply(samples[sample_indices[i-1]].loc);
- mpoint u1 = iso.apply(samples[sample_indices[i]].loc);
-
- data_->x[bid].push_back(p0, p1, rat_element<1, 0>(u0.x, u1.x));
- data_->y[bid].push_back(p0, p1, rat_element<1, 0>(u0.y, u1.y));
- data_->z[bid].push_back(p0, p1, rat_element<1, 0>(u0.z, u1.z));
- data_->r[bid].push_back(p0, p1, rat_element<1, 0>(u0.radius, u1.radius));
- }
+ for (auto i: util::count_along(segments)) {
+ auto& seg = segments[i];
+ double p0 = seg_pos[i];
+ double p1 = seg_pos[i+1];
+ if (p0==p1) continue;
+
+ mpoint u0 = iso.apply(seg.prox);
+ mpoint u1 = iso.apply(seg.dist);
+
+ data_->x[bid].push_back(p0, p1, rat_element<1, 0>(u0.x, u1.x));
+ data_->y[bid].push_back(p0, p1, rat_element<1, 0>(u0.y, u1.y));
+ data_->z[bid].push_back(p0, p1, rat_element<1, 0>(u0.z, u1.z));
+ data_->r[bid].push_back(p0, p1, rat_element<1, 0>(u0.radius, u1.radius));
}
}
}
diff --git a/arbor/morph/primitives.cpp b/arbor/morph/primitives.cpp
index 40d9061b4ad430697f8236017c611095a9d265a7..ced745abf2010645ac9b9b83a3b50a710c61afbc 100644
--- a/arbor/morph/primitives.cpp
+++ b/arbor/morph/primitives.cpp
@@ -57,14 +57,6 @@ double distance(const mpoint& a, const mpoint& b) {
return std::sqrt(dx*dx + dy*dy + dz*dz);
}
-bool is_collocated(const msample& a, const msample& b) {
- return is_collocated(a.loc, b.loc);
-}
-
-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;
}
@@ -140,12 +132,16 @@ bool test_invariants(const mcable_list& l) {
&& l.end()==std::find_if(l.begin(), l.end(), [](const mcable& c) {return !test_invariants(c);});
}
+bool operator==(const msegment& l, const msegment& r) {
+ return l.prox==r.prox && l.dist==r.dist && l.tag==r.tag;
+}
+
std::ostream& operator<<(std::ostream& o, const mpoint& p) {
return o << "(point " << p.x << " " << p.y << " " << p.z << " " << p.radius << ")";
}
-std::ostream& operator<<(std::ostream& o, const msample& s) {
- return o << "(sample " << s.loc << " " << s.tag << ")";
+std::ostream& operator<<(std::ostream& o, const msegment& s) {
+ return o << "(segment " << s.prox << " " << s.dist << " " << s.tag << ")";
}
std::ostream& operator<<(std::ostream& o, const mlocation& l) {
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index 89bfce428cf09eaf90e274edbc60c544f99dafd2..6658ed453a6104fb90443fcbc2f75bb34f66f796 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -146,36 +146,24 @@ mextent thingify_(const tagged_& reg, const mprovider& p) {
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);};
+ auto matches = [reg](auto& seg){return seg.tag==reg.tag;};
+ auto not_matches = [&matches](auto& seg) {return !matches(seg);};
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.
+ auto& segs = m.branch_segments(i); // Range of segments in the branch.
+ auto locs = util::make_range(e.branch_segment_locations(i));
- 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);
+ auto beg = std::cbegin(segs);
+ auto end = std::cend(segs);
- // Find the next sample that matches reg.tag.
- auto start = std::find_if(beg+1, end, matches);
+ // Find the next section that matches reg.tag.
+ auto start = std::find_if(beg, end, [reg](auto& seg){return seg.tag==reg.tag;});
while (start!=end) {
// find the next sample that does not match reg.tag
- auto first = start-1;
auto last = std::find_if(start, end, not_matches);
- auto l = first==beg? 0.: e.sample_location(*first).pos;
- auto r = last==end? 1.: e.sample_location(*(last-1)).pos;
+ auto l = locs[start-beg].pos;
+ auto r = last==end? 1: locs[last-beg].pos;
L.push_back({i, l, r});
// Find the next sample in the branch that matches reg.tag.
diff --git a/arbor/morph/sample_tree.cpp b/arbor/morph/sample_tree.cpp
deleted file mode 100644
index 4ba6d256ac42697f4978f05a4c97922376480ea4..0000000000000000000000000000000000000000
--- a/arbor/morph/sample_tree.cpp
+++ /dev/null
@@ -1,127 +0,0 @@
-#include <iostream>
-#include <stdexcept>
-#include <vector>
-
-#include <arbor/morph/morphexcept.hpp>
-#include <arbor/morph/sample_tree.hpp>
-
-#include "io/sepval.hpp"
-#include "util/span.hpp"
-#include "util/transform.hpp"
-
-namespace arb {
-
-sample_tree::sample_tree(std::vector<msample> samples, std::vector<msize_t> parents) {
- if (samples.size()!=parents.size()) {
- throw std::invalid_argument("sample and parent vectors differ in size");
- }
- reserve(samples.size());
- for (auto i: util::make_span(samples.size())) {
- append(parents[i], samples[i]);
- }
-}
-
-void sample_tree::reserve(msize_t n) {
- samples_.reserve(n);
- parents_.reserve(n);
- props_.reserve(n);
-}
-
-msize_t sample_tree::append(msize_t p, const msample& s) {
- if ((empty() && p!=mnpos) || p>=size()) {
- if (p!=mnpos) throw invalid_sample_parent(p, size());
- }
-
- auto id = size();
- samples_.push_back(s);
- parents_.push_back(p);
-
- // Set the point properties for the new point, and update those of its parent as needed.
- point_prop prop = point_prop_mask_none;
- if (!id) {
- // if adding the first sample mark it as root
- set_root(prop);
- }
- else {
- // Mark the new node as terminal, and unset the parent sample's terminal bit.
- set_terminal(prop);
- const bool term_parent = is_terminal(props_[p]); // track if the parent was terminal.
- unset_terminal(props_[p]);
-
- // Mark if the new sample is collocated with its parent.
- if (is_collocated(s, samples_[p])) {
- set_collocated(prop);
- }
-
- // Set parent to be a fork if it was not a terminal point before the
- // new sample was added (and if it isn't the root).
- if (p && !term_parent) {
- set_fork(props_[p]);
- }
- }
- props_.push_back(prop);
-
- return id;
-}
-
-msize_t sample_tree::append(const msample& s) {
- return append(empty()? mnpos: size()-1, s);
-}
-
-msize_t sample_tree::append(msize_t p, const std::vector<msample>& slist) {
- if (!slist.size()) return size();
-
- for (auto& s: slist) {
- p = append(p, s);
- }
-
- return p;
-}
-
-msize_t sample_tree::append(const std::vector<msample>& slist) {
- return append(empty()? mnpos: size()-1, slist);
-}
-
-msize_t sample_tree::size() const {
- return samples_.size();
-}
-
-bool sample_tree::empty() const {
- return samples_.empty();
-}
-
-const std::vector<msample>& sample_tree::samples() const {
- return samples_;
-}
-
-const std::vector<msize_t>& sample_tree::parents() const {
- return parents_;
-}
-
-const std::vector<point_prop>& sample_tree::properties() const {
- return props_;
-}
-
-std::ostream& operator<<(std::ostream& o, const sample_tree& m) {
- auto tstr = util::transform_view(m.parents_,
- [](msize_t i) -> std::string {
- return i==mnpos? "npos": std::to_string(i);
- });
- return o << "(sample_tree (\n " << io::sepval(m.samples_, "\n ") << ")\n"
- << " (" << io::sepval(tstr, ' ') << "))";
-}
-
-sample_tree swc_as_sample_tree(const std::vector<swc_record>& swc_records) {
- sample_tree m;
- m.reserve(swc_records.size());
-
- for (auto i: util::count_along(swc_records)) {
- auto& r = swc_records[i];
- // The parent of soma must be mnpos, while in SWC files is -1
- msize_t p = i==0? mnpos: r.parent_id;
- m.append(p, msample{{r.x, r.y, r.z, r.r}, r.tag});
- }
- return m;
-}
-
-} // namespace arb
diff --git a/arbor/morph/segment_tree.cpp b/arbor/morph/segment_tree.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8629cef536287a475ceea878eb4737b84894f1fe
--- /dev/null
+++ b/arbor/morph/segment_tree.cpp
@@ -0,0 +1,126 @@
+#include <iostream>
+#include <stdexcept>
+#include <unordered_map>
+#include <vector>
+
+#include <arbor/morph/morphexcept.hpp>
+#include <arbor/morph/segment_tree.hpp>
+
+#include "io/sepval.hpp"
+#include "util/span.hpp"
+#include "util/transform.hpp"
+
+namespace arb {
+
+void segment_tree::reserve(msize_t n) {
+ segments_.reserve(n);
+ parents_.reserve(n);
+ seg_children_.reserve(n);
+}
+
+msize_t segment_tree::append(msize_t p, const mpoint& prox, const mpoint& dist, int tag) {
+ if (p>=size() && p!=mnpos) {
+ throw invalid_segment_parent(p, size());
+ }
+
+ auto id = size();
+ segments_.push_back(msegment{prox, dist, tag});
+ parents_.push_back(p);
+
+ // Set the point properties for the new point, and update those of the parent.
+ seg_children_.push_back({});
+ if (p!=mnpos) {
+ seg_children_[p].increment();
+ }
+
+ return id;
+}
+
+msize_t segment_tree::append(msize_t p, const mpoint& dist, int tag) {
+ // If attaching to the root both prox and dist ends must be specified.
+ if (p==mnpos) {
+ throw invalid_segment_parent(p, size());
+ }
+ if (p>=size()) {
+ throw invalid_segment_parent(p, size());
+ }
+ return append(p, segments_[p].dist, dist, tag);
+}
+
+msize_t segment_tree::size() const {
+ return segments_.size();
+}
+
+bool segment_tree::empty() const {
+ return segments_.empty();
+}
+
+const std::vector<msegment>& segment_tree::segments() const {
+ return segments_;
+}
+
+const std::vector<msize_t>& segment_tree::parents() const {
+ return parents_;
+}
+
+bool segment_tree::is_fork(msize_t i) const {
+ if (i>=size()) throw no_such_segment(i);
+ return seg_children_[i].is_fork();
+}
+bool segment_tree::is_terminal(msize_t i) const {
+ if (i>=size()) throw no_such_segment(i);
+ return seg_children_[i].is_terminal();
+}
+bool segment_tree::is_root(msize_t i) const {
+ if (i>=size()) throw no_such_segment(i);
+ return parents_[i]==mnpos;
+}
+
+std::ostream& operator<<(std::ostream& o, const segment_tree& m) {
+ auto tstr = util::transform_view(m.parents_,
+ [](msize_t i) -> std::string {
+ return i==mnpos? "npos": std::to_string(i);
+ });
+ bool one_line = m.size()<2u;
+ return o << "(segment_tree (" << (one_line? "": "\n ") << io::sepval(m.segments_, "\n ")
+ << (one_line? ") (": ")\n (")
+ << io::sepval(tstr, ' ') << "))";
+}
+
+segment_tree swc_as_segment_tree(const std::vector<swc_record>& swc_records) {
+ if (swc_records.size()<2) {
+ throw "At least two swc records are required to construct a segment tree.";
+ }
+
+ auto point = [&swc_records] (msize_t i) {
+ auto& r = swc_records[i];
+ return mpoint{r.x, r.y, r.z, r.r};
+ };
+ auto tag = [&swc_records] (msize_t i) {
+ return swc_records[i].tag;
+ };
+
+ segment_tree tree;
+ tree.reserve(swc_records.size());
+
+ std::unordered_map<msize_t, msize_t> segmap;
+
+ tree.append(mnpos, point(0), point(1), tag(1));
+ segmap[0] = mnpos;
+ segmap[1] = 0;
+ for (unsigned i=2; i<swc_records.size(); ++i) {
+ msize_t p = segmap[swc_records[i].parent_id];
+ if (p==mnpos) {
+ tree.append(p, point(0), point(i), tag(i));
+ }
+ else {
+ tree.append(p, point(i), tag(i));
+ }
+ segmap[i] = i-1;
+ }
+
+ return tree;
+}
+
+} // namespace arb
+
diff --git a/example/dryrun/branch_cell.hpp b/example/dryrun/branch_cell.hpp
index b4181c868811803391b61d8a9bf650534da68a8e..51d674f243a7dfd1160d11eec9f2b7682b6c9827 100644
--- a/example/dryrun/branch_cell.hpp
+++ b/example/dryrun/branch_cell.hpp
@@ -7,7 +7,9 @@
#include <arbor/common_types.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/morph/segment_tree.hpp>
+#include <string>
#include <sup/json_params.hpp>
// Parameters used to generate the random cell morphologies.
@@ -49,11 +51,12 @@ double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
}
arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
- arb::sample_tree tree;
+ arb::segment_tree tree;
// Add soma.
- double soma_radius = 12.6157/2.0;
- tree.append(arb::mnpos, {{0,0,0,soma_radius}, 1}); // For area of 500 μm².
+ double srad = 12.6157/2.0; // soma radius
+ int stag = 1; // soma tag
+ tree.append(arb::mnpos, {0, 0,-srad, srad}, {0, 0, srad, srad}, stag); // For area of 500 μm².
std::vector<std::vector<unsigned>> levels;
levels.push_back({0});
@@ -62,9 +65,10 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
std::mt19937 gen(gid);
std::uniform_real_distribution<double> dis(0, 1);
- double dend_radius = 0.5; // Diameter of 1 μm for each cable.
+ double drad = 0.5; // Diameter of 1 μm for each dendrite cable.
+ int dtag = 3; // Dendrite tag.
- double dist_from_soma = soma_radius;
+ double dist_from_soma = srad; // Start dendrite at the edge of the soma.
for (unsigned i=0; i<params.max_depth; ++i) {
// Branch prob at this level.
double bp = interp(params.branch_probs, i, params.max_depth);
@@ -78,14 +82,12 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
for (unsigned j=0; j<2; ++j) {
if (dis(gen)<bp) {
auto z = dist_from_soma;
- auto p = tree.append(sec, {{0,0,z,dend_radius}, 3});
- if (nc>1) {
- auto dz = l/nc;
- for (unsigned k=1; k<nc; ++k) {
- p = tree.append(p, {{0,0,z+k*dz, dend_radius}, 3});
- }
+ auto dz = l/nc;
+ auto p = sec;
+ for (unsigned k=0; k<nc; ++k) {
+ p = tree.append(p, {0,0,z+(k+1)*dz, drad}, dtag);
}
- sec_ids.push_back(tree.append(p, {{0,0,z+l,dend_radius}, 3}));
+ sec_ids.push_back(p);
}
}
}
@@ -100,20 +102,20 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
arb::label_dict d;
using arb::reg::tagged;
- d.set("soma", tagged(1));
- d.set("dendrites", join(tagged(3), tagged(4)));
+ d.set("soma", tagged(stag));
+ d.set("dend", tagged(dtag));
- arb::cable_cell cell(arb::morphology(tree, true), d);
+ arb::cable_cell cell(arb::morphology(tree), d);
cell.paint("soma", "hh");
- cell.paint("dendrites", "pas");
+ cell.paint("dend", "pas");
cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
// Add spike threshold detector at the soma.
cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
// Add a synapse to the mid point of the first dendrite.
- cell.place(arb::mlocation{1, 0.5}, "expsyn");
+ cell.place(arb::mlocation{0, 0.5}, "expsyn");
// Add additional synapses that will not be connected to anything.
for (unsigned i=1u; i<params.synapses; ++i) {
@@ -121,7 +123,7 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
}
// Make a CV between every sample in the sample tree.
- cell.default_parameters.discretization = arb::cv_policy_every_sample();
+ cell.default_parameters.discretization = arb::cv_policy_every_segment();
return cell;
}
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index a2db630be5aef7d7ccee9c9c1764fff3ca070b49..69c95b882a7cc7be017f5b243dbb1d2cd0384add 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -279,12 +279,11 @@ void write_trace_json(const std::vector<arb::trace_vector<double>>& traces, unsi
arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration) {
// Create the sample tree that defines the morphology.
- arb::sample_tree tree;
+ arb::segment_tree tree;
double soma_rad = 22.360679775/2.0; // convert diameter to radius in μm
- tree.append({{0,0,0,soma_rad}, 1}); // soma is a single sample point
- double dend_rad = 3./2;
- tree.append(0, {{0,0,soma_rad, dend_rad}, 3}); // proximal point of the dendrite
- tree.append(1, {{0,0,soma_rad+300, dend_rad}, 3}); // distal end of the dendrite
+ tree.append(arb::mnpos, {0,0,0,soma_rad}, {0,0,2*soma_rad,soma_rad}, 1); // soma
+ double dend_rad = 3./2; // μm
+ tree.append(0, {0,0,2*soma_rad, dend_rad}, {0,0,2*soma_rad+300, dend_rad}, 3); // dendrite
// Create a label dictionary that creates a single region that covers the whole cell.
arb::label_dict d;
diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp
index 5a0aa4d23f1bf40cdd4d9191162b84d01b8f8e79..66121c07b473e981bc528ce7e01dc7ae73eea2f8 100644
--- a/example/generators/generators.cpp
+++ b/example/generators/generators.cpp
@@ -47,9 +47,9 @@ public:
// capacitance: 0.01 F/m² [default]
// synapses: 1 * expsyn
arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
- arb::sample_tree tree;
+ arb::segment_tree tree;
double r = 18.8/2.0; // convert 18.8 μm diameter to radius
- tree.append({{0,0,0,r}, 1});
+ tree.append(arb::mnpos, {0,0,-r,r}, {0,0,r,r}, 1);
arb::label_dict d;
d.set("soma", arb::reg::tagged(1));
diff --git a/example/lfp/lfp.cpp b/example/lfp/lfp.cpp
index 268fd286f283d8238d67d3e873e2a8fa4b60688f..6a0066b641d4baf7f75d8f5ced9cc583b77143cd 100644
--- a/example/lfp/lfp.cpp
+++ b/example/lfp/lfp.cpp
@@ -7,9 +7,10 @@
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/place_pwlin.hpp>
#include <arbor/morph/region.hpp>
+#include <arbor/morph/segment_tree.hpp>
+#include <arbor/sampling.hpp>
#include <arbor/simple_sampler.hpp>
#include <arbor/simulation.hpp>
-#include <arbor/sampling.hpp>
#include <arbor/util/any.hpp>
#include <arbor/util/any_ptr.hpp>
@@ -51,11 +52,11 @@ struct lfp_demo_recipe: public arb::recipe {
return cell_;
}
- virtual std::vector<arb::event_generator> event_generators(cell_gid_type) const {
+ virtual std::vector<arb::event_generator> event_generators(cell_gid_type) const override {
return {events_};
}
- any get_global_properties(arb::cell_kind) const {
+ any get_global_properties(arb::cell_kind) const override {
arb::cable_cell_global_properties gprop;
gprop.default_parameters = arb::neuron_parameter_defaults;
return gprop;
@@ -70,13 +71,11 @@ private:
using namespace arb;
// Set up morphology as two branches:
+ segment_tree tree;
// * soma, length 20 μm radius 10 μm, with SWC tag 1.
+ tree.append(arb::mnpos, {0, 0, 10, 10}, {0, 0, -10, 10}, 1);
// * apical dendrite, length 490 μm, radius 1 μm, with SWC tag 4.
- sample_tree tree;
- tree.append({{0, 0, +10, 10}, 1}); // (root point)
- tree.append({{0, 0, -10, 10}, 1});
- tree.append(0, {{0, 0, 10, 1}, 4}); // attach to root point.
- tree.append({{0, 0, 500, 1}, 4});
+ tree.append(arb::mnpos, {0, 0, 10, 1}, {0, 0, 500, 1}, 4);
cell_ = cable_cell(tree);
@@ -84,8 +83,8 @@ private:
cell_.default_parameters.axial_resistivity = 100; // [Ω·cm]
cell_.default_parameters.membrane_capacitance = 0.01; // [F/m²]
- // Twenty CVs per branch, except for the soma.
- cell_.default_parameters.discretization = cv_policy_fixed_per_branch(20, cv_policy_flag::single_root_cv);
+ // Twenty CVs per branch on the dendrites (tag 4).
+ cell_.default_parameters.discretization = cv_policy_fixed_per_branch(20, arb::reg::tagged(4));
// Add pas and hh mechanisms:
cell_.paint(reg::tagged(1), "hh"); // (default parameters)
@@ -260,10 +259,9 @@ int main(int argc, char** argv) {
std::vector<std::vector<std::array<double, 3>>> samples;
for (unsigned branch = 0; branch<cell_morphology.num_branches(); ++branch) {
samples.push_back({});
- auto branch_range = cell_morphology.branch_indexes(branch);
- for (auto i_ptr = branch_range.first; i_ptr!=branch_range.second; ++i_ptr) {
- arb::msample s = cell_morphology.samples()[*i_ptr];
- samples.back().push_back(std::array<double, 3>{s.loc.x, s.loc.z, s.loc.radius});
+ for (auto& seg: cell_morphology.branch_segments(branch)) {
+ samples.back().push_back(std::array<double, 3>{seg.prox.x, seg.prox.z, seg.prox.radius});
+ samples.back().push_back(std::array<double, 3>{seg.dist.x, seg.dist.z, seg.dist.radius});
}
}
diff --git a/example/probe-demo/probe-demo.cpp b/example/probe-demo/probe-demo.cpp
index 895dc659995a2e35ba2f228cd4c7cc036c848578..0598f02888843cbadde9649f529c7efd1e73c68e 100644
--- a/example/probe-demo/probe-demo.cpp
+++ b/example/probe-demo/probe-demo.cpp
@@ -117,8 +117,9 @@ struct cable_recipe: public arb::recipe {
const double length = 1000; // [µm]
const double diam = 1; // [µm]
- sample_tree samples({msample{0, 0, 0, 0.5*diam}, msample{length, 0, 0, 0.5*diam}}, {mnpos, 0u});
- cable_cell c(samples);
+ segment_tree tree;
+ tree.append(arb::mnpos, {0, 0, 0, 0.5*diam}, {length, 0, 0, 0.5*diam}, 1);
+ cable_cell c(tree);
c.paint(reg::all(), "hh"); // HH mechanism over whole cell.
c.place(mlocation{0, 0.}, i_clamp{0., INFINITY, 1.}); // Inject a 1 nA current indefinitely.
diff --git a/example/ring/branch_cell.hpp b/example/ring/branch_cell.hpp
index b4181c868811803391b61d8a9bf650534da68a8e..51d674f243a7dfd1160d11eec9f2b7682b6c9827 100644
--- a/example/ring/branch_cell.hpp
+++ b/example/ring/branch_cell.hpp
@@ -7,7 +7,9 @@
#include <arbor/common_types.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/morph/segment_tree.hpp>
+#include <string>
#include <sup/json_params.hpp>
// Parameters used to generate the random cell morphologies.
@@ -49,11 +51,12 @@ double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
}
arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
- arb::sample_tree tree;
+ arb::segment_tree tree;
// Add soma.
- double soma_radius = 12.6157/2.0;
- tree.append(arb::mnpos, {{0,0,0,soma_radius}, 1}); // For area of 500 μm².
+ double srad = 12.6157/2.0; // soma radius
+ int stag = 1; // soma tag
+ tree.append(arb::mnpos, {0, 0,-srad, srad}, {0, 0, srad, srad}, stag); // For area of 500 μm².
std::vector<std::vector<unsigned>> levels;
levels.push_back({0});
@@ -62,9 +65,10 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
std::mt19937 gen(gid);
std::uniform_real_distribution<double> dis(0, 1);
- double dend_radius = 0.5; // Diameter of 1 μm for each cable.
+ double drad = 0.5; // Diameter of 1 μm for each dendrite cable.
+ int dtag = 3; // Dendrite tag.
- double dist_from_soma = soma_radius;
+ double dist_from_soma = srad; // Start dendrite at the edge of the soma.
for (unsigned i=0; i<params.max_depth; ++i) {
// Branch prob at this level.
double bp = interp(params.branch_probs, i, params.max_depth);
@@ -78,14 +82,12 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
for (unsigned j=0; j<2; ++j) {
if (dis(gen)<bp) {
auto z = dist_from_soma;
- auto p = tree.append(sec, {{0,0,z,dend_radius}, 3});
- if (nc>1) {
- auto dz = l/nc;
- for (unsigned k=1; k<nc; ++k) {
- p = tree.append(p, {{0,0,z+k*dz, dend_radius}, 3});
- }
+ auto dz = l/nc;
+ auto p = sec;
+ for (unsigned k=0; k<nc; ++k) {
+ p = tree.append(p, {0,0,z+(k+1)*dz, drad}, dtag);
}
- sec_ids.push_back(tree.append(p, {{0,0,z+l,dend_radius}, 3}));
+ sec_ids.push_back(p);
}
}
}
@@ -100,20 +102,20 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
arb::label_dict d;
using arb::reg::tagged;
- d.set("soma", tagged(1));
- d.set("dendrites", join(tagged(3), tagged(4)));
+ d.set("soma", tagged(stag));
+ d.set("dend", tagged(dtag));
- arb::cable_cell cell(arb::morphology(tree, true), d);
+ arb::cable_cell cell(arb::morphology(tree), d);
cell.paint("soma", "hh");
- cell.paint("dendrites", "pas");
+ cell.paint("dend", "pas");
cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
// Add spike threshold detector at the soma.
cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
// Add a synapse to the mid point of the first dendrite.
- cell.place(arb::mlocation{1, 0.5}, "expsyn");
+ cell.place(arb::mlocation{0, 0.5}, "expsyn");
// Add additional synapses that will not be connected to anything.
for (unsigned i=1u; i<params.synapses; ++i) {
@@ -121,7 +123,7 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
}
// Make a CV between every sample in the sample tree.
- cell.default_parameters.discretization = arb::cv_policy_every_sample();
+ cell.default_parameters.discretization = arb::cv_policy_every_segment();
return cell;
}
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
index dc623b8a44ec91aafd12ba0ad535e9e6c8d6741b..0a2137988ad56c9dca71c828e5827af25b5858e9 100644
--- a/example/ring/ring.cpp
+++ b/example/ring/ring.cpp
@@ -43,7 +43,7 @@ struct ring_params {
std::string name = "default";
unsigned num_cells = 10;
double min_delay = 10;
- double duration = 100;
+ double duration = 200;
cell_parameters cell;
};
diff --git a/example/single/single.cpp b/example/single/single.cpp
index ac8e67e52aea49a5eed8741136b54c7c7a04ff58..920849e3c8d87336811def28d0752ad1ee7e5dfd 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -8,7 +8,7 @@
#include <arbor/load_balance.hpp>
#include <arbor/cable_cell.hpp>
#include <arbor/morph/morphology.hpp>
-#include <arbor/morph/sample_tree.hpp>
+#include <arbor/morph/segment_tree.hpp>
#include <arbor/swcio.hpp>
#include <arbor/simulation.hpp>
#include <arbor/simple_sampler.hpp>
@@ -128,7 +128,7 @@ options parse_options(int argc, char** argv) {
opt.swc_file = swc.value();
}
else if (auto nseg = parse<unsigned>(arg, 'n', "cv-per-branch")) {
- opt.policy = arb::cv_policy_fixed_per_branch(nseg.value(), arb::cv_policy_flag::single_root_cv);
+ opt.policy = arb::cv_policy_fixed_per_branch(nseg.value());
}
else {
usage(argv[0], "[-m|--morphology SWCFILE] [-d|--dt TIME] [-t|--t-end TIME] [-w|--weight WEIGHT] [-n|--cv-per-branch N]");
@@ -144,18 +144,17 @@ options parse_options(int argc, char** argv) {
// to 0.2 µm.
arb::morphology default_morphology() {
- arb::sample_tree samples;
+ arb::segment_tree tree;
- auto p = samples.append(arb::msample{{ 0.0, 0.0, 0.0, 6.3}, 1});
- 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});
+ tree.append(arb::mnpos, { -6.3, 0.0, 0.0, 6.3}, { 6.3, 0.0, 0.0, 6.3}, 1);
+ tree.append( 0, { 6.3, 0.0, 0.0, 0.5}, {206.3, 0.0, 0.0, 0.2}, 3);
- return arb::morphology(std::move(samples));
+ return arb::morphology(tree);
}
arb::morphology read_swc(const std::string& path) {
std::ifstream f(path);
if (!f) throw std::runtime_error("unable to open SWC file: "+path);
- return arb::morphology(arb::swc_as_sample_tree(arb::parse_swc_file(f)));
+ return arb::morphology(arb::swc_as_segment_tree(arb::parse_swc_file(f)));
}
diff --git a/python/cells.cpp b/python/cells.cpp
index 89fe1fc696267858268f66ab3b10a30f576327b7..30d72bf3a78ef7178d12b8c148782be40655f07f 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -7,6 +7,7 @@
#include <arbor/morph/label_dict.hpp>
#include <arbor/morph/locset.hpp>
#include <arbor/morph/region.hpp>
+#include <arbor/morph/segment_tree.hpp>
#include <arbor/schedule.hpp>
#include <arbor/spike_source_cell.hpp>
#include <arbor/util/any.hpp>
@@ -470,12 +471,11 @@ void register_cells(pybind11::module& m) {
return arb::cable_cell(m, labels.dict);
}), "morphology"_a, "labels"_a)
.def(pybind11::init(
- [](const arb::sample_tree& t, const label_dict_proxy& labels) {
+ [](const arb::segment_tree& t, const label_dict_proxy& labels) {
return arb::cable_cell(arb::morphology(t), labels.dict);
}),
- "morphology"_a, "labels"_a,
- "Construct with a morphology derived from a sample_tree, with automatic detection of whether\n"
- "the morphology has a spherical root/soma.")
+ "segment_tree"_a, "labels"_a,
+ "Construct with a morphology derived from a segment tree.")
.def_property_readonly("num_branches",
[](const arb::cable_cell& c) {return c.morphology().num_branches();},
"The number of unbranched cable sections in the morphology.")
@@ -603,9 +603,9 @@ void register_cells(pybind11::module& m) {
[](arb::cable_cell& c, const char* label) {return c.concrete_region(label).cables();},
"label"_a, "The cable segments of the cell morphology for a region label.")
// Discretization control.
- .def("compartments_on_samples",
- [](arb::cable_cell& c) {c.default_parameters.discretization = arb::cv_policy_every_sample{};},
- "Decompose each branch into compartments defined by sample locations.")
+ .def("compartments_on_segments",
+ [](arb::cable_cell& c) {c.default_parameters.discretization = arb::cv_policy_every_segment{};},
+ "Decompose each branch into compartments defined by segments.")
.def("compartments_length",
[](arb::cable_cell& c, double len) {
c.default_parameters.discretization = arb::cv_policy_max_extent{len};
diff --git a/python/example/ring.py b/python/example/ring.py
index ee58e1e9c64e42d2a28e2454278a7eb27c2e843f..8ddaaf54e31d0efcc2d48f76bb43aee6f31c8d8c 100644
--- a/python/example/ring.py
+++ b/python/example/ring.py
@@ -15,8 +15,8 @@ import matplotlib.pyplot as plt
def make_cable_cell(gid):
b = arbor.flat_cell_builder()
- # Soma with radius 6 μm.
- s = b.add_sphere(6, "soma")
+ # Soma with radius 6 μm, modelled as cylinder of length 2*radius
+ s = b.add_cable(length=12, radius=6, name="soma", ncomp=1)
# Single dendrite of length 100 μm and radius 2 μm attached to soma.
b1 = b.add_cable(parent=s, length=100, radius=2, name="dend", ncomp=1)
# Attach two dendrites of length 50 μm to the end of the first dendrite.
diff --git a/python/example/single_cell_builder.py b/python/example/single_cell_builder.py
index 486b328a248e1a593b2117ab30f0325f2bd2f3b8..ad38d656215b51df2a1738efef7548aba265c22d 100644
--- a/python/example/single_cell_builder.py
+++ b/python/example/single_cell_builder.py
@@ -10,29 +10,30 @@ b = arbor.flat_cell_builder()
# The soma (at the root of the tree) is marked 's', and
# the end of each branch i is marked 'bi'.
#
-# b5
+# b4
# /
# /
-# b2---b4
+# b1---b3
# /
# /
-# s-------b1
+# s-------b0
# \
# \
-# b3
+# b2
-# Add a spherical soma with radius 6 um
-s = b.add_sphere(6, "soma")
+# Start with a spherical soma with radius 6 μm,
+# approximated with a cylinder of: length = diameter = 12 μm.
+s = b.add_cable(length=12, radius=6, name="soma", ncomp=1)
# Add the dendrite cables, labelling those closest to the soma "dendn",
# and those furthest with "dendx" because we will set different electrical
# properties for the two regions.
-b1 = b.add_cable(parent=s, length=100, radius=2, name="dendn", ncomp=100)
+b0 = b.add_cable(parent=s, length=100, radius=2, name="dendn", ncomp=100)
# Radius tapers from 2 to 0.5 over the length of the branch.
-b2 = b.add_cable(parent=b1, length= 50, radius=(2,0.5), name="dendn", ncomp=50)
-b3 = b.add_cable(parent=b1, length= 50, radius=1, name="dendn", ncomp=50)
-b4 = b.add_cable(parent=b2, length= 50, radius=1, name="dendx", ncomp=50)
-b5 = b.add_cable(parent=b2, length= 50, radius=1, name="dendx", ncomp=50)
+b1 = b.add_cable(parent=b0, length= 50, radius=(2,0.5), name="dendn", ncomp=50)
+b2 = b.add_cable(parent=b0, length= 50, radius=1, name="dendn", ncomp=50)
+b3 = b.add_cable(parent=b1, length= 50, radius=1, name="dendx", ncomp=50)
+b4 = b.add_cable(parent=b1, length= 50, radius=1, name="dendx", ncomp=50)
# Combine the "dendn" and "dendx" regions into a single "dend" region.
# The dendrites were labelled as such so that we can set different
@@ -40,7 +41,7 @@ b5 = b.add_cable(parent=b2, length= 50, radius=1, name="dendx", ncomp=50)
# set other properties on the whole dendrites.
b.add_label('dend', '(join (region "dendn") (region "dendx"))')
# Location of stimuli, in the middle of branch 2.
-b.add_label('stim_site', '(location 2 0.5)')
+b.add_label('stim_site', '(location 1 0.5)')
# The root of the tree (equivalent to '(location 0 0)')
b.add_label('root', '(root)')
# The tips of the dendrites (3 locations at b4, b3, b5).
diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py
index 993f6d987aef0ae5543e08133d346ab93483ed47..1cf996233d140f144fa119a4a79d230dd6cd66ea 100644
--- a/python/example/single_cell_swc.py
+++ b/python/example/single_cell_swc.py
@@ -1,3 +1,13 @@
+# NOTE: deprecating spherical roots changes the behavior of this model.
+# There is no soma, because only the root sample has tag 1, which will be
+# ignored as it is always the proximal end of any cable segment.
+# The fix is to:
+# - Write an swc interpreter that inserts a cylinder with the
+# appropriate properties.
+# - Extend the cable-only descriptions to handle detached cables, to
+# preserve surface area and correct starting locations of cables
+# attached to the soma.
+
import arbor
from arbor import mechanism as mech
from arbor import location as loc
@@ -5,19 +15,23 @@ import matplotlib.pyplot as plt
# Load a cell morphology from an swc file.
# The model has 31 branches, including soma, dendrites and axon.
-tree = arbor.load_swc('../../test/unit/swc/example.swc')
+#tree = arbor.load_swc('../../test/unit/swc/example.swc')
+tree = arbor.load_swc('example.swc')
# Define the regions and locsets in the model.
defs = {'soma': '(tag 1)', # soma has tag 1 in swc files.
'axon': '(tag 2)', # axon has tag 2 in swc files.
'dend': '(tag 3)', # dendrites have tag 3 in swc files.
'root': '(root)', # the start of the soma in this morphology is at the root of the cell.
- 'stim_site': '(location 1 0.5)'} # site for the stimulus, in the middle of branch 1.
+ 'stim_site': '(location 0 0.5)', # site for the stimulus, in the middle of branch 1.
+ 'axon_end': '(restrict (terminal) (region "axon"))'} # end of the axon.
labels = arbor.label_dict(defs)
# Combine morphology with region and locset definitions to make a cable cell.
cell = arbor.cable_cell(tree, labels)
+print(cell.locations('axon_end'))
+
# Set initial membrane potential to -55 mV
cell.set_properties(Vm=-55)
# Use Nernst to calculate reversal potential for calcium.
@@ -30,8 +44,8 @@ cell.paint('dend', 'pas')
# Increase resistivity on dendrites.
cell.paint('dend', rL=500)
# Attach stimuli that inject 0.8 nA currents for 1 ms, starting at 3 and 8 ms.
-cell.place('stim_site', arbor.iclamp(3, 1, current=0.8))
-cell.place('stim_site', arbor.iclamp(8, 1, 0.8))
+cell.place('stim_site', arbor.iclamp(3, 1, current=0.5))
+cell.place('stim_site', arbor.iclamp(8, 1, current=1))
# Detect spikes at the soma with a voltage threshold of -10 mV.
cell.place('root', arbor.spike_detector(-10))
@@ -44,8 +58,8 @@ m = arbor.single_cell_model(cell)
# Attach voltage probes that sample at 50 kHz.
m.probe('voltage', where='root', frequency=50000)
m.probe('voltage', where=loc(2,1), frequency=50000)
-m.probe('voltage', where=loc(4,1), frequency=50000)
-m.probe('voltage', where=loc(30,1), frequency=50000)
+m.probe('voltage', where='stim_site', frequency=50000)
+m.probe('voltage', where='axon_end', frequency=50000)
# Simulate the cell for 15 ms.
tfinal=15
@@ -68,7 +82,7 @@ legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces]
ax.legend(legend_labels)
ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='swc morphology demo')
plt.xlim(0,tfinal)
-plt.ylim(-80,50)
+plt.ylim(-80,80)
ax.grid()
plot_to_file=False
diff --git a/python/flat_cell_builder.cpp b/python/flat_cell_builder.cpp
index 53662f86b25a948062f5548d4d1380cff52c1a99..e201677389c864465f5f7c8e3d6dfb06ee5db345 100644
--- a/python/flat_cell_builder.cpp
+++ b/python/flat_cell_builder.cpp
@@ -4,7 +4,7 @@
#include <arbor/cable_cell.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
-#include <arbor/morph/sample_tree.hpp>
+#include <arbor/morph/segment_tree.hpp>
#include "conversion.hpp"
#include "error.hpp"
@@ -15,12 +15,9 @@
namespace pyarb {
class flat_cell_builder {
- // The sample tree describing the morphology: constructed on the fly as
- // the cables/spheres are added with add_cable/add_sphere.
- arb::sample_tree tree_;
-
- // The distal sample id of each cable.
- std::vector<arb::msize_t> cable_distal_id_;
+ // The segment tree describing the morphology, constructed additatively with
+ // segments that are attached to existing segments using add_cable.
+ arb::segment_tree tree_;
// The number of unique region names used to label cables as they
// are added to the cell.
@@ -28,6 +25,8 @@ class flat_cell_builder {
// Map from region names to the tag used to identify them.
std::unordered_map<std::string, int> tag_map_;
+ std::vector<arb::msize_t> cable_distal_segs_;
+
arb::label_dict dict_;
// The morphology is cached, and only updated on request when it is out of date.
@@ -35,26 +34,20 @@ class flat_cell_builder {
mutable arb::morphology morpho_;
mutable std::mutex mutex_;
- // Set on construction and unchanged thereafter.
- // Indicates whether
- bool spherical_ = false;
-
public:
flat_cell_builder() = default;
- arb::msize_t add_sphere(double radius, const char* name) {
- cached_morpho_ = false;
- spherical_ = true;
- if (size()) {
- throw pyarb_error("Add soma to non-empty cell.");
- }
- tree_.append({{0,0,0,radius}, get_tag(name)});
- cable_distal_id_.push_back(0);
- return 0;
+
+ // Add a new cable that is attached to the last cable added to the cell.
+ // Returns the id of the new cable.
+ arb::msize_t add_cable(double len,
+ double r1, double r2, const char* region, int ncomp)
+ {
+ return add_cable(size()? size()-1: arb::mnpos, len, r1, r2, region, ncomp);
}
- // Add a new branch that is attached to parent.
+ // Add a new cable that is attached to the parent cable.
// Returns the id of the new cable.
arb::msize_t add_cable(arb::msize_t parent, double len,
double r1, double r2, const char* region, int ncomp)
@@ -71,42 +64,23 @@ public:
int tag = get_tag(region);
const bool at_root = parent==mnpos;
- // Can't attach a cable to the root on a cell with spherical soma.
- if (at_root && spherical_) {
- throw pyarb_error("Invalid parent id.");
- }
// Parent id must be in the range [0, size())
if (!at_root && parent>=size()) {
throw pyarb_error("Invalid parent id.");
}
- // Calculating the sample that is the parent of this branch is
- // neccesarily complicated to handle cable segments that are attached
- // to the root.
- arb::msize_t p = at_root? (size()? 0: mnpos): cable_distal_id_[parent];
-
- double z = at_root? 0: // attach to root of non-spherical cell
- spherical_&&!parent? soma_rad(): // attach to spherical root
- tree_.samples()[p].loc.z; // attach to end of a cable
-
- // Only add a first point at the very beginning of the cable if
- // the cable is not attached to another
- const bool add_first_point = p==arb::mnpos // attached to the "root"
- || (!p && spherical_) // attached to a spherical root
- || (r1!=tree_.samples()[p].loc.radius);
- // proximal radius does not match r1
- if (add_first_point) {
- p = tree_.append(p, {{0,0,z,r1}, tag});
- }
- if (ncomp>1) {
- double dz = len/ncomp;
- double dr = (r2-r1)/ncomp;
- for (auto i=1; i<ncomp; ++i) {
- p = tree_.append(p, {{0,0,z+i*dz, r1+i*dr}, tag});
- }
+ arb::msize_t p = at_root? mnpos: cable_distal_segs_[parent];
+
+ double z = at_root? 0: // attach to root
+ tree_.segments()[p].dist.z; // attach to end of a cable
+
+ double dz = len/ncomp;
+ double dr = (r2-r1)/ncomp;
+ for (auto i=0; i<ncomp; ++i) {
+ p = tree_.append(p, {0,0,z+i*dz, r1+i*dr}, {0,0,z+(i+1)*dz, r1+(i+1)*dr}, tag);
}
- p = tree_.append(p, {{0,0,z+len,r2}, tag});
- cable_distal_id_.push_back(p);
+
+ cable_distal_segs_.push_back(p);
return size()-1;
}
@@ -151,7 +125,7 @@ public:
}
}
- const arb::sample_tree& samples() const {
+ const arb::segment_tree& segments() const {
return tree_;
}
@@ -170,16 +144,15 @@ public:
const arb::morphology& morphology() const {
const std::lock_guard<std::mutex> guard(mutex_);
if (!cached_morpho_) {
- morpho_ = arb::morphology(tree_, spherical_);
+ morpho_ = arb::morphology(tree_);
cached_morpho_ = true;
}
return morpho_;
}
arb::cable_cell build() const {
- // Make cable_cell from sample tree and dictionary.
auto c = arb::cable_cell(morphology(), dict_);
- c.default_parameters.discretization = arb::cv_policy_every_sample{};
+ c.default_parameters.discretization = arb::cv_policy_every_segment{};
return c;
}
@@ -217,14 +190,9 @@ public:
}
}
- // Only valid if called on a non-empty tree with spherical soma.
- double soma_rad() const {
- return tree_.samples()[0].loc.radius;
- }
-
- // The number of cable segements (plus one optional soma) used to construct the cell.
+ // The number of cable segements used in the cell.
std::size_t size() const {
- return cable_distal_id_.size();
+ return tree_.size();
}
};
@@ -234,8 +202,24 @@ void register_flat_builder(pybind11::module& m) {
pybind11::class_<flat_cell_builder> builder(m, "flat_cell_builder");
builder
.def(pybind11::init<>())
- .def("add_sphere", &flat_cell_builder::add_sphere,
- "radius"_a, "name"_a)
+ .def("add_cable",
+ [](flat_cell_builder& b, double len, pybind11::object rad, const char* name, int ncomp) {
+ using pybind11::isinstance;
+ using pybind11::cast;
+ if (auto radius = try_cast<double>(rad) ) {
+ return b.add_cable(len, *radius, *radius, name, ncomp);
+ }
+
+ if (auto radii = try_cast<std::pair<double, double>>(rad)) {
+ return b.add_cable(len, radii->first, radii->second, name, ncomp);
+ }
+ else {
+ throw pyarb_error(
+ "Radius parameter is not a scalar (constant branch radius) or "
+ "a tuple (radius at proximal and distal ends respectively).");
+ }
+ },
+ "length"_a, "radius"_a, "name"_a, "ncomp"_a=1)
.def("add_cable",
[](flat_cell_builder& b, arb::msize_t p, double len, pybind11::object rad, const char* name, int ncomp) {
using pybind11::isinstance;
@@ -256,8 +240,8 @@ void register_flat_builder(pybind11::module& m) {
"parent"_a, "length"_a, "radius"_a, "name"_a, "ncomp"_a=1)
.def("add_label", &flat_cell_builder::add_label,
"name"_a, "description"_a)
- .def_property_readonly("samples",
- [](const flat_cell_builder& b) { return b.samples(); })
+ .def_property_readonly("segments",
+ [](const flat_cell_builder& b) { return b.segments(); })
.def_property_readonly("labels",
[](const flat_cell_builder& b) { return b.labels(); })
.def_property_readonly("morphology",
diff --git a/python/morph_parse.cpp b/python/morph_parse.cpp
index ba35b9a0c4a3bffc56a76be49745b28a0cdac8ab..2df63fc59bde458b682b8f923fa2a5194b69db47 100644
--- a/python/morph_parse.cpp
+++ b/python/morph_parse.cpp
@@ -227,8 +227,6 @@ std::unordered_multimap<std::string, evaluator> eval_map {
"'location' with 2 arguments: (branch_id:integer position:real)")},
{"terminal", make_call<>(arb::ls::terminal,
"'terminal' with 0 arguments")},
- {"sample", make_call<int>(arb::ls::sample,
- "'sample' with 1 argument: (sample_id:integer)")},
{"distal", make_call<arb::region>(arb::ls::most_distal,
"'distal' with 1 argument: (reg:region)")},
{"proximal",make_call<arb::region>(arb::ls::most_proximal,
diff --git a/python/morphology.cpp b/python/morphology.cpp
index ea2ea8101b9ea64c632782f16444dab0fb5cd1fd..ec634edd20673fe187320fc428bd01934f6f5d20 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -5,7 +5,7 @@
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
-#include <arbor/morph/sample_tree.hpp>
+#include <arbor/morph/segment_tree.hpp>
#include <arbor/swcio.hpp>
#include "error.hpp"
@@ -17,7 +17,7 @@ void register_morphology(pybind11::module& m) {
using namespace pybind11::literals;
//
- // primitives: points, samples, locations, cables... etc.
+ // primitives: points, segments, locations, cables... etc.
//
m.attr("mnpos") = arb::mnpos;
@@ -45,37 +45,32 @@ void register_morphology(pybind11::module& m) {
.def("__repr__",
[](arb::mlocation l) { return util::pprintf("(location {} {})", l.branch, l.pos); });
- // arb::msample
- pybind11::class_<arb::msample> msample(m, "sample");
- msample
+ // arb::mpoint
+ pybind11::class_<arb::mpoint> mpoint(m, "mpoint");
+ mpoint
.def(pybind11::init(
- [](arb::mpoint loc, int tag){
- return arb::msample{loc, tag};}),
- "location"_a, "tag"_a)
- .def(pybind11::init(
- [](double x, double y, double z, double r, int tag){
- return arb::msample{{x,y,z,r}, tag};}),
- "x"_a, "y"_a, "z"_a, "radius"_a, "tag"_a,
- "spatial values {x, y, z, radius} are in μm.")
- .def_readonly("x", &arb::msample::loc)
- .def_property_readonly("x",
- [](const arb::msample& s) {return s.loc.x;},
- "X coordinate [μm].")
- .def_property_readonly("y",
- [](const arb::msample& s) {return s.loc.y;},
- "Y coordinate [μm].")
- .def_property_readonly("z",
- [](const arb::msample& s) {return s.loc.z;},
- "Z coordinate [μm].")
- .def_readonly("tag", &arb::msample::tag,
- "Property tag of sample. "
- "If loaded from standard SWC file the following tags are used: soma=1, axon=2, dendrite=3, apical dendrite=4, however arbitrary tags can be used.")
+ [](double x, double y, double z, double r) {
+ return arb::mpoint{x,y,z,r};
+ }),
+ "x"_a, "y"_a, "z"_a, "radius"_a, "All values in μm.")
+ .def_readonly("x", &arb::mpoint::x, "X coordinate [μm].")
+ .def_readonly("y", &arb::mpoint::y, "Y coordinate [μm].")
+ .def_readonly("z", &arb::mpoint::z, "Z coordinate [μm].")
+ .def_readonly("radius", &arb::mpoint::radius,
+ "Radius of cable at sample location centered at coordinates [μm].")
.def("__str__",
- [](const arb::msample& s) {
- return util::pprintf("{}", s);})
+ [](const arb::mpoint& p) {
+ return util::pprintf("<arbor.mpoint: x {}, y {}, z {}, radius {}>", p.x, p.y, p.z, p.radius);
+ })
.def("__repr__",
- [](const arb::msample& s) {
- return util::pprintf("{}", s);});
+ [](const arb::mpoint& p) {return util::pprintf("{}>", p);});
+
+ // arb::msegment
+ pybind11::class_<arb::msegment> msegment(m, "msegment");
+ msegment
+ .def_readonly("prox", &arb::msegment::prox, "the location and radius of the proximal end.")
+ .def_readonly("dist", &arb::msegment::dist, "the location and radius of the distal end.")
+ .def_readonly("tag", &arb::msegment::tag, "tag meta-data.");
// arb::mcable
pybind11::class_<arb::mcable> cable(m, "cable");
@@ -99,47 +94,46 @@ void register_morphology(pybind11::module& m) {
.def("__repr__", [](const arb::mcable& c) { return util::pprintf("{}", c); });
//
- // Higher-level data structures (sample_tree, morphology)
+ // Higher-level data structures (segment_tree, morphology)
//
- // arb::sample_tree
- pybind11::class_<arb::sample_tree> sample_tree(m, "sample_tree");
- sample_tree
+ // arb::segment_tree
+ pybind11::class_<arb::segment_tree> segment_tree(m, "segment_tree");
+ segment_tree
// constructors
.def(pybind11::init<>())
// modifiers
- .def("reserve", &arb::sample_tree::reserve)
- .def("append", [](arb::sample_tree& t, arb::msample s){return t.append(s);},
- "Append a sample whose parent is the last sample added to the tree.")
- .def("append", [](arb::sample_tree& t, arb::msize_t p, arb::msample s){return t.append(p, s);},
- "parent"_a, "sample"_a,
- "Append a sample.")
+ .def("reserve", &arb::segment_tree::reserve)
+ .def("append", [](arb::segment_tree& t, arb::msize_t parent, arb::mpoint prox, arb::mpoint dist, int tag) {
+ return t.append(parent, prox, dist, tag);
+ },
+ "parent"_a, "prox"_a, "dist"_a, "tag"_a,
+ "Append a segment to the tree.")
+ .def("append", [](arb::segment_tree& t, arb::msize_t parent, arb::mpoint dist, int tag) {
+ return t.append(parent, dist, tag);
+ },
+ "parent"_a, "dist"_a, "tag"_a,
+ "Append a segment to the tree.")
.def("append",
- [](arb::sample_tree& t, double x, double y, double z, double radius, int tag) {
- return t.append(arb::msample{{x,y,z,radius}, tag});
- },
- "x"_a, "y"_a, "z"_a, "radius"_a, "tag"_a,
- "Append a sample whose parent is the last sample added to the tree.")
- .def("append",
- [](arb::sample_tree& t, arb::msize_t p, double x, double y, double z, double radius, int tag) {
- return t.append(p, arb::msample{{x,y,z,radius}, tag});
+ [](arb::segment_tree& t, arb::msize_t p, double x, double y, double z, double radius, int tag) {
+ return t.append(p, arb::mpoint{x,y,z,radius}, tag);
},
"parent"_a, "x"_a, "y"_a, "z"_a, "radius"_a, "tag"_a,
- "Append a sample.")
+ "Append a segment to the tree, using the distal location of the parent segment as the proximal end.")
// properties
- .def_property_readonly("empty", [](const arb::sample_tree& st){return st.empty();},
+ .def_property_readonly("empty", [](const arb::segment_tree& st){return st.empty();},
"Indicates whether the sample tree is empty (i.e. whether it has size 0)")
- .def_property_readonly("size", [](const arb::sample_tree& st){return st.size();},
+ .def_property_readonly("size", [](const arb::segment_tree& st){return st.size();},
"The number of samples in the sample tree.")
- .def_property_readonly("parents", [](const arb::sample_tree& st){return st.parents();},
+ .def_property_readonly("parents", [](const arb::segment_tree& st){return st.parents();},
"A list with the parent index of each sample.")
- .def_property_readonly("samples", [](const arb::sample_tree& st){return st.samples();},
+ .def_property_readonly("segments", [](const arb::segment_tree& st){return st.segments();},
"A list of the samples.")
- .def("__str__", [](const arb::sample_tree& s) {
- return util::pprintf("<arbor.sample_tree:\n{}>", s);});
+ .def("__str__", [](const arb::segment_tree& s) {
+ return util::pprintf("<arbor.segment_tree:\n{}>", s);});
- // Function that creates a sample_tree from an swc file.
- // Wraps calls to C++ functions arb::parse_swc_file() and arb::swc_as_sample_tree().
+ // Function that creates a segment_tree from an swc file.
+ // Wraps calls to C++ functions arb::parse_swc_file() and arb::swc_as_segment_tree().
m.def("load_swc",
[](std::string fname) {
std::ifstream fid{fname};
@@ -149,7 +143,7 @@ void register_morphology(pybind11::module& m) {
try {
auto records = arb::parse_swc_file(fid);
arb::swc_canonicalize(records);
- return arb::swc_as_sample_tree(records);
+ return arb::swc_as_segment_tree(records);
}
catch (arb::swc_error& e) {
// Try to produce helpful error messages for SWC parsing errors.
@@ -158,7 +152,7 @@ void register_morphology(pybind11::module& m) {
e.line_number, fname, e.what()));
}
},
- "Load an swc file and convert to a sample_tree.");
+ "Load an swc file and convert to a segment_tree.");
// arb::morphology
@@ -166,43 +160,26 @@ void register_morphology(pybind11::module& m) {
morph
// constructors
.def(pybind11::init(
- [](arb::sample_tree t){
+ [](arb::segment_tree t){
return arb::morphology(std::move(t));
}))
- .def(pybind11::init(
- [](arb::sample_tree t, bool spherical_root) {
- return arb::morphology(std::move(t), spherical_root);
- }), "sample_tree"_a, "spherical_root"_a)
// morphology's interface is read-only by design, so most of it can
// be implemented as read-only properties.
.def_property_readonly("empty",
[](const arb::morphology& m){return m.empty();},
"Whether the morphology is empty.")
- .def_property_readonly("spherical_root",
- [](const arb::morphology& m){return m.spherical_root();},
- "Whether the root of the morphology is spherical.")
.def_property_readonly("num_branches",
[](const arb::morphology& m){return m.num_branches();},
"The number of branches in the morphology.")
- .def_property_readonly("num_samples",
- [](const arb::morphology& m){return m.num_samples();},
- "The number of samples in the morphology.")
- .def_property_readonly("samples",
- [](const arb::morphology& m){return m.samples();},
- "All of the samples in the morphology.")
- .def_property_readonly("sample_parents",
- [](const arb::morphology& m){return m.sample_parents();},
- "The parent indexes of each sample.")
.def("branch_parent", &arb::morphology::branch_parent,
"i"_a, "The parent branch of branch i.")
.def("branch_children", &arb::morphology::branch_children,
"i"_a, "The child branches of branch i.")
- .def("branch_indexes",
+ .def("branch_segments",
[](const arb::morphology& m, arb::msize_t i) {
- auto p = m.branch_indexes(i);
- return std::vector<arb::msize_t>(p.first, p.second);
+ return m.branch_segments(i);
},
- "i"_a, "Range of indexes into the sample points in branch i.")
+ "i"_a, "A list of the segments in branch i, ordered from proximal to distal ends of the branch.")
.def("__str__",
[](const arb::morphology& m) {
return util::pprintf("<arbor.morphology:\n{}>", m);
diff --git a/test/common_cells.cpp b/test/common_cells.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..43cbcc92da557c80afc1a31ec0705fff6c744e02
--- /dev/null
+++ b/test/common_cells.cpp
@@ -0,0 +1,262 @@
+#include "common_cells.hpp"
+
+namespace arb {
+
+// Generate a segment tree from a sequence of points and parent index.
+arb::segment_tree segments_from_points(
+ std::vector<arb::mpoint> points,
+ std::vector<arb::msize_t> parents,
+ std::vector<int> tags)
+{
+ using arb::mnpos;
+ using arb::msize_t;
+
+ const auto np = points.size();
+
+ // Sanity check the input.
+ if (parents.size()!=np || np<2) {
+ throw std::runtime_error("segments_from_points: every point must have a parent.");
+ }
+ if (tags.size() && tags.size()!=np) {
+ throw std::runtime_error("segments_from_points: every point must have a tag.");
+ }
+ auto tag = [&tags](msize_t i) {return tags.size()? tags[i]: 1;};
+ arb::segment_tree tree;
+ std::unordered_map<msize_t, msize_t> segmap;
+
+ tree.append(mnpos, points[0], points[1], tag(1));
+ segmap[0] = mnpos;
+ segmap[1] = 0;
+ for (unsigned i=2; i<np; ++i) {
+ auto p = segmap[parents[i]];
+ if (p==mnpos) {
+ tree.append(p, points[0], points[i], tag(i));
+ }
+ else {
+ tree.append(p, points[i], tag(i));
+ }
+ segmap[i] = i-1;
+ }
+
+ return tree;
+}
+
+
+int soma_cell_builder::get_tag(const std::string& name) {
+ auto it = tag_map.find(name);
+ // If the name is not in the map, make a unique tag.
+ // Tags start from 1.
+ if (it==tag_map.end()) {
+ tag_map[name] = ++tag_count;
+ return tag_count;
+ }
+ return it->second;
+}
+
+soma_cell_builder::soma_cell_builder(double r) {
+ auto tag = get_tag("soma");
+ tree.append(arb::mnpos, {0,0,0,r}, {0,0,2*r,r}, tag);
+ cv_boundaries.push_back({0, 1.});
+ branch_distal_id.push_back(0);
+ n_soma_children = 0;
+}
+
+mlocation soma_cell_builder::location(mlocation loc) const {
+ if (loc.branch>=branch_distal_id.size()) {
+ throw cable_cell_error("location not on cable_cell.");
+ }
+ if (n_soma_children==0) {
+ return loc;
+ }
+ if (n_soma_children==1) {
+ if (loc.branch<2) {
+ double soma_len = tree.segments()[branch_distal_id[0]].dist.z;
+ double total_len = tree.segments()[branch_distal_id[1]].dist.z;
+ // relative position of the end of the soma on the first branch
+ double split = soma_len/total_len;
+ double pos = loc.branch==0?
+ split*loc.pos:
+ split + (1-split)*loc.pos;
+ return {0, pos};
+ }
+ return {loc.branch-1, loc.pos};
+ }
+ return loc;
+}
+
+mcable soma_cell_builder::cable(mcable cab) const {
+ if (cab.branch>=branch_distal_id.size()) {
+ throw cable_cell_error("cable not on cable_cell.");
+ }
+ auto beg = location({cab.branch, cab.prox_pos});
+ auto end = location({cab.branch, cab.dist_pos});
+ return {beg.branch, beg.pos, end.pos};
+}
+
+// Add a new branch that is attached to parent_branch.
+// Returns the id of the new branch.
+msize_t soma_cell_builder::add_branch(
+ msize_t parent_branch,
+ double len, double r1, double r2, int ncomp,
+ const std::string& region)
+{
+ // Get tag id of region (add a new tag if region does not already exist).
+ int tag = get_tag(region);
+
+ msize_t p = branch_distal_id[parent_branch];
+ auto& ploc = tree.segments()[p].dist;
+
+ double z = ploc.z;
+ if (ploc.radius!=r1) {
+ p = tree.append(p, {0,0,z,r1}, tag);
+ }
+ if (ncomp>1) {
+ double dz = len/ncomp;
+ double dr = (r2-r1)/ncomp;
+ for (auto i=1; i<ncomp; ++i) {
+ p = tree.append(p, {0,0,z+i*dz, r1+i*dr}, tag);
+ }
+ }
+ p = tree.append(p, {0,0,z+len,r2}, tag);
+ branch_distal_id.push_back(p);
+
+ msize_t bid = branch_distal_id.size()-1;
+ for (int i = 0; i<ncomp; ++i) {
+ cv_boundaries.push_back(mlocation{bid, (2*i+1.)/(2.*ncomp)});
+ }
+
+ if (!parent_branch) ++n_soma_children;
+
+ return bid;
+}
+
+cable_cell soma_cell_builder::make_cell() const {
+ // Test that a valid tree was generated, that is, every branch has
+ // either 0 children, or at least 2 children.
+ for (auto i: branch_distal_id) {
+ // skip soma
+ if (i<2) continue;
+ if (!tree.is_fork(i) && !tree.is_terminal(i)) {
+ throw cable_cell_error(
+ "attempt to construct a cable_cell from a soma_cell_builder "
+ "where a non soma branch has only one child branch.");
+ }
+ }
+
+ // Make label dictionary with one entry for each tag.
+ label_dict dict;
+ for (auto& tag: tag_map) {
+ dict.set(tag.first, reg::tagged(tag.second));
+ }
+
+ // Make cable_cell from sample tree and dictionary.
+ cable_cell c(tree, dict);
+ auto boundaries = cv_boundaries;
+ for (auto& b: boundaries) {
+ b = location(b);
+ }
+ c.default_parameters.discretization = cv_policy_explicit(boundaries);
+ return c;
+}
+
+/*
+ * Create cell with just a soma:
+ *
+ * Soma:
+ * diameter: 18.8 µm
+ * mechanisms: HH (default params)
+ * bulk resistivitiy: 100 Ω·cm [default]
+ * capacitance: 0.01 F/m² [default]
+ *
+ * Stimuli:
+ * soma centre, t=[10 ms, 110 ms), 0.1 nA
+ */
+
+cable_cell make_cell_soma_only(bool with_stim) {
+ soma_cell_builder builder(18.8/2.0);
+
+ auto c = builder.make_cell();
+ c.paint("soma", "hh");
+ if (with_stim) {
+ c.place(builder.location({0,0.5}), i_clamp{10., 100., 0.1});
+ }
+
+ return c;
+}
+
+/*
+ * Create cell with a soma and unbranched dendrite:
+ *
+ * Common properties:
+ * bulk resistivity: 100 Ω·cm [default]
+ * capacitance: 0.01 F/m² [default]
+ *
+ * Soma:
+ * mechanisms: HH (default params)
+ * diameter: 12.6157 µm
+ *
+ * Dendrite:
+ * mechanisms: passive (default params)
+ * diameter: 1 µm
+ * length: 200 µm
+ * compartments: 4
+ *
+ * Stimulus:
+ * end of dendrite, t=[5 ms, 85 ms), 0.3 nA
+ */
+
+cable_cell make_cell_ball_and_stick(bool with_stim) {
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
+
+ auto c = builder.make_cell();
+ c.paint("soma", "hh");
+ c.paint("dend", "pas");
+ if (with_stim) {
+ c.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
+ }
+
+ return c;
+}
+
+/*
+ * Create cell with a soma and three-branch dendrite with single branch point:
+ *
+ * O----======
+ *
+ * Common properties:
+ * bulk resistivity: 100 Ω·cm [default]
+ * capacitance: 0.01 F/m² [default]
+ *
+ * Soma:
+ * mechanisms: HH (default params)
+ * diameter: 12.6157 µm
+ *
+ * Dendrites:
+ * mechanisms: passive (default params)
+ * diameter: 1 µm
+ * length: 100 µm
+ * compartments: 4
+ *
+ * Stimulus:
+ * end of first terminal branch, t=[5 ms, 85 ms), 0.45 nA
+ * end of second terminal branch, t=[40 ms, 50 ms), -0.2 nA
+ */
+
+cable_cell make_cell_ball_and_3stick(bool with_stim) {
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
+ builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
+ builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
+
+ auto c = builder.make_cell();
+ c.paint("soma", "hh");
+ c.paint("dend", "pas");
+ if (with_stim) {
+ c.place(builder.location({2,1}), i_clamp{5., 80., 0.45});
+ c.place(builder.location({3,1}), i_clamp{40., 10.,-0.2});
+ }
+
+ return c;
+}
+} // namespace arb
diff --git a/test/common_cells.hpp b/test/common_cells.hpp
index 83d972ae9ad7f953df0a4f14dc662196488ceebe..bad660a9b0d9434b4e4d89dd97e5637e9026c0a4 100644
--- a/test/common_cells.hpp
+++ b/test/common_cells.hpp
@@ -7,86 +7,35 @@
namespace arb {
+// Generate a segment tree from a sequence of points and parent index.
+arb::segment_tree segments_from_points(std::vector<arb::mpoint> points,
+ std::vector<arb::msize_t> parents,
+ std::vector<int> tags={});
+
class soma_cell_builder {
- double soma_rad;
- sample_tree tree;
+ segment_tree tree;
std::vector<msize_t> branch_distal_id;
std::unordered_map<std::string, int> tag_map;
- locset cv_boundaries = mlocation{0, 1.};
+ mlocation_list cv_boundaries;
int tag_count = 0;
+ int n_soma_children = 0;
// Get tag id of region.
// Add a new tag if region with that name has not already had a tag associated with it.
- int get_tag(const std::string& name) {
- auto it = tag_map.find(name);
- // If the name is not in the map, make a unique tag.
- // Tags start from 1.
- if (it==tag_map.end()) {
- tag_map[name] = ++tag_count;
- return tag_count;
- }
- return it->second;
- }
+ int get_tag(const std::string& name);
public:
- soma_cell_builder(double r): soma_rad(r) {
- tree.append({{0,0,0,r}, get_tag("soma")});
- branch_distal_id.push_back(0);
- }
+ soma_cell_builder(double r);
// Add a new branch that is attached to parent_branch.
// Returns the id of the new branch.
msize_t add_branch(msize_t parent_branch, double len, double r1, double r2, int ncomp,
- const std::string& region)
- {
- // Get tag id of region (add a new tag if region does not already exist).
- int tag = get_tag(region);
-
- msize_t p = branch_distal_id[parent_branch];
- double z = parent_branch? tree.samples()[p].loc.z: soma_rad;
-
- p = tree.append(p, {{0,0,z,r1}, tag});
- if (ncomp>1) {
- double dz = len/ncomp;
- double dr = (r2-r1)/ncomp;
- for (auto i=1; i<ncomp; ++i) {
- p = tree.append(p, {{0,0,z+i*dz, r1+i*dr}, tag});
- }
- }
- p = tree.append(p, {{0,0,z+len,r2}, tag});
- branch_distal_id.push_back(p);
-
- msize_t bid = branch_distal_id.size()-1;
- for (int i = 0; i<ncomp; ++i) {
- cv_boundaries = sum(cv_boundaries, mlocation{bid, (2*i+1.)/(2.*ncomp)});
- }
- return bid;
- }
-
- cable_cell make_cell() const {
- // Test that a valid tree was generated, that is, every branch has
- // either 0 children, or at least 2 children.
- for (auto i: branch_distal_id) {
- if (i==0) continue;
- auto prop = tree.properties()[i];
- if (!is_fork(prop) && !is_terminal(prop)) {
- throw cable_cell_error(
- "attempt to construct a cable_cell from a soma_cell_builder "
- "where a branch has only one child branch.");
- }
- }
-
- // Make label dictionary with one entry for each tag.
- label_dict dict;
- for (auto& tag: tag_map) {
- dict.set(tag.first, reg::tagged(tag.second));
- }
-
- // Make cable_cell from sample tree and dictionary.
- cable_cell c(tree, dict);
- c.default_parameters.discretization = cv_policy_explicit(cv_boundaries);
- return c;
- }
+ const std::string& region);
+
+ mlocation location(mlocation) const;
+ mcable cable(mcable) const;
+
+ cable_cell make_cell() const;
};
/*
@@ -102,17 +51,7 @@ public:
* soma centre, t=[10 ms, 110 ms), 0.1 nA
*/
-inline cable_cell make_cell_soma_only(bool with_stim = true) {
- soma_cell_builder builder(18.8/2.0);
-
- auto c = builder.make_cell();
- c.paint("soma", "hh");
- if (with_stim) {
- c.place(mlocation{0,0.5}, i_clamp{10., 100., 0.1});
- }
-
- return c;
-}
+cable_cell make_cell_soma_only(bool with_stim = true);
/*
* Create cell with a soma and unbranched dendrite:
@@ -135,19 +74,7 @@ inline cable_cell make_cell_soma_only(bool with_stim = true) {
* end of dendrite, t=[5 ms, 85 ms), 0.3 nA
*/
-inline cable_cell make_cell_ball_and_stick(bool with_stim = true) {
- soma_cell_builder builder(12.6157/2.0);
- builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
-
- auto c = builder.make_cell();
- c.paint("soma", "hh");
- c.paint("dend", "pas");
- if (with_stim) {
- c.place(mlocation{1,1}, i_clamp{5, 80, 0.3});
- }
-
- return c;
-}
+cable_cell make_cell_ball_and_stick(bool with_stim = true);
/*
* Create cell with a soma and three-branch dendrite with single branch point:
@@ -173,21 +100,6 @@ inline cable_cell make_cell_ball_and_stick(bool with_stim = true) {
* end of second terminal branch, t=[40 ms, 50 ms), -0.2 nA
*/
-inline cable_cell make_cell_ball_and_3stick(bool with_stim = true) {
- soma_cell_builder builder(12.6157/2.0);
- builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
- builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
- builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
-
- auto c = builder.make_cell();
- c.paint("soma", "hh");
- c.paint("dend", "pas");
- if (with_stim) {
- c.place(mlocation{2,1}, i_clamp{5., 80., 0.45});
- c.place(mlocation{3,1}, i_clamp{40., 10.,-0.2});
- }
-
- return c;
-}
+cable_cell make_cell_ball_and_3stick(bool with_stim = true);
} // namespace arb
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index c13aa40bc82116bf3c43122b6cc00068e6acad18..ee680570c7fef37f59ead714d18047c889efd6bf 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -68,7 +68,7 @@ endforeach()
# set(proto_mechanisms pas hh expsyn exp2syn test_kin1 test_kinlva test_ca)
# set(mech_proto_dir "${CMAKE_CURRENT_BINARY_DIR}/mech_proto")
# file(MAKE_DIRECTORY "${mech_proto_dir}")
-#
+#
# build_modules(
# ${proto_mechanisms}
# SOURCE_DIR "${PROJECT_SOURCE_DIR}/mechanisms/mod"
@@ -82,6 +82,7 @@ endforeach()
# Unit test sources
set(unit_sources
+ ../common_cells.cpp
test_algorithms.cpp
test_any.cpp
test_any_visitor.cpp
@@ -138,10 +139,11 @@ set(unit_sources
test_range.cpp
test_ratelem.cpp
test_schedule.cpp
- test_spike_source.cpp
test_scope_exit.cpp
+ test_segment_tree.cpp
test_simd.cpp
test_span.cpp
+ test_spike_source.cpp
test_spikes.cpp
test_spike_store.cpp
test_stats.cpp
diff --git a/test/unit/common.hpp b/test/unit/common.hpp
index 54db33d3aa75436913a37044efd1f244badd8019..00f569caa4f39a513f551b6b8ee79f8522d75d4c 100644
--- a/test/unit/common.hpp
+++ b/test/unit/common.hpp
@@ -8,7 +8,9 @@
#include <cmath>
#include <string>
#include <type_traits>
+#include <unordered_map>
#include <utility>
+#include <vector>
#include "../gtest.h"
@@ -16,7 +18,7 @@
namespace std {
template <typename A, typename B>
- std::ostream& operator<<(std::ostream& out, const std::pair<A, B>& p) {
+ ::std::ostream& operator<<(::std::ostream& out, const ::std::pair<A, B>& p) {
return out << '(' << p.first << ',' << p.second << ')';
}
}
diff --git a/test/unit/common_morphologies.hpp b/test/unit/common_morphologies.hpp
index de00aac7e494db36fb21f6fadf3da04c810a4929..46a7dd4e5216e9e740dc1c98a3ed9987b1100b34 100644
--- a/test/unit/common_morphologies.hpp
+++ b/test/unit/common_morphologies.hpp
@@ -2,16 +2,35 @@
// A set of morphologies for testing discretization.
+#include "util/span.hpp"
+
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/segment_tree.hpp>
+
+#include <cstring>
#include <utility>
#include <vector>
-#include <arbor/morph/morphology.hpp>
+
+#include "io/sepval.hpp"
namespace common_morphology {
-inline std::vector<arb::msample> make_morph_samples(unsigned n) {
- std::vector<arb::msample> ms;
- for (unsigned i = 0; i<n; ++i) ms.push_back({{0., 0., (double)i, 0.5}, 5});
- return ms;
+inline arb::morphology make_morph(std::vector<arb::msize_t> parents, const char* name="") {
+ auto nseg = parents.size();;
+ auto point = [](int i) {return arb::mpoint{0, 0, (double)i, 0.5};};
+ int tag = 1;
+
+ if (!parents.size()) return {};
+
+ arb::segment_tree tree;
+ tree.append(arb::mnpos, point(-1), point(0), tag);
+ for (arb::msize_t i=1; i<nseg; ++i) {
+ int p = parents[i]==arb::mnpos? -1: parents[i];
+ tree.append(parents[i], point(p), point(i), tag);
+ }
+
+ return arb::morphology(tree);
}
// Test morphologies for CV determination:
@@ -20,16 +39,16 @@ inline std::vector<arb::msample> make_morph_samples(unsigned n) {
static const arb::morphology m_empty;
-// One branch.
-static const arb::morphology m_reg_b1{arb::sample_tree(make_morph_samples(2), {arb::mnpos, 0u}), false};
+// regular root, one branch
+static const arb::morphology m_reg_b1 = make_morph({arb::mnpos});
// Six branches:
// branch 0 has child branches 1 and 2; branch 2 has child branches 3, 4 and 5.
-static const arb::morphology m_reg_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 2u, 2u, 2u}), false};
+static const arb::morphology m_reg_b6 = make_morph({arb::mnpos, 0u, 0u, 1u, 1u, 1u});
// Six branches, mutiple top level branches:
// branch 0 has child branches 1 and 2; branch 3 has child branches 4 and 5.
-static const arb::morphology m_mlt_b6{arb::sample_tree(make_morph_samples(7), {arb::mnpos, 0u, 1u, 1u, 0u, 4u, 4u}), false};
+static const arb::morphology m_mlt_b6 = make_morph({arb::mnpos, 0u, 0u, arb::mnpos, 3u, 3u});
static std::pair<const char*, arb::morphology> test_morphologies[] = {
{"m_empty", m_empty},
diff --git a/test/unit/test_cv_policy.cpp b/test/unit/test_cv_policy.cpp
index bba0fa1531736602163a966dfae4eb70365375ce..2ffbfbc3afe4281020b6ccbbe2e12a9998d7031f 100644
--- a/test/unit/test_cv_policy.cpp
+++ b/test/unit/test_cv_policy.cpp
@@ -15,6 +15,7 @@
#include "util/rangeutil.hpp"
#include "util/span.hpp"
+#include "common_cells.hpp"
#include "common_morphologies.hpp"
#include "morph_pred.hpp"
@@ -266,26 +267,27 @@ TEST(cv_policy, max_extent) {
}
}
-TEST(cv_policy, every_sample) {
+TEST(cv_policy, every_segment) {
using namespace cv_policy_flag;
// Cell with root branch and two child branches, with multiple samples per branch.
// Fork is at (0., 0., 4.0).
- std::vector<msample> ms;
+ std::vector<mpoint> points;
- ms.push_back({{ 0., 0., 0., 0.5}, 5});
- for (auto i: make_span(4)) ms.push_back({{ 0., 0., i+1., 0.5}, 5});
- for (auto i: make_span(4)) ms.push_back({{ 0., i+1., 4.0, 0.5}, 5});
- for (auto i: make_span(4)) ms.push_back({{i+1., 0, 4.0, 0.5}, 5});
+ points.push_back({ 0., 0., 0., 0.5});
+ for (auto i: make_span(4)) points.push_back({ 0., 0., i+1., 0.5});
+ for (auto i: make_span(4)) points.push_back({ 0., i+1., 4.0, 0.5});
+ for (auto i: make_span(4)) points.push_back({i+1., 0, 4.0, 0.5});
std::vector<msize_t> parents = {mnpos, 0, 1, 2, 3, 4, 5, 6, 7, 4, 9, 10, 11 };
- morphology m{sample_tree(ms, parents), false};
+ auto tree = segments_from_points(points, parents);
+ morphology m{tree};
// Including all samples:
{
cable_cell cell(m);
- cv_policy pol = cv_policy_every_sample();
+ cv_policy pol = cv_policy_every_segment();
mlocation_list expected = {
{0, 0}, {0, 0.25}, {0, 0.5}, {0, 0.75}, {0, 1},
@@ -299,7 +301,7 @@ TEST(cv_policy, every_sample) {
{
cable_cell cell(m);
region reg = join(reg::branch(1), reg::branch(2));
- cv_policy pol = cv_policy_every_sample(reg);
+ cv_policy pol = cv_policy_every_segment(reg);
// Get samples from branches 1 and 2, plus boundary points from completions of each
// branch, viz. (0, 1), (2, 0), (1, 1) from branch 1 and (0, 1), (1, 0), (2, 1) from
@@ -327,7 +329,7 @@ TEST(cv_policy, domain) {
EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_fixed_per_branch(3, reg1, interior_forks).domain()));
EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3, reg1).domain()));
EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3, reg1, interior_forks).domain()));
- EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_every_sample(reg1).domain()));
+ EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_every_segment(reg1).domain()));
EXPECT_TRUE(region_eq(cell.provider(), join(reg1, reg2), (cv_policy_single(reg1)+cv_policy_single(reg2)).domain()));
EXPECT_TRUE(region_eq(cell.provider(), join(reg1, reg2), (cv_policy_single(reg1)|cv_policy_single(reg2)).domain()));
diff --git a/test/unit/test_event_delivery.cpp b/test/unit/test_event_delivery.cpp
index dea55b57179b2c92159a9d96190f3e672af19f47..d21fa4a5ec7f578bb35f33a14fadcc3b85bb36e7 100644
--- a/test/unit/test_event_delivery.cpp
+++ b/test/unit/test_event_delivery.cpp
@@ -6,10 +6,11 @@
// * Inject events one per cell in a given order, and confirm generated spikes
// are in the same order.
+#include <arbor/cable_cell.hpp>
#include <arbor/common_types.hpp>
#include <arbor/domain_decomposition.hpp>
+#include <arbor/morph/segment_tree.hpp>
#include <arbor/simulation.hpp>
-#include <arbor/cable_cell.hpp>
#include <arbor/spike.hpp>
#include <arbor/spike_event.hpp>
@@ -27,8 +28,8 @@ struct test_recipe: public n_cable_cell_recipe {
explicit test_recipe(int n): n_cable_cell_recipe(n, test_cell()) {}
static cable_cell test_cell() {
- sample_tree st;
- st.append({0,0,0,10,1});
+ segment_tree st;
+ st.append(mnpos, {0,0, 0,10}, {0,0,20,10}, 1);
label_dict d;
d.set("soma", arb::reg::tagged(1));
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index 4d4be869b9f96e9273367e4975e64b358d803f7c..dc5f9527e8cd3764f632a0dd4abae9672e9929e0 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -7,6 +7,8 @@
#include <arbor/mechcat.hpp>
#include <arbor/util/optional.hpp>
+#include "arbor/morph/morphology.hpp"
+#include "arbor/morph/segment_tree.hpp"
#include "fvm_layout.hpp"
#include "util/maputil.hpp"
#include "util/rangeutil.hpp"
@@ -26,18 +28,35 @@ using util::count_along;
using util::value_by_key;
namespace {
- std::vector<cable_cell> two_cell_system() {
+ struct system {
+ std::vector<soma_cell_builder> builders;
std::vector<cable_cell> cells;
+ };
+
+ system two_cell_system() {
+ system s;
+ auto& cells = s.cells;
+
+ // Cell 0: simple ball and stick
+ {
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
- // Cell 0: simple ball and stick (see common_cells.hpp)
- cells.push_back(make_cell_ball_and_stick());
+ cells.push_back(builder.make_cell());
+ auto& cell = cells.back();
+ cell.paint("soma", "hh");
+ cell.paint("dend", "pas");
+ cell.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
+
+ s.builders.push_back(std::move(builder));
+ }
// Cell 1: ball and 3-stick, but with uneven dendrite
// length and heterogeneous electrical properties:
//
// Bulk resistivity: 90 Ω·cm
// capacitance:
- // soma: 0.01 F/m² [default]
+ // soma: 0.01 F/m² [default]
// branch 1: 0.017 F/m²
// branch 2: 0.013 F/m²
// branch 3: 0.018 F/m²
@@ -60,46 +79,54 @@ namespace {
//
// All dendrite branches with 4 compartments.
- soma_cell_builder builder(7.);
- auto b1 = builder.add_branch(0, 200, 0.5, 0.5, 4, "dend");
- auto b2 = builder.add_branch(1, 300, 0.4, 0.4, 4, "dend");
- auto b3 = builder.add_branch(1, 180, 0.35, 0.35, 4, "dend");
- auto cell = builder.make_cell();
+ {
+ soma_cell_builder b(7.);
+ auto b1 = b.add_branch(0, 200, 0.5, 0.5, 4, "dend");
+ auto b2 = b.add_branch(1, 300, 0.4, 0.4, 4, "dend");
+ auto b3 = b.add_branch(1, 180, 0.35, 0.35, 4, "dend");
+ cells.push_back(b.make_cell());
+ auto& cell = cells.back();
- cell.paint("soma", "hh");
- cell.paint("dend", "pas");
+ cell.paint("soma", "hh");
+ cell.paint("dend", "pas");
- using ::arb::reg::branch;
- cell.paint(branch(b1), membrane_capacitance{0.017});
- cell.paint(branch(b2), membrane_capacitance{0.013});
- cell.paint(branch(b3), membrane_capacitance{0.018});
+ using ::arb::reg::branch;
+ auto c1 = reg::cable(b1-1, b.location({b1, 0}).pos, 1);
+ auto c2 = reg::cable(b2-1, b.location({b2, 0}).pos, 1);
+ auto c3 = reg::cable(b3-1, b.location({b3, 0}).pos, 1);
+ cell.paint(c1, membrane_capacitance{0.017});
+ cell.paint(c2, membrane_capacitance{0.013});
+ cell.paint(c3, membrane_capacitance{0.018});
- cell.place(mlocation{2,1}, i_clamp{5., 80., 0.45});
- cell.place(mlocation{3,1}, i_clamp{40., 10.,-0.2});
+ cell.place(b.location({2,1}), i_clamp{5., 80., 0.45});
+ cell.place(b.location({3,1}), i_clamp{40., 10.,-0.2});
- cell.default_parameters.axial_resistivity = 90;
+ cell.default_parameters.axial_resistivity = 90;
- cells.push_back(std::move(cell));
+ s.builders.push_back(std::move(b));
+ }
- return cells;
+ return s;
}
void check_two_cell_system(std::vector<cable_cell>& cells) {
ASSERT_EQ(2u, cells.size());
- ASSERT_EQ(2u, cells[0].morphology().num_branches());
- ASSERT_EQ(4u, cells[1].morphology().num_branches());
+ ASSERT_EQ(1u, cells[0].morphology().num_branches());
+ ASSERT_EQ(3u, cells[1].morphology().num_branches());
}
} // namespace
TEST(fvm_layout, mech_index) {
- std::vector<cable_cell> cells = two_cell_system();
+ auto system = two_cell_system();
+ auto& cells = system.cells;
+ auto& builders = system.builders;
check_two_cell_system(cells);
// Add four synapses of two varieties across the cells.
- cells[0].place(mlocation{1, 0.4}, "expsyn");
- cells[0].place(mlocation{1, 0.4}, "expsyn");
- cells[1].place(mlocation{2, 0.4}, "exp2syn");
- cells[1].place(mlocation{3, 0.4}, "expsyn");
+ cells[0].place(builders[0].location({1, 0.4}), "expsyn");
+ cells[0].place(builders[0].location({1, 0.4}), "expsyn");
+ cells[1].place(builders[1].location({2, 0.4}), "exp2syn");
+ cells[1].place(builders[1].location({3, 0.4}), "expsyn");
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
@@ -208,13 +235,16 @@ TEST(fvm_layout, coalescing_synapses) {
using L=std::initializer_list<unsigned>;
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
+
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
- cell.place(mlocation{1, 0.3}, "expsyn");
- cell.place(mlocation{1, 0.5}, "expsyn");
- cell.place(mlocation{1, 0.7}, "expsyn");
- cell.place(mlocation{1, 0.9}, "expsyn");
+ cell.place(builder.location({1, 0.3}), "expsyn");
+ cell.place(builder.location({1, 0.5}), "expsyn");
+ cell.place(builder.location({1, 0.7}), "expsyn");
+ cell.place(builder.location({1, 0.9}), "expsyn");
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
@@ -224,13 +254,13 @@ TEST(fvm_layout, coalescing_synapses) {
EXPECT_EQ(ivec({1, 1, 1, 1}), expsyn_config.multiplicity);
}
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
// Add synapses of two varieties.
- cell.place(mlocation{1, 0.3}, "expsyn");
- cell.place(mlocation{1, 0.5}, "exp2syn");
- cell.place(mlocation{1, 0.7}, "expsyn");
- cell.place(mlocation{1, 0.9}, "exp2syn");
+ cell.place(builder.location({1, 0.3}), "expsyn");
+ cell.place(builder.location({1, 0.5}), "exp2syn");
+ cell.place(builder.location({1, 0.7}), "expsyn");
+ cell.place(builder.location({1, 0.9}), "exp2syn");
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
@@ -244,12 +274,12 @@ TEST(fvm_layout, coalescing_synapses) {
EXPECT_EQ(ivec({1, 1}), exp2syn_config.multiplicity);
}
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
- cell.place(mlocation{1, 0.3}, "expsyn");
- cell.place(mlocation{1, 0.5}, "expsyn");
- cell.place(mlocation{1, 0.7}, "expsyn");
- cell.place(mlocation{1, 0.9}, "expsyn");
+ cell.place(builder.location({1, 0.3}), "expsyn");
+ cell.place(builder.location({1, 0.5}), "expsyn");
+ cell.place(builder.location({1, 0.7}), "expsyn");
+ cell.place(builder.location({1, 0.9}), "expsyn");
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
@@ -259,13 +289,13 @@ TEST(fvm_layout, coalescing_synapses) {
EXPECT_TRUE(expsyn_config.multiplicity.empty());
}
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
// Add synapses of two varieties.
- cell.place(mlocation{1, 0.3}, "expsyn");
- cell.place(mlocation{1, 0.5}, "exp2syn");
- cell.place(mlocation{1, 0.7}, "expsyn");
- cell.place(mlocation{1, 0.9}, "exp2syn");
+ cell.place(builder.location({1, 0.3}), "expsyn");
+ cell.place(builder.location({1, 0.5}), "exp2syn");
+ cell.place(builder.location({1, 0.7}), "expsyn");
+ cell.place(builder.location({1, 0.9}), "exp2syn");
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
@@ -279,13 +309,13 @@ TEST(fvm_layout, coalescing_synapses) {
EXPECT_TRUE(exp2syn_config.multiplicity.empty());
}
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
// Add synapses of two varieties.
- cell.place(mlocation{1, 0.3}, "expsyn");
- cell.place(mlocation{1, 0.3}, "expsyn");
- cell.place(mlocation{1, 0.7}, "expsyn");
- cell.place(mlocation{1, 0.7}, "expsyn");
+ cell.place(builder.location({1, 0.3}), "expsyn");
+ cell.place(builder.location({1, 0.3}), "expsyn");
+ cell.place(builder.location({1, 0.7}), "expsyn");
+ cell.place(builder.location({1, 0.7}), "expsyn");
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
@@ -295,13 +325,13 @@ TEST(fvm_layout, coalescing_synapses) {
EXPECT_EQ(ivec({2, 2}), expsyn_config.multiplicity);
}
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
// Add synapses of two varieties.
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 0, 0.2));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 0, 0.2));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 0.1, 0.2));
- cell.place(mlocation{1, 0.7}, syn_desc("expsyn", 0.1, 0.2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 0.2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 0.2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 0.1, 0.2));
+ cell.place(builder.location({1, 0.7}), syn_desc("expsyn", 0.1, 0.2));
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
@@ -317,17 +347,17 @@ TEST(fvm_layout, coalescing_synapses) {
}
}
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
// Add synapses of two varieties.
- cell.place(mlocation{1, 0.7}, syn_desc("expsyn", 0, 3));
- cell.place(mlocation{1, 0.7}, syn_desc("expsyn", 1, 3));
- cell.place(mlocation{1, 0.7}, syn_desc("expsyn", 0, 3));
- cell.place(mlocation{1, 0.7}, syn_desc("expsyn", 1, 3));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 0, 2));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 1, 2));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 0, 2));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 1, 2));
+ cell.place(builder.location({1, 0.7}), syn_desc("expsyn", 0, 3));
+ cell.place(builder.location({1, 0.7}), syn_desc("expsyn", 1, 3));
+ cell.place(builder.location({1, 0.7}), syn_desc("expsyn", 0, 3));
+ cell.place(builder.location({1, 0.7}), syn_desc("expsyn", 1, 3));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 1, 2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 1, 2));
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
@@ -344,19 +374,19 @@ TEST(fvm_layout, coalescing_synapses) {
}
}
{
- cable_cell cell = make_cell_ball_and_stick();
+ auto cell = builder.make_cell();
// Add synapses of two varieties.
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 1, 2));
- cell.place(mlocation{1, 0.3}, syn_desc_2("exp2syn", 4, 1));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 1, 2));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 5, 1));
- cell.place(mlocation{1, 0.3}, syn_desc_2("exp2syn", 1, 3));
- cell.place(mlocation{1, 0.3}, syn_desc("expsyn", 1, 2));
- cell.place(mlocation{1, 0.7}, syn_desc_2("exp2syn", 2, 2));
- cell.place(mlocation{1, 0.7}, syn_desc_2("exp2syn", 2, 1));
- cell.place(mlocation{1, 0.7}, syn_desc_2("exp2syn", 2, 1));
- cell.place(mlocation{1, 0.7}, syn_desc_2("exp2syn", 2, 2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 1, 2));
+ cell.place(builder.location({1, 0.3}), syn_desc_2("exp2syn", 4, 1));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 1, 2));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 5, 1));
+ cell.place(builder.location({1, 0.3}), syn_desc_2("exp2syn", 1, 3));
+ cell.place(builder.location({1, 0.3}), syn_desc("expsyn", 1, 2));
+ cell.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 2));
+ cell.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 1));
+ cell.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 1));
+ cell.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 2));
fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
@@ -376,7 +406,9 @@ TEST(fvm_layout, coalescing_synapses) {
}
TEST(fvm_layout, synapse_targets) {
- std::vector<cable_cell> cells = two_cell_system();
+ auto system = two_cell_system();
+ auto& cells = system.cells;
+ auto& builders = system.builders;
// Add synapses with different parameter values so that we can
// ensure: 1) CVs for each synapse mechanism are sorted while
@@ -393,14 +425,14 @@ TEST(fvm_layout, synapse_targets) {
return mechanism_desc(name).set("e", syn_e.at(idx));
};
- cells[0].place(mlocation{1, 0.9}, syn_desc("expsyn", 0));
- cells[0].place(mlocation{0, 0.5}, syn_desc("expsyn", 1));
- cells[0].place(mlocation{1, 0.4}, syn_desc("expsyn", 2));
+ cells[0].place(builders[0].location({1, 0.9}), syn_desc("expsyn", 0));
+ cells[0].place(builders[0].location({0, 0.5}), syn_desc("expsyn", 1));
+ cells[0].place(builders[0].location({1, 0.4}), syn_desc("expsyn", 2));
- cells[1].place(mlocation{2, 0.4}, syn_desc("exp2syn", 3));
- cells[1].place(mlocation{1, 0.4}, syn_desc("exp2syn", 4));
- cells[1].place(mlocation{3, 0.4}, syn_desc("expsyn", 5));
- cells[1].place(mlocation{3, 0.7}, syn_desc("exp2syn", 6));
+ cells[1].place(builders[1].location({2, 0.4}), syn_desc("exp2syn", 3));
+ cells[1].place(builders[1].location({1, 0.4}), syn_desc("exp2syn", 4));
+ cells[1].place(builders[1].location({3, 0.4}), syn_desc("expsyn", 5));
+ cells[1].place(builders[1].location({3, 0.7}), syn_desc("exp2syn", 6));
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
@@ -522,11 +554,11 @@ TEST(fvm_layout, density_norm_area) {
std::vector<double> expected_gl(ncv, dflt_gl);
// Last 1/6 of branch 1
- double seg1_area_right = cells[0].embedding().integrate_area(mcable{1, 5./6., 1.});
+ double seg1_area_right = cells[0].embedding().integrate_area(builder.cable({1, 5./6., 1.}));
// First 1/6 of branch 2
- double seg2_area_left = cells[0].embedding().integrate_area(mcable{2, 0., 1./6.});
+ double seg2_area_left = cells[0].embedding().integrate_area(builder.cable({2, 0., 1./6.}));
// First 1/6 of branch 3
- double seg3_area_left = cells[0].embedding().integrate_area(mcable{3, 0., 1./6.});
+ double seg3_area_left = cells[0].embedding().integrate_area(builder.cable({3, 0., 1./6.}));
// CV 0: soma
// CV1: left of branch 1
@@ -610,19 +642,19 @@ TEST(fvm_layout, density_norm_area_partial) {
auto cell = builder.make_cell();
cell.default_parameters.discretization = cv_policy_fixed_per_branch(1);
- cell.paint(mcable{1, 0., 0.3}, hh_begin);
- cell.paint(mcable{1, 0.4, 1.}, hh_end);
+ cell.paint(builder.cable({1, 0., 0.3}), hh_begin);
+ cell.paint(builder.cable({1, 0.4, 1.}), hh_end);
std::vector<cable_cell> cells{std::move(cell)};
+ // Area of whole cell (which is area of the 1 branch)
+ double area = cells[0].embedding().integrate_area({0, 0., 1});
// First 30% of branch 1.
- double b1_area_begin = cells[0].embedding().integrate_area(mcable{1, 0., 0.3});
+ double b1_area_begin = cells[0].embedding().integrate_area(builder.cable({1, 0., 0.3}));
// Last 60% of branch 1.
- double b1_area_end = cells[0].embedding().integrate_area(mcable{1, 0.4, 1.});
- // All of branch 1.
- double b1_area = cells[0].embedding().integrate_area(mcable{1, 0., 1.});
+ double b1_area_end = cells[0].embedding().integrate_area(builder.cable({1, 0.4, 1.}));
- double expected_norm_area = (b1_area_begin+b1_area_end)/b1_area;
+ double expected_norm_area = (b1_area_begin+b1_area_end)/area;
double expected_gnabar = dflt_gnabar;
double expected_gkbar = (dflt_gkbar*b1_area_begin + end_gkbar*b1_area_end)/(b1_area_begin + b1_area_end);
double expected_gl = (dflt_gl*b1_area_begin + end_gl*b1_area_end)/(b1_area_begin + b1_area_end);
@@ -699,29 +731,29 @@ TEST(fvm_layout, ion_weights) {
// 1/2 of branch 1 and the initial 1/2 of branches 2 and 3.
//
// Geometry:
- // soma 0: radius 5 µm
- // dend 1: 100 µm long, 1 µm diameter cylinder, tag 2
- // dend 2: 200 µm long, 1 µm diameter cylinder, tag 3
- // dend 3: 100 µm long, 1 µm diameter cylinder, tag 4
+ // soma 0: radius 5 µm, area 100π μm²
+ // dend 1: 100 µm long, 1 µm diameter cylinder, area 100π μm²
+ // dend 2: 200 µm long, 1 µm diameter cylinder, area 200π μm²
+ // dend 3: 100 µm long, 1 µm diameter cylinder, area 100π μm²
//
// The radius of the soma is chosen such that the surface area of soma is
// the same as a 100µm dendrite, which makes it easier to describe the
// expected weights.
- auto construct_cell = []() {
- soma_cell_builder builder(5);
- builder.add_branch(0, 100, 0.5, 0.5, 1, "dend");
- builder.add_branch(1, 200, 0.5, 0.5, 1, "dend");
- builder.add_branch(1, 100, 0.5, 0.5, 1, "dend");
- return builder.make_cell();
- };
+ soma_cell_builder builder(5);
+ builder.add_branch(0, 100, 0.5, 0.5, 1, "dend");
+ builder.add_branch(1, 200, 0.5, 0.5, 1, "dend");
+ builder.add_branch(1, 100, 0.5, 0.5, 1, "dend");
using uvec = std::vector<fvm_size_type>;
using ivec = std::vector<fvm_index_type>;
using fvec = std::vector<fvm_value_type>;
+ //uvec mech_branches[] = {
+ //{0}, {0,2}, {2, 3}, {0, 1, 2, 3}, {3}
+ //};
uvec mech_branches[] = {
- {0}, {0,2}, {2, 3}, {0, 1, 2, 3}, {3}
+ {0}, {0,2}
};
ivec expected_ion_cv[] = {
@@ -748,10 +780,11 @@ TEST(fvm_layout, ion_weights) {
for (auto run: count_along(mech_branches)) {
SCOPED_TRACE("run "+std::to_string(run));
- auto c = construct_cell();
+ auto c = builder.make_cell();
for (auto i: mech_branches[run]) {
- c.paint(reg::branch(i), "test_ca");
+ auto cab = builder.cable({i, 0, 1});
+ c.paint(reg::cable(cab.branch, cab.prox_pos, cab.dist_pos), "test_ca");
}
std::vector<cable_cell> cells{std::move(c)};
@@ -847,8 +880,10 @@ TEST(fvm_layout, vinterp_cable) {
// should interpolate between that and the parent.
// Cable cell with just one branch, non-spherical root.
- morphology morph(sample_tree({msample{0., 0., 0., 1.}, msample{10., 0., 0., 1.}}, {mnpos, 0u}));
- cable_cell cell(morph);
+ arb::segment_tree tree;
+ tree.append(mnpos, { 0,0,0,1}, {10,0,0,1}, 1);
+ arb::morphology m(tree);
+ cable_cell cell(m);
// CV midpoints at branch pos 0.1, 0.3, 0.5, 0.7, 0.9.
// Expect voltage reference locations to be CV modpoints.
@@ -901,10 +936,12 @@ TEST(fvm_layout, vinterp_forked) {
// interpolation.
// Cable cell with three branchses; branches 0 has child branches 1 and 2.
- morphology morph(sample_tree(
- {{0., 0., 0., 1.}, {10., 0., 0., 1}, {10., 20., 0., 1}, {10., -20., 0., 1}},
- {mnpos, 0u, 1u, 1u}));
- cable_cell cell(morph);
+ segment_tree tree;
+ tree.append(mnpos, {0., 0., 0., 1.}, {10., 0., 0., 1}, 1);
+ tree.append( 0, {10., 20., 0., 1}, 1);
+ tree.append( 0, {10.,-20., 0., 1}, 1);
+ morphology m(tree);
+ cable_cell cell(m);
// CV 0 contains branch 0 and the fork point; CV 1 and CV 2 have CV 0 as parent,
// and contain branches 1 and 2 respectively, excluding the fork point.
@@ -1003,10 +1040,12 @@ TEST(fvm_layout, iinterp) {
// Use the same cell/discretiazation as in vinterp_forked test:
// Cable cell with three branchses; branches 0 has child branches 1 and 2.
- morphology morph(sample_tree(
- {{0., 0., 0., 1.}, {10., 0., 0., 1}, {10., 20., 0., 1}, {10., -20., 0., 1}},
- {mnpos, 0u, 1u, 1u}));
- cable_cell cell(morph);
+ segment_tree tree;
+ tree.append(mnpos, {0., 0., 0., 1.}, {10., 0., 0., 1}, 1);
+ tree.append( 0, {10., 20., 0., 1}, 1);
+ tree.append( 0, {10.,-20., 0., 1}, 1);
+ morphology m(tree);
+ cable_cell cell(m);
// CV 0 contains branch 0 and the fork point; CV 1 and CV 2 have CV 0 as parent,
// and contain branches 1 and 2 respectively, excluding the fork point.
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 4110db15708690fc5c7ca3bcafba487a96501ff5..1771eda6085b34217a4a19f9b3922286cf2a854e 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -266,12 +266,12 @@ TEST(fvm_lowered, target_handles) {
make_cell_ball_and_3stick()
};
- EXPECT_EQ(cells[0].morphology().num_branches(), 2u);
- EXPECT_EQ(cells[1].morphology().num_branches(), 4u);
+ EXPECT_EQ(cells[0].morphology().num_branches(), 1u);
+ EXPECT_EQ(cells[1].morphology().num_branches(), 3u);
// (in increasing target order)
- cells[0].place(mlocation{1, 0.4}, "expsyn");
- cells[0].place(mlocation{0, 0.5}, "expsyn");
+ cells[0].place(mlocation{0, 0.7}, "expsyn");
+ cells[0].place(mlocation{0, 0.3}, "expsyn");
cells[1].place(mlocation{2, 0.2}, "exp2syn");
cells[1].place(mlocation{2, 0.8}, "expsyn");
@@ -282,9 +282,9 @@ TEST(fvm_lowered, target_handles) {
probe_association_map probe_map;
auto test_target_handles = [&](fvm_cell& cell) {
- mechanism *expsyn = find_mechanism(cell, "expsyn");
+ mechanism* expsyn = find_mechanism(cell, "expsyn");
ASSERT_TRUE(expsyn);
- mechanism *exp2syn = find_mechanism(cell, "exp2syn");
+ mechanism* exp2syn = find_mechanism(cell, "exp2syn");
ASSERT_TRUE(exp2syn);
unsigned expsyn_id = expsyn->mechanism_id();
@@ -341,8 +341,10 @@ TEST(fvm_lowered, stimulus) {
std::vector<cable_cell> cells;
cells.push_back(make_cell_ball_and_stick(false));
- cells[0].place(mlocation{1,1}, i_clamp{5., 80., 0.3});
- cells[0].place(mlocation{0,0.5}, i_clamp{1., 2., 0.1});
+ // At end of stick
+ cells[0].place(mlocation{0,1}, i_clamp{5., 80., 0.3});
+ // On the soma CV, which is over the approximate interval: (cable 0 0 0.1)
+ cells[0].place(mlocation{0,0.05}, i_clamp{1., 2., 0.1});
const fvm_size_type soma_cv = 0u;
const fvm_size_type tip_cv = 5u;
@@ -427,9 +429,9 @@ TEST(fvm_lowered, derived_mechs) {
std::vector<cable_cell> cells;
cells.reserve(3);
+ soma_cell_builder builder(6);
+ builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
for (int i = 0; i<3; ++i) {
- soma_cell_builder builder(6);
- builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
auto cell = builder.make_cell();
switch (i) {
@@ -451,7 +453,7 @@ TEST(fvm_lowered, derived_mechs) {
rec.catalogue() = make_unit_test_catalogue();
rec.catalogue().derive("custom_kin1", "test_kin1", {{"tau", 20.0}});
- cable_probe_total_ion_current_density where{mlocation{1, 0.3}};
+ cable_probe_total_ion_current_density where{builder.location({1, 0.3})};
rec.add_probe(0, 0, where);
rec.add_probe(1, 0, where);
rec.add_probe(2, 0, where);
@@ -758,8 +760,8 @@ TEST(fvm_lowered, point_ionic_current) {
// Test area-weighted linear combination of ion species concentrations
TEST(fvm_lowered, weighted_write_ion) {
- // Create a cell with 4 branches (same morphopology as in fvm_layout.ion_weights test):
- // - Soma (branch 0) plus three dendrites (1, 2, 3) meeting at a branch point.
+ // Create a cell with 3 branches (same morphopology as in fvm_layout.ion_weights test):
+ // - Soma (part of branch 0) plus three dendrites (d1, d2, d3) meeting at a branch point.
// - Dendritic segments are given 1 compartments each.
//
// /
@@ -773,10 +775,10 @@ TEST(fvm_lowered, weighted_write_ion) {
// 1/2 of branch 1 and the initial 1/2 of branches 2 and 3.
//
// Geometry:
- // soma 0: radius 5 µm
- // dend 1: 100 µm long, 1 µm diameter cylinder
- // dend 2: 200 µm long, 1 µm diameter cylinder
- // dend 3: 100 µm long, 1 µm diameter cylinder
+ // soma 0: 10 µm long, 10 µm diameter cylinder: area = 100π μm²
+ // dend 1: 100 µm long, 1 µm diameter cylinder: area = 100π μm²
+ // dend 2: 200 µm long, 1 µm diameter cylinder: area = 200π μm²
+ // dend 3: 100 µm long, 1 µm diameter cylinder: area = 100π μm²
//
// The radius of the soma is chosen such that the surface area of soma is
// the same as a 100 µm dendrite, which makes it easier to describe the
@@ -802,10 +804,10 @@ TEST(fvm_lowered, weighted_write_ion) {
const double con_ext = 120;
// Ca ion reader test_kinlva on CV 2 and 3 via branch 2:
- c.paint(reg::branch(2), "test_kinlva");
+ c.paint(reg::branch(1), "test_kinlva");
// Ca ion writer test_ca on CV 2 and 4 via branch 3:
- c.paint(reg::branch(3), "test_ca");
+ c.paint(reg::branch(2), "test_ca");
cable1d_recipe rec(c);
rec.catalogue() = make_unit_test_catalogue();
@@ -897,7 +899,7 @@ TEST(fvm_lowered, gj_coords_simple) {
soma_cell_builder b(2.1);
b.add_branch(0, 10, 0.3, 0.2, 5, "dend");
auto c = b.make_cell();
- c.place(mlocation{1, 0.8}, gap_junction_site{});
+ c.place(b.location({1, 0.8}), gap_junction_site{});
cells.push_back(std::move(c));
}
@@ -905,7 +907,7 @@ TEST(fvm_lowered, gj_coords_simple) {
soma_cell_builder b(2.4);
b.add_branch(0, 10, 0.3, 0.2, 2, "dend");
auto c = b.make_cell();
- c.place(mlocation{1, 1}, gap_junction_site{});
+ c.place(b.location({1, 1}), gap_junction_site{});
cells.push_back(std::move(c));
}
@@ -980,7 +982,7 @@ TEST(fvm_lowered, gj_coords_complex) {
b0.add_branch(0, 8, 0.3, 0.2, 4, "dend");
auto c0 = b0.make_cell();
- mlocation c0_gj[2] = {{1, 1}, {1, 0.5}};
+ mlocation c0_gj[2] = {b0.location({1, 1}), b0.location({1, 0.5})};
c0.place(c0_gj[0], gap_junction_site{});
c0.place(c0_gj[1], gap_junction_site{});
@@ -991,7 +993,7 @@ TEST(fvm_lowered, gj_coords_complex) {
b1.add_branch(1, 5, 0.2, 0.2, 5, "dend");
auto c1 = b1.make_cell();
- mlocation c1_gj[4] = {{2, 1}, {1, 1}, {1, 0.45}, {1, 0.1}};
+ mlocation c1_gj[4] = {b1.location({2, 1}), b1.location({1, 1}), b1.location({1, 0.45}), b1.location({1, 0.1})};
c1.place(c1_gj[0], gap_junction_site{});
c1.place(c1_gj[1], gap_junction_site{});
@@ -1007,7 +1009,7 @@ TEST(fvm_lowered, gj_coords_complex) {
b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
auto c2 = b2.make_cell();
- mlocation c2_gj[3] = {{1, 0.5}, {4, 1}, {2, 1}};
+ mlocation c2_gj[3] = {b2.location({1, 0.5}), b2.location({4, 1}), b2.location({2, 1})};
c2.place(c2_gj[0], gap_junction_site{});
c2.place(c2_gj[1], gap_junction_site{});
diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_mc_cell_group.cpp
index 1be195302555b8cf917c8bc1c78031c92f489d42..2971974752102f622f80b9cca10aa7e3f4bc0d89 100644
--- a/test/unit/test_mc_cell_group.cpp
+++ b/test/unit/test_mc_cell_group.cpp
@@ -26,8 +26,8 @@ namespace {
cable_cell c = builder.make_cell();
c.paint("soma", "hh");
c.paint("dend", "pas");
- c.place(mlocation{1,1}, i_clamp{5, 80, 0.3});
- c.place(mlocation{0, 0}, threshold_detector{0});
+ c.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
+ c.place(builder.location({0, 0}), threshold_detector{0});
return c;
}
}
@@ -66,7 +66,7 @@ TEST(mc_cell_group, sources) {
for (int i=0; i<20; ++i) {
cells.push_back(make_cell());
if (i==0 || i==3 || i==17) {
- cells.back().place(mlocation{1, 0.3}, threshold_detector{2.3});
+ cells.back().place(mlocation{0, 0.3}, threshold_detector{2.3});
}
EXPECT_EQ(1u + (i==0 || i==3 || i==17), cells.back().detectors().size());
diff --git a/test/unit/test_mc_cell_group_gpu.cpp b/test/unit/test_mc_cell_group_gpu.cpp
index e92e8926b3a217858893a51cc4817b134acfeabb..5a6d11cccfa414dab265b74d8e4123db15978daf 100644
--- a/test/unit/test_mc_cell_group_gpu.cpp
+++ b/test/unit/test_mc_cell_group_gpu.cpp
@@ -27,8 +27,8 @@ namespace {
cable_cell c = builder.make_cell();
c.paint("soma", "hh");
c.paint("dend", "pas");
- c.place(mlocation{1,1}, i_clamp{5, 80, 0.3});
- c.place(mlocation{0, 0}, threshold_detector{0});
+ c.place(builder.location({1, 1}), i_clamp{5, 80, 0.3});
+ c.place(builder.location({0, 0}), threshold_detector{0});
return c;
}
}
diff --git a/test/unit/test_morph_embedding.cpp b/test/unit/test_morph_embedding.cpp
index b079c602331c93a0354c54c9f83ae625ecce2b95..f8de6f4a8048b0fbfb1b525cc90fb1374a26c9a2 100644
--- a/test/unit/test_morph_embedding.cpp
+++ b/test/unit/test_morph_embedding.cpp
@@ -1,16 +1,18 @@
#include <cmath>
+#include <unordered_map>
#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 <arbor/morph/segment_tree.hpp>
#include "util/piecewise.hpp"
#include "../test/gtest.h"
#include "common.hpp"
+#include "common_cells.hpp"
using namespace arb;
using embedding = embed_pwlin;
@@ -34,9 +36,9 @@ using embedding = embed_pwlin;
}
}
-TEST(embedding, samples_and_branch_length) {
+TEST(embedding, segments_and_branch_length) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
using loc = mlocation;
// A single unbranched cable with 5 sample points.
@@ -44,23 +46,26 @@ TEST(embedding, samples_and_branch_length) {
// 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},
+ svec points = {
+ { 0, 0, 0, 2},
+ { 1, 0, 0, 2},
+ { 3, 0, 0, 2},
+ { 7, 0, 0, 2},
+ {10, 0, 0, 2},
};
- sample_tree sm(samples, parents);
- morphology m(sm, false);
+ morphology m(segments_from_points(points, parents));
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})));
+ auto nloc = 5u;
+ EXPECT_EQ(nloc, em.segment_locations().size());
+ auto locs = arb::util::make_range(em.branch_segment_locations(0));
+ EXPECT_EQ(nloc, locs.size());
+ EXPECT_TRUE(location_eq(m, locs[0], (loc{0,0})));
+ EXPECT_TRUE(location_eq(m, locs[1], (loc{0,0.1})));
+ EXPECT_TRUE(location_eq(m, locs[2], (loc{0,0.3})));
+ EXPECT_TRUE(location_eq(m, locs[3], (loc{0,0.7})));
+ EXPECT_TRUE(location_eq(m, locs[4], (loc{0,1})));
EXPECT_EQ(10., em.branch_length(0));
}
@@ -73,93 +78,66 @@ TEST(embedding, samples_and_branch_length) {
// 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);
+ pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
+
+ svec points = {
+ { 0, 0, 0, 2},
+ { 10, 0, 0, 2},
+ {100, 0, 0, 2},
+ { 0, 10, 0, 2},
+ { 0,100, 0, 2},
+ {100,100, 0, 2},
+ { 0,130, 0, 2},
+ { 0,300, 0, 2},
+ };
+ morphology m(segments_from_points(points, parents));
- 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());
+ ASSERT_EQ(4u, m.num_branches());
- embedding em(m);
+ 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));
- }
+ auto locs = arb::util::make_range(em.branch_segment_locations(0));
+ EXPECT_TRUE(location_eq(m, locs[0], (loc{0,0})));
+ EXPECT_TRUE(location_eq(m, locs[1], (loc{0,0.1})));
+ EXPECT_TRUE(location_eq(m, locs[2], (loc{0,1})));
+
+ locs = em.branch_segment_locations(1);
+ EXPECT_TRUE(location_eq(m, locs[0], (loc{1,0})));
+ EXPECT_TRUE(location_eq(m, locs[1], (loc{1,0.1})));
+ EXPECT_TRUE(location_eq(m, locs[2], (loc{1,1})));
+
+ locs = em.branch_segment_locations(2);
+ EXPECT_TRUE(location_eq(m, locs[0], (loc{2,0})));
+ EXPECT_TRUE(location_eq(m, locs[1], (loc{2,1})));
+
+ locs = em.branch_segment_locations(3);
+ EXPECT_TRUE(location_eq(m, locs[0], (loc{3,0})));
+ EXPECT_TRUE(location_eq(m, locs[1], (loc{3,0.15})));
+ EXPECT_TRUE(location_eq(m, locs[2], (loc{3,1})));
+
+ 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 svec = std::vector<mpoint>;
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}
+ svec points = {
+ { 0, 0, 0, 10},
+ {10, 0, 0, 20},
+ {30, 0, 0, 10},
+ {30, 10, 0, 5},
+ {30, 0, 50, 5}
};
- morphology m(sample_tree(samples, parents), false);
+ morphology m(segments_from_points(points, parents));
embedding em(m);
EXPECT_DOUBLE_EQ(30., em.branch_length(0));
@@ -183,20 +161,20 @@ TEST(embedding, partial_branch_length) {
TEST(embedding, partial_area) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
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}
+ svec points = {
+ { 0, 0, 0, 10},
+ {10, 0, 0, 20},
+ {30, 0, 0, 10},
+ {30, 10, 0, 5},
+ {30, 0, 50, 5}
};
- morphology m(sample_tree(samples, parents), false);
+ morphology m(segments_from_points(points, parents));
embedding em(m);
// Cable 1: single truncated cone, length L = 10,
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
index 40307aa0b8630b3266ce70afc84c1edd195e3535..fd5ddbafa14f9b92eaeab240ee1bc39711bca54f 100644
--- a/test/unit/test_morph_expr.cpp
+++ b/test/unit/test_morph_expr.cpp
@@ -9,11 +9,12 @@
#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"
+#include "common.hpp"
+#include "common_cells.hpp"
#include "morph_pred.hpp"
using namespace arb;
@@ -66,16 +67,13 @@ TEST(locset, expr_repn) {
auto root = ls::root();
auto term = ls::terminal();
- auto samp = ls::sample(42);
auto loc = ls::location(2, 0.5);
auto bdy = ls::boundary(reg::tagged(1));
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(sum(root, term, loc)), "(sum (sum (root) (terminal)) (location 2 0.5))");
EXPECT_EQ(to_string(loc), "(location 2 0.5)");
EXPECT_EQ(to_string(bdy), "(boundary (tag 1))");
}
@@ -94,18 +92,18 @@ TEST(locset, invalid_mlocation) {
TEST(locset, thingify_named) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
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});
+ auto sm = segments_from_points(svec{ {0,0,0,1}, {10,0,0,1} }, pvec{mnpos, 0});
{
label_dict dict;
dict.set("banana", banana);
dict.set("cake", cake);
- mprovider mp(morphology(sm, false), dict);
+ mprovider mp(morphology(sm), dict);
EXPECT_EQ(thingify(locset("cake"), mp), thingify(cake, mp));
EXPECT_EQ(thingify(locset("banana"), mp), thingify(banana, mp));
@@ -118,7 +116,7 @@ TEST(locset, thingify_named) {
dict.set("topping", locset("fruit"));
dict.set("fruit", locset("strawberry"));
- EXPECT_THROW(mprovider(morphology(sm, false), dict), unbound_name);
+ EXPECT_THROW(mprovider(morphology(sm), dict), unbound_name);
}
{
label_dict dict;
@@ -127,26 +125,26 @@ TEST(locset, thingify_named) {
dict.set("topping", locset("fruit"));
dict.set("fruit", sum(locset("banana"), locset("topping")));
- EXPECT_THROW(mprovider(morphology(sm, false), dict), circular_definition);
+ EXPECT_THROW(mprovider(morphology(sm), dict), circular_definition);
}
}
TEST(region, thingify_named) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
region banana = reg::branch(0);
region cake = reg::cable(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});
+ auto sm = segments_from_points(svec{ {0,0,0,1}, {10,0,0,1} }, pvec{mnpos, 0});
{
label_dict dict;
dict.set("banana", banana);
dict.set("cake", cake);
- mprovider mp(morphology(sm, false), dict);
+ mprovider mp(morphology(sm), dict);
EXPECT_EQ(thingify(region("cake"), mp), thingify(cake, mp));
EXPECT_EQ(thingify(region("banana"), mp), thingify(banana, mp));
@@ -159,7 +157,7 @@ TEST(region, thingify_named) {
dict.set("topping", region("fruit"));
dict.set("fruit", region("strawberry"));
- EXPECT_THROW(mprovider(morphology(sm, false), dict), unbound_name);
+ EXPECT_THROW(mprovider(morphology(sm), dict), unbound_name);
}
{
label_dict dict;
@@ -168,7 +166,7 @@ TEST(region, thingify_named) {
dict.set("topping", region("fruit"));
dict.set("fruit", join(region("cake"), region("topping")));
- EXPECT_THROW(mprovider(morphology(sm, false), dict), circular_definition);
+ EXPECT_THROW(mprovider(morphology(sm), dict), circular_definition);
}
}
@@ -176,11 +174,10 @@ TEST(region, thingify_named) {
TEST(locset, thingify) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
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);
@@ -198,43 +195,23 @@ TEST(locset, thingify) {
// 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},
+ svec points = {
+ { 0, 0, 0, 2},
+ { 10, 0, 0, 2},
+ {100, 0, 0, 2},
+ { 0, 10, 0, 2},
+ { 0,100, 0, 2},
+ {100,100, 0, 2},
+ { 0,200, 0, 2},
+ { 0,300, 0, 2},
};
- sample_tree sm(samples, parents);
+ auto sm = segments_from_points(points, 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}}));
-
- // Check round-trip of implicit locset conversions.
- // (Use a locset which is non-trivially a multiset in order to
- // test the fold in the constructor.)
- EXPECT_EQ(thingify(multi, mp), thingify(locset(thingify(multi, mp)), mp));
- }
- {
- mprovider mp(morphology(sm, false));
+ auto mp = mprovider(morphology(sm));
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}}));
@@ -246,7 +223,7 @@ TEST(locset, thingify) {
EXPECT_THROW(thingify(begb4, mp), no_such_branch);
}
{
- mprovider mp(morphology(sm, false));
+ auto mp = mprovider(morphology(sm));
auto all = reg::all();
auto ls0 = thingify(ls::uniform(all, 0, 9, 12), mp);
@@ -301,7 +278,7 @@ TEST(locset, thingify) {
EXPECT_TRUE(found == 2);
}
{
- mprovider mp(morphology(sm, false));
+ auto mp = mprovider(morphology(sm));
auto sub_reg = join(reg::cable(0, 0.2, 0.7), reg::cable(1, 0.1, 1), reg::cable(3, 0.5, 0.6));
auto ls0 = thingify(ls::uniform(sub_reg, 0, 10000, 72), mp);
@@ -331,7 +308,8 @@ TEST(locset, thingify) {
TEST(region, thingify_simple_morphologies) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
+ using tvec = std::vector<int>;
using cl = mcable_list;
// A single unbranched cable with 5 sample points.
@@ -339,15 +317,17 @@ TEST(region, thingify_simple_morphologies) {
// 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},
+ svec points = {
+ { 0, 0, 0, 2},
+ { 1, 0, 0, 2},
+ { 3, 0, 0, 2},
+ { 7, 0, 0, 2},
+ {10, 0, 0, 2},
};
- sample_tree sm(samples, parents);
- mprovider mp(morphology(sm, false));
+ tvec tags = {1, 1, 2, 1, 2};
+
+ auto sm = segments_from_points(points, parents, tags);
+ auto mp = mprovider(morphology(sm));
auto h1 = reg::cable(0, 0, 0.5);
auto h2 = reg::cable(0, 0.5, 1);
@@ -398,25 +378,28 @@ TEST(region, thingify_simple_morphologies) {
//
// sample ids:
// 0 |
- // 1 3 |
- // 2 4 |
+ // 1 |
+ // 2 4 |
+ // 3 5 |
// tags:
// 1 |
+ // 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},
+ pvec parents = {mnpos, 0, 1, 2, 1, 4};
+ svec points = {
+ { 0, 0, 0, 5},
+ { 10, 0, 0, 5},
+ { 10, 0, 0, 2},
+ {100, 0, 0, 2},
+ { 10, 0, 0, 2},
+ {100, 0, 0, 2},
};
+ tvec tags = {1, 1, 3, 3, 2, 2};
- // Build morphology with a spherical root
- sample_tree sm(samples, parents);
- mprovider mp(morphology(sm, true));
+ auto sm = segments_from_points(points, parents, tags);
+ auto mp = mprovider(morphology(sm));
using ls::location;
using reg::tagged;
@@ -431,9 +414,9 @@ TEST(region, thingify_simple_morphologies) {
locset end1_ = location(1,1);
auto reg0_ = distal_interval(start1_, 45);
- auto reg1_ = distal_interval(mid0_, 74);
+ auto reg1_ = distal_interval(mid0_, 50);
auto reg2_ = proximal_interval(end1_, 45);
- auto reg3_ = proximal_interval(end1_, 91);
+ auto reg3_ = proximal_interval(end1_, 95);
auto reg4_ = distal_interval(end1_, 100);
auto reg5_ = distal_interval(start1_, 0);
auto reg6_ = proximal_interval(start1_, 0);
@@ -444,9 +427,9 @@ TEST(region, thingify_simple_morphologies) {
EXPECT_TRUE(region_eq(mp, join(tagged(1), tagged(2), tagged(3)), mcable_list{{0,0,1}, {1,0,1}, {2,0,1}}));
EXPECT_TRUE(region_eq(mp, join(tagged(1), tagged(2), tagged(3)), all()));
EXPECT_TRUE(region_eq(mp, reg0_, mcable_list{{1,0,0.5}}));
- EXPECT_TRUE(region_eq(mp, reg1_, mcable_list{{0,0.5,1}, {1,0,0.8}, {2,0,0.8}}));
+ EXPECT_TRUE(region_eq(mp, reg1_, mcable_list{{0,0.5,1}, {1,0,0.5}, {2,0,0.5}}));
EXPECT_TRUE(region_eq(mp, reg2_, mcable_list{{1,0.5,1}}));
- EXPECT_TRUE(region_eq(mp, reg3_, mcable_list{{0, 0.75, 1}, {1,0,1}}));
+ EXPECT_TRUE(region_eq(mp, reg3_, mcable_list{{0, 0.5, 1}, {1,0,1}}));
EXPECT_TRUE(region_eq(mp, reg4_, mcable_list{{1,1,1}}));
EXPECT_TRUE(region_eq(mp, reg5_, mcable_list{{1,0,0}}));
EXPECT_TRUE(region_eq(mp, reg6_, mcable_list{{1,0,0}}));
@@ -455,7 +438,8 @@ TEST(region, thingify_simple_morphologies) {
TEST(region, thingify_moderate_morphologies) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
+ using tvec = std::vector<int>;
using cl = mcable_list;
// Test multi-level morphologies.
@@ -476,20 +460,21 @@ TEST(region, thingify_moderate_morphologies) {
// 3
{
pvec parents = {mnpos, 0, 1, 0, 3, 4, 4, 6};
- svec samples = {
- {{ 0, 0, 0, 1}, 1},
- {{ 10, 0, 0, 1}, 3},
- {{100, 0, 0, 3}, 3},
- {{ 0, 10, 0, 1}, 2},
- {{ 0,100, 0, 5}, 2},
- {{100,100, 0, 2}, 4},
- {{ 0,200, 0, 1}, 3},
- {{ 0,300, 0, 3}, 3},
+ svec points = {
+ { 0, 0, 0, 1},
+ { 10, 0, 0, 1},
+ {100, 0, 0, 3},
+ { 0, 10, 0, 1},
+ { 0,100, 0, 5},
+ {100,100, 0, 2},
+ { 0,200, 0, 1},
+ { 0,300, 0, 3},
};
- sample_tree sm(samples, parents);
+ tvec tags = {1, 3, 3, 2, 2, 4, 3, 3};
+ auto sm = segments_from_points(points, parents, tags);
// Without spherical root
- mprovider mp(morphology(sm, false));
+ auto mp = mprovider(morphology(sm));
using ls::location;
using reg::tagged;
@@ -764,25 +749,26 @@ TEST(region, thingify_moderate_morphologies) {
TEST(region, thingify_complex_morphologies) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
{
pvec parents = {mnpos, 0, 1, 0, 3, 4, 5, 5, 7, 7, 4, 10};
- svec samples = {
- {{ 0, 0, 0, 2}, 3}, //0
- {{ 10, 0, 0, 2}, 3}, //1
- {{100, 0, 0, 2}, 3}, //2
- {{ 0, 10, 0, 2}, 3}, //3
- {{ 0,100, 0, 2}, 3}, //4
- {{100,100, 0, 2}, 3}, //5
- {{200,100, 0, 2}, 3}, //6
- {{100,200, 0, 2}, 3}, //7
- {{200,200, 0, 2}, 3}, //8
- {{100,300, 0, 2}, 3}, //9
- {{ 0,200, 0, 2}, 3}, //10
- {{ 0,300, 0, 2}, 3}, //11
+ svec points = {
+ { 0, 0, 0, 2}, //0
+ { 10, 0, 0, 2}, //1
+ {100, 0, 0, 2}, //2
+ { 0, 10, 0, 2}, //3
+ { 0,100, 0, 2}, //4
+ {100,100, 0, 2}, //5
+ {200,100, 0, 2}, //6
+ {100,200, 0, 2}, //7
+ {200,200, 0, 2}, //8
+ {100,300, 0, 2}, //9
+ { 0,200, 0, 2}, //10
+ { 0,300, 0, 2}, //11
};
- sample_tree sm(samples, parents);
- auto m = morphology(sm, false);
+
+ auto sm = segments_from_points(points, parents);
+ auto m = morphology(sm);
{
mprovider mp(m);
using reg::cable;
@@ -818,23 +804,23 @@ TEST(region, thingify_complex_morphologies) {
}
{
pvec parents = {mnpos, 0, 1, 1, 2, 3, 0, 6, 7, 8, 7};
- svec samples = {
- {{ 0, 10, 10, 1}, 1},
- {{ 0, 30, 30, 1}, 2},
- {{ 0, 60,-20, 1}, 2},
- {{ 0, 90, 70, 1}, 2},
- {{ 0, 80,-10, 1}, 2},
- {{ 0,100,-40, 1}, 2},
- {{ 0,-50,-50, 1}, 2},
- {{ 0, 20,-30, 2}, 2},
- {{ 0, 40,-80, 2}, 2},
- {{ 0,-30,-80, 3}, 2},
- {{ 0, 90,-70, 5}, 2}
+ svec points = {
+ {0, 10, 10, 1},
+ {0, 30, 30, 1},
+ {0, 60,-20, 1},
+ {0, 90, 70, 1},
+ {0, 80,-10, 1},
+ {0,100,-40, 1},
+ {0,-50,-50, 1},
+ {0, 20,-30, 2},
+ {0, 40,-80, 2},
+ {0,-30,-80, 3},
+ {0, 90,-70, 5},
};
- sample_tree sm(samples, parents);
+ auto sm = segments_from_points(points, parents);
// Without spherical root
- mprovider mp(morphology(sm, false));
+ auto mp = mprovider(morphology(sm));
using reg::all;
using reg::nil;
diff --git a/test/unit/test_morph_place.cpp b/test/unit/test_morph_place.cpp
index 7ba2308efc8a54301f58b375b3ceb3ca3c94f502..dc7bebb5eb3a743fc05519a9e1e75cbc68fe4490 100644
--- a/test/unit/test_morph_place.cpp
+++ b/test/unit/test_morph_place.cpp
@@ -5,12 +5,12 @@
#include <arbor/morph/place_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"
+#include "common_cells.hpp"
using namespace arb;
@@ -183,21 +183,21 @@ TEST(isometry, compose) {
TEST(place_pwlin, cable) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
// L-shaped simple cable.
// 0.25 of the length in the z-direction,
// 0.75 of the length in the x-direction.
pvec parents = {mnpos, 0, 1};
- svec samples = {
- {{ 0, 0, 0, 2}, 1},
- {{ 0, 0, 1, 2}, 1},
- {{ 3, 0, 1, 2}, 1}
+ svec points = {
+ {0, 0, 0, 2},
+ {0, 0, 1, 2},
+ {3, 0, 1, 2}
};
- sample_tree sm(samples, parents);
- morphology m(sm, false);
+ auto sm = segments_from_points(points, parents);
+ morphology m(sm);
{
// With no transformation:
@@ -251,28 +251,28 @@ TEST(place_pwlin, cable) {
TEST(place_pwlin, branched) {
using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
+ using svec = std::vector<mpoint>;
// Y-shaped branched morphology.
// Second branch (branch 1) tapers radius.
pvec parents = {mnpos, 0, 1, 2, 3, 4, 2, 6};
- svec samples = {
+ svec points = {
// branch 0
- {{ 0, 0, 0, 2}, 1},
- {{ 0, 0, 1, 2}, 1},
- {{ 3, 0, 1, 2}, 1},
+ { 0, 0, 0, 2},
+ { 0, 0, 1, 2},
+ { 3, 0, 1, 2},
// branch 1
- {{ 3, 0, 1, 1.0}, 1},
- {{ 3, 1, 1, 1.0}, 1},
- {{ 3, 2, 1, 0.0}, 1},
+ { 3, 0, 1, 1.0},
+ { 3, 1, 1, 1.0},
+ { 3, 2, 1, 0.0},
// branch 2
- {{ 3, 0, 1, 2} , 1},
- {{ 3, -1, 1, 2} , 1}
+ { 3, 0, 1, 2},
+ { 3, -1, 1, 2}
};
- sample_tree sm(samples, parents);
- morphology m(sm, false);
+ auto sm = segments_from_points(points, parents);
+ morphology m(sm);
isometry iso = isometry::translate(2, 3, 4);
place_pwlin pl(m, iso);
diff --git a/test/unit/test_morph_primitives.cpp b/test/unit/test_morph_primitives.cpp
index c627e589b9fe31bb10c7f25f2b6eea7012ef3a6b..c077cdde39e78ffa3ecfaf914664accd913c92fe 100644
--- a/test/unit/test_morph_primitives.cpp
+++ b/test/unit/test_morph_primitives.cpp
@@ -34,12 +34,6 @@ TEST(morph_primitives, is_collocated) {
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) {
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
index 48d987488bac96473b15ec5fb95f09c3b47d7fa1..5cc929b611b6837b0156fa19c85a8876a2b76858 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -7,10 +7,10 @@
#include <arbor/morph/morphexcept.hpp>
#include <arbor/morph/morphology.hpp>
-#include <arbor/morph/sample_tree.hpp>
#include <arbor/cable_cell.hpp>
-#include "morph/mbranch.hpp"
+#include "arbor/morph/primitives.hpp"
+#include "arbor/morph/segment_tree.hpp"
#include "util/span.hpp"
#include "morph_pred.hpp"
@@ -36,348 +36,111 @@ TEST(morphology, mpoint) {
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;
-
- EXPECT_FALSE(arb::is_terminal(p));
- EXPECT_FALSE(arb::is_fork(p));
- EXPECT_FALSE(arb::is_root(p));
- EXPECT_FALSE(arb::is_collocated(p));
-
- arb::set_root(p);
- EXPECT_FALSE(arb::is_terminal(p));
- EXPECT_FALSE(arb::is_fork(p));
- EXPECT_TRUE(arb::is_root(p));
- EXPECT_FALSE(arb::is_collocated(p));
-
- arb::set_terminal(p);
- EXPECT_TRUE(arb::is_terminal(p));
- EXPECT_FALSE(arb::is_fork(p));
- EXPECT_TRUE(arb::is_root(p));
- EXPECT_FALSE(arb::is_collocated(p));
-
- arb::unset_root(p);
- EXPECT_TRUE(arb::is_terminal(p));
- EXPECT_FALSE(arb::is_fork(p));
- EXPECT_FALSE(arb::is_root(p));
- EXPECT_FALSE(arb::is_collocated(p));
-
- arb::set_collocated(p);
- EXPECT_TRUE(arb::is_terminal(p));
- EXPECT_FALSE(arb::is_fork(p));
- EXPECT_FALSE(arb::is_root(p));
- EXPECT_TRUE(arb::is_collocated(p));
-
- arb::set_fork(p);
- EXPECT_TRUE(arb::is_terminal(p));
- EXPECT_TRUE(arb::is_fork(p));
- EXPECT_FALSE(arb::is_root(p));
- EXPECT_TRUE(arb::is_collocated(p));
-
- arb::unset_fork(p);
- arb::unset_terminal(p);
- arb::unset_collocated(p);
- EXPECT_FALSE(arb::is_terminal(p));
- EXPECT_FALSE(arb::is_fork(p));
- EXPECT_FALSE(arb::is_root(p));
- EXPECT_FALSE(arb::is_collocated(p));
-}
-
-// TODO: test sample_tree marks properties correctly
-
-// Test internal function that parses a parent list, and marks
-// each node as either root, sequential, fork or terminal.
-TEST(sample_tree, properties) {
- const auto npos = arb::mnpos;
- using arb::sample_tree;
- using pp = arb::point_prop;
- using pvec = std::vector<arb::msize_t>;
-
- pp c = arb::point_prop_mask_collocated;
- pp r = arb::point_prop_mask_root;
- pp t = arb::point_prop_mask_terminal;
- pp s = arb::point_prop_mask_none;
- pp f = arb::point_prop_mask_fork;
- pp tc = t+c;
- pp sc = s+c;
- pp fc = f+c;
-
- // make a sample tree from a parent vector with non-collocated points.
- auto make_tree = [] (const pvec& parents) {
- sample_tree st;
- for (auto p: parents) st.append(p, {{0.,0.,double(st.size()),1.}, 1});
- return st;
- };
- // make a sample tree from a parent vector with collocated points.
- auto make_colloc_tree = [] (const pvec& parents) {
- sample_tree st;
- for (auto p: parents) st.append(p, {{0.,0.,0.,1.}, 1});
- return st;
- };
-
- {
- EXPECT_EQ(make_tree({npos}).properties(), std::vector<pp>{r});
- EXPECT_EQ(make_colloc_tree({npos}).properties(), std::vector<pp>{r});
- }
-
- {
- EXPECT_EQ(make_tree({npos,0}).properties(), (std::vector<pp>{r,t}));
- EXPECT_EQ(make_colloc_tree({npos,0}).properties(), (std::vector<pp>{r,tc}));
- }
-
- {
- EXPECT_EQ(make_tree({npos,0,1,2}).properties(), (std::vector<pp>{r,s,s,t}));
- EXPECT_EQ(make_colloc_tree({npos,0,1,2}).properties(), (std::vector<pp>{r,sc,sc,tc}));
- }
-
- {
- EXPECT_EQ(make_tree({npos,0,1,2,0,4,5}).properties(), (std::vector<pp>{r,s,s,t,s,s,t}));
- EXPECT_EQ(make_colloc_tree({npos,0,1,2,0,4,5}).properties(), (std::vector<pp>{r,sc,sc,tc,sc,sc,tc}));
- }
-
- {
- EXPECT_EQ(make_tree({npos,0,1,2,3,2,4,4,7}).properties(), (std::vector<pp>{r,s,f,s,f,t,t,s,t}));
- EXPECT_EQ(make_colloc_tree({npos,0,1,2,3,2,4,4,7}).properties(), (std::vector<pp>{r,sc,fc,sc,fc,tc,tc,sc,tc}));
- }
-}
-
-TEST(morphology, branches_from_parent_index) {
- const auto npos = arb::mnpos;
- using pvec = std::vector<arb::msize_t>;
- using mb = arb::impl::mbranch;
-
- // make a sample tree from a parent vector with non-collocated points.
- auto make_tree = [] (const pvec& parents) {
- arb::sample_tree st;
- for (auto p: parents) st.append(p, {{0.,0.,double(st.size()),1.}, 1});
- return st;
- };
-
- { // single sample: can only be used to build a morphology with one spherical branch
- pvec parents{npos};
- auto tree = make_tree(parents);
- auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
- EXPECT_EQ(1u, bc.size());
- EXPECT_EQ(mb({0}, npos), bc[0]);
-
- // A cable morphology can't be constructed from a single sample.
- EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), false),
- arb::incomplete_branch);
- }
- {
- pvec parents = {npos, 0};
- auto tree = make_tree(parents);
- auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
- EXPECT_EQ(1u, bc.size());
- EXPECT_EQ(mb({0,1}, npos), bc[0]);
-
- // 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::incomplete_branch);
- }
-
- {
- pvec parents{npos, 0, 1};
-
- // With cable soma: one cable with 3 samples.
- auto tree = make_tree(parents);
- auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
- EXPECT_EQ(1u, bc.size());
- EXPECT_EQ(mb({0,1,2},npos), bc[0]);
-
- // With spherical soma: one sphere and a 2-segment cable.
- // The cable branch is attached to the sphere (i.e. the sphere is the parent branch).
- auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
- EXPECT_EQ(2u, bs.size());
- EXPECT_EQ(mb({0},npos), bs[0]);
- EXPECT_EQ(mb({1,2},0), bs[1]);
- }
-
- {
- pvec parents{npos, 0, 0};
- 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::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);
- EXPECT_EQ(2u, bc.size());
- EXPECT_EQ(mb({0,1},npos), bc[0]);
- EXPECT_EQ(mb({0,2},npos), bc[1]);
- }
-
- {
- pvec parents{npos, 0, 1, 2};
- auto tree = make_tree(parents);
-
- // With cable soma: one cable with 4 samples.
- auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
- EXPECT_EQ(1u, bc.size());
- EXPECT_EQ(mb({0,1,2,3},npos), bc[0]);
-
- // With spherical soma: one sphere and on 3-segment cable.
- // The cable branch is attached to the sphere (i.e. the sphere is the parent branch).
- auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
- EXPECT_EQ(2u, bs.size());
- EXPECT_EQ(mb({0},npos), bs[0]);
- EXPECT_EQ(mb({1,2,3},0), bs[1]);
- }
-
- {
- pvec parents{npos, 0, 1, 0};
- auto tree = make_tree(parents);
-
- // With cable soma: two cables with 3 and 2 samples respectively.
- auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
- EXPECT_EQ(2u, bc.size());
- EXPECT_EQ(mb({0,1,2},npos), bc[0]);
- 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::incomplete_branch);
- }
-
- {
- pvec parents{npos, 0, 1, 0, 3};
- auto tree = make_tree(parents);
-
- // With cable soma: two cables with 3 samples each [0,1,2] and [0,3,4]
- auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
- EXPECT_EQ(2u, bc.size());
- EXPECT_EQ(mb({0,1,2},npos), bc[0]);
- EXPECT_EQ(mb({0,3,4},npos), bc[1]);
-
- // With spherical soma: one sphere and 2 2-sample cables.
- // each of the cable branches is attached to the sphere (i.e. the sphere is the parent branch).
- auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
- EXPECT_EQ(3u, bs.size());
- EXPECT_EQ(mb({0},npos), bs[0]);
- EXPECT_EQ(mb({1,2},0), bs[1]);
- EXPECT_EQ(mb({3,4},0), bs[2]);
- }
-
- {
- pvec parents{npos, 0, 1, 0, 3, 4, 4, 6};
- auto tree = make_tree(parents);
-
- // With cable soma: 4 cables: [0,1,2] [0,3,4] [4,5] [4,6,7]
- auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
- EXPECT_EQ(4u, bc.size());
- EXPECT_EQ(mb({0,1,2},npos), bc[0]);
- EXPECT_EQ(mb({0,3,4},npos), bc[1]);
- EXPECT_EQ(mb({4,5}, 1), bc[2]);
- EXPECT_EQ(mb({4,6,7},1), bc[3]);
-
- // With spherical soma: 1 sphere and 4 cables: [1,2] [3,4] [4,5] [4,6,7]
- auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
- EXPECT_EQ(5u, bs.size());
- EXPECT_EQ(mb({0}, npos), bs[0]);
- EXPECT_EQ(mb({1,2}, 0), bs[1]);
- EXPECT_EQ(mb({3,4}, 0), bs[2]);
- EXPECT_EQ(mb({4,5}, 2), bs[3]);
- EXPECT_EQ(mb({4,6,7}, 2), bs[4]);
- }
-}
-
// For different parent index vectors, attempt multiple valid and invalid sample sets.
TEST(morphology, construction) {
constexpr auto npos = arb::mnpos;
using arb::util::make_span;
- using ms = arb::msample;
+ using mp = arb::mpoint;
using pvec = std::vector<arb::msize_t>;
{
pvec p = {npos, 0};
- std::vector<ms> s = {
- {{0.0, 0.0, 0.0, 1.0}, 1},
- {{0.0, 0.0, 1.0, 1.0}, 1} };
+ std::vector<mp> s = {
+ {0.0, 0.0, 0.0, 1.0},
+ {0.0, 0.0, 1.0, 1.0}};
- arb::sample_tree sm(s, p);
- auto m = arb::morphology(sm);
+ arb::segment_tree tree;
+ tree.append(npos, s[0], s[1], 1);
+ auto m = arb::morphology(tree);
EXPECT_EQ(1u, m.num_branches());
}
{
pvec p = {npos, 0, 1};
- { // 2-segment cable (1 seg soma, 1 seg dendrite)
- std::vector<ms> s = {
- {{0.0, 0.0, 0.0, 5.0}, 1},
- {{0.0, 0.0, 5.0, 1.0}, 1},
- {{0.0, 0.0, 8.0, 1.0}, 2} };
-
- arb::sample_tree sm(s, p);
- auto m = arb::morphology(sm);
-
- EXPECT_EQ(1u, m.num_branches());
- }
- { // spherical soma and single-segment cable
- std::vector<ms> s = {
- {{0.0, 0.0, 0.0, 5.0}, 1},
- {{0.0, 0.0, 1.0, 1.0}, 2},
- {{0.0, 0.0, 8.0, 1.0}, 2} };
-
- arb::sample_tree sm(s, p);
- auto m = arb::morphology(sm);
+ // 2-segment cable
+ std::vector<mp> s = {
+ {0.0, 0.0, 0.0, 5.0},
+ {0.0, 0.0, 1.0, 1.0},
+ {0.0, 0.0, 8.0, 1.0} };
+
+ arb::segment_tree sm;
+ sm.append(npos, s[0], s[1], 2);
+ sm.append( 0, s[1], s[2], 2);
+ auto m = arb::morphology(sm);
- EXPECT_EQ(2u, m.num_branches());
- }
+ EXPECT_EQ(1u, m.num_branches());
}
{
// 0 |
// 1 3 |
// 2 |
pvec p = {npos, 0, 1, 0};
- {
- // two cables: 1x2 segments, 1x1 segment.
- std::vector<ms> s = {
- {{0.0, 0.0, 0.0, 5.0}, 1},
- {{0.0, 0.0, 5.0, 1.0}, 1},
- {{0.0, 0.0, 6.0, 1.0}, 2},
- {{0.0, 4.0, 0.0, 1.0}, 1}};
-
- arb::sample_tree sm(s, p);
- auto m = arb::morphology(sm);
-
- EXPECT_EQ(2u, m.num_branches());
- }
- {
- // error: spherical soma with a single point cable attached via sample 3
- std::vector<ms> s = {
- {{0.0, 0.0, 0.0, 5.0}, 1},
- {{0.0, 0.0, 5.0, 1.0}, 2},
- {{0.0, 0.0, 8.0, 1.0}, 2},
- {{0.0, 5.0, 0.0, 1.0}, 2}};
-
- arb::sample_tree sm(s, p);
- EXPECT_THROW((arb::morphology(sm)), arb::incomplete_branch);
- }
+ // two cables: 1x2 segments, 1x1 segment.
+ std::vector<mp> s = {
+ {0.0, 0.0, 0.0, 5.0},
+ {0.0, 0.0, 5.0, 1.0},
+ {0.0, 0.0, 6.0, 1.0},
+ {0.0, 4.0, 0.0, 1.0}};
+
+ arb::segment_tree tree;
+ tree.append(npos, s[0], s[1], 1);
+ tree.append( 0, s[1], s[2], 2);
+ tree.append(npos, s[0], s[3], 1);
+ auto m = arb::morphology(tree);
+
+ EXPECT_EQ(2u, m.num_branches());
}
{
// 0 |
// 1 3 |
// 2 4 |
- pvec p = {npos, 0, 1, 0, 3};
- {
- // two cables: 1x2 segments, 1x1 segment.
- std::vector<ms> s = {
- {{0.0, 0.0, 0.0, 5.0}, 1},
- {{0.0, 0.0, 5.0, 1.0}, 2},
- {{0.0, 0.0, 8.0, 1.0}, 2},
- {{0.0, 5.0, 0.0, 1.0}, 2},
- {{0.0, 8.0, 0.0, 1.0}, 2}};
-
- arb::sample_tree sm(s, p);
- auto m = arb::morphology(sm);
-
- EXPECT_EQ(3u, m.num_branches());
- }
+
+ // two cables: 2 segments each
+ std::vector<mp> s = {
+ {0.0, 0.0, 0.0, 5.0},
+ {0.0, 0.0, 5.0, 1.0},
+ {0.0, 0.0, 8.0, 1.0},
+ {0.0, 5.0, 0.0, 1.0},
+ {0.0, 8.0, 0.0, 1.0}};
+
+ arb::segment_tree tree;
+ tree.append(npos, s[0], s[1], 1);
+ tree.append( 0, s[1], s[2], 1);
+ tree.append(npos, s[0], s[3], 1);
+ tree.append( 2, s[3], s[4], 1);
+ auto m = arb::morphology(tree);
+
+ EXPECT_EQ(2u, m.num_branches());
+ }
+ {
+ // 0 |
+ // 1 3 |
+ // 2 4 5 |
+ pvec p = {npos, 0, 1, 0};
+ // 4 cables
+ std::vector<mp> s = {
+ {0.0, 0.0, 0.0, 5.0},
+ {0.0, 0.0, 5.0, 1.0},
+ {0.0, 0.0, 6.0, 1.0},
+ {0.0, 4.0, 0.0, 1.0},
+ {0.0, 4.0, 1.0, 1.0},
+ {0.0, 4.0, 2.0, 1.0}};
+
+ arb::segment_tree tree;
+ tree.append(npos, s[0], s[1], 1);
+ tree.append( 0, s[1], s[2], 2);
+ tree.append(npos, s[0], s[3], 1);
+ tree.append( 2, s[3], s[4], 1);
+ tree.append( 2, s[4], s[5], 1);
+ auto m = arb::morphology(tree);
+
+ EXPECT_EQ(4u, m.num_branches());
}
}
// test that morphology generates branch child-parent structure correctly.
TEST(morphology, branches) {
using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
+ using svec = std::vector<arb::mpoint>;
auto npos = arb::mnpos;
auto check_terminal_branches = [](const arb::morphology& m) {
@@ -391,13 +154,14 @@ TEST(morphology, branches) {
};
{
- // 0
- pvec parents = {npos};
- svec samples = {
- {{ 0,0,0,3}, 1},
+ // 0 - 1
+ svec p = {
+ {0.,0.,0.,3.},
+ {10.,0.,0.,3.},
};
- arb::sample_tree sm(samples, parents);
- arb::morphology m(sm);
+ arb::segment_tree tree;
+ tree.append(npos, p[0], p[1], 1);
+ arb::morphology m(tree);
EXPECT_EQ(1u, m.num_branches());
EXPECT_EQ(npos, m.branch_parent(0));
@@ -406,93 +170,97 @@ TEST(morphology, branches) {
check_terminal_branches(m);
}
{
- // 0 - 1
- pvec parents = {npos, 0};
- svec samples = {
- {{ 0,0,0,3}, 1},
- {{ 10,0,0,3}, 1},
+ // 0 - 1 - 2
+
+ // the morphology is a single unbranched cable.
+ svec p = {
+ {0.,0.,0.,3.},
+ {10.,0.,0.,3.},
+ {100,0,0,3},
};
- arb::sample_tree sm(samples, parents);
- arb::morphology m(sm);
+ arb::segment_tree tree;
+ tree.append(npos, p[0], p[1], 1);
+ tree.append( 0, p[1], p[2], 1);
+ arb::morphology m(tree);
EXPECT_EQ(1u, m.num_branches());
EXPECT_EQ(npos, m.branch_parent(0));
EXPECT_EQ(pvec{}, m.branch_children(0));
+ EXPECT_EQ(pvec{0}, m.terminal_branches());
check_terminal_branches(m);
}
{
- // 0 - 1 - 2
- pvec parents = {npos, 0, 1};
- {
- // All samples have same tag -> the morphology is a single unbranched cable.
- svec samples = {
- {{ 0,0,0,3}, 1},
- {{10,0,0,3}, 1},
- {{100,0,0,3}, 1},
- };
- arb::sample_tree sm(samples, parents);
- arb::morphology m(sm);
-
- 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.
- svec samples = {
- {{ 0,0,0,10}, 1},
- {{ 10,0,0, 3}, 3},
- {{100,0,0, 3}, 3},
- };
- arb::sample_tree sm(samples, parents);
- arb::morphology m(sm);
-
- EXPECT_EQ(2u, m.num_branches());
- EXPECT_EQ(npos, m.branch_parent(0));
- 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);
- }
+ // 6 segemnts and six branches.
+ // A single branch at the root that bifurcates, and a further 3 branches
+ // attached to the first bifurcation.
+ // 0 |
+ // / \ |
+ // 1 2 |
+ // / | \ |
+ // 3 4 5 |
+
+ svec p = {
+ {0., 0.,0.,3.},
+ {1., 0.,0.,3.},
+ {2., 1.,0.,3.},
+ {2.,-1.,0.,3.},
+ {3., 2.,0.,3.},
+ {3., 1.,0.,3.},
+ {3., 0.,0.,3.},
+ };
+ arb::segment_tree tree;
+ tree.append(npos, p[0], p[1], 1);
+ tree.append( 0, p[1], p[2], 1);
+ tree.append( 0, p[1], p[3], 1);
+ tree.append( 1, p[2], p[4], 1);
+ tree.append( 1, p[2], p[5], 1);
+ tree.append( 1, p[2], p[6], 1);
+ arb::morphology m(tree);
+
+ EXPECT_EQ(6u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ( 0u, m.branch_parent(1));
+ EXPECT_EQ( 0u, m.branch_parent(2));
+ EXPECT_EQ( 1u, m.branch_parent(3));
+ EXPECT_EQ( 1u, m.branch_parent(4));
+ EXPECT_EQ( 1u, m.branch_parent(5));
+ EXPECT_EQ((pvec{1,2}), m.branch_children(0));
+ EXPECT_EQ((pvec{3,4,5}), m.branch_children(1));
+ EXPECT_EQ((pvec{}), m.branch_children(2));
+ EXPECT_EQ((pvec{}), m.branch_children(3));
+ EXPECT_EQ((pvec{}), m.branch_children(4));
+ EXPECT_EQ((pvec{}), m.branch_children(5));
+ EXPECT_EQ((pvec{2,3,4,5}), m.terminal_branches());
+
+ check_terminal_branches(m);
}
{
- // 2 - 0 - 1
- pvec parents = {npos, 0, 0};
-
- svec samples = {
- {{ 0, 0,0, 5}, 3},
- {{ 10, 0,0, 5}, 3},
- {{ 0,10,0, 5}, 3},
+ // 0 |
+ // / \ |
+ // 1 2 |
+
+ svec p = {
+ { 0, 0,0, 5},
+ {10, 0,0, 5},
+ { 0,10,0, 5},
};
- arb::sample_tree sm(samples, parents);
- arb::morphology m(sm);
+ arb::segment_tree tree;
+ tree.append(npos, p[0], p[1], 3);
+ tree.append(npos, p[0], p[2], 3);
+ arb::morphology m(tree);
EXPECT_EQ(2u, m.num_branches());
EXPECT_EQ(npos, m.branch_parent(0));
EXPECT_EQ(npos, m.branch_parent(1));
EXPECT_EQ(pvec{}, m.branch_children(0));
EXPECT_EQ(pvec{}, m.branch_children(1));
+ EXPECT_EQ((pvec{0,1}), m.terminal_branches());
check_terminal_branches(m);
}
{
- // 0 - 1 - 2 - 3
- pvec parents = {npos, 0, 1, 2};
- }
- {
- // 0 |
- // / \ |
- // 1 3 |
- // / |
- // 2 |
- pvec parents = {npos, 0, 1, 0};
- }
- {
- // Eight samples
+ // 7 segments, 4 branches
//
// 0 |
// / \ |
@@ -503,61 +271,37 @@ TEST(morphology, branches) {
// 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);
- arb::morphology m(sm);
-
- EXPECT_EQ(5u, m.num_branches());
- EXPECT_EQ(npos, m.branch_parent(0));
- EXPECT_EQ(0u, m.branch_parent(1));
- EXPECT_EQ(0u, m.branch_parent(2));
- EXPECT_EQ(2u, m.branch_parent(3));
- EXPECT_EQ(2u, m.branch_parent(4));
- EXPECT_EQ((pvec{1,2}), m.branch_children(0));
- EXPECT_EQ((pvec{}), m.branch_children(1));
- 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 = {
- {{ 0, 0, 0, 10}, 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);
-
- EXPECT_EQ(4u, m.num_branches());
- EXPECT_EQ(npos, m.branch_parent(0));
- EXPECT_EQ(npos, m.branch_parent(1));
- EXPECT_EQ(1u, m.branch_parent(2));
- EXPECT_EQ(1u, m.branch_parent(3));
- EXPECT_EQ((pvec{}), m.branch_children(0));
- 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);
- }
+ svec p = {
+ { 0, 0, 0, 10},
+ { 10, 0, 0, 2},
+ {100, 0, 0, 2},
+ { 0, 10, 0, 2},
+ { 0,100, 0, 2},
+ {100,100, 0, 2},
+ { 0,200, 0, 2},
+ { 0,300, 0, 2},
+ };
+ arb::segment_tree tree;
+ tree.append(npos, p[0], p[1], 1);
+ tree.append( 0, p[1], p[2], 1);
+ tree.append(npos, p[0], p[3], 1);
+ tree.append( 2, p[3], p[4], 1);
+ tree.append( 3, p[4], p[5], 1);
+ tree.append( 3, p[4], p[6], 1);
+ tree.append( 5, p[6], p[7], 1);
+ arb::morphology m(tree);
+
+ EXPECT_EQ(4u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(npos, m.branch_parent(1));
+ EXPECT_EQ(1u, m.branch_parent(2));
+ EXPECT_EQ(1u, m.branch_parent(3));
+ EXPECT_EQ((pvec{}), m.branch_children(0));
+ 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);
}
}
@@ -575,116 +319,79 @@ TEST(morphology, swc) {
// Load swc samples from file.
auto swc_samples = arb::parse_swc_file(fid);
- // Build a sample_tree from swc samples.
- auto sm = arb::swc_as_sample_tree(swc_samples);
- EXPECT_EQ(1058u, sm.size()); // file contains 195 samples
+ // Build a segmewnt_tree from swc samples.
+ auto sm = arb::swc_as_segment_tree(swc_samples);
+ EXPECT_EQ(1057u, sm.size()); // file contains 195 samples
// Test that the morphology contains the expected number of branches.
auto m = arb::morphology(sm);
- EXPECT_EQ(31u, m.num_branches());
+ EXPECT_EQ(30u, m.num_branches());
}
#endif
-TEST(morphology, minset) {
- using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
- using ll = arb::mlocation_list;
+arb::morphology make_4_branch_morph() {
+ using svec = std::vector<arb::mpoint>;
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},
+ // Eight points, 7 segments, 3 branches.
+ // sample branch
+ // 0 0
+ // 1 3 0 1
+ // 2 4 0 1
+ // 5 6 2 3
+ // 7 3
+ svec p = {
+ { 0, 0, 0, 2},
+ { 10, 0, 0, 2},
+ {100, 0, 0, 2},
+ { 0, 10, 0, 2},
+ { 0,100, 0, 2},
+ {100,100, 0, 2},
+ { 0,200, 0, 2},
+ { 0,300, 0, 2},
};
- 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}}));
- }
+ arb::segment_tree tree;
+ tree.append(npos, p[0], p[1], 1);
+ tree.append( 0, p[1], p[2], 1);
+ tree.append(npos, p[0], p[3], 1);
+ tree.append( 2, p[3], p[4], 1);
+ tree.append( 3, p[4], p[5], 1);
+ tree.append( 3, p[4], p[6], 1);
+ tree.append( 4, p[6], p[7], 1);
+
+ return arb::morphology(tree);
+}
+
+TEST(morphology, minset) {
+ auto m = make_4_branch_morph();
+
+ using ll = arb::mlocation_list;
+
+ 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}}));
}
TEST(morphology, maxset) {
- using pvec = std::vector<arb::msize_t>;
- using svec = std::vector<arb::msample>;
+ auto m = make_4_branch_morph();
+
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{}), maxset(m, ll{}));
- EXPECT_EQ((ll{{2,0.1}}), maxset(m, ll{{2,0.1}}));
- EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), maxset(m, ll{{0,0.5}, {1,0.5}}));
- EXPECT_EQ((ll{{0,0.5}}), maxset(m, ll{{0,0.5}}));
- EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
- EXPECT_EQ((ll{{0,0.5}, {2,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0.5}, {2,0.5}}));
- EXPECT_EQ((ll{{0,0.5}, {2,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {2,0.5}}));
- EXPECT_EQ((ll{{0,0.5}, {2,0.5}, {3,1}}), maxset(m, ll{{0,0}, {0,0.5}, {2,0.5}, {3,0}, {3,1}}));
- EXPECT_EQ((ll{{2,0.5}, {3,0.5}}), maxset(m, ll{{1,0.5}, {2,0.5}, {3,0}, {3,0.2}, {3,0.5}}));
- }
- {
- arb::morphology m(sm, true);
-
- EXPECT_EQ((ll{}), maxset(m, ll{}));
- EXPECT_EQ((ll{{2,0.1}}), maxset(m, ll{{2,0.1}}));
- EXPECT_EQ((ll{{1,0.5}}), maxset(m, ll{{0,0.5}, {1,0.5}}));
- EXPECT_EQ((ll{{1,0}, {2,0}}), maxset(m, ll{{2,0}, {1,0}, {0,0}}));
- EXPECT_EQ((ll{{1,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
- EXPECT_EQ((ll{{1,1}, {3,0.1}, {4,0.7}}), maxset(m, ll{{1,0.5}, {1,1}, {3,0.1}, {4,0.5}, {4,0.7}}));
- }
+ EXPECT_EQ((ll{}), maxset(m, ll{}));
+ EXPECT_EQ((ll{{2,0.1}}), maxset(m, ll{{2,0.1}}));
+ EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), maxset(m, ll{{0,0.5}, {1,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}}), maxset(m, ll{{0,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {1,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0}, {1,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {2,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {1,0.5}, {2,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {2,0.5}}), maxset(m, ll{{0,0}, {0,0.5}, {2,0.5}}));
+ EXPECT_EQ((ll{{0,0.5}, {2,0.5}, {3,1}}), maxset(m, ll{{0,0}, {0,0.5}, {2,0.5}, {3,0}, {3,1}}));
+ EXPECT_EQ((ll{{2,0.5}, {3,0.5}}), maxset(m, ll{{1,0.5}, {2,0.5}, {3,0}, {3,0.2}, {3,0.5}}));
}
// Testing mextent; intersection and join operations are
@@ -694,32 +401,9 @@ TEST(mextent, invariants) {
using namespace arb;
using testing::cablelist_eq;
- using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
- using cl = mcable_list;
+ auto m = make_4_branch_morph();
- // Eight samples
- // no-sphere
- // sample branch
- // 0 0
- // 1 3 0 1
- // 2 4 0 1
- // 5 6 2 3
- // 7 3
- 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},
- };
-
- // Instantiate morphology with spherical_root=false.
- morphology m(sample_tree(samples, parents), false);
+ using cl = mcable_list;
mextent x1(cl{{1, 0.25, 0.75}});
ASSERT_TRUE(x1.test_invariants(m));
@@ -749,31 +433,9 @@ TEST(mextent, intersects) {
using namespace arb;
using testing::cablelist_eq;
- using pvec = std::vector<msize_t>;
- using svec = std::vector<msample>;
- using cl = mcable_list;
-
- // Eight samples
- // no-sphere
- // sample branch
- // 0 0
- // 1 3 0 1
- // 2 4 0 1
- // 5 6 2 3
- // 7 3
- 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},
- };
+ auto m = make_4_branch_morph();
- morphology m(sample_tree(samples, parents));
+ using cl = mcable_list;
mextent x1(cl{{1, 0.25, 0.75}});
EXPECT_TRUE(x1.intersects(mlocation{1, 0.25}));
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index ee17642978dd83596772b83ff5a8016884742f14..0d6a48c2826a09d2a99cdb2858968925da8db4ec 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -81,12 +81,17 @@ struct backend_access<gpu::backend> {
// linearly tapered branches.
static morphology make_y_morphology() {
- return morphology(sample_tree(
- {msample{{0., 0., 0., 1.}, 0},
- msample{{100., 0., 0., 0.8}, 0},
- msample{{100., 100., 0., 0.5}, 0},
- msample{{100., 0, 100., 0.4}, 0}},
- {mnpos, 0u, 1u, 1u}));
+ segment_tree tree;
+ tree.append(mnpos, {0., 0., 0., 1.}, {100., 0., 0., 0.8}, 1);
+ tree.append(0, {100., 100., 0., 0.5}, 0);
+ tree.append(0, {100., 0, 100., 0.4}, 0);
+ return {tree};
+}
+
+static morphology make_stick_morphology() {
+ segment_tree tree;
+ tree.append(mnpos, {0., 0., 0., 1.}, {100., 0., 0., 1.0}, 1);
+ return {tree};
}
template <typename Backend>
@@ -94,7 +99,11 @@ void run_v_i_probe_test(const context& ctx) {
using fvm_cell = typename backend_access<Backend>::fvm_cell;
auto deref = [](const fvm_value_type* p) { return backend_access<Backend>::deref(p); };
- cable_cell bs = make_cell_ball_and_stick(false);
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+ cable_cell bs = builder.make_cell();
+
bs.default_parameters.discretization = cv_policy_fixed_per_branch(1);
i_clamp stim(0, 100, 0.3);
@@ -256,7 +265,10 @@ void run_expsyn_g_probe_test(const context& ctx) {
mlocation loc0{1, 0.8};
mlocation loc1{1, 1.0};
- cable_cell bs = make_cell_ball_and_stick(false);
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+ cable_cell bs = builder.make_cell();
bs.place(loc0, "expsyn");
bs.place(loc1, "expsyn");
bs.default_parameters.discretization = cv_policy_fixed_per_branch(2);
@@ -485,7 +497,7 @@ void run_ion_density_probe_test(const context& ctx) {
// Simple constant diameter cable, 3 CVs.
- cable_cell cable(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ cable_cell cable(make_stick_morphology());
cable.default_parameters.discretization = cv_policy_fixed_per_branch(3);
// Calcium ions everywhere, half written by write_ca1, half by write_ca2.
@@ -647,7 +659,8 @@ void run_partial_density_probe_test(const context& ctx) {
// Each cell is a simple constant diameter cable, with 3 CVs each.
- morphology m(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 1.}, 0}}, {mnpos, 0u}));
+ auto m = make_stick_morphology();
+
cells[0] = cable_cell(m);
cells[0].default_parameters.discretization = cv_policy_fixed_per_branch(3);
@@ -760,7 +773,7 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
// Cell is a tapered cable with 3 CVs.
- morphology m(sample_tree({msample{{0., 0., 0., 1.}, 0}, msample{{100., 0., 0., 0.8}, 0}}, {mnpos, 0u}));
+ auto m = make_stick_morphology();
cable_cell cell(m);
const unsigned n_cv = 3;
@@ -947,7 +960,12 @@ void run_multi_probe_test(const context& ctx) {
template <typename Backend>
void run_v_sampled_probe_test(const context& ctx) {
- cable_cell bs = make_cell_ball_and_stick(false);
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+
+ cable_cell bs = builder.make_cell();
+
bs.default_parameters.discretization = cv_policy_fixed_per_branch(1);
std::vector<cable_cell> cells = {bs, bs};
@@ -1116,7 +1134,11 @@ void run_exact_sampling_probe_test(const context& ctx) {
adhoc_recipe() {
gprop_.default_parameters = neuron_parameter_defaults;
- cells_.assign(4, make_cell_ball_and_stick(false));
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
+
+ cells_.assign(4, builder.make_cell());
cells_[0].place(mlocation{1, 0.1}, "expsyn");
cells_[1].place(mlocation{1, 0.1}, "exp2syn");
cells_[2].place(mlocation{1, 0.9}, "expsyn");
diff --git a/test/unit/test_segment_tree.cpp b/test/unit/test_segment_tree.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4f5baff1d7aba5163dabbd0f2fd293b138e78098
--- /dev/null
+++ b/test/unit/test_segment_tree.cpp
@@ -0,0 +1,113 @@
+#include <algorithm>
+#include <fstream>
+#include <cmath>
+#include <random>
+#include <string>
+#include <vector>
+
+#include "../test/gtest.h"
+
+#include <arbor/morph/morphexcept.hpp>
+
+#include "arbor/morph/primitives.hpp"
+#include "util/span.hpp"
+
+#include "morph_pred.hpp"
+
+// Test basic functions on properties of mpoint
+TEST(segment_tree, 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(segment_tree, empty) {
+ using mp = arb::mpoint;
+ arb::segment_tree tree;
+ EXPECT_TRUE(tree.empty());
+ tree.append(arb::mnpos, mp{0,0,0,1}, mp{0,0,1,1}, 1);
+ EXPECT_FALSE(tree.empty());
+}
+
+TEST(segment_tree, invalid_append) {
+ using mp = arb::mpoint;
+ using arb::mnpos;
+ arb::segment_tree tree;
+
+ // Test that appropriate exceptions are thrown when attempting to create
+ // a segment by implicitly extending from the root.
+ EXPECT_THROW(tree.append(mnpos, mp{0,0,0,1}, 1), arb::invalid_segment_parent);
+
+ tree.append(mnpos, mp{0,0,0,1}, mp{0,0,1,1}, 1);
+
+ EXPECT_THROW(tree.append(mnpos, mp{0,0,0,1}, 1), arb::invalid_segment_parent);
+
+ // Test that an exception is thrown when attempting to append to
+ // a segment that is not in the tree (parent>=num_segs)
+ EXPECT_THROW(tree.append(1, mp{0,0,0,1}, mp{0,0,1,1}, 1), arb::invalid_segment_parent);
+ EXPECT_THROW(tree.append(2, mp{0,0,0,1}, mp{0,0,1,1}, 1), arb::invalid_segment_parent);
+ EXPECT_THROW(tree.append(2, mp{0,0,1,1}, 1), arb::invalid_segment_parent);
+}
+
+// Generate some random morphologies of different sizes, and verify that
+// the correct tree is constructed.
+TEST(segment_tree, fuzz) {
+ using mp = arb::mpoint;
+ using arb::mnpos;
+
+ int nrun = 10;
+ int max_size = 1<<12;
+
+ std::mt19937 gen(0);
+
+ auto make_point = [](arb::msize_t i) {return mp{0,0,double(i),1};};
+
+ for (auto nseg=2; nseg<=max_size; nseg*=2) {
+ for (int run=0; run<nrun; ++run) {
+ arb::segment_tree tree;
+ std::vector<arb::msize_t> parents;
+
+ parents.reserve(nseg);
+ tree.reserve(nseg);
+ for (int i=0; i<nseg; ++i) {
+ std::uniform_int_distribution<int> distrib(-1, i-1);
+
+ // Draw random segment to attach to.
+ arb::msize_t p = distrib(gen);
+ parents.push_back(p);
+
+ // Attach segment.
+ // Use an implicit append every third time when possible.
+ if (p!=mnpos && i%3)
+ tree.append(p, make_point(i), i);
+ else
+ tree.append(p, make_point(p), make_point(i), i);
+
+ // Validate that the correct number of segments created.
+ EXPECT_EQ(tree.size(), std::size_t(i+1));
+ }
+ EXPECT_EQ(parents, tree.parents());
+ for (int i=0; i<nseg; ++i) {
+ arb::msegment seg = tree.segments()[i];
+ EXPECT_EQ(seg.prox, make_point(parents[i]));
+ EXPECT_EQ(seg.dist, make_point(i));
+ EXPECT_EQ(seg.tag, i);
+ }
+ }
+ }
+}
+