From 717f7cb50a8d1a1685868d6ba08e891e5447ce9c Mon Sep 17 00:00:00 2001
From: Sam Yates <halfflat@gmail.com>
Date: Tue, 28 Jan 2020 17:38:16 +0100
Subject: [PATCH] Fix error in pw_ratpoly building in embed_pwlin. (#945)
* Fix '==' vs '=' error in assignment.
* Use last element (index 2) not index 1 of parent area and ixa tile to start from on new branch.
* Expose cross-branch consistency via point-to-point integration functions for branch length and area.
* Extend unit tests to suit.
---
arbor/include/arbor/morph/embed_pwlin.hpp | 3 +++
arbor/morph/embed_pwlin.cpp | 18 ++++++++++++++++--
test/unit/test_morph_embedding.cpp | 15 ++++++++++++++-
3 files changed, 33 insertions(+), 3 deletions(-)
diff --git a/arbor/include/arbor/morph/embed_pwlin.hpp b/arbor/include/arbor/morph/embed_pwlin.hpp
index d6947ac8..0338186b 100644
--- a/arbor/include/arbor/morph/embed_pwlin.hpp
+++ b/arbor/include/arbor/morph/embed_pwlin.hpp
@@ -35,10 +35,12 @@ struct embed_pwlin {
// Computed length of mcable.
double integrate_length(mcable c) const;
+ double integrate_length(mlocation proxmal, mlocation distal) const;
double integrate_length(msize_t bid, const pw_constant_fn&) const;
// Membrane surface area of given mcable.
double integrate_area(mcable c) const;
+ double integrate_area(mlocation proxmal, mlocation distal) const;
double integrate_area(msize_t bid, const pw_constant_fn&) const;
// Integrated inverse cross-sectional area of given mcable.
@@ -50,6 +52,7 @@ struct embed_pwlin {
return integrate_length(mcable{bid, 0, 1});
}
+
private:
std::vector<mlocation> sample_locations_;
std::shared_ptr<embed_pwlin_data> data_;
diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp
index d81e9e4c..eaf430a6 100644
--- a/arbor/morph/embed_pwlin.cpp
+++ b/arbor/morph/embed_pwlin.cpp
@@ -135,6 +135,20 @@ double embed_pwlin::integrate_ixa(mcable c) const {
return integrate_ixa(c.branch, pw_constant_fn{{c.prox_pos, c.dist_pos}, {1.}});
}
+// Point to point integration:
+
+double embed_pwlin::integrate_length(mlocation proximal, mlocation distal) const {
+ return interpolate(data_->length, distal.branch, distal.pos) -
+ interpolate(data_->length, proximal.branch, proximal.pos);
+}
+
+double embed_pwlin::integrate_area(mlocation proximal, mlocation distal) const {
+ return interpolate(data_->area, distal.branch, distal.pos) -
+ interpolate(data_->area, proximal.branch, proximal.pos);
+}
+
+// Subregions defined by geometric inequalities:
+
mcable_list embed_pwlin::radius_cmp(msize_t bid, double val, comp_op op) const {
switch (op) {
case comp_op::lt: return data_cmp(data_->radius, bid, val, [](auto l, auto r){return l < r;});
@@ -215,8 +229,8 @@ embed_pwlin::embed_pwlin(const arb::morphology& m) {
double length_0 = parent==mnpos? 0: data_->length[parent].back().second[1];
data_->length[bid].push_back(0., 1, rat_element<1, 0>(length_0, length_0+branch_length));
- double area_0 = parent==mnpos? 0: data_->area[parent].back().second[1];
- double ixa_0 = parent==mnpos? 0: data_->ixa[parent].back().second[1];
+ 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.
diff --git a/test/unit/test_morph_embedding.cpp b/test/unit/test_morph_embedding.cpp
index 0c72b62f..b079c602 100644
--- a/test/unit/test_morph_embedding.cpp
+++ b/test/unit/test_morph_embedding.cpp
@@ -170,9 +170,15 @@ TEST(embedding, partial_branch_length) {
EXPECT_DOUBLE_EQ(10., em.integrate_length(mcable{1, 0., 1.}));
EXPECT_DOUBLE_EQ(7.5, em.integrate_length(mcable{1, 0.25, 1.0}));
- // expect 2*0.25+3*0.5 = 2.0 times corresponding cable length.
+ // Expect 2*0.25+3*0.5 = 2.0 times corresponding cable length.
pw_elements<double> pw({0.25, 0.5, 1.}, {2., 3.});
EXPECT_DOUBLE_EQ(20., em.integrate_length(1, pw));
+
+ // Distamce between points on different branches:
+ ASSERT_EQ(3u, m.num_branches());
+ ASSERT_EQ(0u, m.branch_parent(2));
+ EXPECT_DOUBLE_EQ(em.integrate_length(mcable{0, 0.75, 1.})+em.integrate_length(mcable{2, 0, 0.5}),
+ em.integrate_length(mlocation{0, 0.75}, mlocation{2, 0.5}));
}
TEST(embedding, partial_area) {
@@ -233,6 +239,13 @@ TEST(embedding, partial_area) {
double expected_pw_area = 5.*sub_area1+7.*(sub_area2+sub_area3);
EXPECT_TRUE(near_relative(expected_pw_area, em.integrate_area(0, pw), reltol));
+ // Area between points on different branches:
+ ASSERT_EQ(3u, m.num_branches());
+ ASSERT_EQ(0u, m.branch_parent(2));
+ EXPECT_TRUE(near_relative(
+ em.integrate_area(mcable{0, 0.8, 1.})+em.integrate_area(mcable{2, 0, 0.3}),
+ em.integrate_area(mlocation{0, 0.8}, mlocation{2, 0.3}), reltol));
+
// Integrated inverse cross-sectional area in cable 1 from 0.1 to 0.4:
// radius râ‚€ = 9.5, râ‚ = 8, length = 3.
--
GitLab