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 0157edd8b889d01de253f8545b06cbf91b9b8416..55d6af4d257348526713c1fdc381c4c2792651fd 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -40,13 +40,13 @@ flags = [
'-x',
'c++',
'-I',
+ 'src',
+ '-I',
'.',
- '-isystem',
- '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1'
- '-isystem',
- '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1/bits'
'-I',
- 'src',
+ 'include',
+# '-isystem',
+# '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
# '-I',
# '/usr/include/c++/4.9.2',
# '-isystem',
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
new file mode 100644
index 0000000000000000000000000000000000000000..975d601d09491f6b6418200b795368dbd9b14f7e
--- /dev/null
+++ b/data/hh.mod
@@ -0,0 +1,124 @@
+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 = 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/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/formulation.tex b/docs/formulation.tex
index 88e8e0a55f9a5bff24865cdc316c8f49ae814d78..0676ff4f9b03d1e1b3df4aa16963c858a4bb14e4 100644
--- a/docs/formulation.tex
+++ b/docs/formulation.tex
@@ -193,7 +193,7 @@ where $a_{i,\ell}$ and $a_{i,r}$ are the radii of at the left and right end of t
%-------------------------------------------------------------------------------
\subsubsection{Putting it all together}
%-------------------------------------------------------------------------------
-By substituting the volume averaging of the temporal derivative in~\eq{eq:dvdt} approximations for the flux over the surfaces in~\eq{eq:q_ij} and~\eq{eq:cv_volume} respectively into the conservation equation~\eq{eq:cable_balance} we get the following ODE defined for each node in the cell
+By substituting the volume averaging of the temporal derivative in~\eq{eq:dvdt} approximations for the flux over the surfaces in~\eq{eq:J_ij} and~\eq{eq:cv_volume} respectively into the conservation equation~\eq{eq:cable_balance} we get the following ODE defined for each node in the cell
\begin{equation}
\sigma_i c_m \dder{V_i}{t}
= -\sum_{j\in\mathcal{N}_i} {\frac{\sigma_{i,j}}{r_L \Delta x_{i,j}} (V_i-V_j)} - \sigma_i\cdot(i_m(V_i) - i_e(x_i)),
@@ -228,18 +228,25 @@ The current $i_m$ is often a nonlinear function of voltage, so if it was formula
The equations can be rearranged to have all unknown voltage values on the lhs, and values that can be calculated directly on the rhs:
\begin{align}
- & V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\frac{\alpha_{ij}}{\sigma_i} (V_i^{k+1}-V_j^{k+1})}
+ & \sigma_i V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\Delta t \alpha_{ij} (V_i^{k+1}-V_j^{k+1})}
\nonumber \\
- = & V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e),
+ = & \sigma_i \left( V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e) \right),
\label{eq:ode_linsys}
\end{align}
where the value
\begin{equation}
- \alpha_{ij} = \alpha_{ji} = \frac{\Delta t \sigma_{ij}}{ c_m r_L \Delta x_{ij}}
+ \alpha_{ij} = \alpha_{ji} = \frac{\sigma_{ij}}{ c_m r_L \Delta x_{ij}}
\label{eq:alpha_linsys}
\end{equation}
is a constant that can be computed for each interface between adjacent compartments during set up.
+The left hand side of \eq{eq:ode_linsys} can be rearranged
+\begin{equation}
+ \left[ \sigma_i + \sum_{j\in\mathcal{N}_i} {\Delta t \alpha_{ij}} \right] V_i^{k+1}
+ - \sum_{j\in\mathcal{N}_i} { \Delta t \alpha_{ij} V_j^{k+1}},
+\end{equation}
+which gives the coefficients for the linear system.
+
%-------------------------------------------------------------------------------
\subsubsection{Example: unbranched uniform cable}
%-------------------------------------------------------------------------------
@@ -247,18 +254,22 @@ For an unrbanched uniform cable of constant radius $a$, with length $L$ and $n$
\begin{align}
\Delta x_{ij} &= \Delta x = \frac{L}{n-1}, \nonumber \\
\sigma_{ij} &= \pi a^2, \nonumber \\
- \sigma_{i} &= 2 \pi a \Delta x, \nonumber \\
- \alpha_{ij} &= \frac{\pi a^2\Delta t}{c_m r_L\Delta x}, \nonumber \\
- \frac{\alpha_{ij}}{\sigma_i}
- &= \frac{a\Delta t}{2c_m r_L\Delta x^2}. \nonumber
+ \sigma_{i} &= \sigma_s = 2 \pi a \Delta x, \nonumber \\
+ \alpha_{ij} &= \alpha = \frac{\pi a^2}{c_m r_L\Delta x}, \nonumber
\end{align}
-With these simplifications, the lhs of the linear system is
+With these simplifications, the LHS of the linear system is
\begin{align}
- & V_i^{k+1} + \beta (V_i^{k+1}-V_{i+1}^{k+1}) + \beta (V_i^{k+1}-V_{i-1}^{k+1})
- \nonumber \\
- = & (1+2\beta)V_i^{k+1} - \beta V_{i+1}^{k+1} - \beta V_{i-1}^{k+1}.
+ \left[\sigma_s + 2\Delta t\alpha\right] & V_{i}^{k+1}
+ - \Delta t \alpha V_{i+1}^{k+1}
+ - \Delta t \alpha V_{i-1}^{k+1}
+ \nonumber \\
+ \left[\sigma_s/2 + \Delta t\alpha\right] & V_{1}^{k+1}
+ - \Delta t \alpha V_{2}^{k+1}
+ \nonumber \\
+ \left[\sigma_s/2 + \Delta t\alpha\right] & V_{n}^{k+1}
+ - \Delta t \alpha V_{n-1}^{k+1}
+ \nonumber
\end{align}
-where $\beta=\frac{a\Delta t}{2c_m r_L\Delta x^2}$.
The end points of the cable, i.e. the compartments for $x_1$ and $x_n$, have to be handled differently.
If we assume that a no-flux boundary condition, i.e. $\vv{J}\cdot\vv{n}=0$, is imposed at the end of the cable, the lhs of the linear system are
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 90c9d8f8c4446615bfd3e045cd64a082a0d1d016..8a8286bda8e314d47f5e89b2b454bd6b84709e67 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,7 +1,10 @@
set(HEADERS
swcio.hpp
- )
+)
set(BASE_SOURCES
+ cell.cpp
+ mechanism_interface.cpp
+ parameter_list.cpp
swcio.cpp
)
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b3943c688b33b7a2114f2dc88d2f3d22905d273c
--- /dev/null
+++ b/src/algorithms.hpp
@@ -0,0 +1,250 @@
+#pragma once
+
+#include <iostream>
+
+#include <algorithm>
+#include <numeric>
+#include <type_traits>
+#include <vector>
+
+#include "util.hpp"
+
+/*
+ * Some simple wrappers around stl algorithms to improve readability of code
+ * that uses them.
+ *
+ * For example, a simple sum() wrapper for finding the sum of a container,
+ * is much more readable than using std::accumulate()
+ *
+ */
+
+namespace nest {
+namespace mc {
+namespace algorithms{
+
+ template <typename C>
+ typename C::value_type
+ sum(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::accumulate(c.begin(), c.end(), value_type{0});
+ }
+
+ template <typename C>
+ typename C::value_type
+ mean(C const& c)
+ {
+ return sum(c)/c.size();
+ }
+
+ template <typename C>
+ C make_index(C const& c)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "make_index only applies to integral types"
+ );
+
+ C out(c.size()+1);
+ out[0] = 0;
+ std::partial_sum(c.begin(), c.end(), out.begin()+1);
+ return out;
+ }
+
+ /// works like std::is_sorted(), but with stronger condition that succesive
+ /// elements must be greater than those before them
+ template <typename C>
+ bool is_strictly_monotonic_increasing(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::is_sorted(
+ c.begin(),
+ c.end(),
+ [] (value_type const& lhs, value_type const& rhs) {
+ return lhs <= rhs;
+ }
+ );
+ }
+
+ template <typename C>
+ bool is_strictly_monotonic_decreasing(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::is_sorted(
+ c.begin(),
+ c.end(),
+ [] (value_type const& lhs, value_type const& rhs) {
+ return lhs >= rhs;
+ }
+ );
+ }
+
+ template <
+ typename C,
+ typename = typename std::enable_if<std::is_integral<typename C::value_type>::value>
+ >
+ bool is_minimal_degree(C const& c)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "is_minimal_degree only applies to integral types"
+ );
+
+ if(c.size()==0u) {
+ return true;
+ }
+
+ using value_type = typename C::value_type;
+ if(c[0] != value_type(0)) {
+ return false;
+ }
+ auto i = value_type(1);
+ auto it = std::find_if(
+ c.begin()+1, c.end(), [&i](value_type v) { return v>=(i++); }
+ );
+ return it==c.end();
+ }
+
+ template <typename C>
+ bool is_positive(C const& c)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "is_positive only applies to integral types"
+ );
+ for(auto v : c) {
+ if(v<1) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ template<typename C>
+ bool has_contiguous_segments(const C &parent_index)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ if (!is_minimal_degree(parent_index)) {
+ return false;
+ }
+
+ int n = parent_index.size();
+ std::vector<bool> is_leaf(n, false);
+
+ for(auto i=1; i<n; ++i) {
+ auto p = parent_index[i];
+ if(is_leaf[p]) {
+ return false;
+ }
+
+ if(p != i-1) {
+ // we have a branch and i-1 is a leaf node
+ is_leaf[i-1] = true;
+ }
+ }
+
+ return true;
+ }
+
+ template<typename C>
+ std::vector<typename C::value_type> child_count(const C &parent_index)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ std::vector<typename C::value_type> count(parent_index.size(), 0);
+ for (std::size_t i = 1; i < parent_index.size(); ++i) {
+ ++count[parent_index[i]];
+ }
+
+ return count;
+ }
+
+ template<typename C>
+ std::vector<typename C::value_type> branches(const C &parent_index)
+ {
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ EXPECTS(has_contiguous_segments(parent_index));
+
+ auto num_child = child_count(parent_index);
+ std::vector<typename C::value_type> branch_runs(
+ parent_index.size(), 0
+ );
+
+ std::size_t num_branches = (num_child[0] == 1) ? 1 : 0;
+ for (std::size_t i = 1; i < parent_index.size(); ++i) {
+ auto p = parent_index[i];
+ if (num_child[p] > 1) {
+ ++num_branches;
+ }
+
+ branch_runs[i] = num_branches;
+ }
+
+ return branch_runs;
+ }
+
+ template<typename C>
+ 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)
+ {
+ //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
+} // namespace mc
+} // namespace nest
diff --git a/src/cell.cpp b/src/cell.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6986f0673ca05827ee89493aca2230ec16cec30
--- /dev/null
+++ b/src/cell.cpp
@@ -0,0 +1,199 @@
+#include "cell.hpp"
+#include "tree.hpp"
+
+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
+ segments_.push_back(make_segment<placeholder_segment>());
+ parents_.push_back(0);
+}
+
+int cell::num_segments() const
+{
+ return segments_.size();
+}
+
+//
+// note: I think that we have to enforce that the soma is the first
+// segment that is added
+//
+soma_segment* cell::add_soma(value_type radius, point_type center)
+{
+ if(has_soma()) {
+ throw std::domain_error(
+ "attempt to add a soma to a cell that already has one"
+ );
+ }
+
+ // add segment for the soma
+ if(center.is_set()) {
+ segments_[0] = make_segment<soma_segment>(radius, center);
+ }
+ else {
+ segments_[0] = make_segment<soma_segment>(radius);
+ }
+
+ return segments_[0]->as_soma();
+}
+
+cable_segment* cell::add_cable(cell::index_type parent, segment_ptr&& cable)
+{
+ // check for a valid parent id
+ if(cable->is_soma()) {
+ throw std::domain_error(
+ "attempt to add a soma as a segment"
+ );
+ }
+
+ // check for a valid parent id
+ if(parent>num_segments()) {
+ throw std::out_of_range(
+ "parent index of cell segment is out of range"
+ );
+ }
+ segments_.push_back(std::move(cable));
+ parents_.push_back(parent);
+
+ return segments_.back()->as_cable();
+}
+
+segment* cell::segment(int index)
+{
+ if(index<0 || index>=num_segments()) {
+ throw std::out_of_range(
+ "attempt to access a segment with invalid index"
+ );
+ }
+ return segments_[index].get();
+}
+
+segment const* cell::segment(int index) const
+{
+ if(index<0 || index>=num_segments()) {
+ throw std::out_of_range(
+ "attempt to access a segment with invalid index"
+ );
+ }
+ return segments_[index].get();
+}
+
+
+bool cell::has_soma() const
+{
+ return !segment(0)->is_placeholder();
+}
+
+soma_segment* cell::soma()
+{
+ if(has_soma()) {
+ return segment(0)->as_soma();
+ }
+ return nullptr;
+}
+
+cable_segment* cell::cable(int index)
+{
+ if(index>0 && index<num_segments()) {
+ return segment(index)->as_cable();
+ }
+ return nullptr;
+}
+
+cell::value_type cell::volume() const
+{
+ return
+ std::accumulate(
+ segments_.begin(), segments_.end(),
+ 0.,
+ [](double value, segment_ptr const& seg) {
+ return seg->volume() + value;
+ }
+ );
+}
+
+cell::value_type cell::area() const
+{
+ return
+ std::accumulate(
+ segments_.begin(), segments_.end(),
+ 0.,
+ [](double value, segment_ptr const& seg) {
+ return seg->area() + value;
+ }
+ );
+}
+
+std::vector<segment_ptr> const& cell::segments() const
+{
+ return segments_;
+}
+
+std::vector<int> cell::compartment_counts() const
+{
+ std::vector<int> comp_count;
+ comp_count.reserve(num_segments());
+ for(auto const& s : segments()) {
+ comp_count.push_back(s->num_compartments());
+ }
+ return comp_count;
+}
+
+size_t cell::num_compartments() const
+{
+ auto n = 0u;
+ for(auto &s : segments_) {
+ n += s->num_compartments();
+ }
+ return n;
+}
+
+compartment_model cell::model() const
+{
+ compartment_model m;
+
+ m.tree = cell_tree(parents_);
+ auto counts = compartment_counts();
+ m.parent_index = make_parent_index(m.tree.graph(), counts);
+ m.segment_index = algorithms::make_index(counts);
+
+ 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_;
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/src/cell.hpp b/src/cell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ba5f9b6bff227231cbb4f46e27dd2af4ff2943ce
--- /dev/null
+++ b/src/cell.hpp
@@ -0,0 +1,144 @@
+#pragma once
+
+#include <mutex>
+#include <stdexcept>
+#include <thread>
+#include <vector>
+
+#include <segment.hpp>
+#include <cell_tree.hpp>
+#include <stimulus.hpp>
+
+namespace nest {
+namespace mc {
+
+/// wrapper around compartment layout information derived from a high level cell
+/// description
+struct compartment_model {
+ cell_tree tree;
+ std::vector<int> parent_index;
+ 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:
+
+ // types
+ using index_type = int;
+ using value_type = double;
+ using point_type = point<value_type>;
+
+ // constructor
+ cell();
+
+ /// add a soma to the cell
+ /// radius must be specified
+ 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
+ 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>
+ cable_segment* add_cable(index_type parent, Args ...args);
+
+ /// the number of segments in the cell
+ int num_segments() const;
+
+ bool has_soma() const;
+
+ class segment* segment(int index);
+ class segment const* segment(int index) const;
+
+ /// access pointer to the soma
+ /// returns nullptr if the cell has no soma
+ soma_segment* soma();
+
+ /// access pointer to a cable segment
+ /// will throw an std::out_of_range exception if
+ /// the cable index is not valid
+ cable_segment* cable(int index);
+
+ /// the volume of the cell
+ value_type volume() const;
+
+ /// the surface area of the cell
+ value_type area() const;
+
+ /// the total number of compartments over all segments
+ size_t num_compartments() const;
+
+ std::vector<segment_ptr> const& segments() const;
+
+ /// return reference to array that enumerates the index of the parent of
+ /// each segment
+ std::vector<int> const& segment_parents() const;
+
+ /// return a vector with the compartment count for each segment in the cell
+ std::vector<int> compartment_counts() const;
+
+ 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>
+cable_segment* cell::add_cable(cell::index_type parent, Args ...args)
+{
+ // check for a valid parent id
+ if(parent>=num_segments()) {
+ throw std::out_of_range(
+ "parent index of cell segment is out of range"
+ );
+ }
+ segments_.push_back(make_segment<cable_segment>(std::forward<Args>(args)...));
+ parents_.push_back(parent);
+
+ return segments_.back()->as_cable();
+}
+
+} // namespace mc
+} // namespace nest
+
diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp
index 0d4c496f459821c59caa80319d1536bba71b3671..43fac2429a657545251159a3fa97002b7243fe86 100644
--- a/src/cell_tree.hpp
+++ b/src/cell_tree.hpp
@@ -13,6 +13,9 @@
#include "tree.hpp"
#include "util.hpp"
+namespace nest {
+namespace mc {
+
/// The tree data structure that describes the segments of a cell tree.
/// A cell is represented as a tree where each node may have any number of
/// children. Typically in a cell only the soma has more than two segments,
@@ -25,38 +28,40 @@
/// sets it appears that the soma was always index 0, however we need more
/// flexibility in choosing the root.
class cell_tree {
- public :
-
+ using range = memory::Range;
+public :
// use a signed 16-bit integer for storage of indexes, which is reasonable given
// that typical cells have at most 1000-2000 segments
using int_type = int16_t;
using index_type = memory::HostVector<int_type>;
using index_view = index_type::view_type;
+ /// default empty constructor
+ cell_tree() = default;
+
/// construct from a parent index
cell_tree(std::vector<int> const& parent_index)
{
// handle the case of an empty parent list, which implies a single-compartment model
- std::vector<int> segment_index;
if(parent_index.size()>0) {
- segment_index = tree_.init_from_parent_index(parent_index);
+ tree_ = tree(parent_index);
}
else {
- segment_index = tree_.init_from_parent_index(std::vector<int>({0}));
+ tree_ = tree(std::vector<int>({0}));
}
-
- // if needed, calculate meta-data like length[] and end[] arrays for data
}
/// construct from a tree
// copy constructor
- cell_tree(tree const& t)
- : tree_(t)
+ cell_tree(tree const& t, int s)
+ : tree_(t),
+ soma_(s)
{ }
// move constructor
- cell_tree(tree&& t)
- : tree_(std::move(t))
+ cell_tree(tree&& t, int s)
+ : tree_(std::move(t)),
+ soma_(s)
{ }
/// construct from a cell tree
@@ -66,11 +71,31 @@ class cell_tree {
soma_(other.soma())
{ }
+ // assignment from rvalue
+ cell_tree& operator=(cell_tree&& other)
+ {
+ std::swap(other.tree_, tree_);
+ std::swap(other.soma_, soma_);
+ return *this;
+ }
+
+ // assignment
+ cell_tree& operator=(cell_tree const& other)
+ {
+ tree_ = other.tree_;
+ soma_ = other.soma_;
+ return *this;
+ }
+
// move constructor
cell_tree(cell_tree&& other)
- : tree_(std::move(other.tree_)),
- soma_(other.soma())
- { }
+ {
+ *this = std::move(other);
+ }
+
+ tree const& graph() const {
+ return tree_;
+ }
int_type soma() const {
return soma_;
@@ -148,8 +173,7 @@ class cell_tree {
return depth;
}
- private :
-
+private :
/// helper type for sub-tree computation
/// use in balance()
@@ -166,12 +190,11 @@ class cell_tree {
}
std::string to_string() const {
- std::string s;
-
- s += "[" + std::to_string(root) + ","
- + std::to_string(diameter) + "," + std::to_string(depth) + "]";
-
- return s;
+ return
+ "[" + std::to_string(root) + ","
+ + std::to_string(diameter) + ","
+ + std::to_string(depth) +
+ "]";
}
int root;
@@ -261,8 +284,17 @@ class cell_tree {
}
}
- // storage for the tree structure of cell segments
+ //////////////////////////////////////////////////
+ // state
+ //////////////////////////////////////////////////
+
+ /// storage for the tree structure of cell segments
tree tree_;
+ /// index of the soma
int_type soma_ = 0;
};
+
+} // namespace mc
+} // namespace nest
+
diff --git a/src/compartment.hpp b/src/compartment.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b92360a560ba5de30a036505d50d457b7b5a5a9c
--- /dev/null
+++ b/src/compartment.hpp
@@ -0,0 +1,166 @@
+#pragma once
+
+#include <iterator>
+#include <utility>
+
+namespace nest {
+namespace mc {
+
+/// Defines the simplest type of compartment
+/// The compartment is a conic frustrum
+struct compartment {
+ using value_type = double;
+ using size_type = int;
+ using value_pair = std::pair<value_type, value_type>;
+
+ compartment() = delete;
+
+ compartment(
+ size_type idx,
+ value_type len,
+ value_type r1,
+ value_type r2
+ )
+ : index{idx},
+ radius{r1, r2},
+ length{len}
+ {}
+
+
+ size_type index;
+ std::pair<value_type, value_type> radius;
+ value_type length;
+};
+
+/// The simplest type of compartment iterator :
+/// - divide a segment into n compartments of equal length
+/// - assume that the radius varies linearly from one end of the segment
+/// to the other
+class compartment_iterator :
+ public std::iterator<std::forward_iterator_tag, compartment>
+{
+
+ public:
+
+ using base = std::iterator<std::forward_iterator_tag, compartment>;
+ using size_type = base::value_type::size_type;
+ using real_type = base::value_type::value_type;
+
+ compartment_iterator() = delete;
+
+ compartment_iterator(
+ size_type idx,
+ real_type len,
+ real_type rad,
+ real_type delta_rad
+ )
+ : index_(idx),
+ radius_(rad),
+ delta_radius_(delta_rad),
+ length_(len)
+ { }
+
+ compartment_iterator(compartment_iterator const& other) = default;
+
+ compartment_iterator& operator++()
+ {
+ index_++;
+ radius_ += delta_radius_;
+ return *this;
+ }
+
+ compartment_iterator operator++(int)
+ {
+ compartment_iterator ret(*this);
+ operator++();
+ return ret;
+ }
+
+ compartment operator*() const
+ {
+ return
+ compartment(
+ index_, length_, radius_, radius_ + delta_radius_
+ );
+ }
+
+ bool operator==(compartment_iterator const& other) const
+ {
+ return other.index_ == index_;
+ }
+
+ bool operator!=(compartment_iterator const& other) const
+ {
+ return !operator==(other);
+ }
+
+ private :
+
+ size_type index_;
+ real_type radius_;
+ const real_type delta_radius_;
+ const real_type length_;
+};
+
+class compartment_range {
+ public:
+
+ using size_type = compartment_iterator::size_type;
+ using real_type = compartment_iterator::real_type;
+
+ compartment_range(
+ size_type num_compartments,
+ real_type radius_L,
+ real_type radius_R,
+ real_type length
+ )
+ : num_compartments_(num_compartments),
+ radius_L_(radius_L),
+ radius_R_(radius_R),
+ length_(length)
+ {}
+
+ compartment_iterator begin() const
+ {
+ return {0, compartment_length(), radius_L_, delta_radius()};
+ }
+
+ compartment_iterator cbegin() const
+ {
+ return begin();
+ }
+
+ /// With 0-based indexing compartment number "num_compartments_" is
+ /// one past the end
+ compartment_iterator end() const
+ {
+ return {num_compartments_, 0, 0, 0};
+ }
+
+ compartment_iterator cend() const
+ {
+ return end();
+ }
+
+ real_type delta_radius() const
+ {
+ return (radius_R_ - radius_L_) / num_compartments_;
+ }
+
+ real_type compartment_length() const
+ {
+ return length_ / num_compartments_;
+ }
+
+ private:
+
+ size_type num_compartments_;
+ real_type radius_L_;
+ real_type radius_R_;
+ real_type length_;
+};
+
+} // namespace mc
+} // namespace nest
+
+
diff --git a/src/fvm.hpp b/src/fvm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..69d23d0d1fe97e341ad8de3dc25ce8c0771ef78c
--- /dev/null
+++ b/src/fvm.hpp
@@ -0,0 +1,493 @@
+#pragma once
+
+#include <algorithm>
+#include <map>
+#include <string>
+#include <vector>
+
+#include <algorithms.hpp>
+#include <cell.hpp>
+#include <ion.hpp>
+#include <math.hpp>
+#include <matrix.hpp>
+#include <mechanism.hpp>
+#include <mechanism_interface.hpp>
+#include <util.hpp>
+#include <segment.hpp>
+#include <stimulus.hpp>
+
+#include <vector/include/Vector.hpp>
+
+namespace nest {
+namespace mc {
+namespace fvm {
+
+template <typename T, typename I>
+class fvm_cell {
+ public :
+
+ /// the real number type
+ using value_type = T;
+ /// 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
+ using index_view = typename index_type::view_type;
+
+ /// the container used for values
+ using vector_type = memory::HostVector<value_type>;
+ /// view into value container
+ using vector_view = typename vector_type::view_type;
+
+ /// constructor
+ fvm_cell(nest::mc::cell const& cell);
+
+ /// build the matrix for a given time step
+ void setup_matrix(value_type dt);
+
+ /// TODO this should be const
+ /// which requires const_view in the vector library
+ matrix_type& jacobian() {
+ 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() {
+ return cv_areas_;
+ }
+
+ /// TODO this should be const
+ /// return the capacitance of each CV surface
+ /// this is the total capacitance, not per unit area,
+ /// i.e. equivalent to sigma_i * c_m
+ vector_view cv_capacitance() {
+ return cv_capacitance_;
+ }
+
+ /// 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_;
+
+ /// cv_areas_[i] is the surface area of CV i
+ vector_type cv_areas_;
+
+ /// alpha_[i] is the following value at the CV face between
+ /// CV i and its parent, required when constructing linear system
+ /// face_alpha_[i] = area_face / (c_m * r_L * delta_x);
+ vector_type face_alpha_;
+
+ /// cv_capacitance_[i] is the capacitance of CV i per unit area (i.e. c_m)
+ vector_type cv_capacitance_;
+
+ /// the average current over the surface of each CV
+ /// current_ = i_m - i_e
+ /// so the total current over the surface of CV i is
+ /// current_[i] * cv_areas_
+ vector_type current_;
+
+ /// 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 ////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+template <typename T, typename I>
+fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
+: cv_areas_ {cell.num_compartments(), T(0)}
+, face_alpha_ {cell.num_compartments(), T(0)}
+, cv_capacitance_{cell.num_compartments(), T(0)}
+, current_ {cell.num_compartments(), T(0)}
+, voltage_ {cell.num_compartments(), T(0)}
+{
+ 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;
+
+ auto seg_idx = 0;
+ for(auto const& s : cell.segments()) {
+ if(auto soma = s->as_soma()) {
+ // assert the assumption that the soma is at 0
+ if(seg_idx!=0) {
+ throw std::domain_error(
+ "FVM lowering encountered soma with non-zero index"
+ );
+ }
+ auto area = math::area_sphere(soma->radius());
+ cv_areas_[0] += area;
+ cv_capacitance_[0] += area * soma->mechanism("membrane").get("c_m").value;
+ }
+ else if(auto cable = s->as_cable()) {
+ // loop over each compartment in the cable
+ // each compartment has the face between two CVs at its centre
+ // the centers of the CVs are the end points of the compartment
+ //
+ // __________________________________
+ // | ........ | .cvleft. | cv |
+ // | ........ L ........ C R
+ // |__________|__________|__________|
+ //
+ // The compartment has end points marked L and R (left and right).
+ // The left compartment is assumed to be closer to the soma
+ // (i.e. it follows the minimal degree ordering)
+ // The face is at the center, marked C.
+ // The full control volume to the left (marked with .)
+ auto c_m = cable->mechanism("membrane").get("c_m").value;
+ auto r_L = cable->mechanism("membrane").get("r_L").value;
+ for(auto c : cable->compartments()) {
+ auto i = segment_index[seg_idx] + c.index;
+ auto j = parent_index[i];
+
+ 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);
+
+ auto halflen = c.length/2;
+
+ auto al = math::area_frustrum(halflen, left(c.radius), radius_center);
+ auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
+ cv_areas_[j] += al;
+ cv_areas_[i] += ar;
+ cv_capacitance_[j] += al * c_m;
+ cv_capacitance_[i] += ar * c_m;
+ }
+ }
+ else {
+ throw std::domain_error("FVM lowering encountered unsuported segment type");
+ }
+ ++seg_idx;
+ }
+
+ // normalize the capacitance by cv_area
+ for(auto i=0u; i<size(); ++i) {
+ cv_capacitance_[i] /= cv_areas_[i];
+ }
+
+ /////////////////////////////////////////////
+ // create mechanisms
+ /////////////////////////////////////////////
+
+ // FIXME : candidate for a private member function
+
+ // 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>
+void fvm_cell<T, I>::setup_matrix(T dt)
+{
+ using memory::all;
+
+ // convenience accesors to matrix storage
+ auto l = matrix_.l();
+ auto d = matrix_.d();
+ auto u = matrix_.u();
+ auto p = matrix_.p();
+ auto rhs = matrix_.rhs();
+
+ // The matrix has the following layout in memory
+ // where j is the parent index of i, i.e. i<j
+ //
+ // d[i] is the diagonal entry at a_ii
+ // u[i] is the upper triangle entry at a_ji
+ // l[i] is the lower triangle entry at a_ij
+ //
+ // d[j] . . u[i]
+ // . . .
+ // . . .
+ // l[i] . . d[i]
+ //
+ //d(all) = 1.0;
+ d(all) = cv_areas_;
+ for(auto i=1u; i<d.size(); ++i) {
+ // TODO get this right
+ // probably requires scaling a by cv_areas_[i] and cv_areas_[p[i]]
+ auto a = 1e5*dt * face_alpha_[i];
+
+ d[i] += a;
+ l[i] = -a;
+ u[i] = -a;
+
+ // add contribution to the diagonal of parent
+ d[p[i]] += a;
+ }
+ //std::cout << "d " << d << " l " << l << " u " << u << "\n";
+
+ // the RHS of the linear system is
+ // V[i] - dt/cm*(im - ie)
+ auto factor = 10.*dt;
+ for(auto i=0u; i<d.size(); ++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/helpers.hpp b/src/helpers.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cae1e30a5b7b43835cced5292d7f06e11b1c5d7f
--- /dev/null
+++ b/src/helpers.hpp
@@ -0,0 +1,15 @@
+#include <algorithm>
+
+namespace nestmc {
+namespace range{
+
+ template <typename C>
+ typename C::value_type
+ sum(C const& c)
+ {
+ using value_type = typename C::value_type;
+ return std::accumulate(c.begin(), c.end(), value_type{0});
+ }
+
+}
+}
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/math.hpp b/src/math.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..86eccfce67e7f53fe6645b15808c201c88634e19
--- /dev/null
+++ b/src/math.hpp
@@ -0,0 +1,83 @@
+#pragma once
+
+#include <utility>
+#include <cmath>
+
+namespace nest {
+namespace mc {
+namespace math {
+ template <typename T>
+ T constexpr pi()
+ {
+ return T(3.1415926535897932384626433832795);
+ }
+
+ template <typename T>
+ T constexpr mean(T a, T b)
+ {
+ return (a+b) / T(2);
+ }
+
+ template <typename T>
+ T constexpr mean(std::pair<T,T> const& p)
+ {
+ return (p.first+p.second) / T(2);
+ }
+
+ template <typename T>
+ T constexpr square(T a)
+ {
+ return a*a;
+ }
+
+ template <typename T>
+ T constexpr cube(T a)
+ {
+ return a*a*a;
+ }
+
+ // area of a circle
+ // pi r-squared
+ template <typename T>
+ T constexpr area_circle(T r)
+ {
+ return pi<T>() * square(r);
+ }
+
+ // volume of a conic frustrum
+ template <typename T>
+ T constexpr area_frustrum(T L, T r1, T r2)
+ {
+ return pi<T>() * (r1+r2) * sqrt(square(L) + square(r1-r2));
+ }
+
+ // surface area of conic frustrum, not including the area of the
+ // circles at either end (i.e. just the area of the surface of rotation)
+ template <typename T>
+ T constexpr volume_frustrum(T L, T r1, T r2)
+ {
+ // this formulation uses one less multiplication
+ return pi<T>()/T(3) * (square(r1+r2) - r1*r2) * L;
+ //return pi<T>()/T(3) * (square(r1) + r1*r2 + square(r2)) * L;
+ }
+
+ // volume of a sphere
+ // four-thirds pi r-cubed
+ template <typename T>
+ T constexpr volume_sphere(T r)
+ {
+ return T(4)/T(3) * pi<T>() * cube(r);
+ }
+
+ // surface area of a sphere
+ // four pi r-squared
+ template <typename T>
+ T constexpr area_sphere(T r)
+ {
+ return T(4) * pi<T>() * square(r);
+ }
+
+} // namespace math
+} // namespace mc
+} // namespace nest
+
diff --git a/src/matrix.hpp b/src/matrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..db9d2df7f56e8bfffb73051e996acb13d9fa0369
--- /dev/null
+++ b/src/matrix.hpp
@@ -0,0 +1,209 @@
+#pragma once
+
+#include <type_traits>
+
+#include "util.hpp"
+
+#include <vector/include/Vector.hpp>
+
+namespace nest {
+namespace mc {
+
+/// matrix storage structure
+template<typename T, typename I>
+class matrix {
+ public:
+
+ // define basic types
+ using value_type = T;
+ using size_type = I;
+
+ // define storage types
+ using vector_type = memory::HostVector<value_type>;
+ using vector_view_type = typename vector_type::view_type;
+ using index_type = memory::HostVector<size_type>;
+ using index_view_type = typename index_type::view_type;
+
+ matrix() = default;
+
+ /// construct matrix for one or more cells, with combined parent index
+ /// and a cell index
+ template <
+ typename LHS, typename RHS,
+ typename = typename
+ std::enable_if<
+ util::is_container<LHS>::value &&
+ util::is_container<RHS>::value
+ >
+ >
+ matrix(LHS&& pi, RHS&& ci)
+ : parent_index_(std::forward<LHS>(pi))
+ , cell_index_(std::forward<RHS>(ci))
+ {
+ setup();
+ }
+
+ /// construct matrix for a single cell described by a parent index
+ template <
+ typename IDX,
+ typename = typename
+ std::enable_if< util::is_container<IDX>::value >
+ >
+ matrix(IDX&& pi)
+ : parent_index_(std::forward<IDX>(pi))
+ , cell_index_(2)
+ {
+ cell_index_[0] = 0;
+ cell_index_[1] = size();
+ setup();
+ }
+
+ /// the dimension of the matrix (i.e. the number of rows or colums)
+ std::size_t size() const
+ {
+ return parent_index_.size();
+ }
+
+ /// the total memory used to store the matrix
+ std::size_t memory() const
+ {
+ 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);
+ return s;
+ }
+
+ /// the number of cell matrices that have been packed together
+ size_type num_cells() const
+ {
+ 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
+ vector_view_type l()
+ {
+ return l_;
+ }
+
+ /// the vector holding the diagonal of the matrix
+ vector_view_type d()
+ {
+ return d_;
+ }
+
+ /// the vector holding the upper part of the matrix
+ vector_view_type u()
+ {
+ return u_;
+ }
+
+ /// the vector holding the right hand side of the linear equation system
+ 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_type p()
+ {
+ return parent_index_;
+ }
+
+ /// solve the linear system
+ /// upon completion the solution is stored in the RHS storage
+ /// and can be accessed via rhs()
+ void solve()
+ {
+ /*
+ 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
+ for(auto m=0; m<ncells; ++m) {
+ auto first = cell_index_[m];
+ auto last = cell_index_[m+1];
+
+ // backward sweep
+ for(auto i=last-1; i>first; --i) {
+ auto factor = l_[i] / d_[i];
+ d_[p[i]] -= factor * u_[i];
+ rhs_[p[i]] -= factor * rhs_[i];
+ }
+ rhs_[first] /= d_[first];
+
+ // forward sweep
+ for(auto i=first+1; i<last; ++i) {
+ rhs_[i] -= u_[i] * rhs_[p[i]];
+ 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:
+
+ /// allocate memory for storing matrix and right hand side vector
+ void setup()
+ {
+ const auto n = size();
+ constexpr auto default_value
+ = std::numeric_limits<value_type>::quiet_NaN();
+
+ l_ = vector_type(n, default_value);
+ d_ = vector_type(n, default_value);
+ u_ = vector_type(n, default_value);
+ rhs_ = vector_type(n, default_value);
+ }
+
+ /// the parent indice that describe matrix structure
+ index_type parent_index_;
+
+ /// indexes that point to the start of each cell in the index
+ index_type cell_index_;
+
+ /// storage for lower, diagonal, upper and rhs
+ 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
+} // namespace mc
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea69c292e1f5c50efe0fd0e6e7dbf092c32b16fb
--- /dev/null
+++ b/src/mechanism.hpp
@@ -0,0 +1,93 @@
+#pragma once
+
+#pragma once
+
+#include <memory>
+#include <string>
+
+#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;
+
+ // 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 indexed_view_type = indexed_view<value_type, size_type>;
+
+ using ion_type = ion<value_type, size_type>;
+
+ 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_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 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;
+
+ view_type vec_v_;
+ view_type vec_i_;
+ index_type node_index_;
+};
+
+template <typename T, typename I>
+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::view_type vec_v,
+ typename M::view_type vec_i,
+ typename M::index_view 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
new file mode 100644
index 0000000000000000000000000000000000000000..10fad45fa0d9776c4f45302e6bf196f995765e99
--- /dev/null
+++ b/src/mechanism_interface.cpp
@@ -0,0 +1,44 @@
+#include "mechanism_interface.hpp"
+
+//
+// include the mechanisms
+//
+
+#include <mechanisms/hh.hpp>
+#include <mechanisms/pas.hpp>
+
+namespace nest {
+namespace mc {
+namespace mechanisms {
+
+std::map<std::string, mechanism_helper_ptr<value_type, index_type>> mechanism_helpers;
+
+void setup_mechanism_helpers() {
+ 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
+} // namespace nest
+} // namespace mc
+
diff --git a/src/mechanism_interface.hpp b/src/mechanism_interface.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c59caaad6c655d6a69c5206b956b1cb42959cfdf
--- /dev/null
+++ b/src/mechanism_interface.hpp
@@ -0,0 +1,58 @@
+#pragma once
+
+#include <map>
+#include <string>
+
+#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 value_type = T;
+ using size_type = I;
+ using index_type = memory::HostVector<I>;
+ using index_view = typename index_type::view_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_ptr<T,I> new_mechanism(view_type, view_type, index_view) const = 0;
+
+ virtual void set_parameters(mechanism_ptr_type&, parameter_list const&) const = 0;
+};
+
+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();
+
+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.cpp b/src/parameter_list.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..86c641aa809e6efe1c8f03247cc64f9f1b3a01f0
--- /dev/null
+++ b/src/parameter_list.cpp
@@ -0,0 +1,107 @@
+#include <algorithm>
+#include <iostream>
+
+#include "parameter_list.hpp"
+
+namespace nest {
+namespace mc {
+
+ bool parameter_list::add_parameter(parameter p)
+ {
+ if(has_parameter(p.name)) {
+ return false;
+ }
+
+ parameters_.push_back(std::move(p));
+
+ return true;
+ }
+
+ bool parameter_list::has_parameter(std::string const& n) const
+ {
+ return find_by_name(n) != parameters_.end();
+ }
+
+ int parameter_list::num_parameters() const
+ {
+ return parameters_.size();
+ }
+
+ // returns true if parameter was succesfully updated
+ // returns false if parameter was not updated, i.e. if
+ // - no parameter with name n exists
+ // - value is not in the valid range
+ bool parameter_list::set(std::string const& n, parameter_list::value_type v)
+ {
+ auto p = find_by_name(n);
+ if(p!=parameters_.end()) {
+ if(p->is_in_range(v)) {
+ p->value = v;
+ return true;
+ }
+ }
+ return false;
+ }
+
+ parameter& parameter_list::get(std::string const& n)
+ {
+ auto it = find_by_name(n);
+ if(it==parameters_.end()) {
+ throw std::domain_error(
+ "parameter list does not contain parameter"
+ );
+ }
+ return *it;
+ }
+
+ std::string const& parameter_list::name() const {
+ return mechanism_name_;
+ }
+
+ std::vector<parameter> const& parameter_list::parameters() const
+ {
+ return parameters_;
+ }
+
+ auto parameter_list::find_by_name(std::string const& n)
+ -> decltype(parameters_.begin())
+ {
+ return
+ std::find_if(
+ parameters_.begin(), parameters_.end(),
+ [&n](parameter const& p) {return p.name == n;}
+ );
+ }
+
+ auto parameter_list::find_by_name(std::string const& n) const
+ -> decltype(parameters_.begin())
+ {
+ return
+ std::find_if(
+ parameters_.begin(), parameters_.end(),
+ [&n](parameter const& p) {return p.name == n;}
+ );
+ }
+} // namespace mc
+} // namespace nest
+
+std::ostream&
+operator<<(std::ostream& o, nest::mc::parameter const& p)
+{
+ return o
+ << "parameter("
+ << "name " << p.name
+ << " : value " << p.value
+ << " : range " << p.range
+ << ")";
+}
+
+std::ostream&
+operator<<(std::ostream& o, nest::mc::parameter_list const& l)
+{
+ o << "parameters \"" << l.name() << "\" :\n";
+ for(nest::mc::parameter const& p : l.parameters()) {
+ o << " " << p << "\n";
+ }
+ return o;
+}
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c06988756a077fdbb860fc8f42043b178b1b7e7
--- /dev/null
+++ b/src/parameter_list.hpp
@@ -0,0 +1,230 @@
+#pragma once
+
+#include <limits>
+#include <ostream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace nest {
+namespace mc {
+
+ template <typename T>
+ struct value_range {
+ using value_type = T;
+
+ value_range()
+ : min(std::numeric_limits<value_type>::quiet_NaN()),
+ max(std::numeric_limits<value_type>::quiet_NaN())
+ { }
+
+ value_range(value_type const& left, value_type const& right)
+ : min(left),
+ max(right)
+ {
+ if(left>right) {
+ throw std::out_of_range(
+ "parameter range must have left <= right"
+ );
+ }
+ }
+
+ bool has_lower_bound() const
+ {
+ return min==min;
+ }
+
+ bool has_upper_bound() const
+ {
+ return max==max;
+ }
+
+ bool is_in_range(value_type const& v) const
+ {
+ if(has_lower_bound()) {
+ if(v<min) return false;
+ }
+ if(has_upper_bound()) {
+ if(v>max) return false;
+ }
+ return true;
+ }
+
+ value_type min;
+ value_type max;
+ };
+
+ struct parameter {
+ using value_type = double;
+ using range_type = value_range<value_type>;
+
+ parameter() = delete;
+
+ parameter(std::string n, value_type v, range_type r=range_type())
+ : name(std::move(n)),
+ value(v),
+ range(r)
+ {
+ if(!is_in_range(v)) {
+ throw std::out_of_range(
+ "parameter value is out of permitted value range"
+ );
+ }
+ }
+
+ bool is_in_range(value_type v) const
+ {
+ return range.is_in_range(v);
+ }
+
+ std::string name;
+ value_type value;
+ range_type range;
+ };
+
+ // Use a dumb container class for now
+ // might have to use a more sophisticated interface in the future if need be
+ class parameter_list {
+ public :
+
+ using value_type = double;
+
+ parameter_list(std::string n)
+ : mechanism_name_(std::move(n))
+ { }
+
+ bool has_parameter(std::string const& n) const;
+ bool add_parameter(parameter p);
+
+ // returns true if parameter was succesfully updated
+ // returns false if parameter was not updated, i.e. if
+ // - no parameter with name n exists
+ // - value is not in the valid range
+ bool set(std::string const& n, value_type v);
+ parameter& get(std::string const& n);
+
+ std::string const& name() const;
+
+ std::vector<parameter> const& parameters() const;
+
+ int num_parameters() const;
+
+ private:
+
+ std::vector<parameter> parameters_;
+ std::string mechanism_name_;
+
+ // need two versions for const and non-const applications
+ auto find_by_name(std::string const& n)
+ -> decltype(parameters_.begin());
+
+ auto find_by_name(std::string const& n) const
+ -> decltype(parameters_.begin());
+
+ };
+
+ ///////////////////////////////////////////////////////////////////////////
+ // predefined parameter sets
+ ///////////////////////////////////////////////////////////////////////////
+
+ /// default set of parameters for the cell membrane that are added to every
+ /// segment when it is created.
+ class membrane_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;
+
+ membrane_parameters()
+ : base("membrane")
+ {
+ // 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
+ }
+ };
+
+ /// parameters for the classic Hodgkin & Huxley model (1952)
+ class hh_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;
+
+ 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});
+ }
+ };
+
+ /// 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});
+ }
+ };
+
+} // namespace mc
+} // namespace nest
+
+template <typename T>
+std::ostream& operator<<(std::ostream& o, nest::mc::value_range<T> const& r)
+{
+ o << "[";
+ if(r.has_lower_bound())
+ o << r.min;
+ else
+ o<< "-inf";
+ o << ", ";
+ if(r.has_upper_bound())
+ o << r.max;
+ else
+ o<< "inf";
+ return o << "]";
+}
+
+std::ostream& operator<<(std::ostream& o, nest::mc::parameter const& p);
+std::ostream& operator<<(std::ostream& o, nest::mc::parameter_list const& l);
+
diff --git a/src/point.hpp b/src/point.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3b048bb3db7fb580443e019c98292821f715ffb2
--- /dev/null
+++ b/src/point.hpp
@@ -0,0 +1,97 @@
+#pragma once
+
+#include <limits>
+#include <ostream>
+
+namespace nest {
+namespace mc {
+
+template <typename T>
+struct point {
+ using value_type = T;
+ value_type x;
+ value_type y;
+ value_type z;
+
+ constexpr point(T a, T b, T c)
+ : x{a}, y{b}, z{c}
+ {}
+
+ // initialize to NaN by default
+ constexpr point()
+ : x(std::numeric_limits<T>::quiet_NaN()),
+ y(std::numeric_limits<T>::quiet_NaN()),
+ z(std::numeric_limits<T>::quiet_NaN())
+ {}
+
+ constexpr bool is_set() const
+ {
+ return (x==x && y==y && z==z);
+ }
+};
+
+template <typename T>
+constexpr point<T>
+operator+ (
+ point<T> const& lhs,
+ point<T> const& rhs
+) {
+ return point<T>(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
+}
+
+template <typename T>
+constexpr point<T>
+operator- (
+ point<T> const& lhs,
+ point<T> const& rhs
+) {
+ return point<T>(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z);
+}
+
+template <typename T>
+constexpr point<T>
+operator* (
+ T alpha,
+ point<T> const& p
+) {
+ return point<T>(alpha*p.x, alpha*p.y, alpha*p.z);
+}
+
+template <typename T>
+T
+norm(point<T> const& p)
+{
+ return sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
+}
+
+template <typename T>
+constexpr T
+dot(
+ point<T> const& lhs,
+ point<T> const& rhs
+) {
+ return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
+}
+
+template<typename T>
+bool operator==(const point<T> &rhs, const point<T> &lhs)
+{
+ return (rhs.x == lhs.x) &&
+ (rhs.y == lhs.y) &&
+ (rhs.z == lhs.z);
+}
+
+template<typename T>
+bool operator!=(const point<T> &rhs, const point<T> &lhs)
+{
+ return !(rhs == lhs);
+}
+
+} // namespace mc
+} // namespace nest
+
+template <typename T>
+std::ostream& operator << (std::ostream& o, nest::mc::point<T> const& p)
+{
+ return o << "[" << p.x << ", " << p.y << ", " << p.z << "]";
+}
diff --git a/src/segment.hpp b/src/segment.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fd1d8ded6101cc8e9834b9062df3da97de5cbc8c
--- /dev/null
+++ b/src/segment.hpp
@@ -0,0 +1,441 @@
+#pragma once
+
+#include <cmath>
+
+#include <vector>
+
+#include "compartment.hpp"
+#include "math.hpp"
+#include "parameter_list.hpp"
+#include "point.hpp"
+#include "algorithms.hpp"
+#include "util.hpp"
+
+namespace nest {
+namespace mc {
+
+template <typename T,
+ typename valid = typename std::is_floating_point<T>::type>
+struct segment_properties {
+ T rL = 180.0; // resistivity [Ohm.cm]
+ T cm = 0.01; // capacitance [F/m^2] : 10 nF/mm^2 = 0.01 F/m^2
+};
+
+enum class segmentKind {
+ soma,
+ dendrite,
+ axon,
+ none
+};
+
+// forward declarations of segment specializations
+class soma_segment;
+class cable_segment;
+
+// abstract base class for a cell segment
+class segment {
+ public:
+
+ using value_type = double;
+ using point_type = point<value_type>;
+
+ segmentKind kind() const {
+ return kind_;
+ }
+
+ bool is_soma() const
+ {
+ return kind_==segmentKind::soma;
+ }
+
+ bool is_dendrite() const
+ {
+ return kind_==segmentKind::dendrite;
+ }
+
+ bool is_axon() const
+ {
+ return kind_==segmentKind::axon;
+ }
+
+ virtual int num_compartments() const = 0;
+ virtual void set_compartments(int) = 0;
+
+ virtual value_type volume() const = 0;
+ virtual value_type area() const = 0;
+
+ virtual ~segment() = default;
+
+ virtual cable_segment* as_cable()
+ {
+ return nullptr;
+ }
+
+ virtual soma_segment* as_soma()
+ {
+ return nullptr;
+ }
+
+ virtual bool is_placeholder() const
+ {
+ return false;
+ }
+
+ segment_properties<value_type> properties;
+
+ void add_mechanism(parameter_list p)
+ {
+ auto it = std::find_if(
+ mechanisms_.begin(), mechanisms_.end(),
+ [&p](parameter_list const& l){return l.name() == p.name();}
+ );
+ if(it!=mechanisms_.end()) {
+ throw std::out_of_range(
+ "attempt to add a mechanism parameter set to a segment which has an existing mechanism with the same name"
+ );
+ }
+
+ mechanisms_.push_back(std::move(p));
+ }
+
+ parameter_list& mechanism(std::string const& n)
+ {
+ 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;
+ }
+
+ 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_;
+ std::vector<parameter_list> mechanisms_;
+};
+
+class placeholder_segment : public segment
+{
+ public:
+
+ using base = segment;
+ using base::kind_;
+ using base::value_type;
+
+ placeholder_segment()
+ {
+ kind_ = segmentKind::none;
+ }
+
+ value_type volume() const override
+ {
+ return std::numeric_limits<value_type>::quiet_NaN();
+ }
+
+ value_type area() const override
+ {
+ return std::numeric_limits<value_type>::quiet_NaN();
+ }
+
+ bool is_placeholder() const override
+ {
+ return true;
+ }
+
+ int num_compartments() const override
+ {
+ return 0;
+ }
+
+ virtual void set_compartments(int) override
+ { }
+};
+
+class soma_segment : public segment
+{
+ public :
+
+ using base = segment;
+ using base::kind_;
+ using base::value_type;
+ using base::point_type;
+
+ soma_segment() = delete;
+
+ soma_segment(value_type r)
+ : radius_{r}
+ {
+ kind_ = segmentKind::soma;
+ mechanisms_.push_back(membrane_parameters());
+ }
+
+ soma_segment(value_type r, point_type const &c)
+ : soma_segment(r)
+ {
+ center_ = c;
+ }
+
+ value_type volume() const override
+ {
+ return math::volume_sphere(radius_);
+ }
+
+ value_type area() const override
+ {
+ return math::area_sphere(radius_);
+ }
+
+ value_type radius() const
+ {
+ return radius_;
+ }
+
+ point_type const& center() const
+ {
+ return center_;
+ }
+
+ soma_segment* as_soma() override
+ {
+ return this;
+ }
+
+ /// soma has one and one only compartments
+ int num_compartments() const override
+ {
+ return 1;
+ }
+
+ void set_compartments(int n) override
+ { }
+
+ private :
+
+ // store the center and radius of the soma
+ // note that the center may be undefined
+ value_type radius_;
+ point_type center_;
+};
+
+class cable_segment : public segment
+{
+ public :
+
+ using base = segment;
+ using base::kind_;
+ using base::value_type;
+ using base::point_type;
+
+ // delete the default constructor
+ cable_segment() = delete;
+
+ // constructors for a cable with no location information
+ cable_segment(
+ segmentKind k,
+ std::vector<value_type> r,
+ std::vector<value_type> lens
+ ) {
+ kind_ = k;
+ assert(k==segmentKind::dendrite || k==segmentKind::axon);
+
+ radii_ = std::move(r);
+ lengths_ = std::move(lens);
+
+ // add default membrane parameters
+ mechanisms_.push_back(membrane_parameters());
+ }
+
+ cable_segment(
+ segmentKind k,
+ value_type r1,
+ value_type r2,
+ value_type len
+ )
+ : cable_segment{k, std::vector<value_type>{r1, r2}, std::vector<value_type>{len}}
+ { }
+
+ // constructor that lets the user describe the cable as a
+ // seriew of radii and locations
+ cable_segment(
+ segmentKind k,
+ std::vector<value_type> r,
+ std::vector<point_type> p
+ ) {
+ kind_ = k;
+ assert(k==segmentKind::dendrite || k==segmentKind::axon);
+
+ radii_ = std::move(r);
+ locations_ = std::move(p);
+ update_lengths();
+
+ // add default membrane parameters
+ mechanisms_.push_back(membrane_parameters());
+ }
+
+ // helper that lets user specify with the radius and location
+ // of just the end points of the cable
+ // i.e. describing the cable as a single frustrum
+ cable_segment(
+ segmentKind k,
+ value_type r1,
+ value_type r2,
+ point_type const& p1,
+ point_type const& p2
+ )
+ : cable_segment(k, {r1, r2}, {p1, p2})
+ { }
+
+ value_type volume() const override
+ {
+ auto sum = value_type{0};
+ for(auto i=0; i<num_sub_segments(); ++i) {
+ sum += math::volume_frustrum(lengths_[i], radii_[i], radii_[i+1]);
+ }
+ return sum;
+ }
+
+ value_type area() const override
+ {
+ auto sum = value_type{0};
+ for(auto i=0; i<num_sub_segments(); ++i) {
+ sum += math::area_frustrum(lengths_[i], radii_[i], radii_[i+1]);
+ }
+ return sum;
+ }
+
+ value_type length() const
+ {
+ return algorithms::sum(lengths_);
+ }
+
+ bool has_locations() const
+ {
+ return locations_.size() > 0;
+ }
+
+ // the number sub-segments that define the cable segment
+ int num_sub_segments() const
+ {
+ return radii_.size()-1;
+ }
+
+ std::vector<value_type> const& lengths() const
+ {
+ return lengths_;
+ }
+
+ cable_segment* as_cable() override
+ {
+ return this;
+ }
+
+ int num_compartments() const override
+ {
+ return num_compartments_;
+ }
+
+ void set_compartments(int n) override
+ {
+ if(n<1) {
+ throw std::out_of_range(
+ "number of compartments in a segment must be one or more"
+ );
+ }
+ num_compartments_ = n;
+ }
+
+ value_type radius(value_type loc) const
+ {
+ if(loc>=1.) return radii_.back();
+ if(loc<=0.) return radii_.front();
+
+ auto len = length();
+ value_type pos = loc*len;
+
+ // This could be cached using a partial sum.
+ // In fact a lot of this stuff can be cached if
+ // we find ourselves having to do it over and over again.
+ // The time to cache it might be when update_lengths() is called.
+ auto sum = value_type(0);
+ auto i=0;
+ for(i=0; i<num_sub_segments(); ++i) {
+ if(sum+lengths_[i]>pos) {
+ break;
+ }
+ sum += lengths_[i];
+ }
+
+ auto rel = (len - sum)/lengths_[i];
+
+ return rel*radii_[i] + (1.-rel)*radii_[i+1];
+ }
+
+ /// iterable range type for simple compartment representation
+ compartment_range compartments() const
+ {
+ return {num_compartments(), radii_.front(), radii_.back(), length()};
+ }
+
+ private :
+
+ void update_lengths()
+ {
+ if(locations_.size()) {
+ lengths_.resize(num_sub_segments());
+ for(auto i=0; i<num_sub_segments(); ++i) {
+ lengths_[i] = norm(locations_[i] - locations_[i+1]);
+ }
+ }
+ }
+
+ int num_compartments_ = 1;
+ std::vector<value_type> lengths_;
+ std::vector<value_type> radii_;
+ std::vector<point_type> locations_;
+};
+
+/// Unique pointer wrapper for abstract segment base class
+using segment_ptr = std::unique_ptr<segment>;
+
+/// Helper for constructing segments in a segment_ptr unique pointer wrapper.
+/// Forwards the supplied arguments to construct a segment of type SegmentType.
+/// e.g. auto my_cable = make_segment<cable>(segmentKind::dendrite, ... );
+template <typename SegmentType, typename... Args>
+segment_ptr make_segment(Args&&... args)
+{
+ return segment_ptr(new SegmentType(std::forward<Args>(args)...));
+}
+
+} // namespace mc
+} // namespace nest
+
diff --git a/src/stimulus.hpp b/src/stimulus.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f7c7538cf398a7154588126620d97b7c054953a4
--- /dev/null
+++ b/src/stimulus.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+namespace nest {
+namespace mc {
+
+class i_clamp {
+ public:
+
+ using value_type = double;
+
+ i_clamp(value_type del, value_type dur, value_type amp)
+ : delay_(del),
+ duration_(dur),
+ amplitude_(amp)
+ {}
+
+ value_type delay() const {
+ return delay_;
+ }
+ value_type duration() const {
+ return duration_;
+ }
+ value_type amplitude() const {
+ return amplitude_;
+ }
+
+ void set_delay(value_type d) {
+ delay_ = d;
+ }
+ void set_duration(value_type d) {
+ duration_ = d;
+ }
+ void set_amplitude(value_type a) {
+ amplitude_ = a;
+ }
+
+ // current is set to amplitude for time in the half open interval:
+ // t \in [delay, delay+duration)
+ value_type amplitude(double t) {
+ if(t>=delay_ && t<(delay_+duration_)) {
+ return amplitude_;
+ }
+ return 0;
+ }
+
+ private:
+
+ value_type delay_ = 0;
+ value_type duration_ = 0;
+ value_type amplitude_ = 0;
+};
+
+} // namespace mc
+} // namespace nest
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 28e189a6f5a19df31c9ee961ffe0fe79be1a7676..f96676449d9526154d4f7d0cd3531da62bba715e 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -4,17 +4,19 @@
#include <sstream>
#include <unordered_set>
+#include <algorithms.hpp>
+#include <point.hpp>
#include <swcio.hpp>
+#include <util.hpp>
-namespace nestmc
-{
-namespace io
-{
+namespace nest {
+namespace mc {
+namespace io {
//
-// cell_record implementation
+// swc_record implementation
//
-void cell_record::renumber(id_type new_id, std::map<id_type, id_type> &idmap)
+void swc_record::renumber(id_type new_id, std::map<id_type, id_type> &idmap)
{
auto old_id = id_;
id_ = new_id;
@@ -29,12 +31,13 @@ void cell_record::renumber(id_type new_id, std::map<id_type, id_type> &idmap)
idmap.insert(std::make_pair(old_id, new_id));
}
-void cell_record::check_consistency() const
+void swc_record::check_consistency() const
{
- // Check cell type as well; enum's do not offer complete type safety,
+ // Check record type as well; enum's do not offer complete type safety,
// since you can cast anything that fits to its underlying type
- if (type_ < 0 || type_ > custom) {
- throw std::invalid_argument("unknown cell type");
+ if (static_cast<int>(type_) < 0 ||
+ static_cast<int>(type_) > static_cast<int>(kind::custom)) {
+ throw std::invalid_argument("unknown record type");
}
if (id_ < 0) {
@@ -54,24 +57,24 @@ void cell_record::check_consistency() const
}
}
-std::istream &operator>>(std::istream &is, cell_record &cell)
+std::istream &operator>>(std::istream &is, swc_record &record)
{
swc_parser parser;
- parser.parse_record(is, cell);
+ parser.parse_record(is, record);
return is;
}
-std::ostream &operator<<(std::ostream &os, const cell_record &cell)
+std::ostream &operator<<(std::ostream &os, const swc_record &record)
{
// output in one-based indexing
- os << cell.id_+1 << " "
- << cell.type_ << " "
- << std::setprecision(7) << cell.x_ << " "
- << std::setprecision(7) << cell.y_ << " "
- << std::setprecision(7) << cell.z_ << " "
- << std::setprecision(7) << cell.r_ << " "
- << ((cell.parent_id_ == -1) ? cell.parent_id_ : cell.parent_id_+1);
+ os << record.id_+1 << " "
+ << static_cast<int>(record.type_) << " "
+ << std::setprecision(7) << record.x_ << " "
+ << std::setprecision(7) << record.y_ << " "
+ << std::setprecision(7) << record.z_ << " "
+ << std::setprecision(7) << record.r_ << " "
+ << ((record.parent_id_ == -1) ? record.parent_id_ : record.parent_id_+1);
return os;
}
@@ -104,22 +107,22 @@ T parse_value_strict(std::istream &is, const swc_parser &parser)
return val;
}
-// specialize parsing for cell types
+// specialize parsing for record types
template<>
-cell_record::kind parse_value_strict(std::istream &is, const swc_parser &parser)
+swc_record::kind parse_value_strict(std::istream &is, const swc_parser &parser)
{
- cell_record::id_type val;
+ swc_record::id_type val;
check_parse_status(is >> val, parser);
- // Let cell_record's constructor check for the type validity
- return static_cast<cell_record::kind>(val);
+ // Let swc_record's constructor check for the type validity
+ return static_cast<swc_record::kind>(val);
}
//
// swc_parser implementation
//
-std::istream &swc_parser::parse_record(std::istream &is, cell_record &cell)
+std::istream &swc_parser::parse_record(std::istream &is, swc_record &record)
{
while (!is.eof() && !is.bad()) {
// consume empty and comment lines first
@@ -146,7 +149,7 @@ std::istream &swc_parser::parse_record(std::istream &is, cell_record &cell)
std::istringstream line(linebuff_);
try {
- cell = parse_record(line);
+ record = parse_record(line);
} catch (std::invalid_argument &e) {
// Rethrow as a parse error
throw swc_parse_error(e.what(), lineno_);
@@ -155,67 +158,165 @@ std::istream &swc_parser::parse_record(std::istream &is, cell_record &cell)
return is;
}
-cell_record swc_parser::parse_record(std::istringstream &is)
+swc_record swc_parser::parse_record(std::istringstream &is)
{
auto id = parse_value_strict<int>(is, *this);
- auto type = parse_value_strict<cell_record::kind>(is, *this);
- auto x = parse_value_strict<float>(is, *this);
- auto y = parse_value_strict<float>(is, *this);
- auto z = parse_value_strict<float>(is, *this);
- auto r = parse_value_strict<float>(is, *this);
- auto parent_id = parse_value_strict<cell_record::id_type>(is, *this);
+ auto type = parse_value_strict<swc_record::kind>(is, *this);
+ auto x = parse_value_strict<swc_record::coord_type>(is, *this);
+ auto y = parse_value_strict<swc_record::coord_type>(is, *this);
+ auto z = parse_value_strict<swc_record::coord_type>(is, *this);
+ auto r = parse_value_strict<swc_record::coord_type>(is, *this);
+ auto parent_id = parse_value_strict<swc_record::id_type>(is, *this);
// Convert to zero-based, leaving parent_id as-is if -1
if (parent_id != -1) {
parent_id--;
}
- return cell_record(type, id-1, x, y, z, r, parent_id);
+ return swc_record(type, id-1, x, y, z, r, parent_id);
}
-cell_record_range_clean::cell_record_range_clean(std::istream &is)
+swc_record_range_clean::swc_record_range_clean(std::istream &is)
{
- std::unordered_set<cell_record::id_type> ids;
+ std::unordered_set<swc_record::id_type> ids;
- std::size_t num_trees = 0;
- cell_record::id_type last_id = -1;
- bool needsort = false;
+ std::size_t num_trees = 0;
+ swc_record::id_type last_id = -1;
+ bool needsort = false;
- cell_record curr_cell;
- for (auto c : swc_get_records<swc_io_raw>(is)) {
- if (c.parent() == -1 && ++num_trees > 1) {
+ swc_record curr_record;
+ for (auto r : swc_get_records<swc_io_raw>(is)) {
+ if (r.parent() == -1 && ++num_trees > 1) {
// only a single tree is allowed
break;
}
- auto inserted = ids.insert(c.id());
+ auto inserted = ids.insert(r.id());
if (inserted.second) {
- // not a duplicate; insert cell
- cells_.push_back(c);
- if (!needsort && c.id() < last_id) {
+ // not a duplicate; insert record
+ records_.push_back(r);
+ if (!needsort && r.id() < last_id) {
needsort = true;
}
- last_id = c.id();
+ last_id = r.id();
}
}
if (needsort) {
- std::sort(cells_.begin(), cells_.end());
+ std::sort(records_.begin(), records_.end());
}
- // Renumber cells if necessary
- std::map<cell_record::id_type, cell_record::id_type> idmap;
- cell_record::id_type next_id = 0;
- for (auto &c : cells_) {
- if (c.id() != next_id) {
- c.renumber(next_id, idmap);
+ // Renumber records if necessary
+ std::map<swc_record::id_type, swc_record::id_type> idmap;
+ swc_record::id_type next_id = 0;
+ for (auto &r : records_) {
+ if (r.id() != next_id) {
+ r.renumber(next_id, idmap);
}
++next_id;
}
+
+ // Reject if branches are not contiguously numbered
+ std::vector<swc_record::id_type> parent_list = { 0 };
+ for (std::size_t i = 1; i < records_.size(); ++i) {
+ parent_list.push_back(records_[i].parent());
+ }
+
+ if (!nest::mc::algorithms::has_contiguous_segments(parent_list)) {
+ throw swc_parse_error("branches are not contiguously numbered", 0);
+ }
+}
+
+//
+// Convenience functions for returning the radii and the coordinates of a series
+// of swc records
+//
+static std::vector<swc_record::coord_type>
+swc_radii(const std::vector<swc_record> &records)
+{
+ std::vector<swc_record::coord_type> radii;
+ for (const auto &r : records) {
+ radii.push_back(r.radius());
+ }
+
+ return radii;
+}
+
+static std::vector<nest::mc::point<swc_record::coord_type> >
+swc_points(const std::vector<swc_record> &records)
+{
+ std::vector<nest::mc::point<swc_record::coord_type> > points;
+ for (const auto &r : records) {
+ points.push_back(r.coord());
+ }
+
+ return points;
+}
+
+static void make_cable(cell &cell,
+ const std::vector<swc_record::id_type> &branch_index,
+ const std::vector<swc_record> &branch_run)
+{
+ auto new_parent = branch_index[branch_run.back().id()] - 1;
+ cell.add_cable(new_parent, nest::mc::segmentKind::dendrite,
+ swc_radii(branch_run), swc_points(branch_run));
+}
+
+cell swc_read_cell(std::istream &is)
+{
+ cell newcell;
+ std::vector<swc_record::id_type> parent_list;
+ std::vector<swc_record> swc_records;
+ for (const auto &r : swc_get_records<swc_io_clean>(is)) {
+ swc_records.push_back(r);
+ parent_list.push_back(r.parent());
+ }
+
+ // The parent of soma must be 0
+ if (!parent_list.empty()) {
+ parent_list[0] = 0;
+ }
+
+ auto branch_index = nest::mc::algorithms::branches(parent_list);
+ std::vector<swc_record> branch_run;
+
+ branch_run.reserve(parent_list.size());
+ auto last_branch_point = branch_index[0];
+ for (auto i = 0u; i < swc_records.size(); ++i) {
+ if (branch_index[i] != last_branch_point) {
+ // New branch encountered; add to cell the current one
+ const auto &p = branch_run.back();
+ if (p.parent() == -1) {
+ // This is a soma
+ newcell.add_soma(p.radius(), p.coord());
+ last_branch_point = i;
+ } else {
+ last_branch_point = i - 1;
+ make_cable(newcell, branch_index, branch_run);
+ }
+
+ // Reset the branch run
+ branch_run.clear();
+ if (p.parent() != -1) {
+ // Add parent of the current cell to the branch,
+ // if not branching from soma
+ branch_run.push_back(swc_records[parent_list[i]]);
+ }
+ }
+
+ branch_run.push_back(swc_records[i]);
+ }
+
+ if (!branch_run.empty()) {
+ make_cable(newcell, branch_index, branch_run);
+ }
+
+ return newcell;
}
-} // end of nestmc::io
-} // end of nestmc
+} // namespace io
+} // namespace mc
+} // namespace nest
diff --git a/src/swcio.hpp b/src/swcio.hpp
index e267c184b4e462b2015fb8c38866e2133d60b8b1..a58bbac9aee9a32c03ae9bb33eb9289782a537b2 100644
--- a/src/swcio.hpp
+++ b/src/swcio.hpp
@@ -6,23 +6,21 @@
#include <string>
#include <vector>
-namespace nestmc
-{
-
-namespace io
-{
+#include <cell.hpp>
+#include <point.hpp>
+namespace nest {
+namespace mc {
+namespace io {
-class cell_record
+class swc_record
{
public:
using id_type = int;
+ using coord_type = double;
- // FIXME: enum's are not completely type-safe, since they can accept
- // anything that can be casted to their underlying type.
- //
// More on SWC files: http://research.mssm.edu/cnic/swc.html
- enum kind {
+ enum class kind {
undefined = 0,
soma,
axon,
@@ -33,10 +31,10 @@ public:
custom
};
- // cell records assume zero-based indexing; root's parent remains -1
- cell_record(kind type, int id,
- float x, float y, float z, float r,
- int parent_id)
+ // swc records assume zero-based indexing; root's parent remains -1
+ swc_record(swc_record::kind type, int id,
+ coord_type x, coord_type y, coord_type z, coord_type r,
+ int parent_id)
: type_(type)
, id_(id)
, x_(x)
@@ -48,8 +46,8 @@ public:
check_consistency();
}
- cell_record()
- : type_(cell_record::undefined)
+ swc_record()
+ : type_(swc_record::kind::undefined)
, id_(0)
, x_(0)
, y_(0)
@@ -58,10 +56,10 @@ public:
, parent_id_(-1)
{ }
- cell_record(const cell_record &other) = default;
- cell_record &operator=(const cell_record &other) = default;
+ swc_record(const swc_record &other) = default;
+ swc_record &operator=(const swc_record &other) = default;
- bool strict_equals(const cell_record &other) const
+ bool strict_equals(const swc_record &other) const
{
return id_ == other.id_ &&
x_ == other.x_ &&
@@ -72,43 +70,43 @@ public:
}
// Equality and comparison operators
- friend bool operator==(const cell_record &lhs,
- const cell_record &rhs)
+ friend bool operator==(const swc_record &lhs,
+ const swc_record &rhs)
{
return lhs.id_ == rhs.id_;
}
- friend bool operator<(const cell_record &lhs,
- const cell_record &rhs)
+ friend bool operator<(const swc_record &lhs,
+ const swc_record &rhs)
{
return lhs.id_ < rhs.id_;
}
- friend bool operator<=(const cell_record &lhs,
- const cell_record &rhs)
+ friend bool operator<=(const swc_record &lhs,
+ const swc_record &rhs)
{
return (lhs < rhs) || (lhs == rhs);
}
- friend bool operator!=(const cell_record &lhs,
- const cell_record &rhs)
+ friend bool operator!=(const swc_record &lhs,
+ const swc_record &rhs)
{
return !(lhs == rhs);
}
- friend bool operator>(const cell_record &lhs,
- const cell_record &rhs)
+ friend bool operator>(const swc_record &lhs,
+ const swc_record &rhs)
{
return !(lhs < rhs) && (lhs != rhs);
}
- friend bool operator>=(const cell_record &lhs,
- const cell_record &rhs)
+ friend bool operator>=(const swc_record &lhs,
+ const swc_record &rhs)
{
return !(lhs < rhs);
}
- friend std::ostream &operator<<(std::ostream &os, const cell_record &cell);
+ friend std::ostream &operator<<(std::ostream &os, const swc_record &record);
kind type() const
{
@@ -125,41 +123,46 @@ public:
return parent_id_;
}
- float x() const
+ coord_type x() const
{
return x_;
}
- float y() const
+ coord_type y() const
{
return y_;
}
- float z() const
+ coord_type z() const
{
return z_;
}
- float radius() const
+ coord_type radius() const
{
return r_;
}
- float diameter() const
+ coord_type diameter() const
{
return 2*r_;
}
+ nest::mc::point<coord_type> coord() const
+ {
+ return nest::mc::point<coord_type>(x_, y_, z_);
+ }
+
void renumber(id_type new_id, std::map<id_type, id_type> &idmap);
private:
void check_consistency() const;
- kind type_; // cell type
- id_type id_; // cell id
- float x_, y_, z_; // cell coordinates
- float r_; // cell radius
- id_type parent_id_; // cell parent's id
+ kind type_; // record type
+ id_type id_; // record id
+ coord_type x_, y_, z_; // record coordinates
+ coord_type r_; // record radius
+ id_type parent_id_; // record parent's id
};
@@ -206,11 +209,11 @@ public:
return lineno_;
}
- std::istream &parse_record(std::istream &is, cell_record &cell);
+ std::istream &parse_record(std::istream &is, swc_record &record);
private:
// Read the record from a string stream; will be treated like a single line
- cell_record parse_record(std::istringstream &is);
+ swc_record parse_record(std::istringstream &is);
std::string delim_;
std::string comment_prefix_;
@@ -219,15 +222,15 @@ private:
};
-std::istream &operator>>(std::istream &is, cell_record &cell);
+std::istream &operator>>(std::istream &is, swc_record &record);
-class cell_record_stream_iterator :
- public std::iterator<std::forward_iterator_tag, cell_record>
+class swc_record_stream_iterator :
+ public std::iterator<std::forward_iterator_tag, swc_record>
{
public:
struct eof_tag { };
- cell_record_stream_iterator(std::istream &is)
+ swc_record_stream_iterator(std::istream &is)
: is_(is)
, eof_(false)
{
@@ -236,19 +239,19 @@ public:
read_next_record();
}
- cell_record_stream_iterator(std::istream &is, eof_tag)
+ swc_record_stream_iterator(std::istream &is, eof_tag)
: is_(is)
, eof_(true)
{ }
- cell_record_stream_iterator(const cell_record_stream_iterator &other)
+ swc_record_stream_iterator(const swc_record_stream_iterator &other)
: is_(other.is_)
, parser_(other.parser_)
, curr_record_(other.curr_record_)
, eof_(other.eof_)
{ }
- cell_record_stream_iterator &operator++()
+ swc_record_stream_iterator &operator++()
{
if (eof_) {
throw std::out_of_range("attempt to read past eof");
@@ -258,14 +261,14 @@ public:
return *this;
}
- cell_record_stream_iterator operator++(int)
+ swc_record_stream_iterator operator++(int)
{
- cell_record_stream_iterator ret(*this);
+ swc_record_stream_iterator ret(*this);
operator++();
return ret;
}
- value_type operator*()
+ value_type operator*() const
{
if (eof_) {
throw std::out_of_range("attempt to read past eof");
@@ -274,7 +277,7 @@ public:
return curr_record_;
}
- bool operator==(const cell_record_stream_iterator &other) const
+ bool operator==(const swc_record_stream_iterator &other) const
{
if (eof_ && other.eof_) {
return true;
@@ -283,13 +286,13 @@ public:
}
}
- bool operator!=(const cell_record_stream_iterator &other)
+ bool operator!=(const swc_record_stream_iterator &other) const
{
return !(*this == other);
}
friend std::ostream &operator<<(std::ostream &os,
- const cell_record_stream_iterator &iter)
+ const swc_record_stream_iterator &iter)
{
os << "{ is_.tellg(): " << iter.is_.tellg() << ", "
<< "curr_record_: " << iter.curr_record_ << ", "
@@ -309,7 +312,7 @@ private:
std::istream &is_;
swc_parser parser_;
- cell_record curr_record_;
+ swc_record curr_record_;
// indicator of eof; we need a way to define an end() iterator without
// seeking to the end of file
@@ -317,28 +320,33 @@ private:
};
-class cell_record_range_raw
+class swc_record_range_raw
{
public:
- using value_type = cell_record;
- using reference = value_type &;
- using const_referene = const value_type &;
- using iterator = cell_record_stream_iterator;
- using const_iterator = const cell_record_stream_iterator;
+ using value_type = swc_record;
+ using reference = value_type &;
+ using const_reference = const value_type &;
+ using iterator = swc_record_stream_iterator;
+ using const_iterator = const swc_record_stream_iterator;
- cell_record_range_raw(std::istream &is)
+ swc_record_range_raw(std::istream &is)
: is_(is)
{ }
- iterator begin()
+ iterator begin() const
{
- return cell_record_stream_iterator(is_);
+ return swc_record_stream_iterator(is_);
}
- iterator end()
+ iterator end() const
{
iterator::eof_tag eof;
- return cell_record_stream_iterator(is_, eof);
+ return swc_record_stream_iterator(is_, eof);
+ }
+
+ bool empty() const
+ {
+ return begin() == end();
}
private:
@@ -346,58 +354,66 @@ private:
};
//
-// Reads cells from an input stream until an eof is encountered and returns a
-// cleaned sequence of cell records.
+// Reads records from an input stream until an eof is encountered and returns a
+// cleaned sequence of swc records.
//
// For more information check here:
// https://github.com/eth-cscs/cell_algorithms/wiki/SWC-file-parsing
//
-class cell_record_range_clean
+class swc_record_range_clean
{
public:
- using value_type = cell_record;
+ using value_type = swc_record;
using reference = value_type &;
using const_referene = const value_type &;
- using iterator = std::vector<cell_record>::iterator;
- using const_iterator = std::vector<cell_record>::const_iterator;
+ using iterator = std::vector<swc_record>::iterator;
+ using const_iterator = std::vector<swc_record>::const_iterator;
- cell_record_range_clean(std::istream &is);
+ swc_record_range_clean(std::istream &is);
iterator begin()
{
- return cells_.begin();
+ return records_.begin();
}
iterator end()
{
- return cells_.end();
+ return records_.end();
}
std::size_t size()
{
- return cells_.size();
+ return records_.size();
+ }
+
+ bool empty() const
+ {
+ return records_.empty();
}
private:
- std::vector<cell_record> cells_;
+ std::vector<swc_record> records_;
};
struct swc_io_raw
{
- using cell_range_type = cell_record_range_raw;
+ using record_range_type = swc_record_range_raw;
};
struct swc_io_clean
{
- using cell_range_type = cell_record_range_clean;
+ using record_range_type = swc_record_range_clean;
};
template<typename T = swc_io_clean>
- typename T::cell_range_type swc_get_records(std::istream &is)
+typename T::record_range_type swc_get_records(std::istream &is)
{
- return typename T::cell_range_type(is);
+ return typename T::record_range_type(is);
}
-} // end of nestmc::io
-} // end of nestmc
+cell swc_read_cell(std::istream &is);
+
+} // namespace io
+} // namespace mc
+} // namespace nest
diff --git a/src/tree.hpp b/src/tree.hpp
index 6b53aa6f329490a3c431817445f9b2be4dda25b1..d517f7e4570fc9278bab22ff974a3baaa964ea50 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -7,11 +7,15 @@
#include <cassert>
#include "vector/include/Vector.hpp"
+#include "algorithms.hpp"
#include "util.hpp"
-using range = memory::Range;
+namespace nest {
+namespace mc {
class tree {
+ using range = memory::Range;
+
public :
using int_type = int16_t;
@@ -19,30 +23,51 @@ class tree {
using index_type = memory::HostVector<int_type>;
using index_view = index_type::view_type;
+ tree() = default;
+
+ tree& operator=(tree&& other) {
+ std::swap(data_, other.data_);
+ std::swap(child_index_, other.child_index_);
+ std::swap(children_, other.children_);
+ std::swap(parents_, other.parents_);
+ return *this;
+ }
+
+ tree& operator=(tree const& other) {
+ data_ = other.data_;
+ set_ranges(other.num_nodes());
+ return *this;
+ }
+
+ // copy constructors take advantage of the assignment operators
+ // defined above
tree(tree const& other)
- : data_(other.data_)
{
- set_ranges(other.num_nodes());
+ *this = other;
}
tree(tree&& other)
- : data_(std::move(other.data_))
{
- set_ranges(other.num_nodes());
+ *this = std::move(other);
}
- tree() = default;
-
/// create the tree from a parent_index
template <typename I>
- std::vector<I>
- init_from_parent_index(std::vector<I> const& parent_index)
+ tree(std::vector<I> const& parent_index)
{
+ // validate the inputs
+ if(!algorithms::is_minimal_degree(parent_index)) {
+ throw std::domain_error(
+ "parent index used to build a tree did not satisfy minimal degree ordering"
+ );
+ }
+
// n = number of compartment in cell
auto n = parent_index.size();
- // On completion of this loop child_count[i] is the number of children of compartment i
- // compensate count for compartment 0, which has itself as its own parent
+ // On completion of this loop child_count[i] is the number of children
+ // of compartment i compensate count for compartment 0, which has itself
+ // as its own parent
index_type child_count(n, 0);
child_count[0] = -1;
for(auto i : parent_index) {
@@ -123,9 +148,6 @@ class tree {
child_index_[i+1] = child_index_[i];
}
child_index_[0] = 0;
-
- // return the branch index to the caller for later use
- return branch_index;
}
size_t num_children() const {
@@ -135,7 +157,10 @@ class tree {
return child_index_[b+1] - child_index_[b];
}
size_t num_nodes() const {
- return child_index_.size() - 1;
+ // the number of nodes is the size of the child index minus 1
+ // ... except for the case of an empty tree
+ auto sz = child_index_.size();
+ return sz ? sz - 1 : 0;
}
/// return the child index
@@ -219,19 +244,26 @@ class tree {
}
void set_ranges(int nnode) {
- auto nchild = nnode - 1;
- // data_ is partitioned as follows:
- // data_ = [children_[nchild], child_index_[nnode+1], parents_[nnode]]
- assert(data_.size() == unsigned(nchild + (nnode+1) + nnode));
- children_ = data_(0, nchild);
- child_index_ = data_(nchild, nchild+nnode+1);
- parents_ = data_(nchild+nnode+1, memory::end);
-
- // check that arrays have appropriate size
- // this should be moved into a unit test
- assert(children_.size() == unsigned(nchild));
- assert(child_index_.size() == unsigned(nnode+1));
- assert(parents_.size() == unsigned(nnode));
+ if(nnode) {
+ auto nchild = nnode - 1;
+ // data_ is partitioned as follows:
+ // data_ = [children_[nchild], child_index_[nnode+1], parents_[nnode]]
+ assert(data_.size() == unsigned(nchild + (nnode+1) + nnode));
+ children_ = data_(0, nchild);
+ child_index_ = data_(nchild, nchild+nnode+1);
+ parents_ = data_(nchild+nnode+1, memory::end);
+
+ // check that arrays have appropriate size
+ // this should be moved into a unit test
+ assert(children_.size() == unsigned(nchild));
+ assert(child_index_.size() == unsigned(nnode+1));
+ assert(parents_.size() == unsigned(nnode));
+ }
+ else {
+ children_ = data_(0, 0);
+ child_index_ = data_(0, 0);
+ parents_ = data_(0, 0);
+ }
}
/// Renumber the sub-tree with old_node as its root with new_node as
@@ -264,7 +296,6 @@ class tree {
// check for the senitel that indicates that the old root has
// been processed
if(old_node==-1) {
- //assert(parent_node<0); // sanity check
return new_node;
}
@@ -304,15 +335,71 @@ class tree {
}
}
if(add_parent_as_child) {
- new_node = add_children(new_node, old_tree.parent(old_node), old_node, p, old_tree);
+ new_node =
+ add_children(
+ new_node, old_tree.parent(old_node), old_node, p, old_tree
+ );
}
return new_node;
}
+ //////////////////////////////////////////////////
+ // state
+ //////////////////////////////////////////////////
index_type data_;
- index_view children_;
- index_view child_index_;
- index_view parents_;
+ // provide default parameters so that tree type can
+ // be default constructed
+ index_view children_ = data_(0, 0);
+ index_view child_index_= data_(0, 0);
+ index_view parents_ = data_(0, 0);
};
+
+template <typename C>
+std::vector<int> make_parent_index(tree const& t, C const& counts)
+{
+ using range = memory::Range;
+
+ if( !algorithms::is_positive(counts)
+ || counts.size() != t.num_nodes() )
+ {
+ throw std::domain_error(
+ "make_parent_index requires one non-zero count per segment"
+ );
+ }
+ auto index = algorithms::make_index(counts);
+ auto num_compartments = index.back();
+ std::vector<int> parent_index(num_compartments);
+ auto pos = 0;
+ for(int i : range(0, t.num_nodes())) {
+ // get the parent of this segment
+ // taking care for the case where the root node has -1 as its parent
+ auto parent = t.parent(i);
+ parent = parent>=0 ? parent : 0;
+
+ // the index of the first compartment in the segment
+ // is calculated differently for the root (i.e when i==parent)
+ if(i!=parent) {
+ parent_index[pos++] = index[parent+1]-1;
+ }
+ else {
+ parent_index[pos++] = parent;
+ }
+ // number the remaining compartments in the segment consecutively
+ while(pos<index[i+1]) {
+ parent_index[pos] = pos-1;
+ pos++;
+ }
+ }
+
+ // if one of these assertions is tripped, we have to improve
+ // the input validation above
+ assert(pos==num_compartments);
+ assert(algorithms::is_minimal_degree(parent_index));
+
+ return parent_index;
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/src/util.hpp b/src/util.hpp
index 17a3a49b0a19ee867229f71e08ab27b2c6962d57..706ab7da5a17d493823b8c74d1f18fa6b4a9efaf 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -2,6 +2,12 @@
#include "vector/include/Vector.hpp"
+#ifdef DEBUG
+#define EXPECTS(expression) assert(expression)
+#else
+#define EXPECTS(expression)
+#endif
+
/*
using memory::util::red;
using memory::util::yellow;
@@ -11,7 +17,9 @@ using memory::util::blue;
using memory::util::cyan;
*/
+#include <memory>
#include <ostream>
+#include <utility>
#include <vector>
template <typename T>
@@ -20,9 +28,103 @@ operator << (std::ostream &o, std::vector<T>const& v)
{
o << "[";
for(auto const& i: v) {
- o << i << ", ";
+ o << i << ",";
}
o << "]";
return o;
}
+template <typename T>
+std::ostream& print(std::ostream &o, std::vector<T>const& v)
+{
+ o << "[";
+ for(auto const& i: v) {
+ o << i << ",";
+ }
+ o << "]";
+ return o;
+}
+
+namespace nest {
+namespace mc {
+namespace util {
+
+ // just because we aren't using C++14, doesn't mean we shouldn't go
+ // without make_unique
+ template <typename T, typename... Args>
+ std::unique_ptr<T> make_unique(Args&&... args) {
+ return std::unique_ptr<T>(new T(std::forward<Args>(args) ...));
+ }
+
+ /// helper for taking the first value in a std::pair
+ template <typename L, typename R>
+ L const& left(std::pair<L, R> const& p)
+ {
+ return p.first;
+ }
+
+ /// helper for taking the second value in a std::pair
+ template <typename L, typename R>
+ R const& right(std::pair<L, R> const& p)
+ {
+ return p.second;
+ }
+
+ template <typename T>
+ struct is_std_vector : std::false_type {};
+
+ template <typename T, typename C>
+ struct is_std_vector<std::vector<T,C>> : std::true_type {};
+
+ template <typename T>
+ struct is_container :
+ std::conditional<
+ is_std_vector<typename std::decay<T>::type>::value ||
+ memory::is_array<T>::value,
+ std::true_type,
+ 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 7e787ba46b50a48d29ca3d5dd8d70f8116354d5d..cb33af4d972dd398952b1879cca764c502d8f42c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,26 +1,56 @@
set(HEADERS
+ ../src/cell.hpp
+ ../src/cell_tree.hpp
+ ../src/math.hpp
+ ../src/point.hpp
+ ../src/segment.hpp
../src/swcio.hpp
../src/tree.hpp
- ../src/cell_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_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
new file mode 100644
index 0000000000000000000000000000000000000000..f69569c9611068250cc3265524cfe64e1e10bd55
--- /dev/null
+++ b/tests/test_algorithms.cpp
@@ -0,0 +1,482 @@
+#include "gtest.h"
+
+#include <vector>
+
+#include <algorithms.hpp>
+#include <util.hpp>
+
+TEST(algorithms, sum)
+{
+ // sum of 10 times 2 is 20
+ std::vector<int> v1(10, 2);
+ EXPECT_EQ(10*2, nest::mc::algorithms::sum(v1));
+
+ // make an array 1:20 and sum it up using formula for arithmetic sequence
+ std::vector<int> v2(20);
+ std::iota(v2.begin(), v2.end(), 1);
+ auto n = 20;
+ EXPECT_EQ((n+1)*n/2, nest::mc::algorithms::sum(v2));
+}
+
+TEST(algorithms, make_index)
+{
+ {
+ std::vector<int> v(10, 1);
+ auto index = nest::mc::algorithms::make_index(v);
+
+ EXPECT_EQ(index.size(), 11u);
+ EXPECT_EQ(index.front(), 0);
+ EXPECT_EQ(index.back(), nest::mc::algorithms::sum(v));
+ }
+
+ {
+ std::vector<int> v(10);
+ std::iota(v.begin(), v.end(), 1);
+ auto index = nest::mc::algorithms::make_index(v);
+
+ EXPECT_EQ(index.size(), 11u);
+ EXPECT_EQ(index.front(), 0);
+ EXPECT_EQ(index.back(), nest::mc::algorithms::sum(v));
+ }
+}
+
+TEST(algorithms, minimal_degree)
+{
+ {
+ std::vector<int> v = {0};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 0, 1, 2, 3, 4};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 0, 1, 2, 0, 4};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 0, 1, 2, 0, 4, 5, 4};
+ EXPECT_TRUE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {1};
+ EXPECT_FALSE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 2};
+ EXPECT_FALSE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+
+ {
+ std::vector<int> v = {0, 1, 2};
+ EXPECT_FALSE(nest::mc::algorithms::is_minimal_degree(v));
+ }
+}
+
+TEST(algorithms, is_strictly_monotonic_increasing)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{0}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{0, 1, 2, 3}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{8, 20, 42, 89}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{0, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{8, 20, 20, 89}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_increasing(
+ std::vector<int>{3, 2, 1, 0}
+ )
+ );
+}
+
+TEST(algorithms, is_strictly_monotonic_decreasing)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{0}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{3, 2, 1, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{0, 1, 2, 3}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{8, 20, 42, 89}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{0, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_strictly_monotonic_decreasing(
+ std::vector<int>{8, 20, 20, 89}
+ )
+ );
+}
+
+TEST(algorithms, is_positive)
+{
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{}
+ )
+ );
+ EXPECT_TRUE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{3, 2, 1}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{3, 2, 1, 0}
+ )
+ );
+ EXPECT_FALSE(
+ nest::mc::algorithms::is_positive(
+ std::vector<int>{-1}
+ )
+ );
+}
+
+TEST(algorithms, has_contiguous_segments)
+{
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // /|\.
+ // 3 7 4
+ // / \.
+ // 5 6
+ //
+ EXPECT_FALSE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 2, 3, 4, 2}
+ )
+ );
+
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // /|\.
+ // 3 6 5
+ // / \.
+ // 4 7
+ //
+ EXPECT_FALSE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 3, 2, 2, 5}
+ )
+ );
+
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // /|\.
+ // 3 7 5
+ // / \.
+ // 4 6
+ //
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 3, 2, 5, 2}
+ )
+ );
+
+ //
+ // 0
+ // |
+ // 1
+ // / \.
+ // 2 7
+ // / \.
+ // 3 5
+ // / \.
+ // 4 6
+ //
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0, 0, 1, 2, 3, 2, 5, 1}
+ )
+ );
+
+ // Soma-only list
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{0}
+ )
+ );
+
+ // Empty list
+ EXPECT_TRUE(
+ nest::mc::algorithms::has_contiguous_segments(
+ std::vector<int>{}
+ )
+ );
+}
+
+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)
+{
+ {
+ //
+ // 0
+ // /|\.
+ // 1 4 6
+ // / | \.
+ // 2 5 7
+ // / \.
+ // 3 8
+ // / \.
+ // 9 11
+ // / \.
+ // 10 12
+ // \.
+ // 13
+ //
+ std::vector<int> parent_index =
+ { 0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12 };
+ std::vector<int> expected_child_count =
+ { 3, 1, 1, 0, 1, 0, 1, 1, 2, 1, 0, 1, 1, 0 };
+
+ // auto count = nest::mc::algorithms::child_count(parent_index);
+ EXPECT_EQ(expected_child_count,
+ nest::mc::algorithms::child_count(parent_index));
+ }
+
+}
+
+TEST(algorithms, branches)
+{
+ using namespace nest::mc;
+
+ {
+ //
+ // 0
+ // /|\.
+ // 1 4 6
+ // / | \.
+ // 2 5 7
+ // / \.
+ // 3 8
+ // / \.
+ // 9 11
+ // / \.
+ // 10 12
+ // \.
+ // 13
+ //
+ std::vector<int> parent_index =
+ { 0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12 };
+ std::vector<int> expected_branches =
+ { 0, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5 };
+
+ auto actual_branches = algorithms::branches(parent_index);
+ EXPECT_EQ(expected_branches, actual_branches);
+ }
+
+ {
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // |
+ // 3
+ //
+ std::vector<int> parent_index =
+ { 0, 0, 1, 2 };
+ std::vector<int> expected_branches =
+ { 0, 1, 1, 1 };
+
+ auto actual_branches = algorithms::branches(parent_index);
+ EXPECT_EQ(expected_branches, actual_branches);
+ }
+
+ {
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // / \.
+ // 3 4
+ // \.
+ // 5
+ //
+ std::vector<int> parent_index =
+ { 0, 0, 1, 2, 2, 4 };
+ std::vector<int> expected_branches =
+ { 0, 1, 1, 2, 3, 3 };
+
+ auto actual_branches = algorithms::branches(parent_index);
+ EXPECT_EQ(expected_branches, actual_branches);
+ }
+
+ {
+ std::vector<int> parent_index = { 0 };
+ std::vector<int> expected_branches = { 0 };
+
+ 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_cell.cpp b/tests/test_cell.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..88c8dcaced167d4030529adb8a348a116158c17f
--- /dev/null
+++ b/tests/test_cell.cpp
@@ -0,0 +1,165 @@
+#include "gtest.h"
+
+#include "../src/cell.hpp"
+
+TEST(cell_type, soma)
+{
+ // test that insertion of a soma works
+ // define with no centre point
+ {
+ nest::mc::cell c;
+ auto soma_radius = 2.1;
+
+ EXPECT_EQ(c.has_soma(), false);
+ c.add_soma(soma_radius);
+ EXPECT_EQ(c.has_soma(), true);
+
+ auto s = c.soma();
+ EXPECT_EQ(s->radius(), soma_radius);
+ EXPECT_EQ(s->center().is_set(), false);
+ }
+
+ // test that insertion of a soma works
+ // define with centre point @ (0,0,1)
+ {
+ nest::mc::cell c;
+ auto soma_radius = 3.2;
+
+ EXPECT_EQ(c.has_soma(), false);
+ c.add_soma(soma_radius, {0,0,1});
+ EXPECT_EQ(c.has_soma(), true);
+
+ auto s = c.soma();
+ EXPECT_EQ(s->radius(), soma_radius);
+ EXPECT_EQ(s->center().is_set(), true);
+
+ // add expression template library for points
+ //EXPECT_EQ(s->center(), point<double>(0,0,1));
+ }
+}
+
+TEST(cell_type, add_segment)
+{
+ using namespace nest::mc;
+ // add a pre-defined segment
+ {
+ cell c;
+
+ auto soma_radius = 2.1;
+ auto cable_radius = 0.1;
+ auto cable_length = 8.3;
+
+ // add a soma because we need something to attach the first dendrite to
+ c.add_soma(soma_radius, {0,0,1});
+
+ auto seg =
+ make_segment<cable_segment>(
+ segmentKind::dendrite,
+ cable_radius, cable_radius, cable_length
+ );
+ c.add_cable(0, std::move(seg));
+
+ EXPECT_EQ(c.num_segments(), 2);
+ }
+
+ // add segment on the fly
+ {
+ cell c;
+
+ auto soma_radius = 2.1;
+ auto cable_radius = 0.1;
+ auto cable_length = 8.3;
+
+ // add a soma because we need something to attach the first dendrite to
+ c.add_soma(soma_radius, {0,0,1});
+
+ c.add_cable(
+ 0,
+ segmentKind::dendrite, cable_radius, cable_radius, cable_length
+ );
+
+ EXPECT_EQ(c.num_segments(), 2);
+ }
+ {
+ cell c;
+
+ auto soma_radius = 2.1;
+ auto cable_radius = 0.1;
+ auto cable_length = 8.3;
+
+ // add a soma because we need something to attach the first dendrite to
+ c.add_soma(soma_radius, {0,0,1});
+
+ c.add_cable(
+ 0,
+ segmentKind::dendrite,
+ std::vector<double>{cable_radius, cable_radius, cable_radius, cable_radius},
+ std::vector<double>{cable_length, cable_length, cable_length}
+ );
+
+ EXPECT_EQ(c.num_segments(), 2);
+ }
+}
+
+TEST(cell_type, multiple_cables)
+{
+ using namespace nest::mc;
+
+ // generate a cylindrical cable segment of length 1/pi and radius 1
+ // volume = 1
+ // area = 2
+ auto seg = [](segmentKind k) {
+ return make_segment<cable_segment>( k, 1.0, 1.0, 1./math::pi<double>() );
+ };
+
+ // add a pre-defined segment
+ {
+ cell c;
+
+ auto soma_radius = std::pow(3./(4.*math::pi<double>()), 1./3.);
+
+ // cell strucure as follows
+ // left : segment numbering
+ // right : segment type (soma, axon, dendrite)
+ //
+ // 0 s
+ // / \ / \.
+ // 1 2 d a
+ // / \ / \.
+ // 3 4 d d
+
+ // add a soma
+ c.add_soma(soma_radius, {0,0,1});
+
+ // hook the dendrite and axons
+ c.add_cable(0, seg(segmentKind::dendrite));
+ c.add_cable(0, seg(segmentKind::axon));
+ c.add_cable(1, seg(segmentKind::dendrite));
+ c.add_cable(1, seg(segmentKind::dendrite));
+
+ EXPECT_EQ(c.num_segments(), 5);
+ // each of the 5 segments has volume 1 by design
+ EXPECT_EQ(c.volume(), 5.);
+ // each of the 4 cables has volume 2., and the soma has an awkward area
+ // that isn't a round number
+ EXPECT_EQ(c.area(), 8. + math::area_sphere(soma_radius));
+
+ // construct the graph
+ const auto model = c.model();
+ auto const& con = model.tree;
+
+ EXPECT_EQ(con.num_segments(), 5u);
+ EXPECT_EQ(con.parent(0), -1);
+ EXPECT_EQ(con.parent(1), 0);
+ EXPECT_EQ(con.parent(2), 0);
+ EXPECT_EQ(con.parent(3), 1);
+ EXPECT_EQ(con.parent(4), 1);
+ EXPECT_EQ(con.num_children(0), 2u);
+ EXPECT_EQ(con.num_children(1), 2u);
+ EXPECT_EQ(con.num_children(2), 0u);
+ EXPECT_EQ(con.num_children(3), 0u);
+ EXPECT_EQ(con.num_children(4), 0u);
+
+ }
+}
+
diff --git a/tests/test_compartments.cpp b/tests/test_compartments.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c109f9cda9a335ce32f695f85ccd0295446a5d92
--- /dev/null
+++ b/tests/test_compartments.cpp
@@ -0,0 +1,133 @@
+#include <cmath>
+
+#include "gtest.h"
+
+#include "../src/compartment.hpp"
+#include "../src/util.hpp"
+
+using nest::mc::util::left;
+using nest::mc::util::right;
+
+// not much to test here: just test that values passed into the constructor
+// are correctly stored in members
+TEST(compartments, compartment)
+{
+
+ {
+ nest::mc::compartment c(100, 1.2, 2.1, 2.2);
+ EXPECT_EQ(c.index, 100);
+ EXPECT_EQ(c.length, 1.2);
+ EXPECT_EQ(left(c.radius), 2.1);
+ EXPECT_EQ(right(c.radius), 2.2);
+
+ auto c2 = c;
+ EXPECT_EQ(c2.index, 100);
+ EXPECT_EQ(c2.length, 1.2);
+ EXPECT_EQ(left(c2.radius), 2.1);
+ EXPECT_EQ(right(c2.radius), 2.2);
+ }
+
+ {
+ nest::mc::compartment c{100, 1, 2, 3};
+ EXPECT_EQ(c.index, 100);
+ EXPECT_EQ(c.length, 1.);
+ EXPECT_EQ(left(c.radius), 2.);
+ EXPECT_EQ(right(c.radius), 3.);
+ }
+}
+
+TEST(compartments, compartment_iterator)
+{
+ // iterator with arguments
+ // idx = 0
+ // len = 2.5
+ // rad = 1
+ // delta_rad = 2
+ nest::mc::compartment_iterator it(0, 2.5, 1, 2);
+
+ // so each time the iterator is incremented
+ // idx is incremented by 1
+ // len is unchanged
+ // rad is incremented by 2
+
+ // check the prefix increment
+ ++it;
+ {
+ auto c = *it;
+ EXPECT_EQ(c.index, 1);
+ EXPECT_EQ(left(c.radius), 3.0);
+ EXPECT_EQ(right(c.radius), 5.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+
+ // check postfix increment
+
+ // returned iterator should be unchanged
+ {
+ auto c = *(it++);
+ EXPECT_EQ(c.index, 1);
+ EXPECT_EQ(left(c.radius), 3.0);
+ EXPECT_EQ(right(c.radius), 5.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+ // while the iterator itself was updated
+ {
+ auto c = *it;
+ EXPECT_EQ(c.index, 2);
+ EXPECT_EQ(left(c.radius), 5.0);
+ EXPECT_EQ(right(c.radius), 7.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+
+ // check that it can be copied and compared
+ {
+ // copy iterator
+ auto it2 = it;
+ auto c = *it2;
+ EXPECT_EQ(c.index, 2);
+ EXPECT_EQ(left(c.radius), 5.0);
+ EXPECT_EQ(right(c.radius), 7.0);
+ EXPECT_EQ(c.length, 2.5);
+
+ // comparison
+ EXPECT_EQ(it2, it);
+ it2++;
+ EXPECT_NE(it2, it);
+
+ // check the copy has updated correctly when incremented
+ c= *it2;
+ EXPECT_EQ(c.index, 3);
+ EXPECT_EQ(left(c.radius), 7.0);
+ EXPECT_EQ(right(c.radius), 9.0);
+ EXPECT_EQ(c.length, 2.5);
+ }
+}
+
+TEST(compartments, compartment_range)
+{
+ {
+ nest::mc::compartment_range rng(10, 1.0, 2.0, 10.);
+
+ EXPECT_EQ((*rng.begin()).index, 0);
+ EXPECT_EQ((*rng.end()).index, 10);
+ EXPECT_NE(rng.begin(), rng.end());
+
+ auto count = 0;
+ for(auto c : rng) {
+ EXPECT_EQ(c.index, count);
+ auto er = 1.0 + double(count)/10.;
+ EXPECT_TRUE(std::fabs(left(c.radius)-er)<1e-15);
+ EXPECT_TRUE(std::fabs(right(c.radius)-(er+0.1))<1e-15);
+ EXPECT_EQ(c.length, 1.0);
+ ++count;
+ }
+ EXPECT_EQ(count, 10);
+ }
+
+ // test case of zero length range
+ {
+ nest::mc::compartment_range rng(0, 1.0, 1.0, 0.);
+
+ EXPECT_EQ(rng.begin(), rng.end());
+ }
+}
diff --git a/tests/test_fvm.cpp b/tests/test_fvm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..85a25f3579aa25e2e84860b3144294269ff0895a
--- /dev/null
+++ 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
new file mode 100644
index 0000000000000000000000000000000000000000..73663d6815d876e1740da4cee084c46255a74e3f
--- /dev/null
+++ b/tests/test_matrix.cpp
@@ -0,0 +1,121 @@
+#include <numeric>
+
+#include "gtest.h"
+
+#include <matrix.hpp>
+#include <math.hpp>
+#include <vector/include/Vector.hpp>
+
+TEST(matrix, construct_from_parent_only)
+{
+ using matrix_type = nest::mc::matrix<double, int>;
+
+ // pass parent index as a std::vector by reference
+ // the input should not be moved from
+ {
+ std::vector<int> p = {0,0,1};
+ matrix_type m{p};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 3u);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+
+ // pass parent index as a std::vector by rvalue reference
+ // the input should not be moved from
+ {
+ std::vector<int> p = {0,0,1};
+ matrix_type m{std::move(p)};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 3u);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+
+ // pass parent index as a HostVector by reference
+ // expect that the input is not moved from
+ {
+ memory::HostVector<int> p(3, 0);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{p};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 3u);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+ // pass parent index as a HostVector by rvalue reference
+ // expect that the input is moved from
+ {
+ memory::HostVector<int> p(3, 0);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{std::move(p)};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3u);
+ EXPECT_EQ(p.size(), 0u); // 0 implies moved from
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+}
+
+TEST(matrix, solve)
+{
+ using matrix_type = nest::mc::matrix<double, int>;
+
+ // trivial case : 1x1 matrix
+ {
+ matrix_type m{std::vector<int>{0}};
+
+ m.d()(memory::all) = 2;
+ m.l()(memory::all) = -1;
+ m.u()(memory::all) = -1;
+ m.rhs()(memory::all) = 1;
+
+ m.solve();
+
+ EXPECT_EQ(m.rhs()[0], 0.5);
+ }
+ // matrices in the range of 2x2 to 1000x1000
+ {
+ using namespace nest::mc::math;
+ for(auto n : memory::Range(2,1001)) {
+ auto p = std::vector<int>(n);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{p};
+
+ EXPECT_EQ(m.size(), n);
+ EXPECT_EQ(m.num_cells(), 1);
+
+ m.d()(memory::all) = 2;
+ m.l()(memory::all) = -1;
+ m.u()(memory::all) = -1;
+ m.rhs()(memory::all) = 1;
+
+ m.solve();
+
+ auto x = m.rhs();
+ auto err = square(std::fabs(2.*x[0] - x[1] - 1.));
+ for(auto i : memory::Range(1,n-1)) {
+ err += square(std::fabs(2.*x[i] - x[i-1] - x[i+1] - 1.));
+ }
+ err += square(std::fabs(2.*x[n-1] - x[n-2] - 1.));
+
+ 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_point.cpp b/tests/test_point.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b5bfeb9c1783a6c921ca97f3eb6ae4afd535c311
--- /dev/null
+++ b/tests/test_point.cpp
@@ -0,0 +1,110 @@
+#include <limits>
+#include <cmath>
+
+#include "gtest.h"
+
+#include "../src/point.hpp"
+
+using namespace nest::mc;
+
+TEST(point, construction)
+{
+ {
+ // default constructor
+ point<float> p;
+
+ EXPECT_FALSE(p.is_set());
+
+ // expect NaN, which returns false when comparing for equality
+ EXPECT_NE(p.x, p.x);
+ EXPECT_NE(p.y, p.y);
+ EXPECT_NE(p.z, p.z);
+ }
+
+ {
+ // initializer list
+ point<float> p{1, 2, 3};
+ EXPECT_EQ(p.x, 1.);
+ EXPECT_EQ(p.y, 2.);
+ EXPECT_EQ(p.z, 3.);
+ }
+
+ {
+ // explicit call to constructor
+ point<double> p1(1, 2, 3);
+
+ EXPECT_EQ(p1.x, 1.);
+ EXPECT_EQ(p1.y, 2.);
+ EXPECT_EQ(p1.z, 3.);
+
+ // copy constructor
+ auto p2 = p1;
+ EXPECT_EQ(p1, p2);
+ }
+}
+
+TEST(point, constexpr)
+{
+ // perform test using constexpr
+ constexpr point<double> p1(1, 2, 3);
+ constexpr point<double> p2(1, 2, 3);
+
+ constexpr auto p = p1 + p2;
+
+ static_assert(p.x == p1.x+p2.x, "failed x-component");
+ static_assert(p.y == p1.y+p2.y, "failed y-component");
+ static_assert(p.z == p1.z+p2.z, "failed z-component");
+}
+
+TEST(point, addition)
+{
+ // perform test using constexpr
+ point<double> p1(1, 2, 3);
+ point<double> p2(1, 2, 3);
+
+ auto p = p1 + p2;
+
+ EXPECT_EQ(point<double>(2, 4, 6), p);
+}
+
+TEST(point, subtraction)
+{
+ // perform test using constexpr
+ point<double> p1(1, 2, 3);
+ point<double> p2(1, 2, 3);
+
+ auto p = p1 - p2;
+
+ EXPECT_EQ(point<double>(0, 0, 0), p);
+}
+
+TEST(point, scalar_prod)
+{
+ // perform test using constexpr
+ point<double> p1(1, 2, 3);
+
+ auto p = 0.5 * p1;
+
+ EXPECT_EQ(point<double>(0.5, 1.0, 1.5), p);
+}
+
+TEST(point, norm)
+{
+ // don't use constexpr, because the sqrt is not constexpr
+ point<double> p1(1, 1, 1);
+ point<double> p2(1, 2, 3);
+
+ EXPECT_EQ(norm(p1), std::sqrt(3.));
+ EXPECT_EQ(norm(p2), std::sqrt(1.+4.+9.));
+}
+
+TEST(point, dot)
+{
+ // perform test using constexpr
+ constexpr point<double> p1(1, -1, 1);
+ constexpr point<double> p2(1, 2, 3);
+
+ static_assert(dot(p1,p2)==2., "unable to perform constexpr dot product");
+
+ EXPECT_EQ(dot(p1,p2), 2.);
+}
diff --git a/tests/test_segment.cpp b/tests/test_segment.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9c3df625db9ef87110474d35a5ce88f3f4e07f09
--- /dev/null
+++ b/tests/test_segment.cpp
@@ -0,0 +1,98 @@
+#include <limits>
+
+#include "gtest.h"
+
+#include "../src/segment.hpp"
+
+
+TEST(segments, soma)
+{
+ using namespace nest::mc;
+ using nest::mc::math::pi;
+
+ {
+ auto s = make_segment<soma_segment>(1.0);
+
+ EXPECT_EQ(s->volume(), pi<double>()*4./3.);
+ EXPECT_EQ(s->area(), pi<double>()*4.);
+ EXPECT_EQ(s->kind(), segmentKind::soma);
+ }
+
+ {
+ auto s = make_segment<soma_segment>(1.0, point<double>(0., 1., 2.));
+
+ EXPECT_EQ(s->volume(), pi<double>()*4./3.);
+ EXPECT_EQ(s->area(), pi<double>()*4.);
+ EXPECT_EQ(s->kind(), segmentKind::soma);
+ }
+}
+
+TEST(segments, cable)
+{
+ using namespace nest::mc;
+ using nest::mc::math::pi;
+
+ // take advantage of fact that a cable segment with constant radius 1 and
+ // length 1 has volume=1. and area=2
+ auto length = 1./pi<double>();
+ auto radius = 1.;
+
+ // single cylindrical frustrum
+ {
+ auto s = make_segment<cable_segment>(segmentKind::dendrite, radius, radius, length);
+
+ EXPECT_EQ(s->volume(), 1.0);
+ EXPECT_EQ(s->area(), 2.0);
+ EXPECT_EQ(s->kind(), segmentKind::dendrite);
+ }
+
+ // cable made up of three identical cylindrical frustrums
+ {
+ auto s =
+ make_segment<cable_segment>(
+ segmentKind::axon,
+ std::vector<double>{radius, radius, radius, radius},
+ std::vector<double>{length, length, length}
+ );
+
+ EXPECT_EQ(s->volume(), 3.0);
+ EXPECT_EQ(s->area(), 6.0);
+ EXPECT_EQ(s->kind(), segmentKind::axon);
+ }
+}
+
+TEST(segments, cable_positions)
+{
+ using namespace nest::mc;
+ using nest::mc::math::pi;
+
+ // single frustrum of length 1 and radii 1 and 2
+ // the centre of each end are at the origin (0,0,0) and (0,1,0)
+ {
+ auto s =
+ make_segment<cable_segment>(
+ segmentKind::dendrite,
+ 1, 2,
+ point<double>(0,0,0), point<double>(0,1,0)
+ );
+
+ EXPECT_EQ(s->volume(), math::volume_frustrum(1., 1., 2.));
+ EXPECT_EQ(s->area(), math::area_frustrum (1., 1., 2.));
+ EXPECT_EQ(s->kind(), segmentKind::dendrite);
+ }
+
+ // cable made up of three frustrums
+ // that emulate the frustrum from the previous single-frustrum case
+ {
+ auto s =
+ make_segment<cable_segment>(
+ segmentKind::axon,
+ std::vector<double>{1, 1.5, 2},
+ std::vector<point<double>>{ {0,0,0}, {0,0.5,0}, {0,1,0} }
+ );
+
+ EXPECT_EQ(s->volume(), math::volume_frustrum(1., 1., 2.));
+ EXPECT_EQ(s->area(), math::area_frustrum(1., 1., 2.));
+ EXPECT_EQ(s->kind(), segmentKind::axon);
+ }
+}
diff --git a/tests/test_stimulus.cpp b/tests/test_stimulus.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8fd9307b61a780ad93f430338f7d76c941e20bc0
--- /dev/null
+++ b/tests/test_stimulus.cpp
@@ -0,0 +1,33 @@
+#include "gtest.h"
+
+#include "../src/stimulus.hpp"
+
+TEST(stimulus, i_clamp)
+{
+ using namespace nest::mc;
+
+ // stimulus with delay 2, duration 0.5, amplitude 6.0
+ i_clamp stim(2.0, 0.5, 6.0);
+
+ EXPECT_EQ(stim.delay(), 2.0);
+ EXPECT_EQ(stim.duration(), 0.5);
+ EXPECT_EQ(stim.amplitude(), 6.0);
+
+ // test that current only turned on in the half open interval
+ // t \in [2, 2.5)
+ EXPECT_EQ(stim.amplitude(0.0), 0.0);
+ EXPECT_EQ(stim.amplitude(1.0), 0.0);
+ EXPECT_EQ(stim.amplitude(2.0), 6.0);
+ EXPECT_EQ(stim.amplitude(2.4999), 6.0);
+ EXPECT_EQ(stim.amplitude(2.5), 0.0);
+
+ // update: delay 1.0, duration 1.5, amplitude 3.0
+ stim.set_delay(1.0);
+ stim.set_duration(1.5);
+ stim.set_amplitude(3.0);
+
+ EXPECT_EQ(stim.delay(), 1.0);
+ EXPECT_EQ(stim.duration(), 1.5);
+ EXPECT_EQ(stim.amplitude(), 3.0);
+}
+
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 7ecd5fc013d98e66428d7a85cb658360e4ca6f2d..679260d656c9c7964f0eea599951c371886123f4 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -3,15 +3,17 @@
#include <iostream>
#include <fstream>
#include <numeric>
+#include <type_traits>
#include <vector>
#include "gtest.h"
+#include <cell.hpp>
#include <swcio.hpp>
// SWC tests
-void expect_cell_equals(const nestmc::io::cell_record &expected,
- const nestmc::io::cell_record &actual)
+void expect_record_equals(const nest::mc::io::swc_record &expected,
+ const nest::mc::io::swc_record &actual)
{
EXPECT_EQ(expected.id(), actual.id());
EXPECT_EQ(expected.type(), actual.type());
@@ -22,175 +24,200 @@ void expect_cell_equals(const nestmc::io::cell_record &expected,
EXPECT_EQ(expected.parent(), actual.parent());
}
-TEST(cell_record, construction)
+TEST(swc_record, construction)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// force an invalid type
- cell_record::kind invalid_type = static_cast<cell_record::kind>(100);
- EXPECT_THROW(cell_record cell(invalid_type, 7, 1., 1., 1., 1., 5),
+ swc_record::kind invalid_type = static_cast<swc_record::kind>(100);
+ EXPECT_THROW(swc_record record(invalid_type, 7, 1., 1., 1., 1., 5),
std::invalid_argument);
}
{
// invalid id
- EXPECT_THROW(cell_record cell(
- cell_record::custom, -3, 1., 1., 1., 1., 5),
+ EXPECT_THROW(swc_record record(
+ swc_record::kind::custom, -3, 1., 1., 1., 1., 5),
std::invalid_argument);
}
{
// invalid parent id
- EXPECT_THROW(cell_record cell(
- cell_record::custom, 0, 1., 1., 1., 1., -5),
+ EXPECT_THROW(swc_record record(
+ swc_record::kind::custom, 0, 1., 1., 1., 1., -5),
std::invalid_argument);
}
{
// invalid radius
- EXPECT_THROW(cell_record cell(
- cell_record::custom, 0, 1., 1., 1., -1., -1),
+ EXPECT_THROW(swc_record record(
+ swc_record::kind::custom, 0, 1., 1., 1., -1., -1),
std::invalid_argument);
}
{
// parent_id > id
- EXPECT_THROW(cell_record cell(
- cell_record::custom, 0, 1., 1., 1., 1., 2),
+ EXPECT_THROW(swc_record record(
+ swc_record::kind::custom, 0, 1., 1., 1., 1., 2),
std::invalid_argument);
}
{
// parent_id == id
- EXPECT_THROW(cell_record cell(
- cell_record::custom, 0, 1., 1., 1., 1., 0),
+ EXPECT_THROW(swc_record record(
+ swc_record::kind::custom, 0, 1., 1., 1., 1., 0),
std::invalid_argument);
}
{
// check standard construction by value
- cell_record cell(cell_record::custom, 0, 1., 1., 1., 1., -1);
- EXPECT_EQ(cell.id(), 0);
- EXPECT_EQ(cell.type(), cell_record::custom);
- EXPECT_EQ(cell.x(), 1.);
- EXPECT_EQ(cell.y(), 1.);
- EXPECT_EQ(cell.z(), 1.);
- EXPECT_EQ(cell.radius(), 1.);
- EXPECT_EQ(cell.diameter(), 2*1.);
- EXPECT_EQ(cell.parent(), -1);
+ swc_record record(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
+ EXPECT_EQ(record.id(), 0);
+ EXPECT_EQ(record.type(), swc_record::kind::custom);
+ EXPECT_EQ(record.x(), 1.);
+ EXPECT_EQ(record.y(), 1.);
+ EXPECT_EQ(record.z(), 1.);
+ EXPECT_EQ(record.radius(), 1.);
+ EXPECT_EQ(record.diameter(), 2*1.);
+ EXPECT_EQ(record.parent(), -1);
}
{
// check copy constructor
- cell_record cell_orig(cell_record::custom, 0, 1., 1., 1., 1., -1);
- cell_record cell(cell_orig);
- expect_cell_equals(cell_orig, cell);
+ swc_record record_orig(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
+ swc_record record(record_orig);
+ expect_record_equals(record_orig, record);
}
}
-TEST(cell_record, comparison)
+TEST(swc_record, comparison)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// check comparison operators
- cell_record cell0(cell_record::custom, 0, 1., 1., 1., 1., -1);
- cell_record cell1(cell_record::custom, 0, 2., 3., 4., 5., -1);
- cell_record cell2(cell_record::custom, 1, 2., 3., 4., 5., -1);
- EXPECT_EQ(cell0, cell1);
- EXPECT_LT(cell0, cell2);
- EXPECT_GT(cell2, cell1);
+ swc_record record0(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
+ swc_record record1(swc_record::kind::custom, 0, 2., 3., 4., 5., -1);
+ swc_record record2(swc_record::kind::custom, 1, 2., 3., 4., 5., -1);
+ EXPECT_EQ(record0, record1);
+ EXPECT_LT(record0, record2);
+ EXPECT_GT(record2, record1);
}
}
TEST(swc_parser, invalid_input)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// check incomplete lines; missing parent
std::istringstream is("1 1 14.566132 34.873772 7.857000 0.717830\n");
- cell_record cell;
- EXPECT_THROW(is >> cell, swc_parse_error);
+ swc_record record;
+ EXPECT_THROW(is >> record, swc_parse_error);
}
{
// Check non-parsable values
std::istringstream is(
"1a 1 14.566132 34.873772 7.857000 0.717830 -1\n");
- cell_record cell;
- EXPECT_THROW(is >> cell, swc_parse_error);
+ swc_record record;
+ EXPECT_THROW(is >> record, swc_parse_error);
}
{
- // Check invalid cell type
+ // Check invalid record type
std::istringstream is(
"1 10 14.566132 34.873772 7.857000 0.717830 -1\n");
- cell_record cell;
- EXPECT_THROW(is >> cell, swc_parse_error);
+ swc_record record;
+ EXPECT_THROW(is >> record, swc_parse_error);
+ }
+
+ {
+ // Non-contiguous numbering in branches is considered invalid
+ // 1
+ // / \.
+ // 2 3
+ // /
+ // 4
+ std::stringstream is;
+ is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
+ is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "3 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 1 14.566132 34.873772 7.857000 0.717830 2\n";
+
+ std::vector<swc_record> records;
+ try {
+ for (auto c : swc_get_records<swc_io_clean>(is)) {
+ records.push_back(std::move(c));
+ }
+
+ FAIL() << "expected swc_parse_error, none was thrown\n";
+ } catch (const swc_parse_error &e) {
+ SUCCEED();
+ }
}
}
TEST(swc_parser, valid_input)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
- // check empty file; no cell may be parsed
- cell_record cell, cell_orig;
+ // check empty file; no record may be parsed
+ swc_record record, record_orig;
std::istringstream is("");
- EXPECT_NO_THROW(is >> cell);
- expect_cell_equals(cell_orig, cell);
+ EXPECT_NO_THROW(is >> record);
+ expect_record_equals(record_orig, record);
}
{
// check comment-only file not ending with a newline;
- // no cell may be parsed
- cell_record cell, cell_orig;
+ // no record may be parsed
+ swc_record record, record_orig;
std::istringstream is("#comment\n#comment");
- EXPECT_NO_THROW(is >> cell);
- expect_cell_equals(cell_orig, cell);
+ EXPECT_NO_THROW(is >> record);
+ expect_record_equals(record_orig, record);
}
{
// check last line case (no newline at the end)
std::istringstream is("1 1 14.566132 34.873772 7.857000 0.717830 -1");
- cell_record cell;
- EXPECT_NO_THROW(is >> cell);
- EXPECT_EQ(0, cell.id()); // zero-based indexing
- EXPECT_EQ(cell_record::soma, cell.type());
- EXPECT_FLOAT_EQ(14.566132, cell.x());
- EXPECT_FLOAT_EQ(34.873772, cell.y());
- EXPECT_FLOAT_EQ( 7.857000, cell.z());
- EXPECT_FLOAT_EQ( 0.717830, cell.radius());
- EXPECT_FLOAT_EQ( -1, cell.parent());
+ swc_record record;
+ EXPECT_NO_THROW(is >> record);
+ EXPECT_EQ(0, record.id()); // zero-based indexing
+ EXPECT_EQ(swc_record::kind::soma, record.type());
+ EXPECT_FLOAT_EQ(14.566132, record.x());
+ EXPECT_FLOAT_EQ(34.873772, record.y());
+ EXPECT_FLOAT_EQ( 7.857000, record.z());
+ EXPECT_FLOAT_EQ( 0.717830, record.radius());
+ EXPECT_FLOAT_EQ( -1, record.parent());
}
{
// check valid input with a series of records
- std::vector<cell_record> cells_orig = {
- cell_record(cell_record::soma, 0,
+ std::vector<swc_record> records_orig = {
+ swc_record(swc_record::kind::soma, 0,
14.566132, 34.873772, 7.857000, 0.717830, -1),
- cell_record(cell_record::dendrite, 1,
+ swc_record(swc_record::kind::dendrite, 1,
14.566132+1, 34.873772+1, 7.857000+1, 0.717830+1, -1)
};
std::stringstream swc_input;
swc_input << "# this is a comment\n";
swc_input << "# this is a comment\n";
- for (auto c : cells_orig)
+ for (auto c : records_orig)
swc_input << c << "\n";
swc_input << "# this is a final comment\n";
std::size_t nr_records = 0;
- for (auto cell : swc_get_records<swc_io_raw>(swc_input)) {
- ASSERT_LT(nr_records, cells_orig.size());
- expect_cell_equals(cells_orig[nr_records], cell);
+ for (auto record : swc_get_records<swc_io_raw>(swc_input)) {
+ ASSERT_LT(nr_records, records_orig.size());
+ expect_record_equals(records_orig[nr_records], record);
++nr_records;
}
}
@@ -198,7 +225,7 @@ TEST(swc_parser, valid_input)
TEST(swc_parser, from_allen_db)
{
- using namespace nestmc;
+ using namespace nest::mc;
auto fname = "../data/example.swc";
std::ifstream fid(fname);
@@ -207,8 +234,8 @@ TEST(swc_parser, from_allen_db)
return;
}
- // load the cell records into a std::vector
- std::vector<io::cell_record> nodes;
+ // load the record records into a std::vector
+ std::vector<io::swc_record> nodes;
for (auto node : io::swc_get_records<io::swc_io_raw>(fid)) {
nodes.push_back(std::move(node));
}
@@ -219,43 +246,43 @@ TEST(swc_parser, from_allen_db)
TEST(swc_parser, input_cleaning)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// Check duplicates
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
- EXPECT_EQ(2u, swc_get_records(is).size());
+ EXPECT_EQ(2u, swc_get_records<swc_io_clean>(is).size());
}
{
// Check multiple trees
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "3 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
- auto cells = swc_get_records(is);
- EXPECT_EQ(2u, cells.size());
+ auto records = swc_get_records<swc_io_clean>(is);
+ EXPECT_EQ(2u, records.size());
}
{
// Check unsorted input
std::stringstream is;
- is << "3 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "3 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- std::array<cell_record::id_type, 4> expected_id_list = {{ 0, 1, 2, 3 }};
+ std::array<swc_record::id_type, 4> expected_id_list = {{ 0, 1, 2, 3 }};
auto expected_id = expected_id_list.cbegin();
- for (auto c : swc_get_records(is)) {
+ for (auto c : swc_get_records<swc_io_clean>(is)) {
EXPECT_EQ(*expected_id, c.id());
++expected_id;
}
@@ -268,20 +295,20 @@ TEST(swc_parser, input_cleaning)
// Check holes in numbering
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "21 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "31 1 14.566132 34.873772 7.857000 0.717830 21\n";
- is << "41 1 14.566132 34.873772 7.857000 0.717830 21\n";
- is << "51 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "61 1 14.566132 34.873772 7.857000 0.717830 51\n";
+ is << "21 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "31 2 14.566132 34.873772 7.857000 0.717830 21\n";
+ is << "41 2 14.566132 34.873772 7.857000 0.717830 21\n";
+ is << "51 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "61 2 14.566132 34.873772 7.857000 0.717830 51\n";
- std::array<cell_record::id_type, 6> expected_id_list =
+ std::array<swc_record::id_type, 6> expected_id_list =
{{ 0, 1, 2, 3, 4, 5 }};
- std::array<cell_record::id_type, 6> expected_parent_list =
+ std::array<swc_record::id_type, 6> expected_parent_list =
{{ -1, 0, 1, 1, 0, 4 }};
auto expected_id = expected_id_list.cbegin();
auto expected_parent = expected_parent_list.cbegin();
- for (auto c : swc_get_records(is)) {
+ for (auto c : swc_get_records<swc_io_clean>(is)) {
EXPECT_EQ(*expected_id, c.id());
EXPECT_EQ(*expected_parent, c.parent());
++expected_id;
@@ -294,29 +321,29 @@ TEST(swc_parser, input_cleaning)
}
}
-TEST(cell_record_ranges, raw)
+TEST(swc_record_ranges, raw)
{
- using namespace nestmc::io;
+ using namespace nest::mc::io;
{
// Check valid usage
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "3 1 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "3 2 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
- std::vector<cell_record> cells;
+ std::vector<swc_record> records;
for (auto c : swc_get_records<swc_io_raw>(is)) {
- cells.push_back(c);
+ records.push_back(c);
}
- EXPECT_EQ(4u, cells.size());
+ EXPECT_EQ(4u, records.size());
bool entered = false;
- auto citer = cells.begin();
+ auto citer = records.begin();
for (auto c : swc_get_records<swc_io_raw>(is)) {
- expect_cell_equals(c, *citer++);
+ expect_record_equals(c, *citer++);
entered = true;
}
@@ -343,7 +370,7 @@ TEST(cell_record_ranges, raw)
auto iter = swc_get_records<swc_io_raw>(is).begin();
auto iend = swc_get_records<swc_io_raw>(is).end();
- cell_record c;
+ swc_record c;
EXPECT_NO_THROW(c = *iter++);
EXPECT_EQ(-1, c.parent());
EXPECT_EQ(iend, iter);
@@ -356,14 +383,14 @@ TEST(cell_record_ranges, raw)
// Check parse error context
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- is << "2 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "3 10 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "4 1 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
- std::vector<cell_record> cells;
+ std::vector<swc_record> records;
try {
for (auto c : swc_get_records<swc_io_raw>(is)) {
- cells.push_back(c);
+ records.push_back(c);
}
ADD_FAILURE() << "expected an exception\n";
@@ -371,4 +398,92 @@ TEST(cell_record_ranges, raw)
EXPECT_EQ(3u, e.lineno());
}
}
+
+ {
+ // Test empty range
+ std::stringstream is("");
+ EXPECT_TRUE(swc_get_records<swc_io_raw>(is).empty());
+ EXPECT_TRUE(swc_get_records<swc_io_clean>(is).empty());
+ }
+}
+
+TEST(swc_io, cell_construction)
+{
+ using namespace nest::mc;
+
+ {
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // / \.
+ // 3 4
+ // \.
+ // 5
+ //
+
+ std::stringstream is;
+ is << "1 1 0 0 0 2.1 -1\n";
+ is << "2 2 0.1 1.2 1.2 1.3 1\n";
+ is << "3 2 1.0 2.0 2.2 1.1 2\n";
+ is << "4 2 1.5 3.3 1.3 2.2 3\n";
+ is << "5 2 2.5 5.3 2.5 0.7 3\n";
+ is << "6 2 3.5 2.3 3.7 3.4 5\n";
+
+ using point_type = point<double>;
+ std::vector<point_type> points = {
+ { 0.0, 0.0, 0.0 },
+ { 0.1, 1.2, 1.2 },
+ { 1.0, 2.0, 2.2 },
+ { 1.5, 3.3, 1.3 },
+ { 2.5, 5.3, 2.5 },
+ { 3.5, 2.3, 3.7 },
+ };
+
+ cell cell = io::swc_read_cell(is);
+ EXPECT_TRUE(cell.has_soma());
+ EXPECT_EQ(4, cell.num_segments());
+
+ EXPECT_EQ(norm(points[1]-points[2]), cell.cable(1)->length());
+ EXPECT_EQ(norm(points[2]-points[3]), cell.cable(2)->length());
+ EXPECT_EQ(norm(points[2]-points[4]) + norm(points[4]-points[5]),
+ cell.cable(3)->length());
+
+
+ // Check each segment separately
+ EXPECT_TRUE(cell.segment(0)->is_soma());
+ EXPECT_EQ(2.1, cell.soma()->radius());
+ EXPECT_EQ(point_type(0, 0, 0), cell.soma()->center());
+
+ for (auto i = 1; i < cell.num_segments(); ++i) {
+ EXPECT_TRUE(cell.segment(i)->is_dendrite());
+ }
+
+ EXPECT_EQ(1, cell.cable(1)->num_sub_segments());
+ EXPECT_EQ(1, cell.cable(2)->num_sub_segments());
+ EXPECT_EQ(2, cell.cable(3)->num_sub_segments());
+
+
+ // Check the radii
+ EXPECT_EQ(1.3, cell.cable(1)->radius(0));
+ EXPECT_EQ(1.1, cell.cable(1)->radius(1));
+
+ EXPECT_EQ(1.1, cell.cable(2)->radius(0));
+ EXPECT_EQ(2.2, cell.cable(2)->radius(1));
+
+ EXPECT_EQ(1.1, cell.cable(3)->radius(0));
+ EXPECT_EQ(3.4, cell.cable(3)->radius(1));
+
+ auto len_ratio = norm(points[2]-points[4]) / cell.cable(3)->length();
+ EXPECT_NEAR(.7, cell.cable(3)->radius(len_ratio), 1e-15);
+
+ // Double-check radii at joins are equal
+ EXPECT_EQ(cell.cable(1)->radius(1),
+ cell.cable(2)->radius(0));
+
+ EXPECT_EQ(cell.cable(1)->radius(1),
+ cell.cable(3)->radius(0));
+ }
}
diff --git a/tests/test_tree.cpp b/tests/test_tree.cpp
index 58435e15240344665921799dbf07ff8eca211f5c..a468d80ca8d4a61804d4fd0e242cad466bff8f8c 100644
--- a/tests/test_tree.cpp
+++ b/tests/test_tree.cpp
@@ -11,6 +11,8 @@
using json = nlohmann::json;
using range = memory::Range;
+using namespace nest::mc;
+
TEST(cell_tree, from_parent_index) {
// tree with single branch corresponding to the root node
// this is equivalent to a single compartment model
@@ -28,8 +30,17 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_segments(), 1u);
EXPECT_EQ(tree.num_children(0), 0u);
}
- // tree with two segments off the root node
+
{
+ //
+ // 0 0
+ // / \ / \.
+ // 1 4 => 1 2
+ // / \.
+ // 2 5
+ // /
+ // 3
+ //
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4};
cell_tree tree(parent_index);
@@ -41,9 +52,18 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(2), 0u);
}
{
- // tree with three segments off the root node
+ //
+ // 0 0
+ // /|\ /|\.
+ // 1 4 6 => 1 2 3
+ // / | \.
+ // 2 5 7
+ // / \.
+ // 3 8
+ //
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8};
+
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 4u);
// the root has 3 children
@@ -52,9 +72,29 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(1), 0u);
EXPECT_EQ(tree.num_children(2), 0u);
EXPECT_EQ(tree.num_children(3), 0u);
+
+ // Check new structure
+ EXPECT_EQ(-1, tree.parent(0));
+ EXPECT_EQ(0, tree.parent(1));
+ EXPECT_EQ(0, tree.parent(2));
+ EXPECT_EQ(0, tree.parent(3));
}
{
- // tree with three segments off the root node, and another 2 segments off of the third branch from the root node
+ //
+ // 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 =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12};
cell_tree tree(parent_index);
@@ -68,6 +108,14 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(2), 0u);
EXPECT_EQ(tree.num_children(4), 0u);
EXPECT_EQ(tree.num_children(5), 0u);
+
+ // Check new structure
+ EXPECT_EQ(-1, tree.parent(0));
+ EXPECT_EQ(0, tree.parent(1));
+ EXPECT_EQ(0, tree.parent(2));
+ EXPECT_EQ(0, tree.parent(3));
+ EXPECT_EQ(3, tree.parent(4));
+ EXPECT_EQ(3, tree.parent(5));
}
{
//
@@ -137,8 +185,7 @@ TEST(tree, change_root) {
// |
// 2
std::vector<int> parent_index = {0,0,0};
- tree t;
- t.init_from_parent_index(parent_index);
+ tree t(parent_index);
t.change_root(1);
EXPECT_EQ(t.num_nodes(), 3u);
@@ -156,8 +203,7 @@ TEST(tree, change_root) {
// / \ |
// 3 4 4
std::vector<int> parent_index = {0,0,0,1,1};
- tree t;
- t.init_from_parent_index(parent_index);
+ tree t(parent_index);
t.change_root(1u);
EXPECT_EQ(t.num_nodes(), 5u);
@@ -181,8 +227,7 @@ TEST(tree, change_root) {
// / \.
// 5 6
std::vector<int> parent_index = {0,0,0,1,1,4,4};
- tree t;
- t.init_from_parent_index(parent_index);
+ tree t(parent_index);
t.change_root(1);
@@ -233,21 +278,89 @@ 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");
}
}
// this test doesn't test anything yet... it just loads each cell in turn
// from a json file and creates a .dot file for it
-TEST(cell_tree, json_load) {
+TEST(cell_tree, json_load)
+{
json cell_data;
std::ifstream("../data/cells_small.json") >> cell_data;
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");
}
}
+TEST(tree, make_parent_index)
+{
+ // just the soma
+ {
+ std::vector<int> parent_index = {0};
+ std::vector<int> counts = {1};
+ nest::mc::tree t(parent_index);
+ auto new_parent_index = make_parent_index(t, counts);
+ EXPECT_EQ(parent_index.size(), new_parent_index.size());
+ }
+ // just a soma with 5 compartments
+ {
+ std::vector<int> parent_index = {0};
+ std::vector<int> counts = {5};
+ nest::mc::tree t(parent_index);
+ auto new_parent_index = make_parent_index(t, counts);
+ EXPECT_EQ(new_parent_index.size(), (unsigned)counts[0]);
+ EXPECT_EQ(new_parent_index[0], 0);
+ for(auto i=1u; i<new_parent_index.size(); ++i) {
+ EXPECT_EQ((unsigned)new_parent_index[i], i-1);
+ }
+ }
+ // some trees with single compartment per segment
+ {
+ auto trees = {
+ // 0
+ // |
+ // 1
+ std::vector<int>{0,0},
+ // 0
+ // / \.
+ // 1 2
+ std::vector<int>{0,0,0},
+ // 0
+ // / \.
+ // 1 4
+ // / \ |\.
+ // 2 3 5 6
+ std::vector<int>{0,0,0,1,1,2,2}
+ };
+ for(auto &parent_index : trees) {
+ std::vector<int> counts(parent_index.size(), 1);
+ nest::mc::tree t(parent_index);
+ auto new_parent_index = make_parent_index(t, counts);
+ EXPECT_EQ(parent_index, new_parent_index);
+ }
+ }
+ // a tree with multiple compartments per segment
+ //
+ // 0
+ // / \.
+ // 1 8
+ // / \.
+ // 2 9
+ // /.
+ // 3
+ // / \.
+ // 4 6
+ // / \.
+ // 5 7
+ {
+ std::vector<int> parent_index = {0,0,1,2,3,4,3,6,0,8};
+ std::vector<int> counts = {1,3,2,2,2};
+ nest::mc::tree t(parent_index);
+ auto new_parent_index = make_parent_index(t, counts);
+ EXPECT_EQ(parent_index, new_parent_index);
+ }
+}
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 7fc001ef98c01383d384541de08b861328c29405..8153d030abba20ddd174da9804406139c393a808 160000
--- a/vector
+++ b/vector
@@ -1 +1 @@
-Subproject commit 7fc001ef98c01383d384541de08b861328c29405
+Subproject commit 8153d030abba20ddd174da9804406139c393a808