diff --git a/nrn/ball_and_stick.py b/nrn/ball_and_stick.py
index ce857fc1a103ab20e00e2b2a60ade74ecffa7de9..5fe977fe594ada556ebce1da09c9a15c3f3b03e3 100644
--- a/nrn/ball_and_stick.py
+++ b/nrn/ball_and_stick.py
@@ -1,13 +1,20 @@
-from neuron import h, gui
+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='dend')
dend.connect(soma(1))
-h.psection(sec=soma)
-
# 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
@@ -16,7 +23,6 @@ h.psection(sec=soma)
soma.L = soma.diam = 12.6157 # Makes a soma of 500 microns squared.
dend.L = 200 # microns
dend.diam = 1 # microns
-print "Surface area of soma =", h.area(0.5, sec=soma)
for sec in h.allsec():
sec.Ra = 100 # Axial resistance in Ohm * cm
@@ -34,51 +40,99 @@ dend.insert('pas')
dend.g_pas = 0.001 # Passive conductance in S/cm2
dend.e_pas = -65 # Leak reversal potential mV
-#for sec in h.allsec():
-# h.psection(sec=sec)
-
stim = h.IClamp(dend(1))
stim.delay = 5
-stim.dur = 1
-stim.amp = 5*0.1
+stim.dur = 80
+stim.amp = 3*0.1
-somastyles = ['b-', 'r-']
-dendstyles = ['b-.', 'r-.']
-dt = [0.02, 0.001]
-nsegs = [10, 100]
-from matplotlib import pyplot
-pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
-pyplot.grid()
-pos = 0
-for pos in [0] :
- dend.nseg=nsegs[pos]
-
- v_soma = h.Vector() # Membrane potential vector
- v_dend = h.Vector() # Membrane potential vector
- t_vec = h.Vector() # Time stamp vector
- v_soma.record(soma(0.5)._ref_v)
- v_dend.record(dend(1.0)._ref_v)
- t_vec.record(h._ref_t)
- simdur = 25.0
+if args.plot :
+ pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
+ pyplot.grid()
- h.tstop = simdur
- h.dt = dt[pos]
- h.run()
+simdur = 100.0
+h.tstop = simdur
+h.dt = 0.001
+
+start = timer()
+results = []
+for nseg in [5, 11, 21, 41, 81, 161] :
- pyplot.plot(t_vec, v_soma, somastyles[pos], linewidth=2)
- pyplot.plot(t_vec, v_dend, dendstyles[pos], linewidth=2)
+ print 'simulation with ', nseg, ' compartments in dendrite...'
-pyplot.xlabel('time (ms)')
-pyplot.ylabel('mV')
+ dend.nseg=nseg
-data = np.loadtxt('../tests/v.dat')
-t = data[:,0]
-nfields = data.shape[1]
+ # record voltages
+ v_soma = h.Vector() # soma
+ v_dend = h.Vector() # middle of dendrite
+ v_clamp= h.Vector() # end of dendrite at clamp location
-v = data[:,1]
-pyplot.plot(t, v, 'g', linewidth=2)
-v = data[:,2]
-pyplot.plot(t, v, 'g-.', linewidth=2)
+ v_soma.record( soma(0.5)._ref_v)
+ v_dend.record( dend(0.5)._ref_v)
+ v_clamp.record(dend(1.0)._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.5))
+ spike_counter_dend.thresh = -10
+ spike_counter_clamp = h.APCount(dend(1.0))
+ spike_counter_clamp.thresh = 10
+
+ 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()
-pyplot.show()
+ 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_stick.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/nest_soma.py b/nrn/nest_soma.py
deleted file mode 100644
index 4acec3ed5cb1ac670abc58acc6099dc1dbf696ed..0000000000000000000000000000000000000000
--- a/nrn/nest_soma.py
+++ /dev/null
@@ -1,23 +0,0 @@
-import os.path
-from matplotlib import pyplot
-import numpy as np
-
-step = 0
-while os.path.isfile('../tests/v_' + str(step) + '.dat') :
- fname = '../tests/v_' + str(step) + '.dat'
- print 'loading ' + fname
- data = np.loadtxt(fname)
- t = data[:,0]
- v = data[:,1]
- pyplot.plot(t, v, 'b', label=fname)
- step = step + 1
-
-pyplot.xlabel('time (ms)')
-pyplot.ylabel('mV')
-
-pyplot.xlim([87.5,89.5])
-#pyplot.ylim([0,40])
-
-pyplot.legend()
-pyplot.show()
-
diff --git a/nrn/soma.py b/nrn/soma.py
index 65661b6f4f3739e2e6f5c743e2f8058bd469ed8d..3af4f50f3e6e1ae7761591ab68fd32a1187f9c87 100644
--- a/nrn/soma.py
+++ b/nrn/soma.py
@@ -50,14 +50,14 @@ h.tstop = simdur
# run neuron with multiple dt
start = timer()
results = []
-for dt in [0.02, 0.01, 0.001, 0.0001]:
+for dt in [0.05, 0.02, 0.01, 0.001, 0.0001]:
h.dt = dt
h.run()
results.append({"dt": dt, "spikes": s_vec.to_python()})
if args.plot :
pyplot.plot(t_vec, v_vec, label='neuron ' + str(dt))
-end = timer()
+end = timer()
print "took ", end-start, " seconds"
# save the spike info as in json format
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4cc842363ad1474e408f7b44e9bb717197bb33fb..11b88107784c30865635a48e9386b726505e9618 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -32,6 +32,7 @@ set(TEST_SOURCES
set(VALIDATION_SOURCES
# unit tests
validate_soma.cpp
+ validate_ball_and_stick.cpp
# unit test driver
validate.cpp
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 41f17ab5187c640093352eea88e20dae018c78a3..aa4d9bba21448e6607d1490af61512a6ba25c812 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -48,7 +48,7 @@ TEST(run, single_compartment)
results[1].push_back(model.voltage()[0]);
}
- write_vis_file("v.dat", results);
+ testing::write_vis_file("v.dat", results);
}
TEST(run, ball_and_stick)
@@ -102,7 +102,7 @@ TEST(run, ball_and_stick)
results[2].push_back(model.voltage()[pos]);
}
- write_vis_file("v.dat", results);
+ testing::write_vis_file("v.dat", results);
}
TEST(run, cable)
diff --git a/tests/util.hpp b/tests/util.hpp
index 19b0bfbb92a8a12c041ef05ff897ff7306fb7154..7a0027448e46c860471aeed68ebd8ed8e411dfe0 100644
--- a/tests/util.hpp
+++ b/tests/util.hpp
@@ -2,10 +2,17 @@
#include <iostream>
#include <string>
#include <vector>
+#include <iomanip>
#include <cmath>
#include <util.hpp>
+#include <json/src/json.hpp>
+
+// helpful code for running tests
+// a bit messy: refactor when it gets heavier and obvious patterns emerge...
+
+namespace testing{
[[gnu::unused]] static
void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
@@ -32,6 +39,28 @@ void write_vis_file(const std::string& fname, std::vector<std::vector<double>> v
}
}
+[[gnu::unused]] static
+nlohmann::json
+load_spike_data(const std::string& input_name)
+{
+ nlohmann::json cell_data;
+ auto fid = std::ifstream(input_name);
+ if(!fid.is_open()) {
+ std::cerr << "error : unable to open file " << input_name
+ << " : run the validation generation script first\n";
+ return {};
+ }
+
+ try {
+ fid >> cell_data;
+ }
+ catch (...) {
+ std::cerr << "error : incorrectly formatted json file " << input_name << "\n";
+ return {};
+ }
+ return cell_data;
+}
+
template <typename T>
std::vector<T> find_spikes(std::vector<T> const& v, T threshold, T dt)
{
@@ -72,16 +101,25 @@ struct spike_comparison {
}
};
+[[gnu::unused]] static
std::ostream&
operator<< (std::ostream& o, spike_comparison const& spikes)
{
- return o << "[" << spikes.min << ", " << spikes.max << "], "
- << "mean " << spikes.mean << ", rms " << spikes.rms
- << ", diffs " << spikes.diff;
+ // use snprintf because C++ is just awful for formatting output
+ char buffer[512];
+ snprintf(
+ 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
+ );
+ return o << buffer;
}
template <typename T>
-spike_comparison compare_spikes(std::vector<T> const& spikes, std::vector<T> const& baseline)
+spike_comparison compare_spikes(
+ std::vector<T> const& spikes,
+ std::vector<T> const& baseline)
{
spike_comparison c;
@@ -112,3 +150,4 @@ spike_comparison compare_spikes(std::vector<T> const& spikes, std::vector<T> con
return c;
}
+} // namespace testing
diff --git a/tests/validate_ball_and_stick.cpp b/tests/validate_ball_and_stick.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eac29e2e84f20653acdb7555714f9ed6c6e3e555
--- /dev/null
+++ b/tests/validate_ball_and_stick.cpp
@@ -0,0 +1,162 @@
+#include <fstream>
+
+#include "gtest.h"
+#include "util.hpp"
+
+#include <cell.hpp>
+#include <fvm.hpp>
+
+#include <json/src/json.hpp>
+
+// compares results with those generated by nrn/soma.py
+TEST(ball_and_stick, 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
+ 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
+
+ // add stimulus
+ cell.add_stimulus({1,1}, {5., 80., 0.3});
+
+ // load data from file
+ auto cell_data = testing::load_spike_data("../nrn/ball_and_stick.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"];
+ dendrite->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,0.5}, graph);
+ auto clamp_comp = nest::mc::find_compartment_index({1,1.0}, 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( {run["nseg"], dt, v, measurements} );
+ }
+
+ // print results
+ 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();}
+ );
+ std::cout << std::setw(5) << r.n_comparments
+ << " compartments : " << *m
+ << "\n";
+ }
+
+ // 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()
+ );
+ }
+ }
+
+ // 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
+ // simulation)
+ for(auto j=0; j<3; ++j) {
+ EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<0.1);
+ }
+}
+
diff --git a/tests/validate_soma.cpp b/tests/validate_soma.cpp
index a90e0fa8cd59057afd6ca9ea91d3990c93e9ffcc..b65b4ecb3a0de12d475c6380b24b43b824b65599 100644
--- a/tests/validate_soma.cpp
+++ b/tests/validate_soma.cpp
@@ -8,8 +8,9 @@
#include <json/src/json.hpp>
-// based on hh/Neuron/steps_A.py
-TEST(soma, resolutions)
+// compares results with those generated by nrn/soma.py
+// single compartment model with HH channels
+TEST(soma, neuron_baseline)
{
using namespace nest::mc;
using namespace nlohmann;
@@ -31,24 +32,19 @@ TEST(soma, resolutions)
fvm::fvm_cell<double, int> model(cell);
// load data from file
- std::string input_name = "../nrn/soma.json";
- json cell_data;
- {
- auto fid = std::ifstream(input_name);
- if(!fid.is_open()) {
- std::cerr << "error : unable to open file " << input_name
- << " : run the validation generation script first\n";
- return;
- }
+ auto cell_data = testing::load_spike_data("../nrn/soma.json");
+ EXPECT_TRUE(cell_data.size()>0);
+ if(cell_data.size()==0) return;
- try {
- fid >> cell_data;
- }
- catch (...) {
- std::cerr << "error : incorrectly formatted json file " << input_name << "\n";
- return;
- }
- }
+ json& nrn =
+ *std::min_element(
+ cell_data.begin(), cell_data.end(),
+ [](json& lhs, json& rhs) {return lhs["dt"]<rhs["dt"];}
+ );
+ std::vector<double> nrn_spike_times = nrn["spikes"];
+
+ std::cout << "baseline with time step size " << nrn["dt"] << " ms\n";
+ std::cout << "baseline spikes : " << nrn["spikes"] << "\n";
for(auto& run : cell_data) {
std::vector<double> v;
@@ -69,15 +65,99 @@ TEST(soma, resolutions)
v.push_back(model.voltage()[0]);
}
- // get the spike times from the NEST MC and NEURON simulations respectively
- auto nst_spike_times = find_spikes(v, 0., dt);
- auto nrn_spike_times = run["spikes"].get<std::vector<double>>();
- auto comparison = compare_spikes(nst_spike_times, nrn_spike_times);
+ // get the spike times from the NEST MC simulation
+ auto nst_spike_times = testing::find_spikes(v, 0., dt);
+ // compare NEST MC and baseline NEURON results
+ auto comparison = testing::compare_spikes(nst_spike_times, nrn_spike_times);
// Assert that relative error is less than 1%.
- // For a 100 ms simulation this asserts that the difference between NEST and NEURON
- // at a given time step size is less than 1 ms.
+ // For a 100 ms simulation this asserts that the difference between
+ // NEST and the most accurate NEURON simulation is less than 1 ms.
+ std::cout << "MAX ERROR @ " << dt << " is " << comparison.max_relative_error()*100. << "\n";
EXPECT_TRUE(comparison.max_relative_error()*100. < 1.);
}
}
+// check if solution converges
+TEST(soma, convergence)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ // Soma with diameter 18.8um and HH channel
+ auto soma = cell.add_soma(18.8/2.0);
+ soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
+ soma->add_mechanism(hh_parameters());
+
+ // add stimulus to the soma
+ cell.add_stimulus({0,0.5}, {10., 100., 0.1});
+
+ // make the lowered finite volume cell
+ fvm::fvm_cell<double, int> model(cell);
+
+ // generate baseline solution with small dt=0.0001
+ std::vector<double> baseline_spike_times;
+ {
+ auto dt = 1e-4;
+ std::vector<double> v;
+
+ // 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 tfinal = 120.; // ms
+ int nt = tfinal/dt;
+ v.push_back(model.voltage()[0]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ // save voltage at soma
+ v.push_back(model.voltage()[0]);
+ }
+
+ baseline_spike_times = testing::find_spikes(v, 0., dt);
+ }
+
+ testing::spike_comparison previous;
+ for(auto dt : {0.05, 0.02, 0.01, 0.005, 0.001}) {
+ std::vector<double> v;
+
+ // 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 tfinal = 120.; // ms
+ int nt = tfinal/dt;
+ v.push_back(model.voltage()[0]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ // save voltage at soma
+ v.push_back(model.voltage()[0]);
+ }
+
+ // get the spike times from the NEST MC simulation
+ auto nst_spike_times = testing::find_spikes(v, 0., dt);
+
+ // compare NEST MC and baseline NEURON results
+ auto comparison = testing::compare_spikes(nst_spike_times, baseline_spike_times);
+ std::cout << "dt " << std::setw(8) << dt << " : " << comparison << "\n";
+ if(!previous.is_valid()) {
+ previous = std::move(comparison);
+ }
+ else {
+ // Assert that relative error is less than 1%.
+ // For a 100 ms simulation this asserts that the difference between
+ // NEST and the most accurate NEURON simulation is less than 1 ms.
+ EXPECT_TRUE(comparison.max_relative_error() < previous.max_relative_error());
+ EXPECT_TRUE(comparison.rms < previous.rms);
+ EXPECT_TRUE(comparison.max < previous.max);
+ }
+ }
+}