From e7a8fb6f2d7a80fa6a5d687d731edb6a8a98c803 Mon Sep 17 00:00:00 2001
From: Sam Yates <sam@quux.dropbear.id.au>
Date: Tue, 11 Oct 2016 04:51:37 +0200
Subject: [PATCH] Complex compartments
* Use divided compartments to determine FVM coefficients.
* Pick correct control volume in FVM from sgement position (avoids
off-by-half error.)
* Add colour override functionality to tsplot: `--colour` option.
* Add const accessor for cell soma.
* Source formatting, comments in `math.hpp`
* Fix `range_view`: was using incorrectly named type trait.
* Add unit test for `range_view`.
* Allow points of discontinuity to be omitted from L-infinity norm
calculations.
* Add `-d, --min-dt` option to `validate.exe` to control time
step in validation convergence tests.
* Add validation test: confirm divided compartment policy does
not effect results on simple frustrum dendrites.
* Change default max compartments on validation tests to 100
(ad hoc observed convergence limit at dt circa 0.001 ms;
finder spatial division would required much finer dt.)
* Make NEURON validation data generation scripts use CVODE by
default, and with `secondorder=2` when non-zero `dt` is given.
---
scripts/tsplot | 31 ++++-
src/cell.cpp | 12 +-
src/cell.hpp | 7 +-
src/fvm_multicell.hpp | 130 +++++++++++++-----
src/math.hpp | 55 +++-----
src/util/rangeutil.hpp | 4 +-
tests/unit/test_range.cpp | 9 ++
tests/validation/CMakeLists.txt | 1 +
tests/validation/convergence_test.hpp | 22 ++-
tests/validation/trace_analysis.cpp | 39 ++++++
tests/validation/trace_analysis.hpp | 5 +
tests/validation/validate.cpp | 5 +
tests/validation/validate_ball_and_stick.cpp | 36 +++--
.../validate_compartment_policy.cpp | 95 +++++++++++++
tests/validation/validate_soma.cpp | 17 ++-
tests/validation/validate_synapses.cpp | 5 +-
tests/validation/validation_data.hpp | 6 +-
validation/ref/neuron/ball_and_3stick.py | 1 -
validation/ref/neuron/ball_and_squiggle.py | 1 -
validation/ref/neuron/generate_validation.sh | 2 +-
validation/ref/neuron/nrn_validation.py | 29 +++-
validation/ref/neuron/soma.py | 2 -
22 files changed, 394 insertions(+), 120 deletions(-)
create mode 100644 tests/validation/validate_compartment_policy.cpp
diff --git a/scripts/tsplot b/scripts/tsplot
index 4280a9dc..64186c90 100755
--- a/scripts/tsplot
+++ b/scripts/tsplot
@@ -30,6 +30,11 @@ def parse_clargs():
l, r = (float_or_none(x) for x in s.split(','))
return (l,r)
+ def parse_colour_spec(s):
+ colour, tests = s.split(':',1)
+ tests = tests.split(',')
+ return colour, tests
+
P = argparse.ArgumentParser(description='Plot time series data on one or more graphs.')
P.add_argument('inputs', metavar='FILE', nargs='+',
help='time series data in JSON format')
@@ -46,6 +51,10 @@ def parse_clargs():
type=lambda s: s.split(','),
action='append',
help='select only series matching EXPR')
+ P.add_argument('-c', '--colour', metavar='COLOUR:EXPR,...', dest='colours',
+ type=parse_colour_spec,
+ action='append',
+ help='use colour COLOUR a base for series matching EXPR')
P.add_argument('-o', '--output', metavar='FILE', dest='outfile',
help='save plot to file FILE')
P.add_argument('--dpi', metavar='NUM', dest='dpi',
@@ -143,10 +152,15 @@ def run_select(expr, v):
return test not in str(val)
else:
if isinstance(val, numbers.Number):
- try:
- test=int(test)
- except ValueError:
- test=float(test)
+ if re.match(r'true$', test, re.I):
+ test=True
+ elif re.match(r'false$', test, re.I):
+ test=False
+ else:
+ try:
+ test=int(test)
+ except ValueError:
+ test=float(test)
if op=='=':
return val==test
@@ -311,7 +325,7 @@ def round_numeric_(x):
if not isinstance(x,float): return x
return "{:6g}".format(x)
-def plot_plots(plot_groups, axis='time', save=None, dpi=None, scale=None):
+def plot_plots(plot_groups, axis='time', colour_overrides=[], save=None, dpi=None, scale=None):
nplots = len(plot_groups)
plot_groups = sorted(plot_groups, key=lambda g: g.group_label())
@@ -392,6 +406,10 @@ def plot_plots(plot_groups, axis='time', save=None, dpi=None, scale=None):
line = next(lines)
for s, l in series_by_unit[ui]:
c = next(colours)
+ for colour, tests in colour_overrides:
+ if all([run_select(t, s.meta) for t in tests]):
+ c = colour
+
plot.plot(s.t, s.y, color=c, ls=line, label=l)
# treat exluded points especially
ex = s.excluded()
@@ -442,4 +460,5 @@ plots = gather_ts_plots(tss, groupby)
if not args.outfile:
M.interactive(False)
-plot_plots(plots, axis=axis, save=args.outfile, dpi=args.dpi, scale=args.scale)
+colours = args.colours if args.colours else []
+plot_plots(plots, axis=axis, colour_overrides=colours, save=args.outfile, dpi=args.dpi, scale=args.scale)
diff --git a/src/cell.cpp b/src/cell.cpp
index f247a2a8..12035ad6 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -101,12 +101,12 @@ bool cell::has_soma() const
return !segment(0)->is_placeholder();
}
-soma_segment* cell::soma()
-{
- if(has_soma()) {
- return segment(0)->as_soma();
- }
- return nullptr;
+soma_segment* cell::soma() {
+ return has_soma()? segment(0)->as_soma(): nullptr;
+}
+
+const soma_segment* cell::soma() const {
+ return has_soma()? segment(0)->as_soma(): nullptr;
}
cable_segment* cell::cable(index_type index)
diff --git a/src/cell.hpp b/src/cell.hpp
index c691e684..7fdbc8f7 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -110,7 +110,7 @@ public:
/// parent is the index of the parent segment for the cable section
/// args are the arguments to be used to consruct the new cable
template <typename... Args>
- cable_segment* add_cable(index_type parent, Args ...args);
+ cable_segment* add_cable(index_type parent, Args&&... args);
/// the number of segments in the cell
size_type num_segments() const;
@@ -118,11 +118,12 @@ public:
bool has_soma() const;
class segment* segment(index_type index);
- class segment const* segment(index_type index) const;
+ const class segment* segment(index_type index) const;
/// access pointer to the soma
/// returns nullptr if the cell has no soma
soma_segment* soma();
+ const soma_segment* soma() const;
/// access pointer to a cable segment
/// will throw an std::out_of_range exception if
@@ -229,7 +230,7 @@ bool cell_basic_equality(cell const& lhs, cell const& rhs);
// create a cable by forwarding cable construction parameters provided by the user
template <typename... Args>
-cable_segment* cell::add_cable(cell::index_type parent, Args ...args)
+cable_segment* cell::add_cable(cell::index_type parent, Args&&... args)
{
// check for a valid parent id
if(parent>=num_segments()) {
diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp
index e89ca117..ee3363d4 100644
--- a/src/fvm_multicell.hpp
+++ b/src/fvm_multicell.hpp
@@ -9,6 +9,7 @@
#include <algorithms.hpp>
#include <cell.hpp>
+#include <compartment.hpp>
#include <event_queue.hpp>
#include <ion.hpp>
#include <math.hpp>
@@ -31,7 +32,20 @@ namespace nest {
namespace mc {
namespace fvm {
-template <typename Value, typename Index>
+inline int find_cv_index(const segment_location& loc, const compartment_model& graph) {
+ const auto& si = graph.segment_index;
+ const auto seg = loc.segment;
+
+ auto first = si[seg];
+ auto n = si[seg+1] - first;
+
+ int index = static_cast<int>(n*loc.position+0.5);
+ index = index==0? graph.parent_index[first]: first+(index-1);
+
+ return index;
+};
+
+template <typename Value, typename Index, typename DivPolicy=div_compartment_integrator>
class fvm_multicell {
public:
fvm_multicell() = default;
@@ -233,8 +247,8 @@ private:
//////////////////////////////// Implementation ////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
-template <typename T, typename I>
-void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
+template <typename T, typename I, typename DivPolicy>
+void fvm_multicell<T, I, DivPolicy>::compute_cv_area_unnormalized_capacitance(
std::pair<size_type, size_type> comp_ival,
const segment* seg,
index_type &parent)
@@ -260,38 +274,80 @@ void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
cv_capacitance_[i] += area * soma->mechanism("membrane").get("c_m").value;
}
else if (auto cable = seg->as_cable()) {
- // loop over each compartment in the cable
- // each compartment has the face between two CVs at its centre
- // the centers of the CVs are the end points of the compartment
+ // Loop over each compartment in the cable
+ //
+ // Each compartment i straddles the ith control volume on the right
+ // and the jth control volume on the left, where j is the parent index
+ // of i.
//
- // __________________________________
- // | ........ | .cvleft. | cv |
- // | ........ L ........ C R
- // |__________|__________|__________|
+ // Dividing the comparment into two halves, the centre face C
+ // corresponds to the shared face between the two control volumes,
+ // the surface areas in each half contribute to the surface area of
+ // the respective control volumes, and the volumes and lengths of
+ // each half are used to calculate the flux coefficients that
+ // for the connection between the two control volumes and which
+ // (after scaling by inverse capacitance) is stored in
+ // `face_alpha[i]`.
+ //
//
- // The compartment has end points marked L and R (left and right).
- // The left compartment is assumed to be closer to the soma
- // (i.e. it follows the minimal degree ordering)
- // The face is at the center, marked C.
- // The full control volume to the left (marked with .)
+ // +------- cv j --------+------- cv i -------+
+ // | | |
+ // v v v
+ // ____________________________________________
+ // | ........ | ........ | | |
+ // | ........ L ........ C R |
+ // |__________|__________|__________|_________|
+ // ^ ^
+ // | |
+ // +--- compartment i ---+
+ //
+ // The first control volume of any cell corresponds to the soma
+ // and the first half of the first cable compartment of that cell.
auto c_m = cable->mechanism("membrane").get("c_m").value;
auto r_L = cable->mechanism("membrane").get("r_L").value;
- const auto& compartments = cable->compartments();
- EXPECTS(size(compartments)==ncomp);
+ auto divs = div_compartments<DivPolicy>(cable, ncomp);
for (auto i: util::make_span(comp_ival)) {
- const auto& c = compartments[i-comp_ival.first];
+ const auto& div = divs(i-comp_ival.first);
auto j = parent[i];
- auto radius_center = math::mean(c.radius);
- auto area_face = math::area_circle(radius_center);
- face_alpha_[i] = area_face / (c_m * r_L * c.length);
-
- auto halflen = c.length/2;
- auto al = math::area_frustrum(halflen, left(c.radius), radius_center);
- auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
+ // Conductance approximated by weighted harmonic mean of mean
+ // conductances in each half.
+ //
+ // Mean conductances:
+ // c₠= 1/h₠∫₠A(x)/R dx
+ // c₂ = 1/h₂ ∫₂ A(x)/R dx
+ //
+ // where A(x) is the cross-sectional area, R is the bulk
+ // resistivity, h is the length of the interval and the
+ // integrals are taken over the intervals respectively.
+ // Equivalently, in terms of the semi-compartment volumes
+ // Vâ‚ and Vâ‚‚:
+ //
+ // câ‚ = 1/R·Vâ‚/hâ‚
+ // c₂ = 1/R·V₂/h₂
+ //
+ // Weighted harmonic mean, with h = hâ‚+hâ‚‚:
+ //
+ // c = (hâ‚/h·c₯¹+hâ‚‚/h·c₂¯¹)¯¹
+ // = 1/R · hVâ‚Vâ‚‚/(h₂²Vâ‚+h₲Vâ‚‚)
+
+ auto h1 = div.left.length;
+ auto V1 = div.left.volume;
+ auto h2 = div.right.length;
+ auto V2 = div.right.volume;
+ auto h = h1+h2;
+
+ auto conductance = 1/r_L*h*V1*V2/(h2*h2*V1+h1*h1*V2);
+ face_alpha_[i] = conductance / (c_m * h);
+
+ // auto halflen = c.length/2;
+ // auto al = math::area_frustrum(halflen, left(c.radius), radius_center);
+ // auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
+ auto al = div.left.area;
+ auto ar = div.right.area;
cv_areas_[j] += al;
cv_areas_[i] += ar;
@@ -304,9 +360,9 @@ void fvm_multicell<T, I>::compute_cv_area_unnormalized_capacitance(
}
}
-template <typename T, typename I>
+template <typename T, typename I, typename DivPolicy>
template <typename Cells, typename Detectors, typename Targets, typename Probes>
-void fvm_multicell<T, I>::initialize(
+void fvm_multicell<T, I, DivPolicy>::initialize(
const Cells& cells,
Detectors& detector_handles,
Targets& target_handles,
@@ -401,7 +457,7 @@ void fvm_multicell<T, I>::initialize(
auto& map_entry = syn_mech_map[syn_mech_index];
- size_type syn_comp = comp_ival.first+find_compartment_index(syn.location, graph);
+ size_type syn_comp = comp_ival.first+find_cv_index(syn.location, graph);
map_entry.push_back(syn_comp);
}
@@ -412,7 +468,7 @@ void fvm_multicell<T, I>::initialize(
// add the stimuli
for (const auto& stim: c.stimuli()) {
- auto idx = comp_ival.first+find_compartment_index(stim.location, graph);
+ auto idx = comp_ival.first+find_cv_index(stim.location, graph);
stimuli_.push_back({idx, stim.clamp});
}
@@ -420,7 +476,7 @@ void fvm_multicell<T, I>::initialize(
for (const auto& detector: c.detectors()) {
EXPECTS(detectors_count < detectors_size);
- auto comp = comp_ival.first+find_compartment_index(detector.location, graph);
+ auto comp = comp_ival.first+find_cv_index(detector.location, graph);
*detector_hi++ = comp;
++detectors_count;
}
@@ -429,7 +485,7 @@ void fvm_multicell<T, I>::initialize(
for (const auto& probe: c.probes()) {
EXPECTS(probes_count < probes_size);
- auto comp = comp_ival.first+find_compartment_index(probe.location, graph);
+ auto comp = comp_ival.first+find_cv_index(probe.location, graph);
switch (probe.kind) {
case probeKind::membrane_voltage:
*probe_hi++ = {&fvm_multicell::voltage_, comp};
@@ -555,8 +611,8 @@ void fvm_multicell<T, I>::initialize(
reset();
}
-template <typename T, typename I>
-void fvm_multicell<T, I>::setup_matrix(T dt) {
+template <typename T, typename I, typename DivPolicy>
+void fvm_multicell<T, I, DivPolicy>::setup_matrix(T dt) {
using memory::all;
// convenience accesors to matrix storage
@@ -599,8 +655,8 @@ void fvm_multicell<T, I>::setup_matrix(T dt) {
}
}
-template <typename T, typename I>
-void fvm_multicell<T, I>::reset() {
+template <typename T, typename I, typename DivPolicy>
+void fvm_multicell<T, I, DivPolicy>::reset() {
voltage_(memory::all) = resting_potential_;
t_ = 0.;
for (auto& m : mechanisms_) {
@@ -608,8 +664,8 @@ void fvm_multicell<T, I>::reset() {
}
}
-template <typename T, typename I>
-void fvm_multicell<T, I>::advance(T dt) {
+template <typename T, typename I, typename DivPolicy>
+void fvm_multicell<T, I, DivPolicy>::advance(T dt) {
using memory::all;
PE("current");
diff --git a/src/math.hpp b/src/math.hpp
index 36767bcb..50d760ef 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -9,90 +9,73 @@ namespace mc {
namespace math {
template <typename T>
-T constexpr pi()
-{
+T constexpr pi() {
return T(3.1415926535897932384626433832795l);
}
template <typename T = float>
-T constexpr infinity()
-{
+T constexpr infinity() {
return std::numeric_limits<T>::infinity();
}
template <typename T>
-T constexpr mean(T a, T b)
-{
+T constexpr mean(T a, T b) {
return (a+b) / T(2);
}
template <typename T>
-T constexpr mean(std::pair<T,T> const& p)
-{
+T constexpr mean(std::pair<T,T> const& p) {
return (p.first+p.second) / T(2);
}
template <typename T>
-T constexpr square(T a)
-{
+T constexpr square(T a) {
return a*a;
}
template <typename T>
-T constexpr cube(T a)
-{
+T constexpr cube(T a) {
return a*a*a;
}
-// area of a circle
-// pi r-squared
+// Area of circle radius r.
template <typename T>
-T constexpr area_circle(T r)
-{
+T constexpr area_circle(T r) {
return pi<T>() * square(r);
}
-// volume of a conic frustrum
+// Surface area of conic frustrum excluding the discs at each end,
+// with length L, end radii r1, r2.
template <typename T>
-T constexpr area_frustrum(T L, T r1, T r2)
-{
+T constexpr area_frustrum(T L, T r1, T r2) {
return pi<T>() * (r1+r2) * sqrt(square(L) + square(r1-r2));
}
-// surface area of conic frustrum, not including the area of the
-// circles at either end (i.e. just the area of the surface of rotation)
+// Volume of conic frustrum of length L, end radii r1, r2.
template <typename T>
-T constexpr volume_frustrum(T L, T r1, T r2)
-{
- // this formulation uses one less multiplication
+T constexpr volume_frustrum(T L, T r1, T r2) {
return pi<T>()/T(3) * (square(r1+r2) - r1*r2) * L;
- //return pi<T>()/T(3) * (square(r1) + r1*r2 + square(r2)) * L;
}
-// volume of a sphere
-// four-thirds pi r-cubed
+// Volume of a sphere radius r.
template <typename T>
-T constexpr volume_sphere(T r)
-{
+T constexpr volume_sphere(T r) {
return T(4)/T(3) * pi<T>() * cube(r);
}
-// surface area of a sphere
-// four pi r-squared
+// Surface area of a sphere radius r.
template <typename T>
-T constexpr area_sphere(T r)
-{
+T constexpr area_sphere(T r) {
return T(4) * pi<T>() * square(r);
}
-// linear interpolation in interval [a,b]
+// Linear interpolation by u in interval [a,b]: (1-u)*a + u*b.
template <typename T, typename U>
T constexpr lerp(T a, T b, U u) {
- // (1-u)*a + u*b
return std::fma(u, b, std::fma(-u, a, a));
}
-// -1, 0 or 1 according to sign of parameter
+// Return -1, 0 or 1 according to sign of parameter.
template <typename T>
int signum(T x) {
return (x<T(0)) - (x>T(0));
diff --git a/src/util/rangeutil.hpp b/src/util/rangeutil.hpp
index bc47dec2..855cbbec 100644
--- a/src/util/rangeutil.hpp
+++ b/src/util/rangeutil.hpp
@@ -33,14 +33,14 @@ range<const T*> singleton_view(const T& item) {
// Non-owning views and subviews
template <typename Seq>
-range<typename sequence_traits<Seq>::iterator_type, typename sequence_traits<Seq>::sentinel_type>
+range<typename sequence_traits<Seq>::iterator, typename sequence_traits<Seq>::sentinel>
range_view(Seq& seq) {
return make_range(std::begin(seq), std::end(seq));
}
template <
typename Seq,
- typename Iter = typename sequence_traits<Seq>::iterator_type,
+ typename Iter = typename sequence_traits<Seq>::iterator,
typename Size = typename sequence_traits<Seq>::size_type
>
enable_if_t<is_forward_iterator<Iter>::value, range<Iter>>
diff --git a/tests/unit/test_range.cpp b/tests/unit/test_range.cpp
index 2acf1412..3e014126 100644
--- a/tests/unit/test_range.cpp
+++ b/tests/unit/test_range.cpp
@@ -146,6 +146,15 @@ TEST(range, const_iterator) {
EXPECT_TRUE((std::is_same<const int&, decltype(r_const.front())>::value));
}
+TEST(range, view) {
+ std::vector<int> xs = { 1, 2, 3, 4, 5 };
+ auto r = util::range_view(xs);
+
+ r[3] = 7;
+ std::vector<int> check = { 1, 2, 3, 7, 5 };
+ EXPECT_EQ(check, xs);
+}
+
TEST(range, sentinel) {
const char *cstr = "hello world";
std::string s;
diff --git a/tests/validation/CMakeLists.txt b/tests/validation/CMakeLists.txt
index 783ef003..bbc03f1d 100644
--- a/tests/validation/CMakeLists.txt
+++ b/tests/validation/CMakeLists.txt
@@ -1,6 +1,7 @@
set(VALIDATION_SOURCES
# unit tests
validate_ball_and_stick.cpp
+ validate_compartment_policy.cpp
validate_soma.cpp
validate_synapses.cpp
diff --git a/tests/validation/convergence_test.hpp b/tests/validation/convergence_test.hpp
index 71700ea7..f0980dfd 100644
--- a/tests/validation/convergence_test.hpp
+++ b/tests/validation/convergence_test.hpp
@@ -70,7 +70,7 @@ public:
}
template <typename Model>
- void run(Model& m, Param p, float t_end, float dt) {
+ void run(Model& m, Param p, float t_end, float dt, const std::vector<float>& excl={}) {
// reset samplers and attach to probe locations
for (auto& se: cell_samplers_) {
se.sampler.reset();
@@ -91,7 +91,7 @@ public:
// compute metrics
if (run_validation_) {
- double linf = linf_distance(trace, ref_data_[label]);
+ double linf = linf_distance(trace, ref_data_[label], excl);
auto pd = peak_delta(trace, ref_data_[label]);
conv_tbl_.push_back({label, p, linf, pd});
@@ -116,5 +116,23 @@ public:
}
};
+/*
+ * Extract time points to exclude from current stimulus end-points.
+ */
+
+inline std::vector<float> stimulus_ends(const cell& c) {
+ std::vector<float> ts;
+
+ for (const auto& stimulus: c.stimuli()) {
+ float t0 = stimulus.clamp.delay();
+ float t1 = t0+stimulus.clamp.duration();
+ ts.push_back(t0);
+ ts.push_back(t1);
+ }
+
+ util::sort(ts);
+ return ts;
+}
+
} // namespace mc
} // namespace nest
diff --git a/tests/validation/trace_analysis.cpp b/tests/validation/trace_analysis.cpp
index a754cfe2..ba024ea6 100644
--- a/tests/validation/trace_analysis.cpp
+++ b/tests/validation/trace_analysis.cpp
@@ -48,6 +48,45 @@ double linf_distance(const trace_data& u, const trace_data& ref) {
[&](trace_entry x) { return std::abs(x.v-f(x.t)); }));
}
+/*
+ * Compute linf distance as above, but excluding sample points that lie
+ * near points in `excl`.
+ *
+ * `excl` contains the times to exclude, in ascending order.
+ */
+
+double linf_distance(const trace_data& u, const trace_data& ref, const std::vector<float>& excl) {
+ trace_interpolant f{ref};
+
+ trace_data reduced;
+ unsigned nexcl = excl.size();
+ unsigned ei = 0;
+
+ unsigned nu = u.size();
+ unsigned ui = 0;
+
+ while (ei<nexcl && ui<nu) {
+ float t = excl[ei++];
+
+ unsigned uj = ui;
+ while (uj<nu && u[uj].t<t) ++uj;
+
+ // include points up to and including uj-2, and then proceed from point uj+1,
+ // excluding the two points closest to the discontinuity.
+
+ if (uj>1+ui) {
+ util::append(reduced, util::subrange_view(u, ui, uj-1));
+ }
+ ui = uj+1;
+ }
+
+ if (ui<nu) {
+ util::append(reduced, util::subrange_view(u, ui, nu));
+ }
+
+ return linf_distance(reduced, ref);
+}
+
std::vector<trace_peak> local_maxima(const trace_data& u) {
std::vector<trace_peak> peaks;
if (u.size()<2) return peaks;
diff --git a/tests/validation/trace_analysis.hpp b/tests/validation/trace_analysis.hpp
index 6d9c921d..493a15f7 100644
--- a/tests/validation/trace_analysis.hpp
+++ b/tests/validation/trace_analysis.hpp
@@ -21,6 +21,11 @@ namespace mc {
double linf_distance(const trace_data& u, const trace_data& ref);
+// Compute linf distance as above, excluding samples near
+// times given in `excl`, monotonically increasing.
+
+double linf_distance(const trace_data& u, const trace_data& ref, const std::vector<float>& excl);
+
// Find local maxima (peaks) in a trace, excluding end points.
struct trace_peak {
diff --git a/tests/validation/validate.cpp b/tests/validation/validate.cpp
index de70f9bf..9f071aa2 100644
--- a/tests/validation/validate.cpp
+++ b/tests/validation/validate.cpp
@@ -21,6 +21,8 @@ const char* usage_str =
" -p, --path=DIR Look for validation reference data in DIR\n"
" -m, --max-comp=N Run convergence tests to a maximum of N\n"
" compartments per segment\n"
+" -d, --min-dt=DT Run convergence tests with or with a minimumf\n"
+" timestep DT\n"
" -h, --help Display usage information and exit\n";
int main(int argc, char **argv) {
@@ -41,6 +43,9 @@ int main(int argc, char **argv) {
else if (auto o = parse_opt<int>(arg, 'm', "max-comp")) {
g_trace_io.set_max_ncomp(*o);
}
+ else if (auto o = parse_opt<float>(arg, 'd', "min-dt")) {
+ g_trace_io.set_min_dt(*o);
+ }
else if (auto o = parse_opt<void>(arg, 'v', "verbose")) {
g_trace_io.set_verbose(true);
}
diff --git a/tests/validation/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp
index f3ca4c94..015aad96 100644
--- a/tests/validation/validate_ball_and_stick.cpp
+++ b/tests/validation/validate_ball_and_stick.cpp
@@ -26,17 +26,21 @@ void run_ncomp_convergence_test(
const util::path& ref_data_path,
const cell& c,
SamplerInfoSeq& samplers,
- float t_end=100.f,
- float dt=0.001)
+ float t_end=100.f)
{
auto max_ncomp = g_trace_io.max_ncomp();
+ auto dt = g_trace_io.min_dt();
+
nlohmann::json meta = {
{"name", "membrane voltage"},
{"model", model_name},
+ {"dt", dt},
{"sim", "nestmc"},
{"units", "mV"}
};
+ auto exclude = stimulus_ends(c);
+
convergence_test_runner<int> R("ncomp", samplers, meta);
R.load_reference_data(ref_data_path);
@@ -48,13 +52,13 @@ void run_ncomp_convergence_test(
}
model<lowered_cell> m(singleton_recipe{c});
- R.run(m, ncomp, t_end, dt);
+ R.run(m, ncomp, t_end, dt, exclude);
}
R.report();
R.assert_all_convergence();
}
-TEST(ball_and_stick, neuron_ref) {
+TEST(ball_and_taper, neuron_ref) {
using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
cell c = make_cell_ball_and_stick();
@@ -123,23 +127,29 @@ TEST(rallpack1, numeric_ref) {
250.f);
}
-TEST(ball_and_taper, neuron_ref) {
- using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
+template <typename Policy>
+using lowered_cell_div = fvm::fvm_multicell<double, cell_local_size_type, Policy>;
- cell c = make_cell_ball_and_taper();
+TEST(ball_and_squiggle, neuron_ref) {
+ cell c = make_cell_ball_and_squiggle();
add_common_voltage_probes(c);
float sample_dt = 0.025f;
sampler_info samplers[] = {
{"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
- {"taper.mid", {0u, 1u}, simple_sampler(sample_dt)},
- {"taper.end", {0u, 2u}, simple_sampler(sample_dt)}
+ {"dend.mid", {0u, 1u}, simple_sampler(sample_dt)},
+ {"dend.end", {0u, 2u}, simple_sampler(sample_dt)}
};
- run_ncomp_convergence_test<lowered_cell>(
- "ball_and_taper",
- "neuron_ball_and_taper.json",
+ run_ncomp_convergence_test<lowered_cell_div<div_compartment_sampler>>(
+ "ball_and_squiggle_sampler",
+ "neuron_ball_and_squiggle.json",
c,
samplers);
-}
+ run_ncomp_convergence_test<lowered_cell_div<div_compartment_integrator>>(
+ "ball_and_squiggle_integrator",
+ "neuron_ball_and_squiggle.json",
+ c,
+ samplers);
+}
diff --git a/tests/validation/validate_compartment_policy.cpp b/tests/validation/validate_compartment_policy.cpp
new file mode 100644
index 00000000..6ee21327
--- /dev/null
+++ b/tests/validation/validate_compartment_policy.cpp
@@ -0,0 +1,95 @@
+#include <fstream>
+#include <utility>
+
+#include <json/json.hpp>
+
+#include <common_types.hpp>
+#include <cell.hpp>
+#include <fvm_multicell.hpp>
+#include <model.hpp>
+#include <recipe.hpp>
+#include <simple_sampler.hpp>
+#include <util/rangeutil.hpp>
+
+#include "gtest.h"
+
+#include "../test_common_cells.hpp"
+#include "../test_util.hpp"
+#include "trace_analysis.hpp"
+#include "validation_data.hpp"
+
+using namespace nest::mc;
+
+/*
+ * Expect dendtrites composed of a simple frustrum to give
+ * essentially identical results no matter the compartment
+ * division policy.
+ */
+
+template <typename CompPolicy>
+std::vector<trace_data> run_model(const cell& c, float sample_dt, float t_end, float dt) {
+ model<fvm::fvm_multicell<double, cell_local_size_type, div_compartment_by_ends>> m{singleton_recipe(c)};
+
+ const auto& probes = m.probes();
+ std::size_t n_probes = probes.size();
+ std::vector<simple_sampler> samplers(n_probes, sample_dt);
+
+ for (unsigned i = 0; i<n_probes; ++i) {
+ m.attach_sampler(probes[i].id, samplers[i].sampler<>());
+ }
+
+ m.run(t_end, dt);
+ std::vector<trace_data> traces;
+ for (auto& s: samplers) {
+ traces.push_back(std::move(s.trace));
+ }
+ return traces;
+}
+
+
+void run_test(cell&& c) {
+ add_common_voltage_probes(c);
+
+ float sample_dt = .025;
+ float t_end = 100;
+ float dt = 0.001;
+
+ auto traces_by_ends = run_model<div_compartment_by_ends>(c, sample_dt, t_end, dt);
+ auto traces_sampler = run_model<div_compartment_sampler>(c, sample_dt, t_end, dt);
+ auto traces_integrator = run_model<div_compartment_integrator>(c, sample_dt, t_end, dt);
+
+ auto n_trace = traces_by_ends.size();
+ ASSERT_GT(n_trace, 0);
+ ASSERT_EQ(n_trace, traces_sampler.size());
+ ASSERT_EQ(n_trace, traces_integrator.size());
+
+ for (unsigned i = 0; i<n_trace; ++i) {
+ auto& t1 = traces_by_ends[i];
+ auto& t2 = traces_sampler[i];
+ auto& t3 = traces_integrator[i];
+
+ // expect all traces to be (close to) the same
+ double epsilon = 1e-6;
+ double tol = epsilon*util::max_value(
+ util::transform_view(values(t1), [](double x) { return std::abs(x); }));
+ EXPECT_GE(tol, linf_distance(t1, t2));
+ EXPECT_GE(tol, linf_distance(t2, t3));
+ EXPECT_GE(tol, linf_distance(t3, t1));
+ }
+}
+
+TEST(compartment_policy, validate_ball_and_stick) {
+ SCOPED_TRACE("ball_and_stick");
+ run_test(make_cell_ball_and_stick());
+}
+
+TEST(compartment_policy, validate_ball_and_3stick) {
+ SCOPED_TRACE("ball_and_3stick");
+ run_test(make_cell_ball_and_3stick());
+}
+
+TEST(compartment_policy, validate_ball_and_taper) {
+ SCOPED_TRACE("ball_and_taper");
+ run_test(make_cell_ball_and_taper());
+}
+
diff --git a/tests/validation/validate_soma.cpp b/tests/validation/validate_soma.cpp
index 9f341a50..4f619a9f 100644
--- a/tests/validation/validate_soma.cpp
+++ b/tests/validation/validate_soma.cpp
@@ -41,10 +41,21 @@ TEST(soma, numeric_ref) {
R.load_reference_data("numeric_soma.json");
float t_end = 100.f;
- for (auto dt: {0.05f, 0.02f, 0.01f, 0.005f, 0.001f}) {
- m.reset();
- R.run(m, dt, t_end, dt);
+
+ // use dt = 0.05, 0.025, 0.01, 0.005, 0.0025, ...
+ double max_oo_dt = std::round(1.0/g_trace_io.min_dt());
+ for (double base = 100; ; base *= 10) {
+ for (double multiple: {5., 2.5, 1.}) {
+ double oo_dt = base/multiple;
+ if (oo_dt>max_oo_dt) goto end;
+
+ m.reset();
+ float dt = float(1./oo_dt);
+ R.run(m, dt, t_end, dt);
+ }
}
+end:
+
R.report();
R.assert_all_convergence();
}
diff --git a/tests/validation/validate_synapses.cpp b/tests/validation/validate_synapses.cpp
index 9460d0c3..11b050e2 100644
--- a/tests/validation/validate_synapses.cpp
+++ b/tests/validation/validate_synapses.cpp
@@ -44,6 +44,9 @@ void run_synapse_test(
{{0u, 0u}, 40.0, 0.04}
};
+ // exclude points of discontinuity from linf analysis
+ std::vector<float> exclude = {10.f, 20.f, 40.f};
+
float sample_dt = 0.025f;
sampler_info samplers[] = {
{"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
@@ -59,7 +62,7 @@ void run_synapse_test(
model<lowered_cell> m(singleton_recipe{c});
m.group(0).enqueue_events(synthetic_events);
- R.run(m, ncomp, t_end, dt);
+ R.run(m, ncomp, t_end, dt, exclude);
}
R.report();
R.assert_all_convergence();
diff --git a/tests/validation/validation_data.hpp b/tests/validation/validation_data.hpp
index 254d6eee..f6a4b821 100644
--- a/tests/validation/validation_data.hpp
+++ b/tests/validation/validation_data.hpp
@@ -52,6 +52,9 @@ public:
void set_max_ncomp(int ncomp) { max_ncomp_ = ncomp; }
int max_ncomp() const { return max_ncomp_; }
+ void set_min_dt(float dt) { min_dt_ = dt; }
+ float min_dt() const { return min_dt_; }
+
void set_datadir(const util::path& dir) { datadir_ = dir; }
void set_output(const util::path& file) {
@@ -74,7 +77,8 @@ private:
std::ofstream out_;
nlohmann::json jtraces_ = nlohmann::json::array();
bool verbose_flag_ = false;
- int max_ncomp_ = 1000;
+ int max_ncomp_ = 100;
+ float min_dt_ = 0.001f;
};
extern trace_io g_trace_io;
diff --git a/validation/ref/neuron/ball_and_3stick.py b/validation/ref/neuron/ball_and_3stick.py
index a1c31e7a..a2613067 100644
--- a/validation/ref/neuron/ball_and_3stick.py
+++ b/validation/ref/neuron/ball_and_3stick.py
@@ -19,7 +19,6 @@ model.add_iclamp(5, 80, 0.45, to='dend2')
model.add_iclamp(40, 10, -0.2, to='dend3')
simdur = 100.0
-dt = 0.001
data = V.run_nrn_sim(simdur, report_dt=10, model='ball_and_3stick')
print json.dumps(data)
diff --git a/validation/ref/neuron/ball_and_squiggle.py b/validation/ref/neuron/ball_and_squiggle.py
index 42b2916a..4d4e0222 100644
--- a/validation/ref/neuron/ball_and_squiggle.py
+++ b/validation/ref/neuron/ball_and_squiggle.py
@@ -21,7 +21,6 @@ model.add_dendrite('dend', geom)
model.add_iclamp(5, 80, 0.3, to='dend')
simdur = 100.0
-dt = 0.001
data = V.run_nrn_sim(simdur, report_dt=10, model='ball_and_squiggle')
print json.dumps(data)
diff --git a/validation/ref/neuron/generate_validation.sh b/validation/ref/neuron/generate_validation.sh
index 81c2ecfe..d19893dd 100755
--- a/validation/ref/neuron/generate_validation.sh
+++ b/validation/ref/neuron/generate_validation.sh
@@ -24,7 +24,7 @@ if [ ! -d "$destdir" ]; then
usage
fi
-for v in soma ball_and_stick ball_and_taper ball_and_3stick simple_exp_synapse simple_exp2_synapse; do
+for v in soma ball_and_stick ball_and_squiggle ball_and_taper ball_and_3stick simple_exp_synapse simple_exp2_synapse; do
(cd "$scriptdir"; nrniv -nobanner -python ./$v.py) > "$destdir/neuron_$v.json" &
done
wait
diff --git a/validation/ref/neuron/nrn_validation.py b/validation/ref/neuron/nrn_validation.py
index a5ffa2d2..961b2771 100644
--- a/validation/ref/neuron/nrn_validation.py
+++ b/validation/ref/neuron/nrn_validation.py
@@ -37,7 +37,9 @@ default_model_parameters = {
'tau': 2.0, # Exponential synapse time constant
'tau1': 0.5, # Exp2 synapse tau1
'tau2': 2.0, # Exp2 synapse tau2
- 'ncomp': 1001 # Number of compartments (nseg) in dendrites
+ 'ncomp': 1001, # Number of compartments (nseg) in dendrites
+ 'dt': 0.0, # (Simulation parameter) default dt, 0 => use cvode adaptive
+ 'abstol': 1e-6 # (Simulation parameter) abstol for cvode if used
}
def override_defaults_from_args(args=sys.argv):
@@ -160,7 +162,10 @@ class VModel:
# Samples at cable mid- and end-points taken every `sample_dt`;
# Voltage on all compartments per section reported every `report_dt`.
-def run_nrn_sim(tend, sample_dt=0.025, report_t=None, report_dt=None, dt=0.001, **meta):
+def run_nrn_sim(tend, sample_dt=0.025, report_t=None, report_dt=None, dt=None, **meta):
+ if dt is None:
+ dt = default_model_parameters['dt']
+
# Instrument mid-point and ends of each section for traces.
vtraces = []
vtrace_t_hoc = h.Vector()
@@ -204,7 +209,18 @@ def run_nrn_sim(tend, sample_dt=0.025, report_t=None, report_dt=None, dt=0.001,
vreports.append((s.name(), s.L, s.nseg, ps, vs))
# Run sim
- h.dt = dt
+ if dt==0:
+ # Use CVODE instead
+ h.cvode.active(1)
+ abstol = default_model_parameters['abstol']
+ h.cvode.atol(abstol)
+ common_meta = { 'dt': 0, 'cvode': True, 'abstol': abstol }
+ else:
+ h.dt = dt
+ h.steps_per_ms = 1/dt # or else NEURON might noisily fudge dt
+ common_meta = { 'dt': dt, 'cvode': False }
+
+ h.secondorder = 2
h.tstop = tend
h.run()
@@ -212,7 +228,7 @@ def run_nrn_sim(tend, sample_dt=0.025, report_t=None, report_dt=None, dt=0.001,
traces = []
vtrace_t = list(vtrace_t_hoc)
- traces.append(combine(meta, common_ncomp, {
+ traces.append(combine(common_meta, meta, common_ncomp, {
'name': 'membrane voltage',
'sim': 'neuron',
'units': 'mV',
@@ -225,7 +241,10 @@ def run_nrn_sim(tend, sample_dt=0.025, report_t=None, report_dt=None, dt=0.001,
obs = np.column_stack([np.array(v) for v in vs])
xs = [length*p for p in ps]
for i, t in enumerate(report_t):
- traces.append(combine(meta, {
+ if i>=obs.shape[0]:
+ break
+
+ traces.append(combine(common_meta, meta, {
'name': 'membrane voltage',
'sim': 'neuron',
'units': {'x': 'µm', name: 'mV'},
diff --git a/validation/ref/neuron/soma.py b/validation/ref/neuron/soma.py
index 706af0f5..f1380032 100644
--- a/validation/ref/neuron/soma.py
+++ b/validation/ref/neuron/soma.py
@@ -11,8 +11,6 @@ model = V.VModel()
model.add_soma(18.8, Ra=100)
model.add_iclamp(10, 100, 0.1)
-# NB: this doesn't seem to have converged with
-# the default dt of 0.001.
data = V.run_nrn_sim(100, report_dt=None, model='soma')
print json.dumps(data)
V.nrn_stop()
--
GitLab