add validation test for ball and stick

nrn/ball_and_stick.py

-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()
 

nrn/nest_soma.py

-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()
-

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

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

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)

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

tests/validate_ball_and_stick.cpp

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

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