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

@@ -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_;

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.

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.