diff --git a/nrn/ball_and_stick.py b/nrn/ball_and_stick.py
new file mode 100644
index 0000000000000000000000000000000000000000..ce857fc1a103ab20e00e2b2a60ade74ecffa7de9
--- /dev/null
+++ b/nrn/ball_and_stick.py
@@ -0,0 +1,84 @@
+from neuron import h, gui
+import numpy as np
+
+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
+# 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.
+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
+ 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
+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
+
+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
+
+ h.tstop = simdur
+ h.dt = dt[pos]
+ h.run()
+
+ pyplot.plot(t_vec, v_soma, somastyles[pos], linewidth=2)
+ pyplot.plot(t_vec, v_dend, dendstyles[pos], linewidth=2)
+
+pyplot.xlabel('time (ms)')
+pyplot.ylabel('mV')
+
+data = np.loadtxt('../tests/v.dat')
+t = data[:,0]
+nfields = data.shape[1]
+
+v = data[:,1]
+pyplot.plot(t, v, 'g', linewidth=2)
+v = data[:,2]
+pyplot.plot(t, v, 'g-.', linewidth=2)
+
+pyplot.show()
+
diff --git a/nrn/nest_soma.py b/nrn/nest_soma.py
new file mode 100644
index 0000000000000000000000000000000000000000..4acec3ed5cb1ac670abc58acc6099dc1dbf696ed
--- /dev/null
+++ b/nrn/nest_soma.py
@@ -0,0 +1,23 @@
+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()
+
diff --git a/nrn/soma.py b/nrn/soma.py
new file mode 100644
index 0000000000000000000000000000000000000000..f8912c4f7a0ea6bc204695cbef683b01e5172e44
--- /dev/null
+++ b/nrn/soma.py
@@ -0,0 +1,67 @@
+import os.path
+from neuron import h, gui
+from matplotlib import pyplot
+import numpy as np
+
+soma = h.Section(name='soma')
+
+soma.L = soma.diam = 18.8
+soma.Ra = 123
+
+print "Surface area of soma =", h.area(0.5, sec=soma)
+
+soma.insert('hh')
+
+stim = h.IClamp(soma(0.5))
+stim.delay = 10
+stim.dur = 100
+stim.amp = 0.1
+
+v_vec = h.Vector() # Membrane potential vector
+t_vec = h.Vector() # Time stamp vector
+v_vec.record(soma(0.5)._ref_v)
+t_vec.record(h._ref_t)
+simdur = 120
+
+# initialize plot
+pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
+
+pyplot.subplot(2,1,1)
+pyplot.grid()
+
+h.tstop = simdur
+
+# run neuron with multiple dt
+
+for dt in [0.02, 0.01, 0.005, 0.002, 0.001]:
+ 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()
+
diff --git a/src/fvm.hpp b/src/fvm.hpp
index c69f76289a07b6061c2586eb37a84115507575ce..19b27c45268f8be09a4d10e9a79adf2a5d714ac6 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -434,6 +434,8 @@ void fvm_cell<T, I>::setup_matrix(T dt)
template <typename T, typename I>
void fvm_cell<T, I>::initialize()
{
+ t_ = 0.;
+
// initialize mechanism states
for(auto& m : mechanisms_) {
m->nrn_init();
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 918e00d667a8014783d84473fd59776ad81705ac..4cc842363ad1474e408f7b44e9bb717197bb33fb 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -8,10 +8,10 @@ set(HEADERS
../src/tree.hpp
)
-set(TEST_SOURCES
- # google test framework
- gtest-all.cpp
+# google test framework
+add_library(gtest gtest-all.cpp gtest.h)
+set(TEST_SOURCES
# unit tests
test_algorithms.cpp
test_cell.cpp
@@ -26,14 +26,29 @@ set(TEST_SOURCES
test_tree.cpp
# unit test driver
- main.cpp
+ test.cpp
+)
+
+set(VALIDATION_SOURCES
+ # unit tests
+ validate_soma.cpp
+
+ # unit test driver
+ validate.cpp
)
add_executable(test.exe ${TEST_SOURCES} ${HEADERS})
+add_executable(validate.exe ${VALIDATION_SOURCES} ${HEADERS})
-target_link_libraries(test.exe LINK_PUBLIC cellalgo)
+target_link_libraries(test.exe LINK_PUBLIC cellalgo gtest)
+target_link_libraries(validate.exe LINK_PUBLIC cellalgo gtest)
set_target_properties(test.exe
PROPERTIES
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
)
+
+set_target_properties(validate.exe
+ PROPERTIES
+ RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
+)
diff --git a/tests/main.cpp b/tests/test.cpp
similarity index 100%
rename from tests/main.cpp
rename to tests/test.cpp
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 82ac6fd452039bfb7f3d703c4145f5bc8601a931..41f17ab5187c640093352eea88e20dae018c78a3 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -1,40 +1,18 @@
#include <fstream>
#include "gtest.h"
+#include "util.hpp"
#include "../src/cell.hpp"
#include "../src/fvm.hpp"
-void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
-{
- auto m = values.size();
- if(!m) return;
-
- std::ofstream fid(fname);
- if(!fid.is_open()) return;
-
- auto n = values[0].size();
- for(const auto& v : values) {
- std::cout << n << " ----- " << v.size() << "\n";
- if(n!=v.size()) {
- std::cerr << "all output arrays must have the same length\n";
- return;
- }
- }
-
- for(auto i=0u; i<n; ++i) {
- for(auto j=0u; j<m; ++j) {
- fid << " " << values[j][i];
- }
- fid << "\n";
- }
-}
-
// based on hh/Neuron/steps_A.py
TEST(run, single_compartment)
{
using namespace nest::mc;
+ std::vector<std::vector<double>> results(2);
+
nest::mc::cell cell;
// setup global state for the mechanisms
@@ -46,7 +24,8 @@ TEST(run, single_compartment)
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., 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);
@@ -57,31 +36,26 @@ TEST(run, single_compartment)
model.initialize(); // have to do this _after_ initial conditions are set
// run the simulation
- auto dt = 0.025; // ms
- auto tfinal = 25.; // ms
+ auto dt = 0.0005; // ms
+ auto tfinal = 120.; // ms
int nt = tfinal/dt;
- std::vector<double> result;
- result.push_back(model.voltage()[0]);
+ results[0].push_back(0.);
+ results[1].push_back(model.voltage()[0]);
for(auto i=0; i<nt; ++i) {
model.advance(dt);
- result.push_back(model.voltage()[0]);
+ // save voltage at soma
+ results[0].push_back((i+1)*dt);
+ results[1].push_back(model.voltage()[0]);
}
- {
- std::ofstream fid("v.dat");
- auto t = 0.;
- for(auto v:result) {
- fid << t << " " << v << "\n";
- t += dt;
- }
- }
+ write_vis_file("v.dat", results);
}
TEST(run, ball_and_stick)
{
using namespace nest::mc;
- std::vector<std::vector<double>> results(4);
+ std::vector<std::vector<double>> results(3);
nest::mc::cell cell;
@@ -94,12 +68,14 @@ TEST(run, ball_and_stick)
// 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->set_compartments(81); // 5 compartments
+ dendrite->set_compartments(10);
dendrite->add_mechanism(pas_parameters());
+ dendrite->mechanism("membrane").set("r_L", 100); // no effect for single compartment cell
+
// add stimulus
- //cell.add_stimulus({0,0}, {5., 3., 0.1}); // soma
- cell.add_stimulus({1,1}, {5., 30., 0.1});
+ //cell.add_stimulus({0,0}, {5., 1., 0.1}); // soma
+ cell.add_stimulus({1,1}, {5., 1., 0.5});
// make the lowered finite volume cell
fvm::fvm_cell<double, int> model(cell);
@@ -112,19 +88,18 @@ TEST(run, ball_and_stick)
model.initialize(); // have to do this _after_ initial conditions are set
// run the simulation
- auto dt = 0.02; // ms
+ auto pos = model.size()-1;
+ auto dt = 0.001; // ms
auto tfinal = 25.; // ms
int nt = tfinal/dt;
results[0].push_back(0.);
results[1].push_back(model.voltage()[0]);
- results[2].push_back(model.voltage()[10]);
- results[3].push_back(model.voltage()[20]);
+ results[2].push_back(model.voltage()[pos]);
for(auto i=0; i<nt; ++i) {
model.advance(dt);
results[0].push_back((i+1)*dt);
results[1].push_back(model.voltage()[0]);
- results[2].push_back(model.voltage()[10]);
- results[3].push_back(model.voltage()[20]);
+ results[2].push_back(model.voltage()[pos]);
}
write_vis_file("v.dat", results);
diff --git a/tests/test_tree.cpp b/tests/test_tree.cpp
index 9bb320d1e4cf5f1b0e8b1f5e320384f32935d728..2afcb1ed0e24a2bfb0123ba25abab393f59c9ec4 100644
--- a/tests/test_tree.cpp
+++ b/tests/test_tree.cpp
@@ -232,7 +232,7 @@ TEST(cell_tree, balance) {
EXPECT_EQ(t.parent(5), 2);
EXPECT_EQ(t.parent(6), 5);
- t.to_graphviz("cell.dot");
+ //t.to_graphviz("cell.dot");
}
}
@@ -246,8 +246,7 @@ TEST(cell_tree, json_load)
for(auto c : range(0,cell_data.size())) {
std::vector<int> parent_index = cell_data[c]["parent_index"];
cell_tree tree(parent_index);
- //tree.to_graphviz("cell_" + std::to_string(c) + ".dot");
- tree.to_graphviz("cell" + std::to_string(c) + ".dot");
+ //tree.to_graphviz("cell" + std::to_string(c) + ".dot");
}
}
diff --git a/tests/util.hpp b/tests/util.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cf10847fbf0de27657be11eb54b545538c3799b5
--- /dev/null
+++ b/tests/util.hpp
@@ -0,0 +1,30 @@
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+
+static
+void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
+{
+ auto m = values.size();
+ if(!m) return;
+
+ std::ofstream fid(fname);
+ if(!fid.is_open()) return;
+
+ auto n = values[0].size();
+ for(const auto& v : values) {
+ if(n!=v.size()) {
+ std::cerr << "all output arrays must have the same length\n";
+ return;
+ }
+ }
+
+ for(auto i=0u; i<n; ++i) {
+ for(auto j=0u; j<m; ++j) {
+ fid << " " << values[j][i];
+ }
+ fid << "\n";
+ }
+}
+
diff --git a/tests/validate.cpp b/tests/validate.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d814fca1512dd79b22763a4f9a769f2984ca5269
--- /dev/null
+++ b/tests/validate.cpp
@@ -0,0 +1,7 @@
+#include "gtest.h"
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}
+
diff --git a/tests/validate_soma.cpp b/tests/validate_soma.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b88c09f40b7426d54b06047cb290c993cb9083d5
--- /dev/null
+++ b/tests/validate_soma.cpp
@@ -0,0 +1,55 @@
+#include <fstream>
+
+#include "gtest.h"
+#include "util.hpp"
+
+#include "../src/cell.hpp"
+#include "../src/fvm.hpp"
+
+// based on hh/Neuron/steps_A.py
+TEST(soma, resolutions)
+{
+ 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., 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);
+
+ // 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;
+ results[0].push_back(0.);
+ results[1].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]);
+ }
+ write_vis_file("v_" + std::to_string(i) + ".dat", results);
+ ++i;
+ }
+}
+
diff --git a/tests/vis.py b/tests/vis.py
deleted file mode 100644
index 2a1834819da3f6d1c6e0865ec194c15ce637bb35..0000000000000000000000000000000000000000
--- a/tests/vis.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from matplotlib import pyplot
-import numpy as np
-
-data = np.loadtxt('v.dat')
-t = data[:,0]
-nfields = data.shape[1]
-
-for i in range(1,nfields) :
- v = data[:,i]
- pyplot.plot(t, v)
-
-pyplot.xlabel('time (ms)')
-pyplot.ylabel('mV')
-pyplot.grid()
-pyplot.show()
-
-
-#mi = min(v)
-#ma = max(v)
-#rng = (ma-mi)/2 * 1.1
-#mn = 0.5*(mi+ma)
-#pyplot.ylim(mn-rng, mn+rng)