diff --git a/tests/test_common_cells.hpp b/tests/test_common_cells.hpp
index 8186c5ed0c678ba9b3ccd06b130ea62d63a65e5d..41e2fea727b9e217899a0f8717b7085b7fc0f84a 100644
--- a/tests/test_common_cells.hpp
+++ b/tests/test_common_cells.hpp
@@ -196,6 +196,10 @@ inline cell make_cell_simple_cable(bool with_stim = true) {
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);
if (seg->is_dendrite()) {
seg->set_compartments(4);
}
diff --git a/tests/validation/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp
index 0a57fa840e52e70ea7dcaef4dc14c0f54ee86d93..7927a748edc09830b151b68350ce6b35cfb3615a 100644
--- a/tests/validation/validate_ball_and_stick.cpp
+++ b/tests/validation/validate_ball_and_stick.cpp
@@ -21,8 +21,191 @@
using namespace nest::mc;
-// TODO: further consolidate common code
+struct sampler_info {
+ const char* label;
+ cell_member_type probe;
+ simple_sampler sampler;
+};
+
+template <
+ typename lowered_cell,
+ typename SamplerInfoSeq
+>
+void run_ncomp_convergence_test(
+ const char* model_name,
+ const char* refdata_path,
+ const cell& c,
+ SamplerInfoSeq& samplers,
+ float t_end=100.f)
+{
+ using nlohmann::json;
+
+ SCOPED_TRACE(model_name);
+ auto& V = g_trace_io;
+ bool verbose = V.verbose();
+ int max_ncomp = V.max_ncomp();
+
+ auto keys = util::transform_view(samplers,
+ [](const sampler_info& se) { return se.label; });
+
+ bool run_validation = false;
+ std::map<std::string, trace_data> ref_data;
+ try {
+ ref_data = V.load_traces(refdata_path);
+
+ run_validation = std::all_of(keys.begin(), keys.end(),
+ [&](const char* key) { return ref_data.count(key)>0; });
+
+ EXPECT_TRUE(run_validation);
+ }
+ catch (std::runtime_error&) {
+ ADD_FAILURE() << "failure loading reference data: " << refdata_path;
+ }
+
+ std::map<std::string, std::vector<conv_entry<int>>> conv_results;
+
+ for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
+ for (auto& seg: c.segments()) {
+ if (!seg->is_soma()) {
+ seg->set_compartments(ncomp);
+ }
+ }
+ model<lowered_cell> m(singleton_recipe{c});
+
+ // reset samplers and attach to probe locations
+ for (auto& se: samplers) {
+ se.sampler.reset();
+ m.attach_sampler(se.probe, se.sampler.template sampler<>());
+ }
+
+ m.run(100, 0.001);
+
+ for (auto& se: samplers) {
+ std::string key = se.label;
+ const simple_sampler& s = se.sampler;
+
+ // save trace
+ json meta = {
+ {"name", "membrane voltage"},
+ {"model", model_name},
+ {"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: keys) {
+ SCOPED_TRACE(key);
+
+ const auto& results = conv_results[key];
+ assert_convergence(results);
+ }
+}
+
+TEST(ball_and_stick, neuron_ref) {
+ using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
+
+ cell c = make_cell_ball_and_stick();
+ add_common_voltage_probes(c);
+
+ float sample_dt = 0.025f;
+ sampler_info samplers[] = {
+ {"soma.mid", {0u, 0u}, 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_stick",
+ "neuron_ball_and_stick.json",
+ c,
+ samplers);
+}
+
+TEST(ball_and_taper, neuron_ref) {
+ using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
+
+ cell c = make_cell_ball_and_taper();
+ 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)}
+ };
+
+ run_ncomp_convergence_test<lowered_cell>(
+ "ball_and_taper",
+ "neuron_ball_and_taper.json",
+ c,
+ samplers);
+}
+
+TEST(ball_and_3stick, neuron_ref) {
+ using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
+
+ cell c = make_cell_ball_and_3stick();
+ add_common_voltage_probes(c);
+
+ float sample_dt = 0.025f;
+ sampler_info samplers[] = {
+ {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
+ {"dend1.mid", {0u, 1u}, simple_sampler(sample_dt)},
+ {"dend1.end", {0u, 2u}, simple_sampler(sample_dt)},
+ {"dend2.mid", {0u, 3u}, simple_sampler(sample_dt)},
+ {"dend2.end", {0u, 4u}, simple_sampler(sample_dt)},
+ {"dend3.mid", {0u, 5u}, simple_sampler(sample_dt)},
+ {"dend3.end", {0u, 6u}, simple_sampler(sample_dt)}
+ };
+
+ run_ncomp_convergence_test<lowered_cell>(
+ "ball_and_3stick",
+ "neuron_ball_and_3stick.json",
+ c,
+ samplers);
+}
+
+TEST(rallpack1, numeric_ref) {
+ using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
+
+ cell c = make_cell_simple_cable();
+
+ // three probes: left end, 30% along, right end.
+ c.add_probe({{1, 0.0}, probeKind::membrane_voltage});
+ c.add_probe({{1, 0.3}, probeKind::membrane_voltage});
+ c.add_probe({{1, 1.0}, probeKind::membrane_voltage});
+
+ float sample_dt = 0.025f;
+ sampler_info samplers[] = {
+ {"cable.x0.0", {0u, 0u}, simple_sampler(sample_dt)},
+ {"cable.x0.3", {0u, 1u}, simple_sampler(sample_dt)},
+ {"cable.x1.0", {0u, 2u}, simple_sampler(sample_dt)},
+ };
+
+ run_ncomp_convergence_test<lowered_cell>(
+ "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
@@ -293,4 +476,5 @@ TEST(ball_and_3stick, neuron_ref) {
assert_convergence(results);
}
}
+#endif
diff --git a/tests/validation/validate_soma.cpp b/tests/validation/validate_soma.cpp
index 26d919b07c10a5966734c0e12dd7a7296e55f139..3e848b65856691bef455a6ab0b76dd9c7c1dda1e 100644
--- a/tests/validation/validate_soma.cpp
+++ b/tests/validation/validate_soma.cpp
@@ -20,7 +20,7 @@
using namespace nest::mc;
-TEST(soma, neuron_ref) {
+TEST(soma, numeric_ref) {
// compare voltages against reference data produced from
// nrn/ball_and_taper.py
@@ -32,7 +32,8 @@ TEST(soma, neuron_ref) {
bool verbose = V.verbose();
// load validation data
- auto ref_data = V.load_traces("neuron_soma.json");
+ //auto ref_data = V.load_traces("neuron_soma.json");
+ auto ref_data = V.load_traces("numeric_soma.json");
const char* key = "soma.mid";
bool run_validation = ref_data.count(key);
EXPECT_TRUE(run_validation);
diff --git a/tests/validation/validation_data.hpp b/tests/validation/validation_data.hpp
index 70e99a52227b70f9772dd94370ec38fb7301e2d1..254d6eeed77cd1d99c3ba31a924a8db361e7ff21 100644
--- a/tests/validation/validation_data.hpp
+++ b/tests/validation/validation_data.hpp
@@ -70,7 +70,7 @@ public:
}
private:
- util::path datadir_ = DATADIR "/validation";
+ util::path datadir_ = DATADIR;
std::ofstream out_;
nlohmann::json jtraces_ = nlohmann::json::array();
bool verbose_flag_ = false;
diff --git a/validation/ref/neuron/nrn_validation.py b/validation/ref/neuron/nrn_validation.py
index 66563885b05ba9212bd017ea1c90b514c93c6ff2..a5ffa2d2e238148d68dee137dd92daa17ce89d89 100644
--- a/validation/ref/neuron/nrn_validation.py
+++ b/validation/ref/neuron/nrn_validation.py
@@ -33,7 +33,7 @@ default_model_parameters = {
'g_pas': 0.001, # Passive conductance in S/cm^2
'e_pas': -65.0, # Leak reversal potential in mV
'Ra': 100.0, # Intracellular resistivity in Ω·cm
- 'cm': 1.0, # Membrane areal capacitance in µF/cm2
+ 'cm': 1.0, # Membrane areal capacitance in µF/cm^2
'tau': 2.0, # Exponential synapse time constant
'tau1': 0.5, # Exp2 synapse tau1
'tau2': 2.0, # Exp2 synapse tau2
diff --git a/validation/ref/neuron/soma.py b/validation/ref/neuron/soma.py
index f00d6bf012263d45a57209e49a7d094ba07f38d5..706af0f559a3172284bb45d9ca503e47a9f63cc2 100644
--- a/validation/ref/neuron/soma.py
+++ b/validation/ref/neuron/soma.py
@@ -6,9 +6,6 @@ import nrn_validation as V
V.override_defaults_from_args()
-# dendrite geometry: all 100 µm long, 1 µm diameter.
-geom = [(0,1), (100, 1)]
-
model = V.VModel()
model.add_soma(18.8, Ra=100)
diff --git a/validation/ref/numeric/CMakeLists.txt b/validation/ref/numeric/CMakeLists.txt
index a697439f0f8681efa9f0a0016d023e1897e7684c..7aaf6b7ae3157ede5cf66eadf82d73f73cbb107b 100644
--- a/validation/ref/numeric/CMakeLists.txt
+++ b/validation/ref/numeric/CMakeLists.txt
@@ -1,4 +1,7 @@
# note: function add_validation_data defined in validation/CMakeLists.txt
-# placeholder
-add_validation_data(OUTPUT foo.out DEPENDS foo.sh COMMAND bash foo.sh)
+add_validation_data(
+ OUTPUT numeric_soma.json
+ DEPENDS numeric_soma.jl HHChannels.jl
+ COMMAND julia numeric_soma.jl)
+
diff --git a/validation/ref/numeric/HHChannels.jl b/validation/ref/numeric/HHChannels.jl
new file mode 100644
index 0000000000000000000000000000000000000000..8ce1bd993fb8aa2f86f68dc25e885a4e6e45c832
--- /dev/null
+++ b/validation/ref/numeric/HHChannels.jl
@@ -0,0 +1,143 @@
+module HHChannels
+
+export Stim, run_hh
+
+using Sundials
+using SIUnits.ShortUnits
+
+immutable HHParam
+ c_m # membrane spacific capacitance
+ gnabar # Na channel cross-membrane conductivity
+ gkbar # K channel cross-membrane conductivity
+ gl # Leak conductivity
+ ena # Na channel reversal potential
+ ek # K channel reversal potential
+ el # Leak reversal potential
+ q10 # temperature dependent rate coefficient
+ # (= 3^((T-T₀)/10K) with T₀ = 6.3 °C)
+
+ # constructor with default values, corresponding
+ # to a resting potential of -65 mV and temperature 6.3 °C
+ HHParam(;
+ #c_m = 0.01nF*m^-2,
+ c_m = 0.01F*m^-2,
+ gnabar = .12S*cm^-2,
+ ena = 115.0mV + -65.0mV,
+ gkbar = .036S*cm^-2,
+ ek = -12.0mV + -65.0mV,
+ gl = .0003S*cm^-2,
+ el = -54.3mV,
+ q10 = 1
+ ) = new(c_m, gnabar, gkbar, gl, ena, ek, el, q10)
+
+end
+
+immutable Stim
+ t0 # start time of stimulus
+ t1 # stop time of stimulus
+ j # stimulus current density
+
+ Stim() = new(0s, 0s, 0A/m^2)
+ Stim(t0, t1, j) = new(t0, t1, j)
+end
+
+vtrap(x,y) = x/(exp(x/y) - 1.0)
+
+# "m" sodium activation system
+function m_lims(v, q10)
+ alpha = .1mV^-1 * vtrap(-(v+40mV),10mV)
+ beta = 4 * exp(-(v+65mV)/18mV)
+ sum = alpha + beta
+ mtau = 1ms / (q10*sum)
+ minf = alpha/sum
+ return mtau, minf
+end
+
+# "h" sodium inactivation system
+function h_lims(v, q10)
+ alpha = 0.07*exp(-(v+65mV)/20mV)
+ beta = 1 / (exp(-(v+35mV)/10mV) + 1)
+ sum = alpha + beta
+ htau = 1ms / (q10*sum)
+ hinf = alpha/sum
+ return htau, hinf
+end
+
+# "n" potassium activation system
+function n_lims(v, q10)
+ alpha = .01mV^-1 * vtrap(-(v+55mV),10mV)
+ beta = .125*exp(-(v+65mV)/80mV)
+ sum = alpha + beta
+ ntau = 1ms / (q10*sum)
+ ninf = alpha/sum
+ return ntau, ninf
+end
+
+# Choose initial conditions for the system such that the gating variables
+# are at steady state for the user-specified voltage v
+function initial_conditions(v, q10)
+ mtau, minf = m_lims(v, q10)
+ htau, hinf = h_lims(v, q10)
+ ntau, ninf = n_lims(v, q10)
+
+ return (v, minf, hinf, ninf)
+end
+
+# Given time t and state (v, m, h, n),
+# return (vdot, mdot, hdot, ndot)
+function f(t, state; p=HHParam(), stim=Stim())
+ v, m, h, n = state
+
+ # calculate current density due to ion channels
+ gna = p.gnabar*m*m*m*h
+ gk = p.gkbar*n*n*n*n
+
+
+ ina = gna*(v - p.ena)
+ ik = gk*(v - p.ek)
+ il = p.gl*(v - p.el)
+
+ itot = ik + ina + il
+
+ # calculate current density due to stimulus
+ if t>=stim.t0 && t<stim.t1
+ itot -= stim.j
+ end
+
+ # calculate the voltage dependent rates for the gating variables
+ mtau, minf = m_lims(v, p.q10)
+ htau, hinf = h_lims(v, p.q10)
+ ntau, ninf = n_lims(v, p.q10)
+
+ return (-itot/p.c_m, (minf-m)/mtau, (hinf-h)/htau, (ninf-n)/ntau)
+end
+
+function run_hh(t_end; v0=-65mV, stim=Stim(), param=HHParam(), sample_dt=0.01ms)
+ v_scale = 1V
+ t_scale = 1s
+
+ v0, m0, h0, n0 = initial_conditions(v0, param.q10)
+ y0 = [ v0/v_scale, m0, h0, n0 ]
+
+ samples = collect(0s: sample_dt: t_end)
+
+ fbis(t, y, ydot) = begin
+ vdot, mdot, hdot, ndot =
+ f(t*t_scale, (y[1]*v_scale, y[2], y[3], y[4]), stim=stim, p=param)
+
+ ydot[1], ydot[2], ydot[3], ydot[4] =
+ vdot*t_scale/v_scale, mdot*t_scale, hdot*t_scale, ndot*t_scale
+
+ return Sundials.CV_SUCCESS
+ end
+
+ # Ideally would run with vector absolute tolerance to account for v_scale,
+ # but this would prevent us using the nice cvode wrapper.
+
+ res = Sundials.cvode(fbis, y0, map(t->t/t_scale, samples), abstol=1e-6, reltol=5e-10)
+
+ # Use map here because of issues with type deduction with arrays and SIUnits.
+ return samples, map(v->v*v_scale, res[:, 1])
+end
+
+end # module HHChannels
diff --git a/validation/ref/numeric/hh_soma.jl b/validation/ref/numeric/hh_soma.jl
deleted file mode 100644
index 1db5813f0349b45f90ac3db1dcdef18352716058..0000000000000000000000000000000000000000
--- a/validation/ref/numeric/hh_soma.jl
+++ /dev/null
@@ -1,138 +0,0 @@
-using Sundials
-using SIUnits.ShortUnits
-
-c_m = 0.01nF*m^-2
-
-radius = 18.8μm / 2
-surface_area = 4 * pi * radius * radius
-
-gnabar = .12S*cm^-2
-gkbar = .036S*cm^-2
-gl = .0003S*cm^-2
-el = -54.3mV
-#celsius= 6.3
-#q10 = 3^((celsius - 6.3)/10)
-q10 = 1
-
-# define the resting potential for the membrane
-vrest = -65.0mV
-
-# define the resting potentials for ion species
-ena = 115.0mV+vrest
-ek = -12.0mV+vrest
-eca = 12.5mV*log(2.0/5e-5)
-
-vtrap(x,y) = x/(exp(x/y) - 1.0)
-
-function print()
- println("q10 ", q10)
- println("vrest ", vrest)
- println("ena ", ena)
- println("ek ", ek)
- println("eca ", eca)
-end
-
-#"m" sodium activation system
-function m_lims(v)
- alpha = .1mV^-1 * vtrap(-(v+40mV),10mV)
- beta = 4 * exp(-(v+65mV)/18mV)
- sum = alpha + beta
- mtau = 1ms / (q10*sum)
- minf = alpha/sum
- return mtau, minf
-end
-
-#"h" sodium inactivation system
-function h_lims(v)
- alpha = 0.07*exp(-(v+65mV)/20mV)
- beta = 1 / (exp(-(v+35mV)/10mV) + 1)
- sum = alpha + beta
- htau = 1ms / (q10*sum)
- hinf = alpha/sum
- return htau, hinf
-end
-
-#"n" potassium activation system
-function n_lims(v)
- alpha = .01mV^-1 * vtrap(-(v+55mV),10mV)
- beta = .125*exp(-(v+65mV)/80mV)
- sum = alpha + beta
- ntau = 1ms / (q10*sum)
- ninf = alpha/sum
- return ntau, ninf
-end
-
-# v = y[1] V
-# m = y[2]
-# h = y[3]
-# n = y[4]
-
-# dv/dt = ydot[1] V/s
-# dm/dt = ydot[2] /s
-# dh/dt = ydot[3] /s
-# dn/dt = ydot[4] /s
-
-# choose initial conditions for the system such that the gating variables
-# are at steady state for the user-specified voltage v
-function initial_conditions(v)
- mtau, minf = m_lims(v)
- htau, hinf = h_lims(v)
- ntau, ninf = n_lims(v)
-
- return [v/V, minf, hinf, ninf]
-end
-
-# calculate the lhs of the ODE system
-function f(t, y, ydot)
- # copy variables into helper variable
- v = y[1]V
- m, h, n = y[2], y[3], y[4]
-
- # calculate current due to ion channels
- gna = gnabar*m*m*m*h
- gk = gkbar*n*n*n*n
- ina = gna*(v - ena)
- ik = gk*(v - ek)
- il = gl*(v - el)
- imembrane = ik + ina + il
-
- # calculate current due to stimulus
- #c.add_stimulus({0,0.5}, {10., 100., 0.1});
- ielectrode = 0.0nA / surface_area
- time = t*ms
- if time>=10ms && time<100ms
- ielectrode = 0.1nA / surface_area
- end
-
- # calculate the total membrane current
- i = -imembrane + ielectrode
-
- # calculate the voltage dependent rates for the gating variables
- mtau, minf = m_lims(v)
- ntau, ninf = n_lims(v)
- htau, hinf = h_lims(v)
-
- # set the derivatives
- # note values are in SI units, as determined by the scaling factors:
- ydot[1] = i/c_m / (V/s)
- ydot[2] = (minf-m)/mtau / (1/s)
- ydot[3] = (hinf-h)/htau / (1/s)
- ydot[4] = (ninf-n)/ntau / (1/s)
-
- return Sundials.CV_SUCCESS
-end
-
-###########################################################
-# now we actually run the model
-###########################################################
-
-# from 0 to 100 ms in 1000 steps
-t = collect(linspace(0.0, 0.1, 1001));
-
-# these tolerances are as tight as they will go without breaking convergence of the iterative schemes
-y0 = initial_conditions(vrest)
-res = Sundials.cvode(f, y0, t, abstol=1e-6, reltol=5e-10);
-
-#using Plots
-#gr()
-#plot(t, res[:,1])
diff --git a/validation/ref/numeric/numeric_soma.jl b/validation/ref/numeric/numeric_soma.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6c431f845cbf5acbb3a1825e8d52951482d255a1
--- /dev/null
+++ b/validation/ref/numeric/numeric_soma.jl
@@ -0,0 +1,27 @@
+#!/usr/bin/env julia
+
+include("HHChannels.jl")
+
+using JSON
+using SIUnits.ShortUnits
+using HHChannels
+
+radius = 18.8µm/2
+area = 4*pi*radius^2
+
+stim = Stim(10ms, 100ms, 0.1nA/area)
+ts, vs = run_hh(100ms, stim=stim, sample_dt=0.025ms)
+
+trace = Dict(
+ :name => "membrane voltage",
+ :sim => "numeric",
+ :model => "soma",
+ :units => "mV",
+ :data => Dict(
+ :time => map(t->t/ms, ts),
+ symbol("soma.mid") => map(v->v/mV, vs)
+ )
+)
+
+println(JSON.json([trace]))
+