diff --git a/nrn/generate_validation_data.sh b/nrn/generate_validation_data.sh
new file mode 100755
index 0000000000000000000000000000000000000000..8627a0c56ecd1f2137178a193be99207cb06ef46
--- /dev/null
+++ b/nrn/generate_validation_data.sh
@@ -0,0 +1,2 @@
+python2.7 soma.py
+
diff --git a/nrn/soma.py b/nrn/soma.py
index f8912c4f7a0ea6bc204695cbef683b01e5172e44..0e73ad624123ec70aed3e9d6dba4747524cd3c17 100644
--- a/nrn/soma.py
+++ b/nrn/soma.py
@@ -1,7 +1,18 @@
import os.path
-from neuron import h, gui
from matplotlib import pyplot
import numpy as np
+import json
+import argparse
+from neuron import gui, h
+
+parser = argparse.ArgumentParser(description='generate spike train info for a soma with hh channels')
+parser.add_argument('--plot', action='store_true', dest='plot')
+args = parser.parse_args()
+
+if args.plot :
+ print '-- plotting turned on'
+else :
+ print '-- plotting turned off'
soma = h.Section(name='soma')
@@ -17,51 +28,43 @@ stim.delay = 10
stim.dur = 100
stim.amp = 0.1
+spike_counter = h.APCount(soma(0.5))
+spike_counter.thresh = 0
+
v_vec = h.Vector() # Membrane potential vector
t_vec = h.Vector() # Time stamp vector
+s_vec = h.Vector() # Time stamp vector
v_vec.record(soma(0.5)._ref_v)
t_vec.record(h._ref_t)
+spike_counter.record(s_vec)
+
simdur = 120
# initialize plot
-pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
-
-pyplot.subplot(2,1,1)
-pyplot.grid()
+if args.plot :
+ pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
h.tstop = simdur
# run neuron with multiple dt
-
-for dt in [0.02, 0.01, 0.005, 0.002, 0.001]:
+results = []
+#for dt in [0.02, 0.01, 0.005, 0.0005, 0.0001]:
+for dt in [0.02, 0.0001]:
h.dt = dt
h.run()
- pyplot.plot(t_vec, v_vec, label='neuron ' + str(dt))
- #pyplot.plot(t_vec, v_vec, 'k', label='neuron ' + str(dt))
-
-pyplot.xlim([102.5,105])
-pyplot.ylim([0,40])
-
-pyplot.legend()
-pyplot.subplot(2,1,2)
-pyplot.grid()
-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)
- pyplot.plot(t, v, label=fname)
- step = step + 1
-
-pyplot.xlabel('time (ms)')
-pyplot.ylabel('mV')
-
-pyplot.xlim([102.5,105])
-pyplot.ylim([0,40])
-
-pyplot.legend()
-pyplot.show()
+ results.append({"dt": dt, "spikes": s_vec.to_python()})
+ if args.plot :
+ pyplot.plot(t_vec, v_vec, label='neuron ' + str(dt))
+
+# save the spike info as in json format
+fp = open('soma.json', 'w')
+json.dump(results, fp, indent=1)
+
+if args.plot :
+ pyplot.xlabel('time (ms)')
+ pyplot.ylabel('mV')
+ pyplot.xlim([0, 120])
+ pyplot.grid()
+ pyplot.legend()
+ pyplot.show()
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 19b27c45268f8be09a4d10e9a79adf2a5d714ac6..833c8cd15953f93d49f9e9ff9c67c2ba43674396 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -257,8 +257,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
cv_capacitance_[i] /= cv_areas_[i];
}
- std::cout << "capacitance " << cv_capacitance_ << "\n";
-
/////////////////////////////////////////////
// create mechanisms
/////////////////////////////////////////////
diff --git a/tests/util.hpp b/tests/util.hpp
index cf10847fbf0de27657be11eb54b545538c3799b5..873bc2055fb6d41f2fda60887223a37a8ec60556 100644
--- a/tests/util.hpp
+++ b/tests/util.hpp
@@ -3,7 +3,7 @@
#include <string>
#include <vector>
-static
+[[gnu::unused]] static
void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
{
auto m = values.size();
@@ -28,3 +28,30 @@ void write_vis_file(const std::string& fname, std::vector<std::vector<double>> v
}
}
+template <typename T>
+std::vector<T> find_spikes(std::vector<T> const& v, T threshold, T dt)
+{
+ if(v.size()<2) {
+ return {};
+ }
+ bool is_spiking = v[0] >= threshold;
+ auto it = v.begin() + 1;
+
+ std::vector<T> times;
+ for(auto i=1u; i<v.size(); ++i) {
+ if(is_spiking && v[i]<threshold) {
+ is_spiking = false;
+ }
+ else {
+ if(*it>=threshold) {
+ is_spiking = true;
+ auto pos = (threshold-v[i-1]) / (v[i]-v[i-1]);
+ times.push_back((i-1+pos)*dt);
+ }
+ }
+ }
+
+
+ return times;
+}
+
diff --git a/tests/validate_soma.cpp b/tests/validate_soma.cpp
index b88c09f40b7426d54b06047cb290c993cb9083d5..70fe80c384e55a9b88c55aff7e020f7c6100bf65 100644
--- a/tests/validate_soma.cpp
+++ b/tests/validate_soma.cpp
@@ -3,13 +3,16 @@
#include "gtest.h"
#include "util.hpp"
-#include "../src/cell.hpp"
-#include "../src/fvm.hpp"
+#include <cell.hpp>
+#include <fvm.hpp>
+
+#include <json/src/json.hpp>
// based on hh/Neuron/steps_A.py
TEST(soma, resolutions)
{
using namespace nest::mc;
+ using namespace nlohmann;
nest::mc::cell cell;
@@ -22,15 +25,34 @@ TEST(soma, resolutions)
soma->add_mechanism(hh_parameters());
// add stimulus to the soma
- //cell.add_stimulus({0,0}, {5., 3., 0.1});
cell.add_stimulus({0,0.5}, {10., 100., 0.1});
// make the lowered finite volume cell
fvm::fvm_cell<double, int> model(cell);
- auto i =0;
- for(auto dt : {0.02, 0.01, 0.005, 0.002, 0.001}) {
- std::vector<std::vector<double>> results(2);
+ // 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;
+ }
+
+ try {
+ fid >> cell_data;
+ }
+ catch (...) {
+ std::cerr << "error : incorrectly formatted json file " << input_name << "\n";
+ return;
+ }
+ }
+
+ for(auto& run : cell_data) {
+ std::vector<double> v;
+ double dt = run["dt"];
// set initial conditions
using memory::all;
@@ -40,16 +62,16 @@ TEST(soma, resolutions)
// run the simulation
auto tfinal = 120.; // ms
int nt = tfinal/dt;
- results[0].push_back(0.);
- results[1].push_back(model.voltage()[0]);
+ v.push_back(model.voltage()[0]);
for(auto i=0; i<nt; ++i) {
model.advance(dt);
// save voltage at soma
- results[0].push_back((i+1)*dt);
- results[1].push_back(model.voltage()[0]);
+ v.push_back(model.voltage()[0]);
}
- write_vis_file("v_" + std::to_string(i) + ".dat", results);
- ++i;
+
+ std::cout << v << "\n";
+ auto spike_times = find_spikes(v, 0., dt);
+ std::cout << "dt " << dt << " : spikes " << spike_times << "\n";
}
}