diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 5fbfafbe20e42f52e83b4ebc87adaa713d12bd43..7eb69f6397841b1309ceb8d8721484a763b5060a 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -82,8 +82,6 @@ std::vector<V> unique_union(const std::vector<V>& a, const std::vector<V>& b) {
return u;
}
-} // anonymous namespace
-
// Convert mcable_map values to a piecewise function over an mcable.
// The projection gives the map from the values in the mcable_map to the values in the piecewise function.
template <typename T, typename U, typename Proj = get_value>
@@ -124,8 +122,14 @@ pw_elements<U> pw_over_cable(const mcable_map<T>& mm, mcable cable, U dflt_value
}
return pw;
}
+} // anonymous namespace
+
+
+// Building CV geometry
+// --------------------
// Construct cv_geometry for cell from locset describing CV boundary points.
+
cv_geometry cv_geometry_from_ends(const cable_cell& cell, const locset& lset) {
auto pop = [](auto& vec) { auto h = vec.back(); return vec.pop_back(), h; };
@@ -271,6 +275,7 @@ namespace impl {
}
// Merge CV geometry lists in-place.
+
cv_geometry& append(cv_geometry& geom, const cv_geometry& right) {
using impl::tail;
using impl::append_offset;
@@ -311,6 +316,7 @@ cv_geometry& append(cv_geometry& geom, const cv_geometry& right) {
}
// Combine two fvm_cv_geometry groups in-place.
+
fvm_cv_discretization& append(fvm_cv_discretization& dczn, const fvm_cv_discretization& right) {
append(dczn.geometry, right.geometry);
@@ -326,6 +332,10 @@ fvm_cv_discretization& append(fvm_cv_discretization& dczn, const fvm_cv_discreti
return dczn;
}
+
+// FVM discretization
+// ------------------
+
fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global_dflt) {
const auto& dflt = cell.default_parameters;
fvm_cv_discretization D;
@@ -446,8 +456,234 @@ fvm_cv_discretization fvm_cv_discretize(const std::vector<cable_cell>& cells,
return combined;
}
+// Voltage interpolation
+// ---------------------
+//
+// Interpolated voltages and axial current at a given site are determined
+// from 'voltage references'. A voltage reference is a CV from which the
+// membrane voltage is taken, and a location within that CV where the
+// voltage is deemed to be accurate.
+//
+// A CV that includes no fork points has one reference location which is
+// the centre of the CV (by branch length). Otherwise, every fork in a CV
+// is regarded as being a reference location.
+//
+// Voltage references should comprise adjacent CVs, however should the site
+// lie between fork points within the one CV, there is nothing to interpolate
+// and the voltage references will all come from the one CV containing the
+// site.
+
+struct voltage_reference {
+ fvm_index_type cv = -1;
+ mlocation loc;
+};
+
+struct voltage_reference_pair {
+ voltage_reference proximal;
+ voltage_reference distal;
+};
+
+// Collection of other locations that are coincident under projection.
+std::vector<mlocation> coincident_locations(const morphology& m, mlocation x) {
+ std::vector<mlocation> result;
+ if (x.pos==0) {
+ msize_t parent_bid = m.branch_parent(x.branch);
+ if (parent_bid!=mnpos) {
+ result.push_back({parent_bid, 1});
+ }
+ for (msize_t sibling_bid: m.branch_children(parent_bid)) {
+ if (sibling_bid!=x.branch) {
+ result.push_back({sibling_bid, 0});
+ }
+ }
+ }
+ else if (x.pos==1) {
+ for (msize_t child_bid: m.branch_children(x.branch)) {
+ result.push_back({child_bid, 0});
+ }
+ }
+ return result;
+}
+
+// Test if location intersects (sorted) sequence of cables.
+template <typename Seq>
+bool cables_intersect_location(Seq&& cables, mlocation x) {
+ struct cmp_branch {
+ bool operator()(const mcable& c, msize_t bid) const { return c.branch<bid; }
+ bool operator()(msize_t bid, const mcable& c) const { return bid<c.branch; }
+ };
+
+ using std::begin;
+ using std::end;
+ auto eqr = std::equal_range(begin(cables), end(cables), x.branch, cmp_branch{});
+
+ return util::any_of(util::make_range(eqr),
+ [&x](const mcable& c) { return c.prox_pos<=x.pos && x.pos<=c.dist_pos; });
+}
+
+voltage_reference_pair fvm_voltage_reference_points(const morphology& morph, const cv_geometry& geom, fvm_size_type cell_idx, mlocation site) {
+ voltage_reference site_ref, parent_ref, child_ref;
+ bool check_parent = true, check_child = true;
+ msize_t bid = site.branch;
+
+ // 'Simple' CVs contain no fork points, and are represented by a single cable.
+ auto cv_simple = [&geom](auto cv) { return geom.cables(cv).size()==1u; };
+
+ auto cv_midpoint = [&geom](auto cv) {
+ // Under assumption that CV is simple:
+ mcable c = geom.cables(cv).front();
+ return mlocation{c.branch, (c.prox_pos+c.dist_pos)/2};
+ };
+
+ auto cv_contains_fork = [&](auto cv, mlocation x) {
+ // CV contains fork if it intersects any location coincident with x
+ // other than x itselfv.
+
+ if (cv_simple(cv)) return false;
+ auto locs = coincident_locations(morph, x);
+
+ return util::any_of(locs, [&](mlocation y) { return cables_intersect_location(geom.cables(cv), y); });
+ };
+
+ site_ref.cv = geom.location_cv(cell_idx, site, cv_prefer::cv_empty);
+ if (cv_simple(site_ref.cv)) {
+ site_ref.loc = cv_midpoint(site_ref.cv);
+ }
+ else if (cv_contains_fork(site_ref.cv, mlocation{bid, 0})) {
+ site_ref.loc = mlocation{bid, 0};
+ check_parent = false;
+ }
+ else {
+ // CV not simple, and without head of branch as fork point, must contain
+ // tail of branch as a fork point.
+ arb_assert(cv_contains_fork(site_ref.cv, mlocation{bid, 1}));
+
+ site_ref.loc = mlocation{bid, 1};
+ check_child = false;
+ }
+
+ if (check_parent) {
+ parent_ref.cv = geom.cv_parent[site_ref.cv];
+ }
+ if (parent_ref.cv!=-1) {
+ parent_ref.loc = cv_simple(parent_ref.cv)? cv_midpoint(parent_ref.cv): mlocation{bid, 0};
+ arb_assert(parent_ref.loc.branch==bid);
+ }
+
+ if (check_child) {
+ for (auto child_cv: geom.children(site_ref.cv)) {
+ mcable child_prox_cable = geom.cables(child_cv).front();
+ if (child_prox_cable.branch==bid) {
+ child_ref.cv = child_cv;
+ break;
+ }
+ }
+ }
+ if (child_ref.cv!=-1) {
+ child_ref.loc = cv_simple(child_ref.cv)? cv_midpoint(child_ref.cv): mlocation{bid, 1};
+ arb_assert(child_ref.loc.branch==bid);
+ }
+
+ // If both child and parent references are possible, pick based on distallity with respect
+ // to the site_ref location.
+
+ if (child_ref.cv!=-1 && parent_ref.cv!=-1) {
+ if (site.pos<site_ref.loc.pos) child_ref.cv = -1; // i.e. use parent.
+ else parent_ref.cv = -1; // i.e. use child.
+ }
+
+ voltage_reference_pair result;
+ if (child_ref.cv!=-1) {
+ result.proximal = site_ref;
+ result.distal = child_ref;
+ }
+ else if (parent_ref.cv!=-1) {
+ result.proximal = parent_ref;
+ result.distal = site_ref;
+ }
+ else {
+ result.proximal = site_ref;
+ result.distal = site_ref;
+ }
+
+ return result;
+}
+
+// Interpolate membrane voltage from reference points in adjacent CVs.
+
+fvm_voltage_interpolant fvm_interpolate_voltage(const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx, mlocation site) {
+ auto& embedding = cell.embedding();
+ fvm_voltage_interpolant vi;
+
+ auto vrefs = fvm_voltage_reference_points(cell.morphology(), D.geometry, cell_idx, site);
+ vi.proximal_cv = vrefs.proximal.cv;
+ vi.distal_cv = vrefs.distal.cv;
+
+ arb_assert(vrefs.proximal.loc.branch==site.branch);
+ arb_assert(vrefs.distal.loc.branch==site.branch);
+
+ if (vrefs.proximal.cv==vrefs.distal.cv) { // (no interpolation)
+ vi.proximal_coef = 1.0;
+ vi.distal_coef = 0.0;
+ }
+ else {
+ msize_t bid = site.branch;
+
+ arb_assert(vrefs.proximal.loc.pos<vrefs.distal.loc.pos);
+ mcable rr_span = mcable{bid, vrefs.proximal.loc.pos , vrefs.distal.loc.pos};
+ double rr_resistance = embedding.integrate_ixa(rr_span, D.axial_resistivity[0].at(bid));
+
+ // Note: site is not necessarily distal to the most proximal reference point.
+ bool flip_rs = vrefs.proximal.loc.pos>site.pos;
+ mcable rs_span = flip_rs? mcable{bid, site.pos, vrefs.proximal.loc.pos}
+ : mcable{bid, vrefs.proximal.loc.pos, site.pos};
+
+ double rs_resistance = embedding.integrate_ixa(rs_span, D.axial_resistivity[0].at(bid));
+ if (flip_rs) {
+ rs_resistance = -rs_resistance;
+ }
+
+ double p = rs_resistance/rr_resistance;
+ vi.proximal_coef = 1-p;
+ vi.distal_coef = p;
+ }
+ return vi;
+}
+
+// Axial current as linear combination of membrane voltages at reference points in adjacent CVs.
+
+fvm_voltage_interpolant fvm_axial_current(const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx, mlocation site) {
+ auto& embedding = cell.embedding();
+ fvm_voltage_interpolant vi;
+
+ auto vrefs = fvm_voltage_reference_points(cell.morphology(), D.geometry, cell_idx, site);
+ vi.proximal_cv = vrefs.proximal.cv;
+ vi.distal_cv = vrefs.distal.cv;
+
+ if (vi.proximal_cv==vi.distal_cv) {
+ vi.proximal_coef = 0;
+ vi.distal_coef = 0;
+ }
+ else {
+ msize_t bid = site.branch;
+
+ arb_assert(vrefs.proximal.loc.pos<vrefs.distal.loc.pos);
+ mcable rr_span = mcable{bid, vrefs.proximal.loc.pos , vrefs.distal.loc.pos};
+ double rr_conductance = 100/embedding.integrate_ixa(rr_span, D.axial_resistivity[cell_idx].at(bid)); // [µS]
+
+ vi.proximal_coef = -rr_conductance;
+ vi.distal_coef = rr_conductance;
+ }
+
+ return vi;
+}
+
+// FVM mechanism data
+// ------------------
+
// CVs are absolute (taken from combined discretization) so do not need to be shifted.
// Only target numbers need to be shifted.
+
fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& right) {
using impl::append_offset;
@@ -514,6 +750,8 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
return combined;
}
+// Construct FVM mechanism data for a single cell.
+
fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop,
const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx)
{
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index 434ba6f5e80363941fe43f9e064e8480e5a3529d..e3ade27957e539c0a595d0bd37b4b4e1f2c6052f 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -203,6 +203,24 @@ fvm_cv_discretization& append(fvm_cv_discretization&, const fvm_cv_discretizatio
fvm_cv_discretization fvm_cv_discretize(const cable_cell& cell, const cable_cell_parameter_set& global_dflt);
fvm_cv_discretization fvm_cv_discretize(const std::vector<cable_cell>& cells, const cable_cell_parameter_set& global_defaults, const arb::execution_context& ctx={});
+
+// Interpolant data for voltage, axial current probes.
+//
+// The proximal CV will always be either the same as distal CV (trivial
+// interpolation) or the parent of the distal CV.
+
+struct fvm_voltage_interpolant {
+ fvm_index_type proximal_cv, distal_cv;
+ fvm_value_type proximal_coef, distal_coef;
+};
+
+// Interpolated membrane voltage.
+fvm_voltage_interpolant fvm_interpolate_voltage(const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx, mlocation site);
+
+// Axial current as linear combiantion of voltages.
+fvm_voltage_interpolant fvm_axial_current(const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx, mlocation site);
+
+
// Post-discretization data for point and density mechanism instantiation.
struct fvm_mechanism_config {
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index c09f25eef696e26127fb0da910b989f197796028..54403866cde46c114a53eb831ae9b7a00a601412 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -1,3 +1,4 @@
+#include <limits>
#include <string>
#include <vector>
@@ -13,6 +14,7 @@
#include "io/sepval.hpp"
#include "common.hpp"
+#include "common_morphologies.hpp"
#include "unit_test_catalogue.hpp"
#include "../common_cells.hpp"
@@ -835,3 +837,223 @@ TEST(fvm_layout, revpot) {
ASSERT_EQ(1u, M.mechanisms.count(write_eb_ec.name()));
EXPECT_EQ((std::vector<fvm_index_type>(1, soma1_index)), M.mechanisms.at(write_eb_ec.name()).cv);
}
+
+TEST(fvm_layout, vinterp_cable) {
+ // On a simple cable, expect CVs used forinterpolation to change at
+ // the midpoints of interior CVs. Every site in the proximal CV should
+ // interpolate between that and the next; every site in the distal CV
+ // 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);
+
+ // CV midpoints at branch pos 0.1, 0.3, 0.5, 0.7, 0.9.
+ // Expect voltage reference locations to be CV modpoints.
+ cell.default_parameters.discretization = cv_policy_fixed_per_branch(5);
+ fvm_cv_discretization D = fvm_cv_discretize(cell, neuron_parameter_defaults);
+
+ // Test locations, either side of CV midpoints plus extrema, CV boundaries.
+ double site_pos[] = { 0., 0.03, 0.11, 0.2, 0.28, 0.33, 0.4, 0.46, 0.55, 0.6, 0.75, 0.8, 0.83, 0.95, 1.};
+
+ for (auto pos: site_pos) {
+ mlocation site{0, pos};
+
+ fvm_index_type expected_distal;
+ if (pos<0.3) {
+ expected_distal = 1;
+ }
+ else if (pos<0.5) {
+ expected_distal = 2;
+ }
+ else if (pos<0.7) {
+ expected_distal = 3;
+ }
+ else {
+ expected_distal = 4;
+ }
+ fvm_index_type expected_proximal = expected_distal-1;
+
+ fvm_voltage_interpolant I = fvm_interpolate_voltage(cell, D, 0, site);
+
+ EXPECT_EQ(expected_proximal, I.proximal_cv);
+ EXPECT_EQ(expected_distal, I.distal_cv);
+
+ // Cable has constant diameter, so interpolation coefficients should
+ // be simple linear functions of branch position.
+
+ double prox_refpos = I.proximal_cv*0.2+0.1;
+ double dist_refpos = I.distal_cv*0.2+0.1;
+
+ // (Tortuous fp manipulation along the way makes the error greater than 4 ulp).
+ const double relerr = 32*std::numeric_limits<double>::epsilon();
+
+ EXPECT_TRUE(testing::near_relative((dist_refpos-pos)/0.2, I.proximal_coef, relerr));
+ EXPECT_TRUE(testing::near_relative((pos-prox_refpos)/0.2, I.distal_coef, relerr));
+ }
+}
+
+TEST(fvm_layout, vinterp_forked) {
+ // If a CV contains points at both ends of a branch, there will be
+ // no other adjacent CV on the same branch that we can use for
+ // 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);
+
+ // 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.
+ mlocation_list cv_ends{{1, 0.}, {2, 0.}};
+ cell.default_parameters.discretization = cv_policy_explicit(cv_ends);
+ fvm_cv_discretization D = fvm_cv_discretize(cell, neuron_parameter_defaults);
+
+ // Points in branch 0 should only get CV 0 for interpolation.
+ {
+ fvm_voltage_interpolant I = fvm_interpolate_voltage(cell, D, 0, mlocation{0, 0.3});
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_EQ(0, I.distal_cv);
+ EXPECT_EQ(1, I.proximal_coef+I.distal_coef);
+ }
+ // Points in branches 1 and 2 should get CV 0 and CV 1 or 2 respectively.
+ {
+ fvm_voltage_interpolant I = fvm_interpolate_voltage(cell, D, 0, mlocation{1, 0});
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_EQ(1., I.proximal_coef);
+ EXPECT_EQ(1, I.distal_cv);
+ EXPECT_EQ(0., I.distal_coef);
+
+ // Past the midpoint, we're extrapolating.
+ I = fvm_interpolate_voltage(cell, D, 0, mlocation{1, 0.7});
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_LT(I.proximal_coef, 0.);
+ EXPECT_EQ(1, I.distal_cv);
+ EXPECT_GT(I.distal_coef, 1.);
+
+ I = fvm_interpolate_voltage(cell, D, 0, mlocation{2, 0});
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_EQ(1., I.proximal_coef);
+ EXPECT_EQ(2, I.distal_cv);
+ EXPECT_EQ(0., I.distal_coef);
+
+ I = fvm_interpolate_voltage(cell, D, 0, mlocation{2, 0.7});
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_LT(I.proximal_coef, 0.);
+ EXPECT_EQ(2, I.distal_cv);
+ EXPECT_GT(I.distal_coef, 1.);
+ }
+}
+
+TEST(fvm_layout, iinterp) {
+ // If we get two distinct interpolation points back, the coefficients
+ // should match the face-conductance.
+
+ // 1. Vertex-delimited and vertex-centred discretizations.
+ using namespace common_morphology;
+
+ std::vector<cable_cell> cells;
+ std::vector<std::string> label;
+ for (auto& p: test_morphologies) {
+ if (p.second.empty()) continue;
+
+ cells.emplace_back(p.second);
+ cells.back().default_parameters.discretization = cv_policy_fixed_per_branch(3);
+ label.push_back(p.first+": forks-at-end"s);
+
+ cells.emplace_back(p.second);
+ cells.back().default_parameters.discretization = cv_policy_fixed_per_branch(3, cv_policy_flag::interior_forks);
+ label.push_back(p.first+": interior-forks"s);
+ }
+
+ fvm_cv_discretization D = fvm_cv_discretize(cells, neuron_parameter_defaults);
+ for (unsigned cell_idx = 0; cell_idx<cells.size(); ++cell_idx) {
+ SCOPED_TRACE(label[cell_idx]);
+ unsigned n_branch = D.geometry.n_branch(cell_idx);
+ for (msize_t bid = 0; bid<n_branch; ++bid) {
+ for (double pos: {0., 0.3, 0.4, 0.7, 1.}) {
+ mlocation x{bid, pos};
+ SCOPED_TRACE(x);
+
+ fvm_voltage_interpolant I = fvm_axial_current(cells[cell_idx], D, cell_idx, x);
+
+ // With the given discretization policies, should only have no interpolation when
+ // the cell has only the once CV.
+
+ if (D.geometry.cell_cvs(cell_idx).size()==1) {
+ EXPECT_EQ(I.proximal_cv, I.distal_cv);
+ EXPECT_EQ(D.geometry.cell_cvs(cell_idx).front(), I.proximal_cv);
+ }
+ else {
+ EXPECT_EQ(D.geometry.cv_parent.at(I.distal_cv), I.proximal_cv);
+ EXPECT_TRUE(I.proximal_cv>=D.geometry.cell_cv_interval(cell_idx).first);
+
+ double fc = D.face_conductance.at(I.distal_cv);
+ EXPECT_DOUBLE_EQ(-fc, I.proximal_coef);
+ EXPECT_DOUBLE_EQ(+fc, I.distal_coef);
+ }
+ }
+ }
+ }
+
+ // 2. Weird discretization: test points where the interpolated current has to be zero.
+ // 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);
+
+ // 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.
+ mlocation_list cv_ends{{1, 0.}, {2, 0.}};
+ cell.default_parameters.discretization = cv_policy_explicit(cv_ends);
+ D = fvm_cv_discretize(cell, neuron_parameter_defaults);
+
+ // Expect axial current interpolations on branches 1 and 2 to match CV 1 and 2
+ // face-conductances; CV 0 contains the fork point, so there is nothing to
+ // interpolate from on branch 0.
+
+ // Branch 0:
+ for (double pos: {0., 0.1, 0.8, 1.}) {
+ mlocation x{0, pos};
+ SCOPED_TRACE(x);
+
+ fvm_voltage_interpolant I = fvm_axial_current(cell, D, 0, x);
+
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_EQ(0, I.distal_cv);
+ EXPECT_EQ(0., I.proximal_coef);
+ EXPECT_EQ(0., I.distal_coef);
+ }
+
+ // Branch 1:
+ double fc1 = D.face_conductance[1];
+ for (double pos: {0., 0.1, 0.8, 1.}) {
+ mlocation x{1, pos};
+ SCOPED_TRACE(x);
+
+ fvm_voltage_interpolant I = fvm_axial_current(cell, D, 0, x);
+
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_EQ(1, I.distal_cv);
+ EXPECT_EQ(-fc1, I.proximal_coef);
+ EXPECT_EQ(+fc1, I.distal_coef);
+ }
+
+ // Branch 2:
+ double fc2 = D.face_conductance[2];
+ for (double pos: {0., 0.1, 0.8, 1.}) {
+ mlocation x{2, pos};
+ SCOPED_TRACE(x);
+
+ fvm_voltage_interpolant I = fvm_axial_current(cell, D, 0, x);
+
+ EXPECT_EQ(0, I.proximal_cv);
+ EXPECT_EQ(2, I.distal_cv);
+ EXPECT_EQ(-fc2, I.proximal_coef);
+ EXPECT_EQ(+fc2, I.distal_coef);
+ }
+}