diff --git a/CMakeLists.txt b/CMakeLists.txt
index 62c167a698d58d3ed8470ef14d430eebc156e0f1..86bd6f78cf4ba6d419ac0242d088aa54cca6a8cc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -88,16 +88,34 @@ else()
set(use_external_modcc ON BOOL)
endif()
-# whether to attempt to use nrniv to build validation data
+# Validation data generation
+
+# destination directory for generated data
+set(VALIDATION_DATA_DIR "${CMAKE_SOURCE_DIR}/validation/data" CACHE PATH "location of generated validation data")
+
+# Whether to build validation data at all
+set(BUILD_VALIDATION_DATA ON CACHE BOOL "generate validation data")
+
+# Whether to attempt to use julia to build validation data
+find_program(JULIA_BIN julia)
+if(JULIA_BIN STREQUAL "JULIA_BIN-NOTFOUND")
+ message(STATUS "julia not found; will not automatically build validation data sets from julia scripts")
+ set(BUILD_JULIA_VALIDATION_DATA FALSE)
+else()
+ set(BUILD_JULIA_VALIDATION_DATA TRUE)
+endif()
+
+# Whether to attempt to use nrniv to build validation data
# (if we find nrniv, do)
find_program(NRNIV_BIN nrniv)
if(NRNIV_BIN STREQUAL "NRNIV_BIN-NOTFOUND")
- message(STATUS "nrniv not found; will not automatically build validation data sets")
- set(BUILD_VALIDATION_DATA FALSE)
+ message(STATUS "nrniv not found; will not automatically build NEURON validation data sets")
+ set(BUILD_NRN_VALIDATION_DATA FALSE)
else()
- set(BUILD_VALIDATION_DATA TRUE)
+ set(BUILD_NRN_VALIDATION_DATA TRUE)
endif()
+
include_directories(${CMAKE_SOURCE_DIR}/tclap/include)
include_directories(${CMAKE_SOURCE_DIR}/vector)
include_directories(${CMAKE_SOURCE_DIR}/include)
@@ -109,12 +127,17 @@ if( "${WITH_TBB}" STREQUAL "ON" )
include_directories(${TBB_INCLUDE_DIRS})
endif()
+# only include validation data if flag is set
+if(BUILD_VALIDATION_DATA)
+ add_subdirectory(validation)
+endif()
+
# only compile modcc if it is not provided externally
if(use_external_modcc)
add_subdirectory(modcc)
endif()
+
add_subdirectory(mechanisms)
-add_subdirectory(nrn)
add_subdirectory(src)
add_subdirectory(tests)
add_subdirectory(miniapp)
diff --git a/docs/passive_cable/cable_computation.tex b/docs/passive_cable/cable_computation.tex
index 39a372196c006c7c67837f039df88c903cb439b4..888aa8a11890879818e2964ba4644744e434ab79 100644
--- a/docs/passive_cable/cable_computation.tex
+++ b/docs/passive_cable/cable_computation.tex
@@ -54,16 +54,16 @@ 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
+ \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}
@@ -87,9 +87,9 @@ and the linear membrane capacitance $c$. These determine $\lambda$ and $\tau$ by
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,
+ 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.
@@ -98,18 +98,18 @@ 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: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
+ \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
@@ -158,7 +158,7 @@ and thus
Consequently,
\begin{equation}
\begin{aligned}
- G(x, s) &= \frac{1}{ms}\cdot\frac{e^{mx}+e^{2mL-mx}}{1-e^{2mL}}\\
+ 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}
@@ -199,10 +199,10 @@ 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
+ |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}
@@ -218,8 +218,8 @@ 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}.
+ \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}
@@ -235,9 +235,9 @@ Let $G(x,s)=f(x,s)/h(s)$, where
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).
+ 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}
@@ -255,11 +255,11 @@ $m_k$ is non-zero for $k\geq 1$ and
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,
+ \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$.
@@ -272,7 +272,7 @@ In terms of $a_k$,
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\}.
+ 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}
@@ -285,16 +285,16 @@ 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\}.
+ 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.
+ \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}
@@ -313,9 +313,9 @@ For real $\alpha<1$ and $z>0$, \textcite[][Theorem 2.3]{borwein2009} give the up
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} \\
+ \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}
diff --git a/nrn/CMakeLists.txt b/nrn/CMakeLists.txt
deleted file mode 100644
index 28325c091cea5aec8b0666ac46ba6503e82cdc9a..0000000000000000000000000000000000000000
--- a/nrn/CMakeLists.txt
+++ /dev/null
@@ -1,26 +0,0 @@
-# The validation scripts to run (without .py extension)
-
-set(validations
- ball_and_stick
- ball_and_3stick
- ball_and_taper
- simple_exp_synapse
- simple_exp2_synapse
- soma)
-
-# Only try and make validation sets if we can find nrniv
-if(BUILD_VALIDATION_DATA)
- set(common "${CMAKE_CURRENT_SOURCE_DIR}/nrn_validation.py")
- foreach(v ${validations})
- set(out "${CMAKE_SOURCE_DIR}/data/validation/neuron_${v}.json")
- set(src "${CMAKE_CURRENT_SOURCE_DIR}/${v}.py")
- add_custom_command(
- OUTPUT "${out}"
- DEPENDS "${src}" "${common}"
- WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
- COMMAND ${NRNIV_BIN} -nobanner -python ${src} > ${out})
- list(APPEND all_neuron_validation "${out}")
- endforeach()
- add_custom_target(validation_data DEPENDS ${all_neuron_validation})
-endif()
-
diff --git a/scripts/tsplot b/scripts/tsplot
index 8e89cbda96451680bfc9de8087ba558c0282defc..4280a9dc902cfa738e5eaf8bb7043008d96baae0 100755
--- a/scripts/tsplot
+++ b/scripts/tsplot
@@ -44,6 +44,7 @@ def parse_clargs():
help='plot series with same KEYs on the same axes')
P.add_argument('-s', '--select', metavar='EXPR,...', dest='select',
type=lambda s: s.split(','),
+ action='append',
help='select only series matching EXPR')
P.add_argument('-o', '--output', metavar='FILE', dest='outfile',
help='save plot to file FILE')
@@ -212,7 +213,7 @@ def read_json_timeseries(j, axis='time', select=[]):
meta = dict(j.items() + {'label': key, 'data': None, 'units': units(key)}.items())
del meta['data']
- if all([run_select(sel, meta) for sel in select]):
+ if not select or any([all([run_select(s, meta) for s in term]) for term in select]):
ts_list.append(TimeSeries(times, jdata[key], **meta))
return ts_list
@@ -422,7 +423,7 @@ args = parse_clargs()
tss = []
axis = args.axis if args.axis else 'time'
for filename in args.inputs:
- select = args.select if args.select else []
+ select = args.select
with open(filename) as f:
j = json.load(f)
tss.extend(read_json_timeseries(j, axis, select))
diff --git a/src/math.hpp b/src/math.hpp
index 4d8ee51c103cc74107938135b8094276fde10c47..bf964ba38c07487914ffcf5768b7f39aee63a5d1 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -1,9 +1,9 @@
#pragma once
#include <cmath>
+#include <limits>
#include <utility>
-
namespace nest {
namespace mc {
namespace math {
@@ -14,6 +14,12 @@ T constexpr pi()
return T(3.1415926535897932384626433832795);
}
+template <typename T = float>
+T constexpr infinity()
+{
+ return std::numeric_limits<T>::infinity();
+}
+
template <typename T>
T constexpr mean(T a, T b)
{
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index eb3ee8fe2f3c3ce1018bb113c966cb556e12916e..caa752b9489a483e91f718350bd5c456ded0ae66 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -155,7 +155,7 @@ namespace mc {
// 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
+ base::add_parameter({"r_L", 100.00, {0., 1e9}}); // equivalent to Ra in Neuron : Ohm.cm
}
};
diff --git a/tests/test_common_cells.hpp b/tests/test_common_cells.hpp
index 9e302365dbea18cff2b4f82a7e74905c9143c5ed..c6863b71712960ca28ae0f9e8d08c7b2bff138ef 100644
--- a/tests/test_common_cells.hpp
+++ b/tests/test_common_cells.hpp
@@ -1,4 +1,7 @@
+#include <cmath>
+
#include <cell.hpp>
+#include <math.hpp>
#include <parameter_list.hpp>
namespace nest {
@@ -9,8 +12,9 @@ namespace mc {
*
* Soma:
* diameter: 18.8 µm
- * mechanisms: membrane, HH
- * memrane resistance: 123 Ω·cm
+ * mechanisms: HH (default params)
+ * bulk resistivitiy: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
*
* Stimuli:
* soma centre, t=[10 ms, 110 ms), 0.1 nA
@@ -20,7 +24,7 @@ inline cell make_cell_soma_only(bool with_stim = true) {
cell c;
auto soma = c.add_soma(18.8/2.0);
- soma->mechanism("membrane").set("r_L", 123);
+ soma->mechanism("membrane").set("r_L", 100);
soma->add_mechanism(hh_parameters());
if (with_stim) {
@@ -33,15 +37,19 @@ inline cell make_cell_soma_only(bool with_stim = true) {
/*
* Create cell with a soma and unbranched dendrite:
*
+ * Common properties:
+ * bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
+ *
* Soma:
- * mechanisms: HH
+ * mechanisms: HH (default params)
* diameter: 12.6157 µm
*
* Dendrite:
- * mechanisms: none
+ * mechanisms: passive (default params)
* diameter: 1 µm
* length: 200 µm
- * membrane resistance: 100 Ω·cm
+ * bulk resistivity: 100 Ω·cm
* compartments: 4
*
* Stimulus:
@@ -54,10 +62,15 @@ inline cell make_cell_ball_and_stick(bool with_stim = true) {
auto soma = c.add_soma(12.6157/2.0);
soma->add_mechanism(hh_parameters());
- auto dendrite = c.add_cable(0, segmentKind::dendrite, 1.0/2, 1.0/2, 200.0);
- dendrite->add_mechanism(pas_parameters());
- dendrite->mechanism("membrane").set("r_L", 100);
- dendrite->set_compartments(4);
+ c.add_cable(0, segmentKind::dendrite, 1.0/2, 1.0/2, 200.0);
+
+ for (auto& seg: c.segments()) {
+ seg->mechanism("membrane").set("r_L", 100);
+ if (seg->is_dendrite()) {
+ seg->add_mechanism(pas_parameters());
+ seg->set_compartments(4);
+ }
+ }
if (with_stim) {
c.add_stimulus({1,1}, {5., 80., 0.3});
@@ -68,16 +81,20 @@ inline cell make_cell_ball_and_stick(bool with_stim = true) {
/*
* Create cell with a soma and unbranched tapered dendrite:
*
+ * Common properties:
+ * bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
+ *
* Soma:
- * mechanisms: HH
+ * mechanisms: HH (default params)
* diameter: 12.6157 µm
*
* Dendrite:
- * mechanisms: none
+ * mechanisms: passive (default params)
* diameter proximal: 1 µm
* diameter distal: 0.2 µm
* length: 200 µm
- * membrane resistance: 100 Ω·cm
+ * bulk resistivity: 100 Ω·cm
* compartments: 4
*
* Stimulus:
@@ -90,10 +107,15 @@ inline cell make_cell_ball_and_taper(bool with_stim = true) {
auto soma = c.add_soma(12.6157/2.0);
soma->add_mechanism(hh_parameters());
- auto dendrite = c.add_cable(0, segmentKind::dendrite, 1.0/2, 0.2/2, 200.0);
- dendrite->add_mechanism(pas_parameters());
- dendrite->mechanism("membrane").set("r_L", 100);
- dendrite->set_compartments(4);
+ c.add_cable(0, segmentKind::dendrite, 1.0/2, 0.2/2, 200.0);
+
+ for (auto& seg: c.segments()) {
+ seg->mechanism("membrane").set("r_L", 100);
+ if (seg->is_dendrite()) {
+ seg->add_mechanism(pas_parameters());
+ seg->set_compartments(4);
+ }
+ }
if (with_stim) {
c.add_stimulus({1,1}, {5., 80., 0.3});
@@ -106,15 +128,18 @@ inline cell make_cell_ball_and_taper(bool with_stim = true) {
*
* O----======
*
+ * Common properties:
+ * bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
+ *
* Soma:
- * mechanisms: HH
+ * mechanisms: HH (default params)
* diameter: 12.6157 µm
*
* Dendrites:
- * mechanisms: membrane
+ * mechanisms: passive (default params)
* diameter: 1 µm
* length: 100 µm
- * membrane resistance: 100 Ω·cm
* compartments: 4
*
* Stimulus:
@@ -128,15 +153,14 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
auto soma = c.add_soma(12.6157/2.0);
soma->add_mechanism(hh_parameters());
- // add dendrite of length 200 um and diameter 1 um with passive channel
c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100);
c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
for (auto& seg: c.segments()) {
+ seg->mechanism("membrane").set("r_L", 100);
if (seg->is_dendrite()) {
seg->add_mechanism(pas_parameters());
- seg->mechanism("membrane").set("r_L", 100);
seg->set_compartments(4);
}
}
@@ -148,6 +172,76 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
return c;
}
+/*
+ * Create 'soma-less' cell with single cable, with physical
+ * parameters from Rallpack 1 model.
+ *
+ * Common properties:
+ * mechanisms: passive
+ * membrane conductance: 0.000025 S/cm² ( = 1/(4Ω·m²) )
+ * membrane reversal potential: -65 mV (default)
+ * diameter: 1 µm
+ * length: 1000 µm
+ * bulk resistivity: 100 Ω·cm
+ * capacitance: 0.01 F/m² [default]
+ * compartments: 4
+ *
+ * Stimulus:
+ * end of dendrite, t=[0 ms, inf), 0.1 nA
+ *
+ * Note: zero-volume soma added with same mechanisms, as
+ * work-around for some existing fvm modelling issues.
+ *
+ * TODO: Set the correct values when parameters are generally
+ * settable!
+ *
+ * We can't currently change leak parameters
+ * from defaults, so we scale other electrical parameters
+ * proportionally.
+ */
+
+inline cell make_cell_simple_cable(bool with_stim = true) {
+ cell c;
+
+ c.add_soma(0);
+ c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 1000);
+
+ double r_L = 100;
+ double c_m = 0.01;
+ double gbar = 0.000025;
+ double I = 0.1;
+
+ // fudge factor! can't change passive membrane
+ // conductance from gbar0 = 0.001
+
+ double gbar0 = 0.001;
+ double f = gbar/gbar0;
+
+ // scale everything else
+ r_L *= f;
+ c_m /= f;
+ I /= f;
+
+ for (auto& seg: c.segments()) {
+ seg->add_mechanism(pas_parameters());
+ seg->mechanism("membrane").set("r_L", r_L);
+ seg->mechanism("membrane").set("c_m", c_m);
+ // seg->mechanism("pas").set("g", gbar);
+
+ if (seg->is_dendrite()) {
+ seg->set_compartments(4);
+ }
+ }
+
+ if (with_stim) {
+ // stimulus in the middle of our zero-volume 'soma'
+ // corresponds to proximal end of cable.
+ c.add_stimulus({0,0.5}, {0., math::infinity<>(), I});
+ }
+ return c;
+}
+
+
/*
* Attach voltage probes at each cable mid-point and end-point,
* and at soma mid-point.
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index b6a7d811f24a88fdd618ea80e2e2ac17d6e919f5..60ca6ab0572c74a580399364609d16b9c242adaa 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -12,6 +12,7 @@ set(TEST_SOURCES
test_cell_group.cpp
test_lexcmp.cpp
test_mask_stream.cpp
+ test_math.cpp
test_matrix.cpp
test_mechanisms.cpp
test_nop.cpp
diff --git a/tests/unit/test_math.cpp b/tests/unit/test_math.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b32bd194cb75b9df1c31c34c0323fdc0b95bce9e
--- /dev/null
+++ b/tests/unit/test_math.cpp
@@ -0,0 +1,36 @@
+#include <cmath>
+#include <limits>
+
+#include "gtest.h"
+#include <math.hpp>
+
+using namespace nest::mc::math;
+
+TEST(math, infinity) {
+ // check values for float, double, long double
+ auto finf = infinity<float>();
+ EXPECT_TRUE((std::is_same<float, decltype(finf)>::value));
+ EXPECT_TRUE(std::isinf(finf));
+ EXPECT_GT(finf, 0.f);
+
+ auto dinf = infinity<double>();
+ EXPECT_TRUE((std::is_same<double, decltype(dinf)>::value));
+ EXPECT_TRUE(std::isinf(dinf));
+ EXPECT_GT(dinf, 0.0);
+
+ auto ldinf = infinity<long double>();
+ EXPECT_TRUE((std::is_same<long double, decltype(ldinf)>::value));
+ EXPECT_TRUE(std::isinf(ldinf));
+ EXPECT_GT(ldinf, 0.0l);
+
+ // check default value promotes correctly (i.e., acts like INFINITY)
+ struct {
+ float f;
+ double d;
+ long double ld;
+ } check = {infinity<>(), infinity<>(), infinity<>()};
+
+ EXPECT_EQ(std::numeric_limits<float>::infinity(), check.f);
+ EXPECT_EQ(std::numeric_limits<double>::infinity(), check.d);
+ EXPECT_EQ(std::numeric_limits<long double>::infinity(), check.ld);
+}
diff --git a/tests/validation/CMakeLists.txt b/tests/validation/CMakeLists.txt
index 5f61ecc32af9c2083cc69f3945bf8bc8d58417db..783ef00321a34313ea36d90d17d9ec3862afa148 100644
--- a/tests/validation/CMakeLists.txt
+++ b/tests/validation/CMakeLists.txt
@@ -12,8 +12,9 @@ set(VALIDATION_SOURCES
validate.cpp
)
-add_definitions("-DDATADIR=\"${CMAKE_SOURCE_DIR}/data\"")
-
+if(VALIDATION_DATA_DIR)
+ add_definitions("-DDATADIR=\"${VALIDATION_DATA_DIR}\"")
+endif()
add_executable(validate.exe ${VALIDATION_SOURCES} ${HEADERS})
set(TARGETS validate.exe)
@@ -36,8 +37,8 @@ foreach(target ${TARGETS})
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
)
- if (BUILD_VALIDATION_DATA)
- add_dependencies(${target} validation_data)
+ if(BUILD_VALIDATION_DATA)
+ add_dependencies(${target} validation_data)
endif()
endforeach()
diff --git a/tests/validation/hh_soma.jl b/tests/validation/hh_soma.jl
deleted file mode 100644
index 82e7df8b5f6c2c594b8b35f584cbf41a42e58e6b..0000000000000000000000000000000000000000
--- a/tests/validation/hh_soma.jl
+++ /dev/null
@@ -1,132 +0,0 @@
-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/trace_analysis.cpp b/tests/validation/trace_analysis.cpp
index d641e34bcbad20ea1c64aece9d47f9c17a4856f0..a754cfe228b8893a6d557eb40092cbde8de4dfa5 100644
--- a/tests/validation/trace_analysis.cpp
+++ b/tests/validation/trace_analysis.cpp
@@ -83,7 +83,7 @@ util::optional<trace_peak> peak_delta(const trace_data& a, const trace_data& b)
auto p = local_maxima(a);
auto q = local_maxima(b);
- if (p.size()!=q.size()) return util::nothing;
+ if (p.size()!=q.size() || p.empty()) return util::nothing;
auto max_delta = p[0]-q[0];
diff --git a/tests/validation/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp
index 0a57fa840e52e70ea7dcaef4dc14c0f54ee86d93..e01eb4cf9c6f95a2631e1bbf370762967df71487 100644
--- a/tests/validation/validate_ball_and_stick.cpp
+++ b/tests/validation/validate_ball_and_stick.cpp
@@ -21,66 +21,75 @@
using namespace nest::mc;
-// TODO: further consolidate common code
+struct sampler_info {
+ const char* label;
+ cell_member_type probe;
+ simple_sampler sampler;
+};
+
+template <
+ typename lowered_cell,
+ typename SamplerInfoSeq
+>
+void run_ncomp_convergence_test(
+ const char* model_name,
+ const char* ref_data_path,
+ const cell& c,
+ SamplerInfoSeq& samplers,
+ float t_end=100.f,
+ float dt=0.001)
+{
+ using nlohmann::json;
+
+ SCOPED_TRACE(model_name);
-TEST(ball_and_stick, neuron_ref) {
- // compare voltages against reference data produced from
- // nrn/ball_and_stick.py
-
- using namespace nlohmann;
-
- using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
auto& V = g_trace_io;
-
bool verbose = V.verbose();
int max_ncomp = V.max_ncomp();
- // load validation data
- auto ref_data = V.load_traces("neuron_ball_and_stick.json");
- bool run_validation =
- ref_data.count("soma.mid") &&
- ref_data.count("dend.mid") &&
- ref_data.count("dend.end");
+ auto keys = util::transform_view(samplers,
+ [](const sampler_info& se) { return se.label; });
- EXPECT_TRUE(run_validation);
+ bool run_validation = false;
+ std::map<std::string, trace_data> ref_data;
+ try {
+ ref_data = V.load_traces(ref_data_path);
- // generate test data
- cell c = make_cell_ball_and_stick();
- add_common_voltage_probes(c);
+ run_validation = std::all_of(keys.begin(), keys.end(),
+ [&](const char* key) { return ref_data.count(key)>0; });
- float sample_dt = .025;
- std::pair<const char *, simple_sampler> samplers[] = {
- {"soma.mid", simple_sampler(sample_dt)},
- {"dend.mid", simple_sampler(sample_dt)},
- {"dend.end", simple_sampler(sample_dt)}
- };
+ EXPECT_TRUE(run_validation);
+ }
+ catch (std::runtime_error&) {
+ ADD_FAILURE() << "failure loading reference data: " << ref_data_path;
+ }
std::map<std::string, std::vector<conv_entry<int>>> conv_results;
for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- for (auto& se: samplers) {
- se.second.reset();
+ for (auto& seg: c.segments()) {
+ if (!seg->is_soma()) {
+ seg->set_compartments(ncomp);
+ }
}
- c.cable(1)->set_compartments(ncomp);
model<lowered_cell> m(singleton_recipe{c});
- // the ball-and-stick-cell (should) have three voltage probes:
- // centre of soma, centre of dendrite, end of dendrite.
-
- m.attach_sampler({0u, 0u}, samplers[0].second.sampler<>());
- m.attach_sampler({0u, 1u}, samplers[1].second.sampler<>());
- m.attach_sampler({0u, 2u}, samplers[2].second.sampler<>());
+ // reset samplers and attach to probe locations
+ for (auto& se: samplers) {
+ se.sampler.reset();
+ m.attach_sampler(se.probe, se.sampler.template sampler<>());
+ }
- m.run(100, 0.001);
+ m.run(t_end, dt);
for (auto& se: samplers) {
- std::string key = se.first;
- const simple_sampler& s = se.second;
+ std::string key = se.label;
+ const simple_sampler& s = se.sampler;
// save trace
json meta = {
{"name", "membrane voltage"},
- {"model", "ball_and_stick"},
+ {"model", model_name},
{"sim", "nestmc"},
{"ncomp", ncomp},
{"units", "mV"}};
@@ -101,7 +110,7 @@ TEST(ball_and_stick, neuron_ref) {
report_conv_table(std::cout, conv_results, "ncomp");
}
- for (auto key: util::transform_view(samplers, util::first)) {
+ for (auto key: keys) {
SCOPED_TRACE(key);
const auto& results = conv_results[key];
@@ -109,188 +118,92 @@ TEST(ball_and_stick, neuron_ref) {
}
}
-TEST(ball_and_taper, neuron_ref) {
- // compare voltages against reference data produced from
- // nrn/ball_and_taper.py
-
- using namespace nlohmann;
-
+TEST(ball_and_stick, neuron_ref) {
using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
- auto& V = g_trace_io;
-
- bool verbose = V.verbose();
- int max_ncomp = V.max_ncomp();
- // load validation data
- auto ref_data = V.load_traces("neuron_ball_and_taper.json");
- bool run_validation =
- ref_data.count("soma.mid") &&
- ref_data.count("taper.mid") &&
- ref_data.count("taper.end");
-
- EXPECT_TRUE(run_validation);
-
- // generate test data
- cell c = make_cell_ball_and_taper();
+ cell c = make_cell_ball_and_stick();
add_common_voltage_probes(c);
- float sample_dt = .025;
- std::pair<const char *, simple_sampler> samplers[] = {
- {"soma.mid", simple_sampler(sample_dt)},
- {"taper.mid", simple_sampler(sample_dt)},
- {"taper.end", simple_sampler(sample_dt)}
+ float sample_dt = 0.025f;
+ sampler_info samplers[] = {
+ {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
+ {"dend.mid", {0u, 1u}, simple_sampler(sample_dt)},
+ {"dend.end", {0u, 2u}, simple_sampler(sample_dt)}
};
- std::map<std::string, std::vector<conv_entry<int>>> conv_results;
-
- for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- for (auto& se: samplers) {
- se.second.reset();
- }
- c.cable(1)->set_compartments(ncomp);
- model<lowered_cell> m(singleton_recipe{c});
-
- // the ball-and-stick-cell (should) have three voltage probes:
- // centre of soma, centre of dendrite, end of dendrite.
-
- m.attach_sampler({0u, 0u}, samplers[0].second.sampler<>());
- m.attach_sampler({0u, 1u}, samplers[1].second.sampler<>());
- m.attach_sampler({0u, 2u}, samplers[2].second.sampler<>());
-
- m.run(100, 0.001);
-
- for (auto& se: samplers) {
- std::string key = se.first;
- const simple_sampler& s = se.second;
-
- // save trace
- json meta = {
- {"name", "membrane voltage"},
- {"model", "ball_and_taper"},
- {"sim", "nestmc"},
- {"ncomp", ncomp},
- {"units", "mV"}};
-
- V.save_trace(key, s.trace, meta);
-
- // compute metrics
- if (run_validation) {
- double linf = linf_distance(s.trace, ref_data[key]);
- auto pd = peak_delta(s.trace, ref_data[key]);
+ run_ncomp_convergence_test<lowered_cell>(
+ "ball_and_stick",
+ "neuron_ball_and_stick.json",
+ c,
+ samplers);
+}
- conv_results[key].push_back({key, ncomp, linf, pd});
- }
- }
- }
+TEST(ball_and_taper, neuron_ref) {
+ using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
- if (verbose && run_validation) {
- report_conv_table(std::cout, conv_results, "ncomp");
- }
+ cell c = make_cell_ball_and_taper();
+ add_common_voltage_probes(c);
- for (auto key: util::transform_view(samplers, util::first)) {
- SCOPED_TRACE(key);
+ float sample_dt = 0.025f;
+ sampler_info samplers[] = {
+ {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
+ {"taper.mid", {0u, 1u}, simple_sampler(sample_dt)},
+ {"taper.end", {0u, 2u}, simple_sampler(sample_dt)}
+ };
- const auto& results = conv_results[key];
- assert_convergence(results);
- }
+ run_ncomp_convergence_test<lowered_cell>(
+ "ball_and_taper",
+ "neuron_ball_and_taper.json",
+ c,
+ samplers);
}
-
TEST(ball_and_3stick, neuron_ref) {
- // compare voltages against reference data produced from
- // nrn/ball_and_3stick.py
-
- using namespace nlohmann;
-
using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
- auto& V = g_trace_io;
-
- bool verbose = V.verbose();
- int max_ncomp = V.max_ncomp();
-
- // load validation data
- auto ref_data = V.load_traces("neuron_ball_and_3stick.json");
- bool run_validation =
- ref_data.count("soma.mid") &&
- ref_data.count("dend1.mid") &&
- ref_data.count("dend1.end") &&
- ref_data.count("dend2.mid") &&
- ref_data.count("dend2.end") &&
- ref_data.count("dend3.mid") &&
- ref_data.count("dend3.end");
-
- EXPECT_TRUE(run_validation);
-
- // generate test data
cell c = make_cell_ball_and_3stick();
add_common_voltage_probes(c);
- float sample_dt = .025;
- std::pair<const char *, simple_sampler> samplers[] = {
- {"soma.mid", simple_sampler(sample_dt)},
- {"dend1.mid", simple_sampler(sample_dt)},
- {"dend1.end", simple_sampler(sample_dt)},
- {"dend2.mid", simple_sampler(sample_dt)},
- {"dend2.end", simple_sampler(sample_dt)},
- {"dend3.mid", simple_sampler(sample_dt)},
- {"dend3.end", simple_sampler(sample_dt)}
+ float sample_dt = 0.025f;
+ sampler_info samplers[] = {
+ {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
+ {"dend1.mid", {0u, 1u}, simple_sampler(sample_dt)},
+ {"dend1.end", {0u, 2u}, simple_sampler(sample_dt)},
+ {"dend2.mid", {0u, 3u}, simple_sampler(sample_dt)},
+ {"dend2.end", {0u, 4u}, simple_sampler(sample_dt)},
+ {"dend3.mid", {0u, 5u}, simple_sampler(sample_dt)},
+ {"dend3.end", {0u, 6u}, simple_sampler(sample_dt)}
};
- std::map<std::string, std::vector<conv_entry<int>>> conv_results;
-
- for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- for (auto& se: samplers) {
- se.second.reset();
- }
- c.cable(1)->set_compartments(ncomp);
- c.cable(2)->set_compartments(ncomp);
- c.cable(3)->set_compartments(ncomp);
- model<lowered_cell> m(singleton_recipe{c});
-
- // the ball-and-3stick-cell (should) have seven voltage probes:
- // centre of soma, followed by centre of section, end of section
- // for each of the three dendrite sections.
-
- for (unsigned i = 0; i < util::size(samplers); ++i) {
- m.attach_sampler({0u, i}, samplers[i].second.sampler<>());
- }
-
- m.run(100, 0.001);
-
- for (auto& se: samplers) {
- std::string key = se.first;
- const simple_sampler& s = se.second;
-
- // save trace
- json meta = {
- {"name", "membrane voltage"},
- {"model", "ball_and_3stick"},
- {"sim", "nestmc"},
- {"ncomp", ncomp},
- {"units", "mV"}};
-
- V.save_trace(key, s.trace, meta);
+ run_ncomp_convergence_test<lowered_cell>(
+ "ball_and_3stick",
+ "neuron_ball_and_3stick.json",
+ c,
+ samplers);
+}
- // compute metrics
- if (run_validation) {
- double linf = linf_distance(s.trace, ref_data[key]);
- auto pd = peak_delta(s.trace, ref_data[key]);
+TEST(rallpack1, numeric_ref) {
+ using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>;
- conv_results[key].push_back({key, ncomp, linf, pd});
- }
- }
- }
+ cell c = make_cell_simple_cable();
- if (verbose && run_validation) {
- report_conv_table(std::cout, conv_results, "ncomp");
- }
+ // three probes: left end, 30% along, right end.
+ c.add_probe({{1, 0.0}, probeKind::membrane_voltage});
+ c.add_probe({{1, 0.3}, probeKind::membrane_voltage});
+ c.add_probe({{1, 1.0}, probeKind::membrane_voltage});
- for (auto key: util::transform_view(samplers, util::first)) {
- SCOPED_TRACE(key);
+ float sample_dt = 0.025f;
+ sampler_info samplers[] = {
+ {"cable.x0.0", {0u, 0u}, simple_sampler(sample_dt)},
+ {"cable.x0.3", {0u, 1u}, simple_sampler(sample_dt)},
+ {"cable.x1.0", {0u, 2u}, simple_sampler(sample_dt)},
+ };
- const auto& results = conv_results[key];
- assert_convergence(results);
- }
+ run_ncomp_convergence_test<lowered_cell>(
+ "rallpack1",
+ "numeric_rallpack1.json",
+ c,
+ samplers,
+ 250.f);
}
diff --git a/tests/validation/validate_soma.cpp b/tests/validation/validate_soma.cpp
index 26d919b07c10a5966734c0e12dd7a7296e55f139..31eb784add73f5b675e8b2943e8faa453a0afb26 100644
--- a/tests/validation/validate_soma.cpp
+++ b/tests/validation/validate_soma.cpp
@@ -20,7 +20,7 @@
using namespace nest::mc;
-TEST(soma, neuron_ref) {
+TEST(soma, numeric_ref) {
// compare voltages against reference data produced from
// nrn/ball_and_taper.py
@@ -32,10 +32,21 @@ TEST(soma, neuron_ref) {
bool verbose = V.verbose();
// load validation data
- auto ref_data = V.load_traces("neuron_soma.json");
+
+ bool run_validation = false;
+ std::map<std::string, trace_data> ref_data;
const char* key = "soma.mid";
- bool run_validation = ref_data.count(key);
- EXPECT_TRUE(run_validation);
+
+ const char* ref_data_path = "numeric_soma.json";
+ try {
+ ref_data = V.load_traces(ref_data_path);
+ run_validation = ref_data.count(key);
+
+ EXPECT_TRUE(run_validation);
+ }
+ catch (std::runtime_error&) {
+ ADD_FAILURE() << "failure loading reference data: " << ref_data_path;
+ }
// generate test data
cell c = make_cell_soma_only();
diff --git a/tests/validation/validation_data.hpp b/tests/validation/validation_data.hpp
index 70e99a52227b70f9772dd94370ec38fb7301e2d1..254d6eeed77cd1d99c3ba31a924a8db361e7ff21 100644
--- a/tests/validation/validation_data.hpp
+++ b/tests/validation/validation_data.hpp
@@ -70,7 +70,7 @@ public:
}
private:
- util::path datadir_ = DATADIR "/validation";
+ util::path datadir_ = DATADIR;
std::ofstream out_;
nlohmann::json jtraces_ = nlohmann::json::array();
bool verbose_flag_ = false;
diff --git a/validation/CMakeLists.txt b/validation/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..99b4bec969a87219693dfe0bbd05591e1acfd6c2
--- /dev/null
+++ b/validation/CMakeLists.txt
@@ -0,0 +1,51 @@
+# Validation data generation
+
+add_custom_target(validation_data)
+
+# Helper function because ffs CMake.
+
+function(make_unique_target_name name path)
+ # Try and make a broadly human readable target name if possible.
+ string(REGEX REPLACE ".*/" "" leaf "${path}")
+ string(REGEX REPLACE "[^a-zA-Z0-9_.+-]" "_" canon "${leaf}")
+
+ # Check against reserved names, of which of course there is no documented list
+ # except in the sodding CMake source code. Seriously. Look at the CMP0037 policy
+ # text for a laugh.
+ if(canon MATCHES "^(all|ALL_BUILD|help|install|INSTALL|preinstall|clean|edit_cache|rebuild_cache|test|RUN_TESTS|package|PACKAGE|package_source|ZERO_CHECK)$")
+ set(canon "${canon}_")
+ endif()
+ while((TARGET "${canon}"))
+ set(canon "${canon}_")
+ endwhile()
+ set("${name}" "${canon}" PARENT_SCOPE)
+endfunction()
+
+# Helper function to add a data generation script that writes to standard output.
+# e.g.:
+# add_validation_data(OUTPUT foo_model.json DEPENDS foo_model.py common.py COMMAND python foo_model.py)
+
+include(CMakeParseArguments)
+function(add_validation_data)
+ cmake_parse_arguments(ADD_VALIDATION_DATA "" "OUTPUT" "DEPENDS;COMMAND" ${ARGN})
+ set(out "${VALIDATION_DATA_DIR}/${ADD_VALIDATION_DATA_OUTPUT}")
+ string(REGEX REPLACE "([^;]+)" "${CMAKE_CURRENT_SOURCE_DIR}/\\1" deps "${ADD_VALIDATION_DATA_DEPENDS}")
+ add_custom_command(
+ OUTPUT "${out}"
+ DEPENDS ${deps}
+ WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
+ COMMAND ${ADD_VALIDATION_DATA_COMMAND} > "${out}")
+
+ # Cmake, why can't we just write add_dependencies(validation_data "${out}")?!
+ make_unique_target_name(ffs_cmake "${out}")
+ add_custom_target("${ffs_cmake}" DEPENDS "${out}")
+ add_dependencies(validation_data "${ffs_cmake}")
+endfunction()
+
+
+if(BUILD_NRN_VALIDATION_DATA)
+ add_subdirectory(ref/neuron)
+endif()
+
+add_subdirectory(ref/numeric)
+
diff --git a/validation/README.md b/validation/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..2357501d57fbe0101a05b1ef077c60eafde22a16
--- /dev/null
+++ b/validation/README.md
@@ -0,0 +1,22 @@
+# Validation data and generation
+
+## Sub-directory organization
+
+`validation/data`
+ ~ Generated validation data
+
+`validation/ref`
+ ~ Reference models
+
+`validation/ref/neuron`
+ ~ NEURON-based reference models, run with `nrniv -python`
+
+`validation/ref/numeric`
+ ~ Direct numerical and analytic models
+
+## Data generation
+
+Data is generated via the `validation_data` CMake target, which is
+a prerequisite for the `validation.exe` test executable.
+
+
diff --git a/data/validation/.keep b/validation/data/.keep
similarity index 100%
rename from data/validation/.keep
rename to validation/data/.keep
diff --git a/validation/ref/neuron/CMakeLists.txt b/validation/ref/neuron/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4df17dc38d23b59c935dad9e02ca0eb19d2b793f
--- /dev/null
+++ b/validation/ref/neuron/CMakeLists.txt
@@ -0,0 +1,18 @@
+# note: function add_validation_data defined in validation/CMakeLists.txt
+
+set(models
+ ball_and_stick
+ ball_and_3stick
+ ball_and_taper
+ simple_exp_synapse
+ simple_exp2_synapse
+ soma)
+
+foreach(model ${models})
+ set(script "${model}.py")
+ add_validation_data(
+ OUTPUT "neuron_${model}.json"
+ DEPENDS "${script}" "nrn_validation.py"
+ COMMAND ${NRNIV_BIN} -nobanner -python "${script}")
+endforeach()
+
diff --git a/nrn/ball_and_3stick.py b/validation/ref/neuron/ball_and_3stick.py
similarity index 100%
rename from nrn/ball_and_3stick.py
rename to validation/ref/neuron/ball_and_3stick.py
diff --git a/nrn/ball_and_stick.py b/validation/ref/neuron/ball_and_stick.py
similarity index 100%
rename from nrn/ball_and_stick.py
rename to validation/ref/neuron/ball_and_stick.py
diff --git a/nrn/ball_and_taper.py b/validation/ref/neuron/ball_and_taper.py
similarity index 100%
rename from nrn/ball_and_taper.py
rename to validation/ref/neuron/ball_and_taper.py
diff --git a/nrn/generate_validation.sh b/validation/ref/neuron/generate_validation.sh
similarity index 100%
rename from nrn/generate_validation.sh
rename to validation/ref/neuron/generate_validation.sh
diff --git a/nrn/nrn_validation.py b/validation/ref/neuron/nrn_validation.py
similarity index 99%
rename from nrn/nrn_validation.py
rename to validation/ref/neuron/nrn_validation.py
index 66563885b05ba9212bd017ea1c90b514c93c6ff2..a5ffa2d2e238148d68dee137dd92daa17ce89d89 100644
--- a/nrn/nrn_validation.py
+++ b/validation/ref/neuron/nrn_validation.py
@@ -33,7 +33,7 @@ default_model_parameters = {
'g_pas': 0.001, # Passive conductance in S/cm^2
'e_pas': -65.0, # Leak reversal potential in mV
'Ra': 100.0, # Intracellular resistivity in Ω·cm
- 'cm': 1.0, # Membrane areal capacitance in µF/cm2
+ 'cm': 1.0, # Membrane areal capacitance in µF/cm^2
'tau': 2.0, # Exponential synapse time constant
'tau1': 0.5, # Exp2 synapse tau1
'tau2': 2.0, # Exp2 synapse tau2
diff --git a/nrn/simple_exp2_synapse.py b/validation/ref/neuron/simple_exp2_synapse.py
similarity index 100%
rename from nrn/simple_exp2_synapse.py
rename to validation/ref/neuron/simple_exp2_synapse.py
diff --git a/nrn/simple_exp_synapse.py b/validation/ref/neuron/simple_exp_synapse.py
similarity index 100%
rename from nrn/simple_exp_synapse.py
rename to validation/ref/neuron/simple_exp_synapse.py
diff --git a/nrn/soma.py b/validation/ref/neuron/soma.py
similarity index 75%
rename from nrn/soma.py
rename to validation/ref/neuron/soma.py
index 867fa7cc7b82e2dabb28d4caa03ca88ca10be3f2..706af0f559a3172284bb45d9ca503e47a9f63cc2 100644
--- a/nrn/soma.py
+++ b/validation/ref/neuron/soma.py
@@ -6,11 +6,9 @@ import nrn_validation as V
V.override_defaults_from_args()
-# dendrite geometry: all 100 µm long, 1 µm diameter.
-geom = [(0,1), (100, 1)]
-
model = V.VModel()
-model.add_soma(18.8, Ra=123)
+
+model.add_soma(18.8, Ra=100)
model.add_iclamp(10, 100, 0.1)
# NB: this doesn't seem to have converged with
diff --git a/validation/ref/numeric/CMakeLists.txt b/validation/ref/numeric/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ddb5001fb9888b6a207043219db6a9bd11de2c69
--- /dev/null
+++ b/validation/ref/numeric/CMakeLists.txt
@@ -0,0 +1,13 @@
+# note: function add_validation_data defined in validation/CMakeLists.txt
+
+if(BUILD_JULIA_VALIDATION_DATA)
+ add_validation_data(
+ OUTPUT numeric_soma.json
+ DEPENDS numeric_soma.jl HHChannels.jl
+ COMMAND ${JULIA_BIN} numeric_soma.jl)
+
+ add_validation_data(
+ OUTPUT numeric_rallpack1.json
+ DEPENDS numeric_rallpack1.jl PassiveCable.jl
+ COMMAND ${JULIA_BIN} numeric_rallpack1.jl)
+endif()
diff --git a/validation/ref/numeric/HHChannels.jl b/validation/ref/numeric/HHChannels.jl
new file mode 100644
index 0000000000000000000000000000000000000000..bf32c8e84ddafb0f32f23c94086ab23c1408d376
--- /dev/null
+++ b/validation/ref/numeric/HHChannels.jl
@@ -0,0 +1,143 @@
+module HHChannels
+
+export Stim, run_hh
+
+using Sundials
+using SIUnits.ShortUnits
+
+immutable HHParam
+ c_m # membrane spacific capacitance
+ gnabar # Na channel cross-membrane conductivity
+ gkbar # K channel cross-membrane conductivity
+ gl # Leak conductivity
+ ena # Na channel reversal potential
+ ek # K channel reversal potential
+ el # Leak reversal potential
+ q10 # temperature dependent rate coefficient
+ # (= 3^((T-T₀)/10K) with T₀ = 6.3 °C)
+
+ # constructor with default values, corresponding
+ # to a resting potential of -65 mV and temperature 6.3 °C
+ HHParam(;
+ #c_m = 0.01nF*m^-2,
+ c_m = 0.01F*m^-2,
+ gnabar = .12S*cm^-2,
+ ena = 115.0mV + -65.0mV,
+ gkbar = .036S*cm^-2,
+ ek = -12.0mV + -65.0mV,
+ gl = .0003S*cm^-2,
+ el = -54.3mV,
+ q10 = 1
+ ) = new(c_m, gnabar, gkbar, gl, ena, ek, el, q10)
+
+end
+
+immutable Stim
+ t0 # start time of stimulus
+ t1 # stop time of stimulus
+ i_e # stimulus current density
+
+ Stim() = new(0s, 0s, 0A/m^2)
+ Stim(t0, t1, i_e) = new(t0, t1, i_e)
+end
+
+vtrap(x,y) = x/(exp(x/y) - 1.0)
+
+# "m" sodium activation system
+function m_lims(v, q10)
+ 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, q10)
+ 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, q10)
+ 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
+
+# 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, q10)
+ mtau, minf = m_lims(v, q10)
+ htau, hinf = h_lims(v, q10)
+ ntau, ninf = n_lims(v, q10)
+
+ return (v, minf, hinf, ninf)
+end
+
+# Given time t and state (v, m, h, n),
+# return (vdot, mdot, hdot, ndot)
+function f(t, state; p=HHParam(), stim=Stim())
+ v, m, h, n = state
+
+ # calculate current density due to ion channels
+ gna = p.gnabar*m*m*m*h
+ gk = p.gkbar*n*n*n*n
+
+
+ ina = gna*(v - p.ena)
+ ik = gk*(v - p.ek)
+ il = p.gl*(v - p.el)
+
+ itot = ik + ina + il
+
+ # calculate current density due to stimulus
+ if t>=stim.t0 && t<stim.t1
+ itot -= stim.i_e
+ end
+
+ # calculate the voltage dependent rates for the gating variables
+ mtau, minf = m_lims(v, p.q10)
+ htau, hinf = h_lims(v, p.q10)
+ ntau, ninf = n_lims(v, p.q10)
+
+ return (-itot/p.c_m, (minf-m)/mtau, (hinf-h)/htau, (ninf-n)/ntau)
+end
+
+function run_hh(t_end; v0=-65mV, stim=Stim(), param=HHParam(), sample_dt=0.01ms)
+ v_scale = 1V
+ t_scale = 1s
+
+ v0, m0, h0, n0 = initial_conditions(v0, param.q10)
+ y0 = [ v0/v_scale, m0, h0, n0 ]
+
+ samples = collect(0s: sample_dt: t_end)
+
+ fbis(t, y, ydot) = begin
+ vdot, mdot, hdot, ndot =
+ f(t*t_scale, (y[1]*v_scale, y[2], y[3], y[4]), stim=stim, p=param)
+
+ ydot[1], ydot[2], ydot[3], ydot[4] =
+ vdot*t_scale/v_scale, mdot*t_scale, hdot*t_scale, ndot*t_scale
+
+ return Sundials.CV_SUCCESS
+ end
+
+ # Ideally would run with vector absolute tolerance to account for v_scale,
+ # but this would prevent us using the nice cvode wrapper.
+
+ res = Sundials.cvode(fbis, y0, map(t->t/t_scale, samples), abstol=1e-6, reltol=5e-10)
+
+ # Use map here because of issues with type deduction with arrays and SIUnits.
+ return samples, map(v->v*v_scale, res[:, 1])
+end
+
+end # module HHChannels
diff --git a/scripts/PassiveCable.jl b/validation/ref/numeric/PassiveCable.jl
similarity index 91%
rename from scripts/PassiveCable.jl
rename to validation/ref/numeric/PassiveCable.jl
index 76450490ce5c4f265403edd35f2bedad4ad7c3bb..6928b74c33948be7162885040cb89870ae245596 100644
--- a/scripts/PassiveCable.jl
+++ b/validation/ref/numeric/PassiveCable.jl
@@ -18,26 +18,29 @@ export cable_normalize, cable, rallpack1
#
# Return:
# g(x, t)
+#
+# TODO: verify correctness when L≠1
-function cable_normalized(x, t, L; tol=1e-8)
+function cable_normalized(x::Float64, t::Float64, L::Float64; tol=1e-8)
if t<=0
return 0.0
else
ginf = -cosh(L-x)/sinh(L)
sum = exp(-t/L)
+ Ltol = L*tol
for k = countfrom(1)
a = k*pi/L
- e = exp(-t*(1+a^2))
+ b = exp(-t*(1+a^2))
- sum += 2/L*e*cos(a*x)/(1+a^2)
- resid_ub = e/(L*a^3*t)
+ sum += 2*b*cos(a*x)/(1+a^2)
+ resid_ub = b/(a^3*t)
- if resid_ub<tol
+ if resid_ub<Ltol
break
end
end
- return ginf+sum
+ return ginf+sum/L;
end
end
diff --git a/validation/ref/numeric/numeric_rallpack1.jl b/validation/ref/numeric/numeric_rallpack1.jl
new file mode 100644
index 0000000000000000000000000000000000000000..ea9755aa682ced365d2e22fd6a59cd7456c41682
--- /dev/null
+++ b/validation/ref/numeric/numeric_rallpack1.jl
@@ -0,0 +1,65 @@
+#!/usr/bin/env julia
+
+include("PassiveCable.jl")
+
+using JSON
+using SIUnits.ShortUnits
+using PassiveCable
+
+# This should run the same effective model
+# as rallpack1, but with differing
+# electrical parameters (see below).
+
+function run_cable(x_prop, ts)
+ # Physical properties:
+
+ # f is a fudge factor. rM needs to be the same
+ # the same as in nestmc, where we cannot yet set
+ # the membrane conductance parameter. Scaling
+ # other parameters proportionally, however,
+ # gives the same dynamics.
+
+ f = 0.1/4
+
+ diam = 1.0µm # cable diameter
+ L = 1.0mm # cable length
+ I = 0.1nA /f # current injection
+ rL = 1.0Ω*m *f # bulk resistivity
+ erev = -65.0mV # (passive) reversal potential
+ rM = 4Ω*m^2 *f # membrane resistivity
+ cM = 0.01F/m^2 /f # membrane specific capacitance
+
+ # convert to linear resistivity, length and time constants
+ area = pi*diam^2/4
+ r = rL/area
+
+ lambda = sqrt(diam/4 * rM/rL)
+ tau = cM*rM
+
+ # compute solutions
+ tol = 1e-8mV
+ return [cable(L*x_prop, t, L, lambda, tau, r, erev, -I, tol=tol) for t in ts]
+end
+
+function run_rallpack1(x_prop, ts)
+ return [rallpack1(0.001*x_prop, t/s)*V for t in ts]
+end
+
+# Generate traces at x=0, x=0.3L, x=L
+
+ts = collect(0s: 0.025ms: 250ms)
+trace = Dict(
+ :name => "membrane voltage",
+ :sim => "numeric",
+ :model => "rallpack1",
+ :units => "mV",
+ :data => Dict(
+ :time => map(t->t/ms, ts),
+ symbol("cable.x0.0") => map(v->v/mV, run_cable(0, ts)),
+ symbol("cable.x0.3") => map(v->v/mV, run_cable(0.3, ts)),
+ symbol("cable.x1.0") => map(v->v/mV, run_cable(1.0, ts))
+ )
+)
+
+println(JSON.json([trace]))
+
diff --git a/validation/ref/numeric/numeric_soma.jl b/validation/ref/numeric/numeric_soma.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6c431f845cbf5acbb3a1825e8d52951482d255a1
--- /dev/null
+++ b/validation/ref/numeric/numeric_soma.jl
@@ -0,0 +1,27 @@
+#!/usr/bin/env julia
+
+include("HHChannels.jl")
+
+using JSON
+using SIUnits.ShortUnits
+using HHChannels
+
+radius = 18.8µm/2
+area = 4*pi*radius^2
+
+stim = Stim(10ms, 100ms, 0.1nA/area)
+ts, vs = run_hh(100ms, stim=stim, sample_dt=0.025ms)
+
+trace = Dict(
+ :name => "membrane voltage",
+ :sim => "numeric",
+ :model => "soma",
+ :units => "mV",
+ :data => Dict(
+ :time => map(t->t/ms, ts),
+ symbol("soma.mid") => map(v->v/mV, vs)
+ )
+)
+
+println(JSON.json([trace]))
+