diff --git a/.gitignore b/.gitignore
index c23dfde0164c2ce71be3a2113f47b2f96e793753..f60f1b41e5cdbbb2ac445207ed0d2e41176e0104 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,6 +32,7 @@
*.dot
*.pdf
*.jpg
+*.dat
# latex output
*.aux
@@ -47,3 +48,5 @@ CMakeCache.txt
cmake_install.cmake
Makefile
+# mechanism implementations generated my modparser
+include/mechanisms
diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index b8378fb82ab46999a5a31f6f5ff516f0c6cc5af6..55d6af4d257348526713c1fdc381c4c2792651fd 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -43,6 +43,8 @@ flags = [
'src',
'-I',
'.',
+ '-I',
+ 'include',
# '-isystem',
# '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
# '-I',
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 61278f34ff79a0b1dc688835f369df0b144926e0..57991c886fc85752b0d7a791375ba8004ce69012 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,6 +14,7 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS "YES")
set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
+include_directories(${CMAKE_SOURCE_DIR}/include)
include_directories(${CMAKE_SOURCE_DIR}/src)
include_directories(${CMAKE_SOURCE_DIR})
diff --git a/data/hh.mod b/data/hh.mod
old mode 100755
new mode 100644
index 6020be0c1f6d99ab90a3639916eb32cdb1fab155..975d601d09491f6b6418200b795368dbd9b14f7e
--- a/data/hh.mod
+++ b/data/hh.mod
@@ -1,17 +1,17 @@
TITLE hh.mod squid sodium, potassium, and leak channels
-COMMENT
- This is the original Hodgkin-Huxley treatment for the set of sodium,
- potassium, and leakage channels found in the squid giant axon membrane.
- ("A quantitative description of membrane current and its application
- conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).)
- Membrane voltage is in absolute mV and has been reversed in polarity
- from the original HH convention and shifted to reflect a resting potential
- of -65 mV.
- Remember to set celsius=6.3 (or whatever) in your HOC file.
- See squid.hoc for an example of a simulation using this model.
- SW Jaslove 6 March, 1992
-ENDCOMMENT
+: COMMENT
+: This is the original Hodgkin-Huxley treatment for the set of sodium,
+: potassium, and leakage channels found in the squid giant axon membrane.
+: ("A quantitative description of membrane current and its application
+: conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).)
+: Membrane voltage is in absolute mV and has been reversed in polarity
+: from the original HH convention and shifted to reflect a resting potential
+: of -65 mV.
+: Remember to set celsius=6.3 (or whatever) in your HOC file.
+: See squid.hoc for an example of a simulation using this model.
+: SW Jaslove 6 March, 1992
+: ENDCOMMENT
UNITS {
(mA) = (milliamp)
@@ -26,14 +26,14 @@ NEURON {
NONSPECIFIC_CURRENT il
RANGE gnabar, gkbar, gl, el, gna, gk
GLOBAL minf, hinf, ninf, mtau, htau, ntau
- THREADSAFE : assigned GLOBALs will be per thread
}
PARAMETER {
- gnabar = .12 (S/cm2) <0,1e9>
- gkbar = .036 (S/cm2) <0,1e9>
- gl = .0003 (S/cm2) <0,1e9>
+ gnabar = .12 (S/cm2) : <0,1e9>
+ gkbar = .036 (S/cm2) : <0,1e9>
+ gl = .0003 (S/cm2) : <0,1e9>
el = -54.3 (mV)
+ celsius = 37 (degC)
}
STATE {
@@ -42,20 +42,17 @@ STATE {
ASSIGNED {
v (mV)
- celsius (degC)
- ena (mV)
- ek (mV)
gna (S/cm2)
gk (S/cm2)
- ina (mA/cm2)
- ik (mA/cm2)
- il (mA/cm2)
- minf hinf ninf
- mtau (ms) htau (ms) ntau (ms)
+ minf
+ hinf
+ ninf
+ mtau (ms)
+ htau (ms)
+ ntau (ms)
}
-? currents
BREAKPOINT {
SOLVE states METHOD cnexp
gna = gnabar*m*m*m*h
@@ -72,18 +69,15 @@ INITIAL {
n = ninf
}
-DERIVATIVE states
-{
+DERIVATIVE states {
rates(v)
m' = (minf-m)/mtau
h' = (hinf-h)/htau
n' = (ninf-n)/ntau
}
-PROCEDURE rates(v(mV))
+PROCEDURE rates(v) :(mV))
{
- :Computes rate and other constants at current v.
- :Call once from HOC to initialize inf at resting v.
LOCAL alpha, beta, sum, q10
q10 = 3^((celsius - 6.3)/10)
@@ -103,18 +97,28 @@ PROCEDURE rates(v(mV))
hinf = alpha/sum
:"n" potassium activation system
- alpha = .01*vtrap(-(v+55),10)
+ alpha = .01*vtrap(-(v+55),10)
beta = .125*exp(-(v+65)/80)
sum = alpha + beta
ntau = 1/(q10*sum)
ninf = alpha/sum
}
-FUNCTION vtrap(x,y) { :Traps for 0 in denominator of rate eqns.
- if (fabs(x/y) < 1e-6) {
- vtrap = y*(1 - x/y/2)
- }else{
- vtrap = x/(exp(x/y) - 1)
- }
+FUNCTION vtrap(x,y) {
+ : Traps for 0 in denominator of rate eqns.
+
+ : Disabled for function inlining...
+ : there is a very good case for making vtrap a builtin function for the language
+ : - it is ubiquitous: used in so many models of ion channels
+ : - platform specific optimizations written in assembly/intrinsics
+ : required to facilitate vectorization/minimize branch misses
+
+ :if (fabs(x/y) < 1e-6) {
+ : vtrap = y*(1 - x/y/2)
+ :}else{
+ : vtrap = x/(exp(x/y) - 1)
+ :}
+
+ vtrap = x/(exp(x/y) - 1)
}
diff --git a/data/passive.mod b/data/passive.mod
new file mode 100644
index 0000000000000000000000000000000000000000..413bd250c8e9a077287c97b82cb9dd2cd44927cb
--- /dev/null
+++ b/data/passive.mod
@@ -0,0 +1,28 @@
+TITLE passive membrane channel
+
+UNITS {
+ (mV) = (millivolt)
+ (mA) = (milliamp)
+ (S) = (siemens)
+}
+
+NEURON {
+ SUFFIX pas
+ NONSPECIFIC_CURRENT i
+ RANGE g, e
+}
+
+INITIAL {}
+
+PARAMETER {
+ g = .001 (S/cm2) :<0,1e9>
+ e = -70 (mV)
+}
+
+ASSIGNED {
+ v (mV)
+}
+
+BREAKPOINT {
+ i = g*(v - e)
+}
diff --git a/docs/report.tex b/docs/report.tex
index 54d93c8dcfd449fe70b095113e2a1930b12cf4e7..0f816e1c25a7aa22b7199df27a0cd3d8216f9dc2 100644
--- a/docs/report.tex
+++ b/docs/report.tex
@@ -72,6 +72,9 @@
\newcommand{\dder}[2]{\frac{\deriv{#1}}{\deriv{#2}}}
\newcommand{\vv}[1]{\bm{#1}\xspace}
+\newcommand{\unit}[1]{\left[{#1}\right]}
+\newcommand{\txtunit}[1]{$\left[{#1}\right]$}
+
%----------------------------------------------------------------------------------------
% ARTICLE INFORMATION
%----------------------------------------------------------------------------------------
diff --git a/docs/symbols.tex b/docs/symbols.tex
index 19db4a374eb9a3697bcb8971ca2cca9e0b40877e..e12b55820da048ffef27bb790fca4e0b363ba23e 100644
--- a/docs/symbols.tex
+++ b/docs/symbols.tex
@@ -1,3 +1,113 @@
+%-------------------------------------------------------------------------------
+%\subsubsection{Balancing Units}
+%-------------------------------------------------------------------------------
+Ensuring that units are balanced and correct requires care.
+Take the description of the nonlinear system of ODEs that arises from the finite volume discretisation
+\begin{equation}
+ \label{eq:linsys_FV}
+ V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\frac{\Delta t \alpha_{ij}}{\sigma_i} (V_i^{k+1}-V_j^{k+1})}
+ = V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e).
+\end{equation}
+The choice of units for a parameter, e.g. $\mu m^2$ or $m^2$ for the area $\sigma_{ij}$, introduces a constant of proportionality wherever it is used ($10^{-12}$ in the case of $\mu m^2 \rightarrow m^2$).
+Wherever terms are added in \eq{eq:linsys_FV} the units must be checked, and constants of proportionality balanced.
+
+First, appropriate units for each of the parameters and variables are chosen in~\tbl{tbl:units}.
+We try to use the same units as NEURON, except for the specific membrane capacitance $c_m$, for which $F\cdot m^{-2}$ is used in place of $nF\cdot mm^{-2}$.
+In \eq{eq:linsys_FV} we choose units of $mV \equiv 10^{-3}V$ for each term because of the $V_i$ terms on either side of the equation.
+
+\begin{table}[hp!]
+\begin{tabular}{lllr}
+ \hline
+ term & units & normalized units & NEURON \\
+ \hline
+ $t$ & $ms$ & $10^{-3} \cdot s$ & yes \\
+ $V$ & $mV$ & $10^{-3} \cdot V$ & yes \\
+ $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 \\
+ $\overline{g}$ & $S\cdot cm^{-2}$ & $10^{4} \cdot A\cdot V^{-1}\cdot m^{-2}$ & yes \\
+ $I_e$ & $nA$ & $10^{-9} \cdot A$ & yes \\
+ \hline
+\end{tabular}
+\caption{The units chosen for parameters and variables in NEST MC. The NEURON column indicates whether the same units have been used as NEURON.}
+\label{tbl:units}
+\end{table}
+
+%------------------------------------------
+\subsubsection{current terms}
+%------------------------------------------
+Membrane current is calculated as follows $i_m = \overline{g}(E-V)$, with units
+\begin{align}
+ \unit{ i_m } &= \unit{ \overline{g} } \unit{ V } \nonumber \\
+ &= 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 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 \\
+ &= 10^{-9}\cdot A \cdot 10^{12} \cdot m^{-2} \nonumber \\
+ &= 10^{3} \cdot A \cdot m ^{-2}, \label{eq:ie_unit}
+\end{align}
+which must be scaled by $10^2$ to match $i_m$ in \eq{eq:im_unit}.
+
+The units for the flux coefficent can be calculated as follows:
+\begin{align}
+ \unit{ \frac{\Delta t}{c_m} } &= 10^{-3} \cdot s \cdot s^{-1}\cdot A^{-1}\cdot V\cdot m^2 \nonumber \\
+ &= 10^{-3} \cdot A^{-1} \cdot V\cdot m^2. \label{eq:dtcm_unit}
+\end{align}
+From \eq{eq:im_unit} and \eq{eq:dtcm_unit} that the units of the full current term are
+\begin{align}
+ \unit{ \frac{\Delta t}{c_m}\left(i_m - i_e\right) }
+ &= 10^{-3} \cdot A^{-1} \cdot V\cdot m^2 \cdot 10 \cdot A \cdot m^{-2} \nonumber \\
+ &= 10^{-2} \cdot V,
+\end{align}
+which must be scaled by $10$ to match the units of $mV\equiv10^{-3}V$.
+%------------------------------------------
+\subsubsection{flux terms}
+%------------------------------------------
+The coefficients in the linear system have the units
+\begin{equation}
+ \unit{ \frac{\Delta t\alpha_{ij}}{\sigma_i} }
+ =
+ \unit{ \frac{\Delta t \sigma_{ij} } {c_m r_L \Delta x_{ij} \sigma_i} }
+ =
+ \unit{ \frac{\Delta t } {c_m r_L \Delta x_{ij} } },
+\end{equation}
+where we we simplify by noting that $\unit{\sigma_{ij}}=\unit{\sigma_i}$.
+The units of the term $c_m r_L$ on the denominator are calculated as follows
+\begin{align}
+ \unit{c_m r_L}
+ &= s \cdot A \cdot V^{-1} \cdot m^{-2} \cdot 10^{-2} \cdot A^{-1} \cdot V \cdot m \nonumber \\
+ &= 10^{-2} \cdot s \cdot m^{-1},
+\end{align}
+so the units of the denominator are
+\begin{align}
+ \unit{c_m r_L \Delta x_{ij}}
+ &= 10^{-2} \cdot s \cdot m^{-1} \cdot 10^{-6} \cdot m \nonumber \\
+ &= 10^{-8} \cdot s,
+\end{align}
+and hence
+\begin{align}
+ \unit{\frac{\Delta t } {c_m r_L \Delta x_{ij} }}
+ &= 10^{8} \cdot s^{-1} \cdot 10^{-3} \cdot s \nonumber \\
+ &= 10^{5}.
+\end{align}
+
+So, the terms with $\alpha_{ij}$ must be scaled by $10^5$ to match the units of $mV$.
+%------------------------------------------
+\subsubsection{discretization with scaling}
+%------------------------------------------
+Here is something that I wish the NEURON documentation had provided:
+\begin{align}
+& V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {10^5 \cdot \frac{\Delta t \alpha_{ij}}{\sigma_i} (V_i^{k+1}-V_j^{k+1})} \nonumber \\
+& = V_i^k - 10\cdot \frac{\Delta t}{c_m}(i_m^{k} - 10^2\cdot I_e/\sigma_i).
+\end{align}
+%------------------------------------------
+\subsection{Supplementary Unit Information}
+%------------------------------------------
+Here is some information about units scraped from Wikipedia for convenience.
+
\begin{table*}[htp!]
\begin{center}
@@ -15,6 +125,8 @@
\hline
\end{tabular}
+ \vspace{20pt}
+
\begin{tabular}{llll}
\hline
symbol & unit & equivalents & SI base \\
@@ -44,27 +156,4 @@
\end{center}
\caption{Symbols and quantities.}
\end{table*}
-%-------------------------------------------------------------------------------
-\subsubsection{Units}
-%-------------------------------------------------------------------------------
-Reffering to the cable equation first defined in~\eq{eq:cable}
-\begin{equation*}
- c_m \pder{V}{t} = \frac{1}{2\pi a r_{L}} \pder{}{x} \left( a^2 \pder{V}{x} \right) - i_m + i_e,
-\end{equation*}
-
-If the units are taken to be
-\begin{itemize}
- \item $c_m = F\cdot cm^{-2}$
- \item $V = V$
- \item $a = cm$
- \item $r_L = \Omega\cdot cm$
-\end{itemize}
-Then the units of each term in equation are $A\cdot cm^{-2}$.
-In practice, the units above are not used, for example distances are usually measured in $\mu m$ and areas in $cm^2$.
-But if we work term-by-term, scaling for these factors is manageable.
-
-A useful identity to use when performing the dimensional analysis relates capacitance and resistance
-\begin{equation*}
- 1~F = 1~\Omega^{-1} \cdot s
-\end{equation*}
diff --git a/include/mechanisms/hh_verb.hpp b/include/mechanisms/hh_verb.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8b67a69d56a0f61b2db8baaa37d03fc34d7e9687
--- /dev/null
+++ b/include/mechanisms/hh_verb.hpp
@@ -0,0 +1,284 @@
+#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
new file mode 100755
index 0000000000000000000000000000000000000000..9fe83f04d1b2775eae25f0e4ec7030cd89ee6163
--- /dev/null
+++ b/mechanisms/generate.sh
@@ -0,0 +1,4 @@
+for mech in pas hh
+do
+ modcc -t cpu -o ../include/mechanisms/$mech.hpp ./mod/$mech.mod
+done
diff --git a/mechanisms/mod/hh.mod b/mechanisms/mod/hh.mod
new file mode 100644
index 0000000000000000000000000000000000000000..8b23f44a4dc3e27cdbb2d78a161b9504c8651073
--- /dev/null
+++ b/mechanisms/mod/hh.mod
@@ -0,0 +1,125 @@
+TITLE hh.mod squid sodium, potassium, and leak channels
+
+: COMMENT
+: This is the original Hodgkin-Huxley treatment for the set of sodium,
+: potassium, and leakage channels found in the squid giant axon membrane.
+: ("A quantitative description of membrane current and its application
+: conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).)
+: Membrane voltage is in absolute mV and has been reversed in polarity
+: from the original HH convention and shifted to reflect a resting potential
+: of -65 mV.
+: Remember to set celsius=6.3 (or whatever) in your HOC file.
+: See squid.hoc for an example of a simulation using this model.
+: SW Jaslove 6 March, 1992
+: ENDCOMMENT
+
+UNITS {
+ (mA) = (milliamp)
+ (mV) = (millivolt)
+ (S) = (siemens)
+}
+
+NEURON {
+ SUFFIX hh
+ USEION na READ ena WRITE ina
+ USEION k READ ek WRITE ik
+ NONSPECIFIC_CURRENT il
+ RANGE gnabar, gkbar, gl, el, gna, gk
+ GLOBAL minf, hinf, ninf, mtau, htau, ntau
+}
+
+PARAMETER {
+ gnabar = .12 (S/cm2) : <0,1e9>
+ gkbar = .036 (S/cm2) : <0,1e9>
+ gl = .0003 (S/cm2) : <0,1e9>
+ el = -54.3 (mV)
+ celsius = 6.3 (degC)
+ :celsius = 37 (degC)
+}
+
+STATE {
+ m h n
+}
+
+ASSIGNED {
+ v (mV)
+
+ gna (S/cm2)
+ gk (S/cm2)
+ minf
+ hinf
+ ninf
+ mtau (ms)
+ htau (ms)
+ ntau (ms)
+}
+
+BREAKPOINT {
+ SOLVE states METHOD cnexp
+ gna = gnabar*m*m*m*h
+ ina = gna*(v - ena)
+ gk = gkbar*n*n*n*n
+ ik = gk*(v - ek)
+ il = gl*(v - el)
+}
+
+INITIAL {
+ rates(v)
+ m = minf
+ h = hinf
+ n = ninf
+}
+
+DERIVATIVE states {
+ rates(v)
+ m' = (minf-m)/mtau
+ h' = (hinf-h)/htau
+ n' = (ninf-n)/ntau
+}
+
+PROCEDURE rates(v) :(mV))
+{
+ LOCAL alpha, beta, sum, q10
+
+ q10 = 3^((celsius - 6.3)/10)
+
+ :"m" sodium activation system
+ alpha = .1 * vtrap(-(v+40),10)
+ beta = 4 * exp(-(v+65)/18)
+ sum = alpha + beta
+ mtau = 1/(q10*sum)
+ minf = alpha/sum
+
+ :"h" sodium inactivation system
+ alpha = .07 * exp(-(v+65)/20)
+ beta = 1 / (exp(-(v+35)/10) + 1)
+ sum = alpha + beta
+ htau = 1/(q10*sum)
+ hinf = alpha/sum
+
+ :"n" potassium activation system
+ alpha = .01*vtrap(-(v+55),10)
+ beta = .125*exp(-(v+65)/80)
+ sum = alpha + beta
+ ntau = 1/(q10*sum)
+ ninf = alpha/sum
+}
+
+FUNCTION vtrap(x,y) {
+ : Traps for 0 in denominator of rate eqns.
+
+ : Disabled for function inlining...
+ : there is a very good case for making vtrap a builtin function for the language
+ : - it is ubiquitous: used in so many models of ion channels
+ : - platform specific optimizations written in assembly/intrinsics
+ : required to facilitate vectorization/minimize branch misses
+
+ :if (fabs(x/y) < 1e-6) {
+ : vtrap = y*(1 - x/y/2)
+ :}else{
+ : vtrap = x/(exp(x/y) - 1)
+ :}
+
+ vtrap = x/(exp(x/y) - 1)
+}
+
diff --git a/mechanisms/mod/pas.mod b/mechanisms/mod/pas.mod
new file mode 100644
index 0000000000000000000000000000000000000000..b7dc2ad6d2a343bb2333c5b8ecfec17060150865
--- /dev/null
+++ b/mechanisms/mod/pas.mod
@@ -0,0 +1,28 @@
+TITLE passive membrane channel
+
+UNITS {
+ (mV) = (millivolt)
+ (mA) = (milliamp)
+ (S) = (siemens)
+}
+
+NEURON {
+ SUFFIX pas
+ NONSPECIFIC_CURRENT i
+ RANGE g, e
+}
+
+INITIAL {}
+
+PARAMETER {
+ g = .001 (S/cm2) :<0,1e9>
+ e = -65 (mV) : we use -65 for the ball and stick model, instead of Neuron default of -70
+}
+
+ASSIGNED {
+ v (mV)
+}
+
+BREAKPOINT {
+ i = g*(v - e)
+}
diff --git a/nrn/ball_and_stick.py b/nrn/ball_and_stick.py
new file mode 100644
index 0000000000000000000000000000000000000000..5fe977fe594ada556ebce1da09c9a15c3f3b03e3
--- /dev/null
+++ b/nrn/ball_and_stick.py
@@ -0,0 +1,138 @@
+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
+
+stim = h.IClamp(dend(1))
+stim.delay = 5
+stim.dur = 80
+stim.amp = 3*0.1
+
+if args.plot :
+ pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
+ pyplot.grid()
+
+simdur = 100.0
+h.tstop = simdur
+h.dt = 0.001
+
+start = timer()
+results = []
+for nseg in [5, 11, 21, 41, 81, 161] :
+
+ print 'simulation with ', nseg, ' compartments in dendrite...'
+
+ dend.nseg=nseg
+
+ # record voltages
+ v_soma = h.Vector() # soma
+ v_dend = h.Vector() # middle of dendrite
+ v_clamp= h.Vector() # end of dendrite at clamp location
+
+ v_soma.record( soma(0.5)._ref_v)
+ v_dend.record( dend(0.5)._ref_v)
+ v_clamp.record(dend(1.0)._ref_v)
+
+ # record spikes
+ # this is a bit verbose, no?
+ spike_counter_soma = h.APCount(soma(0.5))
+ spike_counter_soma.thresh = 0
+ spike_counter_dend = h.APCount(dend(0.5))
+ spike_counter_dend.thresh = -10
+ spike_counter_clamp = h.APCount(dend(1.0))
+ spike_counter_clamp.thresh = 10
+
+ spikes_soma = h.Vector() # soma
+ spikes_dend = h.Vector() # middle of dendrite
+ spikes_clamp= h.Vector() # end of dendrite at clamp location
+
+ spike_counter_soma.record(spikes_soma)
+ spike_counter_dend.record(spikes_dend)
+ spike_counter_clamp.record(spikes_clamp)
+
+ # record time stamps
+ t_vec = h.Vector()
+ t_vec.record(h._ref_t)
+
+ # finally it's time to run the simulation
+ h.run()
+
+ results.append(
+ {
+ "nseg": nseg,
+ "dt" : h.dt,
+ "measurements": {
+ "soma" : {
+ "thresh" : spike_counter_soma.thresh,
+ "spikes" : spikes_soma.to_python()
+ },
+ "dend" : {
+ "thresh" : spike_counter_dend.thresh,
+ "spikes" : spikes_dend.to_python()
+ },
+ "clamp" : {
+ "thresh" : spike_counter_clamp.thresh,
+ "spikes" : spikes_clamp.to_python()
+ }
+ }
+ }
+ )
+
+ if args.plot :
+ pyplot.plot(t_vec, v_soma, 'k', linewidth=1, label='soma ' + str(nseg))
+ pyplot.plot(t_vec, v_dend, 'b', linewidth=1, label='dend ' + str(nseg))
+ pyplot.plot(t_vec, v_clamp, 'r', linewidth=1, label='clamp ' + str(nseg))
+
+# time the simulations
+end = timer()
+print "took ", end-start, " seconds"
+
+# save the spike info as in json format
+fp = open('ball_and_stick.json', 'w')
+json.dump(results, fp, indent=2)
+
+if args.plot :
+ pyplot.xlabel('time (ms)')
+ pyplot.ylabel('mV')
+ pyplot.xlim([0, simdur])
+ pyplot.legend()
+
+ pyplot.show()
+
diff --git a/nrn/generate_validation_data.sh b/nrn/generate_validation_data.sh
new file mode 100755
index 0000000000000000000000000000000000000000..8627a0c56ecd1f2137178a193be99207cb06ef46
--- /dev/null
+++ b/nrn/generate_validation_data.sh
@@ -0,0 +1,2 @@
+python2.7 soma.py
+
diff --git a/nrn/soma.py b/nrn/soma.py
new file mode 100644
index 0000000000000000000000000000000000000000..3af4f50f3e6e1ae7761591ab68fd32a1187f9c87
--- /dev/null
+++ b/nrn/soma.py
@@ -0,0 +1,74 @@
+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 info for a soma with hh channels')
+parser.add_argument('--plot', action='store_true', dest='plot')
+args = parser.parse_args()
+
+if args.plot :
+ print '-- plotting turned on'
+else :
+ print '-- plotting turned off'
+
+soma = h.Section(name='soma')
+
+soma.L = soma.diam = 18.8
+soma.Ra = 123
+
+print "Surface area of soma =", h.area(0.5, sec=soma)
+
+soma.insert('hh')
+
+stim = h.IClamp(soma(0.5))
+stim.delay = 10
+stim.dur = 100
+stim.amp = 0.1
+
+spike_counter = h.APCount(soma(0.5))
+spike_counter.thresh = 0
+
+v_vec = h.Vector() # Membrane potential vector
+t_vec = h.Vector() # Time stamp vector
+s_vec = h.Vector() # Time stamp vector
+v_vec.record(soma(0.5)._ref_v)
+t_vec.record(h._ref_t)
+spike_counter.record(s_vec)
+
+simdur = 120
+
+# initialize plot
+if args.plot :
+ pyplot.figure(figsize=(8,4)) # Default figsize is (8,6)
+
+h.tstop = simdur
+
+# run neuron with multiple dt
+start = timer()
+results = []
+for dt in [0.05, 0.02, 0.01, 0.001, 0.0001]:
+ h.dt = dt
+ h.run()
+ results.append({"dt": dt, "spikes": s_vec.to_python()})
+ if args.plot :
+ pyplot.plot(t_vec, v_vec, label='neuron ' + str(dt))
+
+end = timer()
+print "took ", end-start, " seconds"
+
+# save the spike info as in json format
+fp = open('soma.json', 'w')
+json.dump(results, fp, indent=1)
+
+if args.plot :
+ pyplot.xlabel('time (ms)')
+ pyplot.ylabel('mV')
+ pyplot.xlim([0, 120])
+ pyplot.grid()
+ pyplot.legend()
+ pyplot.show()
+
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 5696287e9b2c56096146c67fe0d891548a0abfd4..8a8286bda8e314d47f5e89b2b454bd6b84709e67 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,10 +1,11 @@
set(HEADERS
swcio.hpp
- )
+)
set(BASE_SOURCES
+ cell.cpp
+ mechanism_interface.cpp
parameter_list.cpp
swcio.cpp
- cell.cpp
)
add_library(cellalgo ${BASE_SOURCES} ${HEADERS})
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 8ac62055f351212fb75c72657fa6157463450f45..b3943c688b33b7a2114f2dc88d2f3d22905d273c 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -7,6 +7,8 @@
#include <type_traits>
#include <vector>
+#include "util.hpp"
+
/*
* Some simple wrappers around stl algorithms to improve readability of code
* that uses them.
@@ -130,15 +132,16 @@ namespace algorithms{
return false;
}
- std::vector<bool> is_leaf(parent_index.size(), false);
+ int n = parent_index.size();
+ std::vector<bool> is_leaf(n, false);
- for (std::size_t i = 1; i < parent_index.size(); ++i) {
+ for(auto i=1; i<n; ++i) {
auto p = parent_index[i];
- if (is_leaf[p]) {
+ if(is_leaf[p]) {
return false;
}
- if (p != i-1) {
+ if(p != i-1) {
// we have a branch and i-1 is a leaf node
is_leaf[i-1] = true;
}
@@ -163,7 +166,7 @@ namespace algorithms{
return count;
}
- template<typename C, bool CheckStrict = true>
+ template<typename C>
std::vector<typename C::value_type> branches(const C &parent_index)
{
static_assert(
@@ -171,11 +174,7 @@ namespace algorithms{
"integral type required"
);
- if (CheckStrict && !has_contiguous_segments(parent_index)) {
- throw std::invalid_argument(
- "parent_index has not contiguous branch numbering"
- );
- }
+ EXPECTS(has_contiguous_segments(parent_index));
auto num_child = child_count(parent_index);
std::vector<typename C::value_type> branch_runs(
@@ -196,9 +195,54 @@ namespace algorithms{
}
template<typename C>
- std::vector<typename C::value_type> branches_fast(const C &parent_index)
+ bool is_sorted(const C& c)
+ {
+ return std::is_sorted(c.begin(), c.end());
+ }
+
+ template<typename C>
+ bool is_unique(const C& c)
+ {
+ return std::adjacent_find(c.begin(), c.end()) == c.end();
+ }
+
+ /// Return and index that maps entries in sub to their corresponding
+ /// values in super, where sub is a subset of super.
+ ///
+ /// Both sets are sorted and have unique entries.
+ /// Complexity is O(n), where n is size of super
+ template<typename C>
+ // C::iterator models forward_iterator
+ // C::value_type is_integral
+ C index_into(const C& super, const C& sub)
{
- return branches<C,false>(parent_index);
+ //EXPECTS {s \in super : \forall s \in sub};
+ EXPECTS(is_unique(super) && is_unique(sub));
+ EXPECTS(is_sorted(super) && is_sorted(sub));
+ EXPECTS(sub.size() <= super.size());
+
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "index_into only applies to integral types"
+ );
+
+ C out(sub.size()); // out will have one entry for each index in sub
+
+ auto sub_it=sub.begin();
+ auto super_it=super.begin();
+ auto sub_idx=0u, super_idx = 0u;
+
+ while(sub_it!=sub.end() && super_it!=super.end()) {
+ if(*sub_it==*super_it) {
+ out[sub_idx] = super_idx;
+ ++sub_it; ++sub_idx;
+ }
+ ++super_it; ++super_idx;
+ }
+
+ EXPECTS(sub_idx==sub.size());
+
+ return out;
}
} // namespace algorithms
diff --git a/src/cell.cpp b/src/cell.cpp
index ab01007a206ea241befea21fbed83c531298a36d..a6986f0673ca05827ee89493aca2230ec16cec30 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -4,6 +4,20 @@
namespace nest {
namespace mc {
+int find_compartment_index(
+ segment_location const& location,
+ compartment_model const& graph
+) {
+ EXPECTS(location.segment<graph.segment_index.size());
+ const auto& si = graph.segment_index;
+ const auto seg = location.segment;
+
+ auto first = si[seg];
+ auto n = si[seg+1] - first;
+ auto index = std::floor(n*location.position);
+ return index<n ? first+index : first+n-1;
+}
+
cell::cell()
{
// insert a placeholder segment for the soma
@@ -20,7 +34,7 @@ int cell::num_segments() const
// note: I think that we have to enforce that the soma is the first
// segment that is added
//
-void cell::add_soma(value_type radius, point_type center)
+soma_segment* cell::add_soma(value_type radius, point_type center)
{
if(has_soma()) {
throw std::domain_error(
@@ -35,9 +49,11 @@ void cell::add_soma(value_type radius, point_type center)
else {
segments_[0] = make_segment<soma_segment>(radius);
}
+
+ return segments_[0]->as_soma();
}
-void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
+cable_segment* cell::add_cable(cell::index_type parent, segment_ptr&& cable)
{
// check for a valid parent id
if(cable->is_soma()) {
@@ -54,6 +70,8 @@ void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
}
segments_.push_back(std::move(cable));
parents_.push_back(parent);
+
+ return segments_.back()->as_cable();
}
segment* cell::segment(int index)
@@ -137,9 +155,9 @@ std::vector<int> cell::compartment_counts() const
return comp_count;
}
-int cell::num_compartments() const
+size_t cell::num_compartments() const
{
- auto n = 0;
+ auto n = 0u;
for(auto &s : segments_) {
n += s->num_compartments();
}
@@ -158,6 +176,20 @@ compartment_model cell::model() const
return m;
}
+
+void cell::add_stimulus( segment_location loc, i_clamp stim)
+{
+ if(!(loc.segment<num_segments())) {
+ throw std::out_of_range(
+ util::pprintf(
+ "can't insert stimulus in segment % of a cell with % segments",
+ loc.segment, num_segments()
+ )
+ );
+ }
+ stimulii_.push_back({loc, std::move(stim)});
+}
+
std::vector<int> const& cell::segment_parents() const
{
return parents_;
diff --git a/src/cell.hpp b/src/cell.hpp
index 68bc4c5a21f99237b18557f690226c0441e580c6..ba5f9b6bff227231cbb4f46e27dd2af4ff2943ce 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -5,8 +5,9 @@
#include <thread>
#include <vector>
-#include "segment.hpp"
-#include "cell_tree.hpp"
+#include <segment.hpp>
+#include <cell_tree.hpp>
+#include <stimulus.hpp>
namespace nest {
namespace mc {
@@ -19,6 +20,21 @@ struct compartment_model {
std::vector<int> segment_index;
};
+struct segment_location {
+ segment_location(int s, double l)
+ : segment(s), position(l)
+ {
+ EXPECTS(position>=0. && position<=1.);
+ }
+ int segment;
+ double position;
+};
+
+int find_compartment_index(
+ segment_location const& location,
+ compartment_model const& graph
+);
+
/// high-level abstract representation of a cell and its segments
class cell {
public:
@@ -33,18 +49,18 @@ class cell {
/// add a soma to the cell
/// radius must be specified
- void add_soma(value_type radius, point_type center=point_type());
+ soma_segment* add_soma(value_type radius, point_type center=point_type());
/// add a cable
/// parent is the index of the parent segment for the cable section
/// cable is the segment that will be moved into the cell
- void add_cable(index_type parent, segment_ptr&& cable);
+ cable_segment* add_cable(index_type parent, segment_ptr&& cable);
/// add a cable by constructing it in place
/// parent is the index of the parent segment for the cable section
/// args are the arguments to be used to consruct the new cable
template <typename... Args>
- void add_cable(index_type parent, Args ...args);
+ cable_segment* add_cable(index_type parent, Args ...args);
/// the number of segments in the cell
int num_segments() const;
@@ -70,7 +86,7 @@ class cell {
value_type area() const;
/// the total number of compartments over all segments
- int num_compartments() const;
+ size_t num_compartments() const;
std::vector<segment_ptr> const& segments() const;
@@ -83,17 +99,33 @@ class cell {
compartment_model model() const;
+ void add_stimulus(segment_location loc, i_clamp stim);
+
+ std::vector<std::pair<segment_location, i_clamp>>&
+ stimulii() {
+ return stimulii_;
+ }
+
+ const std::vector<std::pair<segment_location, i_clamp>>&
+ stimulii() const {
+ return stimulii_;
+ }
+
private:
// storage for connections
std::vector<index_type> parents_;
+
// the segments
std::vector<segment_ptr> segments_;
+
+ // the stimulii
+ std::vector<std::pair<segment_location, i_clamp>> stimulii_;
};
// create a cable by forwarding cable construction parameters provided by the user
template <typename... Args>
-void cell::add_cable(cell::index_type parent, Args ...args)
+cable_segment* cell::add_cable(cell::index_type parent, Args ...args)
{
// check for a valid parent id
if(parent>=num_segments()) {
@@ -103,6 +135,8 @@ void cell::add_cable(cell::index_type parent, Args ...args)
}
segments_.push_back(make_segment<cable_segment>(std::forward<Args>(args)...));
parents_.push_back(parent);
+
+ return segments_.back()->as_cable();
}
} // namespace mc
diff --git a/src/fvm.hpp b/src/fvm.hpp
index b0eedb412c94aec0e3a6e4b43aab86f14772f21a..69d23d0d1fe97e341ad8de3dc25ce8c0771ef78c 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -1,13 +1,20 @@
#pragma once
#include <algorithm>
+#include <map>
+#include <string>
+#include <vector>
+#include <algorithms.hpp>
#include <cell.hpp>
-#include <segment.hpp>
+#include <ion.hpp>
#include <math.hpp>
#include <matrix.hpp>
+#include <mechanism.hpp>
+#include <mechanism_interface.hpp>
#include <util.hpp>
-#include <algorithms.hpp>
+#include <segment.hpp>
+#include <stimulus.hpp>
#include <vector/include/Vector.hpp>
@@ -24,8 +31,16 @@ class fvm_cell {
/// the integral index type
using size_type = I;
+ /// the type used to store matrix information
using matrix_type = matrix<value_type, size_type>;
+ /// mechanism type
+ using mechanism_type =
+ nest::mc::mechanisms::mechanism_ptr<value_type, size_type>;
+
+ /// ion species storage
+ using ion_type = mechanisms::ion<value_type, size_type>;
+
/// the container used for indexes
using index_type = memory::HostVector<size_type>;
/// view into index container
@@ -44,28 +59,85 @@ class fvm_cell {
/// TODO this should be const
/// which requires const_view in the vector library
- matrix_type& jacobian();
+ 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();
+ vector_view cv_areas() {
+ return cv_areas_;
+ }
/// TODO this should be const
/// return the capacitance of each CV surface
- /// note that this is the total capacitance, not per unit area
- /// which is equivalent to sigma_i * c_m
- vector_view cv_capacitance();
+ /// this is the total capacitance, not per unit area,
+ /// i.e. equivalent to sigma_i * c_m
+ vector_view cv_capacitance() {
+ return cv_capacitance_;
+ }
- std::size_t size() const
- {
+ /// return the voltage in each CV
+ vector_view voltage() {
+ return voltage_;
+ }
+
+ std::size_t size() const {
return matrix_.size();
}
+ /// return reference to in iterable container of the mechanisms
+ std::vector<mechanism_type>& mechanisms() {
+ return mechanisms_;
+ }
+
+ /// return reference to list of ions
+ //std::map<mechanisms::ionKind, ion_type> ions_;
+ std::map<mechanisms::ionKind, ion_type>& ions() {
+ return ions_;
+ }
+ std::map<mechanisms::ionKind, ion_type> const& ions() const {
+ return ions_;
+ }
+
+ /// return reference to sodium ion
+ ion_type& ion_na() {
+ return ions_[mechanisms::ionKind::na];
+ }
+ ion_type const& ion_na() const {
+ return ions_[mechanisms::ionKind::na];
+ }
+
+ /// return reference to calcium ion
+ ion_type& ion_ca() {
+ return ions_[mechanisms::ionKind::ca];
+ }
+ ion_type const& ion_ca() const {
+ return ions_[mechanisms::ionKind::ca];
+ }
+
+ /// return reference to pottasium ion
+ ion_type& ion_k() {
+ return ions_[mechanisms::ionKind::k];
+ }
+ ion_type const& ion_k() const {
+ return ions_[mechanisms::ionKind::k];
+ }
+
+ /// make a time step
+ void advance(value_type dt);
+
+ /// set initial states
+ void initialize();
+
private:
+ /// current time
+ value_type t_ = value_type{0};
+
/// the linear system for implicit time stepping of cell state
matrix_type matrix_;
@@ -89,11 +161,18 @@ class fvm_cell {
/// the potential in mV in each CV
vector_type voltage_;
+ /// the set of mechanisms present in the cell
+ std::vector<mechanism_type> mechanisms_;
+
+ /// the ion species
+ std::map<mechanisms::ionKind, ion_type> ions_;
+
+ std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
};
-////////////////////////////////////////////////////////////
-// Implementation
-////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////// Implementation ////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
template <typename T, typename I>
fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
@@ -106,8 +185,16 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
using util::left;
using util::right;
+ // TODO: potential code stink
+ // matrix_ is a member, but it is not initialized with the other members
+ // above because it requires the parent_index, which is calculated
+ // "on the fly" by cell.model().
+ // cell.model() is quite expensive, and the information it calculates is
+ // used elsewhere, so we defer the intialization to inside the constructor
+ // body.
const auto graph = cell.model();
matrix_ = matrix_type(graph.parent_index);
+
auto parent_index = matrix_.p();
auto const& segment_index = graph.segment_index;
@@ -148,7 +235,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
auto radius_center = math::mean(c.radius);
auto area_face = math::area_circle( radius_center );
face_alpha_[i] = area_face / (c_m * r_L * c.length);
- cv_capacitance_[i] = c_m;
auto halflen = c.length/2;
@@ -170,27 +256,127 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
for(auto i=0u; i<size(); ++i) {
cv_capacitance_[i] /= cv_areas_[i];
}
-}
-template <typename T, typename I>
-typename fvm_cell<T,I>::matrix_type&
-fvm_cell<T,I>::jacobian()
-{
- return matrix_;
-}
+ /////////////////////////////////////////////
+ // create mechanisms
+ /////////////////////////////////////////////
-template <typename T, typename I>
-typename fvm_cell<T,I>::vector_view
-fvm_cell<T,I>::cv_areas()
-{
- return cv_areas_;
-}
+ // FIXME : candidate for a private member function
-template <typename T, typename I>
-typename fvm_cell<T,I>::vector_view
-fvm_cell<T,I>::cv_capacitance()
-{
- return cv_capacitance_;
+ // for each mechanism in the cell record the indexes of the segments that
+ // contain the mechanism
+ std::map<std::string, std::vector<int>> mech_map;
+
+ for(auto i=0; i<cell.num_segments(); ++i) {
+ for(const auto& mech : cell.segment(i)->mechanisms()) {
+ // FIXME : Membrane has to be a proper mechanism,
+ // because it is exposed via the public interface.
+ // This if statement is bad
+ if(mech.name() != "membrane") {
+ mech_map[mech.name()].push_back(i);
+ }
+ }
+ }
+
+ // Create the mechanism implementations with the state for each mechanism
+ // instance.
+ // TODO : this works well for density mechanisms (e.g. ion channels), but
+ // does it work for point processes (e.g. synapses)?
+ for(auto& mech : mech_map) {
+ auto& helper = nest::mc::mechanisms::get_mechanism_helper(mech.first);
+
+ // calculate the number of compartments that contain the mechanism
+ auto num_comp = 0u;
+ for(auto seg : mech.second) {
+ num_comp += segment_index[seg+1] - segment_index[seg];
+ }
+
+ // build a vector of the indexes of the compartments that contain
+ // the mechanism
+ std::vector<int> compartment_index(num_comp);
+ auto pos = 0u;
+ for(auto seg : mech.second) {
+ auto seg_size = segment_index[seg+1] - segment_index[seg];
+ std::iota(
+ compartment_index.data() + pos,
+ compartment_index.data() + pos + seg_size,
+ segment_index[seg]
+ );
+ pos += seg_size;
+ }
+
+ // instantiate the mechanism
+ index_view node_index(compartment_index.data(), compartment_index.size());
+ mechanisms_.push_back(
+ helper->new_mechanism(voltage_, current_, node_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
+ 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);
+ }
+ }
+
+ // create the ion state
+ if(auto n = std::distance(indexes.first.begin(), ends.first)) {
+ ions_.emplace(ion, index_view(indexes.first.data(), n));
+ }
+
+ // join the ion reference in each mechanism into the cell-wide ion state
+ for(auto& mech : mechanisms_) {
+ if(mech->uses_ion(ion)) {
+ mech->set_ion(ion, ions_[ion]);
+ }
+ }
+ }
+
+ // FIXME: Hard code parameters for now.
+ // Take defaults for reversal potential of sodium and potassium from
+ // the default values in Neuron.
+ // Neuron's defaults are defined in the file
+ // nrn/src/nrnoc/membdef.h
+ using memory::all;
+
+ constexpr value_type DEF_vrest = -65.0; // same name as #define in Neuron
+
+ ion_na().reversal_potential()(all) = 115+DEF_vrest; // mV
+ ion_na().internal_concentration()(all) = 10.0; // mM
+ ion_na().external_concentration()(all) = 140.0; // mM
+
+ ion_k().reversal_potential()(all) = -12.0+DEF_vrest;// mV
+ ion_k().internal_concentration()(all) = 54.4; // mM
+ ion_k().external_concentration()(all) = 2.5; // mM
+
+ ion_ca().reversal_potential()(all) = 12.5 * std::log(2.0/5e-5);// mV
+ ion_ca().internal_concentration()(all) = 5e-5; // mM
+ ion_ca().external_concentration()(all) = 2.0; // mM
+
+ // add the stimulii
+ for(const auto& stim : cell.stimulii()) {
+ auto idx = find_compartment_index(stim.first, graph);
+ stimulii_.push_back( {idx, stim.second} );
+ }
}
template <typename T, typename I>
@@ -217,9 +403,12 @@ void fvm_cell<T, I>::setup_matrix(T dt)
// . . .
// l[i] . . d[i]
//
+ //d(all) = 1.0;
d(all) = cv_areas_;
for(auto i=1u; i<d.size(); ++i) {
- auto a = dt * face_alpha_[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;
l[i] = -a;
@@ -228,14 +417,77 @@ 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
- // sigma_i * (V[i] - dt/cm*(im - ie))
+ // V[i] - dt/cm*(im - ie)
+ auto factor = 10.*dt;
for(auto i=0u; i<d.size(); ++i) {
- rhs[i] = cv_areas_[i] * (voltage_[i] - dt/cv_capacitance_[i]*current_[i]);
+ //rhs[i] = voltage_[i] - factor/cv_capacitance_[i]*current_[i];
+ rhs[i] = cv_areas_[i]*(voltage_[i] - factor/cv_capacitance_[i]*current_[i]);
+ }
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::initialize()
+{
+ t_ = 0.;
+
+ // initialize mechanism states
+ for(auto& m : mechanisms_) {
+ m->nrn_init();
+ }
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::advance(T dt)
+{
+ using memory::all;
+
+ current_(all) = 0.;
+
+ // update currents from ion channels
+ for(auto& m : mechanisms_) {
+ m->set_params(t_, dt);
+ m->nrn_current();
+ }
+
+ // add current contributions from stimulii
+ for(auto& stim : stimulii_) {
+ auto ie = stim.second.amplitude(t_);
+ auto loc = stim.first;
+
+ // the factor of 100 scales the injected current to 10^2.nA
+ 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";
}
} // namespace fvm
} // namespace mc
} // namespace nest
+
+
diff --git a/src/indexed_view.hpp b/src/indexed_view.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0cbc22a91b3235f70af6b851050af5e5fdefa3a2
--- /dev/null
+++ b/src/indexed_view.hpp
@@ -0,0 +1,50 @@
+#pragma once
+
+#include <vector/include/Vector.hpp>
+
+namespace nest {
+namespace mc {
+
+template <typename T, typename I>
+struct indexed_view {
+ using value_type = T;
+ using size_type = I;
+ using view_type = typename memory::HostVector<T>::view_type;
+ using index_view_type = typename memory::HostVector<I>::view_type;;
+
+ view_type view;
+ index_view_type index; // TODO make this a const view
+
+ indexed_view(view_type v, index_view_type i)
+ : view(v),
+ index(i)
+ {}
+
+ size_type size() const
+ {
+ return index.size();
+ }
+
+ // TODO
+ //
+ // these should either
+ // - return the internal reference type implementations from the containers
+ // e.g. the GPU reference that does a copy-on-read from GPU memory
+ // - or ensure that the result of dereferencing these is properly handled
+ // i.e. by just returning a value for the const version
+
+ value_type&
+ operator[] (const size_type i)
+ {
+ return view[index[i]];
+ }
+
+ value_type const&
+ operator[] (const size_type i) const
+ {
+ return view[index[i]];
+ }
+};
+
+} // namespace mc
+} // namespace nest
diff --git a/src/ion.hpp b/src/ion.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..311e86003b4cfba2b4cd4ceb569fc97a346daca0
--- /dev/null
+++ b/src/ion.hpp
@@ -0,0 +1,111 @@
+#pragma once
+
+#include <vector/include/Vector.hpp>
+#include "indexed_view.hpp"
+
+namespace nest {
+namespace mc {
+namespace mechanisms {
+
+/*******************************************************************************
+
+ ion channels have the following fields :
+
+ ---------------------------------------------------
+ label Ca Na K name
+ ---------------------------------------------------
+ iX ica ina ik current
+ eX eca ena ek reversal_potential
+ Xi cai nai ki internal_concentration
+ Xo cao nao ko external_concentration
+ gX gca gna gk conductance
+ ---------------------------------------------------
+*******************************************************************************/
+
+/// let's enumerate the ion channel types
+enum class ionKind {ca, na, k};
+
+[[gnu::unused]] static
+std::string to_string(ionKind k)
+{
+ switch(k) {
+ case ionKind::na : return "sodium";
+ case ionKind::ca : return "calcium";
+ case ionKind::k : return "pottasium";
+ }
+ return "unkown";
+}
+
+/// and a little helper to iterate over them
+[[gnu::unused]] static
+std::vector<ionKind> ion_kinds()
+{
+ return {ionKind::ca, ionKind::na, ionKind::k};
+}
+
+/// storage for ion channel information in a cell group
+template<typename T, typename I>
+class ion {
+public :
+ // expose tempalte parameters
+ using value_type = T;
+ using size_type = I;
+
+ // define storage types
+ using vector_type = memory::HostVector<value_type>;
+ using index_type = memory::HostVector<size_type>;
+ using vector_view_type = typename vector_type::view_type;
+ using index_view_type = typename index_type::view_type;
+
+ using indexed_view_type = indexed_view<value_type, size_type>;
+
+ ion() = default;
+
+ ion(index_view_type idx)
+ : node_index_{idx}
+ , iX_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()}
+ , eX_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()}
+ , Xi_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()}
+ , Xo_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()}
+ { }
+
+ std::size_t memory() const {
+ return 4u*size() * sizeof(value_type)
+ + size() * sizeof(index_type)
+ + sizeof(ion);
+ }
+
+ vector_view_type current() {
+ return iX_;
+ }
+ vector_view_type reversal_potential() {
+ return eX_;
+ }
+ vector_view_type internal_concentration() {
+ return Xi_;
+ }
+ vector_view_type external_concentration() {
+ return Xo_;
+ }
+
+ index_type node_index() {
+ return node_index_;
+ }
+
+ size_t size() const {
+ return node_index_.size();
+ }
+
+private :
+
+ index_type node_index_;
+ vector_type iX_;
+ vector_type eX_;
+ vector_type Xi_;
+ vector_type Xo_;
+};
+
+} // namespace mechanisms
+} // namespace mc
+} // namespace nest
+
diff --git a/src/matrix.hpp b/src/matrix.hpp
index 16936fa47c26d3c632d08d008e59c76c2e4b0db0..db9d2df7f56e8bfffb73051e996acb13d9fa0369 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -19,10 +19,10 @@ class matrix {
using size_type = I;
// define storage types
- using vector_type = memory::HostVector<value_type>;
- using view_type = typename vector_type::view_type;
- using index_type = memory::HostVector<size_type>;
- using index_view = typename index_type::view_type;
+ 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;
matrix() = default;
@@ -67,7 +67,7 @@ class matrix {
/// the total memory used to store the matrix
std::size_t memory() const
{
- auto s = 5 * (sizeof(value_type) * size() + sizeof(vector_type));
+ auto s = 6 * (sizeof(value_type) * size() + sizeof(vector_type));
s += sizeof(size_type) * (parent_index_.size() + cell_index_.size())
+ 2*sizeof(index_type);
s += sizeof(matrix);
@@ -80,32 +80,45 @@ class matrix {
return cell_index_.size() - 1;
}
+ /// FIXME : make modparser use the correct accessors (l,d,u,rhs) instead of these
+ vector_view_type vec_rhs() { return rhs(); }
+ vector_view_type vec_d() { return d(); }
+ vector_view_type vec_v() { return v(); }
+
/// the vector holding the lower part of the matrix
- view_type l()
+ vector_view_type l()
{
return l_;
}
/// the vector holding the diagonal of the matrix
- view_type d()
+ vector_view_type d()
{
return d_;
}
/// the vector holding the upper part of the matrix
- view_type u()
+ vector_view_type u()
{
return u_;
}
/// the vector holding the right hand side of the linear equation system
- view_type rhs()
+ vector_view_type rhs()
{
return rhs_;
}
+ /// the vector holding the solution (voltage)
+ vector_view_type v()
+ {
+ EXPECTS(has_voltage_);
+
+ return v_;
+ }
+
/// the vector holding the parent index
- index_view p()
+ index_view_type p()
{
return parent_index_;
}
@@ -115,7 +128,16 @@ class matrix {
/// and can be accessed via rhs()
void solve()
{
- index_view const& p = parent_index_;
+ /*
+ std::cout << "solving matrix :\n";
+ std::cout << " l " << l_ << "\n";
+ std::cout << " d " << d_ << "\n";
+ std::cout << " u " << u_ << "\n";
+ std::cout << " rhs " << rhs_ << "\n";
+ */
+
+ // FIXME make a const view
+ index_view_type const& p = parent_index_;
auto const ncells = num_cells();
// loop over submatrices
@@ -137,6 +159,17 @@ class matrix {
rhs_[i] /= d_[i];
}
}
+ //std::cout << " v " << rhs_ << "\n";
+ }
+
+ void add_voltage(vector_view_type v_ext)
+ {
+ EXPECTS(v_ext.size()==size());
+
+ std::cout << "============ adding voltage" << std::endl;
+
+ v_ = v_ext;
+ has_voltage_ = true;
}
private:
@@ -164,8 +197,12 @@ class matrix {
vector_type l_;
vector_type d_;
vector_type u_;
+
/// after calling solve, the solution is stored in rhs_
vector_type rhs_;
+ vector_view_type v_;
+
+ bool has_voltage_=false;
};
} // namespace nest
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
index 3432225885848ed06f1f97513772709025757a20..ea69c292e1f5c50efe0fd0e6e7dbf092c32b16fb 100644
--- a/src/mechanism.hpp
+++ b/src/mechanism.hpp
@@ -5,17 +5,22 @@
#include <memory>
#include <string>
-#include "matrix.hpp"
+#include "indexed_view.hpp"
#include "parameter_list.hpp"
#include "util.hpp"
+#include "ion.hpp"
namespace nest {
namespace mc {
+namespace mechanisms {
enum class mechanismKind {point, density};
template <typename T, typename I>
class mechanism {
+
+public:
+
using value_type = T;
using size_type = I;
@@ -24,31 +29,50 @@ class mechanism {
using view_type = typename vector_type::view_type;
using index_type = memory::HostVector<size_type>;
using index_view = typename index_type::view_type;
- //using indexed_view= IndexedView<value_type, size_type, target()>;
+ using indexed_view_type = indexed_view<value_type, size_type>;
- using matrix_type = matrix<value_type, size_type>;
+ using ion_type = ion<value_type, size_type>;
- mechanism(matrix_type *matrix, index_view node_indices)
- : matrix_(matrix)
- , node_indices_(node_indices)
+ mechanism(view_type vec_v, view_type vec_i, index_view node_index)
+ : vec_v_(vec_v)
+ , vec_i_(vec_i)
+ , node_index_(node_index)
{}
std::size_t size() const
{
- return node_indices_.size();
+ return node_index_.size();
+ }
+
+ index_view node_index() const
+ {
+ return node_index_;
+ }
+
+ value_type voltage(size_type i) const
+ {
+ return vec_v_[node_index_[i]];
+ }
+
+ value_type current(size_type i) const
+ {
+ return vec_i_[node_index_[i]];
}
virtual void set_params(value_type t_, value_type dt_) = 0;
virtual std::string name() const = 0;
virtual std::size_t memory() const = 0;
- virtual void init() = 0;
- virtual void state() = 0;
- virtual void current() = 0;
+ virtual void nrn_init() = 0;
+ virtual void nrn_state() = 0;
+ virtual void nrn_current() = 0;
+ virtual bool uses_ion(ionKind) const = 0;
+ virtual void set_ion(ionKind k, ion_type& i) = 0;
virtual mechanismKind kind() const = 0;
- matrix_type* matrix_;
- index_type node_indices_;
+ view_type vec_v_;
+ view_type vec_i_;
+ index_type node_index_;
};
template <typename T, typename I>
@@ -57,11 +81,13 @@ using mechanism_ptr = std::unique_ptr<mechanism<T,I>>;
template <typename M>
mechanism_ptr<typename M::value_type, typename M::size_type>
make_mechanism(
- typename M::matrix_type* matrix,
- typename M::index_view node_indices
+ typename M::view_type vec_v,
+ typename M::view_type vec_i,
+ typename M::index_view node_indices
) {
- return util::make_unique<M>(matrix, node_indices);
+ return util::make_unique<M>(vec_v, vec_i, node_indices);
}
+} // namespace mechanisms
} // namespace nest
} // namespace mc
diff --git a/src/mechanism_interface.cpp b/src/mechanism_interface.cpp
index ba24eb06ce13c3452bb5db6d8847e3b1b1ce3081..10fad45fa0d9776c4f45302e6bf196f995765e99 100644
--- a/src/mechanism_interface.cpp
+++ b/src/mechanism_interface.cpp
@@ -1,21 +1,41 @@
#include "mechanism_interface.hpp"
-/* include the mechanisms
+//
+// include the mechanisms
+//
+
#include <mechanisms/hh.hpp>
#include <mechanisms/pas.hpp>
-*/
namespace nest {
namespace mc {
namespace mechanisms {
-std::map<std::string, mechanism_helper<double, int>> mechanism_map;
+std::map<std::string, mechanism_helper_ptr<value_type, index_type>> mechanism_helpers;
void setup_mechanism_helpers() {
- /* manually insert
- mechanism_map["hh"] = mechanisms::hh;
- mechanism_map["pas"] = mechanisms::pas;
- */
+ mechanism_helpers["pas"] =
+ make_mechanism_helper<
+ mechanisms::pas::helper<value_type, index_type>
+ >();
+
+ mechanism_helpers["hh"] =
+ make_mechanism_helper<
+ mechanisms::hh::helper<value_type, index_type>
+ >();
+}
+
+mechanism_helper_ptr<value_type, index_type>&
+get_mechanism_helper(const std::string& name)
+{
+ auto helper = mechanism_helpers.find(name);
+ if(helper==mechanism_helpers.end()) {
+ throw std::out_of_range(
+ nest::mc::util::pprintf("there is no mechanism named \'%\'", name)
+ );
+ }
+
+ return helper->second;
}
} // namespace mechanisms
diff --git a/src/mechanism_interface.hpp b/src/mechanism_interface.hpp
index 389ef05997b90c96782dff8593aa5817ca4b2fde..c59caaad6c655d6a69c5206b956b1cb42959cfdf 100644
--- a/src/mechanism_interface.hpp
+++ b/src/mechanism_interface.hpp
@@ -3,30 +3,56 @@
#include <map>
#include <string>
-#include "matrix.hpp"
#include "mechanism.hpp"
#include "parameter_list.hpp"
namespace nest {
namespace mc {
+namespace mechanisms {
-///
+using value_type = double;
+using index_type = int;
+
+/// helper type for building mechanisms
+/// the use of abstract base classes everywhere is a bit ugly
template <typename T, typename I>
struct mechanism_helper {
- using index_type = memory::HostView<I>;
+ using value_type = T;
+ using size_type = I;
+ using index_type = memory::HostVector<I>;
using index_view = typename index_type::view_type;
- using mechanism_type = mechanism<T, I>;
- using matrix_type = typename mechanism_type::matrix_type;
+ using mechanism_ptr_type = mechanism_ptr<T, I>;
+ using view_type = typename mechanism<T,I>::view_type;
virtual std::string name() const = 0;
- virtual mechanism_type new_mechanism(matrix_type*, index_view) const = 0;
+ virtual mechanism_ptr<T,I> new_mechanism(view_type, view_type, index_view) const = 0;
- virtual void set_parameters(mechanism_type&, parameter_list const&) const = 0;
+ virtual void set_parameters(mechanism_ptr_type&, parameter_list const&) const = 0;
};
-extern std::map<std::string, mechanism_helper<double, int>> mechanism_helpers;
+template <typename T, typename I>
+using mechanism_helper_ptr =
+ std::unique_ptr<mechanism_helper<T,I>>;
+
+template <typename M>
+mechanism_helper_ptr<typename M::value_type, typename M::size_type>
+make_mechanism_helper()
+{
+ return util::make_unique<M>();
+}
+
+// for now use a global variable for the map of mechanism helpers
+extern std::map<
+ std::string,
+ mechanism_helper_ptr<value_type, index_type>
+> mechanism_helpers;
+
void setup_mechanism_helpers();
-} // namespace nest
+mechanism_helper_ptr<value_type, index_type>&
+get_mechanism_helper(const std::string& name);
+
+} // namespace mechanisms
} // namespace mc
+} // namespace nest
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index 56848466f58f36171e44b2d6da2c371a4edd84a9..0c06988756a077fdbb860fc8f42043b178b1b7e7 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -146,8 +146,15 @@ namespace mc {
membrane_parameters()
: base("membrane")
{
- base::add_parameter({"r_L", 0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m2
- base::add_parameter({"c_m", 180.00, {0., 1e9}}); // Ohm.cm
+ // from nrn/src/nrnoc/membdef.h
+ //#define DEF_cm 1. // uF/cm^2
+ //#define DEF_Ra 35.4 // ohm-cm
+ //#define DEF_celsius 6.3 // deg-C
+ //#define DEF_vrest -65. // mV
+ // r_L is called Ra in Neuron
+ //base::add_parameter({"c_m", 10e-6, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m^2 == 10^-6 F/cm^2
+ base::add_parameter({"c_m", 0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m^2 == 10^-6 F/cm^2
+ base::add_parameter({"r_L", 180.00, {0., 1e9}}); // equivalent to Ra in Neuron : Ohm.cm
}
};
@@ -169,10 +176,33 @@ namespace mc {
hh_parameters()
: base("hh")
{
- base::add_parameter({"gnabar", 0.12, {0., 1e9}});
- base::add_parameter({"gkbar", 0.036, {0., 1e9}});
- base::add_parameter({"gl", 0.0003,{0., 1e9}});
- base::add_parameter({"el", -54.3});
+ base::add_parameter({"gnabar", 0.12, {0, 1e9}});
+ base::add_parameter({"gkbar", 0.036, {0, 1e9}});
+ base::add_parameter({"gl", 0.0003,{0, 1e9}});
+ base::add_parameter({"el", -54.3});
+ }
+ };
+
+ /// parameters for passive channel
+ class pas_parameters
+ : public parameter_list
+ {
+ public:
+
+ using base = parameter_list;
+
+ using base::value_type;
+
+ using base::set;
+ using base::get;
+ using base::parameters;
+ using base::has_parameter;
+
+ pas_parameters()
+ : base("pas")
+ {
+ base::add_parameter({"g", 0.001, {0, 1e9}});
+ base::add_parameter({"e", -70});
}
};
diff --git a/src/segment.hpp b/src/segment.hpp
index e339e3ac4bc974477a624fe698d11bf8eb521f26..fd1d8ded6101cc8e9834b9062df3da97de5cbc8c 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -113,6 +113,31 @@ class segment {
return *it;
}
+ const parameter_list& mechanism(std::string const& n) const
+ {
+ auto it = std::find_if(
+ mechanisms_.begin(), mechanisms_.end(),
+ [&n](parameter_list const& l){return l.name() == n;}
+ );
+ if(it==mechanisms_.end()) {
+ throw std::out_of_range(
+ "attempt to access a parameter that is not defined in a segment"
+ );
+ }
+
+ return *it;
+ }
+
+ std::vector<parameter_list>& mechanisms()
+ {
+ return mechanisms_;
+ }
+
+ const std::vector<parameter_list>& mechanisms() const
+ {
+ return mechanisms_;
+ }
+
protected:
segmentKind kind_;
diff --git a/src/stimulus.hpp b/src/stimulus.hpp
index 94a4dddfddf1cab80272de74bd3a94f9130989fd..f7c7538cf398a7154588126620d97b7c054953a4 100644
--- a/src/stimulus.hpp
+++ b/src/stimulus.hpp
@@ -1,6 +1,7 @@
#pragma once
-namespace nestmc {
+namespace nest {
+namespace mc {
class i_clamp {
public:
@@ -49,4 +50,5 @@ class i_clamp {
value_type amplitude_ = 0;
};
-} // namespace nestmc
+} // namespace mc
+} // namespace nest
diff --git a/src/util.hpp b/src/util.hpp
index d92392799cf98da27ee1ca62788686a308c02c1d..706ab7da5a17d493823b8c74d1f18fa6b4a9efaf 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -2,6 +2,11 @@
#include "vector/include/Vector.hpp"
+#ifdef DEBUG
+#define EXPECTS(expression) assert(expression)
+#else
+#define EXPECTS(expression)
+#endif
/*
using memory::util::red;
@@ -23,7 +28,7 @@ operator << (std::ostream &o, std::vector<T>const& v)
{
o << "[";
for(auto const& i: v) {
- o << i << ", ";
+ o << i << ",";
}
o << "]";
return o;
@@ -34,7 +39,7 @@ std::ostream& print(std::ostream &o, std::vector<T>const& v)
{
o << "[";
for(auto const& i: v) {
- o << i << ", ";
+ o << i << ",";
}
o << "]";
return o;
@@ -80,6 +85,45 @@ namespace util {
std::false_type
>::type
{};
+
+ // printf with variadic templates for simplified string creation
+ // mostly used to simplify error string creation
+ [[gnu::unused]] static
+ std::string pprintf(const char *s) {
+ std::string errstring;
+ while(*s) {
+ if(*s == '%' && s[1]!='%') {
+ // instead of throwing an exception, replace with ??
+ //throw std::runtime_error("pprintf: the number of arguments did not match the format ");
+ errstring += "<?>";
+ }
+ else {
+ errstring += *s;
+ }
+ ++s;
+ }
+ return errstring;
+ }
+
+ template <typename T, typename ... Args>
+ std::string pprintf(const char *s, T value, Args... args) {
+ std::string errstring;
+ while(*s) {
+ if(*s == '%' && s[1]!='%') {
+ std::stringstream str;
+ str << value;
+ errstring += str.str();
+ errstring += pprintf(++s, args...);
+ return errstring;
+ }
+ else {
+ errstring += *s;
+ ++s;
+ }
+ }
+ return errstring;
+ }
+
} // namespace util
} // namespace mc
} // namespace nest
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b555fc465f3e3292ae2fa48f11dfad17bbe0e6e5..cb33af4d972dd398952b1879cca764c502d8f42c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -8,31 +8,49 @@ set(HEADERS
../src/tree.hpp
)
-set(TEST_SOURCES
- # google test framework
- gtest-all.cpp
+# google test framework
+add_library(gtest gtest-all.cpp gtest.h)
+set(TEST_SOURCES
# unit tests
test_algorithms.cpp
test_cell.cpp
test_compartments.cpp
+ test_fvm.cpp
test_matrix.cpp
+ test_mechanisms.cpp
+ test_parameters.cpp
test_point.cpp
- test_run.cpp
test_segment.cpp
test_stimulus.cpp
test_swcio.cpp
test_tree.cpp
# unit test driver
- main.cpp
+ test.cpp
+)
+
+set(VALIDATION_SOURCES
+ # unit tests
+ validate_soma.cpp
+ validate_ball_and_stick.cpp
+
+ # unit test driver
+ validate.cpp
)
add_executable(test.exe ${TEST_SOURCES} ${HEADERS})
+add_executable(validate.exe ${VALIDATION_SOURCES} ${HEADERS})
-target_link_libraries(test.exe LINK_PUBLIC cellalgo)
+target_link_libraries(test.exe LINK_PUBLIC cellalgo gtest)
+target_link_libraries(validate.exe LINK_PUBLIC cellalgo gtest)
set_target_properties(test.exe
PROPERTIES
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
)
+
+set_target_properties(validate.exe
+ PROPERTIES
+ RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
+)
diff --git a/tests/main.cpp b/tests/test.cpp
similarity index 100%
rename from tests/main.cpp
rename to tests/test.cpp
diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp
index 6503c07d9a413bf2e90f8cbf151b0df1e3e5949b..f69569c9611068250cc3265524cfe64e1e10bd55 100644
--- a/tests/test_algorithms.cpp
+++ b/tests/test_algorithms.cpp
@@ -255,6 +255,84 @@ TEST(algorithms, has_contiguous_segments)
);
}
+TEST(algorithms, is_unique)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_unique(
+ std::vector<int>{}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_unique(
+ std::vector<int>{0}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_unique(
+ std::vector<int>{0,1,100}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_unique(
+ std::vector<int>{0,0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_unique(
+ std::vector<int>{0,1,2,2,3,4}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_unique(
+ std::vector<int>{0,1,2,3,4,4}
+ )
+ );
+}
+
+TEST(algorithms, is_sorted)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{100}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{0,1,2}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{0,2,100}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{0,0}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{0,1,2,2,2,2,3,4,5,5,5}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{0,1,2,1}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_sorted(
+ std::vector<int>{1,0}
+ )
+ );
+}
+
TEST(algorithms, child_count)
{
{
@@ -310,7 +388,7 @@ TEST(algorithms, branches)
std::vector<int> expected_branches =
{ 0, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5 };
- auto actual_branches = algorithms::branches_fast(parent_index);
+ auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
}
@@ -329,7 +407,7 @@ TEST(algorithms, branches)
std::vector<int> expected_branches =
{ 0, 1, 1, 1 };
- auto actual_branches = algorithms::branches_fast(parent_index);
+ auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
}
@@ -350,7 +428,7 @@ TEST(algorithms, branches)
std::vector<int> expected_branches =
{ 0, 1, 1, 2, 3, 3 };
- auto actual_branches = algorithms::branches_fast(parent_index);
+ auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
}
@@ -358,7 +436,47 @@ TEST(algorithms, branches)
std::vector<int> parent_index = { 0 };
std::vector<int> expected_branches = { 0 };
- auto actual_branches = algorithms::branches_fast(parent_index);
+ auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
}
}
+
+TEST(algorithms, index_into)
+{
+ using C = std::vector<int>;
+
+ // by default index_into assumes that the inputs satisfy
+ // quite a strong set of prerequisites
+ //
+ // TODO: test that the EXPECTS() catch bad inputs when DEBUG mode is enabled
+ // put this in a seperate unit test
+ auto tests = {
+ std::make_pair(C{}, C{}),
+ std::make_pair(C{100}, C{}),
+ std::make_pair(C{0,1,3,4,6,7,10,11}, C{0,4,6,7,11}),
+ std::make_pair(C{0,1,3,4,6,7,10,11}, C{0}),
+ std::make_pair(C{0,1,3,4,6,7,10,11}, C{11}),
+ std::make_pair(C{0,1,3,4,6,7,10,11}, C{4}),
+ std::make_pair(C{0,1,3,4,6,7,10,11}, C{0,11}),
+ std::make_pair(C{0,1,3,4,6,7,10,11}, C{4,11}),
+ std::make_pair(C{0,1,3,4,6,7,10,11}, C{0,1,3,4,6,7,10,11})
+ };
+
+ auto test_result = [] (const C& super, const C& sub, const C& index) {
+ if(sub.size()!=index.size()) return false;
+ for(auto i=0u; i<sub.size(); ++i) {
+ if(index[i]>=C::value_type(super.size())) return false;
+ if(super[index[i]]!=sub[i]) return false;
+ }
+ return true;
+ };
+
+ for(auto& t : tests) {
+ EXPECT_TRUE(
+ test_result(
+ t.first, t.second,
+ nest::mc::algorithms::index_into(t.first, t.second)
+ )
+ );
+ }
+}
diff --git a/tests/test_fvm.cpp b/tests/test_fvm.cpp
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..85a25f3579aa25e2e84860b3144294269ff0895a 100644
--- a/tests/test_fvm.cpp
+++ b/tests/test_fvm.cpp
@@ -0,0 +1,105 @@
+#include <fstream>
+
+#include "gtest.h"
+#include "util.hpp"
+
+#include "../src/cell.hpp"
+#include "../src/fvm.hpp"
+
+TEST(fvm, cable)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ cell.add_soma(6e-4); // 6um in cm
+
+ // 1um radius and 4mm long, all in cm
+ cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1);
+ cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1);
+
+ //std::cout << cell.segment(1)->area() << " is the area\n";
+ EXPECT_EQ(cell.model().tree.num_segments(), 3u);
+
+ // add passive to all 3 segments in the cell
+ for(auto& seg :cell.segments()) {
+ seg->add_mechanism(pas_parameters());
+ }
+
+ cell.soma()->add_mechanism(hh_parameters());
+ cell.segment(2)->add_mechanism(hh_parameters());
+
+ auto& soma_hh = cell.soma()->mechanism("hh");
+
+ soma_hh.set("gnabar", 0.12);
+ soma_hh.set("gkbar", 0.036);
+ soma_hh.set("gl", 0.0003);
+ soma_hh.set("el", -54.387);
+
+ cell.segment(1)->set_compartments(4);
+ cell.segment(2)->set_compartments(4);
+
+ using fvm_cell = fvm::fvm_cell<double, int>;
+ fvm_cell fvcell(cell);
+ auto& J = fvcell.jacobian();
+
+ EXPECT_EQ(cell.num_compartments(), 9u);
+
+ // assert that the matrix has one row for each compartment
+ EXPECT_EQ(J.size(), cell.num_compartments());
+
+ fvcell.setup_matrix(0.02);
+
+ // assert that the number of cv areas is the same as the matrix size
+ // i.e. both should equal the number of compartments
+ EXPECT_EQ(fvcell.cv_areas().size(), J.size());
+}
+
+TEST(fvm, init)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ cell.add_soma(12.6157/2.0);
+ //auto& props = cell.soma()->properties;
+
+ cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
+
+ const auto m = cell.model();
+ EXPECT_EQ(m.tree.num_segments(), 2u);
+
+ // in this context (i.e. attached to a segment on a high-level cell)
+ // a mechanism is essentially a set of parameters
+ // - the only "state" is that used to define parameters
+ cell.soma()->add_mechanism(hh_parameters());
+
+ auto& soma_hh = cell.soma()->mechanism("hh");
+
+ soma_hh.set("gnabar", 0.12);
+ soma_hh.set("gkbar", 0.036);
+ soma_hh.set("gl", 0.0003);
+ soma_hh.set("el", -54.3);
+
+ // check that parameter values were set correctly
+ EXPECT_EQ(cell.soma()->mechanism("hh").get("gnabar").value, 0.12);
+ EXPECT_EQ(cell.soma()->mechanism("hh").get("gkbar").value, 0.036);
+ EXPECT_EQ(cell.soma()->mechanism("hh").get("gl").value, 0.0003);
+ EXPECT_EQ(cell.soma()->mechanism("hh").get("el").value, -54.3);
+
+ cell.segment(1)->set_compartments(10);
+
+ using fvm_cell = fvm::fvm_cell<double, int>;
+ fvm_cell fvcell(cell);
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 11u);
+
+ fvcell.setup_matrix(0.01);
+}
+
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
index f0271bd40b6ee109489dfd1f1edb781dca49cfa3..73663d6815d876e1740da4cee084c46255a74e3f 100644
--- a/tests/test_matrix.cpp
+++ b/tests/test_matrix.cpp
@@ -116,7 +116,6 @@ TEST(matrix, solve)
EXPECT_NEAR(0., std::sqrt(err), 1e-8);
}
-
}
}
diff --git a/tests/test_mechanisms.cpp b/tests/test_mechanisms.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d61c89d20b40136044dfd7556157b68947d8e56f
--- /dev/null
+++ b/tests/test_mechanisms.cpp
@@ -0,0 +1,37 @@
+#include "gtest.h"
+
+#include <mechanism_interface.hpp>
+#include <matrix.hpp>
+
+TEST(mechanisms, helpers) {
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ EXPECT_EQ(nest::mc::mechanisms::mechanism_helpers.size(), 2u);
+
+ // verify that the hh and pas channels are available
+ EXPECT_EQ(nest::mc::mechanisms::get_mechanism_helper("hh")->name(), "hh");
+ EXPECT_EQ(nest::mc::mechanisms::get_mechanism_helper("pas")->name(), "pas");
+
+ // check that an out_of_range exception is thrown if an invalid mechanism is
+ // requested
+ ASSERT_THROW(
+ nest::mc::mechanisms::get_mechanism_helper("dachshund"),
+ std::out_of_range
+ );
+
+ //0 1 2 3 4 5 6 7 8 9
+ std::vector<int> parent_index = {0,0,1,2,3,4,0,6,7,8};
+ memory::HostVector<int> node_indices = std::vector<int>{0,6,7,8,9};
+ auto n = node_indices.size();
+
+ //nest::mc::matrix<double, int> matrix(parent_index);
+ memory::HostVector<double> vec_i(n, 0.);
+ memory::HostVector<double> vec_v(n, 0.);
+
+ auto& helper = nest::mc::mechanisms::get_mechanism_helper("hh");
+ auto mech = helper->new_mechanism(vec_v, vec_i, node_indices);
+
+ EXPECT_EQ(mech->name(), "hh");
+ EXPECT_EQ(mech->size(), 5u);
+ //EXPECT_EQ(mech->matrix_, &matrix);
+}
diff --git a/tests/test_parameters.cpp b/tests/test_parameters.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0dce58b0c5472de6db9d721ad78dece7a9b3de94
--- /dev/null
+++ b/tests/test_parameters.cpp
@@ -0,0 +1,42 @@
+
+#include <fstream>
+
+#include "gtest.h"
+#include "util.hpp"
+
+#include <parameter_list.hpp>
+
+// test out the parameter infrastructure
+TEST(parameters, setting)
+{
+ nest::mc::parameter_list list("test");
+ EXPECT_EQ(list.name(), "test");
+ EXPECT_EQ(list.num_parameters(), 0);
+
+ nest::mc::parameter p("a", 0.12, {0, 10});
+
+ // add_parameter() returns a bool that indicates whether
+ // it was able to successfull add the parameter
+ EXPECT_TRUE(list.add_parameter(std::move(p)));
+ EXPECT_EQ(list.num_parameters(), 1);
+
+ // test in place construction of a parameter
+ EXPECT_TRUE(list.add_parameter({"b", -3.0}));
+ EXPECT_EQ(list.num_parameters(), 2);
+
+ // check that adding a parameter that already exists returns false
+ // and does not increase the number of parameters
+ EXPECT_FALSE(list.add_parameter({"b", -3.0}));
+ EXPECT_EQ(list.num_parameters(), 2);
+
+ auto &parms = list.parameters();
+ EXPECT_EQ(parms[0].name, "a");
+ EXPECT_EQ(parms[0].value, 0.12);
+ EXPECT_EQ(parms[0].range.min, 0);
+ EXPECT_EQ(parms[0].range.max, 10);
+
+ EXPECT_EQ(parms[1].name, "b");
+ EXPECT_EQ(parms[1].value, -3);
+ EXPECT_FALSE(parms[1].range.has_lower_bound());
+ EXPECT_FALSE(parms[1].range.has_upper_bound());
+}
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
deleted file mode 100644
index c4e3b1025e82ca257295fdcf8980fd82240d1956..0000000000000000000000000000000000000000
--- a/tests/test_run.cpp
+++ /dev/null
@@ -1,142 +0,0 @@
-#include "gtest.h"
-
-#include "../src/cell.hpp"
-#include "../src/fvm.hpp"
-
-TEST(run, cable)
-{
- using namespace nest::mc;
-
- nest::mc::cell cell;
-
- cell.add_soma(6e-4); // 6um in cm
-
- // 1um radius and 4mm long, all in cm
- cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1);
-
- std::cout << cell.segment(1)->area() << " is the area\n";
- EXPECT_EQ(cell.model().tree.num_segments(), 2u);
-
- cell.soma()->add_mechanism(hh_parameters());
-
- auto& soma_hh = cell.soma()->mechanism("hh");
-
- soma_hh.set("gnabar", 0.12);
- soma_hh.set("gkbar", 0.036);
- soma_hh.set("gl", 0.0003);
- soma_hh.set("el", -54.387);
-
- cell.segment(1)->set_compartments(4);
-
- using fvm_cell = fvm::fvm_cell<double, int>;
- fvm_cell fvcell(cell);
- auto& J = fvcell.jacobian();
- EXPECT_EQ(J.size(), 5u);
-
- fvcell.setup_matrix(0.02);
- EXPECT_EQ(fvcell.cv_areas().size(), J.size());
-
- auto& cable_parms = cell.segment(1)->mechanism("membrane");
- std::cout << soma_hh << std::endl;
- std::cout << cable_parms << std::endl;
-
- std::cout << "l " << J.l() << "\n";
- std::cout << "d " << J.d() << "\n";
- std::cout << "u " << J.u() << "\n";
- std::cout << "p " << J.p() << "\n";
-
- J.rhs()(memory::all) = 1.;
- J.rhs()[0] = 10.;
-
- J.solve();
-}
-
-TEST(run, init)
-{
- using namespace nest::mc;
-
- nest::mc::cell cell;
-
- cell.add_soma(12.6157/2.0);
- //auto& props = cell.soma()->properties;
-
- cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
-
- const auto m = cell.model();
- EXPECT_EQ(m.tree.num_segments(), 2u);
-
- // in this context (i.e. attached to a segment on a high-level cell)
- // a mechanism is essentially a set of parameters
- // - the only "state" is that used to define parameters
- cell.soma()->add_mechanism(hh_parameters());
-
- auto& soma_hh = cell.soma()->mechanism("hh");
-
- soma_hh.set("gnabar", 0.12);
- soma_hh.set("gkbar", 0.036);
- soma_hh.set("gl", 0.0003);
- soma_hh.set("el", -54.3);
-
- // check that parameter values were set correctly
- EXPECT_EQ(cell.soma()->mechanism("hh").get("gnabar").value, 0.12);
- EXPECT_EQ(cell.soma()->mechanism("hh").get("gkbar").value, 0.036);
- EXPECT_EQ(cell.soma()->mechanism("hh").get("gl").value, 0.0003);
- EXPECT_EQ(cell.soma()->mechanism("hh").get("el").value, -54.3);
-
- cell.segment(1)->set_compartments(10);
-
- using fvm_cell = fvm::fvm_cell<double, int>;
- fvm_cell fvcell(cell);
- auto& J = fvcell.jacobian();
- EXPECT_EQ(J.size(), 11u);
-
- fvcell.setup_matrix(0.01);
- std::cout << "areas " << fvcell.cv_areas() << "\n";
-
- //std::cout << "l" << J.l() << "\n";
- //std::cout << "d" << J.d() << "\n";
- //std::cout << "u" << J.u() << "\n";
-
- J.rhs()(memory::all) = 1.;
- J.rhs()[0] = 10.;
-
- J.solve();
-
- //std::cout << "x" << J.rhs() << "\n";
-}
-
-// test out the parameter infrastructure
-TEST(run, parameters)
-{
- nest::mc::parameter_list list("test");
- EXPECT_EQ(list.name(), "test");
- EXPECT_EQ(list.num_parameters(), 0);
-
- nest::mc::parameter p("a", 0.12, {0, 10});
-
- // add_parameter() returns a bool that indicates whether
- // it was able to successfull add the parameter
- EXPECT_EQ(list.add_parameter(std::move(p)), true);
- EXPECT_EQ(list.num_parameters(), 1);
-
- // test in place construction of a parameter
- EXPECT_EQ(list.add_parameter({"b", -3.0}), true);
- EXPECT_EQ(list.num_parameters(), 2);
-
- // check that adding a parameter that already exists returns false
- // and does not increase the number of parameters
- EXPECT_EQ(list.add_parameter({"b", -3.0}), false);
- EXPECT_EQ(list.num_parameters(), 2);
-
- auto &parms = list.parameters();
- EXPECT_EQ(parms[0].name, "a");
- EXPECT_EQ(parms[0].value, 0.12);
- EXPECT_EQ(parms[0].range.min, 0);
- EXPECT_EQ(parms[0].range.max, 10);
-
- EXPECT_EQ(parms[1].name, "b");
- EXPECT_EQ(parms[1].value, -3);
- EXPECT_EQ(parms[1].range.has_lower_bound(), false);
- EXPECT_EQ(parms[1].range.has_upper_bound(), false);
-}
-
diff --git a/tests/test_stimulus.cpp b/tests/test_stimulus.cpp
index 43be56f158e07e409c94b21cca7c293026b8d830..8fd9307b61a780ad93f430338f7d76c941e20bc0 100644
--- a/tests/test_stimulus.cpp
+++ b/tests/test_stimulus.cpp
@@ -4,7 +4,7 @@
TEST(stimulus, i_clamp)
{
- using namespace nestmc;
+ using namespace nest::mc;
// stimulus with delay 2, duration 0.5, amplitude 6.0
i_clamp stim(2.0, 0.5, 6.0);
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 17a836c2065b0a56954d26fa3e1297cbfc16de45..679260d656c9c7964f0eea599951c371886123f4 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -137,7 +137,7 @@ TEST(swc_parser, invalid_input)
{
// Non-contiguous numbering in branches is considered invalid
// 1
- // / \
+ // / \.
// 2 3
// /
// 4
diff --git a/tests/test_tree.cpp b/tests/test_tree.cpp
index a9b923170e14a42edab6a1deb070cb380e8d2378..a468d80ca8d4a61804d4fd0e242cad466bff8f8c 100644
--- a/tests/test_tree.cpp
+++ b/tests/test_tree.cpp
@@ -34,9 +34,9 @@ TEST(cell_tree, from_parent_index) {
{
//
// 0 0
- // / \ / \
+ // / \ / \.
// 1 4 => 1 2
- // / \
+ // / \.
// 2 5
// /
// 3
@@ -54,11 +54,11 @@ TEST(cell_tree, from_parent_index) {
{
//
// 0 0
- // /|\ /|\
+ // /|\ /|\.
// 1 4 6 => 1 2 3
- // / | \
+ // / | \.
// 2 5 7
- // / \
+ // / \.
// 3 8
//
std::vector<int> parent_index =
@@ -82,17 +82,17 @@ TEST(cell_tree, from_parent_index) {
{
//
// 0 0
- // /|\ /|\
+ // /|\ /|\.
// 1 4 6 => 1 2 3
- // / | \ / \
+ // / | \ / \.
// 2 5 7 4 5
- // / \
+ // / \.
// 3 8
- // / \
+ // / \.
// 9 11
- // / \
+ // / \.
// 10 12
- // \
+ // \.
// 13
//
std::vector<int> parent_index =
@@ -278,7 +278,7 @@ TEST(cell_tree, balance) {
EXPECT_EQ(t.parent(5), 2);
EXPECT_EQ(t.parent(6), 5);
- t.to_graphviz("cell.dot");
+ //t.to_graphviz("cell.dot");
}
}
@@ -292,8 +292,7 @@ TEST(cell_tree, json_load)
for(auto c : range(0,cell_data.size())) {
std::vector<int> parent_index = cell_data[c]["parent_index"];
cell_tree tree(parent_index);
- //tree.to_graphviz("cell_" + std::to_string(c) + ".dot");
- tree.to_graphviz("cell" + std::to_string(c) + ".dot");
+ //tree.to_graphviz("cell" + std::to_string(c) + ".dot");
}
}
diff --git a/tests/util.hpp b/tests/util.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7a0027448e46c860471aeed68ebd8ed8e411dfe0
--- /dev/null
+++ b/tests/util.hpp
@@ -0,0 +1,153 @@
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+#include <iomanip>
+
+#include <cmath>
+
+#include <util.hpp>
+#include <json/src/json.hpp>
+
+// helpful code for running tests
+// a bit messy: refactor when it gets heavier and obvious patterns emerge...
+
+namespace testing{
+
+[[gnu::unused]] static
+void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
+{
+ auto m = values.size();
+ if(!m) return;
+
+ std::ofstream fid(fname);
+ if(!fid.is_open()) return;
+
+ auto n = values[0].size();
+ for(const auto& v : values) {
+ if(n!=v.size()) {
+ std::cerr << "all output arrays must have the same length\n";
+ return;
+ }
+ }
+
+ for(auto i=0u; i<n; ++i) {
+ for(auto j=0u; j<m; ++j) {
+ fid << " " << values[j][i];
+ }
+ fid << "\n";
+ }
+}
+
+[[gnu::unused]] static
+nlohmann::json
+load_spike_data(const std::string& input_name)
+{
+ nlohmann::json cell_data;
+ auto fid = std::ifstream(input_name);
+ if(!fid.is_open()) {
+ std::cerr << "error : unable to open file " << input_name
+ << " : run the validation generation script first\n";
+ return {};
+ }
+
+ try {
+ fid >> cell_data;
+ }
+ catch (...) {
+ std::cerr << "error : incorrectly formatted json file " << input_name << "\n";
+ return {};
+ }
+ return cell_data;
+}
+
+template <typename T>
+std::vector<T> find_spikes(std::vector<T> const& v, T threshold, T dt)
+{
+ if(v.size()<2) {
+ return {};
+ }
+
+ std::vector<T> times;
+ for(auto i=1u; i<v.size(); ++i) {
+ if(v[i]>=threshold && v[i-1]<threshold) {
+ auto pos = (threshold-v[i-1]) / (v[i]-v[i-1]);
+ times.push_back((i-1+pos)*dt);
+ }
+ }
+
+ return times;
+}
+
+struct spike_comparison {
+ double min = std::numeric_limits<double>::quiet_NaN();
+ double max = std::numeric_limits<double>::quiet_NaN();
+ double mean = std::numeric_limits<double>::quiet_NaN();
+ double rms = std::numeric_limits<double>::quiet_NaN();
+ std::vector<double> diff;
+
+ // check whether initialized (i.e. has valid results)
+ bool is_valid() const {
+ return min == min;
+ }
+
+ // return maximum relative error
+ double max_relative_error() const {
+ if(!is_valid()) {
+ return std::numeric_limits<double>::quiet_NaN();
+ }
+
+ return *std::max_element(diff.begin(), diff.end());
+ }
+};
+
+[[gnu::unused]] static
+std::ostream&
+operator<< (std::ostream& o, spike_comparison const& spikes)
+{
+ // use snprintf because C++ is just awful for formatting output
+ char buffer[512];
+ snprintf(
+ buffer, sizeof(buffer),
+ "min,max = %10.8f,%10.8f | mean,rms = %10.8f,%10.8f | max_rel = %10.8f",
+ spikes.min, spikes.max, spikes.mean, spikes.rms,
+ spikes.max_relative_error()*100
+ );
+ return o << buffer;
+}
+
+template <typename T>
+spike_comparison compare_spikes(
+ std::vector<T> const& spikes,
+ std::vector<T> const& baseline)
+{
+ spike_comparison c;
+
+ // return default initialized (all NaN) if number of spikes differs
+ if(spikes.size() != baseline.size()) {
+ return c;
+ }
+
+ c.min = std::numeric_limits<double>::max();
+ c.max = 0.;
+ c.mean = 0.;
+ c.rms = 0.;
+
+ auto n = spikes.size();
+ for(auto i=0u; i<n; ++i) {
+ auto error = std::fabs(spikes[i] - baseline[i]);
+ c.min = std::min(c.min, error);
+ c.max = std::max(c.max, error);
+ c.mean += error;
+ c.rms += error*error;
+ // relative difference
+ c.diff.push_back(error/baseline[i]);
+ }
+
+ c.mean /= n;
+ c.rms = std::sqrt(c.rms/n);
+
+ return c;
+}
+
+} // namespace testing
diff --git a/tests/validate.cpp b/tests/validate.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d814fca1512dd79b22763a4f9a769f2984ca5269
--- /dev/null
+++ b/tests/validate.cpp
@@ -0,0 +1,7 @@
+#include "gtest.h"
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}
+
diff --git a/tests/validate_ball_and_stick.cpp b/tests/validate_ball_and_stick.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eac29e2e84f20653acdb7555714f9ed6c6e3e555
--- /dev/null
+++ b/tests/validate_ball_and_stick.cpp
@@ -0,0 +1,162 @@
+#include <fstream>
+
+#include "gtest.h"
+#include "util.hpp"
+
+#include <cell.hpp>
+#include <fvm.hpp>
+
+#include <json/src/json.hpp>
+
+// compares results with those generated by nrn/soma.py
+TEST(ball_and_stick, neuron_baseline)
+{
+ using namespace nest::mc;
+ using namespace nlohmann;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ // Soma with diameter 12.6157 um and HH channel
+ auto soma = cell.add_soma(12.6157/2.0);
+ soma->add_mechanism(hh_parameters());
+
+ // add dendrite of length 200 um and diameter 1 um with passive channel
+ auto dendrite = cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
+ dendrite->add_mechanism(pas_parameters());
+
+ dendrite->mechanism("membrane").set("r_L", 100); // no effect for single compartment cell
+
+ // add stimulus
+ cell.add_stimulus({1,1}, {5., 80., 0.3});
+
+ // load data from file
+ auto cell_data = testing::load_spike_data("../nrn/ball_and_stick.json");
+ EXPECT_TRUE(cell_data.size()>0);
+ if(cell_data.size()==0) return;
+
+ json& nrn =
+ *std::max_element(
+ cell_data.begin(), cell_data.end(),
+ [](json& lhs, json& rhs) {return lhs["nseg"]<rhs["nseg"];}
+ );
+
+ auto& measurements = nrn["measurements"];
+
+ double dt = nrn["dt"];
+ double tfinal = 100.; // ms
+ int nt = tfinal/dt;
+
+ // inline type for storing the results of a simulation along with
+ // the information required to compare two simulations for accuracy
+ struct result {
+ std::vector<std::vector<double>> spikes;
+ std::vector<std::vector<double>> baseline_spikes;
+ std::vector<testing::spike_comparison> comparisons;
+ std::vector<double> thresh;
+ int n_comparments;
+
+ result(int nseg, double dt,
+ std::vector<std::vector<double>> &v,
+ json& measurements
+ ) {
+ n_comparments = nseg;
+ baseline_spikes = {
+ measurements["soma"]["spikes"],
+ measurements["dend"]["spikes"],
+ measurements["clamp"]["spikes"]
+ };
+ thresh = {
+ measurements["soma"]["thresh"],
+ measurements["dend"]["thresh"],
+ measurements["clamp"]["thresh"]
+ };
+ for(auto i=0; i<3; ++i) {
+ // calculate the NEST MC spike times
+ spikes.push_back
+ (testing::find_spikes(v[i], thresh[i], dt));
+ // compare NEST MC and baseline NEURON spike times
+ comparisons.push_back
+ (testing::compare_spikes(spikes[i], baseline_spikes[i]));
+ }
+ }
+ };
+
+ std::vector<result> results;
+ for(auto run_index=0u; run_index<cell_data.size(); ++run_index) {
+ auto& run = cell_data[run_index];
+ int num_compartments = run["nseg"];
+ dendrite->set_compartments(num_compartments);
+ std::vector<std::vector<double>> v(3);
+
+ // make the lowered finite volume cell
+ fvm::fvm_cell<double, int> model(cell);
+ auto graph = cell.model();
+
+ // set initial conditions
+ using memory::all;
+ model.voltage()(all) = -65.;
+ model.initialize(); // have to do this _after_ initial conditions are set
+
+ // run the simulation
+ auto soma_comp = nest::mc::find_compartment_index({0,0}, graph);
+ auto dend_comp = nest::mc::find_compartment_index({1,0.5}, graph);
+ auto clamp_comp = nest::mc::find_compartment_index({1,1.0}, graph);
+ v[0].push_back(model.voltage()[soma_comp]);
+ v[1].push_back(model.voltage()[dend_comp]);
+ v[2].push_back(model.voltage()[clamp_comp]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ // save voltage at soma
+ v[0].push_back(model.voltage()[soma_comp]);
+ v[1].push_back(model.voltage()[dend_comp]);
+ v[2].push_back(model.voltage()[clamp_comp]);
+ }
+
+ results.push_back( {run["nseg"], dt, v, measurements} );
+ }
+
+ // print results
+ for(auto& r : results){
+ // select the location with the largest error for printing
+ auto m =
+ std::max_element(
+ r.comparisons.begin(), r.comparisons.end(),
+ [](testing::spike_comparison& l, testing::spike_comparison& r)
+ {return l.max_relative_error()<r.max_relative_error();}
+ );
+ std::cout << std::setw(5) << r.n_comparments
+ << " compartments : " << *m
+ << "\n";
+ }
+
+ // sort results in ascending order of compartments
+ std::sort(
+ results.begin(), results.end(),
+ [](const result& l, const result& r)
+ {return l.n_comparments<r.n_comparments;}
+ );
+
+ // the strategy for testing is the following:
+ // 1. test that the solution converges to the finest reference solution as
+ // the number of compartments increases (i.e. spatial resolution is
+ // refined)
+ for(auto i=1u; i<results.size(); ++i) {
+ for(auto j=0; j<3; ++j) {
+ EXPECT_TRUE(
+ results[i].comparisons[j].max_relative_error()
+ < results[i-1].comparisons[j].max_relative_error()
+ );
+ }
+ }
+
+ // 2. test that the best solution (i.e. with most compartments) matches the
+ // reference solution closely (less than 0.1% over the course of 100ms
+ // simulation)
+ for(auto j=0; j<3; ++j) {
+ EXPECT_TRUE(results.back().comparisons[j].max_relative_error()*100<0.1);
+ }
+}
+
diff --git a/tests/validate_soma.cpp b/tests/validate_soma.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b65b4ecb3a0de12d475c6380b24b43b824b65599
--- /dev/null
+++ b/tests/validate_soma.cpp
@@ -0,0 +1,163 @@
+#include <fstream>
+
+#include "gtest.h"
+#include "util.hpp"
+
+#include <cell.hpp>
+#include <fvm.hpp>
+
+#include <json/src/json.hpp>
+
+// compares results with those generated by nrn/soma.py
+// single compartment model with HH channels
+TEST(soma, 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 18.8um and HH channel
+ auto soma = cell.add_soma(18.8/2.0);
+ soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
+ soma->add_mechanism(hh_parameters());
+
+ // add stimulus to the soma
+ cell.add_stimulus({0,0.5}, {10., 100., 0.1});
+
+ // make the lowered finite volume cell
+ fvm::fvm_cell<double, int> model(cell);
+
+ // load data from file
+ auto cell_data = testing::load_spike_data("../nrn/soma.json");
+ EXPECT_TRUE(cell_data.size()>0);
+ if(cell_data.size()==0) return;
+
+ json& nrn =
+ *std::min_element(
+ cell_data.begin(), cell_data.end(),
+ [](json& lhs, json& rhs) {return lhs["dt"]<rhs["dt"];}
+ );
+ std::vector<double> nrn_spike_times = nrn["spikes"];
+
+ std::cout << "baseline with time step size " << nrn["dt"] << " ms\n";
+ std::cout << "baseline spikes : " << nrn["spikes"] << "\n";
+
+ for(auto& run : cell_data) {
+ std::vector<double> v;
+ double dt = run["dt"];
+
+ // set initial conditions
+ using memory::all;
+ model.voltage()(all) = -65.;
+ model.initialize(); // have to do this _after_ initial conditions are set
+
+ // run the simulation
+ auto tfinal = 120.; // ms
+ int nt = tfinal/dt;
+ v.push_back(model.voltage()[0]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ // save voltage at soma
+ v.push_back(model.voltage()[0]);
+ }
+
+ // get the spike times from the NEST MC simulation
+ auto nst_spike_times = testing::find_spikes(v, 0., dt);
+ // compare NEST MC and baseline NEURON results
+ auto comparison = testing::compare_spikes(nst_spike_times, nrn_spike_times);
+
+ // Assert that relative error is less than 1%.
+ // For a 100 ms simulation this asserts that the difference between
+ // NEST and the most accurate NEURON simulation is less than 1 ms.
+ std::cout << "MAX ERROR @ " << dt << " is " << comparison.max_relative_error()*100. << "\n";
+ EXPECT_TRUE(comparison.max_relative_error()*100. < 1.);
+ }
+}
+
+// check if solution converges
+TEST(soma, convergence)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ // Soma with diameter 18.8um and HH channel
+ auto soma = cell.add_soma(18.8/2.0);
+ soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
+ soma->add_mechanism(hh_parameters());
+
+ // add stimulus to the soma
+ cell.add_stimulus({0,0.5}, {10., 100., 0.1});
+
+ // make the lowered finite volume cell
+ fvm::fvm_cell<double, int> model(cell);
+
+ // generate baseline solution with small dt=0.0001
+ std::vector<double> baseline_spike_times;
+ {
+ auto dt = 1e-4;
+ std::vector<double> v;
+
+ // set initial conditions
+ using memory::all;
+ model.voltage()(all) = -65.;
+ model.initialize(); // have to do this _after_ initial conditions are set
+
+ // run the simulation
+ auto tfinal = 120.; // ms
+ int nt = tfinal/dt;
+ v.push_back(model.voltage()[0]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ // save voltage at soma
+ v.push_back(model.voltage()[0]);
+ }
+
+ baseline_spike_times = testing::find_spikes(v, 0., dt);
+ }
+
+ testing::spike_comparison previous;
+ for(auto dt : {0.05, 0.02, 0.01, 0.005, 0.001}) {
+ std::vector<double> v;
+
+ // set initial conditions
+ using memory::all;
+ model.voltage()(all) = -65.;
+ model.initialize(); // have to do this _after_ initial conditions are set
+
+ // run the simulation
+ auto tfinal = 120.; // ms
+ int nt = tfinal/dt;
+ v.push_back(model.voltage()[0]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ // save voltage at soma
+ v.push_back(model.voltage()[0]);
+ }
+
+ // get the spike times from the NEST MC simulation
+ auto nst_spike_times = testing::find_spikes(v, 0., dt);
+
+ // compare NEST MC and baseline NEURON results
+ auto comparison = testing::compare_spikes(nst_spike_times, baseline_spike_times);
+ std::cout << "dt " << std::setw(8) << dt << " : " << comparison << "\n";
+ if(!previous.is_valid()) {
+ previous = std::move(comparison);
+ }
+ else {
+ // Assert that relative error is less than 1%.
+ // For a 100 ms simulation this asserts that the difference between
+ // NEST and the most accurate NEURON simulation is less than 1 ms.
+ EXPECT_TRUE(comparison.max_relative_error() < previous.max_relative_error());
+ EXPECT_TRUE(comparison.rms < previous.rms);
+ EXPECT_TRUE(comparison.max < previous.max);
+ }
+ }
+}
diff --git a/vector b/vector
index e95bcbf1de6f833b5113a82dcc6d47bb7f41c397..8153d030abba20ddd174da9804406139c393a808 160000
--- a/vector
+++ b/vector
@@ -1 +1 @@
-Subproject commit e95bcbf1de6f833b5113a82dcc6d47bb7f41c397
+Subproject commit 8153d030abba20ddd174da9804406139c393a808