From 781481ce06fc4f6795441f2f8270a689a37021a4 Mon Sep 17 00:00:00 2001
From: Ben Cumming <bcumming@cscs.ch>
Date: Mon, 23 May 2016 11:17:40 +0200
Subject: [PATCH] single compartment validation test
---
nrn/soma.py | 8 +++--
src/fvm.hpp | 1 -
src/util.hpp | 4 +--
tests/util.hpp | 79 +++++++++++++++++++++++++++++++++++------
tests/validate_soma.cpp | 12 +++++--
5 files changed, 85 insertions(+), 19 deletions(-)
diff --git a/nrn/soma.py b/nrn/soma.py
index 0e73ad62..65661b6f 100644
--- a/nrn/soma.py
+++ b/nrn/soma.py
@@ -1,3 +1,4 @@
+from timeit import default_timer as timer
import os.path
from matplotlib import pyplot
import numpy as np
@@ -47,14 +48,17 @@ if args.plot :
h.tstop = simdur
# run neuron with multiple dt
+start = timer()
results = []
-#for dt in [0.02, 0.01, 0.005, 0.0005, 0.0001]:
-for dt in [0.02, 0.0001]:
+for dt in [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()
+
+print "took ", end-start, " seconds"
# save the spike info as in json format
fp = open('soma.json', 'w')
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 833c8cd1..69d23d0d 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -375,7 +375,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// add the stimulii
for(const auto& stim : cell.stimulii()) {
auto idx = find_compartment_index(stim.first, graph);
- std::cout << "adding stimulus at compartment " << idx << "\n";
stimulii_.push_back( {idx, stim.second} );
}
}
diff --git a/src/util.hpp b/src/util.hpp
index 3dae43f8..706ab7da 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -28,7 +28,7 @@ operator << (std::ostream &o, std::vector<T>const& v)
{
o << "[";
for(auto const& i: v) {
- o << i << ", ";
+ o << i << ",";
}
o << "]";
return o;
@@ -39,7 +39,7 @@ std::ostream& print(std::ostream &o, std::vector<T>const& v)
{
o << "[";
for(auto const& i: v) {
- o << i << ", ";
+ o << i << ",";
}
o << "]";
return o;
diff --git a/tests/util.hpp b/tests/util.hpp
index 873bc205..19b0bfbb 100644
--- a/tests/util.hpp
+++ b/tests/util.hpp
@@ -3,6 +3,10 @@
#include <string>
#include <vector>
+#include <cmath>
+
+#include <util.hpp>
+
[[gnu::unused]] static
void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
{
@@ -34,24 +38,77 @@ 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;
+ if(v[i]>=threshold && v[i-1]<threshold) {
+ auto pos = (threshold-v[i-1]) / (v[i]-v[i-1]);
+ times.push_back((i-1+pos)*dt);
}
- 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;
+}
+
+struct spike_comparison {
+ double min = std::numeric_limits<double>::quiet_NaN();
+ double max = std::numeric_limits<double>::quiet_NaN();
+ double mean = std::numeric_limits<double>::quiet_NaN();
+ double rms = std::numeric_limits<double>::quiet_NaN();
+ std::vector<double> diff;
+
+ // check whether initialized (i.e. has valid results)
+ bool is_valid() const {
+ return min == min;
+ }
+
+ // return maximum relative error
+ double max_relative_error() const {
+ if(!is_valid()) {
+ return std::numeric_limits<double>::quiet_NaN();
}
+
+ return *std::max_element(diff.begin(), diff.end());
+ }
+};
+
+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;
+}
+
+template <typename T>
+spike_comparison compare_spikes(std::vector<T> const& spikes, std::vector<T> const& baseline)
+{
+ spike_comparison c;
+
+ // return default initialized (all NaN) if number of spikes differs
+ if(spikes.size() != baseline.size()) {
+ return c;
}
+ c.min = std::numeric_limits<double>::max();
+ c.max = 0.;
+ c.mean = 0.;
+ c.rms = 0.;
- return times;
+ auto n = spikes.size();
+ for(auto i=0u; i<n; ++i) {
+ auto error = std::fabs(spikes[i] - baseline[i]);
+ c.min = std::min(c.min, error);
+ c.max = std::max(c.max, error);
+ c.mean += error;
+ c.rms += error*error;
+ // relative difference
+ c.diff.push_back(error/baseline[i]);
+ }
+
+ c.mean /= n;
+ c.rms = std::sqrt(c.rms/n);
+
+ return c;
}
diff --git a/tests/validate_soma.cpp b/tests/validate_soma.cpp
index 70fe80c3..a90e0fa8 100644
--- a/tests/validate_soma.cpp
+++ b/tests/validate_soma.cpp
@@ -69,9 +69,15 @@ TEST(soma, resolutions)
v.push_back(model.voltage()[0]);
}
- std::cout << v << "\n";
- auto spike_times = find_spikes(v, 0., dt);
- std::cout << "dt " << dt << " : spikes " << spike_times << "\n";
+ // 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);
+
+ // 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.
+ EXPECT_TRUE(comparison.max_relative_error()*100. < 1.);
}
}
--
GitLab