single compartment validation test

nrn/soma.py

+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')

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} );
     }
 }

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;

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;
 }
 

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.);
     }
 }