diff --git a/data/ball_and_stick.swc b/data/ball_and_stick.swc
new file mode 100644
index 0000000000000000000000000000000000000000..1737539b1b4d73d202cecab2ecb7e5599d93064e
--- /dev/null
+++ b/data/ball_and_stick.swc
@@ -0,0 +1,8 @@
+# ball and stick model with
+# - soma with radius 12.6157/2
+# - dendrite with length 200 and radius 0.5
+
+1 1 0.0 0.0 0.0 6.30785 -1
+2 2 6.30785 0.0 0.0 0.5 1
+3 2 206.30785 0.0 0.0 0.5 2
+
diff --git a/nrn/ball_and_3stick.py b/nrn/ball_and_3stick.py
new file mode 100644
index 0000000000000000000000000000000000000000..3e777261f88fde9c53407d10fe0a6fe74cec4c5e
--- /dev/null
+++ b/nrn/ball_and_3stick.py
@@ -0,0 +1,154 @@
+from timeit import default_timer as timer
+import os.path
+from matplotlib import pyplot
+import numpy as np
+import json
+import argparse
+from neuron import gui, h
+
+parser = argparse.ArgumentParser(description='generate spike train ball and stick model with hh channels at soma and pas channels in dendrite')
+parser.add_argument('--plot', action='store_true', dest='plot')
+args = parser.parse_args()
+
+soma = h.Section(name='soma')
+dend = [
+ h.Section(name='dend0'),
+ h.Section(name='dend1'),
+ h.Section(name='dend2'),
+ ]
+
+dend[0].connect(soma(1))
+dend[1].connect(dend[0](1))
+dend[2].connect(dend[0](1))
+
+# Surface area of cylinder is 2*pi*r*h (sealed ends are implicit).
+# Here we make a square cylinder in that the diameter
+# is equal to the height, so diam = h. ==> Area = 4*pi*r^2
+# We want a soma of 500 microns squared:
+# r^2 = 500/(4*pi) ==> r = 6.2078, diam = 12.6157
+soma.L = soma.diam = 12.6157 # Makes a soma of 500 microns squared.
+for d in dend :
+ d.L = 100
+ d.diam = 1
+
+for sec in h.allsec():
+ sec.Ra = 100 # Axial resistance in Ohm * cm
+ sec.cm = 1 # Membrane capacitance in micro Farads / cm^2
+
+# Insert active Hodgkin-Huxley current in the soma
+soma.insert('hh')
+soma.gnabar_hh = 0.12 # Sodium conductance in S/cm2
+soma.gkbar_hh = 0.036 # Potassium conductance in S/cm2
+soma.gl_hh = 0.0003 # Leak conductance in S/cm2
+soma.el_hh = -54.3 # Reversal potential in mV
+
+# Insert passive current in the dendrite
+for d in dend :
+ d.insert('pas')
+ d.g_pas = 0.001 # Passive conductance in S/cm2
+ d.e_pas = -65 # Leak reversal potential mV
+
+stim = [
+ h.IClamp(dend[1](1)),
+ h.IClamp(dend[2](1))
+ ]
+stim[0].delay = 5
+stim[0].dur = 80
+stim[0].amp = 4.5*0.1
+
+stim[1].delay = 40
+stim[1].dur = 10
+stim[1].amp = -2*0.1
+
+if args.plot :
+ pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
+ pyplot.grid()
+
+simdur = 100.0
+h.tstop = simdur
+h.dt = 0.001
+
+start = timer()
+results = []
+for nseg in [5, 11, 51, 101] :
+
+ print 'simulation with ', nseg, ' compartments in dendrite...'
+
+ for d in dend :
+ d.nseg=nseg
+
+ # record voltages
+ v_soma = h.Vector() # soma
+ v_dend = h.Vector() # at the dendrite branching point
+ v_clamp= h.Vector() # end of dendrite at clamp location
+
+ v_soma.record( soma(0.5)._ref_v)
+ v_dend.record( dend[0](1)._ref_v)
+ v_clamp.record(dend[1](1)._ref_v)
+
+ # record spikes
+ # this is a bit verbose, no?
+ spike_counter_soma = h.APCount(soma(0.5))
+ spike_counter_soma.thresh = -0
+ spike_counter_dend = h.APCount(dend[0](1))
+ spike_counter_dend.thresh = -20
+ spike_counter_clamp = h.APCount(dend[1](1.0))
+ spike_counter_clamp.thresh = 20
+
+ spikes_soma = h.Vector() # soma
+ spikes_dend = h.Vector() # middle of dendrite
+ spikes_clamp= h.Vector() # end of dendrite at clamp location
+
+ spike_counter_soma.record(spikes_soma)
+ spike_counter_dend.record(spikes_dend)
+ spike_counter_clamp.record(spikes_clamp)
+
+ # record time stamps
+ t_vec = h.Vector()
+ t_vec.record(h._ref_t)
+
+ # finally it's time to run the simulation
+ h.run()
+
+ results.append(
+ {
+ "nseg": nseg,
+ "dt" : h.dt,
+ "measurements": {
+ "soma" : {
+ "thresh" : spike_counter_soma.thresh,
+ "spikes" : spikes_soma.to_python()
+ },
+ "dend" : {
+ "thresh" : spike_counter_dend.thresh,
+ "spikes" : spikes_dend.to_python()
+ },
+ "clamp" : {
+ "thresh" : spike_counter_clamp.thresh,
+ "spikes" : spikes_clamp.to_python()
+ }
+ }
+ }
+ )
+
+ if args.plot :
+ pyplot.plot(t_vec, v_soma, 'k', linewidth=1, label='soma ' + str(nseg))
+ pyplot.plot(t_vec, v_dend, 'b', linewidth=1, label='dend ' + str(nseg))
+ pyplot.plot(t_vec, v_clamp, 'r', linewidth=1, label='clamp ' + str(nseg))
+
+# time the simulations
+end = timer()
+print "took ", end-start, " seconds"
+
+# save the spike info as in json format
+fp = open('ball_and_3stick.json', 'w')
+json.dump(results, fp, indent=2)
+
+if args.plot :
+ pyplot.xlabel('time (ms)')
+ pyplot.ylabel('mV')
+ pyplot.xlim([0, simdur])
+ pyplot.legend()
+
+ pyplot.show()
+
diff --git a/nrn/ball_and_stick.py b/nrn/ball_and_stick.py
index 5fe977fe594ada556ebce1da09c9a15c3f3b03e3..24562337a4baffd7a1ae3a1e7885cb07e06987ec 100644
--- a/nrn/ball_and_stick.py
+++ b/nrn/ball_and_stick.py
@@ -55,7 +55,7 @@ h.dt = 0.001
start = timer()
results = []
-for nseg in [5, 11, 21, 41, 81, 161] :
+for nseg in [5, 11, 51, 101] :
print 'simulation with ', nseg, ' compartments in dendrite...'
diff --git a/nrn/generate_validation.sh b/nrn/generate_validation.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f35d3521e9d2ce785a01bc14cd966d4910ecf7f7
--- /dev/null
+++ b/nrn/generate_validation.sh
@@ -0,0 +1,3 @@
+python2.7 ./soma.py
+python2.7 ./ball_and_stick.py
+python2.7 ./ball_and_3stick.py
diff --git a/nrn/generate_validation_data.sh b/nrn/generate_validation_data.sh
deleted file mode 100755
index 8627a0c56ecd1f2137178a193be99207cb06ef46..0000000000000000000000000000000000000000
--- a/nrn/generate_validation_data.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-python2.7 soma.py
-
diff --git a/tests/util.hpp b/tests/util.hpp
index 7a0027448e46c860471aeed68ebd8ed8e411dfe0..b5e25d72d0ab5d141f581d183771f249e48f438b 100644
--- a/tests/util.hpp
+++ b/tests/util.hpp
@@ -111,7 +111,7 @@ operator<< (std::ostream& o, spike_comparison const& spikes)
buffer, sizeof(buffer),
"min,max = %10.8f,%10.8f | mean,rms = %10.8f,%10.8f | max_rel = %10.8f",
spikes.min, spikes.max, spikes.mean, spikes.rms,
- spikes.max_relative_error()*100
+ spikes.max_relative_error()
);
return o << buffer;
}
diff --git a/tests/validate_ball_and_stick.cpp b/tests/validate_ball_and_stick.cpp
index eac29e2e84f20653acdb7555714f9ed6c6e3e555..024c7fdbde036cadf169744a1cc95269bc2778d4 100644
--- a/tests/validate_ball_and_stick.cpp
+++ b/tests/validate_ball_and_stick.cpp
@@ -8,7 +8,7 @@
#include <json/src/json.hpp>
-// compares results with those generated by nrn/soma.py
+// compares results with those generated by nrn/ball_and_stick.py
TEST(ball_and_stick, neuron_baseline)
{
using namespace nest::mc;
@@ -27,7 +27,7 @@ TEST(ball_and_stick, neuron_baseline)
auto dendrite = cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
dendrite->add_mechanism(pas_parameters());
- dendrite->mechanism("membrane").set("r_L", 100); // no effect for single compartment cell
+ dendrite->mechanism("membrane").set("r_L", 100);
// add stimulus
cell.add_stimulus({1,1}, {5., 80., 0.3});
@@ -119,17 +119,170 @@ TEST(ball_and_stick, neuron_baseline)
}
// print results
+ auto colors = {memory::util::kWhite, memory::util::kGreen, memory::util::kYellow};
for(auto& r : results){
- // select the location with the largest error for printing
- auto m =
- std::max_element(
- r.comparisons.begin(), r.comparisons.end(),
- [](testing::spike_comparison& l, testing::spike_comparison& r)
- {return l.max_relative_error()<r.max_relative_error();}
+ auto color = colors.begin();
+ for(auto const& result : r.comparisons) {
+ std::cout << std::setw(5) << r.n_comparments << " compartments : ";
+ std::cout << memory::util::colorize(util::pprintf("%\n", result), *(color++));
+ }
+ }
+
+ // sort results in ascending order of compartments
+ std::sort(
+ results.begin(), results.end(),
+ [](const result& l, const result& r)
+ {return l.n_comparments<r.n_comparments;}
+ );
+
+ // the strategy for testing is the following:
+ // 1. test that the solution converges to the finest reference solution as
+ // the number of compartments increases (i.e. spatial resolution is
+ // refined)
+ for(auto i=1u; i<results.size(); ++i) {
+ for(auto j=0; j<3; ++j) {
+ EXPECT_TRUE(
+ results[i].comparisons[j].max_relative_error()
+ < results[i-1].comparisons[j].max_relative_error()
);
- std::cout << std::setw(5) << r.n_comparments
- << " compartments : " << *m
- << "\n";
+ }
+ }
+
+ // 2. test that the best solution (i.e. with most compartments) matches the
+ // reference solution closely (less than 0.5% over the course of 100ms
+ // simulation)
+ auto tol = 0.5;
+ for(auto j=0; j<3; ++j) {
+ EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<tol);
+ }
+}
+
+// compares results with those generated by nrn/ball_and_3sticks.py
+TEST(ball_and_3stick, neuron_baseline)
+{
+ using namespace nest::mc;
+ using namespace nlohmann;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ // Soma with diameter 12.6157 um and HH channel
+ auto soma = cell.add_soma(12.6157/2.0);
+ soma->add_mechanism(hh_parameters());
+
+ // add dendrite of length 200 um and diameter 1 um with passive channel
+ std::vector<cable_segment*> dendrites;
+ dendrites.push_back(cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100));
+ dendrites.push_back(cell.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100));
+ dendrites.push_back(cell.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100));
+
+ for(auto dend : dendrites) {
+ dend->add_mechanism(pas_parameters());
+ dend->mechanism("membrane").set("r_L", 100);
+ }
+
+ // add stimulus
+ cell.add_stimulus({2,1}, {5., 80., 0.45});
+ cell.add_stimulus({3,1}, {40., 10.,-0.2});
+
+ // load data from file
+ auto cell_data = testing::load_spike_data("../nrn/ball_and_3stick.json");
+ EXPECT_TRUE(cell_data.size()>0);
+ if(cell_data.size()==0) return;
+
+ json& nrn =
+ *std::max_element(
+ cell_data.begin(), cell_data.end(),
+ [](json& lhs, json& rhs) {return lhs["nseg"]<rhs["nseg"];}
+ );
+
+ auto& measurements = nrn["measurements"];
+
+ double dt = nrn["dt"];
+ double tfinal = 100.; // ms
+ int nt = tfinal/dt;
+
+ // inline type for storing the results of a simulation along with
+ // the information required to compare two simulations for accuracy
+ struct result {
+ std::vector<std::vector<double>> spikes;
+ std::vector<std::vector<double>> baseline_spikes;
+ std::vector<testing::spike_comparison> comparisons;
+ std::vector<double> thresh;
+ int n_comparments;
+
+ result(int nseg, double dt,
+ std::vector<std::vector<double>> &v,
+ json& measurements
+ ) {
+ n_comparments = nseg;
+ baseline_spikes = {
+ measurements["soma"]["spikes"],
+ measurements["dend"]["spikes"],
+ measurements["clamp"]["spikes"]
+ };
+ thresh = {
+ measurements["soma"]["thresh"],
+ measurements["dend"]["thresh"],
+ measurements["clamp"]["thresh"]
+ };
+ for(auto i=0; i<3; ++i) {
+ // calculate the NEST MC spike times
+ spikes.push_back
+ (testing::find_spikes(v[i], thresh[i], dt));
+ // compare NEST MC and baseline NEURON spike times
+ comparisons.push_back
+ (testing::compare_spikes(spikes[i], baseline_spikes[i]));
+ }
+ }
+ };
+
+ std::vector<result> results;
+ for(auto run_index=0u; run_index<cell_data.size(); ++run_index) {
+ auto& run = cell_data[run_index];
+ int num_compartments = run["nseg"];
+ for(auto dend : dendrites) {
+ dend->set_compartments(num_compartments);
+ }
+ std::vector<std::vector<double>> v(3);
+
+ // make the lowered finite volume cell
+ fvm::fvm_cell<double, int> model(cell);
+ auto graph = cell.model();
+
+ // set initial conditions
+ using memory::all;
+ model.voltage()(all) = -65.;
+ model.initialize(); // have to do this _after_ initial conditions are set
+
+ // run the simulation
+ auto soma_comp = nest::mc::find_compartment_index({0,0.}, graph);
+ auto dend_comp = nest::mc::find_compartment_index({1,1.}, graph);
+ auto clamp_comp = nest::mc::find_compartment_index({2,1.}, graph);
+ v[0].push_back(model.voltage()[soma_comp]);
+ v[1].push_back(model.voltage()[dend_comp]);
+ v[2].push_back(model.voltage()[clamp_comp]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ // save voltage at soma
+ v[0].push_back(model.voltage()[soma_comp]);
+ v[1].push_back(model.voltage()[dend_comp]);
+ v[2].push_back(model.voltage()[clamp_comp]);
+ }
+
+ results.push_back( {num_compartments, dt, v, measurements} );
+ }
+
+ // print results
+ auto colors = {memory::util::kWhite, memory::util::kGreen, memory::util::kYellow};
+ for(auto& r : results){
+ auto color = colors.begin();
+ for(auto const& result : r.comparisons) {
+ std::cout << std::setw(5) << r.n_comparments << " compartments : ";
+ std::cout << memory::util::colorize(util::pprintf("%\n", result), *(color++));
+ }
}
// sort results in ascending order of compartments
@@ -153,10 +306,11 @@ TEST(ball_and_stick, neuron_baseline)
}
// 2. test that the best solution (i.e. with most compartments) matches the
- // reference solution closely (less than 0.1% over the course of 100ms
+ // reference solution closely (less than 0.5% over the course of 100ms
// simulation)
+ auto tol = 0.5;
for(auto j=0; j<3; ++j) {
- EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<0.1);
+ EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<tol);
}
}