diff --git a/.gitignore b/.gitignore
index f60f1b41e5cdbbb2ac445207ed0d2e41176e0104..a18a8b7d7af992216fe6d6bd229358fe0aa88b95 100644
--- a/.gitignore
+++ b/.gitignore
@@ -50,3 +50,17 @@ Makefile
# mechanism implementations generated my modparser
include/mechanisms
+
+# external build stuff
+external/bin
+external/modparser-build
+external/modparser-configure
+external/modparser-done
+external/modparser-download
+external/modparser-install
+external/modparser-mkdir
+external/modparser-patch
+external/modparser-update
+external/tmp
+mechanisms/*.hpp
+
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 dcfb69f8cfc435418dfb95348d613c9830e6ad6f..588ca1a5ea28ef04d17b318e754d479e5489eb9a 160000
--- a/external/modparser
+++ b/external/modparser
@@ -1 +1 @@
-Subproject commit dcfb69f8cfc435418dfb95348d613c9830e6ad6f
+Subproject commit 588ca1a5ea28ef04d17b318e754d479e5489eb9a
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/CMakeLists.txt b/mechanisms/CMakeLists.txt
index ac68a5e963e42acf6b67557ca524d61572fc4902..203dedff949356b1a715146baf0ad6d6714aac56 100644
--- a/mechanisms/CMakeLists.txt
+++ b/mechanisms/CMakeLists.txt
@@ -1,4 +1,4 @@
-set(mechanisms pas hh)
+set(mechanisms pas hh expsyn)
set(modcc "${CMAKE_BINARY_DIR}/external/bin/modcc")
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 100644
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 98b871cb193ee09a595416128e6342b84fa028b9..2f336f8454bd8e3e4f4675865665c6bbb42ec654 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 fab1163b1bbf9b8c43dbf954b4a7081c4e9dcc7f..1e04fc2c3bcc4e0fe887c7048997afd272fa4bff 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 d0c6f90272ee5f48fc7727a1aab519e3bebe5fdc..4a908bba49887948e37bfa7fe3409b19df90bbb3 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -2,20 +2,25 @@
#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 <segment.hpp>
+#include <stimulus.hpp>
+#include <util.hpp>
+
#include <vector/include/Vector.hpp>
-#include "algorithms.hpp"
-#include "cell.hpp"
-#include "ion.hpp"
-#include "math.hpp"
-#include "matrix.hpp"
-#include "mechanism.hpp"
-#include "mechanism_interface.hpp"
-#include "segment.hpp"
-#include "stimulus.hpp"
-#include "util.hpp"
+#include <mechanisms/expsyn.hpp>
namespace nest {
namespace mc {
@@ -44,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);
@@ -56,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_;
}
@@ -83,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();
@@ -129,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
@@ -160,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_;
@@ -167,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_;
};
////////////////////////////////////////////////////////////////////////////////
@@ -292,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];
@@ -305,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
@@ -355,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
@@ -376,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>
@@ -405,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;
@@ -416,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)
@@ -460,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 e86c9b897b66a529252721b739ae20f96cfa16ae..9e2c9e7d10d68c2f973cffaaf005ceb939fdcf41 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -18,10 +18,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;
@@ -84,41 +86,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 7a06a7330e06bb14000e7c262c3db55fe4dc7135..c9c801fd198c764d61f3acf653055c8957daf00e 100644
--- a/src/mechanism.hpp
+++ b/src/mechanism.hpp
@@ -35,6 +35,7 @@ public:
: vec_v_(vec_v)
, vec_i_(vec_i)
, node_index_(node_index)
+ , vec_area_(nullptr, 0)
{}
std::size_t size() const
@@ -63,14 +64,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>
@@ -87,5 +94,5 @@ make_mechanism(
}
} // namespace mechanisms
-} // namespace nest
} // namespace mc
+} // namespace nest
diff --git a/src/tree.hpp b/src/tree.hpp
index be96e3619424ec54b143e694108188419802e67b..3d76a9dc1e8171f62026320a98f71078e91f1dc5 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -20,7 +20,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;
@@ -163,22 +164,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_;
}
@@ -288,7 +289,7 @@ class tree {
int new_node,
int old_node,
int parent_node,
- index_view p,
+ view_type p,
tree const& old_tree
)
{
@@ -350,9 +351,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 51d0aa9c97644c3e3ade1d35f6346c0c44bf7c77..9b61220058ee4f83d4102c68102dd85c7bfa9ece 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/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/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);
+ }
+}
+