diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index dd831a11743225f78c866ddc1462b732a90ebb7e..7403dc1e7ab5f7e9aeb3dc066abf3bcc0d87b7be 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -57,6 +57,10 @@ segment* cable_cell::segment(index_type index) {
assert_valid_segment(index);
return segments_[index].get();
}
+segment const* cable_cell::parent(index_type index) const {
+ assert_valid_segment(index);
+ return segments_[parents_[index]].get();
+}
segment const* cable_cell::segment(index_type index) const {
assert_valid_segment(index);
diff --git a/arbor/fvm_compartment.hpp b/arbor/fvm_compartment.hpp
index 0d626328764637a0ef36ba31190c08017d68aa49..4790f595e3f15d07fa192ccd870dda2a8ead0649 100644
--- a/arbor/fvm_compartment.hpp
+++ b/arbor/fvm_compartment.hpp
@@ -34,6 +34,16 @@ struct semi_compartment {
return *this;
}
+ static semi_compartment half_sphere(value_type r1) {
+ using namespace math;
+ return semi_compartment{
+ r1,
+ area_sphere(r1)/2,
+ volume_sphere(r1)/2,
+ {r1, r1}
+ };
+ }
+
static semi_compartment frustrum(value_type l, value_type r1, value_type r2) {
using namespace math;
return semi_compartment{
@@ -111,15 +121,17 @@ public:
// `lengths` must be a sequence of length at least one.
// `radii` must be a sequence of length `size(lengths)+1`.
template <typename Radii, typename Lengths>
- div_compartment_sampler(size_type n, const Radii& radii, const Lengths& lengths) {
+ div_compartment_sampler(size_type n, const Radii& radii, const Lengths& lengths, bool with_soma = false) {
// set up offset-to-subsegment lookup and interpolation
using namespace util;
+ with_soma_ = with_soma;
segs_ = make_partition(offsets_, lengths);
compat::compiler_barrier_if_icc_leq(20160415u);
nseg_ = size(segs_);
- scale_ = segs_.bounds().second/n;
+ scale_ = with_soma_? (segs_.bounds().second - segs_.front().second)/(n-1): segs_.bounds().second/n;
+
assign(radii_, radii);
arb_assert(size(radii_)==size(offsets_));
}
@@ -144,6 +156,8 @@ protected:
size_type i; // index
value_type p; // proportion [0,1] along sub-segment
+ sub_segment_index(): i(0), p(0) {}
+
sub_segment_index(size_type i_, value_type p_): i(i_), p(p_) {}
bool operator<(sub_segment_index x) const {
return i<x.i || (i==x.i && p<x.p);
@@ -169,6 +183,7 @@ protected:
return math::lerp(radii_[s.i], radii_[s.i+1], s.p);
}
+ bool with_soma_ = false; //segment passed in includes the soma information in radii[0] and lengths[0]
size_type nseg_ = 0;
std::vector<value_type> offsets_;
std::vector<value_type> radii_;
@@ -180,15 +195,33 @@ protected:
class div_compartment_integrator: public div_compartment_sampler {
public:
template <typename Radii, typename Lengths>
- div_compartment_integrator(size_type n, const Radii& radii, const Lengths& lengths):
- div_compartment_sampler(n, radii, lengths) {}
+ div_compartment_integrator(size_type n, const Radii& radii, const Lengths& lengths, bool with_soma = false):
+ div_compartment_sampler(n, radii, lengths, with_soma) {}
div_compartment operator()(size_type i) const {
using namespace math;
-
- auto sleft = locate(scale_*i);
- auto scentre = locate(scale_*(i+0.5));
- auto sright = locate(scale_*(i+1));
+ sub_segment_index sleft, scentre, sright;
+ // If segment includes the soma, divisions are specially calculated
+ if (with_soma_) {
+ auto soma_l = segs_.front().second;
+ if (i == 0) {
+ sleft = locate(soma_l/2);
+ scentre = locate(soma_l);
+ sright = locate(soma_l + scale_/4);
+ } else if (i == 1) {
+ sleft = locate(soma_l + scale_/4);
+ scentre = locate(soma_l + scale_/2);
+ sright = locate(soma_l + scale_);
+ } else {
+ sleft = locate(soma_l + scale_ * (i-1));
+ scentre = locate(soma_l + scale_ * (i-0.5));
+ sright = locate(soma_l + scale_ * i);
+ }
+ } else {
+ sleft = locate(scale_ * i);
+ scentre = locate(scale_ * (i + 0.5));
+ sright = locate(scale_ * (i + 1));
+ }
return div_compartment(i, integrate(sleft, scentre), integrate(scentre, sright));
}
@@ -199,6 +232,12 @@ protected:
auto seg = segs_[a.i];
auto l = (b.p-a.p)*(seg.second-seg.first);
+
+ // If segment includes the soma it will be at index 0
+ // The soma is represented as a sphere
+ if (with_soma_ && a.i==0) {
+ return semi_compartment::half_sphere(radii_.front());
+ }
return semi_compartment::frustrum(l, radius_at(a), radius_at(b));
}
@@ -206,6 +245,7 @@ protected:
sub_segment_index x = std::min(b, sub_segment_index(a.i, 1));
auto s = sub_segment_frustrum(a, x);
+
while (a.i<b.i) {
++a.i;
a.p = 0;
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index c350770a8b911b53bc2f2eebfd48b3f4936b2f55..dd4d5b51042516ff61153b06ccbe808ffaf4a07d 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -42,6 +42,11 @@ namespace {
explicit compartment_model(const cable_cell& cell) {
tree = arb::tree(cell.parents());
auto counts = cell.compartment_counts();
+ for (unsigned i = 0; i < cell.segments().size(); i++) {
+ if (!cell.segment(i)->is_soma() && cell.parent(i)->is_soma()) {
+ counts[i]++;
+ }
+ }
parent_index = make_parent_index(tree, counts);
segment_index = algorithms::make_index(counts);
}
@@ -129,41 +134,55 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
util::make_partition(D.cell_segment_bounds,
transform_view(cells, [](const cable_cell& c) { return c.num_segments(); }));
- std::vector<index_type> cell_comp_bounds;
- auto cell_comp_part = make_partition(cell_comp_bounds,
- transform_view(cells, [](const cable_cell& c) { return c.num_compartments(); }));
+ std::vector<index_type> cell_cv_bounds;
+ auto cell_cv_part = make_partition(cell_cv_bounds,
+ transform_view(cells, [](const cable_cell& c) {
+ unsigned ncv = 0;
+ for (unsigned i = 0; i < c.segments().size(); i++) {
+ ncv += c.segment(i)->num_compartments();
+ if (!c.segment(i)->is_soma() && c.parent(i)->is_soma()) {
+ ncv++;
+ }
+ }
+ return ncv;
+ }));
D.ncell = cells.size();
- D.ncomp = cell_comp_part.bounds().second;
+ D.ncv = cell_cv_part.bounds().second;
- D.face_conductance.assign(D.ncomp, 0.);
- D.cv_area.assign(D.ncomp, 0.);
- D.cv_capacitance.assign(D.ncomp, 0.);
- D.parent_cv.assign(D.ncomp, index_type(-1));
- D.cv_to_cell.resize(D.ncomp);
+ D.face_conductance.assign(D.ncv, 0.);
+ D.cv_area.assign(D.ncv, 0.);
+ D.cv_capacitance.assign(D.ncv, 0.);
+ D.parent_cv.assign(D.ncv, index_type(-1));
+ D.cv_to_cell.resize(D.ncv);
for (auto i: make_span(0, D.ncell)) {
- util::fill(subrange_view(D.cv_to_cell, cell_comp_part[i]), static_cast<index_type>(i));
+ util::fill(subrange_view(D.cv_to_cell, cell_cv_part[i]), static_cast<index_type>(i));
}
- std::vector<size_type> seg_comp_bounds;
+ std::vector<size_type> seg_cv_bounds;
for (auto i: make_span(0, D.ncell)) {
const auto& c = cells[i];
- compartment_model cell_graph (c);
- auto cell_comp_ival = cell_comp_part[i];
+ compartment_model cell_graph(c);
+ auto cell_cv_ival = cell_cv_part[i];
- auto cell_comp_base = cell_comp_ival.first;
- for (auto k: make_span(cell_comp_ival)) {
- D.parent_cv[k] = cell_graph.parent_index[k-cell_comp_base]+cell_comp_base;
+ auto cell_cv_base = cell_cv_ival.first;
+ for (auto k: make_span(cell_cv_ival)) {
+ D.parent_cv[k] = cell_graph.parent_index[k-cell_cv_base]+cell_cv_base;
}
// Compartment index range for each segment in this cell.
- seg_comp_bounds.clear();
- auto seg_comp_part = make_partition(
- seg_comp_bounds,
- transform_view(c.segments(), [](const segment_ptr& s) { return s->num_compartments(); }),
- cell_comp_base);
+ seg_cv_bounds.clear();
+ auto seg_cv_part = make_partition(
+ seg_cv_bounds,
+ transform_view(make_span(c.num_segments()), [&c](const unsigned s) {
+ if (!c.segment(s)->is_soma() && c.parent(s)->is_soma()) {
+ return c.segment(s)->num_compartments() + 1;
+ }
+ return c.segment(s)->num_compartments();
+ }),
+ cell_cv_base);
- const auto nseg = seg_comp_part.size();
+ const auto nseg = seg_cv_part.size();
if (nseg==0) {
throw arbor_internal_error("fvm_layout: cannot discretrize cell with no segments");
}
@@ -179,12 +198,9 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
segment_info soma_info;
- size_type soma_cv = cell_comp_base;
+ size_type soma_cv = cell_cv_base;
value_type soma_area = math::area_sphere(soma->radius());
- D.cv_area[soma_cv] = soma_area; // [µm²]
- D.cv_capacitance[soma_cv] = soma_area*soma->cm; // [pF]
-
soma_info.proximal_cv = soma_cv;
soma_info.distal_cv = soma_cv;
soma_info.distal_cv_area = soma_area;
@@ -192,8 +208,8 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
// Other segments must all be cable segments.
for (size_type j = 1; j<nseg; ++j) {
- const auto& seg_comp_ival = seg_comp_part[j];
- const auto ncomp = seg_comp_ival.second-seg_comp_ival.first;
+ const auto& seg_cv_ival = seg_cv_part[j];
+ const auto ncv = seg_cv_ival.second-seg_cv_ival.first;
segment_info seg_info;
@@ -204,19 +220,31 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
auto cm = cable->cm; // [F/m²]
auto rL = cable->rL; // [Ω·cm]
- auto divs = div_compartment_integrator(ncomp, cable->radii(), cable->lengths());
+ bool soma_parent = c.parent(j)->as_soma() ? true : false; //segment's parent is a soma
- seg_info.parent_cv = D.parent_cv[seg_comp_ival.first];
- seg_info.parent_cv_area = divs(0).left.area;
+ auto radii = cable->radii();
+ auto lengths = cable->lengths();
- seg_info.proximal_cv = seg_comp_ival.first;
- seg_info.distal_cv = seg_comp_ival.second-1;
- seg_info.distal_cv_area = divs(ncomp-1).right.area;
+ // If segment has soma parent, send soma information to div_compartment_integrator
+ if(soma_parent) {
+ radii.insert(radii.begin(), soma->radius());
+ lengths.insert(lengths.begin(), soma->radius()*2);
+ }
+
+ auto divs = div_compartment_integrator(ncv, radii, lengths, soma_parent);
+
+ seg_info.parent_cv = soma_parent ? seg_cv_ival.first : D.parent_cv[seg_cv_ival.first];
+ seg_info.parent_cv_area = soma_parent ? divs(0).right.area + divs(1).left.area : divs(0).left.area;
+ seg_info.soma_parent = soma_parent;
+
+ seg_info.proximal_cv = soma_parent ? seg_cv_ival.first + 1 : seg_cv_ival.first;
+ seg_info.distal_cv = seg_cv_ival.second-1;
+ seg_info.distal_cv_area = divs(ncv-1).right.area;
D.segments.push_back(seg_info);
- for (auto i: make_span(seg_comp_ival)) {
- const auto& div = divs(i-seg_comp_ival.first);
+ for (auto i: make_span(seg_cv_ival)) {
+ const auto& div = divs(i-seg_cv_ival.first);
auto j = D.parent_cv[i];
auto h1 = div.left.length; // [µm]
@@ -238,10 +266,13 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
D.cv_capacitance[i] += ar * cm; // [pF]
}
}
+
+ D.cv_area[soma_cv] = soma_area; // [µm²]
+ D.cv_capacitance[soma_cv] = soma_area*soma->cm; // [pF]
}
// Number of CVs per cell is exactly number of compartments.
- D.cell_cv_bounds = std::move(cell_comp_bounds);
+ D.cell_cv_bounds = std::move(cell_cv_bounds);
return D;
}
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index 90c202203e5fa696b781ada08b999fc1f8fc25c7..a87814ef9396dcdc901edc9f6503cde05f895495 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -18,6 +18,7 @@ struct segment_info {
using value_type = fvm_value_type;
using index_type = fvm_index_type;
+ bool soma_parent = false; // segments parent is soma
value_type parent_cv_area = 0;
value_type distal_cv_area = 0;
@@ -56,7 +57,7 @@ struct fvm_discretization {
using index_type = fvm_index_type; // In particular, used for CV indices.
size_type ncell;
- size_type ncomp;
+ size_type ncv;
// Note: if CV j has no parent, parent_cv[j] = j. TODO: confirm!
std::vector<index_type> parent_cv;
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 59400d301ea60b513a88608f90ef9ea5724f7cb5..b1b13406ca835ca02b6fe67a9964927d102979e3 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -375,7 +375,7 @@ void fvm_lowered_cell_impl<B>::initialize(
fvm_discretization D = fvm_discretize(cells);
- std::vector<index_type> cv_to_intdom(D.ncomp);
+ std::vector<index_type> cv_to_intdom(D.ncv);
std::transform(D.cv_to_cell.begin(), D.cv_to_cell.end(), cv_to_intdom.begin(),
[&cell_to_intdom](index_type i){ return cell_to_intdom[i]; });
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index 14be4c2392c863126419f42b372cc60adf23843d..b91a6e11f2434c10535c620a72620c6d375dcbee 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -164,6 +164,7 @@ public:
bool has_soma() const;
class segment* segment(index_type index);
+ const class segment* parent(index_type index) const;
const class segment* segment(index_type index) const;
/// access pointer to the soma
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index a7f471c8638de0326767cd532d4ea8c6a301d87b..ae394dd5f0158f6c7f607bca015da74f44c12cb3 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -151,7 +151,7 @@ TEST(fvm_layout, topology) {
// | 14 | 15 | 16 | 17|
EXPECT_EQ(2u, D.ncell);
- EXPECT_EQ(18u, D.ncomp);
+ EXPECT_EQ(20u, D.ncv);
unsigned nseg = 6;
EXPECT_EQ(nseg, D.segments.size());
@@ -161,11 +161,11 @@ TEST(fvm_layout, topology) {
ASSERT_EQ(D.ncell, D.cell_segment_part().size());
ASSERT_EQ(D.ncell, D.cell_cv_part().size());
- ASSERT_EQ(D.ncomp, D.parent_cv.size());
- ASSERT_EQ(D.ncomp, D.cv_to_cell.size());
- ASSERT_EQ(D.ncomp, D.face_conductance.size());
- ASSERT_EQ(D.ncomp, D.cv_area.size());
- ASSERT_EQ(D.ncomp, D.cv_capacitance.size());
+ ASSERT_EQ(D.ncv, D.parent_cv.size());
+ ASSERT_EQ(D.ncv, D.cv_to_cell.size());
+ ASSERT_EQ(D.ncv, D.face_conductance.size());
+ ASSERT_EQ(D.ncv, D.cv_area.size());
+ ASSERT_EQ(D.ncv, D.cv_capacitance.size());
// Partitions of CVs and segments by cell:
@@ -175,31 +175,32 @@ TEST(fvm_layout, topology) {
EXPECT_EQ(spair(0, 2), D.cell_segment_part()[0]);
EXPECT_EQ(spair(2, nseg), D.cell_segment_part()[1]);
- EXPECT_EQ(ipair(0, 5), D.cell_cv_part()[0]);
- EXPECT_EQ(ipair(5, D.ncomp), D.cell_cv_part()[1]);
+ EXPECT_EQ(ipair(0, 6), D.cell_cv_part()[0]);
+ EXPECT_EQ(ipair(6, D.ncv), D.cell_cv_part()[1]);
// Segment and CV parent relationships:
using ivec = std::vector<fvm_index_type>;
- EXPECT_EQ(ivec({0,0,1,2,3,5,5,6,7,8,9,10,11,12,9,14,15,16}), D.parent_cv);
+ EXPECT_EQ(ivec({0,0,1,2,3,4,6,6,7,8,9,10,11,12,13,14,11,16,17,18}), D.parent_cv);
+
EXPECT_FALSE(D.segments[0].has_parent());
- EXPECT_EQ(0, D.segments[1].parent_cv);
+ EXPECT_EQ(1, D.segments[1].parent_cv);
EXPECT_FALSE(D.segments[2].has_parent());
- EXPECT_EQ(5, D.segments[3].parent_cv);
- EXPECT_EQ(9, D.segments[4].parent_cv);
- EXPECT_EQ(9, D.segments[5].parent_cv);
+ EXPECT_EQ(7, D.segments[3].parent_cv);
+ EXPECT_EQ(11, D.segments[4].parent_cv);
+ EXPECT_EQ(11, D.segments[5].parent_cv);
// Segment CV ranges (half-open, exclusing parent):
EXPECT_EQ(ipair(0,1), D.segments[0].cv_range());
- EXPECT_EQ(ipair(1,5), D.segments[1].cv_range());
- EXPECT_EQ(ipair(5,6), D.segments[2].cv_range());
- EXPECT_EQ(ipair(6,10), D.segments[3].cv_range());
- EXPECT_EQ(ipair(10,14), D.segments[4].cv_range());
- EXPECT_EQ(ipair(14,18), D.segments[5].cv_range());
+ EXPECT_EQ(ipair(2,6), D.segments[1].cv_range());
+ EXPECT_EQ(ipair(6,7), D.segments[2].cv_range());
+ EXPECT_EQ(ipair(8,12), D.segments[3].cv_range());
+ EXPECT_EQ(ipair(12,16), D.segments[4].cv_range());
+ EXPECT_EQ(ipair(16,20), D.segments[5].cv_range());
// CV to cell index:
@@ -227,25 +228,27 @@ TEST(fvm_layout, area) {
}
unsigned n = 4; // compartments per dendritic segment
- EXPECT_FLOAT_EQ(A[0]+A[1]/(2*n), D.cv_area[0]);
- EXPECT_FLOAT_EQ(A[1]/n, D.cv_area[1]);
+ EXPECT_FLOAT_EQ(A[0], D.cv_area[0]);
+ EXPECT_FLOAT_EQ(A[1]/(2*n), D.cv_area[1]);
EXPECT_FLOAT_EQ(A[1]/n, D.cv_area[2]);
EXPECT_FLOAT_EQ(A[1]/n, D.cv_area[3]);
- EXPECT_FLOAT_EQ(A[1]/(2*n), D.cv_area[4]);
+ EXPECT_FLOAT_EQ(A[1]/n, D.cv_area[4]);
+ EXPECT_FLOAT_EQ(A[1]/(2*n), D.cv_area[5]);
- EXPECT_FLOAT_EQ(A[2]+A[3]/(2*n), D.cv_area[5]);
- EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[6]);
- EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[7]);
+ EXPECT_FLOAT_EQ(A[2], D.cv_area[6]);
+ EXPECT_FLOAT_EQ(A[3]/(2*n), D.cv_area[7]);
EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[8]);
- EXPECT_FLOAT_EQ((A[3]+A[4]+A[5])/(2*n), D.cv_area[9]);
- EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[10]);
- EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[11]);
+ EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[9]);
+ EXPECT_FLOAT_EQ(A[3]/n, D.cv_area[10]);
+ EXPECT_FLOAT_EQ((A[3]+A[4]+A[5])/(2*n), D.cv_area[11]);
EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[12]);
- EXPECT_FLOAT_EQ(A[4]/(2*n), D.cv_area[13]);
- EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[14]);
- EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[15]);
+ EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[13]);
+ EXPECT_FLOAT_EQ(A[4]/n, D.cv_area[14]);
+ EXPECT_FLOAT_EQ(A[4]/(2*n), D.cv_area[15]);
EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[16]);
- EXPECT_FLOAT_EQ(A[5]/(2*n), D.cv_area[17]);
+ EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[17]);
+ EXPECT_FLOAT_EQ(A[5]/n, D.cv_area[18]);
+ EXPECT_FLOAT_EQ(A[5]/(2*n), D.cv_area[19]);
// Confirm proportional allocation of surface capacitance:
@@ -258,14 +261,11 @@ TEST(fvm_layout, area) {
double cm3 = cells[1].segment(3)->cm;
double c = A[3]/(2*n)*cm1+A[4]/(2*n)*cm2+A[5]/(2*n)*cm3;
- EXPECT_FLOAT_EQ(c, D.cv_capacitance[9]);
-
- // CV 5 should be a weighted sum of soma and first segment
- // capacitcance from cell 1.
+ EXPECT_FLOAT_EQ(c, D.cv_capacitance[11]);
double cm0 = cells[1].soma()->cm;
- c = A[2]*cm0+A[3]/(2*n)*cm1;
- EXPECT_FLOAT_EQ(c, D.cv_capacitance[5]);
+ c = A[2]*cm0;
+ EXPECT_FLOAT_EQ(c, D.cv_capacitance[6]);
// Confirm face conductance within a constant diameter
// equals a/h·1/rL where a is the cross sectional
@@ -280,7 +280,7 @@ TEST(fvm_layout, area) {
double g = a/h/cable->rL; // [µm·S/cm]
g *= 100; // [µS]
- EXPECT_FLOAT_EQ(g, D.face_conductance[11]);
+ EXPECT_FLOAT_EQ(g, D.face_conductance[13]);
}
TEST(fvm_layout, mech_index) {
@@ -308,20 +308,20 @@ TEST(fvm_layout, mech_index) {
// Proportional area contrib: soma area/CV area.
EXPECT_EQ(mechanismKind::density, hh_config.kind);
- EXPECT_EQ(ivec({0,5}), hh_config.cv);
+ EXPECT_EQ(ivec({0,6}), hh_config.cv);
- fvec norm_area({area(cells[0].soma())/D.cv_area[0], area(cells[1].soma())/D.cv_area[5]});
+ fvec norm_area({area(cells[0].soma())/D.cv_area[0], area(cells[1].soma())/D.cv_area[6]});
EXPECT_TRUE(testing::seq_almost_eq<double>(norm_area, hh_config.norm_area));
// Three expsyn synapses, two 0.4 along segment 1, and one 0.4 along segment 5.
// These two synapses can be coalesced into 1 synapse
// 0.4 along => second (non-parent) CV for segment.
- EXPECT_EQ(ivec({2, 15}), expsyn_config.cv);
+ EXPECT_EQ(ivec({3, 17}), expsyn_config.cv);
// One exp2syn synapse, 0.4 along segment 4.
- EXPECT_EQ(ivec({11}), exp2syn_config.cv);
+ EXPECT_EQ(ivec({13}), exp2syn_config.cv);
// There should be a K and Na ion channel associated with each
// hh mechanism node.
@@ -330,8 +330,8 @@ TEST(fvm_layout, mech_index) {
ASSERT_EQ(1u, M.ions.count("k"s));
EXPECT_EQ(0u, M.ions.count("ca"s));
- EXPECT_EQ(ivec({0,5}), M.ions.at("na"s).cv);
- EXPECT_EQ(ivec({0,5}), M.ions.at("k"s).cv);
+ EXPECT_EQ(ivec({0,6}), M.ions.at("na"s).cv);
+ EXPECT_EQ(ivec({0,6}), M.ions.at("k"s).cv);
}
TEST(fvm_layout, coalescing_synapses) {
@@ -371,7 +371,7 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 2, 3, 4}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 3, 4, 5}), expsyn_config.cv);
EXPECT_EQ(ivec({1, 1, 1, 1}), expsyn_config.multiplicity);
}
{
@@ -387,17 +387,16 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 3}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 4}), expsyn_config.cv);
EXPECT_EQ(ivec({1, 1}), expsyn_config.multiplicity);
auto &exp2syn_config = M.mechanisms.at("exp2syn");
- EXPECT_EQ(ivec({2, 4}), exp2syn_config.cv);
+ EXPECT_EQ(ivec({3, 5}), exp2syn_config.cv);
EXPECT_EQ(ivec({1, 1}), exp2syn_config.multiplicity);
}
{
cable_cell cell = make_cell_ball_and_stick();
- // Add synapses of two varieties.
cell.add_synapse({1, 0.3}, "expsyn");
cell.add_synapse({1, 0.5}, "expsyn");
cell.add_synapse({1, 0.7}, "expsyn");
@@ -407,7 +406,7 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 2, 3, 4}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 3, 4, 5}), expsyn_config.cv);
EXPECT_TRUE(expsyn_config.multiplicity.empty());
}
{
@@ -423,11 +422,11 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 3}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 4}), expsyn_config.cv);
EXPECT_TRUE(expsyn_config.multiplicity.empty());
auto &exp2syn_config = M.mechanisms.at("exp2syn");
- EXPECT_EQ(ivec({2, 4}), exp2syn_config.cv);
+ EXPECT_EQ(ivec({3, 5}), exp2syn_config.cv);
EXPECT_TRUE(exp2syn_config.multiplicity.empty());
}
{
@@ -443,7 +442,7 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 3}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 4}), expsyn_config.cv);
EXPECT_EQ(ivec({2, 2}), expsyn_config.multiplicity);
}
{
@@ -459,7 +458,7 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 1, 3}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 2, 4}), expsyn_config.cv);
EXPECT_EQ(ivec({2, 1, 1}), expsyn_config.multiplicity);
EXPECT_EQ(fvec({0, 0.1, 0.1}), expsyn_config.param_values[0].second);
EXPECT_EQ(fvec({0.2, 0.2, 0.2}), expsyn_config.param_values[1].second);
@@ -481,7 +480,7 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 1, 3, 3}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 2, 4, 4}), expsyn_config.cv);
EXPECT_EQ(ivec({4, 6, 5, 7, 0, 2, 1, 3}), expsyn_config.target);
EXPECT_EQ(ivec({2, 2, 2, 2}), expsyn_config.multiplicity);
EXPECT_EQ(fvec({0, 1, 0, 1}), expsyn_config.param_values[0].second);
@@ -506,14 +505,14 @@ TEST(fvm_layout, coalescing_synapses) {
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
- EXPECT_EQ(ivec({1, 1}), expsyn_config.cv);
+ EXPECT_EQ(ivec({2, 2}), expsyn_config.cv);
EXPECT_EQ(ivec({0, 2, 5, 3}), expsyn_config.target);
EXPECT_EQ(ivec({3, 1}), expsyn_config.multiplicity);
EXPECT_EQ(fvec({1, 5}), expsyn_config.param_values[0].second);
EXPECT_EQ(fvec({2, 1}), expsyn_config.param_values[1].second);
auto &exp2syn_config = M.mechanisms.at("exp2syn");
- EXPECT_EQ(ivec({1, 1, 3, 3}), exp2syn_config.cv);
+ EXPECT_EQ(ivec({2, 2, 4, 4}), exp2syn_config.cv);
EXPECT_EQ(ivec({4, 1, 7, 8, 6, 9}), exp2syn_config.target);
EXPECT_EQ(ivec({1, 1, 2, 2}), exp2syn_config.multiplicity);
EXPECT_EQ(fvec({1, 4, 2, 2}), exp2syn_config.param_values[0].second);
@@ -671,7 +670,7 @@ TEST(fvm_layout, density_norm_area) {
seg.add_mechanism(hh);
}
- int ncv = 10;
+ int ncv = 11; //ncomp + 1
std::vector<double> expected_gkbar(ncv, dflt_gkbar);
std::vector<double> expected_gl(ncv, dflt_gl);
@@ -683,28 +682,30 @@ TEST(fvm_layout, density_norm_area) {
auto seg2_divs = div_by_ends(segs[2]->as_cable());
auto seg3_divs = div_by_ends(segs[3]->as_cable());
- // CV 0: mix of soma and left of segment 1
- expected_gl[0] = wmean(soma_area, dflt_gl, seg1_divs(0).left.area, seg1_gl);
-
+ // CV 0: soma
+ // CV1: left of segment 1
+ expected_gl[0] = dflt_gl;
expected_gl[1] = seg1_gl;
+
expected_gl[2] = seg1_gl;
+ expected_gl[3] = seg1_gl;
- // CV 3: mix of right of segment 1 and left of segments 2 and 3.
- expected_gkbar[3] = wmean(seg1_divs(2).right.area, dflt_gkbar, seg2_divs(0).left.area, seg2_gkbar, seg3_divs(0).left.area, seg3_gkbar);
- expected_gl[3] = wmean(seg1_divs(2).right.area, seg1_gl, seg2_divs(0).left.area, dflt_gl, seg3_divs(0).left.area, seg3_gl);
+ // CV 4: mix of right of segment 1 and left of segments 2 and 3.
+ expected_gkbar[4] = wmean(seg1_divs(2).right.area, dflt_gkbar, seg2_divs(0).left.area, seg2_gkbar, seg3_divs(0).left.area, seg3_gkbar);
+ expected_gl[4] = wmean(seg1_divs(2).right.area, seg1_gl, seg2_divs(0).left.area, dflt_gl, seg3_divs(0).left.area, seg3_gl);
- // CV 4-6: just segment 2
- expected_gkbar[4] = seg2_gkbar;
+ // CV 5-7: just segment 2
expected_gkbar[5] = seg2_gkbar;
expected_gkbar[6] = seg2_gkbar;
+ expected_gkbar[7] = seg2_gkbar;
- // CV 7-9: just segment 3
- expected_gkbar[7] = seg3_gkbar;
+ // CV 8-10: just segment 3
expected_gkbar[8] = seg3_gkbar;
expected_gkbar[9] = seg3_gkbar;
- expected_gl[7] = seg3_gl;
+ expected_gkbar[10] = seg3_gkbar;
expected_gl[8] = seg3_gl;
expected_gl[9] = seg3_gl;
+ expected_gl[10] = seg3_gl;
cable_cell_global_properties gprop;
fvm_discretization D = fvm_discretize(cells);
@@ -717,12 +718,14 @@ TEST(fvm_layout, density_norm_area) {
double area_relerr = 10*std::numeric_limits<double>::epsilon();
EXPECT_TRUE(testing::near_relative(D.cv_area[0],
- soma_area+seg1_divs(0).left.area, area_relerr));
+ soma_area, area_relerr));
EXPECT_TRUE(testing::near_relative(D.cv_area[1],
+ seg1_divs(0).left.area, area_relerr));
+ EXPECT_TRUE(testing::near_relative(D.cv_area[2],
seg1_divs(0).right.area+seg1_divs(1).left.area, area_relerr));
- EXPECT_TRUE(testing::near_relative(D.cv_area[3],
+ EXPECT_TRUE(testing::near_relative(D.cv_area[4],
seg1_divs(2).right.area+seg2_divs(0).left.area+seg3_divs(0).left.area, area_relerr));
- EXPECT_TRUE(testing::near_relative(D.cv_area[6],
+ EXPECT_TRUE(testing::near_relative(D.cv_area[7],
seg2_divs(2).right.area, area_relerr));
// Grab the HH parameters from the mechanism.
@@ -803,11 +806,11 @@ TEST(fvm_layout, ion_weights) {
};
ivec expected_ion_cv[] = {
- {0}, {0, 1, 2}, {1, 2, 3}, {0, 1, 2, 3}, {1, 3}
+ {0}, {0, 2, 3}, {2, 3, 4}, {0, 1, 2, 3, 4}, {2, 4}
};
fvec expected_iconc_norm_area[] = {
- {1./3}, {1./3, 1./2, 0.}, {1./4, 0., 0.}, {0., 0., 0., 0.}, {3./4, 0.}
+ {0.}, {0., 1./2, 0.}, {1./4, 0., 0.}, {0., 0., 0., 0., 0.}, {3./4, 0.}
};
cable_cell_global_properties gprop;
@@ -828,6 +831,7 @@ TEST(fvm_layout, ion_weights) {
auto& ca = M.ions.at("ca"s);
EXPECT_EQ(expected_ion_cv[run], ca.cv);
+
EXPECT_TRUE(testing::seq_almost_eq<fvm_value_type>(expected_iconc_norm_area[run], ca.iconc_norm_area));
EXPECT_TRUE(util::all_of(ca.econc_norm_area, [](fvm_value_type v) { return v==1.; }));
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index e6d574630b1af13dfea47c414c1e2f9c91d77ad0..7bd1b18cea57dddbb1b5d188178ff9bf44244551 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -227,7 +227,7 @@ TEST(fvm_lowered, matrix_init)
fvcell.initialize({0}, cable1d_recipe(cell), cell_to_intdom, targets, probe_map);
auto& J = fvcell.*private_matrix_ptr;
- EXPECT_EQ(J.size(), 11u);
+ EXPECT_EQ(J.size(), 12u);
// Test that the matrix is initialized with sensible values
@@ -315,7 +315,7 @@ TEST(fvm_lowered, stimulus) {
// delay | 5 | 1
// duration | 80 | 2
// amplitude | 0.3 | 0.1
- // CV | 4 | 0
+ // CV | 5 | 0
execution_context context;
@@ -326,7 +326,7 @@ TEST(fvm_lowered, stimulus) {
cells[0].add_stimulus({0,0.5}, {1., 2., 0.1});
const fvm_size_type soma_cv = 0u;
- const fvm_size_type tip_cv = 4u;
+ const fvm_size_type tip_cv = 5u;
// now we have two stims :
//
@@ -683,10 +683,10 @@ TEST(fvm_lowered, weighted_write_ion) {
const double con_int = 80;
const double con_ext = 120;
- // Ca ion reader test_kinlva on CV 1 and 2 via segment 2:
+ // Ca ion reader test_kinlva on CV 2 and 3 via segment 2:
c.segments()[2] ->add_mechanism("test_kinlva");
- // Ca ion writer test_ca on CV 1 and 3 via segment 3:
+ // Ca ion writer test_ca on CV 2 and 4 via segment 3:
c.segments()[3] ->add_mechanism("test_ca");
std::vector<target_handle> targets;
@@ -703,7 +703,7 @@ TEST(fvm_lowered, weighted_write_ion) {
ion.init_concentration();
std::vector<unsigned> ion_nodes = util::assign_from(ion.node_index_);
- std::vector<unsigned> expected_ion_nodes = {1, 2, 3};
+ std::vector<unsigned> expected_ion_nodes = {2, 3, 4};
EXPECT_EQ(expected_ion_nodes, ion_nodes);
std::vector<double> ion_iconc_weights = util::assign_from(ion.weight_Xi_);
@@ -785,11 +785,11 @@ TEST(fvm_lowered, gj_coords_simple) {
return g * 1e3 / D.cv_area[i];
};
- EXPECT_EQ(pair({4,8}), GJ[0].loc);
- EXPECT_EQ(weight(0.5, 4), GJ[0].weight);
+ EXPECT_EQ(pair({5,10}), GJ[0].loc);
+ EXPECT_EQ(weight(0.5, 5), GJ[0].weight);
- EXPECT_EQ(pair({8,4}), GJ[1].loc);
- EXPECT_EQ(weight(0.5, 8), GJ[1].weight);
+ EXPECT_EQ(pair({10,5}), GJ[1].loc);
+ EXPECT_EQ(weight(0.5, 10), GJ[1].weight);
}
TEST(fvm_lowered, gj_coords_complex) {
@@ -895,10 +895,10 @@ TEST(fvm_lowered, gj_coords_complex) {
return g * 1e3 / D.cv_area[i];
};
- std::vector<pair> expected_loc = {{4, 14}, {4,11}, {2,21}, {14, 4}, {11,4} ,{8,28}, {6, 24}, {21,2}, {28,8}, {24, 6}};
+ std::vector<pair> expected_loc = {{5, 16}, {5,13}, {3,24}, {16, 5}, {13,5} ,{10,31}, {8, 27}, {24,3}, {31,10}, {27, 8}};
std::vector<double> expected_weight = {
- weight(0.03, 4), weight(0.04, 4), weight(0.01, 2), weight(0.03, 14), weight(0.04, 11),
- weight(0.02, 8), weight(0.01, 6), weight(0.01, 21), weight(0.02, 28), weight(0.01, 24)
+ weight(0.03, 5), weight(0.04, 5), weight(0.01, 3), weight(0.03, 16), weight(0.04, 13),
+ weight(0.02, 10), weight(0.01, 8), weight(0.01, 24), weight(0.02, 31), weight(0.01, 27)
};
for (unsigned i = 0; i < GJ.size(); i++) {
@@ -912,7 +912,6 @@ TEST(fvm_lowered, gj_coords_complex) {
}
EXPECT_TRUE(found);
}
- std::cout << std::endl;
}
TEST(fvm_lowered, cell_group_gj) {