diff --git a/CMakeLists.txt b/CMakeLists.txt
index 57991c886fc85752b0d7a791375ba8004ce69012..ccc56d8a708f7c3270277f2c54796656655aacfc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,6 +14,7 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS "YES")
set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
+include_directories(${CMAKE_SOURCE_DIR}/external)
include_directories(${CMAKE_SOURCE_DIR}/include)
include_directories(${CMAKE_SOURCE_DIR}/src)
include_directories(${CMAKE_SOURCE_DIR})
diff --git a/data/ball_and_stick.swc b/data/ball_and_stick.swc
new file mode 100644
index 0000000000000000000000000000000000000000..1737539b1b4d73d202cecab2ecb7e5599d93064e
--- /dev/null
+++ b/data/ball_and_stick.swc
@@ -0,0 +1,8 @@
+# ball and stick model with
+# - soma with radius 12.6157/2
+# - dendrite with length 200 and radius 0.5
+
+1 1 0.0 0.0 0.0 6.30785 -1
+2 2 6.30785 0.0 0.0 0.5 1
+3 2 206.30785 0.0 0.0 0.5 2
+
diff --git a/mechanisms/generate.sh b/mechanisms/generate.sh
index d9033f0dbc1cb34e6ef10b1b42d6bdd9e37b0cd7..94755fbb91b912e5eada0148af1c72ea91945ecd 100755
--- a/mechanisms/generate.sh
+++ b/mechanisms/generate.sh
@@ -1,4 +1,4 @@
for mech in pas hh
do
- ../modparser/bin/modcc -t cpu -o ../include/mechanisms/$mech.hpp ./mod/$mech.mod
+ ../external/modparser/bin/modcc -t cpu -o ../include/mechanisms/$mech.hpp ./mod/$mech.mod
done
diff --git a/nrn/ball_and_3stick.py b/nrn/ball_and_3stick.py
new file mode 100644
index 0000000000000000000000000000000000000000..15ece3d3da151daf7bdb90f1af1d6cd368352489
--- /dev/null
+++ b/nrn/ball_and_3stick.py
@@ -0,0 +1,154 @@
+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='dend0'),
+ h.Section(name='dend1'),
+ h.Section(name='dend2'),
+ ]
+
+dend[0].connect(soma(1))
+dend[1].connect(dend[0](1))
+dend[2].connect(dend[0](1))
+
+# 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
+# We want a soma of 500 microns squared:
+# r^2 = 500/(4*pi) ==> r = 6.2078, diam = 12.6157
+soma.L = soma.diam = 12.6157 # Makes a soma of 500 microns squared.
+for d in dend :
+ d.L = 100
+ d.diam = 1
+
+for sec in h.allsec():
+ sec.Ra = 100 # Axial resistance in Ohm * cm
+ sec.cm = 1 # Membrane capacitance in micro Farads / cm^2
+
+# Insert active Hodgkin-Huxley current in the soma
+soma.insert('hh')
+soma.gnabar_hh = 0.12 # Sodium conductance in S/cm2
+soma.gkbar_hh = 0.036 # Potassium conductance in S/cm2
+soma.gl_hh = 0.0003 # Leak conductance in S/cm2
+soma.el_hh = -54.3 # Reversal potential in mV
+
+# Insert passive current in the dendrite
+for d in dend :
+ d.insert('pas')
+ d.g_pas = 0.001 # Passive conductance in S/cm2
+ d.e_pas = -65 # Leak reversal potential mV
+
+stim = [
+ h.IClamp(dend[1](1)),
+ h.IClamp(dend[2](1))
+ ]
+stim[0].delay = 5
+stim[0].dur = 80
+stim[0].amp = 4.5*0.1
+
+stim[1].delay = 40
+stim[1].dur = 10
+stim[1].amp = -2*0.1
+
+if args.plot :
+ pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
+ pyplot.grid()
+
+simdur = 100.0
+h.tstop = simdur
+h.dt = 0.001
+
+start = timer()
+results = []
+for nseg in [5, 11, 51, 101] :
+
+ print 'simulation with ', nseg, ' compartments in dendrite...'
+
+ for d in dend :
+ d.nseg=nseg
+
+ # record voltages
+ v_soma = h.Vector() # soma
+ v_dend = h.Vector() # at the dendrite branching point
+ v_clamp= h.Vector() # end of dendrite at clamp location
+
+ v_soma.record( soma(0.5)._ref_v)
+ v_dend.record( dend[0](1)._ref_v)
+ v_clamp.record(dend[1](1)._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](1))
+ spike_counter_dend.thresh = -20
+ spike_counter_clamp = h.APCount(dend[1](1.0))
+ spike_counter_clamp.thresh = 20
+
+ 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()
+
+ 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_3stick.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/ball_and_stick.py b/nrn/ball_and_stick.py
index 5fe977fe594ada556ebce1da09c9a15c3f3b03e3..24562337a4baffd7a1ae3a1e7885cb07e06987ec 100644
--- a/nrn/ball_and_stick.py
+++ b/nrn/ball_and_stick.py
@@ -55,7 +55,7 @@ h.dt = 0.001
start = timer()
results = []
-for nseg in [5, 11, 21, 41, 81, 161] :
+for nseg in [5, 11, 51, 101] :
print 'simulation with ', nseg, ' compartments in dendrite...'
diff --git a/nrn/generate_validation.sh b/nrn/generate_validation.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f35d3521e9d2ce785a01bc14cd966d4910ecf7f7
--- /dev/null
+++ b/nrn/generate_validation.sh
@@ -0,0 +1,3 @@
+python2.7 ./soma.py
+python2.7 ./ball_and_stick.py
+python2.7 ./ball_and_3stick.py
diff --git a/nrn/generate_validation_data.sh b/nrn/generate_validation_data.sh
deleted file mode 100755
index 8627a0c56ecd1f2137178a193be99207cb06ef46..0000000000000000000000000000000000000000
--- a/nrn/generate_validation_data.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-python2.7 soma.py
-
diff --git a/src/cell.cpp b/src/cell.cpp
index a6986f0673ca05827ee89493aca2230ec16cec30..8fdd1bb6e0f55cdfa2a3c7ea6632e0fdee74bba6 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -195,5 +195,41 @@ std::vector<int> const& cell::segment_parents() const
return parents_;
}
+// Rough and ready comparison of two cells.
+// We don't use an operator== because equality of two cells is open to
+// interpretation. For example, it is possible to have two viable representations
+// of a cell: with and without location information for the cables.
+//
+// Checks that two cells have the same
+// - number and type of segments
+// - volume and area properties of each segment
+// - number of compartments in each segment
+bool cell_basic_equality(cell const& lhs, cell const& rhs)
+{
+ if(lhs.num_segments() != rhs.num_segments()) {
+ return false;
+ }
+ if(lhs.segment_parents() != rhs.segment_parents()) {
+ return false;
+ }
+ for(auto i=0; i<lhs.num_segments(); ++i) {
+ // a quick and dirty test
+ auto& l = *lhs.segment(i);
+ auto& r = *rhs.segment(i);
+
+ if(l.kind() != r.kind()) return false;
+ if(l.area() != r.area()) return false;
+ if(l.volume() != r.volume()) return false;
+ if(l.as_cable()) {
+ if( l.as_cable()->num_compartments()
+ != r.as_cable()->num_compartments())
+ {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
} // namespace mc
} // namespace nest
diff --git a/src/cell.hpp b/src/cell.hpp
index ba5f9b6bff227231cbb4f46e27dd2af4ff2943ce..fece134947041d705564008e54be61ce93d50f60 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -123,6 +123,12 @@ class cell {
std::vector<std::pair<segment_location, i_clamp>> stimulii_;
};
+// Checks that two cells have the same
+// - number and type of segments
+// - volume and area properties of each segment
+// - number of compartments in each segment
+bool cell_basic_equality(cell const& lhs, cell const& rhs);
+
// create a cable by forwarding cable construction parameters provided by the user
template <typename... Args>
cable_segment* cell::add_cable(cell::index_type parent, Args ...args)
diff --git a/src/segment.hpp b/src/segment.hpp
index fd1d8ded6101cc8e9834b9062df3da97de5cbc8c..decb0cdd88665c421cf159b8d7f49e4f91e03f0f 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -76,6 +76,16 @@ class segment {
return nullptr;
}
+ virtual const cable_segment* as_cable() const
+ {
+ return nullptr;
+ }
+
+ virtual const soma_segment* as_soma() const
+ {
+ return nullptr;
+ }
+
virtual bool is_placeholder() const
{
return false;
@@ -230,6 +240,11 @@ class soma_segment : public segment
return this;
}
+ const soma_segment* as_soma() const override
+ {
+ return this;
+ }
+
/// soma has one and one only compartments
int num_compartments() const override
{
@@ -359,6 +374,11 @@ class cable_segment : public segment
return this;
}
+ const cable_segment* as_cable() const override
+ {
+ return this;
+ }
+
int num_compartments() const override
{
return num_compartments_;
diff --git a/src/util.hpp b/src/util.hpp
index 706ab7da5a17d493823b8c74d1f18fa6b4a9efaf..afa029ff9cfcba65258473875003c00f58c3d16b 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -45,6 +45,21 @@ std::ostream& print(std::ostream &o, std::vector<T>const& v)
return o;
}
+template <typename T>
+bool operator ==(const std::vector<T>& lhs, const std::vector<T>& rhs)
+{
+ if(lhs.size() != rhs.size()) {
+ return false;
+ }
+ return std::equal(lhs.begin(), lhs.end(), rhs.begin());
+}
+
+template <typename T>
+bool operator !=(const std::vector<T>& lhs, const std::vector<T>& rhs)
+{
+ return !(lhs==rhs);
+}
+
namespace nest {
namespace mc {
namespace util {
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 679260d656c9c7964f0eea599951c371886123f4..6968fd432a89842023807c8ed0ffb2ca8140b231 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -487,3 +487,31 @@ TEST(swc_io, cell_construction)
cell.cable(3)->radius(0));
}
}
+
+// check that simple ball and stick model with one dendrite attached to a soma
+// which is used in the validation tests can be loaded from file and matches
+// the one generated with the C++ interface
+TEST(swc_parser, from_file_ball_and_stick)
+{
+ auto fname = "../data/ball_and_stick.swc";
+ std::ifstream fid(fname);
+ if(!fid.is_open()) {
+ std::cerr << "unable to open file " << fname << "... skipping test\n";
+ return;
+ }
+
+ // read the file into a cell object
+ auto cell = nest::mc::io::swc_read_cell(fid);
+
+ // verify that the correct number of nodes was read
+ EXPECT_EQ(cell.num_segments(), 2);
+ EXPECT_EQ(cell.num_compartments(), 2u);
+
+ // make an equivalent cell via C++ interface
+ nest::mc::cell local_cell;
+ local_cell.add_soma(6.30785);
+ local_cell.add_cable(0, nest::mc::segmentKind::dendrite, 0.5, 0.5, 200);
+
+ EXPECT_TRUE(nest::mc::cell_basic_equality(local_cell, cell));
+}
+
diff --git a/tests/util.hpp b/tests/util.hpp
index 7a0027448e46c860471aeed68ebd8ed8e411dfe0..aeb2f8d3e585a1603de6c7745b566a8b11aee12c 100644
--- a/tests/util.hpp
+++ b/tests/util.hpp
@@ -3,6 +3,7 @@
#include <string>
#include <vector>
#include <iomanip>
+#include <chrono>
#include <cmath>
@@ -14,6 +15,22 @@
namespace testing{
+using time_point = std::chrono::time_point<std::chrono::system_clock>;
+using duration_type = std::chrono::duration<double>;
+
+static inline
+time_point tic()
+{
+ return std::chrono::system_clock::now();
+}
+
+static inline
+double toc(time_point start)
+{
+ return duration_type(tic() - start).count();
+}
+
+
[[gnu::unused]] static
void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
{
@@ -111,7 +128,7 @@ operator<< (std::ostream& o, spike_comparison const& spikes)
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
+ spikes.max_relative_error()
);
return o << buffer;
}
diff --git a/tests/validate_ball_and_stick.cpp b/tests/validate_ball_and_stick.cpp
index eac29e2e84f20653acdb7555714f9ed6c6e3e555..a2aef0ada111aa7ecb8d150034a437016e987ce1 100644
--- a/tests/validate_ball_and_stick.cpp
+++ b/tests/validate_ball_and_stick.cpp
@@ -8,7 +8,7 @@
#include <json/src/json.hpp>
-// compares results with those generated by nrn/soma.py
+// compares results with those generated by nrn/ball_and_stick.py
TEST(ball_and_stick, neuron_baseline)
{
using namespace nest::mc;
@@ -27,7 +27,7 @@ TEST(ball_and_stick, neuron_baseline)
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
+ dendrite->mechanism("membrane").set("r_L", 100);
// add stimulus
cell.add_stimulus({1,1}, {5., 80., 0.3});
@@ -119,17 +119,173 @@ TEST(ball_and_stick, neuron_baseline)
}
// print results
+ auto colors = {memory::util::kWhite, memory::util::kGreen, memory::util::kYellow};
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();}
+ auto color = colors.begin();
+ for(auto const& result : r.comparisons) {
+ std::cout << std::setw(5) << r.n_comparments << " compartments : ";
+ std::cout << memory::util::colorize(util::pprintf("%\n", result), *(color++));
+ }
+ }
+
+ // 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()
);
- std::cout << std::setw(5) << r.n_comparments
- << " compartments : " << *m
- << "\n";
+ }
+ }
+
+ // 2. test that the best solution (i.e. with most compartments) matches the
+ // reference solution closely (less than 0.5% over the course of 100ms
+ // simulation)
+ auto tol = 0.5;
+ for(auto j=0; j<3; ++j) {
+ EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<tol);
+ }
+}
+
+// compares results with those generated by nrn/ball_and_3sticks.py
+TEST(ball_and_3stick, 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
+ std::vector<cable_segment*> dendrites;
+ dendrites.push_back(cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100));
+ dendrites.push_back(cell.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100));
+ dendrites.push_back(cell.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100));
+
+ for(auto dend : dendrites) {
+ dend->add_mechanism(pas_parameters());
+ dend->mechanism("membrane").set("r_L", 100);
+ }
+
+ // add stimulus
+ cell.add_stimulus({2,1}, {5., 80., 0.45});
+ cell.add_stimulus({3,1}, {40., 10.,-0.2});
+
+ // load data from file
+ auto cell_data = testing::load_spike_data("../nrn/ball_and_3stick.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;
+ auto start = testing::tic();
+ for(auto run_index=0u; run_index<cell_data.size(); ++run_index) {
+ auto& run = cell_data[run_index];
+ int num_compartments = run["nseg"];
+ for(auto dend : dendrites) {
+ dend->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,1.}, graph);
+ auto clamp_comp = nest::mc::find_compartment_index({2,1.}, 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( {num_compartments, dt, v, measurements} );
+ }
+ auto time_taken = testing::toc(start);
+ std::cout << "took " << time_taken << " seconds\n";
+
+ // print results
+ auto colors = {memory::util::kWhite, memory::util::kGreen, memory::util::kYellow};
+ for(auto& r : results){
+ auto color = colors.begin();
+ for(auto const& result : r.comparisons) {
+ std::cout << std::setw(5) << r.n_comparments << " compartments : ";
+ std::cout << memory::util::colorize(util::pprintf("%\n", result), *(color++));
+ }
}
// sort results in ascending order of compartments
@@ -153,10 +309,11 @@ TEST(ball_and_stick, neuron_baseline)
}
// 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
+ // reference solution closely (less than 0.5% over the course of 100ms
// simulation)
+ auto tol = 0.5;
for(auto j=0; j<3; ++j) {
- EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<0.1);
+ EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<tol);
}
}