diff --git a/README.md b/README.md
index 7c76215207038b458f0511ba145680f36dd1d960..ba282f782b70e02f34538839fb79bb16ab047df8 100644
--- a/README.md
+++ b/README.md
@@ -3,16 +3,22 @@
This is the repository for the NestMC prototype code. Unfortunately we do not have thorough documentation of how-to guides.
Below are some guides for how to build the project and run the miniapp.
Contact us or submit a ticket if you have any questions or want help.
-
+https://github.com/eth-cscs/nestmc-proto
+
+1. Basic installation
+2. MPI
+3. TBB
+4. TBB on Cray systems
+5. Targeting KNL
+6. Examples of environment configuration
+ - Julia
+
+## Basic installation
```bash
-# clone repo
+# clone repository
git clone git@github.com:eth-cscs/nestmc-proto.git
cd nestmc-proto/
-# setup sub modules
-git submodule init
-git submodule update
-
# setup environment
# on a desktop system this is probably not required
# on a cluster this is usually required to make sure that an appropriate
@@ -25,7 +31,7 @@ export CXX=`which g++`
# build main project (out-of-tree)
mkdir build
cd build
-cmake ..
+cmake <path to CMakeLists.txt>
make -j
# test
@@ -42,7 +48,6 @@ To ensure that CMake detects MPI correctly, you should specify the MPI wrapper f
export CXX=mpicxx
export CC=mpicc
cmake <path to CMakeLists.txt> -DWITH_MPI=ON
-
```
## TBB
@@ -81,61 +86,55 @@ cmake <path to CMakeLists.txt> -DWITH_TBB=ON -DSYSTEM_CRAY=ON
# multithreading and MPI
cmake <path to CMakeLists.txt> -DWITH_TBB=ON -DWITH_MPI=ON -DSYSTEM_CRAY=ON
-
```
-# targetting KNL
+## targeting KNL
-## build modparser
+#### build modparser without KNL environment
The source to source compiler "modparser" that generates the C++/CUDA kernels for the ion channels and synapses is in a separate repository.
-It is included in our project as a git submodule, and by default it will be built with the same compiler and flags that are used to build the miniapp and tests.
+By default it will be built with the same compiler and flags that are used to build the miniapp and tests.
This can cause problems if we are cross compiling, e.g. for KNL, because the modparser compiler might not be runnable on the compilation node.
-CMake will look for the source to source compiler executable, `modcc`, in the `PATH` environment variable, and will use the version if finds instead of building its own.
+You are probably best of building the software twice: Once without KNL support to create the modcc parser and next the KNL version using
+the now compiled executable
Modparser requires a C++11 compiler, and has been tested on GCC, Intel, and Clang compilers
- if the default compiler on your is some ancient version of gcc you might need to load a module/set the CC and CXX environment variables.
+CMake will look for the source to source compiler executable, `modcc`, in the `PATH` environment variable, and will use the version if finds instead of building its own.
+So add the g++ compiled modcc to your path
+e.g:
+
```bash
-git clone git@github.com:eth-cscs/modparser.git
-cd modparser
+# First build a 'normal' non KNL version of the software
-# example of setting a C++11 compiler
-export CXX=`which gcc-4.8`
+# Load your environment (see section 6 for detailed example)
+export CC=`which gcc`; export CXX=`which g++`
-cmake .
-make -j
+# Make directory , do the configuration and build
+mkdir build
+cd build
+cmake <path to CMakeLists.txt> -DCMAKE_BUILD_TYPE=release
+make -j8
# set path and test that you can see modcc
export PATH=`pwd`/bin:$PATH
which modcc
```
-## set up environment
+#### set up environment
- source the intel compilers
- source the TBB vars
-- I have only tested with the latest stable version from online, not the version that comes installed sometimes with the Intel compilers.
+- I have only tested with the latest stable version from on-line, not the version that comes installed sometimes with the Intel compilers.
-## build miniapp
+#### build miniapp
```bash
-# clone the repo and set up the submodules
+# clone the repository and set up the submodules
git clone https://github.com/eth-cscs/nestmc-proto.git
cd nestmc-proto
-git submodule init
-git submodule update
-
-# make a path for out of source build
-mkdir build_knl
-cd build_knl
-
-## build miniapp
-
-# setup submodules
-git submodule init
-git submodule update
# make a path for out of source build
mkdir build_knl
@@ -144,13 +143,13 @@ cd build_knl
# run cmake with all the magic flags
export CC=`which icc`
export CXX=`which icpc`
-cmake .. -DCMAKE_BUILD_TYPE=release -DWITH_TBB=ON -DWITH_PROFILING=ON -DVECTORIZE_TARGET=KNL -DUSE_OPTIMIZED_KERNELS=ON
+cmake <path to CMakeLists.txt> -DCMAKE_BUILD_TYPE=release -DWITH_TBB=ON -DWITH_PROFILING=ON -DVECTORIZE_TARGET=KNL -DUSE_OPTIMIZED_KERNELS=ON
make -j
```
The flags passed into cmake are described:
- `-DCMAKE_BUILD_TYPE=release` : build in release mode with `-O3`.
- - `-WITH_TBB=ON` : use TBB for threading on multicore
+ - `-WITH_TBB=ON` : use TBB for threading on multi-core
- `-DWITH_PROFILING=ON` : use internal profilers that print profiling report at end
- `-DVECTORIZE_TARGET=KNL` : generate AVX512 instructions, alternatively you can use:
- `AVX2` for Haswell & Broadwell
@@ -158,7 +157,7 @@ The flags passed into cmake are described:
- `-DUSE_OPTIMIZED_KERNELS=ON` : tell the source to source compiler to generate optimized kernels that use Intel extensions
- without these vectorized code will not be generated.
-## run tests
+#### run tests
Run some unit tests
```bash
@@ -187,7 +186,7 @@ cd miniapp
# a small run to check that everything works
./miniapp.exe -n 1000 -s 200 -t 20 -p 0
-# a larger run for generating meaninful benchmarks
+# a larger run for generating meaningful benchmarks
./miniapp.exe -n 2000 -s 2000 -t 100 -p 0
```
@@ -220,3 +219,12 @@ total | 0.791 100.0 | 38.593 100.0 |
-----------------------------------------------------
```
+## Examples of environment configuration
+### Julia (HBP PCP system)
+``` bash
+module load cmake
+module load intel-ics
+module load openmpi_ics/2.0.0
+module load gcc/6.1.0
+```
+
diff --git a/docs/passive_cable/cable.bib b/docs/passive_cable/cable.bib
new file mode 100644
index 0000000000000000000000000000000000000000..1f47e75b738ede826221dbc4669e419c21bfaa30
--- /dev/null
+++ b/docs/passive_cable/cable.bib
@@ -0,0 +1,33 @@
+@article{bhalla1992,
+ title = {Rallpacks: a set of benchmarks for neuronal simulators},
+ author = {Upinder S. Bhalla and David H. Bilitch and James M. Bower},
+ journal = {Trends in Neurosciences},
+ volume = {15},
+ number = {11},
+ pages = {453--458},
+ year = {1992},
+ month = {11},
+ doi = {10.1016/0166-2236(92)90009-W}
+}
+
+@article{churchill1937,
+ title = {The inversion of the Laplace transformation by a direct expansion in series and its application to boundary-value problems},
+ author = {R. V. Churchill},
+ journal = {Mathematische Zeitschrift},
+ volume = {42},
+ number = {1},
+ pages = {567--579},
+ year = {1937},
+ doi = {10.1007/BF01160095}
+}
+
+@article{borwein2009,
+ title = {Uniform bounds for the complementary incomplete gamma function},
+ author = {Jonathan M Borwein and O-Yeat Chan},
+ journal = {Mathematical Inequalities and Applications},
+ volume = {12},
+ number = {1},
+ pages = {115--121},
+ year = {2009},
+ doi = {10.7153/mia-12-10}
+}
diff --git a/docs/passive_cable/cable_computation.tex b/docs/passive_cable/cable_computation.tex
new file mode 100644
index 0000000000000000000000000000000000000000..39a372196c006c7c67837f039df88c903cb439b4
--- /dev/null
+++ b/docs/passive_cable/cable_computation.tex
@@ -0,0 +1,324 @@
+\documentclass[parskip=half]{scrartcl}
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage[citestyle=authoryear-icomp,bibstyle=authoryear]{biblatex}
+\usepackage{booktabs}
+\usepackage[style=british]{csquotes}
+\usepackage{fontspec}
+\usepackage{siunitx}
+
+\setmainfont[Numbers=Lowercase]{TeX Gyre Pagella}
+\newfontfamily\dispfamily{TeX Gyre Pagella}
+\addtokomafont{disposition}{\dispfamily}
+
+\MakeAutoQuote{«}{»}
+\DeclareMathOperator{\Res}{Res}
+
+\addbibresource{cable.bib}
+
+\title{Analytic solution to passive cable model}
+\author{Sam Yates}
+\date{October 20, 2016}
+
+\begin{document}
+
+\maketitle
+
+\section{Background}
+
+The Rallpack suite \parencite{bhalla1992} is a set of three tests for cable-based neuronal
+simulators, for validation and benchmarking. Validation data is provided with the
+distribution of the suite, which at the time of writing can now be found within the
+\textsc{Genesis} simulator simulator distribution.\footnote{see \url{http://genesis-sim.org}.}
+
+The first two of the tests are based on passive cable models which admit analytic
+solutions, and the Rallpack suite includes code that can generate the reference curves
+for these models. The reference code, however, requires hand-tuning of the iterations;
+the authors state in the internal documentation that
+«the use of 10 terms gives at least six figure accuracy», but this is for the
+fixed time increments used in their dataset and is not validated within the code
+or by the supplied material.
+
+For flexibility of testing and cross-checking of the Rallpack reference data,
+it would be useful to have at hand a robust calculator of the passive cable
+potential solution.
+
+\section{The Rallpack 1 model}
+
+The model comprises a cylindrical cable segment with the physical and electrical
+properties as described in Table~\ref{tbl:rallpack1}. A current of \SI{0.1}{\nA}
+is injected at the left end of the cable from $t=0$, and the initial potential
+is the reversal potential.
+
+\begin{table}[ht]
+ \centering
+ \begin{tabular}{lSl}
+ \toprule
+ Term & {Value} & Property\\
+ \midrule
+ $d$ & \SI{1.0}{\um} & cable diameter \\
+ $L$ & \SI{1.0}{\mm} & cable length \\
+ $R_A$ & \SI{1.0}{\ohm\m} & bulk axial resistivity \\
+ $R_M$ & \SI{4.0}{\ohm\m\squared} & areal membrane resistivity \\
+ $C_M$ & \SI{0.01}{\F\per\m\squared} & areal membrane capacitance \\
+ $E_M$ & \SI{-65.0}{\mV} & membrane reversal potential \\
+ \bottomrule
+ \end{tabular}
+ \caption{Cable properties for the Rallpack 1 model.}
+ \label{tbl:rallpack1}
+\end{table}
+
+The potential on a constant-radius passive cable is governed by the PDE
+\begin{equation}
+ \lambda^2 \frac{\partial^2 v}{\partial x^2} =
+ \tau\frac{\partial v}{\partial t} + v - E,
+\end{equation}
+where $E$ is the membrane reversal potential, and $\lambda$ and $\tau$ are the
+electrotonic length and time constants for the cable. It is convenient
+to consider the electrical properties of the cable per unit-length, in terms
+of the linear axial resistivity $r$, the linear membrane capacitance $r$,
+and the linear membrane capacitance $c$. These determine $\lambda$ and $\tau$ by
+\begin{gather*}
+ \lambda = \frac{1}{r\cdot g},\\
+ \tau = r\cdot c.
+\end{gather*}
+
+With the model boundary conditions,
+\begin{subequations}
+ \begin{align}
+ v(x, 0) &= E, \\
+ \left.\frac{\partial v}{\partial x}\right\vert_{x=0} & = -Ir, \\
+ \left.\frac{\partial v}{\partial x}\right\vert_{x=L} & = 0,
+ \end{align}
+\end{subequations}
+where $I$ is the injected current and $L$ is the cable length.
+
+The solution $v(x, t)$ can be expressed in terms of the solution $g(x, t; L)$
+to a normalized version of the cable equation,
+\begin{subequations}
+ \begin{align}
+ \label{eq:normcable}
+ \frac{\partial^2 g}{\partial x^2} & =
+ \frac{\partial g}{\partial t} + g,
+ \\
+ \label{eq:normcableinitial}
+ g(x, 0) &= 0,
+ \\
+ \label{eq:normcableleft}
+ \left.\frac{\partial g}{\partial x}\right\vert_{x=0} & = 1,
+ \\
+ \label{eq:normcableright}
+ \left.\frac{\partial g}{\partial x}\right\vert_{x=L} & = 0
+ \end{align}
+\end{subequations}
+by
+\begin{equation}
+ v(x, t)= E - Ir\lambda \cdot g(\frac{x}{\lambda}, \frac{t}{\tau}; \frac{L}{\lambda}).
+\end{equation}
+
+\section{Solution to the normalized cable equation}
+
+Let $G(x, s)$ be the Laplace transform of $g(x, t)$ with respect to $t$.
+If the inverse transform of $G$ exists, it must agree with $g$ almost everywhere
+for $t>0$, or for all $t>0$ given the continuity of $g$.
+
+From
+\eqref{eq:normcable} and \eqref{eq:normcableinitial},
+\begin{equation}
+ \label{eq:lap}
+ \frac{\partial^2 G}{\partial x^2} = sG + G.
+\end{equation}
+The boundary conditions \eqref{eq:normcableleft} and \eqref{eq:normcableright} give
+\begin{align}
+ \label{eq:lapleft}
+ \frac{\partial G}{\partial x}(0, s) & = \frac{1}{s}
+ \\
+ \label{eq:lapright}
+ \frac{\partial G}{\partial x}(L, s) & = 0.
+\end{align}
+
+Solutions to \eqref{eq:lap} must be of the form
+\begin{equation}
+ G(x, s) = A(s) e^{mx} +B(s) e^{-mx}
+\end{equation}
+where $m=\sqrt{1+s}$. From \eqref{eq:lapleft} and \eqref{eq:lapright},
+\begin{align*}
+ mA(s) - mB(s) & = \frac{1}{s}
+ \\
+ mA(s)e^{mL} -mB(s)e^{-mL} & = 0
+\end{align*}
+and thus
+\begin{align*}
+ A(s) & = \frac{1}{sm (1-e^{2mL})}
+ \\
+ B(s) & = \frac{e^{2mL}}{sm (1-e^{2mL})}.
+\end{align*}
+
+Consequently,
+\begin{equation}
+ \begin{aligned}
+ G(x, s) &= \frac{1}{ms}\cdot\frac{e^{mx}+e^{2mL-mx}}{1-e^{2mL}}\\
+ &= - \frac{1}{ms}\cdot\frac{\cosh m(L-x)}{\sinh mL}.
+ \end{aligned}
+\end{equation}
+
+\subsection*{Inversion}
+
+Sufficient conditions for the inverse transform of $G(x,s)$ to exist
+and be representable in series form are as follows
+\parencite[][Theorem 4]{churchill1937}:
+\begin{enumerate}
+\item
+ $G(x, s)$ is analytic in some right half-plane $H$,
+\item
+ the singularities of $G(x, s)$ are all poles, and
+\item for some $k>1$, $|s^k G(x, s)|$ is bounded in $H$ and
+ on the circles $|s|=\rho_i$, where $\rho_i$ is an unbounded monotonically
+ increasing sequence.
+\end{enumerate}
+If in addition all the poles $\sigma_j$ of $G$ are simple, then
+the series expansion of the inverse transform is given by
+\begin{equation}
+ g(x,t) = \sum_j e^{\sigma_j t}\, \Res(G;\sigma_j),\quad t>0.
+\end{equation}
+
+$G(x,s)$ has poles when $s=0$, $m=0$, or $\sinh mL=0$.
+In the following, let
+\begin{gather*}
+ a_k = k\pi /L\\
+ m_k = a_k i\\
+ s_k = m_k^2 -1 = -k^2\pi^2/L^2-1.
+\end{gather*}
+for $k\in\mathbb{Z}$.
+The poles are then $0$, $-1$, and $s_k$ for $k\geq 1$.
+
+There are no branch points
+arising from $m=\sqrt{1+s}$, as letting $m=-\sqrt{1+s}$ leaves $G$ unchanged.
+
+For $|s+1|>\epsilon$,
+\begin{equation}
+ \begin{aligned}
+ |s^{3/2}G(x,s)|^2
+ & \leq (1+\epsilon)^{-1} \left| \frac{\cosh m(L-x)}{\sinh mL} \right|^2
+ \\
+ & \leq (1+\epsilon)^{-1} (1+|\coth mL|)^2
+ \label{eq:gbounds}
+ \end{aligned}
+\end{equation}
+which is bounded in the half-plane $\Re(s) \geq -1+\epsilon$.
+As $\coth u+iv$ is periodic in $v$, \eqref{eq:gbounds} is also bounded outside
+the regions $\{s: m\in D_k\}$, where $D_k$ are non-overlapping $\delta$-neighborhoods of $m_k$
+for some fixed small $\delta>0$. These are contained within disks of radius
+$2k\pi\delta/L+\delta^2$ about $s_k$ for $k\geq 1$, which then are separated by
+circles $|s|=\rho_k$ about the origin for some $\rho_k$ in $(|s_k|, |s_k+1|)$.
+
+\textbf{Pole at $s=0$}.\\
+Recalling $m=\sqrt{1+s}$, $m$ and $\sinh mL$ are non-zero in
+a neighbourhood of $s=0$, and so the pole is simple and
+\begin{equation}
+ \begin{aligned}
+ \Res(G; 0) & = - \frac{1}{m}\cdot\left.\frac{\cosh m(L-x)}{\sinh mL}\right|_{s=0}\\
+ & = - \frac{\cosh (L-x)}{\sinh L}.
+ \end{aligned}
+\end{equation}
+
+\textbf{Pole at $s=-1$}.\\
+Let $G(x,s)=f(x,s)/h(s)$, where
+\begin{equation}
+ \begin{aligned}
+ f(x, s) &= -\frac{1}{s}\cosh m(L-x)\\
+ h(s) &= m\sinh mL.
+ \end{aligned}
+ \label{eq:hf}
+\end{equation}
+Noting that $dm/ds = \frac{1}{2}m^{-1}$,
+\begin{equation}
+ \begin{aligned}
+ h'(s) &= \frac{1}{2}m^{-1}\sinh mL + \frac{1}{2}L\cosh mL \\
+ &= \frac{1}{2}L + \frac{1}{2}L + O(m^2) \quad(m\to 0) \\
+ &= L + O(s+1) \quad(s\to -1).
+ \label{eq:hprime}
+ \end{aligned}
+\end{equation}
+The pole is therefore simple, and
+\begin{equation}
+ \Res(G; -1) = f(x, -1)/h'(-1) = -\frac{1}{L}.
+\end{equation}
+
+\textbf{Pole at $s=s_k$, $k \geq 1$}.\\
+Taking $h$ and $f$ as above \eqref{eq:hf},
+$m_k$ is non-zero for $k\geq 1$ and
+\[
+ h'(s_k) = \frac{1}{2}L\cosh m_kL.
+\]
+Consequently the pole is simple and
+\begin{equation}
+ \begin{aligned}
+ \Res(G; s_k)
+ & = f(x, s_k)/h'(s_k)\\
+ & = -\frac{2}{s_k L}\frac{\cosh m_k(L-x)}{\cosh m_kL} \\
+ & = -\frac{2}{s_k L}\frac{\cosh m_kL\cosh m_kx-\sinh m_kL\sinh m_kL}{\cosh m_kL} \\
+ & = -\frac{2}{s_k L}\cosh m_k x,
+ \end{aligned}
+\end{equation}
+as $\sinh m_k=0$.
+
+In terms of $a_k$,
+\begin{equation}
+ \Res(G; s_k) = \frac{2}{L}\cdot\frac{1}{1+a_k^2}\cdot\cos a_kx.
+\end{equation}
+
+The series exapnsion for $g(x, t)$ therefore is
+\begin{equation}
+ g(x, t) = -\frac{\cosh(L-x)}{\sinh L} + \frac{1}{L}e^{-t}\left\{
+ 1+2\sum_{k=1}^\infty \frac{e^{-ta_k^2}}{1+a_k^2}\cos a_k x\right\}.
+ \label{eq:theg}
+\end{equation}
+
+\section{Computation of series}
+
+As $t$ approaches zero, the series \eqref{eq:theg} converges increasingly slowly.
+For computational purposes, an estimate on the series residual allows the determination
+of stopping criteria for a given tolerance.
+
+Let $g_n$ be the partial sum
+\begin{equation}
+ g_n(x, t) = -\frac{\cosh(L-x)}{\sinh L} + \frac{1}{L}e^{-t}\left\{
+ 1+2\sum_{k=1}^n \frac{e^{-ta_k^2}}{1+a_k^2}\cos a_k x\right\}.
+\end{equation}
+so that $g(x, t) =\lim_{n\to\infty} g_n(x,t)$. Let $\bar{g}_n = |g-g_n|$ be the
+residual. The $a_k$ form an increasing sequence, so
+\begin{equation}
+ \begin{aligned}
+ \bar{g}_n(x,t)
+ & \leq \frac{2}{L}e^{-t}\sum_{n+1}^\infty\frac{e^{-ta_k^2}}{1+a_k^2}\\
+ & \leq \frac{2}{L}e^{-t}\int_{a_n}^\infty \frac{e^{-tu^2}}{1+u^2}\,du\\
+ & < \frac{2}{L}e^{-t}\int_{a_n}^\infty \frac{e^{-tu^2}}{u^2}\,du.
+ \end{aligned}
+ \label{eq:gbar}
+\end{equation}
+Substituting $u=\sqrt{v/t}$ gives the identity
+\begin{equation}
+ \int_{a}^\infty \frac{e^{-tu^2}}{u^2}\,du
+ = \frac{1}{2}\sqrt{t}\int_{a^2t}^\infty e^{-tv} v^{\frac{-3}{2}}\,dv
+ = \frac{1}{2}\sqrt{t}\,\Gamma(-\frac{1}{2},a^2 t),
+\end{equation}
+where $\Gamma(\alpha, z)$ is the upper incomplete gamma function.
+
+For real $\alpha<1$ and $z>0$, \textcite[][Theorem 2.3]{borwein2009} give the upper bound
+\[
+ \Gamma(\alpha, z) \leq z^{\alpha-1} e^{-z}.
+\]
+Substituting into \eqref{eq:gbar} gives
+\begin{equation}
+ \begin{aligned}
+ \bar{g}_n(x,t)
+ & < \frac{1}{L}e^{-t}\sqrt{t}\,\Gamma(-\frac{1}{2},a_n^2 t) \\
+ & \leq \frac{1}{L}e^{-t}\sqrt{t}\,(a_n^2 t)^{-\frac{3}{2}}e^{-a_n^2t} \\
+ & = \frac{e^{-t(1+a_n^2)}}{L t a_n^3}.
+ \end{aligned}
+\end{equation}
+
+\printbibliography
+\end{document}
diff --git a/docs/passive_cable/makefile b/docs/passive_cable/makefile
new file mode 100644
index 0000000000000000000000000000000000000000..8756af0b9b8f33a1a145a07259d07a335a86de86
--- /dev/null
+++ b/docs/passive_cable/makefile
@@ -0,0 +1,14 @@
+SOURCES := cable_computation.tex
+
+LATEXMK := latexmk -e '$$clean_ext=q/bbl run.xml/' -xelatex -use-make -halt-on-error
+
+all: cable_computation.pdf
+
+cable_computation.pdf: cable_computation.tex cable.bib
+ $(LATEXMK) $<
+
+clean:
+ for s in $(SOURCES); do $(LATEXMK) -c "$$s"; done
+
+realclean:
+ for s in $(SOURCES); do $(LATEXMK) -C "$$s"; done
diff --git a/scripts/PassiveCable.jl b/scripts/PassiveCable.jl
new file mode 100644
index 0000000000000000000000000000000000000000..76450490ce5c4f265403edd35f2bedad4ad7c3bb
--- /dev/null
+++ b/scripts/PassiveCable.jl
@@ -0,0 +1,141 @@
+module PassiveCable
+
+export cable_normalize, cable, rallpack1
+
+# Compute solution g(x, t) to
+#
+# ∂²g/∂x² - g - ∂g/∂t = 0
+#
+# on [0, L] × [0,∞), subject to:
+#
+# g(x, 0) = 0
+# ∂g/∂x (0, t) = 1
+# ∂g/∂x (L, t) = 0
+#
+# Parameters:
+# x, t, L: as described above
+# tol: absolute error tolerance in result
+#
+# Return:
+# g(x, t)
+
+function cable_normalized(x, t, L; tol=1e-8)
+ if t<=0
+ return 0.0
+ else
+ ginf = -cosh(L-x)/sinh(L)
+ sum = exp(-t/L)
+
+ for k = countfrom(1)
+ a = k*pi/L
+ e = exp(-t*(1+a^2))
+
+ sum += 2/L*e*cos(a*x)/(1+a^2)
+ resid_ub = e/(L*a^3*t)
+
+ if resid_ub<tol
+ break
+ end
+ end
+ return ginf+sum
+ end
+end
+
+
+# Compute solution f(x, t) to
+#
+# λ²∂²f/∂x² - f - τ∂f/∂t = 0
+#
+# on [0, L] x [0, ∞), subject to:
+#
+# f(x, 0) = V
+# ∂f/∂x (0, t) = I·r
+# ∂f/∂x (L, t) = 0
+#
+# where:
+#
+# λ² = 1/(r·g) length constant
+# τ = r·c time constant
+#
+# In the physical model, the parameters correspond to the following:
+#
+# L: length of cable
+# r: linear axial resistivity
+# g: linear membrane conductivity
+# c: linear membrane capacitance
+# V: membrane reversal potential
+# I: injected axial current on the left end (x = 0) of the cable.
+#
+# Note that r, g and c are specific 1-d quantities that differ from
+# the cable resistivity r_L, specific membrane conductivity ḡ and
+# specific membrane capacitance c_m as used elsewhere. If the
+# cross-sectional area is A and cable circumference is f, then
+# these quantities are related by:
+#
+# r = r_L/A
+# g = ḡ·f
+# c = c_m·f
+#
+# Parameters:
+# x: displacement along cable
+# t: time
+# L, lambda, tau, r, V, I: as described above
+# tol: absolute error tolerance in result
+#
+# Return:
+# computed potential at (x,t) on cable.
+
+function cable(x, t, L, lambda, tau, r, V, I; tol=1e-8)
+ scale = I*r*lambda;
+ if scale == 0
+ return V
+ else
+ tol_n = abs(tol/scale)
+ return scale*cable_normalized(x/lambda, t/tau, L/lambda, tol=tol_n) + V
+ end
+end
+
+
+# Rallpack 1 test
+#
+# One sided cable model with the following parameters:
+#
+# RA = 1 Ω·m bulk axial resistivity
+# RM = 4 Ω·m² areal membrane resistivity
+# CM = 0.01 F/m² areal membrane capacitance
+# d = 1 µm cable diameter
+# EM = -65 mV reversal potential
+# I = 0.1 nA injected current
+# L = 1 mm cable length.
+#
+# (This notation aligns with that used in the Rallpacks paper.)
+#
+# Note that the injected current as described in the Rallpack model
+# is trans-membrane, not axial. Consequently we need to swap the
+# sign on I when passing to the cable function.
+#
+# Parameters:
+# x: displacement along cable [m]
+# t: time [s]
+# tol: absolute tolerance for reported potential.
+#
+# Return:
+# computed potential at (x,t) on cable.
+
+function rallpack1(x, t; tol=1e-8)
+ RA = 1
+ RM = 4
+ CM = 1e-2
+ d = 1e-6
+ EM = -65e-3
+ I = 0.1e-9
+ L = 1e-3
+
+ r = 4*RA/(pi*d*d)
+ lambda = sqrt(d/4 * RM/RA)
+ tau = CM*RM
+
+ return cable(x, t, L, lambda, tau, r, EM, -I, tol=tol)
+end
+
+end # module PassiveCable
diff --git a/scripts/README.md b/scripts/README.md
index e645a762eaf88c5f7196a995502589da44c2f89c..373f9f1330840c7cddf9fcf264d92f943c55c33b 100644
--- a/scripts/README.md
+++ b/scripts/README.md
@@ -1,5 +1,4 @@
-tsplot
-------
+#tsplot
The `tsplot` script is a wrapper around matplotlib for displaying a collection of
time series plots.
@@ -78,4 +77,58 @@ of the timeseries data.
Use the `-o` or `--output` option to save the plot as an image, instead of
displaying it interactively.
+#profstats
+`profstats` collects the profiling data output from multiple MPI ranks and performs
+a simple statistical summary.
+
+Input files are in the JSON format emitted by the profiling code.
+
+By default, `profstats` reports the quartiles of the times reported for each
+profiling region and subregion. With the `-r` option, the collated raw times
+are reported instead.
+
+Output is in CSV format.
+
+#PassiveCable.jl
+
+Compute analytic solutions to the simple passive cylindrical dendrite cable
+model with step current injection at one end from t=0.
+
+This is used to generate validation data for the first Rallpack test.
+
+Module exports the following functions:
+
+ * `cable_normalized(x, t, L; tol)`
+
+ Compute potential _V_ at position _x_ in [0, _L_] at time _t_ ≥ 0 according
+ to the normalized cable equation with unit length constant and time constant.
+
+ Neumann boundary conditions: _V'_(0) = 1; _V'_(L) = 0.
+ Initial conditions: _V_( _x_, 0) = 0.
+
+ Absolute tolerance `tol` defaults to 1e-8.
+
+ * `cable(x, t, L, lambda, tau, r, V, I; tol)`
+
+ Compute the potential given:
+ * length constant `lambda`
+ * time constant `tau`,
+ * axial linear resistivity `r`
+ * injected current of `I` at the origin
+ * reversal potential `V`
+
+ Implied units must be compatible, e.g. SI units.
+
+ Absolute tolerance `tol` defaults to 1e-8.
+
+ * `rallpack1(x, t; tol)`
+
+ Compute the value of the potential in the Rallpack 1 test model at position
+ `x` and time `t`.
+
+ Parameters for the underlying cable equation calculation are taken from the
+ Rallpack model description in SI units; as the cable length is 1 mm in this
+ model, `x` can take values in [0, 0.001].
+
+ Absolute tolerance `tol` defaults to 1e-8.
diff --git a/scripts/profstats b/scripts/profstats
index 33f37356836104e57758f0407e1da5fd0625a111..88f68c72e6253fa4c0d240d051d7ce1be1960604 100755
--- a/scripts/profstats
+++ b/scripts/profstats
@@ -1,4 +1,4 @@
-#!env python
+#!/usr/bin/env python2
#coding: utf-8
import json
diff --git a/scripts/tsplot b/scripts/tsplot
index 9ad6e2a5c714a5f634e37c31c88411f63923f130..8e89cbda96451680bfc9de8087ba558c0282defc 100755
--- a/scripts/tsplot
+++ b/scripts/tsplot
@@ -1,4 +1,4 @@
-#!env python
+#!/usr/bin/env python2
#coding: utf-8
import argparse
diff --git a/src/cell.hpp b/src/cell.hpp
index 11f938f52f3e4e69a35f9a3eb951ca2fce95cbd6..c691e68422b856183c941d0fd1c4172d644f3289 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -90,8 +90,11 @@ public:
spike_detectors_(other.spike_detectors_),
probes_(other.probes_)
{
- util::assign_by(segments_, other.segments_,
- [](const segment_ptr& p) { return p->clone(); });
+ // unique_ptr's cannot be copy constructed, do a manual assignment
+ segments_.reserve(other.segments_.size());
+ for (const auto& s: other.segments_) {
+ segments_.push_back(s->clone());
+ }
}
/// add a soma to the cell
@@ -242,4 +245,3 @@ cable_segment* cell::add_cable(cell::index_type parent, Args ...args)
} // namespace mc
} // namespace nest
-
diff --git a/src/util/rangeutil.hpp b/src/util/rangeutil.hpp
index 46101e3eef8160124d73d60e730c31d34c85f7c1..839fb4527be34cf04777c65d7e934712aa7c57e3 100644
--- a/src/util/rangeutil.hpp
+++ b/src/util/rangeutil.hpp
@@ -7,10 +7,12 @@
#include <algorithm>
#include <iterator>
+#include <numeric>
#include <util/meta.hpp>
#include <util/range.hpp>
#include <util/transform.hpp>
+#include <util/meta.hpp>
namespace nest {
namespace mc {
@@ -66,12 +68,12 @@ AssignableContainer& assign(AssignableContainer& c, const Seq& seq) {
return c;
}
+
// Assign sequence to a container with transform `proj`
template <typename AssignableContainer, typename Seq, typename Proj>
AssignableContainer& assign_by(AssignableContainer& c, const Seq& seq, const Proj& proj) {
- auto canon = canonical_view(transform_view(seq, proj));
- c.assign(std::begin(canon), std::end(canon));
+ assign(c, transform_view(seq, proj));
return c;
}
diff --git a/tests/unit/test_cell.cpp b/tests/unit/test_cell.cpp
index 3102d97435d3095b4ef91b885b040ac084f37a5f..eb28e42be495a425a7ff68256453c15e027929ab 100644
--- a/tests/unit/test_cell.cpp
+++ b/tests/unit/test_cell.cpp
@@ -165,3 +165,58 @@ TEST(cell_type, multiple_cables)
}
}
+TEST(cell_type, clone)
+{
+ using namespace nest::mc;
+
+ // make simple cell with multiple segments
+
+ cell c;
+ c.add_soma(2.1);
+ c.add_cable(0, segmentKind::dendrite, 0.3, 0.2, 10);
+ c.segment(1)->set_compartments(3);
+ c.add_cable(1, segmentKind::dendrite, 0.2, 0.15, 20);
+ c.segment(2)->set_compartments(5);
+
+ parameter_list exp_default("expsyn");
+ c.add_synapse({1, 0.3}, exp_default);
+
+ c.add_detector({0, 0.5}, 10.0);
+
+ // make clone
+
+ cell d(clone_cell, c);
+
+ // check equality
+
+ ASSERT_EQ(c.num_segments(), d.num_segments());
+ EXPECT_EQ(c.soma()->radius(), d.soma()->radius());
+ EXPECT_EQ(c.segment(1)->as_cable()->length(), d.segment(1)->as_cable()->length());
+ {
+ const auto& csyns = c.synapses();
+ const auto& dsyns = d.synapses();
+
+ ASSERT_EQ(csyns.size(), dsyns.size());
+ for (unsigned i=0; i<csyns.size(); ++i) {
+ ASSERT_EQ(csyns[i].location, dsyns[i].location);
+ }
+ }
+
+ ASSERT_EQ(1u, c.detectors().size());
+ ASSERT_EQ(1u, d.detectors().size());
+ EXPECT_EQ(c.detectors()[0].threshold, d.detectors()[0].threshold);
+
+ // check clone is independent
+
+ c.add_cable(2, segmentKind::dendrite, 0.15, 0.1, 20);
+ EXPECT_NE(c.num_segments(), d.num_segments());
+
+ d.detectors()[0].threshold = 13.0;
+ ASSERT_EQ(1u, c.detectors().size());
+ ASSERT_EQ(1u, d.detectors().size());
+ EXPECT_NE(c.detectors()[0].threshold, d.detectors()[0].threshold);
+
+ c.segment(1)->set_compartments(7);
+ EXPECT_NE(c.segment(1)->num_compartments(), d.segment(1)->num_compartments());
+ EXPECT_EQ(c.segment(2)->num_compartments(), d.segment(2)->num_compartments());
+}
diff --git a/tests/validation/hh_soma.jl b/tests/validation/hh_soma.jl
new file mode 100644
index 0000000000000000000000000000000000000000..82e7df8b5f6c2c594b8b35f584cbf41a42e58e6b
--- /dev/null
+++ b/tests/validation/hh_soma.jl
@@ -0,0 +1,132 @@
+using Sundials
+using SIUnits.ShortUnits
+
+c_m = 0.01nF*m^-2
+
+radius = 18.8μm / 2
+surface_area = 4 * pi * radius * radius
+
+gnabar = .12S*cm^-2
+gkbar = .036S*cm^-2
+gl = .0003S*cm^-2
+el = -54.3mV
+#celsius= 6.3
+#q10 = 3^((celsius - 6.3)/10)
+q10 = 1
+
+# define the resting potential for the membrane
+vrest = -65.0mV
+
+# define the resting potentials for ion species
+ena = 115.0mV+vrest
+ek = -12.0mV+vrest
+eca = 12.5mV*log(2.0/5e-5)
+
+vtrap(x,y) = x/(exp(x/y) - 1.0)
+
+function print()
+ println("q10 ", q10)
+ println("vrest ", vrest)
+ println("ena ", ena)
+ println("ek ", ek)
+ println("eca ", eca)
+end
+
+#"m" sodium activation system
+function m_lims(v)
+ alpha = .1mV^-1 * vtrap(-(v+40mV),10mV)
+ beta = 4 * exp(-(v+65mV)/18mV)
+ sum = alpha + beta
+ mtau = 1ms / (q10*sum)
+ minf = alpha/sum
+ return mtau, minf
+end
+
+#"h" sodium inactivation system
+function h_lims(v)
+ alpha = 0.07*exp(-(v+65mV)/20mV)
+ beta = 1 / (exp(-(v+35mV)/10mV) + 1)
+ sum = alpha + beta
+ htau = 1ms / (q10*sum)
+ hinf = alpha/sum
+ return htau, hinf
+end
+
+#"n" potassium activation system
+function n_lims(v)
+ alpha = .01mV^-1 * vtrap(-(v+55mV),10mV)
+ beta = .125*exp(-(v+65mV)/80mV)
+ sum = alpha + beta
+ ntau = 1ms / (q10*sum)
+ ninf = alpha/sum
+ return ntau, ninf
+end
+
+# v = y[1]
+# m = y[2]
+# h = y[3]
+# n = y[4]
+
+# choose initial conditions for the system such that the gating variables
+# are at steady state for the user-specified voltage v
+function initial_conditions(v)
+ mtau, minf = m_lims(v)
+ htau, hinf = h_lims(v)
+ ntau, ninf = n_lims(v)
+
+ return [Float64(v), minf, hinf, ninf]
+end
+
+# calculate the lhs of the ODE system
+function f(t, y, ydot)
+ # copy variables into helper variable
+ v = y[1]mV
+ m, h, n = y[2], y[3], y[4]
+
+ # calculate current due to ion channels
+ gna = gnabar*m*m*m*h
+ gk = gkbar*n*n*n*n
+ ina = gna*(v - ena)
+ ik = gk*(v - ek)
+ il = gl*(v - el)
+ imembrane = ik + ina + il
+
+ # calculate current due to stimulus
+ #c.add_stimulus({0,0.5}, {10., 100., 0.1});
+ ielectrode = 0.0nA / surface_area
+ if t>=Float64(10ms) && t<Float64(100ms)
+ ielectrode = 0.1nA / surface_area
+ end
+
+ # calculate the total membrane current
+ i = -imembrane + ielectrode
+
+ # calculate the voltage dependent rates for the gating variables
+ mtau, minf = m_lims(v)
+ ntau, ninf = n_lims(v)
+ htau, hinf = h_lims(v)
+
+ # set the derivatives
+ # note tha these are in SI units, which are indicated in comments
+ ydot[1] = Float64(i/c_m) # V*s^-1
+ ydot[2] = Float64((minf-m)/mtau) # s^-1
+ ydot[3] = Float64((hinf-h)/htau) # s^-1
+ ydot[4] = Float64((ninf-n)/ntau) # s^-1
+
+ return Sundials.CV_SUCCESS
+end
+
+###########################################################
+# now we actually run the model
+###########################################################
+
+# from 0 to 100 ms in 1000 steps
+t = collect(linspace(0.0, 0.1, 1001));
+
+# these tolerances are as tight as they will go without breaking convergence of the iterative schemes
+y0 = initial_conditions(vrest)
+res = Sundials.cvode(f, y0, t, abstol=1e-6, reltol=5e-10);
+
+#using Plots
+#gr()
+#plot(t, res[:,1])
diff --git a/tests/validation/validate.cpp b/tests/validation/validate.cpp
index e3d457b939cac1c584588c8c1f94672a2d655587..de70f9bf477e1ef2cb7c2fb67b29968c29a2f18a 100644
--- a/tests/validation/validate.cpp
+++ b/tests/validation/validate.cpp
@@ -29,7 +29,6 @@ int main(int argc, char **argv) {
communication::global_policy_guard global_guard(argc, argv);
::testing::InitGoogleTest(&argc, argv);
- int rv = 0;
try {
auto arg = argv+1;
while (*arg) {
@@ -54,8 +53,9 @@ int main(int argc, char **argv) {
}
}
- rv = RUN_ALL_TESTS();
+ return RUN_ALL_TESTS();
}
+
catch (to::parse_opt_error& e) {
to::usage(argv[0], usage_str, e.what());
return 1;