diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 55d6af4d257348526713c1fdc381c4c2792651fd..c356012c50dfe37e3a849f16809e390433511306 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -45,6 +45,8 @@ flags = [
'.',
'-I',
'include',
+ '-I',
+ 'external',
# '-isystem',
# '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
# '-I',
diff --git a/docs/symbols.tex b/docs/symbols.tex
index e12b55820da048ffef27bb790fca4e0b363ba23e..baa371ad8efe42d483d2f35325b6d04ef730ca43 100644
--- a/docs/symbols.tex
+++ b/docs/symbols.tex
@@ -25,8 +25,9 @@ In \eq{eq:linsys_FV} we choose units of $mV \equiv 10^{-3}V$ for each term becau
$a,~\Delta x$ & $\mu m$ & $10^{-6} \cdot m$ & yes \\
$\sigma_{i},~\sigma_{ij}$ & $\mu m^2$ & $10^{-12} \cdot m^2$ & yes \\
$c_m$ & $F\cdot m^{-2}$ & $s\cdot A\cdot V^{-1}\cdot m^{-2}$ & no \\
- $r_L$ & $\Omega\cdot cm$ & $ 10^{-2} \cdot A^{-1}\cdot V\cdot m$ & yes \\
+ $r_L$ & $\Omega\cdot cm$ & $10^{-2} \cdot A^{-1}\cdot V\cdot m$ & yes \\
$\overline{g}$ & $S\cdot cm^{-2}$ & $10^{4} \cdot A\cdot V^{-1}\cdot m^{-2}$ & yes \\
+ $g_s$ & $\mu S$ & $10^{-6} \cdot A\cdot V^{-1}$ & yes \\
$I_e$ & $nA$ & $10^{-9} \cdot A$ & yes \\
\hline
\end{tabular}
@@ -43,6 +44,16 @@ Membrane current is calculated as follows $i_m = \overline{g}(E-V)$, with units
&= 10^{4} \cdot A\cdot V^{-1}\cdot m^{-2} \cdot 10^{-3} \cdot V \nonumber \\
&= 10 \cdot A \cdot m^{-2}. \label{eq:im_unit}
\end{align}
+The point process currents are calculated as point sources which must be turned into current densities as follows $i_m = g_s(E-V)/\sigma_i$.
+The units for the synaptic conductance $g_s$ are $\mu S$, so the units are calculated as follows
+\begin{align}
+ \unit{ i_m } &= \unit{ g_s } \unit{ V } \unit{\sigma_i}^-1 \nonumber \\
+ &= 10^{-6} \cdot A\cdot V^{-1} \cdot 10^{-3} \cdot 10^{12} \cdot m^{-2} \nonumber \\
+ &= 10^{3} \cdot A \cdot m^{-2}, \label{eq:ims_unit}
+\end{align}
+which must be scaled by $10^{2}$ to match that of of the density channels in \eq{eq:im_unit}.
+
+
The injected current $I_e$ has units $nA$, which has to be expressed in terms of current per unit area $i_e=I_e / \sigma_i$ with units
\begin{align}
\unit{ i_e } &= \unit{ I_e } \unit{ \sigma_i }^{-1} \nonumber \\
diff --git a/external/modparser b/external/modparser
index d1dd8fa868c0dc7e348093ed4a4b5c6c94e29e34..b48972977a590b90cef6a142cda570984d79cb9e 160000
--- a/external/modparser
+++ b/external/modparser
@@ -1 +1 @@
-Subproject commit d1dd8fa868c0dc7e348093ed4a4b5c6c94e29e34
+Subproject commit b48972977a590b90cef6a142cda570984d79cb9e
diff --git a/external/vector b/external/vector
index 8153d030abba20ddd174da9804406139c393a808..193df3683a3eaa2bf203bbcfa935d7e7019b1719 160000
--- a/external/vector
+++ b/external/vector
@@ -1 +1 @@
-Subproject commit 8153d030abba20ddd174da9804406139c393a808
+Subproject commit 193df3683a3eaa2bf203bbcfa935d7e7019b1719
diff --git a/include/mechanisms/hh_verb.hpp b/include/mechanisms/hh_verb.hpp
deleted file mode 100644
index 8b67a69d56a0f61b2db8baaa37d03fc34d7e9687..0000000000000000000000000000000000000000
--- a/include/mechanisms/hh_verb.hpp
+++ /dev/null
@@ -1,284 +0,0 @@
-#pragma once
-
-#include <cmath>
-#include <limits>
-
-#include <mechanism.hpp>
-#include <mechanism_interface.hpp>
-#include <algorithms.hpp>
-
-namespace nest{ namespace mc{ namespace mechanisms{ namespace hh{
-
-template<typename T, typename I>
-class mechanism_hh : public mechanism<T, I> {
-public:
- using base = mechanism<T, I>;
- using value_type = typename base::value_type;
- using size_type = typename base::size_type;
- using vector_type = typename base::vector_type;
- using view_type = typename base::view_type;
- using index_type = typename base::index_type;
- using index_view = typename index_type::view_type;
- using indexed_view_type= typename base::indexed_view_type;
- using ion_type = typename base::ion_type;
-
- struct Ionna {
- view_type ena;
- view_type ina;
- index_type index;
- std::size_t memory() const { return sizeof(size_type)*index.size(); }
- std::size_t size() const { return index.size(); }
- };
- Ionna ion_na;
- struct Ionk {
- view_type ek;
- view_type ik;
- index_type index;
- std::size_t memory() const { return sizeof(size_type)*index.size(); }
- std::size_t size() const { return index.size(); }
- };
- Ionk ion_k;
-
- mechanism_hh(view_type vec_v, view_type vec_i, index_view node_index)
- : base(vec_v, vec_i, node_index)
- {
- size_type num_fields = 15;
-
- // calculate the padding required to maintain proper alignment of sub arrays
- auto alignment = data_.alignment();
- auto field_size_in_bytes = sizeof(value_type)*size();
- auto remainder = field_size_in_bytes % alignment;
- auto padding = remainder ? (alignment - remainder)/sizeof(value_type) : 0;
- auto field_size = size()+padding;
-
- // allocate memory
- data_ = vector_type(field_size * num_fields);
- data_(memory::all) = std::numeric_limits<value_type>::quiet_NaN();
-
- // asign the sub-arrays
- gnabar = data_(0*field_size, 1*size());
- minf = data_(1*field_size, 2*size());
- h = data_(2*field_size, 3*size());
- m = data_(3*field_size, 4*size());
- gl = data_(4*field_size, 5*size());
- gkbar = data_(5*field_size, 6*size());
- el = data_(6*field_size, 7*size());
- ninf = data_(7*field_size, 8*size());
- mtau = data_(8*field_size, 9*size());
- gna = data_(9*field_size, 10*size());
- gk = data_(10*field_size, 11*size());
- n = data_(11*field_size, 12*size());
- hinf = data_(12*field_size, 13*size());
- ntau = data_(13*field_size, 14*size());
- htau = data_(14*field_size, 15*size());
-
- // set initial values for variables and parameters
- std::fill(gnabar.data(), gnabar.data()+size(), 0.12);
- std::fill(gkbar.data(), gkbar.data()+size(), 0.036);
- std::fill(gl.data(), gl.data()+size(), 0.0003);
- std::fill(el.data(), el.data()+size(), -54.3);
- }
-
- using base::size;
-
- std::size_t memory() const override {
- auto s = std::size_t{0};
- s += data_.size()*sizeof(value_type);
- s += ion_na.memory();
- s += ion_k.memory();
- return s;
- }
-
- void set_params(value_type t_, value_type dt_) override {
- t = t_;
- dt = dt_;
- }
-
- std::string name() const override {
- return "hh";
- }
-
- mechanismKind kind() const override {
- return mechanismKind::density;
- }
-
- bool uses_ion(ionKind k) const override {
- switch(k) {
- case ionKind::na : return true;
- case ionKind::ca : return false;
- case ionKind::k : return true;
- }
- return false;
- }
-
- void set_ion(ionKind k, ion_type& i) override {
- using nest::mc::algorithms::index_into;
- if(k==ionKind::na) {
- ion_na.index = index_into(i.node_index(), node_index_);
- ion_na.ina = i.current();
- ion_na.ena = i.reversal_potential();
- return;
- }
- if(k==ionKind::k) {
- ion_k.index = index_into(i.node_index(), node_index_);
- ion_k.ik = i.current();
- ion_k.ek = i.reversal_potential();
- return;
- }
- throw std::domain_error(nest::mc::util::pprintf("mechanism % does not support ion type\n", name()));
- }
-
- void nrn_state() {
- const indexed_view_type vec_v(vec_v_, node_index_);
- int n_ = node_index_.size();
- for(int i_=0; i_<n_; ++i_) {
- value_type v = vec_v[i_];
- rates(i_, v);
- m[i_] = minf[i_]+(m[i_]-minf[i_])*exp( -dt/mtau[i_]);
- h[i_] = hinf[i_]+(h[i_]-hinf[i_])*exp( -dt/htau[i_]);
- n[i_] = ninf[i_]+(n[i_]-ninf[i_])*exp( -dt/ntau[i_]);
- printf("m h n : %18.14f %18.14f %18.14f\n", m[i_], h[i_], n[i_]);
- }
- }
-
- void rates(const int i_, value_type v) {
- value_type ll3_, ll1_, ll0_, alpha, beta, ll2_, sum, q10;
- q10 = std::pow( 3, (celsius- 6.2999999999999998)/ 10);
- printf("q10 %18.14f\n", q10);
- ll2_ = -(v+ 40);
- ll0_ = ll2_/(exp(ll2_/ 10)- 1);
- alpha = 0.10000000000000001*ll0_;
- beta = 4*exp( -(v+ 65)/ 18);
-
- //std::cout << "v " << v << " dt " << dt << "\n";
- //std::cout << "m (alpha, beta) : " << alpha << " " << beta << "\n";
-
- sum = alpha+beta;
- mtau[i_] = 1/(q10*sum);
- minf[i_] = alpha/sum;
-
- //////////////////////////////////////////////////////////
-
- alpha = 0.070000000000000007*exp( -(v+ 65)/ 20);
- beta = 1/(exp( -(v+ 35)/ 10)+ 1);
-
- //std::cout << "h (alpha, beta) : " << alpha << " " << beta << "\n";
-
- sum = alpha+beta;
- htau[i_] = 1/(q10*sum);
- hinf[i_] = alpha/sum;
-
- //////////////////////////////////////////////////////////
-
- //ll3_ = -v+ 55; // TODO : inlining is breaking
- ll3_ = -(v+ 55);
- ll1_ = ll3_/(exp(ll3_/ 10)- 1);
- alpha = 0.01*ll1_;
- beta = 0.125*exp( -(v+ 65)/ 80);
-
- //std::cout << "n (alpha, beta) : " << alpha << " " << beta << "\n";
-
- //printf("m_inf, h_inf, n_inf : %18.16f %18.16f %18.16f\n",
- //minf[i_], hinf[i_], ninf[i_]);
-
- sum = alpha+beta;
- ntau[i_] = 1/(q10*sum); // TODO : make modparser insert parenthesis
- ninf[i_] = alpha/sum;
- }
-
- void nrn_current() {
- const indexed_view_type ion_ek(ion_k.ek, ion_k.index);
- indexed_view_type ion_ina(ion_na.ina, ion_na.index);
- const indexed_view_type ion_ena(ion_na.ena, ion_na.index);
- indexed_view_type ion_ik(ion_k.ik, ion_k.index);
- indexed_view_type vec_i(vec_i_, node_index_);
- const indexed_view_type vec_v(vec_v_, node_index_);
- int n_ = node_index_.size();
- for(int i_=0; i_<n_; ++i_) {
- value_type ek = ion_ek[i_];
- value_type ena = ion_ena[i_];
- value_type v = vec_v[i_];
- value_type il, ina, ik, current_;
- gna[i_] = gnabar[i_]*m[i_]*m[i_]*m[i_]*h[i_];
- ina = gna[i_]*(v-ena);
- current_ = ina;
- gk[i_] = gkbar[i_]*n[i_]*n[i_]*n[i_]*n[i_];
- ik = gk[i_]*(v-ek);
- current_ = current_+ik;
- il = gl[i_]*(v-el[i_]);
- current_ = current_+il;
- ion_ina[i_] += ina;
- ion_ik[i_] += ik;
- vec_i[i_] += current_;
- printf("i = (l+k+na) %18.16f = %18.16f %18.16f %18.16f\n",
- current_, il, ik, ina);
- }
- }
-
- void nrn_init() {
- const indexed_view_type vec_v(vec_v_, node_index_);
- int n_ = node_index_.size();
- for(int i_=0; i_<n_; ++i_) {
- value_type v = vec_v[i_];
- rates(i_, v);
- m[i_] = minf[i_];
- h[i_] = hinf[i_];
- n[i_] = ninf[i_];
- printf("initial conditions for m, h, n : %16.14f %16.14f %16.14f\n", m[0], h[0], n[0]);
- }
- }
-
- vector_type data_;
- view_type gnabar;
- view_type minf;
- view_type h;
- view_type m;
- view_type gl;
- view_type gkbar;
- view_type el;
- view_type ninf;
- view_type mtau;
- view_type gna;
- view_type gk;
- view_type n;
- view_type hinf;
- view_type ntau;
- view_type htau;
- value_type t = std::numeric_limits<value_type>::quiet_NaN();
- value_type dt = std::numeric_limits<value_type>::quiet_NaN();
- value_type celsius = 6.3; // TODO change from 37
-
- using base::vec_v_;
- using base::vec_i_;
- using base::node_index_;
-
-};
-
-template<typename T, typename I>
-struct helper : public mechanism_helper<T, I> {
- using base = mechanism_helper<T, I>;
- using index_view = typename base::index_view;
- using view_type = typename base::view_type;
- using mechanism_ptr_type = typename base::mechanism_ptr_type;
- using mechanism_type = mechanism_hh<T, I>;
-
- std::string
- name() const override
- {
- return "hh";
- }
-
- mechanism_ptr<T,I>
- new_mechanism(view_type vec_v, view_type vec_i, index_view node_index) const override
- {
- return nest::mc::mechanisms::make_mechanism<mechanism_type>(vec_v, vec_i, node_index);
- }
-
- void
- set_parameters(mechanism_ptr_type&, parameter_list const&) const override
- {
- }
-
-};
-
-}}}} // namespaces
diff --git a/mechanisms/generate.sh b/mechanisms/generate.sh
index 94755fbb91b912e5eada0148af1c72ea91945ecd..1fcd49beb79603e6d79cf35e1e4b9f699fcb2f49 100755
--- a/mechanisms/generate.sh
+++ b/mechanisms/generate.sh
@@ -1,4 +1,8 @@
-for mech in pas hh
+#flags="-t cpu -O"
+flags="-t cpu"
+
+for mech in pas hh expsyn
do
- ../external/modparser/bin/modcc -t cpu -o ../include/mechanisms/$mech.hpp ./mod/$mech.mod
+ echo ../external/modparser/bin/modcc ${flags} -o ../include/mechanisms/$mech.hpp ./mod/$mech.mod
+ ../external/modparser/bin/modcc ${flags} -o ../include/mechanisms/$mech.hpp ./mod/$mech.mod
done
diff --git a/mechanisms/mod/expsyn.mod b/mechanisms/mod/expsyn.mod
new file mode 100755
index 0000000000000000000000000000000000000000..db1f6d6ed61a8b563141ea2516658b8f87a35637
--- /dev/null
+++ b/mechanisms/mod/expsyn.mod
@@ -0,0 +1,41 @@
+NEURON {
+ POINT_PROCESS ExpSyn
+ RANGE tau, e
+ NONSPECIFIC_CURRENT i
+}
+
+UNITS {
+ (nA) = (nanoamp)
+ (mV) = (millivolt)
+ (uS) = (microsiemens)
+}
+
+PARAMETER {
+ : the default for Neuron is 0.1
+ :tau = 0.1 (ms) : <1e-9,1e9>
+ tau = 2.0 (ms) : <1e-9,1e9>
+ e = 0 (mV)
+}
+
+ASSIGNED {}
+
+STATE {
+ g : (uS)
+}
+
+INITIAL {
+ g=0
+}
+
+BREAKPOINT {
+ SOLVE state METHOD cnexp
+ i = g*(v - e)
+}
+
+DERIVATIVE state {
+ g' = -g/tau
+}
+
+NET_RECEIVE(weight) {
+ g = g + weight
+}
diff --git a/nrn/generate_validation.sh b/nrn/generate_validation.sh
index f35d3521e9d2ce785a01bc14cd966d4910ecf7f7..c42b1d5ce66512b8dc5ff46bf13e85f4707bd8ef 100755
--- a/nrn/generate_validation.sh
+++ b/nrn/generate_validation.sh
@@ -1,3 +1,4 @@
python2.7 ./soma.py
python2.7 ./ball_and_stick.py
python2.7 ./ball_and_3stick.py
+python2.7 ./simple_synapse.py
diff --git a/nrn/simple_synapse.py b/nrn/simple_synapse.py
new file mode 100644
index 0000000000000000000000000000000000000000..1779799f716157df995f2321346ad6ff52f1a608
--- /dev/null
+++ b/nrn/simple_synapse.py
@@ -0,0 +1,140 @@
+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))
+
+# 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
+
+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
+
+# create a stimulator that stimulates once at 9 ms
+stim = h.NetStim()
+stim.number = 1
+stim.start = 9
+
+if args.plot :
+ pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
+ pyplot.grid()
+
+simdur = 50.0
+h.tstop = simdur
+h.dt = 0.01
+
+start = timer()
+results = []
+for nseg in [5, 11, 51, 101] :
+
+ print 'simulation with ', nseg, ' compartments in dendrite...'
+
+ dend.nseg=nseg
+
+ # add a synapse
+ syn_ = h.ExpSyn(dend(0.5))
+ syn_.tau = 2
+ ncstim = h.NetCon(stim, syn_)
+ ncstim.delay = 1 # 1 ms delay (arrives at 10ms)
+ ncstim.weight[0] = 0.04 # NetCon weight is a vector
+
+ ncstim1 = h.NetCon(stim, syn_)
+ ncstim1.delay = 11 # 1 ms delay (arrives at 10ms)
+ ncstim1.weight[0] = 0.04 # NetCon weight is a vector
+ ncstim2 = h.NetCon(stim, syn_)
+ ncstim2.delay = 31 # 1 ms delay (arrives at 10ms)
+ ncstim2.weight[0] = 0.04 # NetCon weight is a vector
+
+ # record voltages
+ v_soma = h.Vector() # soma
+ v_dend = h.Vector() # middle of dendrite
+
+ v_soma.record( soma(0.5)._ref_v)
+ v_dend.record( dend(0.5)._ref_v)
+
+ # record spikes
+ spike_counter_soma = h.APCount(soma(0.5))
+ spike_counter_soma.thresh = 20
+ spike_counter_dend = h.APCount(dend(0.5))
+ spike_counter_dend.thresh = -10
+
+ spikes_soma = h.Vector() # soma
+ spikes_dend = h.Vector() # middle of dendrite
+
+ spike_counter_soma.record(spikes_soma)
+ spike_counter_dend.record(spikes_dend)
+
+ # 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()
+ },
+ }
+ }
+ )
+
+ 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))
+
+# time the simulations
+end = timer()
+print "took ", end-start, " seconds"
+
+# save the spike info as in json format
+fp = open('simple_synapse.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/src/cell.cpp b/src/cell.cpp
index 8fdd1bb6e0f55cdfa2a3c7ea6632e0fdee74bba6..687831570cd0b5d332d983982305c8ec127cae47 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -195,6 +195,16 @@ std::vector<int> const& cell::segment_parents() const
return parents_;
}
+void cell::add_synapse(segment_location loc)
+{
+ synapses_.push_back(loc);
+}
+
+const std::vector<segment_location>& cell::synapses() const
+{
+ return synapses_;
+}
+
// 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
diff --git a/src/cell.hpp b/src/cell.hpp
index fece134947041d705564008e54be61ce93d50f60..7b21909166a6e8d3dd4826b8dcf7ec8a125706de 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -111,6 +111,11 @@ class cell {
return stimulii_;
}
+ void add_synapse(segment_location loc);
+
+ const std::vector<segment_location>& synapses() const;
+
+
private:
// storage for connections
@@ -121,6 +126,9 @@ class cell {
// the stimulii
std::vector<std::pair<segment_location, i_clamp>> stimulii_;
+
+ // the synapses
+ std::vector<segment_location> synapses_;
};
// Checks that two cells have the same
diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp
index 43fac2429a657545251159a3fa97002b7243fe86..785a4404e104e0f58db1e361b25ec13c81295c5c 100644
--- a/src/cell_tree.hpp
+++ b/src/cell_tree.hpp
@@ -32,9 +32,10 @@ class cell_tree {
public :
// use a signed 16-bit integer for storage of indexes, which is reasonable given
// that typical cells have at most 1000-2000 segments
- using int_type = int16_t;
- using index_type = memory::HostVector<int_type>;
- using index_view = index_type::view_type;
+ using int_type = int16_t;
+ using index_type = memory::HostVector<int_type>;
+ using view_type = index_type::view_type;
+ using const_view_type = index_type::const_view_type;
/// default empty constructor
cell_tree() = default;
@@ -133,7 +134,7 @@ public :
}
/// returns a list of the children of segment b
- const index_view children(size_t b) const {
+ const_view_type children(size_t b) const {
return tree_.children(b);
}
diff --git a/src/event_queue.hpp b/src/event_queue.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a8b33707a2d82e949717fd243479ce33776120e
--- /dev/null
+++ b/src/event_queue.hpp
@@ -0,0 +1,74 @@
+#pragma once
+
+#include <ostream>
+#include <queue>
+
+#include <cstdint>
+
+namespace nest {
+namespace mc {
+
+struct local_event {
+ uint32_t target;
+ float time;
+ float weight;
+};
+
+inline bool operator < (local_event const& l, local_event const& r) {
+ return l.time < r.time;
+}
+
+inline bool operator > (local_event const& l, local_event const& r) {
+ return l.time > r.time;
+}
+
+class event_queue {
+public :
+ // create
+ event_queue() {}
+
+ // push stuff
+ template <typename Iter>
+ void push(Iter b, Iter e) {
+ for (; b!=e; ++b) {
+ queue_.push(*b);
+ }
+ }
+
+ // push thing
+ void push(local_event e) {
+ queue_.push(e);
+ }
+
+ std::size_t size() const {
+ return queue_.size();
+ }
+
+ // pop until
+ std::pair<bool, local_event> pop_if_before(float t_until) {
+ if (!queue_.empty() && queue_.top().time < t_until) {
+ auto ev = queue_.top();
+ queue_.pop();
+ return {true, ev};
+ }
+ else {
+ return {false, {}};
+ }
+ }
+
+private:
+ std::priority_queue<
+ local_event,
+ std::vector<local_event>,
+ std::greater<local_event>
+ > queue_;
+};
+
+} // namespace nest
+} // namespace mc
+
+inline
+std::ostream& operator<< (std::ostream& o, nest::mc::local_event e)
+{
+ return o << "event[" << e.target << "," << e.time << "," << e.weight << "]";
+}
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 69d23d0d1fe97e341ad8de3dc25ce8c0771ef78c..600ff04d4faecd6fc2a854825b5663652623e8eb 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -2,22 +2,26 @@
#include <algorithm>
#include <map>
+#include <set>
#include <string>
#include <vector>
#include <algorithms.hpp>
#include <cell.hpp>
+#include <event_queue.hpp>
#include <ion.hpp>
#include <math.hpp>
#include <matrix.hpp>
#include <mechanism.hpp>
#include <mechanism_interface.hpp>
-#include <util.hpp>
#include <segment.hpp>
#include <stimulus.hpp>
+#include <util.hpp>
#include <vector/include/Vector.hpp>
+#include <include/mechanisms/expsyn.hpp>
+
namespace nest {
namespace mc {
namespace fvm {
@@ -45,11 +49,13 @@ class fvm_cell {
using index_type = memory::HostVector<size_type>;
/// view into index container
using index_view = typename index_type::view_type;
+ using const_index_view = typename index_type::const_view_type;
/// the container used for values
using vector_type = memory::HostVector<value_type>;
/// view into value container
using vector_view = typename vector_type::view_type;
+ using const_vector_view = typename vector_type::const_view_type;
/// constructor
fvm_cell(nest::mc::cell const& cell);
@@ -57,26 +63,23 @@ class fvm_cell {
/// build the matrix for a given time step
void setup_matrix(value_type dt);
- /// TODO this should be const
/// which requires const_view in the vector library
- matrix_type& jacobian() {
+ const matrix_type& jacobian() {
return matrix_;
}
- /// TODO this should be const
/// return list of CV areas in :
/// um^2
/// 1e-6.mm^2
/// 1e-8.cm^2
- vector_view cv_areas() {
+ const_vector_view cv_areas() const {
return cv_areas_;
}
- /// TODO this should be const
/// return the capacitance of each CV surface
/// this is the total capacitance, not per unit area,
/// i.e. equivalent to sigma_i * c_m
- vector_view cv_capacitance() {
+ const_vector_view cv_capacitance() const {
return cv_capacitance_;
}
@@ -84,6 +87,9 @@ class fvm_cell {
vector_view voltage() {
return voltage_;
}
+ const_vector_view voltage() const {
+ return voltage_;
+ }
std::size_t size() const {
return matrix_.size();
@@ -130,9 +136,16 @@ class fvm_cell {
/// make a time step
void advance(value_type dt);
+ /// advance solution to target time tfinal with maximum step size dt
+ void advance_to(value_type tfinal, value_type dt);
+
/// set initial states
void initialize();
+ event_queue& queue() {
+ return events_;
+ }
+
private:
/// current time
@@ -161,6 +174,11 @@ class fvm_cell {
/// the potential in mV in each CV
vector_type voltage_;
+ /// synapses
+ using synapse_type =
+ mechanisms::ExpSyn::mechanism_ExpSyn<value_type, size_type>;
+ std::size_t synapse_index_;
+
/// the set of mechanisms present in the cell
std::vector<mechanism_type> mechanisms_;
@@ -168,6 +186,9 @@ class fvm_cell {
std::map<mechanisms::ionKind, ion_type> ions_;
std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
+
+ /// event queue
+ event_queue events_;
};
////////////////////////////////////////////////////////////////////////////////
@@ -293,7 +314,7 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// build a vector of the indexes of the compartments that contain
// the mechanism
- std::vector<int> compartment_index(num_comp);
+ index_type compartment_index(num_comp);
auto pos = 0u;
for(auto seg : mech.second) {
auto seg_size = segment_index[seg+1] - segment_index[seg];
@@ -306,41 +327,30 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
}
// instantiate the mechanism
- index_view node_index(compartment_index.data(), compartment_index.size());
mechanisms_.push_back(
- helper->new_mechanism(voltage_, current_, node_index)
+ helper->new_mechanism(voltage_, current_, compartment_index)
);
}
/////////////////////////////////////////////
// build the ion species
- // FIXME : this should be a private member function
/////////////////////////////////////////////
for(auto ion : mechanisms::ion_kinds()) {
- auto indexes =
- std::make_pair(std::vector<int>(size()), std::vector<int>(size()));
- auto ends =
- std::make_pair(indexes.first.begin(), indexes.second.begin());
-
- // after the loop the range
- // [indexes.first.begin(), ends.first)
- // will hold the indexes of the compartments that require ion
+ // find the compartment indexes of all compartments that have a
+ // mechanism that depends on/influences ion
+ std::set<int> index_set;
for(auto& mech : mechanisms_) {
if(mech->uses_ion(ion)) {
- ends.second =
- std::set_union(
- mech->node_index().begin(), mech->node_index().end(),
- indexes.first.begin(), ends.first,
- indexes.second.begin()
- );
- std::swap(indexes.first, indexes.second);
- std::swap(ends.first, ends.second);
+ for(auto idx : mech->node_index()) {
+ index_set.insert(idx);
+ }
}
}
+ std::vector<int> indexes(index_set.begin(), index_set.end());
// create the ion state
- if(auto n = std::distance(indexes.first.begin(), ends.first)) {
- ions_.emplace(ion, index_view(indexes.first.data(), n));
+ if(indexes.size()) {
+ ions_.emplace(ion, index_type(indexes));
}
// join the ion reference in each mechanism into the cell-wide ion state
@@ -356,7 +366,7 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// the default values in Neuron.
// Neuron's defaults are defined in the file
// nrn/src/nrnoc/membdef.h
- using memory::all;
+ auto all = memory::all;
constexpr value_type DEF_vrest = -65.0; // same name as #define in Neuron
@@ -377,6 +387,24 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
auto idx = find_compartment_index(stim.first, graph);
stimulii_.push_back( {idx, stim.second} );
}
+
+ // add the synapses
+ std::vector<size_type> synapse_indexes;
+ synapse_indexes.reserve(cell.synapses().size());
+ for(auto loc : cell.synapses()) {
+ synapse_indexes.push_back(
+ find_compartment_index(loc, graph)
+ );
+ }
+
+ mechanisms_.push_back(
+ mechanisms::make_mechanism<synapse_type>(
+ voltage_, current_, index_view(synapse_indexes)
+ )
+ );
+ synapse_index_ = mechanisms_.size()-1;
+ // don't forget to give point processes access to cv_areas_
+ mechanisms_[synapse_index_]->set_areas(cv_areas_);
}
template <typename T, typename I>
@@ -406,8 +434,6 @@ void fvm_cell<T, I>::setup_matrix(T dt)
//d(all) = 1.0;
d(all) = cv_areas_;
for(auto i=1u; i<d.size(); ++i) {
- // TODO get this right
- // probably requires scaling a by cv_areas_[i] and cv_areas_[p[i]]
auto a = 1e5*dt * face_alpha_[i];
d[i] += a;
@@ -417,7 +443,6 @@ void fvm_cell<T, I>::setup_matrix(T dt)
// add contribution to the diagonal of parent
d[p[i]] += a;
}
- //std::cout << "d " << d << " l " << l << " u " << u << "\n";
// the RHS of the linear system is
// V[i] - dt/cm*(im - ie)
@@ -461,29 +486,44 @@ void fvm_cell<T, I>::advance(T dt)
current_[loc] -= 100.*ie/cv_areas_[loc];
}
- //std::cout << "t " << t_ << " current " << current_;
-
// set matrix diagonals and rhs
setup_matrix(dt);
- //printf("rhs %18.14f d %18.14f\n", matrix_.rhs()[0], matrix_.d()[0]);
-
// solve the linear system
matrix_.solve();
voltage_(all) = matrix_.rhs();
- //printf("v solve %18.14f\n", voltage_[0]);
-
- //std::cout << " v " << voltage_ << "\n";
-
// update states
for(auto& m : mechanisms_) {
m->nrn_state();
}
t_ += dt;
- //std::cout << "******************\n";
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::advance_to(T tfinal, T dt)
+{
+ if(t_>=tfinal) {
+ return;
+ }
+
+ do {
+ auto tnext = std::min(tfinal, t_+dt);
+ auto next = events_.pop_if_before(tnext);
+ // if there is an event before tnext...
+ if(next.first) {
+ tnext = next.second.time;
+ }
+ advance(tnext-t_);
+ t_ = tnext;
+ if(next.first) { // handle event
+ auto &e = next.second;
+
+ mechanisms_[synapse_index_]->net_receive(e.target, e.weight);
+ }
+ } while(t_<tfinal);
}
} // namespace fvm
diff --git a/src/matrix.hpp b/src/matrix.hpp
index db9d2df7f56e8bfffb73051e996acb13d9fa0369..615ace9d1c3f06a7c223c35d98f5a484b2a0868c 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -19,10 +19,12 @@ class matrix {
using size_type = I;
// define storage types
- using vector_type = memory::HostVector<value_type>;
- using vector_view_type = typename vector_type::view_type;
- using index_type = memory::HostVector<size_type>;
- using index_view_type = typename index_type::view_type;
+ using vector_type = memory::HostVector<value_type>;
+ using vector_view_type = typename vector_type::view_type;
+ using const_vector_view_type = typename vector_type::const_view_type;
+ using index_type = memory::HostVector<size_type>;
+ using index_view_type = typename index_type::view_type;
+ using const_index_view_type = typename index_type::const_view_type;
matrix() = default;
@@ -85,41 +87,62 @@ class matrix {
vector_view_type vec_d() { return d(); }
vector_view_type vec_v() { return v(); }
+ const_vector_view_type vec_rhs() const { return rhs(); }
+ const_vector_view_type vec_d() const { return d(); }
+ const_vector_view_type vec_v() const { return v(); }
+
/// the vector holding the lower part of the matrix
- vector_view_type l()
- {
+ vector_view_type l() {
+ return l_;
+ }
+
+ const_vector_view_type l() const {
return l_;
}
/// the vector holding the diagonal of the matrix
- vector_view_type d()
- {
+ vector_view_type d() {
+ return d_;
+ }
+
+ const vector_view_type d() const {
return d_;
}
/// the vector holding the upper part of the matrix
- vector_view_type u()
- {
+ vector_view_type u() {
+ return u_;
+ }
+ const vector_view_type u() const {
return u_;
}
/// the vector holding the right hand side of the linear equation system
- vector_view_type rhs()
- {
+ vector_view_type rhs() {
+ return rhs_;
+ }
+
+ const_vector_view_type rhs() const {
return rhs_;
}
/// the vector holding the solution (voltage)
- vector_view_type v()
- {
+ vector_view_type v() {
EXPECTS(has_voltage_);
+ return v_;
+ }
+ const_vector_view_type v() const {
+ EXPECTS(has_voltage_);
return v_;
}
/// the vector holding the parent index
- index_view_type p()
- {
+ index_view_type p() {
+ return parent_index_;
+ }
+
+ const_index_view_type p() const {
return parent_index_;
}
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
index ea69c292e1f5c50efe0fd0e6e7dbf092c32b16fb..1afd8c31630faa017e95d8a7f48d444a5da93f2e 100644
--- a/src/mechanism.hpp
+++ b/src/mechanism.hpp
@@ -37,6 +37,7 @@ public:
: vec_v_(vec_v)
, vec_i_(vec_i)
, node_index_(node_index)
+ , vec_area_(nullptr, 0)
{}
std::size_t size() const
@@ -65,14 +66,20 @@ public:
virtual void nrn_init() = 0;
virtual void nrn_state() = 0;
virtual void nrn_current() = 0;
+ virtual void net_receive(int, value_type) {};
virtual bool uses_ion(ionKind) const = 0;
virtual void set_ion(ionKind k, ion_type& i) = 0;
+ void set_areas(view_type area) {
+ vec_area_ = area;
+ }
+
virtual mechanismKind kind() const = 0;
view_type vec_v_;
view_type vec_i_;
index_type node_index_;
+ view_type vec_area_;
};
template <typename T, typename I>
@@ -89,5 +96,5 @@ make_mechanism(
}
} // namespace mechanisms
-} // namespace nest
} // namespace mc
+} // namespace nest
diff --git a/src/tree.hpp b/src/tree.hpp
index d517f7e4570fc9278bab22ff974a3baaa964ea50..122d12e5951effa03fd3b97377da20c6c8f88453 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -21,7 +21,8 @@ class tree {
using int_type = int16_t;
using index_type = memory::HostVector<int_type>;
- using index_view = index_type::view_type;
+ using view_type = index_type::view_type;
+ using const_view_type = index_type::const_view_type;
tree() = default;
@@ -164,22 +165,22 @@ class tree {
}
/// return the child index
- const index_view child_index() const {
+ const_view_type child_index() {
return child_index_;
}
/// return the list of all children
- const index_view children() const {
+ const_view_type children() const {
return children_;
}
/// return the list of all children of branch b
- const index_view children(size_t b) const {
+ const_view_type children(size_t b) const {
return children_(child_index_[b], child_index_[b+1]);
}
/// return the list of parents
- const index_view parents() const {
+ const_view_type parents() const {
return parents_;
}
@@ -289,7 +290,7 @@ class tree {
int new_node,
int old_node,
int parent_node,
- index_view p,
+ view_type p,
tree const& old_tree
)
{
@@ -351,9 +352,9 @@ class tree {
// provide default parameters so that tree type can
// be default constructed
- index_view children_ = data_(0, 0);
- index_view child_index_= data_(0, 0);
- index_view parents_ = data_(0, 0);
+ view_type children_ = data_(0, 0);
+ view_type child_index_= data_(0, 0);
+ view_type parents_ = data_(0, 0);
};
template <typename C>
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index cb33af4d972dd398952b1879cca764c502d8f42c..88862afa769fd9165bcf1d25a0096a867ffd66f3 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -16,6 +16,7 @@ set(TEST_SOURCES
test_algorithms.cpp
test_cell.cpp
test_compartments.cpp
+ test_event_queue.cpp
test_fvm.cpp
test_matrix.cpp
test_mechanisms.cpp
@@ -24,6 +25,7 @@ set(TEST_SOURCES
test_segment.cpp
test_stimulus.cpp
test_swcio.cpp
+ test_synapses.cpp
test_tree.cpp
# unit test driver
@@ -32,8 +34,9 @@ set(TEST_SOURCES
set(VALIDATION_SOURCES
# unit tests
- validate_soma.cpp
validate_ball_and_stick.cpp
+ validate_soma.cpp
+ validate_synapses.cpp
# unit test driver
validate.cpp
diff --git a/tests/den_syn.txt b/tests/den_syn.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4629823aa311380d840911e3266f1749edfaa418
--- /dev/null
+++ b/tests/den_syn.txt
@@ -0,0 +1,1002 @@
+ 0 -65 -65
+ 0.05 -64.9987 -65
+ 0.1 -64.9976 -64.9999
+ 0.15 -64.9966 -64.9996
+ 0.2 -64.9957 -64.9994
+ 0.25 -64.9949 -64.999
+ 0.3 -64.9941 -64.9987
+ 0.35 -64.9934 -64.9983
+ 0.4 -64.9927 -64.998
+ 0.45 -64.992 -64.9976
+ 0.5 -64.9914 -64.9972
+ 0.55 -64.9908 -64.9969
+ 0.6 -64.9902 -64.9965
+ 0.65 -64.9897 -64.9961
+ 0.7 -64.9891 -64.9958
+ 0.75 -64.9886 -64.9954
+ 0.8 -64.9882 -64.9951
+ 0.85 -64.9877 -64.9947
+ 0.9 -64.9873 -64.9944
+ 0.95 -64.9868 -64.9941
+ 1 -64.9864 -64.9937
+ 1.05 -64.986 -64.9934
+ 1.1 -64.9857 -64.9931
+ 1.15 -64.9853 -64.9928
+ 1.2 -64.985 -64.9926
+ 1.25 -64.9846 -64.9923
+ 1.3 -64.9843 -64.992
+ 1.35 -64.984 -64.9918
+ 1.4 -64.9837 -64.9915
+ 1.45 -64.9834 -64.9913
+ 1.5 -64.9832 -64.9911
+ 1.55 -64.9829 -64.9908
+ 1.6 -64.9827 -64.9906
+ 1.65 -64.9824 -64.9904
+ 1.7 -64.9822 -64.9902
+ 1.75 -64.982 -64.99
+ 1.8 -64.9818 -64.9899
+ 1.85 -64.9816 -64.9897
+ 1.9 -64.9814 -64.9895
+ 1.95 -64.9812 -64.9894
+ 2 -64.981 -64.9892
+ 2.05 -64.9809 -64.9891
+ 2.1 -64.9807 -64.9889
+ 2.15 -64.9806 -64.9888
+ 2.2 -64.9805 -64.9886
+ 2.25 -64.9803 -64.9885
+ 2.3 -64.9802 -64.9884
+ 2.35 -64.9801 -64.9883
+ 2.4 -64.98 -64.9882
+ 2.45 -64.9799 -64.9881
+ 2.5 -64.9798 -64.988
+ 2.55 -64.9797 -64.9879
+ 2.6 -64.9796 -64.9878
+ 2.65 -64.9795 -64.9877
+ 2.7 -64.9795 -64.9876
+ 2.75 -64.9794 -64.9876
+ 2.8 -64.9793 -64.9875
+ 2.85 -64.9793 -64.9874
+ 2.9 -64.9792 -64.9874
+ 2.95 -64.9792 -64.9873
+ 3 -64.9791 -64.9873
+ 3.05 -64.9791 -64.9872
+ 3.1 -64.9791 -64.9872
+ 3.15 -64.979 -64.9871
+ 3.2 -64.979 -64.9871
+ 3.25 -64.979 -64.987
+ 3.3 -64.979 -64.987
+ 3.35 -64.9789 -64.987
+ 3.4 -64.9789 -64.987
+ 3.45 -64.9789 -64.9869
+ 3.5 -64.9789 -64.9869
+ 3.55 -64.9789 -64.9869
+ 3.6 -64.9789 -64.9869
+ 3.65 -64.9789 -64.9868
+ 3.7 -64.9789 -64.9868
+ 3.75 -64.9789 -64.9868
+ 3.8 -64.9789 -64.9868
+ 3.85 -64.979 -64.9868
+ 3.9 -64.979 -64.9868
+ 3.95 -64.979 -64.9868
+ 4 -64.979 -64.9868
+ 4.05 -64.979 -64.9868
+ 4.1 -64.9791 -64.9868
+ 4.15 -64.9791 -64.9868
+ 4.2 -64.9791 -64.9868
+ 4.25 -64.9791 -64.9868
+ 4.3 -64.9792 -64.9868
+ 4.35 -64.9792 -64.9868
+ 4.4 -64.9792 -64.9869
+ 4.45 -64.9793 -64.9869
+ 4.5 -64.9793 -64.9869
+ 4.55 -64.9793 -64.9869
+ 4.6 -64.9794 -64.9869
+ 4.65 -64.9794 -64.9869
+ 4.7 -64.9795 -64.987
+ 4.75 -64.9795 -64.987
+ 4.8 -64.9796 -64.987
+ 4.85 -64.9796 -64.987
+ 4.9 -64.9796 -64.987
+ 4.95 -64.9797 -64.9871
+ 5 -64.9797 -64.9871
+ 5.05 -64.9798 -64.9871
+ 5.1 -64.9798 -64.9871
+ 5.15 -64.9799 -64.9872
+ 5.2 -64.9799 -64.9872
+ 5.25 -64.98 -64.9872
+ 5.3 -64.98 -64.9872
+ 5.35 -64.9801 -64.9873
+ 5.4 -64.9801 -64.9873
+ 5.45 -64.9802 -64.9873
+ 5.5 -64.9802 -64.9874
+ 5.55 -64.9803 -64.9874
+ 5.6 -64.9803 -64.9874
+ 5.65 -64.9804 -64.9874
+ 5.7 -64.9804 -64.9875
+ 5.75 -64.9805 -64.9875
+ 5.8 -64.9805 -64.9875
+ 5.85 -64.9806 -64.9876
+ 5.9 -64.9806 -64.9876
+ 5.95 -64.9807 -64.9876
+ 6 -64.9808 -64.9877
+ 6.05 -64.9808 -64.9877
+ 6.1 -64.9809 -64.9877
+ 6.15 -64.9809 -64.9878
+ 6.2 -64.981 -64.9878
+ 6.25 -64.981 -64.9878
+ 6.3 -64.9811 -64.9879
+ 6.35 -64.9811 -64.9879
+ 6.4 -64.9812 -64.9879
+ 6.45 -64.9812 -64.988
+ 6.5 -64.9813 -64.988
+ 6.55 -64.9813 -64.988
+ 6.6 -64.9814 -64.9881
+ 6.65 -64.9814 -64.9881
+ 6.7 -64.9815 -64.9881
+ 6.75 -64.9815 -64.9881
+ 6.8 -64.9816 -64.9882
+ 6.85 -64.9816 -64.9882
+ 6.9 -64.9817 -64.9882
+ 6.95 -64.9817 -64.9883
+ 7 -64.9818 -64.9883
+ 7.05 -64.9818 -64.9883
+ 7.1 -64.9819 -64.9884
+ 7.15 -64.9819 -64.9884
+ 7.2 -64.9819 -64.9884
+ 7.25 -64.982 -64.9885
+ 7.3 -64.982 -64.9885
+ 7.35 -64.9821 -64.9885
+ 7.4 -64.9821 -64.9885
+ 7.45 -64.9822 -64.9886
+ 7.5 -64.9822 -64.9886
+ 7.55 -64.9823 -64.9886
+ 7.6 -64.9823 -64.9887
+ 7.65 -64.9823 -64.9887
+ 7.7 -64.9824 -64.9887
+ 7.75 -64.9824 -64.9887
+ 7.8 -64.9825 -64.9888
+ 7.85 -64.9825 -64.9888
+ 7.9 -64.9826 -64.9888
+ 7.95 -64.9826 -64.9889
+ 8 -64.9826 -64.9889
+ 8.05 -64.9827 -64.9889
+ 8.1 -64.9827 -64.9889
+ 8.15 -64.9827 -64.989
+ 8.2 -64.9828 -64.989
+ 8.25 -64.9828 -64.989
+ 8.3 -64.9829 -64.989
+ 8.35 -64.9829 -64.9891
+ 8.4 -64.9829 -64.9891
+ 8.45 -64.983 -64.9891
+ 8.5 -64.983 -64.9891
+ 8.55 -64.983 -64.9892
+ 8.6 -64.9831 -64.9892
+ 8.65 -64.9831 -64.9892
+ 8.7 -64.9831 -64.9892
+ 8.75 -64.9832 -64.9893
+ 8.8 -64.9832 -64.9893
+ 8.85 -64.9832 -64.9893
+ 8.9 -64.9833 -64.9893
+ 8.95 -64.9833 -64.9893
+ 9 -64.9833 -64.9894
+ 9.05 -64.9834 -64.9894
+ 9.1 -64.9834 -64.9894
+ 9.15 -64.9834 -64.9894
+ 9.2 -64.9834 -64.9894
+ 9.25 -64.9835 -64.9895
+ 9.3 -64.9835 -64.9895
+ 9.35 -64.9835 -64.9895
+ 9.4 -64.9835 -64.9895
+ 9.45 -64.9836 -64.9895
+ 9.5 -64.9836 -64.9896
+ 9.55 -64.9836 -64.9896
+ 9.6 -64.9837 -64.9896
+ 9.65 -64.9837 -64.9896
+ 9.7 -64.9837 -64.9896
+ 9.75 -64.9837 -64.9896
+ 9.8 -64.9837 -64.9897
+ 9.85 -64.9838 -64.9897
+ 9.9 -64.9838 -64.9897
+ 9.95 -64.9838 -64.9897
+ 10 -64.9838 -64.9897
+ 10.05 -64.7901 -30.2938
+ 10.1 -63.5961 -24.971
+ 10.15 -61.7946 -22.2475
+ 10.2 -59.8162 -20.5105
+ 10.25 -57.8346 -19.2806
+ 10.3 -55.9105 -18.3638
+ 10.35 -54.0572 -17.6636
+ 10.4 -52.2655 -17.1227
+ 10.45 -50.5111 -16.7017
+ 10.5 -48.7559 -16.3712
+ 10.55 -46.9476 -16.107
+ 10.6 -45.0152 -15.8878
+ 10.65 -42.8632 -15.6933
+ 10.7 -40.3611 -15.5024
+ 10.75 -37.3264 -15.2912
+ 10.8 -33.4992 -15.0302
+ 10.85 -28.5098 -14.6807
+ 10.9 -21.8567 -14.1898
+ 10.95 -12.9706 -13.4856
+ 11 -1.57469 -12.4785
+ 11.05 11.4421 -11.0884
+ 11.1 23.4366 -9.32351
+ 11.15 31.6923 -7.37006
+ 11.2 35.8502 -5.5413
+ 11.25 37.2894 -4.07525
+ 11.3 37.2952 -3.02811
+ 11.35 36.5481 -2.33932
+ 11.4 35.3514 -1.92365
+ 11.45 33.8453 -1.71228
+ 11.5 32.1057 -1.658
+ 11.55 30.1812 -1.72945
+ 11.6 28.1082 -1.90517
+ 11.65 25.9159 -2.16959
+ 11.7 23.6296 -2.51074
+ 11.75 21.2709 -2.91889
+ 11.8 18.8589 -3.38589
+ 11.85 16.4103 -3.90472
+ 11.9 13.9394 -4.46923
+ 11.95 11.4585 -5.07404
+ 12 8.97788 -5.71437
+ 12.05 6.50609 -6.38602
+ 12.1 4.05009 -7.08527
+ 12.15 1.61541 -7.80886
+ 12.2 -0.793676 -8.55394
+ 12.25 -3.17399 -9.31802
+ 12.3 -5.52329 -10.099
+ 12.35 -7.84019 -10.8949
+ 12.4 -10.124 -11.7042
+ 12.45 -12.3747 -12.5256
+ 12.5 -14.5928 -13.358
+ 12.55 -16.7795 -14.2004
+ 12.6 -18.9365 -15.0522
+ 12.65 -21.0664 -15.9129
+ 12.7 -23.1726 -16.7822
+ 12.75 -25.2592 -17.66
+ 12.8 -27.3319 -18.5465
+ 12.85 -29.3975 -19.4423
+ 12.9 -31.4644 -20.3483
+ 12.95 -33.5428 -21.2656
+ 13 -35.6445 -22.1961
+ 13.05 -37.7828 -23.142
+ 13.1 -39.9724 -24.1061
+ 13.15 -42.2282 -25.0919
+ 13.2 -44.5635 -26.1032
+ 13.25 -46.9883 -27.1439
+ 13.3 -49.5042 -28.2183
+ 13.35 -52.1002 -29.3294
+ 13.4 -54.7462 -30.4787
+ 13.45 -57.3884 -31.665
+ 13.5 -59.9493 -32.8823
+ 13.55 -62.3364 -34.1192
+ 13.6 -64.4601 -35.3591
+ 13.65 -66.2559 -36.5808
+ 13.7 -67.6994 -37.7626
+ 13.75 -68.8077 -38.8855
+ 13.8 -69.6265 -39.9364
+ 13.85 -70.2137 -40.9089
+ 13.9 -70.6256 -41.8032
+ 13.95 -70.9098 -42.6236
+ 14 -71.1029 -43.3767
+ 14.05 -71.2319 -44.0704
+ 14.1 -71.3158 -44.7121
+ 14.15 -71.3679 -45.3088
+ 14.2 -71.3974 -45.8666
+ 14.25 -71.4104 -46.3908
+ 14.3 -71.4115 -46.8858
+ 14.35 -71.4037 -47.3555
+ 14.4 -71.3892 -47.8029
+ 14.45 -71.3695 -48.2309
+ 14.5 -71.3458 -48.6415
+ 14.55 -71.319 -49.0369
+ 14.6 -71.2896 -49.4185
+ 14.65 -71.2582 -49.7879
+ 14.7 -71.2253 -50.1461
+ 14.75 -71.191 -50.4941
+ 14.8 -71.1557 -50.8329
+ 14.85 -71.1195 -51.1631
+ 14.9 -71.0827 -51.4855
+ 14.95 -71.0453 -51.8005
+ 15 -71.0076 -52.1085
+ 15.05 -70.9695 -52.4101
+ 15.1 -70.9313 -52.7056
+ 15.15 -70.893 -52.9951
+ 15.2 -70.8545 -53.2791
+ 15.25 -70.8161 -53.5578
+ 15.3 -70.7778 -53.8312
+ 15.35 -70.7395 -54.0997
+ 15.4 -70.7014 -54.3633
+ 15.45 -70.6634 -54.6222
+ 15.5 -70.6256 -54.8765
+ 15.55 -70.588 -55.1264
+ 15.6 -70.5507 -55.3719
+ 15.65 -70.5136 -55.613
+ 15.7 -70.4767 -55.8499
+ 15.75 -70.4402 -56.0827
+ 15.8 -70.4039 -56.3114
+ 15.85 -70.3679 -56.5361
+ 15.9 -70.3323 -56.7568
+ 15.95 -70.2969 -56.9737
+ 16 -70.2619 -57.1866
+ 16.05 -70.2271 -57.3957
+ 16.1 -70.1927 -57.6011
+ 16.15 -70.1586 -57.8028
+ 16.2 -70.1248 -58.0008
+ 16.25 -70.0914 -58.1951
+ 16.3 -70.0582 -58.3859
+ 16.35 -70.0254 -58.5731
+ 16.4 -69.9928 -58.7568
+ 16.45 -69.9606 -58.9371
+ 16.5 -69.9287 -59.114
+ 16.55 -69.8971 -59.2874
+ 16.6 -69.8657 -59.4575
+ 16.65 -69.8347 -59.6243
+ 16.7 -69.8039 -59.7879
+ 16.75 -69.7734 -59.9482
+ 16.8 -69.7432 -60.1054
+ 16.85 -69.7133 -60.2594
+ 16.9 -69.6836 -60.4103
+ 16.95 -69.6541 -60.5581
+ 17 -69.625 -60.7029
+ 17.05 -69.596 -60.8448
+ 17.1 -69.5673 -60.9837
+ 17.15 -69.5389 -61.1197
+ 17.2 -69.5106 -61.2528
+ 17.25 -69.4826 -61.3831
+ 17.3 -69.4548 -61.5106
+ 17.35 -69.4272 -61.6354
+ 17.4 -69.3998 -61.7574
+ 17.45 -69.3726 -61.8768
+ 17.5 -69.3455 -61.9936
+ 17.55 -69.3187 -62.1078
+ 17.6 -69.2921 -62.2194
+ 17.65 -69.2656 -62.3286
+ 17.7 -69.2393 -62.4352
+ 17.75 -69.2131 -62.5395
+ 17.8 -69.1871 -62.6413
+ 17.85 -69.1612 -62.7408
+ 17.9 -69.1355 -62.838
+ 17.95 -69.11 -62.9328
+ 18 -69.0845 -63.0255
+ 18.05 -69.0592 -63.1159
+ 18.1 -69.034 -63.2042
+ 18.15 -69.009 -63.2903
+ 18.2 -68.984 -63.3744
+ 18.25 -68.9592 -63.4564
+ 18.3 -68.9345 -63.5363
+ 18.35 -68.9099 -63.6143
+ 18.4 -68.8854 -63.6903
+ 18.45 -68.8609 -63.7644
+ 18.5 -68.8366 -63.8367
+ 18.55 -68.8124 -63.907
+ 18.6 -68.7882 -63.9756
+ 18.65 -68.7641 -64.0424
+ 18.7 -68.7402 -64.1074
+ 18.75 -68.7162 -64.1707
+ 18.8 -68.6924 -64.2323
+ 18.85 -68.6686 -64.2923
+ 18.9 -68.6449 -64.3507
+ 18.95 -68.6213 -64.4074
+ 19 -68.5977 -64.4626
+ 19.05 -68.5742 -64.5163
+ 19.1 -68.5508 -64.5685
+ 19.15 -68.5274 -64.6192
+ 19.2 -68.504 -64.6685
+ 19.25 -68.4807 -64.7164
+ 19.3 -68.4575 -64.7628
+ 19.35 -68.4343 -64.808
+ 19.4 -68.4112 -64.8518
+ 19.45 -68.3881 -64.8943
+ 19.5 -68.3651 -64.9355
+ 19.55 -68.3421 -64.9755
+ 19.6 -68.3192 -65.0142
+ 19.65 -68.2963 -65.0518
+ 19.7 -68.2735 -65.0882
+ 19.75 -68.2507 -65.1234
+ 19.8 -68.2279 -65.1575
+ 19.85 -68.2052 -65.1906
+ 19.9 -68.1825 -65.2225
+ 19.95 -68.1599 -65.2534
+ 20 -68.1373 -65.2833
+ 20.05 -67.9216 -30.561
+ 20.1 -66.7223 -25.2682
+ 20.15 -64.9549 -22.568
+ 20.2 -63.0569 -20.853
+ 20.25 -61.2029 -19.6469
+ 20.3 -59.4542 -18.758
+ 20.35 -57.8284 -18.0909
+ 20.4 -56.325 -17.5895
+ 20.45 -54.9346 -17.2162
+ 20.5 -53.6429 -16.9436
+ 20.55 -52.4318 -16.7506
+ 20.6 -51.2801 -16.6205
+ 20.65 -50.1639 -16.5396
+ 20.7 -49.0567 -16.4961
+ 20.75 -47.9286 -16.4799
+ 20.8 -46.7461 -16.4815
+ 20.85 -45.4705 -16.492
+ 20.9 -44.0554 -16.5022
+ 20.95 -42.4434 -16.5024
+ 21 -40.5617 -16.4817
+ 21.05 -38.3147 -16.427
+ 21.1 -35.5758 -16.3222
+ 21.15 -32.1776 -16.1465
+ 21.2 -27.9055 -15.8728
+ 21.25 -22.5081 -15.466
+ 21.3 -15.7565 -14.8824
+ 21.35 -7.60396 -14.0738
+ 21.4 1.53283 -13.0023
+ 21.45 10.5759 -11.6711
+ 21.5 18.1225 -10.1598
+ 21.55 23.2744 -8.62532
+ 21.6 26.0558 -7.24189
+ 21.65 27.0384 -6.12681
+ 21.7 26.8268 -5.31593
+ 21.75 25.85 -4.78853
+ 21.8 24.3724 -4.50201
+ 21.85 22.5543 -4.41309
+ 21.9 20.4973 -4.48556
+ 21.95 18.2708 -4.69092
+ 22 15.9247 -5.00698
+ 22.05 13.4968 -5.41609
+ 22.1 11.0165 -5.90396
+ 22.15 8.50691 -6.45874
+ 22.2 5.98607 -7.07046
+ 22.25 3.46822 -7.73069
+ 22.3 0.964352 -8.43221
+ 22.35 -1.51719 -9.1689
+ 22.4 -3.97021 -9.93552
+ 22.45 -6.39031 -10.7276
+ 22.5 -8.77454 -11.5415
+ 22.55 -11.1212 -12.3739
+ 22.6 -13.4295 -13.2222
+ 22.65 -15.6998 -14.0842
+ 22.7 -17.933 -14.9582
+ 22.75 -20.131 -15.8427
+ 22.8 -22.2966 -16.7366
+ 22.85 -24.4335 -17.6393
+ 22.9 -26.5467 -18.5503
+ 22.95 -28.6423 -19.4695
+ 23 -30.7281 -20.3975
+ 23.05 -32.8132 -21.3349
+ 23.1 -34.9085 -22.283
+ 23.15 -37.0262 -23.2438
+ 23.2 -39.1797 -24.2194
+ 23.25 -41.3829 -25.2128
+ 23.3 -43.6488 -26.2273
+ 23.35 -45.9876 -27.2666
+ 23.4 -48.404 -28.3344
+ 23.45 -50.8927 -29.4339
+ 23.5 -53.4335 -30.5671
+ 23.55 -55.9875 -31.7338
+ 23.6 -58.4947 -32.9303
+ 23.65 -60.8783 -34.1486
+ 23.7 -63.0563 -35.3755
+ 23.75 -64.959 -36.5937
+ 23.8 -66.5454 -37.7834
+ 23.85 -67.8104 -38.9257
+ 23.9 -68.7805 -40.0056
+ 23.95 -69.5014 -41.0137
+ 24 -70.0245 -41.9465
+ 24.05 -70.3974 -42.8055
+ 24.1 -70.6598 -43.5954
+ 24.15 -70.8421 -44.3228
+ 24.2 -70.9672 -44.9945
+ 24.25 -71.051 -45.6175
+ 24.3 -71.1052 -46.1979
+ 24.35 -71.1377 -46.7414
+ 24.4 -71.1544 -47.2526
+ 24.45 -71.1593 -47.7357
+ 24.5 -71.1552 -48.1942
+ 24.55 -71.1442 -48.631
+ 24.6 -71.128 -49.0487
+ 24.65 -71.1075 -49.4494
+ 24.7 -71.0837 -49.835
+ 24.75 -71.0572 -50.2069
+ 24.8 -71.0286 -50.5666
+ 24.85 -70.9983 -50.9152
+ 24.9 -70.9665 -51.2536
+ 24.95 -70.9336 -51.5828
+ 25 -70.8998 -51.9033
+ 25.05 -70.8652 -52.216
+ 25.1 -70.8301 -52.5212
+ 25.15 -70.7945 -52.8194
+ 25.2 -70.7585 -53.1112
+ 25.25 -70.7224 -53.3967
+ 25.3 -70.686 -53.6764
+ 25.35 -70.6496 -53.9504
+ 25.4 -70.6131 -54.219
+ 25.45 -70.5766 -54.4825
+ 25.5 -70.5402 -54.7409
+ 25.55 -70.504 -54.9946
+ 25.6 -70.4678 -55.2435
+ 25.65 -70.4318 -55.4878
+ 25.7 -70.396 -55.7276
+ 25.75 -70.3605 -55.9631
+ 25.8 -70.3251 -56.1943
+ 25.85 -70.29 -56.4214
+ 25.9 -70.2552 -56.6443
+ 25.95 -70.2206 -56.8631
+ 26 -70.1863 -57.078
+ 26.05 -70.1523 -57.289
+ 26.1 -70.1186 -57.4962
+ 26.15 -70.0852 -57.6995
+ 26.2 -70.0521 -57.8991
+ 26.25 -70.0193 -58.095
+ 26.3 -69.9867 -58.2873
+ 26.35 -69.9545 -58.476
+ 26.4 -69.9226 -58.6611
+ 26.45 -69.891 -58.8427
+ 26.5 -69.8597 -59.0209
+ 26.55 -69.8286 -59.1957
+ 26.6 -69.7979 -59.367
+ 26.65 -69.7674 -59.5351
+ 26.7 -69.7373 -59.6999
+ 26.75 -69.7074 -59.8614
+ 26.8 -69.6777 -60.0197
+ 26.85 -69.6484 -60.1749
+ 26.9 -69.6192 -60.3269
+ 26.95 -69.5904 -60.4759
+ 27 -69.5618 -60.6218
+ 27.05 -69.5334 -60.7647
+ 27.1 -69.5053 -60.9047
+ 27.15 -69.4774 -61.0418
+ 27.2 -69.4498 -61.1759
+ 27.25 -69.4223 -61.3073
+ 27.3 -69.3951 -61.4358
+ 27.35 -69.3681 -61.5616
+ 27.4 -69.3413 -61.6847
+ 27.45 -69.3146 -61.8051
+ 27.5 -69.2882 -61.9229
+ 27.55 -69.262 -62.0381
+ 27.6 -69.2359 -62.1507
+ 27.65 -69.21 -62.2608
+ 27.7 -69.1843 -62.3684
+ 27.75 -69.1587 -62.4736
+ 27.8 -69.1333 -62.5763
+ 27.85 -69.108 -62.6767
+ 27.9 -69.0829 -62.7748
+ 27.95 -69.0579 -62.8706
+ 28 -69.0331 -62.9642
+ 28.05 -69.0084 -63.0555
+ 28.1 -68.9838 -63.1447
+ 28.15 -68.9593 -63.2317
+ 28.2 -68.9349 -63.3166
+ 28.25 -68.9107 -63.3994
+ 28.3 -68.8865 -63.4803
+ 28.35 -68.8625 -63.5591
+ 28.4 -68.8386 -63.6359
+ 28.45 -68.8147 -63.7109
+ 28.5 -68.791 -63.7839
+ 28.55 -68.7673 -63.8551
+ 28.6 -68.7437 -63.9245
+ 28.65 -68.7202 -63.992
+ 28.7 -68.6968 -64.0579
+ 28.75 -68.6735 -64.1219
+ 28.8 -68.6502 -64.1843
+ 28.85 -68.627 -64.2451
+ 28.9 -68.6038 -64.3042
+ 28.95 -68.5807 -64.3617
+ 29 -68.5577 -64.4176
+ 29.05 -68.5348 -64.472
+ 29.1 -68.5119 -64.5249
+ 29.15 -68.489 -64.5764
+ 29.2 -68.4662 -64.6263
+ 29.25 -68.4435 -64.6749
+ 29.3 -68.4208 -64.7221
+ 29.35 -68.3981 -64.7679
+ 29.4 -68.3755 -64.8123
+ 29.45 -68.353 -64.8555
+ 29.5 -68.3305 -64.8974
+ 29.55 -68.308 -64.938
+ 29.6 -68.2856 -64.9774
+ 29.65 -68.2632 -65.0156
+ 29.7 -68.2408 -65.0526
+ 29.75 -68.2185 -65.0884
+ 29.8 -68.1963 -65.1232
+ 29.85 -68.1741 -65.1568
+ 29.9 -68.1519 -65.1893
+ 29.95 -68.1297 -65.2208
+ 30 -68.1076 -65.2513
+ 30.05 -68.0856 -65.2807
+ 30.1 -68.0635 -65.3092
+ 30.15 -68.0415 -65.3366
+ 30.2 -68.0196 -65.3632
+ 30.25 -67.9977 -65.3888
+ 30.3 -67.9758 -65.4135
+ 30.35 -67.954 -65.4373
+ 30.4 -67.9322 -65.4603
+ 30.45 -67.9104 -65.4824
+ 30.5 -67.8887 -65.5037
+ 30.55 -67.867 -65.5242
+ 30.6 -67.8454 -65.5439
+ 30.65 -67.8238 -65.5628
+ 30.7 -67.8023 -65.581
+ 30.75 -67.7808 -65.5985
+ 30.8 -67.7593 -65.6152
+ 30.85 -67.7379 -65.6312
+ 30.9 -67.7165 -65.6466
+ 30.95 -67.6952 -65.6613
+ 31 -67.6739 -65.6753
+ 31.05 -67.6527 -65.6887
+ 31.1 -67.6315 -65.7015
+ 31.15 -67.6104 -65.7137
+ 31.2 -67.5893 -65.7253
+ 31.25 -67.5683 -65.7363
+ 31.3 -67.5474 -65.7468
+ 31.35 -67.5264 -65.7568
+ 31.4 -67.5056 -65.7662
+ 31.45 -67.4848 -65.775
+ 31.5 -67.464 -65.7834
+ 31.55 -67.4434 -65.7913
+ 31.6 -67.4227 -65.7987
+ 31.65 -67.4022 -65.8057
+ 31.7 -67.3817 -65.8122
+ 31.75 -67.3612 -65.8182
+ 31.8 -67.3409 -65.8239
+ 31.85 -67.3206 -65.8291
+ 31.9 -67.3003 -65.8339
+ 31.95 -67.2802 -65.8383
+ 32 -67.2601 -65.8423
+ 32.05 -67.24 -65.846
+ 32.1 -67.2201 -65.8493
+ 32.15 -67.2002 -65.8522
+ 32.2 -67.1804 -65.8548
+ 32.25 -67.1606 -65.8571
+ 32.3 -67.141 -65.859
+ 32.35 -67.1214 -65.8607
+ 32.4 -67.1019 -65.862
+ 32.45 -67.0824 -65.863
+ 32.5 -67.0631 -65.8638
+ 32.55 -67.0438 -65.8643
+ 32.6 -67.0246 -65.8645
+ 32.65 -67.0055 -65.8644
+ 32.7 -66.9865 -65.8641
+ 32.75 -66.9676 -65.8635
+ 32.8 -66.9487 -65.8627
+ 32.85 -66.9299 -65.8617
+ 32.9 -66.9113 -65.8605
+ 32.95 -66.8927 -65.859
+ 33 -66.8742 -65.8573
+ 33.05 -66.8558 -65.8555
+ 33.1 -66.8374 -65.8534
+ 33.15 -66.8192 -65.8512
+ 33.2 -66.8011 -65.8487
+ 33.25 -66.783 -65.8461
+ 33.3 -66.7651 -65.8433
+ 33.35 -66.7472 -65.8404
+ 33.4 -66.7295 -65.8373
+ 33.45 -66.7118 -65.834
+ 33.5 -66.6943 -65.8306
+ 33.55 -66.6768 -65.8271
+ 33.6 -66.6594 -65.8234
+ 33.65 -66.6422 -65.8196
+ 33.7 -66.625 -65.8157
+ 33.75 -66.6079 -65.8116
+ 33.8 -66.591 -65.8075
+ 33.85 -66.5741 -65.8032
+ 33.9 -66.5573 -65.7988
+ 33.95 -66.5407 -65.7943
+ 34 -66.5241 -65.7897
+ 34.05 -66.5076 -65.7851
+ 34.1 -66.4913 -65.7803
+ 34.15 -66.475 -65.7755
+ 34.2 -66.4589 -65.7705
+ 34.25 -66.4429 -65.7655
+ 34.3 -66.4269 -65.7605
+ 34.35 -66.4111 -65.7553
+ 34.4 -66.3954 -65.7501
+ 34.45 -66.3798 -65.7448
+ 34.5 -66.3643 -65.7395
+ 34.55 -66.3489 -65.7341
+ 34.6 -66.3336 -65.7287
+ 34.65 -66.3184 -65.7232
+ 34.7 -66.3033 -65.7176
+ 34.75 -66.2883 -65.7121
+ 34.8 -66.2735 -65.7065
+ 34.85 -66.2587 -65.7008
+ 34.9 -66.2441 -65.6951
+ 34.95 -66.2295 -65.6894
+ 35 -66.2151 -65.6836
+ 35.05 -66.2008 -65.6779
+ 35.1 -66.1866 -65.6721
+ 35.15 -66.1725 -65.6662
+ 35.2 -66.1585 -65.6604
+ 35.25 -66.1446 -65.6545
+ 35.3 -66.1308 -65.6486
+ 35.35 -66.1172 -65.6428
+ 35.4 -66.1036 -65.6368
+ 35.45 -66.0902 -65.6309
+ 35.5 -66.0768 -65.625
+ 35.55 -66.0636 -65.6191
+ 35.6 -66.0505 -65.6132
+ 35.65 -66.0375 -65.6072
+ 35.7 -66.0246 -65.6013
+ 35.75 -66.0118 -65.5954
+ 35.8 -65.9991 -65.5895
+ 35.85 -65.9865 -65.5836
+ 35.9 -65.9741 -65.5776
+ 35.95 -65.9617 -65.5717
+ 36 -65.9494 -65.5658
+ 36.05 -65.9373 -65.56
+ 36.1 -65.9253 -65.5541
+ 36.15 -65.9133 -65.5482
+ 36.2 -65.9015 -65.5424
+ 36.25 -65.8898 -65.5366
+ 36.3 -65.8782 -65.5308
+ 36.35 -65.8667 -65.525
+ 36.4 -65.8553 -65.5192
+ 36.45 -65.844 -65.5135
+ 36.5 -65.8328 -65.5077
+ 36.55 -65.8217 -65.502
+ 36.6 -65.8107 -65.4963
+ 36.65 -65.7999 -65.4907
+ 36.7 -65.7891 -65.485
+ 36.75 -65.7784 -65.4794
+ 36.8 -65.7679 -65.4739
+ 36.85 -65.7574 -65.4683
+ 36.9 -65.747 -65.4628
+ 36.95 -65.7368 -65.4573
+ 37 -65.7266 -65.4518
+ 37.05 -65.7166 -65.4464
+ 37.1 -65.7066 -65.441
+ 37.15 -65.6968 -65.4356
+ 37.2 -65.687 -65.4303
+ 37.25 -65.6774 -65.425
+ 37.3 -65.6678 -65.4197
+ 37.35 -65.6584 -65.4145
+ 37.4 -65.649 -65.4093
+ 37.45 -65.6398 -65.4041
+ 37.5 -65.6306 -65.399
+ 37.55 -65.6215 -65.3939
+ 37.6 -65.6126 -65.3888
+ 37.65 -65.6037 -65.3838
+ 37.7 -65.5949 -65.3788
+ 37.75 -65.5862 -65.3738
+ 37.8 -65.5776 -65.3689
+ 37.85 -65.5692 -65.3641
+ 37.9 -65.5608 -65.3592
+ 37.95 -65.5524 -65.3544
+ 38 -65.5442 -65.3497
+ 38.05 -65.5361 -65.3449
+ 38.1 -65.5281 -65.3403
+ 38.15 -65.5201 -65.3356
+ 38.2 -65.5123 -65.331
+ 38.25 -65.5045 -65.3265
+ 38.3 -65.4968 -65.3219
+ 38.35 -65.4892 -65.3175
+ 38.4 -65.4817 -65.313
+ 38.45 -65.4743 -65.3086
+ 38.5 -65.467 -65.3042
+ 38.55 -65.4598 -65.2999
+ 38.6 -65.4526 -65.2956
+ 38.65 -65.4455 -65.2914
+ 38.7 -65.4385 -65.2872
+ 38.75 -65.4316 -65.283
+ 38.8 -65.4248 -65.2789
+ 38.85 -65.4181 -65.2748
+ 38.9 -65.4114 -65.2708
+ 38.95 -65.4049 -65.2668
+ 39 -65.3984 -65.2628
+ 39.05 -65.3919 -65.2589
+ 39.1 -65.3856 -65.255
+ 39.15 -65.3794 -65.2512
+ 39.2 -65.3732 -65.2474
+ 39.25 -65.3671 -65.2436
+ 39.3 -65.3611 -65.2399
+ 39.35 -65.3551 -65.2362
+ 39.4 -65.3492 -65.2326
+ 39.45 -65.3434 -65.229
+ 39.5 -65.3377 -65.2254
+ 39.55 -65.3321 -65.2219
+ 39.6 -65.3265 -65.2184
+ 39.65 -65.321 -65.215
+ 39.7 -65.3156 -65.2116
+ 39.75 -65.3102 -65.2082
+ 39.8 -65.3049 -65.2049
+ 39.85 -65.2997 -65.2016
+ 39.9 -65.2945 -65.1984
+ 39.95 -65.2895 -65.1951
+ 40 -65.2844 -65.192
+ 40.05 -65.0852 -30.387
+ 40.1 -63.883 -25.0466
+ 40.15 -62.0728 -22.3136
+ 40.2 -60.0871 -20.5704
+ 40.25 -58.1011 -19.336
+ 40.3 -56.1765 -18.4161
+ 40.35 -54.3282 -17.7138
+ 40.4 -52.5488 -17.1719
+ 40.45 -50.8165 -16.7512
+ 40.5 -49.0972 -16.4224
+ 40.55 -47.3436 -16.1618
+ 40.6 -45.4924 -15.9486
+ 40.65 -43.4595 -15.7636
+ 40.7 -41.1317 -15.5871
+ 40.75 -38.353 -15.3973
+ 40.8 -34.905 -15.168
+ 40.85 -30.4782 -14.8658
+ 40.9 -24.6469 -14.4456
+ 40.95 -16.8869 -13.8458
+ 41 -6.7721 -12.9871
+ 41.05 5.41828 -11.7849
+ 41.1 17.8741 -10.2018
+ 41.15 27.7647 -8.33735
+ 41.2 33.6114 -6.45343
+ 41.25 36.137 -4.83166
+ 41.3 36.7126 -3.60883
+ 41.35 36.2579 -2.77157
+ 41.4 35.221 -2.24412
+ 41.45 33.8085 -1.95123
+ 41.5 32.125 -1.83711
+ 41.55 30.2328 -1.86387
+ 41.6 28.1755 -2.00575
+ 41.65 25.9871 -2.24442
+ 41.7 23.6956 -2.5661
+ 41.75 21.3252 -2.9598
+ 41.8 18.8967 -3.41641
+ 41.85 16.4283 -3.92818
+ 41.9 13.9354 -4.48836
+ 41.95 11.4313 -5.09106
+ 42 8.92706 -5.73107
+ 42.05 6.43183 -6.40384
+ 42.1 3.95303 -7.10534
+ 42.15 1.49652 -7.83205
+ 42.2 -0.93318 -8.58089
+ 42.25 -3.33273 -9.34919
+ 42.3 -5.6998 -10.1346
+ 42.35 -8.03293 -10.9353
+ 42.4 -10.3314 -11.7494
+ 42.45 -12.5953 -12.5756
+ 42.5 -14.8251 -13.4126
+ 42.55 -17.0221 -14.2597
+ 42.6 -19.1881 -15.1159
+ 42.65 -21.3258 -15.9807
+ 42.7 -23.4386 -16.854
+ 42.75 -25.531 -17.7356
+ 42.8 -27.6087 -18.6257
+ 42.85 -29.6787 -19.5248
+ 42.9 -31.7497 -20.4339
+ 42.95 -33.8318 -21.3542
+ 43 -35.9372 -22.2875
+ 43.05 -38.0792 -23.2361
+ 43.1 -40.2724 -24.203
+ 43.15 -42.5312 -25.1914
+ 43.2 -44.8685 -26.2052
+ 43.25 -47.293 -27.2483
+ 43.3 -49.8051 -28.3247
+ 43.35 -52.3917 -29.4373
+ 43.4 -55.0204 -30.5872
+ 43.45 -57.6358 -31.7726
+ 43.5 -60.1599 -32.9871
+ 43.55 -62.502 -34.2189
+ 43.6 -64.5765 -35.451
+ 43.65 -66.3242 -36.6627
+ 43.7 -67.7257 -37.8327
+ 43.75 -68.8005 -38.9431
+ 43.8 -69.595 -39.9815
+ 43.85 -70.1656 -40.9425
+ 43.9 -70.5669 -41.8265
+ 43.95 -70.8445 -42.6379
+ 44 -71.0337 -43.3836
+ 44.05 -71.1606 -44.0709
+ 44.1 -71.2433 -44.7073
+ 44.15 -71.2948 -45.2996
+ 44.2 -71.324 -45.8537
+ 44.25 -71.337 -46.3748
+ 44.3 -71.3381 -46.8671
+ 44.35 -71.3304 -47.3345
+ 44.4 -71.3159 -47.78
+ 44.45 -71.2964 -48.2063
+ 44.5 -71.2728 -48.6155
+ 44.55 -71.246 -49.0096
+ 44.6 -71.2168 -49.3902
+ 44.65 -71.1855 -49.7585
+ 44.7 -71.1526 -50.1159
+ 44.75 -71.1185 -50.4632
+ 44.8 -71.0833 -50.8013
+ 44.85 -71.0472 -51.131
+ 44.9 -71.0105 -51.4529
+ 44.95 -70.9733 -51.7674
+ 45 -70.9358 -52.0751
+ 45.05 -70.8979 -52.3763
+ 45.1 -70.8598 -52.6714
+ 45.15 -70.8217 -52.9607
+ 45.2 -70.7834 -53.2444
+ 45.25 -70.7452 -53.5228
+ 45.3 -70.7071 -53.7961
+ 45.35 -70.6691 -54.0644
+ 45.4 -70.6311 -54.3278
+ 45.45 -70.5934 -54.5866
+ 45.5 -70.5559 -54.8408
+ 45.55 -70.5186 -55.0905
+ 45.6 -70.4815 -55.3359
+ 45.65 -70.4447 -55.577
+ 45.7 -70.4082 -55.8138
+ 45.75 -70.3719 -56.0466
+ 45.8 -70.336 -56.2752
+ 45.85 -70.3003 -56.4999
+ 45.9 -70.265 -56.7206
+ 45.95 -70.23 -56.9374
+ 46 -70.1953 -57.1503
+ 46.05 -70.1609 -57.3595
+ 46.1 -70.1269 -57.5649
+ 46.15 -70.0931 -57.7666
+ 46.2 -70.0597 -57.9646
+ 46.25 -70.0267 -58.159
+ 46.3 -69.9939 -58.3498
+ 46.35 -69.9615 -58.5371
+ 46.4 -69.9293 -58.7209
+ 46.45 -69.8975 -58.9013
+ 46.5 -69.866 -59.0782
+ 46.55 -69.8348 -59.2517
+ 46.6 -69.8039 -59.422
+ 46.65 -69.7733 -59.5889
+ 46.7 -69.743 -59.7525
+ 46.75 -69.7129 -59.913
+ 46.8 -69.6832 -60.0703
+ 46.85 -69.6537 -60.2244
+ 46.9 -69.6245 -60.3755
+ 46.95 -69.5955 -60.5235
+ 47 -69.5668 -60.6684
+ 47.05 -69.5383 -60.8104
+ 47.1 -69.5101 -60.9495
+ 47.15 -69.4821 -61.0856
+ 47.2 -69.4543 -61.2189
+ 47.25 -69.4268 -61.3494
+ 47.3 -69.3994 -61.4771
+ 47.35 -69.3723 -61.6021
+ 47.4 -69.3454 -61.7243
+ 47.45 -69.3187 -61.8439
+ 47.5 -69.2921 -61.9609
+ 47.55 -69.2658 -62.0753
+ 47.6 -69.2396 -62.1871
+ 47.65 -69.2136 -62.2965
+ 47.7 -69.1878 -62.4034
+ 47.75 -69.1621 -62.5078
+ 47.8 -69.1366 -62.6099
+ 47.85 -69.1113 -62.7096
+ 47.9 -69.086 -62.807
+ 47.95 -69.061 -62.9021
+ 48 -69.036 -62.995
+ 48.05 -69.0112 -63.0856
+ 48.1 -68.9865 -63.1741
+ 48.15 -68.9619 -63.2605
+ 48.2 -68.9375 -63.3448
+ 48.25 -68.9131 -63.427
+ 48.3 -68.8889 -63.5072
+ 48.35 -68.8647 -63.5855
+ 48.4 -68.8407 -63.6617
+ 48.45 -68.8168 -63.7361
+ 48.5 -68.7929 -63.8085
+ 48.55 -68.7691 -63.8792
+ 48.6 -68.7454 -63.948
+ 48.65 -68.7218 -64.015
+ 48.7 -68.6983 -64.0803
+ 48.75 -68.6749 -64.1439
+ 48.8 -68.6515 -64.2057
+ 48.85 -68.6282 -64.266
+ 48.9 -68.6049 -64.3246
+ 48.95 -68.5817 -64.3816
+ 49 -68.5586 -64.4371
+ 49.05 -68.5356 -64.491
+ 49.1 -68.5126 -64.5434
+ 49.15 -68.4896 -64.5944
+ 49.2 -68.4667 -64.6439
+ 49.25 -68.4439 -64.692
+ 49.3 -68.4211 -64.7388
+ 49.35 -68.3983 -64.7841
+ 49.4 -68.3756 -64.8282
+ 49.45 -68.353 -64.8709
+ 49.5 -68.3304 -64.9124
+ 49.55 -68.3078 -64.9526
+ 49.6 -68.2853 -64.9916
+ 49.65 -68.2628 -65.0294
+ 49.7 -68.2404 -65.0661
+ 49.75 -68.218 -65.1016
+ 49.8 -68.1956 -65.1359
+ 49.85 -68.1733 -65.1692
+ 49.9 -68.151 -65.2014
+ 49.95 -68.1288 -65.2325
+ 50 -68.1066 -65.2627
+ 50.05 -68.0844 -65.2918
diff --git a/tests/plot.py b/tests/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..754eb6f2cfcf336215d5b1c1ba440138833ad3c8
--- /dev/null
+++ b/tests/plot.py
@@ -0,0 +1,20 @@
+from matplotlib import pyplot
+import numpy as np
+
+raw = np.fromfile("den_syn.txt", sep=" ")
+n = raw.size/3
+data = raw.reshape(n,3)
+
+t = data[:, 0]
+soma = data[:, 1]
+dend = data[:, 2]
+
+pyplot.plot(t, soma)
+pyplot.plot(t, dend)
+
+pyplot.xlabel('time (ms)')
+pyplot.ylabel('mV')
+pyplot.xlim([t[0], t[n-1]])
+pyplot.grid()
+pyplot.show()
+
diff --git a/tests/soma_syn.txt b/tests/soma_syn.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1dd4d8fc78d6493ee8a474a7c047b0b5b05d2e6a
--- /dev/null
+++ b/tests/soma_syn.txt
@@ -0,0 +1,401 @@
+ -65
+ -64.9985
+ -64.997
+ -64.9956
+ -64.9942
+ -64.9929
+ -64.9916
+ -64.9903
+ -64.989
+ -64.9877
+ -64.9865
+ -64.9853
+ -64.9841
+ -64.983
+ -64.9818
+ -64.9807
+ -64.9796
+ -64.9785
+ -64.9774
+ -64.9764
+ -64.9754
+ -64.9744
+ -64.9734
+ -64.9724
+ -64.9715
+ -64.9706
+ -64.9697
+ -64.9688
+ -64.9679
+ -64.9671
+ -64.9663
+ -64.9655
+ -64.9647
+ -64.9639
+ -64.9632
+ -64.9625
+ -64.9618
+ -64.9611
+ -64.9604
+ -64.9598
+ -64.9592
+ -64.9586
+ -64.958
+ -64.9574
+ -64.9569
+ -64.9563
+ -64.9558
+ -64.9553
+ -64.9549
+ -64.9544
+ -64.954
+ -64.9535
+ -64.9531
+ -64.9528
+ -64.9524
+ -64.9521
+ -64.9517
+ -64.9514
+ -64.9511
+ -64.9508
+ -64.9506
+ -64.9503
+ -64.9501
+ -64.9499
+ -64.9497
+ -64.9495
+ -64.9493
+ -64.9492
+ -64.949
+ -64.9489
+ -64.9488
+ -64.9487
+ -64.9486
+ -64.9485
+ -64.9485
+ -64.9484
+ -64.9484
+ -64.9484
+ -64.9484
+ -64.9484
+ -64.9484
+ -64.9484
+ -64.9485
+ -64.9485
+ -64.9486
+ -64.9487
+ -64.9488
+ -64.9489
+ -64.949
+ -64.9491
+ -64.9492
+ -64.9493
+ -64.9495
+ -64.9496
+ -64.9498
+ -64.9499
+ -64.9501
+ -64.9503
+ -64.9505
+ -64.9507
+ -64.9509
+ -64.9511
+ -64.9513
+ -64.9515
+ -64.9518
+ -64.952
+ -64.9522
+ -64.9525
+ -64.9527
+ -64.953
+ -64.9533
+ -64.9535
+ -64.9538
+ -64.9541
+ -64.9543
+ -64.9546
+ -64.9549
+ -64.9552
+ -64.9555
+ -64.9558
+ -64.9561
+ -64.9564
+ -64.9567
+ -64.957
+ -64.9573
+ -64.9576
+ -64.9579
+ -64.9582
+ -64.9585
+ -64.9588
+ -64.9592
+ -64.9595
+ -64.9598
+ -64.9601
+ -64.9604
+ -64.9607
+ -64.961
+ -64.9614
+ -64.9617
+ -64.962
+ -64.9623
+ -64.9626
+ -64.9629
+ -64.9632
+ -64.9635
+ -64.9639
+ -64.9642
+ -64.9645
+ -64.9648
+ -64.9651
+ -64.9654
+ -64.9657
+ -64.966
+ -64.9663
+ -64.9665
+ -64.9668
+ -64.9671
+ -64.9674
+ -64.9677
+ -64.968
+ -64.9682
+ -64.9685
+ -64.9688
+ -64.9691
+ -64.9693
+ -64.9696
+ -64.9698
+ -64.9701
+ -64.9703
+ -64.9706
+ -64.9708
+ -64.9711
+ -64.9713
+ -64.9715
+ -64.9718
+ -64.972
+ -64.9722
+ -64.9724
+ -64.9727
+ -64.9729
+ -64.9731
+ -64.9733
+ -64.9735
+ -64.9737
+ -64.9739
+ -64.9741
+ -64.9743
+ -64.9744
+ -64.9746
+ -64.9748
+ -64.975
+ -64.9751
+ -64.9753
+ -64.9754
+ -64.9756
+ -64.9757
+ -64.9759
+ -64.976
+ -64.9762
+ -64.9763
+ -64.9764
+ -58.7317
+ -53.4025
+ -48.8062
+ -44.7515
+ -41.0083
+ -37.2743
+ -33.138
+ -28.0298
+ -21.1595
+ -11.5092
+ 1.77581
+ 17.7566
+ 31.6112
+ 38.5044
+ 40.1376
+ 39.9792
+ 39.2077
+ 38.0736
+ 36.6678
+ 35.0415
+ 33.2307
+ 31.2648
+ 29.1692
+ 26.9668
+ 24.6785
+ 22.3232
+ 19.9178
+ 17.4773
+ 15.0147
+ 12.5414
+ 10.0669
+ 7.59928
+ 5.14504
+ 2.70941
+ 0.296431
+ -2.09088
+ -4.45044
+ -6.78095
+ -9.08183
+ -11.3532
+ -13.5957
+ -15.8107
+ -18.0002
+ -20.167
+ -22.315
+ -24.449
+ -26.5755
+ -28.7028
+ -30.8412
+ -33.0035
+ -35.2054
+ -37.4655
+ -39.8054
+ -42.2486
+ -44.8193
+ -47.5387
+ -50.4185
+ -53.4512
+ -56.5966
+ -59.7684
+ -62.8311
+ -65.6181
+ -67.9773
+ -69.8235
+ -71.1623
+ -72.0707
+ -72.6551
+ -73.0161
+ -73.2319
+ -73.3571
+ -73.4271
+ -73.4638
+ -73.4804
+ -73.485
+ -73.4823
+ -73.4751
+ -73.4651
+ -73.4532
+ -73.4402
+ -73.4264
+ -73.4121
+ -73.3973
+ -73.3822
+ -73.3668
+ -73.3512
+ -73.3354
+ -73.3193
+ -73.3031
+ -73.2868
+ -73.2702
+ -73.2535
+ -73.2366
+ -73.2195
+ -73.2022
+ -73.1847
+ -73.1671
+ -73.1492
+ -73.1312
+ -73.113
+ -73.0946
+ -73.076
+ -73.0572
+ -73.0382
+ -73.019
+ -72.9995
+ -72.9799
+ -72.96
+ -72.9399
+ -72.9197
+ -72.8991
+ -72.8784
+ -72.8574
+ -72.8362
+ -72.8148
+ -72.7931
+ -72.7712
+ -72.7491
+ -72.7267
+ -72.7041
+ -72.6812
+ -72.6581
+ -72.6347
+ -72.6111
+ -72.5873
+ -72.5632
+ -72.5388
+ -72.5142
+ -72.4894
+ -72.4642
+ -72.4389
+ -72.4132
+ -72.3873
+ -72.3612
+ -72.3348
+ -72.3081
+ -72.2812
+ -72.2541
+ -72.2266
+ -72.1989
+ -72.171
+ -72.1428
+ -72.1143
+ -72.0856
+ -72.0566
+ -72.0274
+ -71.9979
+ -71.9682
+ -71.9382
+ -71.908
+ -71.8775
+ -71.8468
+ -71.8159
+ -71.7846
+ -71.7532
+ -71.7215
+ -71.6896
+ -71.6574
+ -71.625
+ -71.5924
+ -71.5595
+ -71.5265
+ -71.4932
+ -71.4596
+ -71.4259
+ -71.3919
+ -71.3578
+ -71.3234
+ -71.2888
+ -71.254
+ -71.2191
+ -71.1839
+ -71.1485
+ -71.113
+ -71.0772
+ -71.0413
+ -71.0052
+ -70.9689
+ -70.9325
+ -70.8958
+ -70.8591
+ -70.8221
+ -70.785
+ -70.7478
+ -70.7104
+ -70.6729
+ -70.6352
+ -70.5974
+ -70.5595
+ -70.5214
+ -70.4832
+ -70.4449
+ -70.4065
+ -70.368
+ -70.3293
+ -70.2906
+ -70.2518
+ -70.2129
+ -70.1739
+ -70.1348
+ -70.0956
diff --git a/tests/test_event_queue.cpp b/tests/test_event_queue.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3297784d506f6eff10078e249388942b85d2240f
--- /dev/null
+++ b/tests/test_event_queue.cpp
@@ -0,0 +1,99 @@
+#include "gtest.h"
+
+#include <vector>
+
+#include <event_queue.hpp>
+
+TEST(event_queue, push)
+{
+ using namespace nest::mc;
+
+ event_queue q;
+
+ q.push({1u, 2.f, 2.f});
+ q.push({4u, 1.f, 2.f});
+ q.push({8u, 20.f, 2.f});
+ q.push({2u, 8.f, 2.f});
+
+ std::vector<float> times;
+ while(q.size()) {
+ times.push_back(
+ q.pop_if_before(std::numeric_limits<float>::max()).second.time
+ );
+ }
+
+ //std::copy(times.begin(), times.end(), std::ostream_iterator<float>(std::cout, ","));
+ //std::cout << "\n";
+ EXPECT_TRUE(std::is_sorted(times.begin(), times.end()));
+}
+
+TEST(event_queue, push_range)
+{
+ using namespace nest::mc;
+
+ local_event events[] = {
+ {1u, 2.f, 2.f},
+ {4u, 1.f, 2.f},
+ {8u, 20.f, 2.f},
+ {2u, 8.f, 2.f}
+ };
+
+ event_queue q;
+ q.push(std::begin(events), std::end(events));
+
+ std::vector<float> times;
+ while(q.size()) {
+ times.push_back(
+ q.pop_if_before(std::numeric_limits<float>::max()).second.time
+ );
+ }
+
+ EXPECT_TRUE(std::is_sorted(times.begin(), times.end()));
+}
+
+TEST(event_queue, pop_if_before)
+{
+ using namespace nest::mc;
+
+ local_event events[] = {
+ {1u, 1.f, 2.f},
+ {2u, 2.f, 2.f},
+ {3u, 3.f, 2.f},
+ {4u, 4.f, 2.f}
+ };
+
+ event_queue q;
+ q.push(std::begin(events), std::end(events));
+
+ EXPECT_EQ(q.size(), 4u);
+
+ auto e1 = q.pop_if_before(0.);
+ EXPECT_FALSE(e1.first);
+ EXPECT_EQ(q.size(), 4u);
+
+ auto e2 = q.pop_if_before(5.);
+ EXPECT_TRUE(e2.first);
+ EXPECT_EQ(e2.second.target, 1u);
+ EXPECT_EQ(q.size(), 3u);
+
+ auto e3 = q.pop_if_before(5.);
+ EXPECT_TRUE(e3.first);
+ EXPECT_EQ(e3.second.target, 2u);
+ EXPECT_EQ(q.size(), 2u);
+
+ auto e4 = q.pop_if_before(2.5);
+ EXPECT_FALSE(e4.first);
+ EXPECT_EQ(q.size(), 2u);
+
+ auto e5 = q.pop_if_before(5.);
+ EXPECT_TRUE(e5.first);
+ EXPECT_EQ(e5.second.target, 3u);
+ EXPECT_EQ(q.size(), 1u);
+
+ q.pop_if_before(5.);
+ EXPECT_EQ(q.size(), 0u);
+
+ // empty queue should always return "false"
+ auto e6 = q.pop_if_before(100.);
+ EXPECT_FALSE(e6.first);
+}
diff --git a/tests/test_synapses.cpp b/tests/test_synapses.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5a8f2ffac903289663dc21e105523a593bd695a2
--- /dev/null
+++ b/tests/test_synapses.cpp
@@ -0,0 +1,83 @@
+#include "gtest.h"
+#include "util.hpp"
+
+#include <cell.hpp>
+#include <fvm.hpp>
+
+// compares results with those generated by nrn/ball_and_stick.py
+TEST(synapses, add_to_cell)
+{
+ using namespace nest::mc;
+
+ 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());
+
+ cell.add_synapse({0, 0.1});
+ cell.add_synapse({1, 0.2});
+ cell.add_synapse({0, 0.3});
+
+ EXPECT_EQ(3u, cell.synapses().size());
+ EXPECT_EQ(cell.synapses()[0].segment, 0);
+ EXPECT_EQ(cell.synapses()[0].position, 0.1);
+ EXPECT_EQ(cell.synapses()[1].segment, 1);
+ EXPECT_EQ(cell.synapses()[1].position, 0.2);
+ EXPECT_EQ(cell.synapses()[2].segment, 0);
+ EXPECT_EQ(cell.synapses()[2].position, 0.3);
+}
+
+// compares results with those generated by nrn/ball_and_stick.py
+TEST(synapses, basic_state)
+{
+ using namespace nest::mc;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ using synapse_type = mechanisms::ExpSyn::mechanism_ExpSyn<double, int>;
+ auto num_syn = 4;
+ synapse_type::index_type indexes(num_syn);
+ synapse_type::vector_type voltage(num_syn, -65.0);
+ synapse_type::vector_type current(num_syn, 1.0);
+ auto mech = mechanisms::make_mechanism<synapse_type>( voltage, current, indexes );
+
+ auto ptr = dynamic_cast<synapse_type*>(mech.get());
+
+ // parameters initialized to default values
+ for(auto e : ptr->e) {
+ EXPECT_EQ(e, 0.);
+ }
+ for(auto tau : ptr->tau) {
+ EXPECT_EQ(tau, 2.0);
+ }
+
+ // current and voltage vectors correctly hooked up
+ for(auto v : ptr->vec_v_) {
+ EXPECT_EQ(v, -65.);
+ }
+ for(auto i : ptr->vec_i_) {
+ EXPECT_EQ(i, 1.0);
+ }
+
+ // should be initialized to NaN
+ for(auto g : ptr->g) {
+ EXPECT_NE(g, g);
+ }
+
+ // initialize state then check g has been set to zero
+ ptr->nrn_init();
+ for(auto g : ptr->g) {
+ EXPECT_EQ(g, 0.);
+ }
+
+ // call net_receive on two of the synapses
+ ptr->net_receive(1, 3.14);
+ ptr->net_receive(3, 1.04);
+ EXPECT_EQ(ptr->g[1], 3.14);
+ EXPECT_EQ(ptr->g[3], 1.04);
+}
diff --git a/tests/validate_synapses.cpp b/tests/validate_synapses.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7b882a689c94cb5e6b70e8f89115917fc7a6ba0d
--- /dev/null
+++ b/tests/validate_synapses.cpp
@@ -0,0 +1,160 @@
+#include <fstream>
+
+#include "gtest.h"
+#include "util.hpp"
+
+#include <cell.hpp>
+#include <fvm.hpp>
+
+#include <json/src/json.hpp>
+
+
+// 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,
+ nlohmann::json& measurements
+ ) {
+ n_comparments = nseg;
+ baseline_spikes = {
+ measurements["soma"]["spikes"],
+ measurements["dend"]["spikes"]
+ };
+ thresh = {
+ measurements["soma"]["thresh"],
+ measurements["dend"]["thresh"]
+ };
+ for(auto i=0u; i<v.size(); ++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]));
+ }
+ }
+};
+
+// compares results with those generated by nrn/simple_synapse.py
+TEST(simple_synapse, 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);
+ soma->mechanism("membrane").set("r_L", 100);
+
+ // add stimulus
+ cell.add_synapse({1, 0.5});
+
+ // load data from file
+ auto cell_data = testing::load_spike_data("../nrn/simple_synapse.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 = 50.; // ms
+
+ 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(2);
+
+ // 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
+
+ // add the 3 spike events to the queue
+ model.queue().push({0u, 10.0, 0.04});
+ model.queue().push({0u, 20.0, 0.04});
+ model.queue().push({0u, 40.0, 0.04});
+
+ // 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);
+ v[0].push_back(model.voltage()[soma_comp]);
+ v[1].push_back(model.voltage()[dend_comp]);
+ double t = 0.;
+ while(t < tfinal) {
+ t += dt;
+ model.advance_to(t, dt);
+ // save voltage at soma
+ v[0].push_back(model.voltage()[soma_comp]);
+ v[1].push_back(model.voltage()[dend_comp]);
+ }
+
+ results.push_back( {num_compartments, dt, v, measurements} );
+ }
+
+ // 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
+ 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 j=0; j<2; ++j) {
+ EXPECT_TRUE(
+ results.back().comparisons[j].max_relative_error()
+ < results.front().comparisons[j].max_relative_error()
+ );
+ }
+
+ // 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<2; ++j) {
+ EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<tol);
+ }
+}
+