diff --git a/scripts/tsplot b/scripts/tsplot
index 8e89cbda96451680bfc9de8087ba558c0282defc..4280a9dc902cfa738e5eaf8bb7043008d96baae0 100755
--- a/scripts/tsplot
+++ b/scripts/tsplot
@@ -44,6 +44,7 @@ def parse_clargs():
help='plot series with same KEYs on the same axes')
P.add_argument('-s', '--select', metavar='EXPR,...', dest='select',
type=lambda s: s.split(','),
+ action='append',
help='select only series matching EXPR')
P.add_argument('-o', '--output', metavar='FILE', dest='outfile',
help='save plot to file FILE')
@@ -212,7 +213,7 @@ def read_json_timeseries(j, axis='time', select=[]):
meta = dict(j.items() + {'label': key, 'data': None, 'units': units(key)}.items())
del meta['data']
- if all([run_select(sel, meta) for sel in select]):
+ if not select or any([all([run_select(s, meta) for s in term]) for term in select]):
ts_list.append(TimeSeries(times, jdata[key], **meta))
return ts_list
@@ -422,7 +423,7 @@ args = parse_clargs()
tss = []
axis = args.axis if args.axis else 'time'
for filename in args.inputs:
- select = args.select if args.select else []
+ select = args.select
with open(filename) as f:
j = json.load(f)
tss.extend(read_json_timeseries(j, axis, select))
diff --git a/tests/test_common_cells.hpp b/tests/test_common_cells.hpp
index 41e2fea727b9e217899a0f8717b7085b7fc0f84a..b74a18c384a52301c45585d978c380be179c1e01 100644
--- a/tests/test_common_cells.hpp
+++ b/tests/test_common_cells.hpp
@@ -13,6 +13,7 @@ namespace mc {
* diameter: 18.8 µm
* mechanisms: HH (default params)
* bulk resistivitiy: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
*
* Stimuli:
* soma centre, t=[10 ms, 110 ms), 0.1 nA
@@ -37,6 +38,7 @@ inline cell make_cell_soma_only(bool with_stim = true) {
*
* Common properties:
* bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
*
* Soma:
* mechanisms: HH (default params)
@@ -80,6 +82,7 @@ inline cell make_cell_ball_and_stick(bool with_stim = true) {
*
* Common properties:
* bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
*
* Soma:
* mechanisms: HH (default params)
@@ -126,6 +129,7 @@ inline cell make_cell_ball_and_taper(bool with_stim = true) {
*
* Common properties:
* bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
*
* Soma:
* mechanisms: HH (default params)
@@ -173,11 +177,12 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
*
* Common properties:
* mechanisms: passive
- * membrane conductance: 0.000025 S·cm¯² ( = 1/(4Ω·m²) )
+ * membrane conductance: 0.000025 S/cm² ( = 1/(4Ω·m²) )
* membrane reversal potential: -65 mV (default)
* diameter: 1 µm
* length: 1000 µm
* bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
* compartments: 4
*
* Stimulus:
@@ -185,21 +190,43 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
*
* Note: zero-volume soma added with same mechanisms, as
* work-around for some existing fvm modelling issues.
+ *
+ * TODO: Set the correct values when parameters are generally
+ * settable!
+ *
+ * We can't currently change leak parameters
+ * from defaults, so we scale other electrical parameters
+ * proportionally.
*/
inline cell make_cell_simple_cable(bool with_stim = true) {
cell c;
c.add_soma(0);
- c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 1000);
+
+ double r_L = 100;
+ double c_m = 0.01;
+ double gbar = 0.000025;
+ double I = 0.1;
+
+ // fudge factor! can't change passive membrane
+ // conductance from gbar0 = 0.001
+
+ double gbar0 = 0.001;
+ double f = gbar/gbar0;
+
+ // scale everything else
+ r_L *= f;
+ c_m /= f;
+ I /= f;
for (auto& seg: c.segments()) {
seg->add_mechanism(pas_parameters());
- seg->mechanism("membrane").set("r_L", 100);
-
- // parameter support not there yet, so leave
- // as default and modify reference data for now.
- // seg->mechanism("pas").set("g", 0.000025);
+ seg->mechanism("membrane").set("r_L", r_L);
+ seg->mechanism("membrane").set("c_m", c_m);
+ // seg->mechanism("pas").set("g", gbar);
+
if (seg->is_dendrite()) {
seg->set_compartments(4);
}
@@ -208,7 +235,7 @@ inline cell make_cell_simple_cable(bool with_stim = true) {
if (with_stim) {
// stimulus in the middle of our zero-volume 'soma'
// corresponds to proximal end of cable.
- c.add_stimulus({0,0.5}, {0., HUGE_VAL, 0.1});
+ c.add_stimulus({0,0.5}, {0., HUGE_VAL, I});
}
return c;
}
diff --git a/tests/validation/trace_analysis.cpp b/tests/validation/trace_analysis.cpp
index d641e34bcbad20ea1c64aece9d47f9c17a4856f0..a754cfe228b8893a6d557eb40092cbde8de4dfa5 100644
--- a/tests/validation/trace_analysis.cpp
+++ b/tests/validation/trace_analysis.cpp
@@ -83,7 +83,7 @@ util::optional<trace_peak> peak_delta(const trace_data& a, const trace_data& b)
auto p = local_maxima(a);
auto q = local_maxima(b);
- if (p.size()!=q.size()) return util::nothing;
+ if (p.size()!=q.size() || p.empty()) return util::nothing;
auto max_delta = p[0]-q[0];
diff --git a/tests/validation/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp
index 7927a748edc09830b151b68350ce6b35cfb3615a..e01eb4cf9c6f95a2631e1bbf370762967df71487 100644
--- a/tests/validation/validate_ball_and_stick.cpp
+++ b/tests/validation/validate_ball_and_stick.cpp
@@ -33,10 +33,11 @@ template <
>
void run_ncomp_convergence_test(
const char* model_name,
- const char* refdata_path,
+ const char* ref_data_path,
const cell& c,
SamplerInfoSeq& samplers,
- float t_end=100.f)
+ float t_end=100.f,
+ float dt=0.001)
{
using nlohmann::json;
@@ -52,7 +53,7 @@ void run_ncomp_convergence_test(
bool run_validation = false;
std::map<std::string, trace_data> ref_data;
try {
- ref_data = V.load_traces(refdata_path);
+ ref_data = V.load_traces(ref_data_path);
run_validation = std::all_of(keys.begin(), keys.end(),
[&](const char* key) { return ref_data.count(key)>0; });
@@ -60,7 +61,7 @@ void run_ncomp_convergence_test(
EXPECT_TRUE(run_validation);
}
catch (std::runtime_error&) {
- ADD_FAILURE() << "failure loading reference data: " << refdata_path;
+ ADD_FAILURE() << "failure loading reference data: " << ref_data_path;
}
std::map<std::string, std::vector<conv_entry<int>>> conv_results;
@@ -79,7 +80,7 @@ void run_ncomp_convergence_test(
m.attach_sampler(se.probe, se.sampler.template sampler<>());
}
- m.run(100, 0.001);
+ m.run(t_end, dt);
for (auto& se: samplers) {
std::string key = se.label;
@@ -202,279 +203,7 @@ TEST(rallpack1, numeric_ref) {
"rallpack1",
"numeric_rallpack1.json",
c,
- samplers);
-}
-
-#if 0
-TEST(ball_and_stick, neuron_ref) {
- // compare voltages against reference data produced from
- // nrn/ball_and_stick.py
-
- using namespace nlohmann;
-
- using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
- auto& V = g_trace_io;
-
- bool verbose = V.verbose();
- int max_ncomp = V.max_ncomp();
-
- // load validation data
- auto ref_data = V.load_traces("neuron_ball_and_stick.json");
- bool run_validation =
- ref_data.count("soma.mid") &&
- ref_data.count("dend.mid") &&
- ref_data.count("dend.end");
-
- EXPECT_TRUE(run_validation);
-
- // generate test data
- cell c = make_cell_ball_and_stick();
- add_common_voltage_probes(c);
-
- float sample_dt = .025;
- std::pair<const char *, simple_sampler> samplers[] = {
- {"soma.mid", simple_sampler(sample_dt)},
- {"dend.mid", simple_sampler(sample_dt)},
- {"dend.end", simple_sampler(sample_dt)}
- };
-
- std::map<std::string, std::vector<conv_entry<int>>> conv_results;
-
- for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- for (auto& se: samplers) {
- se.second.reset();
- }
- c.cable(1)->set_compartments(ncomp);
- model<lowered_cell> m(singleton_recipe{c});
-
- // the ball-and-stick-cell (should) have three voltage probes:
- // centre of soma, centre of dendrite, end of dendrite.
-
- m.attach_sampler({0u, 0u}, samplers[0].second.sampler<>());
- m.attach_sampler({0u, 1u}, samplers[1].second.sampler<>());
- m.attach_sampler({0u, 2u}, samplers[2].second.sampler<>());
-
- m.run(100, 0.001);
-
- for (auto& se: samplers) {
- std::string key = se.first;
- const simple_sampler& s = se.second;
-
- // save trace
- json meta = {
- {"name", "membrane voltage"},
- {"model", "ball_and_stick"},
- {"sim", "nestmc"},
- {"ncomp", ncomp},
- {"units", "mV"}};
-
- V.save_trace(key, s.trace, meta);
-
- // compute metrics
- if (run_validation) {
- double linf = linf_distance(s.trace, ref_data[key]);
- auto pd = peak_delta(s.trace, ref_data[key]);
-
- conv_results[key].push_back({key, ncomp, linf, pd});
- }
- }
- }
-
- if (verbose && run_validation) {
- report_conv_table(std::cout, conv_results, "ncomp");
- }
-
- for (auto key: util::transform_view(samplers, util::first)) {
- SCOPED_TRACE(key);
-
- const auto& results = conv_results[key];
- assert_convergence(results);
- }
-}
-
-TEST(ball_and_taper, neuron_ref) {
- // compare voltages against reference data produced from
- // nrn/ball_and_taper.py
-
- using namespace nlohmann;
-
- using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
- auto& V = g_trace_io;
-
- bool verbose = V.verbose();
- int max_ncomp = V.max_ncomp();
-
- // load validation data
- auto ref_data = V.load_traces("neuron_ball_and_taper.json");
- bool run_validation =
- ref_data.count("soma.mid") &&
- ref_data.count("taper.mid") &&
- ref_data.count("taper.end");
-
- EXPECT_TRUE(run_validation);
-
- // generate test data
- cell c = make_cell_ball_and_taper();
- add_common_voltage_probes(c);
-
- float sample_dt = .025;
- std::pair<const char *, simple_sampler> samplers[] = {
- {"soma.mid", simple_sampler(sample_dt)},
- {"taper.mid", simple_sampler(sample_dt)},
- {"taper.end", simple_sampler(sample_dt)}
- };
-
- std::map<std::string, std::vector<conv_entry<int>>> conv_results;
-
- for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- for (auto& se: samplers) {
- se.second.reset();
- }
- c.cable(1)->set_compartments(ncomp);
- model<lowered_cell> m(singleton_recipe{c});
-
- // the ball-and-stick-cell (should) have three voltage probes:
- // centre of soma, centre of dendrite, end of dendrite.
-
- m.attach_sampler({0u, 0u}, samplers[0].second.sampler<>());
- m.attach_sampler({0u, 1u}, samplers[1].second.sampler<>());
- m.attach_sampler({0u, 2u}, samplers[2].second.sampler<>());
-
- m.run(100, 0.001);
-
- for (auto& se: samplers) {
- std::string key = se.first;
- const simple_sampler& s = se.second;
-
- // save trace
- json meta = {
- {"name", "membrane voltage"},
- {"model", "ball_and_taper"},
- {"sim", "nestmc"},
- {"ncomp", ncomp},
- {"units", "mV"}};
-
- V.save_trace(key, s.trace, meta);
-
- // compute metrics
- if (run_validation) {
- double linf = linf_distance(s.trace, ref_data[key]);
- auto pd = peak_delta(s.trace, ref_data[key]);
-
- conv_results[key].push_back({key, ncomp, linf, pd});
- }
- }
- }
-
- if (verbose && run_validation) {
- report_conv_table(std::cout, conv_results, "ncomp");
- }
-
- for (auto key: util::transform_view(samplers, util::first)) {
- SCOPED_TRACE(key);
-
- const auto& results = conv_results[key];
- assert_convergence(results);
- }
-}
-
-
-TEST(ball_and_3stick, neuron_ref) {
- // compare voltages against reference data produced from
- // nrn/ball_and_3stick.py
-
- using namespace nlohmann;
-
- using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
- auto& V = g_trace_io;
-
- bool verbose = V.verbose();
- int max_ncomp = V.max_ncomp();
-
- // load validation data
- auto ref_data = V.load_traces("neuron_ball_and_3stick.json");
- bool run_validation =
- ref_data.count("soma.mid") &&
- ref_data.count("dend1.mid") &&
- ref_data.count("dend1.end") &&
- ref_data.count("dend2.mid") &&
- ref_data.count("dend2.end") &&
- ref_data.count("dend3.mid") &&
- ref_data.count("dend3.end");
-
-
- EXPECT_TRUE(run_validation);
-
- // generate test data
- cell c = make_cell_ball_and_3stick();
- add_common_voltage_probes(c);
-
- float sample_dt = .025;
- std::pair<const char *, simple_sampler> samplers[] = {
- {"soma.mid", simple_sampler(sample_dt)},
- {"dend1.mid", simple_sampler(sample_dt)},
- {"dend1.end", simple_sampler(sample_dt)},
- {"dend2.mid", simple_sampler(sample_dt)},
- {"dend2.end", simple_sampler(sample_dt)},
- {"dend3.mid", simple_sampler(sample_dt)},
- {"dend3.end", simple_sampler(sample_dt)}
- };
-
- std::map<std::string, std::vector<conv_entry<int>>> conv_results;
-
- for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- for (auto& se: samplers) {
- se.second.reset();
- }
- c.cable(1)->set_compartments(ncomp);
- c.cable(2)->set_compartments(ncomp);
- c.cable(3)->set_compartments(ncomp);
- model<lowered_cell> m(singleton_recipe{c});
-
- // the ball-and-3stick-cell (should) have seven voltage probes:
- // centre of soma, followed by centre of section, end of section
- // for each of the three dendrite sections.
-
- for (unsigned i = 0; i < util::size(samplers); ++i) {
- m.attach_sampler({0u, i}, samplers[i].second.sampler<>());
- }
-
- m.run(100, 0.001);
-
- for (auto& se: samplers) {
- std::string key = se.first;
- const simple_sampler& s = se.second;
-
- // save trace
- json meta = {
- {"name", "membrane voltage"},
- {"model", "ball_and_3stick"},
- {"sim", "nestmc"},
- {"ncomp", ncomp},
- {"units", "mV"}};
-
- V.save_trace(key, s.trace, meta);
-
- // compute metrics
- if (run_validation) {
- double linf = linf_distance(s.trace, ref_data[key]);
- auto pd = peak_delta(s.trace, ref_data[key]);
-
- conv_results[key].push_back({key, ncomp, linf, pd});
- }
- }
- }
-
- if (verbose && run_validation) {
- report_conv_table(std::cout, conv_results, "ncomp");
- }
-
- for (auto key: util::transform_view(samplers, util::first)) {
- SCOPED_TRACE(key);
-
- const auto& results = conv_results[key];
- assert_convergence(results);
- }
+ samplers,
+ 250.f);
}
-#endif
diff --git a/tests/validation/validate_soma.cpp b/tests/validation/validate_soma.cpp
index 3e848b65856691bef455a6ab0b76dd9c7c1dda1e..31eb784add73f5b675e8b2943e8faa453a0afb26 100644
--- a/tests/validation/validate_soma.cpp
+++ b/tests/validation/validate_soma.cpp
@@ -32,11 +32,21 @@ TEST(soma, numeric_ref) {
bool verbose = V.verbose();
// load validation data
- //auto ref_data = V.load_traces("neuron_soma.json");
- auto ref_data = V.load_traces("numeric_soma.json");
+
+ bool run_validation = false;
+ std::map<std::string, trace_data> ref_data;
const char* key = "soma.mid";
- bool run_validation = ref_data.count(key);
- EXPECT_TRUE(run_validation);
+
+ const char* ref_data_path = "numeric_soma.json";
+ try {
+ ref_data = V.load_traces(ref_data_path);
+ run_validation = ref_data.count(key);
+
+ EXPECT_TRUE(run_validation);
+ }
+ catch (std::runtime_error&) {
+ ADD_FAILURE() << "failure loading reference data: " << ref_data_path;
+ }
// generate test data
cell c = make_cell_soma_only();
diff --git a/validation/ref/numeric/CMakeLists.txt b/validation/ref/numeric/CMakeLists.txt
index 7aaf6b7ae3157ede5cf66eadf82d73f73cbb107b..68dc379216190d25ecf29aefb352a189d3913803 100644
--- a/validation/ref/numeric/CMakeLists.txt
+++ b/validation/ref/numeric/CMakeLists.txt
@@ -5,3 +5,7 @@ add_validation_data(
DEPENDS numeric_soma.jl HHChannels.jl
COMMAND julia numeric_soma.jl)
+add_validation_data(
+ OUTPUT numeric_rallpack1.json
+ DEPENDS numeric_rallpack1.jl PassiveCable.jl
+ COMMAND julia numeric_rallpack1.jl)
diff --git a/validation/ref/numeric/PassiveCable.jl b/validation/ref/numeric/PassiveCable.jl
index 76450490ce5c4f265403edd35f2bedad4ad7c3bb..6928b74c33948be7162885040cb89870ae245596 100644
--- a/validation/ref/numeric/PassiveCable.jl
+++ b/validation/ref/numeric/PassiveCable.jl
@@ -18,26 +18,29 @@ export cable_normalize, cable, rallpack1
#
# Return:
# g(x, t)
+#
+# TODO: verify correctness when L≠1
-function cable_normalized(x, t, L; tol=1e-8)
+function cable_normalized(x::Float64, t::Float64, L::Float64; tol=1e-8)
if t<=0
return 0.0
else
ginf = -cosh(L-x)/sinh(L)
sum = exp(-t/L)
+ Ltol = L*tol
for k = countfrom(1)
a = k*pi/L
- e = exp(-t*(1+a^2))
+ b = exp(-t*(1+a^2))
- sum += 2/L*e*cos(a*x)/(1+a^2)
- resid_ub = e/(L*a^3*t)
+ sum += 2*b*cos(a*x)/(1+a^2)
+ resid_ub = b/(a^3*t)
- if resid_ub<tol
+ if resid_ub<Ltol
break
end
end
- return ginf+sum
+ return ginf+sum/L;
end
end
diff --git a/validation/ref/numeric/numeric_rallpack1.jl b/validation/ref/numeric/numeric_rallpack1.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d0f5997e66310cada8b7258249768b19c3850012
--- /dev/null
+++ b/validation/ref/numeric/numeric_rallpack1.jl
@@ -0,0 +1,60 @@
+#!/usr/bin/env julia
+
+include("PassiveCable.jl")
+
+using JSON
+using SIUnits.ShortUnits
+using PassiveCable
+
+function run_cable(x_prop, ts)
+ # physical properties
+
+ # f is a fudge factor. rM needs to be the same
+ # the same as in nestmc, where we cannot yet set
+ # the membrane conductance parameter. Scaling
+ # other parameters proportionally, however,
+ # gives the same dynamics.
+
+ f = 0.1/4
+
+ diam = 1.0µm # cable diameter
+ L = 1.0mm # cable length
+ I = 0.1nA /f # current injection
+ rL = 1.0Ω*m *f # bulk resistivity
+ erev = -65.0mV # (passive) reversal potential
+ rM = 4Ω*m^2 *f # membrane resistivity
+ cM = 0.01F/m^2 /f # membrane specific capacitance
+
+ # convert to linear resistivity, length and time constants
+ area = pi*diam^2/4
+ r = rL/area
+
+ lambda = sqrt(diam/4 * rM/rL)
+ tau = cM*rM
+
+ # compute solutions
+ tol = 1e-8mV
+ return [cable(L*x_prop, t, L, lambda, tau, r, erev, -I, tol=tol) for t in ts]
+end
+
+function run_rallpack1(x_prop, ts)
+ return [rallpack1(0.001*x_prop, t/s)*V for t in ts]
+end
+# Generate traces at x=0, x=0.3L, x=L
+
+ts = collect(0s: 0.025ms: 250ms)
+trace = Dict(
+ :name => "membrane voltage",
+ :sim => "numeric",
+ :model => "rallpack1",
+ :units => "mV",
+ :data => Dict(
+ :time => map(t->t/ms, ts),
+ symbol("cable.x0.0") => map(v->v/mV, run_cable(0, ts)),
+ symbol("cable.x0.3") => map(v->v/mV, run_rallpack1(0.3, ts)),
+ symbol("cable.x1.0") => map(v->v/mV, run_cable(1.0, ts))
+ )
+)
+
+println(JSON.json([trace]))
+