diff --git a/arbor/include/arbor/morph/embed_pwlin.hpp b/arbor/include/arbor/morph/embed_pwlin.hpp index d6947ac87ed594ea52918aaf5e59296120b7c9fc..0338186bb31fa1b45cde63100ec2f7522c019cb1 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 d81e9e4c6145c40c4223ed7149b7200abc3c33fc..eaf430a6f42bafc9fde9f7ae8d0c28eb30f059ad 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 0c72b62fec69a9819ae8e83fa7082c3f740ebd61..b079c602331c93a0354c54c9f83ae625ecce2b95 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.