diff --git a/.gitignore b/.gitignore
index e9d54cb796031fdcbba7935bec51511be3f6de6b..b710bef2b72c6b9f6eea87de726bf180c4f3a107 100644
--- a/.gitignore
+++ b/.gitignore
@@ -66,3 +66,5 @@ mechanisms/*.hpp
# build path
build
+
+commit.msg
diff --git a/.gitmodules b/.gitmodules
index ae780729fd3c2e22a3a6f5029fb0b611977ccfeb..86e0c20fe68579e7b0f1ccdd1dba498cf16b13fe 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -7,3 +7,6 @@
[submodule "modparser"]
path = external/modparser
url = git@github.com:eth-cscs/modparser.git
+[submodule "external/tclap"]
+ path = external/tclap
+ url = https://github.com/eile/tclap.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a432f4e8d3cadf2899f90b13ff687018e4fecd0d..44e40b639c2328b5512ad4c557a877d4d9cbec8f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,6 +25,12 @@ if(WITH_ASSERTIONS)
add_definitions("-DWITH_ASSERTIONS")
endif()
+# enable traces?
+set(WITH_TRACE OFF CACHE BOOL "enable TRACE() macros in code")
+if(WITH_TRACE)
+ add_definitions("-DWITH_TRACE")
+endif()
+
# TBB support
set(WITH_TBB OFF CACHE BOOL "use TBB for on-node threading" )
if(WITH_TBB)
@@ -56,6 +62,12 @@ if(SYSTEM_CRAY)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -dynamic")
endif()
+# Cray systems
+set(TARGET_KNL OFF CACHE BOOL "target Intel KNL architecture")
+if(TARGET_KNL)
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXXOPT_KNL}")
+endif()
+
# targets for extermal dependencies
include(ExternalProject)
externalproject_add(modparser
@@ -71,6 +83,7 @@ externalproject_add(modparser
include_directories(${CMAKE_SOURCE_DIR}/external)
+include_directories(${CMAKE_SOURCE_DIR}/external/tclap/include)
include_directories(${CMAKE_SOURCE_DIR}/include)
include_directories(${CMAKE_SOURCE_DIR}/src)
include_directories(${CMAKE_SOURCE_DIR}/miniapp)
diff --git a/cmake/CompilerOptions.cmake b/cmake/CompilerOptions.cmake
index e5447d88e7d0009f8ec697bb1463e3a7b737b0d8..6dc09b961cc6a71ce69b1bf9cba5d6821e6949e8 100644
--- a/cmake/CompilerOptions.cmake
+++ b/cmake/CompilerOptions.cmake
@@ -12,3 +12,22 @@ if(${CMAKE_CXX_COMPILER_ID} MATCHES "Clang")
set(CXXOPT_WALL "${CXXOPT_WALL} -Wno-missing-braces")
endif()
+
+if(${CMAKE_CXX_COMPILER_ID} MATCHES "GNU")
+
+ # compiler flags for generating KNL-specific AVX512 instructions
+ # supported in gcc 4.9.x and later
+ #set(CXXOPT_KNL "-mavx512f")
+ set(CXXOPT_KNL "-march=knl")
+endif()
+
+if(${CMAKE_CXX_COMPILER_ID} MATCHES "Intel")
+ # Disable warning for unused template parameter
+ # this is raised by a templated function in the json library
+
+ set(CXXOPT_WALL "${CXXOPT_WALL} -wd488")
+
+ # compiler flags for generating KNL-specific AVX512 instructions
+ set(CXXOPT_KNL "-xMIC-AVX512")
+endif()
+
diff --git a/cmake/FindTBB.cmake b/cmake/FindTBB.cmake
index 7cb10d64935022bf78e4b5d8ffcac3328ebe4ac5..61caa6644d2ff56eee9502819034790601a7f265 100644
--- a/cmake/FindTBB.cmake
+++ b/cmake/FindTBB.cmake
@@ -178,7 +178,7 @@ if(NOT TBB_FOUND)
# Search for the libraries
find_library(TBB_${_comp}_LIBRARY_RELEASE ${_comp}
HINTS ${TBB_LIBRARY} ${TBB_SEARCH_DIR}
- PATHS ${TBB_DEFAULT_SEARCH_DIR}
+ PATHS ${TBB_DEFAULT_SEARCH_DIR} ENV LIBRARY_PATH
PATH_SUFFIXES "${TBB_LIB_PATH_SUFFIX}")
find_library(TBB_${_comp}_LIBRARY_DEBUG ${_comp}_debug
diff --git a/commit.msg b/commit.msg
deleted file mode 100644
index 3abb1935698a6efe63d5ea651d8c7d5c017541fe..0000000000000000000000000000000000000000
--- a/commit.msg
+++ /dev/null
@@ -1,9 +0,0 @@
-* change `probe_sort` enum to scoped enum
- - renamed to `probeKind`
- - refactored out of class
- - updated coding guidelines wiki with enum rules
-* refactor `std::pair` into structs with meaningfull name types in `cell.hpp`
- - not just probes: stimulii and detectors too.
-* add profiling region around sampling in `cell_group`
-* change output data format for traces to json
-* remove white space at end of lines (looking at you Sam)
diff --git a/data/validation/ball_and_3stick.json b/data/validation/ball_and_3stick.json
new file mode 100644
index 0000000000000000000000000000000000000000..e4121f4383e44dd7a7c15c36598b2a85fba728f2
--- /dev/null
+++ b/data/validation/ball_and_3stick.json
@@ -0,0 +1,158 @@
+[
+ {
+ "nseg": 5,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 6.943999999998871,
+ 20.48399999999761,
+ 33.862999999979664,
+ 52.52400000006877,
+ 65.9560000001329,
+ 79.3600000001969
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.1399999999987624,
+ 20.448999999997692,
+ 33.722999999978995,
+ 52.61200000006919,
+ 65.85800000013244,
+ 79.18600000019607
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 7.055999999998809,
+ 20.44199999999771,
+ 33.73499999997905,
+ 52.56900000006898,
+ 65.86400000013246,
+ 79.20200000019615
+ ],
+ "thresh": -20.0
+ }
+ }
+ },
+ {
+ "nseg": 11,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 6.948999999998868,
+ 20.506999999997557,
+ 33.91399999997991,
+ 52.528000000068786,
+ 65.98500000013304,
+ 79.42000000019719
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.147999999998758,
+ 20.478999999997622,
+ 33.77899999997926,
+ 52.62300000006924,
+ 65.8930000001326,
+ 79.24900000019638
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 7.061999999998806,
+ 20.467999999997648,
+ 33.785999999979296,
+ 52.576000000069016,
+ 65.8940000001326,
+ 79.26000000019643
+ ],
+ "thresh": -20.0
+ }
+ }
+ },
+ {
+ "nseg": 51,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 6.949999999998868,
+ 20.512999999997543,
+ 33.92699999997997,
+ 52.530000000068796,
+ 65.99200000013307,
+ 79.43500000019726
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.150999999998756,
+ 20.486999999997604,
+ 33.793999999979334,
+ 52.626000000069254,
+ 65.90200000013265,
+ 79.26500000019645
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 7.062999999998805,
+ 20.473999999997634,
+ 33.79899999997936,
+ 52.578000000069025,
+ 65.90200000013265,
+ 79.2740000001965
+ ],
+ "thresh": -20.0
+ }
+ }
+ },
+ {
+ "nseg": 101,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 6.949999999998868,
+ 20.512999999997543,
+ 33.927999999979974,
+ 52.530000000068796,
+ 65.99300000013308,
+ 79.43500000019726
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.150999999998756,
+ 20.486999999997604,
+ 33.793999999979334,
+ 52.626000000069254,
+ 65.90300000013265,
+ 79.26600000019646
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 7.062999999998805,
+ 20.47499999999763,
+ 33.79999999997936,
+ 52.578000000069025,
+ 65.90200000013265,
+ 79.2750000001965
+ ],
+ "thresh": -20.0
+ }
+ }
+ }
+]
\ No newline at end of file
diff --git a/data/validation/ball_and_stick.json b/data/validation/ball_and_stick.json
new file mode 100644
index 0000000000000000000000000000000000000000..5b459a3a287e012def1539866531efa458f66cc0
--- /dev/null
+++ b/data/validation/ball_and_stick.json
@@ -0,0 +1,158 @@
+[
+ {
+ "nseg": 5,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 6.980999999998851,
+ 20.98699999999644,
+ 34.792999999984104,
+ 48.61100000005008,
+ 62.43100000011607,
+ 76.25100000018206
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.2499999999987015,
+ 21.193999999995956,
+ 34.97599999998498,
+ 48.78900000005093,
+ 62.60700000011691,
+ 76.4270000001829
+ ],
+ "thresh": 10.0
+ },
+ "dend": {
+ "spikes": [
+ 7.169999999998746,
+ 21.177999999995993,
+ 34.97799999998499,
+ 48.79400000005096,
+ 62.614000000116945,
+ 76.43400000018293
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 11,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 6.998999999998841,
+ 21.078999999996224,
+ 34.983999999985016,
+ 48.9080000000515,
+ 62.835000000118,
+ 76.7630000001845
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.277999999998686,
+ 21.29999999999571,
+ 35.17999999998595,
+ 49.09800000005241,
+ 63.0240000001189,
+ 76.9510000001854
+ ],
+ "thresh": 10.0
+ },
+ "dend": {
+ "spikes": [
+ 7.181999999998739,
+ 21.2609999999958,
+ 35.15899999998585,
+ 49.08100000005233,
+ 63.00700000011882,
+ 76.93400000018532
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 51,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 7.003999999998838,
+ 21.10199999999617,
+ 35.03299999998525,
+ 48.98500000005187,
+ 62.9400000001185,
+ 76.89500000018514
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.284999999998682,
+ 21.32599999999565,
+ 35.232999999986205,
+ 49.17800000005279,
+ 63.13200000011942,
+ 77.08700000018605
+ ],
+ "thresh": 10.0
+ },
+ "dend": {
+ "spikes": [
+ 7.185999999998737,
+ 21.28199999999575,
+ 35.20499999998607,
+ 49.15500000005268,
+ 63.10900000011931,
+ 77.06400000018594
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 101,
+ "dt": 0.001,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 7.003999999998838,
+ 21.102999999996168,
+ 35.03499999998526,
+ 48.98800000005188,
+ 62.94400000011852,
+ 76.90000000018516
+ ],
+ "thresh": 0.0
+ },
+ "clamp": {
+ "spikes": [
+ 7.284999999998682,
+ 21.326999999995646,
+ 35.234999999986215,
+ 49.181000000052805,
+ 63.13500000011943,
+ 77.09100000018607
+ ],
+ "thresh": 10.0
+ },
+ "dend": {
+ "spikes": [
+ 7.185999999998737,
+ 21.28299999999575,
+ 35.20699999998608,
+ 49.15700000005269,
+ 63.11300000011933,
+ 77.06900000018597
+ ],
+ "thresh": -10.0
+ }
+ }
+ }
+]
\ No newline at end of file
diff --git a/data/validation/simple_exp2_synapse.json b/data/validation/simple_exp2_synapse.json
new file mode 100644
index 0000000000000000000000000000000000000000..137d9e6f86916e945b9f8c87d02369912fa03df0
--- /dev/null
+++ b/data/validation/simple_exp2_synapse.json
@@ -0,0 +1,90 @@
+[
+ {
+ "nseg": 5,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.23749999999963,
+ 21.72500000000066,
+ 41.2625000000051
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.487500000000605,
+ 41.11250000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 11,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.23749999999963,
+ 21.73750000000066,
+ 41.2750000000051
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.487500000000605,
+ 41.11250000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 51,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.23749999999963,
+ 21.73750000000066,
+ 41.2750000000051
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.487500000000605,
+ 41.11250000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 101,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.23749999999963,
+ 21.73750000000066,
+ 41.2750000000051
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.487500000000605,
+ 41.11250000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ }
+]
\ No newline at end of file
diff --git a/data/validation/simple_exp_synapse.json b/data/validation/simple_exp_synapse.json
new file mode 100644
index 0000000000000000000000000000000000000000..675f8accd5753b720efe80020f3665df9ba47318
--- /dev/null
+++ b/data/validation/simple_exp_synapse.json
@@ -0,0 +1,90 @@
+[
+ {
+ "nseg": 5,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.087499999999638,
+ 21.525000000000613,
+ 41.11250000000506
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.51250000000061,
+ 41.10000000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 11,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.099999999999637,
+ 21.537500000000616,
+ 41.125000000005066
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.525000000000613,
+ 41.11250000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 51,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.099999999999637,
+ 21.537500000000616,
+ 41.125000000005066
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.525000000000613,
+ 41.11250000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ },
+ {
+ "nseg": 101,
+ "dt": 0.0125,
+ "measurements": {
+ "soma": {
+ "spikes": [
+ 11.099999999999637,
+ 21.537500000000616,
+ 41.125000000005066
+ ],
+ "thresh": 20.0
+ },
+ "dend": {
+ "spikes": [
+ 11.087499999999638,
+ 21.525000000000613,
+ 41.11250000000506
+ ],
+ "thresh": -10.0
+ }
+ }
+ }
+]
\ No newline at end of file
diff --git a/data/validation/soma.json b/data/validation/soma.json
new file mode 100644
index 0000000000000000000000000000000000000000..b0de340ce5891e8cc95fdfb5cc92d956751f37da
--- /dev/null
+++ b/data/validation/soma.json
@@ -0,0 +1,62 @@
+[
+ {
+ "dt": 0.05,
+ "spikes": [
+ 12.04999999999985,
+ 27.57499999999897,
+ 42.85000000000118,
+ 58.10000000000465,
+ 73.37500000000811,
+ 88.62500000001158,
+ 103.90000000001505
+ ]
+ },
+ {
+ "dt": 0.02,
+ "spikes": [
+ 12.04999999999985,
+ 27.57499999999897,
+ 42.85000000000118,
+ 58.10000000000465,
+ 73.37500000000811,
+ 88.62500000001158,
+ 103.90000000001505
+ ]
+ },
+ {
+ "dt": 0.01,
+ "spikes": [
+ 12.037499999999584,
+ 27.525000000001977,
+ 42.76250000000544,
+ 57.9875000000089,
+ 73.21250000000188,
+ 88.43749999998803,
+ 103.66249999997419
+ ]
+ },
+ {
+ "dt": 0.001,
+ "spikes": [
+ 12.026000000003206,
+ 27.476999999981313,
+ 42.68400000002178,
+ 57.881000000094346,
+ 73.0770000001669,
+ 88.27300000023946,
+ 103.46900000031202
+ ]
+ },
+ {
+ "dt": 0.0001,
+ "spikes": [
+ 12.024499999978646,
+ 27.473100000350247,
+ 42.67780000085499,
+ 57.87230000135939,
+ 73.06600000186377,
+ 88.25970000236815,
+ 103.45330000287252
+ ]
+ }
+]
\ No newline at end of file
diff --git a/external/modparser b/external/modparser
index 588ca1a5ea28ef04d17b318e754d479e5489eb9a..b200bf6376a2dc30edea98fcc2375fc9be095135 160000
--- a/external/modparser
+++ b/external/modparser
@@ -1 +1 @@
-Subproject commit 588ca1a5ea28ef04d17b318e754d479e5489eb9a
+Subproject commit b200bf6376a2dc30edea98fcc2375fc9be095135
diff --git a/external/tclap b/external/tclap
new file mode 160000
index 0000000000000000000000000000000000000000..f41dcb5ce3d063c9fe95623193bba693338f3edb
--- /dev/null
+++ b/external/tclap
@@ -0,0 +1 @@
+Subproject commit f41dcb5ce3d063c9fe95623193bba693338f3edb
diff --git a/mechanisms/CMakeLists.txt b/mechanisms/CMakeLists.txt
index 203dedff949356b1a715146baf0ad6d6714aac56..cf183d75545c011dbb5fbb3436221ab7ab4d0933 100644
--- a/mechanisms/CMakeLists.txt
+++ b/mechanisms/CMakeLists.txt
@@ -1,4 +1,4 @@
-set(mechanisms pas hh expsyn)
+set(mechanisms pas hh expsyn exp2syn)
set(modcc "${CMAKE_BINARY_DIR}/external/bin/modcc")
diff --git a/mechanisms/mod/exp2syn.mod b/mechanisms/mod/exp2syn.mod
new file mode 100755
index 0000000000000000000000000000000000000000..ea0cf025e9228216b4699161a2e574ca1d1dda6e
--- /dev/null
+++ b/mechanisms/mod/exp2syn.mod
@@ -0,0 +1,50 @@
+NEURON {
+ POINT_PROCESS exp2syn
+ RANGE tau1, tau2, e
+ NONSPECIFIC_CURRENT i
+}
+
+UNITS {
+ (nA) = (nanoamp)
+ (mV) = (millivolt)
+ (uS) = (microsiemens)
+}
+
+PARAMETER {
+ tau1 = .5 (ms) : <1e-9,1e9>
+ tau2 = 2 (ms) : <1e-9,1e9>
+ e=0 (mV)
+}
+
+ASSIGNED {
+ factor
+}
+
+STATE {
+ A : (uS)
+ B : (uS)
+}
+
+INITIAL {
+ LOCAL tp
+ A = 0
+ B = 0
+ tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1)
+ factor = -exp(-tp/tau1) + exp(-tp/tau2)
+ factor = 1/factor
+}
+
+BREAKPOINT {
+ SOLVE state METHOD cnexp
+ i = (B - A)*(v - e)
+}
+
+DERIVATIVE state {
+ A' = -A/tau1
+ B' = -B/tau2
+}
+
+NET_RECEIVE(weight) {
+ A = A + weight*factor
+ B = B + weight*factor
+}
diff --git a/mechanisms/mod/expsyn.mod b/mechanisms/mod/expsyn.mod
index db1f6d6ed61a8b563141ea2516658b8f87a35637..16d10c5cf856cd87e50252dca2163d844588e2b9 100644
--- a/mechanisms/mod/expsyn.mod
+++ b/mechanisms/mod/expsyn.mod
@@ -1,5 +1,5 @@
NEURON {
- POINT_PROCESS ExpSyn
+ POINT_PROCESS expsyn
RANGE tau, e
NONSPECIFIC_CURRENT i
}
diff --git a/miniapp/CMakeLists.txt b/miniapp/CMakeLists.txt
index a5b9561198875a6515b4e8a69f7685ea75c19803..682adefe1d04b057cb330dc9468485b276b797b2 100644
--- a/miniapp/CMakeLists.txt
+++ b/miniapp/CMakeLists.txt
@@ -1,9 +1,9 @@
set(HEADERS
)
set(MINIAPP_SOURCES
- # mpi.cpp
io.cpp
miniapp.cpp
+ miniapp_recipes.cpp
)
add_executable(miniapp.exe ${MINIAPP_SOURCES} ${HEADERS})
diff --git a/miniapp/README.md b/miniapp/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..9f8bd3822fdd8cca08442412100e5676564bdf17
--- /dev/null
+++ b/miniapp/README.md
@@ -0,0 +1,90 @@
+# NEST MC miniapp
+
+The miniapp is a simple benchmark of the NEST MC framework.
+
+## the model
+
+The model is very simple, with three parameters that describe the network, and two for controlling time stepping.
+
+### network/model parameters
+
+The following parameters are used to describe the size, connectivity and resolution of the model:
+
+- `cells` is the total number of cells in the network.
+- `synapses_per_cell` the number of synapses per cell, must be in the range `[0,cells-1]`
+- `compartments` the number of compartments per segment.
+- `all_to_all` toggle whether to make an all to all network
+
+All cells have identical morphology, a soma with a dendrite attached. The dendrite branches as illustrated (roughly) below
+
+
+```
+s = soma
+. = branching point of dendrite
+
+ /
+ /
+s------.
+ \
+ \
+```
+
+The disrcetization of each cell is controlled with the __compartments__ parameter.
+For example, when `compartments=100`, the total number of compartments in each cell is 301: 1 for the soma, and 100 for each of the dendrite segments.
+
+The `synapses_per_cell` parameter is in the range `[0,cells-1]`.
+If it is zero, then there are no connections between the cells (not much of a network).
+By default, the source gid of each synapse is chosen at random from the global set of cells (excluding the cell of the synapse).
+If the `all_to_all` flag is set, `synapses_per_cell` is set to `cells-1`, i.e. one connection for each cell (excluding the cell of the synapse)
+
+Note that the to avoid numerical instability, the number of synapses per cell should be greater than 200 (or zero!).
+The number of synapses per cell required for stability is dependent on the number of compartments per segment (fewer compartments is more stable) and the time step size (smaller time step sizes increase stability).
+If there are numeric instabilities the simulation will print a warning
+```
+warning: solution out of bounds
+```
+
+### time stepping parameters
+
+The time stepping can be controlled via two parameters
+
+- `dt` the time step size in ms (default `dt=0.025`)
+- `tfinal` the length of the simulation in ms (default `tfinal=200`)
+
+## configuration
+
+There are two ways to specify the model properties, i.e. the number of cells, connections per cell, etc.
+
+### command line arguments
+
+- `-n integer` : `ncells`
+- `-s integer` : `synapses_per_cell`
+- `-c integer` : `compartments`
+- `-d float` : `dt`
+- `-t float` : `tfinal`
+- `-i filename` : name of json file with parameters
+
+For example
+
+```
+> ./miniapp.exe -n 1000 -s 500 -c 50 -t 100 -d 0.02
+```
+
+will run a simulation with 1000 cells, with 500 synapses per cell, 50 compartments per segment for a total of 100 ms with 0.02 ms time steps.
+
+### input parameter file
+
+To run the same simulation that we ran with the command line arguments with an input file:
+
+```
+> cat input.json
+{
+ "cells": 1000,
+ "synapses": 500,
+ "compartments": 50,
+ "dt": 0.02,
+ "tfinal": 100.0,
+ "all_to_all": false
+}
+> ./miniapp.exe -i input.json
+```
diff --git a/miniapp/io.cpp b/miniapp/io.cpp
index 887848b549b7bad46574b92ac9fcde74faae5d5c..6664cf47e4ddbf1fcd33cd4589553a29cb967255 100644
--- a/miniapp/io.cpp
+++ b/miniapp/io.cpp
@@ -1,21 +1,106 @@
+#include <fstream>
+#include <exception>
+
+#include <tclap/CmdLine.h>
+#include <json/src/json.hpp>
+
#include "io.hpp"
namespace nest {
namespace mc {
namespace io {
-// read simulation options from json file with name fname
-// for now this is just a placeholder
-options read_options(std::string fname) {
- // 10 cells, 1 synapses per cell, 10 compartments per segment
- return {200, 1, 100};
+/// read simulation options from json file with name fname
+/// if file name is empty or if file is not a valid json file a default
+/// set of parameters is returned :
+/// 1000 cells, 500 synapses per cell, 100 compartments per segment
+cl_options read_options(int argc, char** argv) {
+
+ // set default options
+ const cl_options defopts{"", 1000, 500, "expsyn", 100, 100., 0.025, false};
+
+ cl_options options;
+ // parse command line arguments
+ try {
+ TCLAP::CmdLine cmd("mod2c performance benchmark harness", ' ', "0.1");
+
+ TCLAP::ValueArg<uint32_t> ncells_arg(
+ "n", "ncells", "total number of cells in the model",
+ false, defopts.cells, "non negative integer", cmd);
+ TCLAP::ValueArg<uint32_t> nsynapses_arg(
+ "s", "nsynapses", "number of synapses per cell",
+ false, defopts.synapses_per_cell, "non negative integer", cmd);
+ TCLAP::ValueArg<std::string> syntype_arg(
+ "S", "syntype", "type of synapse (expsyn or exp2syn)",
+ false, defopts.syn_type, "synapse type", cmd);
+ TCLAP::ValueArg<uint32_t> ncompartments_arg(
+ "c", "ncompartments", "number of compartments per segment",
+ false, defopts.compartments_per_segment, "non negative integer", cmd);
+ TCLAP::ValueArg<std::string> ifile_arg(
+ "i", "ifile", "json file with model parameters",
+ false, "","file name string", cmd);
+ TCLAP::ValueArg<double> tfinal_arg(
+ "t", "tfinal", "time to simulate in ms",
+ false, defopts.tfinal, "positive real number", cmd);
+ TCLAP::ValueArg<double> dt_arg(
+ "d", "dt", "time step size in ms",
+ false, defopts.dt, "positive real number", cmd);
+ TCLAP::SwitchArg all_to_all_arg(
+ "m","alltoall","all to all network", cmd, false);
+
+ cmd.parse(argc, argv);
+
+ options.cells = ncells_arg.getValue();
+ options.synapses_per_cell = nsynapses_arg.getValue();
+ options.syn_type = syntype_arg.getValue();
+ options.compartments_per_segment = ncompartments_arg.getValue();
+ options.ifname = ifile_arg.getValue();
+ options.tfinal = tfinal_arg.getValue();
+ options.dt = dt_arg.getValue();
+ options.all_to_all = all_to_all_arg.getValue();
+ }
+ // catch any exceptions in command line handling
+ catch (TCLAP::ArgException& e) {
+ throw usage_error("error parsing command line argument "+e.argId()+": "+e.error());
+ }
+
+ if (options.ifname != "") {
+ std::ifstream fid(options.ifname);
+ if (fid) {
+ // read json data in input file
+ nlohmann::json fopts;
+ fid >> fopts;
+
+ try {
+ options.cells = fopts["cells"];
+ options.synapses_per_cell = fopts["synapses"];
+ options.compartments_per_segment = fopts["compartments"];
+ options.dt = fopts["dt"];
+ options.tfinal = fopts["tfinal"];
+ options.all_to_all = fopts["all_to_all"];
+ }
+ catch (std::exception& e) {
+ throw model_description_error(
+ "unable to parse parameters in "+options.ifname+": "+e.what());
+ }
+ }
+ else {
+ throw usage_error("unable to open model parameter file "+options.ifname);
+ }
+ }
+
+ return options;
}
-std::ostream& operator<<(std::ostream& o, const options& opt) {
+std::ostream& operator<<(std::ostream& o, const cl_options& options) {
o << "simultion options:\n";
- o << " cells : " << opt.cells << "\n";
- o << " compartments/segment : " << opt.compartments_per_segment << "\n";
- o << " synapses/cell : " << opt.synapses_per_cell << "\n";
+ o << " cells : " << options.cells << "\n";
+ o << " compartments/segment : " << options.compartments_per_segment << "\n";
+ o << " synapses/cell : " << options.synapses_per_cell << "\n";
+ o << " simulation time : " << options.tfinal << "\n";
+ o << " dt : " << options.dt << "\n";
+ o << " all to all network : " << (options.all_to_all ? "yes" : "no") << "\n";
+ o << " input file name : " << options.ifname << "\n";
return o;
}
diff --git a/miniapp/io.hpp b/miniapp/io.hpp
index e8a8762ceb34b7d7148dfb7eb831217ac633d3d3..dab584fe2a132f289bb9a2f1ae02711aa29ef5a4 100644
--- a/miniapp/io.hpp
+++ b/miniapp/io.hpp
@@ -1,21 +1,43 @@
#pragma once
-#include <json/src/json.hpp>
+#include <string>
+#include <cstdint>
+#include <iosfwd>
+#include <stdexcept>
+#include <utility>
namespace nest {
namespace mc {
namespace io {
// holds the options for a simulation run
-struct options {
- int cells;
- int synapses_per_cell;
- int compartments_per_segment;
+struct cl_options {
+ std::string ifname;
+ uint32_t cells;
+ uint32_t synapses_per_cell;
+ std::string syn_type;
+ uint32_t compartments_per_segment;
+ double tfinal;
+ double dt;
+ bool all_to_all;
};
-std::ostream& operator<<(std::ostream& o, const options& opt);
+class usage_error: public std::runtime_error {
+public:
+ template <typename S>
+ usage_error(S&& whatmsg): std::runtime_error(std::forward<S>(whatmsg)) {}
+};
+
+class model_description_error: public std::runtime_error {
+public:
+ template <typename S>
+ model_description_error(S&& whatmsg): std::runtime_error(std::forward<S>(whatmsg)) {}
+};
+
+std::ostream& operator<<(std::ostream& o, const cl_options& opt);
+
+cl_options read_options(int argc, char** argv);
-options read_options(std::string fname);
} // namespace io
} // namespace mc
diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index 74b551ff4f07fe795514caafde901de1721108dc..9b26ff86ef37f829667161afca27a129791de3ed 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -1,425 +1,178 @@
+#include <cmath>
+#include <exception>
#include <iostream>
#include <fstream>
-#include <sstream>
+#include <memory>
+#include <json/src/json.hpp>
+
+#include <common_types.hpp>
#include <cell.hpp>
#include <cell_group.hpp>
#include <fvm_cell.hpp>
-#include <mechanism_interface.hpp>
+#include <mechanism_catalogue.hpp>
+#include <model.hpp>
+#include <threading/threading.hpp>
+#include <profiling/profiler.hpp>
+#include <communication/communicator.hpp>
+#include <communication/global_policy.hpp>
+#include <util/ioutil.hpp>
+#include <util/optional.hpp>
#include "io.hpp"
-#include "threading/threading.hpp"
-#include "profiling/profiler.hpp"
-#include "communication/communicator.hpp"
-#include "communication/global_policy.hpp"
-#include "util/optional.hpp"
+#include "miniapp_recipes.hpp"
+#include "trace_sampler.hpp"
-using namespace nest;
+using namespace nest::mc;
-using real_type = double;
-using index_type = int;
-using id_type = uint32_t;
-using numeric_cell = mc::fvm::fvm_cell<real_type, index_type>;
-using cell_group = mc::cell_group<numeric_cell>;
+using global_policy = communication::global_policy;
-using global_policy = nest::mc::communication::global_policy;
-using communicator_type =
- mc::communication::communicator<global_policy>;
+using lowered_cell = fvm::fvm_cell<double, cell_local_size_type>;
+using model_type = model<lowered_cell>;
+using sample_trace_type = sample_trace<model_type::time_type, model_type::value_type>;
-using nest::mc::util::optional;
+void banner();
+std::unique_ptr<recipe> make_recipe(const io::cl_options&);
+std::unique_ptr<sample_trace_type> make_trace(cell_member_type probe_id, probe_spec probe);
+std::pair<cell_gid_type, cell_gid_type> distribute_cells(cell_size_type ncells);
-struct model {
- communicator_type communicator;
- std::vector<cell_group> cell_groups;
+void write_trace_json(const sample_trace_type& trace, const std::string& prefix = "trace_");
- int num_groups() const {
- return cell_groups.size();
- }
+int main(int argc, char** argv) {
+ nest::mc::communication::global_policy_guard global_guard(argc, argv);
- void run(double tfinal, double dt) {
- auto t = 0.;
- auto delta = communicator.min_delay();
- while(t<tfinal) {
- mc::threading::parallel_for::apply(
- 0, num_groups(),
- [&](int i) {
- mc::util::profiler_enter("stepping","events");
- cell_groups[i].enqueue_events(communicator.queue(i));
- mc::util::profiler_leave();
- cell_groups[i].advance(t+delta, dt);
- mc::util::profiler_enter("events");
- communicator.add_spikes(cell_groups[i].spikes());
- cell_groups[i].clear_spikes();
- mc::util::profiler_leave(2);
- }
- );
-
- mc::util::profiler_enter("stepping", "exchange");
- communicator.exchange();
- mc::util::profiler_leave(2);
-
- t += delta;
+ try {
+ std::cout << util::mask_stream(global_policy::id()==0);
+ banner();
+
+ // read parameters
+ io::cl_options options = io::read_options(argc, argv);
+ std::cout << options << "\n";
+ std::cout << "\n";
+ std::cout << ":: simulation to " << options.tfinal << " ms in "
+ << std::ceil(options.tfinal / options.dt) << " steps of "
+ << options.dt << " ms" << std::endl;
+
+ auto recipe = make_recipe(options);
+ auto cell_range = distribute_cells(recipe->num_cells());
+
+ // build model from recipe
+ model_type m(*recipe, cell_range.first, cell_range.second);
+
+ // inject some artificial spikes, 1 per 20 neurons.
+ cell_gid_type spike_cell = 20*((cell_range.first+19)/20);
+ for (; spike_cell<cell_range.second; spike_cell+=20) {
+ m.add_artificial_spike({spike_cell,0u});
}
- }
-
- void init_communicator() {
- mc::util::profiler_enter("setup", "communicator");
- // calculate the source and synapse distribution serially
- std::vector<id_type> target_counts(num_groups());
- std::vector<id_type> source_counts(num_groups());
- for (auto i=0; i<num_groups(); ++i) {
- target_counts[i] = cell_groups[i].cell().synapses()->size();
- source_counts[i] = cell_groups[i].spike_sources().size();
+ // attach samplers to all probes
+ std::vector<std::unique_ptr<sample_trace_type>> traces;
+ const model_type::time_type sample_dt = 0.1;
+ for (auto probe: m.probes()) {
+ traces.push_back(make_trace(probe.id, probe.probe));
+ m.attach_sampler(probe.id, make_trace_sampler(traces.back().get(),sample_dt));
}
- target_map = mc::algorithms::make_index(target_counts);
- source_map = mc::algorithms::make_index(source_counts);
+ // run model
+ m.run(options.tfinal, options.dt);
+ util::profiler_output(0.001);
- // create connections
- communicator = communicator_type(num_groups(), target_counts);
-
- mc::util::profiler_leave(2);
- }
+ std::cout << "there were " << m.num_spikes() << " spikes\n";
- void update_gids() {
- mc::util::profiler_enter("setup", "globalize");
- auto com_policy = communicator.communication_policy();
- auto global_source_map = com_policy.make_map(source_map.back());
- auto domain_idx = communicator.domain_id();
- for (auto i=0; i<num_groups(); ++i) {
- cell_groups[i].set_source_gids(source_map[i]+global_source_map[domain_idx]);
- cell_groups[i].set_target_gids(target_map[i]+communicator.target_gid_from_group_lid(0));
- }
- mc::util::profiler_leave(2);
- }
-
- // TODO : only stored here because init_communicator() and update_gids() are split
- std::vector<id_type> source_map;
- std::vector<id_type> target_map;
-
- // traces from probes
- struct trace_data {
- struct sample_type {
- float time;
- double value;
- };
- std::string name;
- index_type id;
- std::vector<sample_type> samples;
- };
-
- // different traces may be written to by different threads;
- // during simulation, each trace_sampler will be responsible for its
- // corresponding element in the traces vector.
-
- std::vector<trace_data> traces;
-
- // make a sampler that records to traces
- struct simple_sampler_functor {
- std::vector<trace_data> &traces_;
- size_t trace_index_ = 0;
- float requested_sample_time_ = 0;
- float dt_ = 0;
-
- simple_sampler_functor(std::vector<trace_data> &traces, size_t index, float dt) :
- traces_(traces), trace_index_(index), dt_(dt)
- {}
-
- optional<float> operator()(float t, double v) {
- traces_[trace_index_].samples.push_back({t,v});
- return requested_sample_time_ += dt_;
- }
- };
-
- mc::sampler make_simple_sampler(
- index_type probe_gid, const std::string name, index_type id, float dt)
- {
- traces.push_back(trace_data{name, id});
- return {probe_gid, simple_sampler_functor(traces, traces.size()-1, dt)};
- }
-
- void reset_traces() {
- // do not call during simulation: thread-unsafe access to traces.
- traces.clear();
- }
-
- void dump_traces() {
- // do not call during simulation: thread-unsafe access to traces.
+ // save traces
for (const auto& trace: traces) {
- auto path = "trace_" + std::to_string(trace.id)
- + "_" + trace.name + ".json";
-
- nlohmann::json json;
- json["name"] = trace.name;
- for (const auto& sample: trace.samples) {
- json["time"].push_back(sample.time);
- json["value"].push_back(sample.value);
- }
- std::ofstream file(path);
- file << std::setw(1) << json << std::endl;
+ write_trace_json(*trace.get());
}
}
-};
-
-// define some global model parameters
-namespace parameters {
-namespace synapses {
- // synapse delay
- constexpr double delay = 5.0; // ms
-
- // connection weight
- constexpr double weight = 0.005; // uS
-}
-}
-
-///////////////////////////////////////
-// prototypes
-///////////////////////////////////////
-
-/// make a single abstract cell
-mc::cell make_cell(int compartments_per_segment, int num_synapses);
-
-/// do basic setup (initialize global state, print banner, etc)
-void setup();
-
-/// helper function for initializing cells
-cell_group make_lowered_cell(int cell_index, const mc::cell& c);
-
-/// models
-void ring_model(nest::mc::io::options& opt, model& m);
-void all_to_all_model(nest::mc::io::options& opt, model& m);
-
-
-///////////////////////////////////////
-// main
-///////////////////////////////////////
-int main(int argc, char** argv) {
- nest::mc::communication::global_policy_guard global_guard(argc, argv);
-
- setup();
-
- // read parameters
- mc::io::options opt;
- try {
- opt = mc::io::read_options("");
- if (!global_policy::id()) {
- std::cout << opt << "\n";
- }
+ catch (io::usage_error& e) {
+ // only print usage/startup errors on master
+ std::cerr << util::mask_stream(global_policy::id()==0);
+ std::cerr << e.what() << "\n";
+ return 1;
}
catch (std::exception& e) {
- std::cerr << e.what() << std::endl;
- exit(1);
+ std::cerr << e.what() << "\n";
+ return 2;
}
+ return 0;
+}
- model m;
- all_to_all_model(opt, m);
-
- //
- // time stepping
- //
- auto tfinal = 20.;
- auto dt = 0.01;
-
- auto id = m.communicator.domain_id();
-
- if (!id) {
- m.communicator.add_spike({0, 5});
- }
+std::pair<cell_gid_type, cell_gid_type> distribute_cells(cell_size_type num_cells) {
+ // Crude load balancing:
+ // divide [0, num_cells) into num_domains non-overlapping, contiguous blocks
+ // of size as close to equal as possible.
- m.run(tfinal, dt);
+ auto num_domains = communication::global_policy::size();
+ auto domain_id = communication::global_policy::id();
- mc::util::profiler_output(0.001);
+ cell_gid_type cell_from = (cell_gid_type)(num_cells*(domain_id/(double)num_domains));
+ cell_gid_type cell_to = (cell_gid_type)(num_cells*((domain_id+1)/(double)num_domains));
- if (!id) {
- std::cout << "there were " << m.communicator.num_spikes() << " spikes\n";
- }
- m.dump_traces();
-
-#ifdef SPLAT
- if (!global_policy::id()) {
- //for (auto i=0u; i<m.cell_groups.size(); ++i) {
- m.cell_groups[0].splat("cell0.txt");
- m.cell_groups[1].splat("cell1.txt");
- m.cell_groups[2].splat("cell2.txt");
- //}
- }
-#endif
+ return {cell_from, cell_to};
}
-///////////////////////////////////////
-// models
-///////////////////////////////////////
-
-void ring_model(nest::mc::io::options& opt, model& m) {
- //
- // make cells
- //
-
- // make a basic cell
- auto basic_cell = make_cell(opt.compartments_per_segment, 1);
-
- // make a vector for storing all of the cells
- m.cell_groups = std::vector<cell_group>(opt.cells);
-
- // initialize the cells in parallel
- mc::threading::parallel_for::apply(
- 0, opt.cells,
- [&](int i) {
- // initialize cell
- mc::util::profiler_enter("setup");
- mc::util::profiler_enter("make cell");
- m.cell_groups[i] = make_lowered_cell(i, basic_cell);
- mc::util::profiler_leave();
- mc::util::profiler_leave();
- }
- );
-
- //
- // network creation
- //
- m.init_communicator();
-
- for (auto i=0u; i<(id_type)opt.cells; ++i) {
- m.communicator.add_connection({
- i, (i+1)%opt.cells,
- parameters::synapses::weight, parameters::synapses::delay
- });
- }
-
- m.update_gids();
+void banner() {
+ std::cout << "====================\n";
+ std::cout << " starting miniapp\n";
+ std::cout << " - " << threading::description() << " threading support\n";
+ std::cout << " - communication policy: " << global_policy::name() << "\n";
+ std::cout << "====================\n";
}
-void all_to_all_model(nest::mc::io::options& opt, model& m) {
- //
- // make cells
- //
+std::unique_ptr<recipe> make_recipe(const io::cl_options& options) {
+ basic_recipe_param p;
- // make a basic cell
- auto basic_cell = make_cell(opt.compartments_per_segment, opt.cells-1);
+ p.num_compartments = options.compartments_per_segment;
+ p.num_synapses = options.all_to_all? options.cells-1: options.synapses_per_cell;
+ p.synapse_type = options.syn_type;
- // make a vector for storing all of the cells
- id_type ncell_global = opt.cells;
- id_type ncell_local = ncell_global / m.communicator.num_domains();
- int remainder = ncell_global - (ncell_local*m.communicator.num_domains());
- if (m.communicator.domain_id()<remainder) {
- ncell_local++;
+ if (options.all_to_all) {
+ return make_basic_kgraph_recipe(options.cells, p);
}
-
- m.cell_groups = std::vector<cell_group>(ncell_local);
-
- // initialize the cells in parallel
- mc::threading::parallel_for::apply(
- 0, ncell_local,
- [&](int i) {
- mc::util::profiler_enter("setup", "cells");
- m.cell_groups[i] = make_lowered_cell(i, basic_cell);
- mc::util::profiler_leave(2);
- }
- );
-
- //
- // network creation
- //
- m.init_communicator();
-
- // monitor soma and dendrite on a few cells
- float sample_dt = 0.1;
- index_type monitor_group_gids[] = { 0, 1, 2 };
- for (auto gid : monitor_group_gids) {
- if (!m.communicator.is_local_group(gid)) {
- continue;
- }
-
- auto lid = m.communicator.group_lid(gid);
- auto probe_soma = m.cell_groups[lid].probe_gid_range().first;
- auto probe_dend = probe_soma+1;
-
- m.cell_groups[lid].add_sampler(m.make_simple_sampler(probe_soma, "vsoma", gid, sample_dt));
- m.cell_groups[lid].add_sampler(m.make_simple_sampler(probe_dend, "vdend", gid, sample_dt));
- }
-
- mc::util::profiler_enter("setup", "connections");
- // lid is local cell/group id
- for (auto lid=0u; lid<ncell_local; ++lid) {
- auto target = m.communicator.target_gid_from_group_lid(lid);
- auto gid = m.communicator.group_gid_from_group_lid(lid);
- // tid is global cell/group id
- for (auto tid=0u; tid<ncell_global; ++tid) {
- if (gid!=tid) {
- m.communicator.add_connection({
- tid, target++,
- parameters::synapses::weight, parameters::synapses::delay
- });
- }
- }
+ else {
+ return make_basic_rgraph_recipe(options.cells, p);
}
-
- m.communicator.construct();
- mc::util::profiler_leave(2);
-
- m.update_gids();
}
-///////////////////////////////////////
-// function definitions
-///////////////////////////////////////
-
-void setup() {
- // print banner
- if (!global_policy::id()) {
- std::cout << "====================\n";
- std::cout << " starting miniapp\n";
- std::cout << " - " << mc::threading::description() << " threading support\n";
- std::cout << " - communication policy: " << global_policy::name() << "\n";
- std::cout << "====================\n";
- }
-
- // setup global state for the mechanisms
- mc::mechanisms::setup_mechanism_helpers();
+std::unique_ptr<sample_trace_type> make_trace(cell_member_type probe_id, probe_spec probe) {
+ std::string name = "";
+ std::string units = "";
+
+ switch (probe.kind) {
+ case probeKind::membrane_voltage:
+ name = "v";
+ units = "mV";
+ break;
+ case probeKind::membrane_current:
+ name = "i";
+ units = "mA/cm²";
+ break;
+ default: ;
+ }
+ name += probe.location.segment? "dend" : "soma";
+
+ return util::make_unique<sample_trace_type>(probe_id, name, units);
}
-// make a high level cell description for use in simulation
-mc::cell make_cell(int compartments_per_segment, int num_synapses) {
- nest::mc::cell cell;
+void write_trace_json(const sample_trace_type& trace, const std::string& prefix) {
+ auto path = prefix + std::to_string(trace.probe_id.gid) +
+ "." + std::to_string(trace.probe_id.index) + "_" + trace.name + ".json";
- // Soma with diameter 12.6157 um and HH channel
- auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism(mc::hh_parameters());
+ nlohmann::json jrep;
+ jrep["name"] = trace.name;
+ jrep["units"] = trace.units;
+ jrep["cell"] = trace.probe_id.gid;
+ jrep["probe"] = trace.probe_id.index;
- // add dendrite of length 200 um and diameter 1 um with passive channel
- std::vector<mc::cable_segment*> dendrites;
- dendrites.push_back(cell.add_cable(0, mc::segmentKind::dendrite, 0.5, 0.5, 200));
- //dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25, 100));
- //dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25, 100));
+ auto& jt = jrep["data"]["time"];
+ auto& jy = jrep["data"][trace.name];
- for (auto d : dendrites) {
- d->add_mechanism(mc::pas_parameters());
- d->set_compartments(compartments_per_segment);
- d->mechanism("membrane").set("r_L", 100);
+ for (const auto& sample: trace.samples) {
+ jt.push_back(sample.time);
+ jy.push_back(sample.value);
}
-
- // add stimulus
- //cell.add_stimulus({1,1}, {5., 80., 0.3});
-
- cell.add_detector({0,0}, 30);
-
- for (auto i=0; i<num_synapses; ++i) {
- cell.add_synapse({1, 0.5});
- }
-
- // add probes:
- auto probe_soma = cell.add_probe({0, 0}, mc::probeKind::membrane_voltage);
- auto probe_dendrite = cell.add_probe({1, 0.5}, mc::probeKind::membrane_voltage);
-
- EXPECTS(probe_soma==0);
- EXPECTS(probe_dendrite==1);
- (void)probe_soma, (void)probe_dendrite;
-
- return cell;
+ std::ofstream file(path);
+ file << std::setw(1) << jrep << std::endl;
}
-cell_group make_lowered_cell(int cell_index, const mc::cell& c) {
- return cell_group(c);
-}
diff --git a/miniapp/miniapp_recipes.cpp b/miniapp/miniapp_recipes.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..22c59bd851a049e448784b9cfdd635d7e7a96970
--- /dev/null
+++ b/miniapp/miniapp_recipes.cpp
@@ -0,0 +1,234 @@
+#include <cmath>
+#include <random>
+#include <vector>
+#include <utility>
+
+#include <cell.hpp>
+#include <util/debug.hpp>
+
+#include "miniapp_recipes.hpp"
+
+namespace nest {
+namespace mc {
+
+// TODO: split cell description into separate morphology, stimulus, probes, mechanisms etc.
+// description for greater data reuse.
+
+template <typename RNG>
+cell make_basic_cell(
+ unsigned compartments_per_segment,
+ unsigned num_synapses,
+ const std::string& syn_type,
+ RNG& rng)
+{
+ nest::mc::cell cell;
+
+ // Soma with diameter 12.6157 um and HH channel
+ auto soma = cell.add_soma(12.6157/2.0);
+ soma->add_mechanism(mc::hh_parameters());
+
+ // add dendrite of length 200 um and diameter 1 um with passive channel
+ std::vector<mc::cable_segment*> dendrites;
+ dendrites.push_back(cell.add_cable(0, mc::segmentKind::dendrite, 0.5, 0.5, 200));
+ dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
+ dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
+
+ for (auto d : dendrites) {
+ d->add_mechanism(mc::pas_parameters());
+ d->set_compartments(compartments_per_segment);
+ d->mechanism("membrane").set("r_L", 100);
+ }
+
+ cell.add_detector({0,0}, 20);
+
+ auto distribution = std::uniform_real_distribution<float>(0.f, 1.0f);
+ // distribute the synapses at random locations the terminal dendrites in a
+ // round robin manner; the terminal dendrites in this cell have indices 2 and 3.
+ nest::mc::parameter_list syn_default(syn_type);
+ for (unsigned i=0; i<num_synapses; ++i) {
+ cell.add_synapse({2+(i%2), distribution(rng)}, syn_default);
+ }
+
+ return cell;
+}
+
+class basic_cell_recipe: public recipe {
+public:
+ basic_cell_recipe(cell_gid_type ncell, basic_recipe_param param, probe_distribution pdist):
+ ncell_(ncell), param_(std::move(param)), pdist_(std::move(pdist))
+ {
+ delay_distribution_param = exp_param{param_.mean_connection_delay_ms
+ - param_.min_connection_delay_ms};
+ }
+
+ cell_size_type num_cells() const override { return ncell_; }
+
+ cell get_cell(cell_gid_type i) const override {
+ auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
+
+ auto cc = get_cell_count_info(i);
+ auto cell = make_basic_cell(param_.num_compartments, cc.num_targets,
+ param_.synapse_type, gen);
+
+ EXPECTS(cell.num_segments()==basic_cell_segments);
+ EXPECTS(cell.probes().size()==0);
+ EXPECTS(cell.synapses().size()==cc.num_targets);
+ EXPECTS(cell.detectors().size()==cc.num_sources);
+
+ // add probes
+ unsigned n_probe_segs = pdist_.all_segments? basic_cell_segments: 1u;
+ for (unsigned i = 0; i<n_probe_segs; ++i) {
+ if (pdist_.membrane_voltage) {
+ cell.add_probe({{i, i? 0.5: 0.0}, mc::probeKind::membrane_voltage});
+ }
+ if (pdist_.membrane_current) {
+ cell.add_probe({{i, i? 0.5: 0.0}, mc::probeKind::membrane_current});
+ }
+ }
+ EXPECTS(cell.probes().size()==cc.num_probes);
+ return cell;
+ }
+
+ cell_count_info get_cell_count_info(cell_gid_type i) const override {
+ cell_count_info cc = {1, param_.num_synapses, 0 };
+
+ // probe this cell?
+ if (std::floor(i*pdist_.proportion)!=std::floor((i-1.0)*pdist_.proportion)) {
+ std::size_t np = pdist_.membrane_voltage + pdist_.membrane_current;
+ if (pdist_.all_segments) {
+ np *= basic_cell_segments;
+ }
+
+ cc.num_probes = np;
+ }
+
+ return cc;
+ }
+
+protected:
+ template <typename RNG>
+ cell_connection draw_connection_params(RNG& rng) const {
+ std::exponential_distribution<float> delay_dist(delay_distribution_param);
+ float delay = param_.min_connection_delay_ms + delay_dist(rng);
+ float weight = param_.syn_weight_per_cell/param_.num_synapses;
+ return cell_connection{{0, 0}, {0, 0}, weight, delay};
+ }
+
+ cell_gid_type ncell_;
+ basic_recipe_param param_;
+ probe_distribution pdist_;
+ static constexpr int basic_cell_segments = 4;
+
+ using exp_param = std::exponential_distribution<float>::param_type;
+ exp_param delay_distribution_param;
+};
+
+class basic_ring_recipe: public basic_cell_recipe {
+public:
+ basic_ring_recipe(cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist = probe_distribution{}):
+ basic_cell_recipe(ncell, std::move(param), std::move(pdist)) {}
+
+ std::vector<cell_connection> connections_on(cell_gid_type i) const override {
+ std::vector<cell_connection> conns;
+ auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
+
+ cell_gid_type prev = i==0? ncell_-1: i-1;
+ for (unsigned t=0; t<param_.num_synapses; ++t) {
+ cell_connection cc = draw_connection_params(gen);
+ cc.source = {prev, 0};
+ cc.dest = {i, t};
+ conns.push_back(cc);
+ }
+
+ return conns;
+ }
+};
+
+std::unique_ptr<recipe> make_basic_ring_recipe(
+ cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist)
+{
+ return std::unique_ptr<recipe>(new basic_ring_recipe(ncell, param, pdist));
+}
+
+
+class basic_rgraph_recipe: public basic_cell_recipe {
+public:
+ basic_rgraph_recipe(cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist = probe_distribution{}):
+ basic_cell_recipe(ncell, std::move(param), std::move(pdist)) {}
+
+ std::vector<cell_connection> connections_on(cell_gid_type i) const override {
+ std::vector<cell_connection> conns;
+ auto conn_param_gen = std::mt19937(i); // TODO: replace this with hashing generator...
+ auto source_gen = std::mt19937(i*123+457); // ditto
+
+ std::uniform_int_distribution<cell_gid_type> source_distribution(0, ncell_-2);
+
+ for (unsigned t=0; t<param_.num_synapses; ++t) {
+ auto source = source_distribution(source_gen);
+ if (source>=i) ++source;
+
+ cell_connection cc = draw_connection_params(conn_param_gen);
+ cc.source = {source, 0};
+ cc.dest = {i, t};
+ conns.push_back(cc);
+ }
+
+ return conns;
+ }
+};
+
+std::unique_ptr<recipe> make_basic_rgraph_recipe(
+ cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist)
+{
+ return std::unique_ptr<recipe>(new basic_rgraph_recipe(ncell, param, pdist));
+}
+
+class basic_kgraph_recipe: public basic_cell_recipe {
+public:
+ basic_kgraph_recipe(cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist = probe_distribution{}):
+ basic_cell_recipe(ncell, std::move(param), std::move(pdist))
+ {
+ if (std::size_t(param.num_synapses) != ncell-1) {
+ throw invalid_recipe_error("number of synapses per cell must equal number "
+ "of cells minus one in complete graph model");
+ }
+ }
+
+ std::vector<cell_connection> connections_on(cell_gid_type i) const override {
+ std::vector<cell_connection> conns;
+ auto conn_param_gen = std::mt19937(i); // TODO: replace this with hashing generator...
+
+ for (unsigned t=0; t<param_.num_synapses; ++t) {
+ cell_gid_type source = t>=i? t+1: t;
+ EXPECTS(source<ncell_);
+
+ cell_connection cc = draw_connection_params(conn_param_gen);
+ cc.source = {source, 0};
+ cc.dest = {i, t};
+ conns.push_back(cc);
+ }
+
+ return conns;
+ }
+};
+
+std::unique_ptr<recipe> make_basic_kgraph_recipe(
+ cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist)
+{
+ return std::unique_ptr<recipe>(new basic_kgraph_recipe(ncell, param, pdist));
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/miniapp/miniapp_recipes.hpp b/miniapp/miniapp_recipes.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2e519bd1b492e532db86175acbd11520d1bd7a8b
--- /dev/null
+++ b/miniapp/miniapp_recipes.hpp
@@ -0,0 +1,46 @@
+#pragma once
+
+#include <cstddef>
+#include <memory>
+#include <stdexcept>
+
+#include "recipe.hpp"
+
+// miniapp-specific recipes
+
+namespace nest {
+namespace mc {
+
+struct probe_distribution {
+ float proportion = 1.f; // what proportion of cells should get probes?
+ bool all_segments = true; // false => soma only
+ bool membrane_voltage = true;
+ bool membrane_current = true;
+};
+
+struct basic_recipe_param {
+ unsigned num_compartments = 1;
+ unsigned num_synapses = 1;
+ std::string synapse_type = "expsyn";
+ float min_connection_delay_ms = 20.0;
+ float mean_connection_delay_ms = 20.75;
+ float syn_weight_per_cell = 0.3;
+};
+
+std::unique_ptr<recipe> make_basic_ring_recipe(
+ cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist = probe_distribution{});
+
+std::unique_ptr<recipe> make_basic_kgraph_recipe(
+ cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist = probe_distribution{});
+
+std::unique_ptr<recipe> make_basic_rgraph_recipe(
+ cell_gid_type ncell,
+ basic_recipe_param param,
+ probe_distribution pdist = probe_distribution{});
+
+} // namespace mc
+} // namespace nest
diff --git a/miniapp/trace_sampler.hpp b/miniapp/trace_sampler.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..26db0a491444db32a49653f861b18dce464c103d
--- /dev/null
+++ b/miniapp/trace_sampler.hpp
@@ -0,0 +1,70 @@
+#pragma once
+
+#include <cstdlib>
+#include <vector>
+
+#include <common_types.hpp>
+#include <cell.hpp>
+#include <util/optional.hpp>
+
+#include <iostream>
+
+namespace nest {
+namespace mc {
+
+template <typename Time=float, typename Value=double>
+struct sample_trace {
+ using time_type = Time;
+ using value_type = Value;
+
+ struct sample_type {
+ time_type time;
+ value_type value;
+ };
+
+ std::string name;
+ std::string units;
+ cell_member_type probe_id;
+ std::vector<sample_type> samples;
+
+ sample_trace() =default;
+ sample_trace(cell_member_type probe_id, const std::string& name, const std::string& units):
+ name(name), units(units), probe_id(probe_id)
+ {}
+};
+
+template <typename Time=float, typename Value=double>
+struct trace_sampler {
+ using time_type = Time;
+ using value_type = Value;
+
+ float next_sample_t() const { return t_next_sample_; }
+
+ util::optional<time_type> operator()(time_type t, value_type v) {
+ if (t<t_next_sample_) {
+ return t_next_sample_;
+ }
+
+ trace_->samples.push_back({t,v});
+ return t_next_sample_+=sample_dt_;
+ }
+
+ trace_sampler(sample_trace<time_type, value_type> *trace, time_type sample_dt, time_type tfrom=0):
+ trace_(trace), sample_dt_(sample_dt), t_next_sample_(tfrom)
+ {}
+
+private:
+ sample_trace<time_type, value_type> *trace_;
+
+ time_type sample_dt_;
+ time_type t_next_sample_;
+};
+
+// with type deduction ...
+template <typename Time, typename Value>
+trace_sampler<Time, Value> make_trace_sampler(sample_trace<Time, Value> *trace, Time sample_dt, Time tfrom=0) {
+ return trace_sampler<Time, Value>(trace, sample_dt, tfrom);
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/nrn/ball_and_3stick.py b/nrn/ball_and_3stick.py
index 15ece3d3da151daf7bdb90f1af1d6cd368352489..4333072216b56d203708430243e28ef008cdfa68 100644
--- a/nrn/ball_and_3stick.py
+++ b/nrn/ball_and_3stick.py
@@ -4,11 +4,21 @@ from matplotlib import pyplot
import numpy as np
import json
import argparse
+import sys
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()
+
+# hack to make things work with nrniv ... -python:
+# throw away args before -python foo.py if -python present.
+
+if '-python' in sys.argv:
+ argv = sys.argv[sys.argv.index('-python')+2:]
+else:
+ argv = sys.argv
+
+args = parser.parse_args(argv)
soma = h.Section(name='soma')
dend = [
diff --git a/nrn/ball_and_stick.py b/nrn/ball_and_stick.py
index 24562337a4baffd7a1ae3a1e7885cb07e06987ec..10932d1adcda27e60e5c8781f491a665b6d8cd8c 100644
--- a/nrn/ball_and_stick.py
+++ b/nrn/ball_and_stick.py
@@ -4,11 +4,21 @@ from matplotlib import pyplot
import numpy as np
import json
import argparse
+import sys
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()
+
+# hack to make things work with nrniv ... -python:
+# throw away args before -python foo.py if -python present.
+
+if '-python' in sys.argv:
+ argv = sys.argv[sys.argv.index('-python')+2:]
+else:
+ argv = sys.argv
+
+args = parser.parse_args(argv)
soma = h.Section(name='soma')
dend = h.Section(name='dend')
diff --git a/nrn/generate_validation.sh b/nrn/generate_validation.sh
index c42b1d5ce66512b8dc5ff46bf13e85f4707bd8ef..532a7cd738d53a1a0f8342d07f09c016b5e24857 100755
--- a/nrn/generate_validation.sh
+++ b/nrn/generate_validation.sh
@@ -1,4 +1,5 @@
python2.7 ./soma.py
python2.7 ./ball_and_stick.py
python2.7 ./ball_and_3stick.py
-python2.7 ./simple_synapse.py
+python2.7 ./simple_synapse.py --synapse exp2
+python2.7 ./simple_synapse.py --synapse exp
diff --git a/nrn/simple_synapse.py b/nrn/simple_synapse.py
index 1779799f716157df995f2321346ad6ff52f1a608..de7bf22d7c84f652148605bdb3581defc249d208 100644
--- a/nrn/simple_synapse.py
+++ b/nrn/simple_synapse.py
@@ -8,7 +8,20 @@ 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()
+parser.add_argument('--synapse', metavar='SYN', dest='synapse',
+ choices=['exp', 'exp2'],
+ default='exp',
+ help='use synapse type SYN (exp or exp2)')
+
+# hack to make things work with nrniv ... -python:
+# throw away args before -python foo.py if -python present.
+
+if '-python' in sys.argv:
+ argv = sys.argv[sys.argv.index('-python')+2:]
+else:
+ argv = sys.argv
+
+args = parser.parse_args(argv)
soma = h.Section(name='soma')
dend = h.Section(name='dend')
@@ -62,8 +75,16 @@ for nseg in [5, 11, 51, 101] :
dend.nseg=nseg
# add a synapse
- syn_ = h.ExpSyn(dend(0.5))
- syn_.tau = 2
+ if args.synapse == 'exp':
+ syn_ = h.ExpSyn(dend(0.5))
+ syn_.tau = 2
+ elif args.synapse == 'exp2':
+ syn_ = h.Exp2Syn(dend(0.5))
+ syn_.tau1 = 0.5 # artificially slow onset
+ syn_.tau2 = 2
+ else:
+ raise RuntimeError('unrecognized synapse type')
+
ncstim = h.NetCon(stim, syn_)
ncstim.delay = 1 # 1 ms delay (arrives at 10ms)
ncstim.weight[0] = 0.04 # NetCon weight is a vector
@@ -127,7 +148,7 @@ end = timer()
print "took ", end-start, " seconds"
# save the spike info as in json format
-fp = open('simple_synapse.json', 'w')
+fp = open('simple_'+args.synapse+'_synapse.json', 'w')
json.dump(results, fp, indent=2)
if args.plot :
diff --git a/nrn/soma.py b/nrn/soma.py
index 3af4f50f3e6e1ae7761591ab68fd32a1187f9c87..9decb2c1700a6fd113530019835dad3e807c2b9b 100644
--- a/nrn/soma.py
+++ b/nrn/soma.py
@@ -4,11 +4,21 @@ from matplotlib import pyplot
import numpy as np
import json
import argparse
+import sys
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()
+
+# hack to make things work with nrniv ... -python:
+# throw away args before -python foo.py if -python present.
+
+if '-python' in sys.argv:
+ argv = sys.argv[sys.argv.index('-python')+2:]
+else:
+ argv = sys.argv
+
+args = parser.parse_args(argv)
if args.plot :
print '-- plotting turned on'
diff --git a/scripts/README.md b/scripts/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e645a762eaf88c5f7196a995502589da44c2f89c
--- /dev/null
+++ b/scripts/README.md
@@ -0,0 +1,81 @@
+tsplot
+------
+
+The `tsplot` script is a wrapper around matplotlib for displaying a collection of
+time series plots.
+
+## Input data
+
+`tsplot` reads timeseries in JSON format, according to the following conventions.
+
+```
+{
+ "units": <units>
+ "name": <name>
+ <other key-value metadata>
+ "data": {
+ "time": [ <time values> ]
+ <trace name>: [ <trace values> ]
+ }
+}
+```
+
+The `data` object must contain numeric arrays, with at least one with the key `time`;
+other members of `data` correspond to traces sampled at the corresponding time values.
+
+The other members of the top level object are regarded as metadata, with some keys
+treated specially:
+ * `units` are used to distinguish different axes for plotting, and the labels for those
+ axes. It's value is either a string, where the specified unit is taken as applying to
+ all included traces, or an object representing a mapping of trace names to their
+ corresponding unit string.
+ * `name` is taken as the title of the corresponding plot, if it is unambiguous.
+ * `label` is ignored: the _label_ for a trace is its name in the data object.
+
+## Operation
+
+The basic usage is simply:
+```
+tsplot data.json ...
+```
+which will produce an interactive plot of the timeseries provided by the provided
+files, with one trace per subplot.
+
+### Grouping
+
+Traces can be gathered on to the same subplot by grouping by metadata with the
+`-g` or `--group` option. To collect all traces with the same value of the key
+'id' and the same units:
+```
+tsplot -g units,id data.json ...
+```
+A subplot can comprise data with to two differint units, and will be plotted
+with two differing vertical axes.
+
+Note that for the purposes of tsplot, the value of the key _label_ is the
+propertu name of the trace in its json representation.
+
+### Restricting data
+
+The `-t` or `--trange` option exlcudes any points that have a time range outside
+that specified. Ranges are given by two numbers separated by a comma, but one or
+the other can be omitted to indicate that there is no bound on that side. For
+example:
+```
+tsplot -t ,100 data.json ...
+```
+will display all points with a time value less than or equal to 100.
+
+Extreme values for data can be automatically excluded and marked on the plot
+with the `-x` or `--exclude` option, taking a parameter _N_. All values in a
+timeseries that lie outside the interval [ _m_ - _Nr_, _m_ + _Nr_ ] are omitted,
+where _m_ is the median of the finite values in the timeseries, and _r_ is
+the 90% interquantile gap, that is, the difference between the 5% and 95% quantile
+of the timeseries data.
+
+### Output to file
+
+Use the `-o` or `--output` option to save the plot as an image, instead of
+displaying it interactively.
+
+
diff --git a/scripts/tsplot b/scripts/tsplot
new file mode 100755
index 0000000000000000000000000000000000000000..81a5da78c06a24d71f727431ef237a18f2a13b0b
--- /dev/null
+++ b/scripts/tsplot
@@ -0,0 +1,363 @@
+#!env python
+#coding: utf-8
+
+import argparse
+import json
+import numpy as np
+import sys
+import re
+import logging
+import matplotlib as M
+import matplotlib.pyplot as P
+from distutils.version import LooseVersion
+from itertools import chain, islice, cycle
+
+# Run-time check for matplot lib version for line style functionality.
+if LooseVersion(M.__version__)<LooseVersion("1.5.0"):
+ logging.critical("require matplotlib version ≥ 1.5.0")
+ sys.exit(1)
+
+# Read timeseries data from multiple files, plot each in one planel, with common
+# time axis, and optionally sharing a vertical axis as governed by the configuration.
+
+def parse_clargs():
+ def float_or_none(s):
+ try: return float(s)
+ except ValueError: return None
+
+ def parse_range_spec(s):
+ l, r = (float_or_none(x) for x in s.split(','))
+ return (l,r)
+
+ P = argparse.ArgumentParser(description='Plot time series data on one or more graphs.')
+ P.add_argument('inputs', metavar='FILE', nargs='+',
+ help='time series data in JSON format')
+ P.add_argument('-t', '--trange', metavar='RANGE', dest='trange',
+ type=parse_range_spec,
+ help='restrict time axis to RANGE (see below)')
+ P.add_argument('-g', '--group', metavar='KEY,...', dest='groupby',
+ type=lambda s: s.split(','),
+ help='plot series with same KEYs on the same axes')
+ P.add_argument('-o', '--output', metavar='FILE', dest='outfile',
+ help='save plot to file FILE')
+ P.add_argument('-x', '--exclude', metavar='NUM', dest='exclude',
+ type=float,
+ help='remove extreme points outside NUM times the 0.9-interquantile range of the median')
+
+ P.epilog = 'A range is specifed by a pair of floating point numbers min,max where '
+ P.epilog += 'either may be omitted to indicate the minimum or maximum of the corresponding '
+ P.epilog += 'values in the data.'
+
+ # modify args to avoid argparse having a fit when it encounters an option
+ # argument of the form '<negative number>,...'
+
+ argsbis = [' '+a if re.match(r'-[\d.]',a) else a for a in sys.argv[1:]]
+ return P.parse_args(argsbis)
+
+def isstring(s):
+ return isinstance(s,str) or isinstance(s,unicode)
+
+def take(n, s):
+ return islice((i for i in s), 0, n)
+
+class TimeSeries:
+ def __init__(self, ts, ys, **kwargs):
+ self.t = np.array(ts)
+ n = self.t.shape[0]
+
+ self.y = np.full_like(self.t, np.nan)
+ ny = min(len(ys), len(self.y))
+ self.y[:ny] = ys[:ny]
+
+ self.meta = dict(kwargs)
+ self.ex_ts = None
+
+ def trestrict(self, bounds):
+ clip = range_meet(self.trange(), bounds)
+ self.t = np.ma.masked_outside(self.t, v1=clip[0], v2=clip[1])
+ self.y = np.ma.masked_array(self.y, mask=self.t.mask)
+
+ def exclude_outliers(self, iqr_factor):
+ yfinite = np.ma.masked_invalid(self.y).compressed()
+ l_, lq, median, uq, u_ = np.percentile(yfinite, [0, 5.0, 50.0, 95.0, 100])
+ lb = median - iqr_factor*(uq-lq)
+ ub = median + iqr_factor*(uq-lq)
+
+ np_err_save = np.seterr(all='ignore')
+ yex = np.ma.masked_where(np.isfinite(self.y)&(self.y<=ub)&(self.y>=lb), self.y)
+ np.seterr(**np_err_save)
+
+ tex = np.ma.masked_array(self.t, mask=yex.mask)
+ self.ex_ts = TimeSeries(tex.compressed(), yex.compressed())
+ self.ex_ts.meta = dict(self.meta)
+
+ self.y = np.ma.filled(np.ma.masked_array(self.y, mask=~yex.mask), np.nan)
+
+ def excluded(self):
+ return self.ex_ts
+
+ def name(self):
+ return self.meta.get('name',"") # value of 'name' attribute in source
+
+ def label(self):
+ return self.meta.get('label',"") # name of column in source
+
+ def units(self):
+ return self.meta.get('units',"")
+
+ def trange(self):
+ return (min(self.t), max(self.t))
+
+ def yrange(self):
+ return (min(self.y), max(self.y))
+
+
+def read_json_timeseries(source):
+ j = json.load(source)
+
+ # Convention:
+ #
+ # Time series data is represented by an object with a subobject 'data' and optionally
+ # other key/value entries comprising metadata for the time series.
+ #
+ # The 'data' object holds one array of numbers 'time' and zero or more other
+ # numeric arrays of sample values corresponding to the values in 'time'. The
+ # names of these other arrays are taken to be the labels for the plots.
+ #
+ # Units can be specified by a top level entry 'units' which is either a string
+ # (units are common for all data series in the object) or by a map that
+ # takes a label to a unit string.
+
+ ts_list = []
+ jdata = j['data']
+ ncol = len(jdata)
+
+ times = jdata['time']
+ nsample = len(times)
+
+ def units(label):
+ try:
+ unitmap = j['units']
+ if isstring(unitmap):
+ return unitmap
+ else:
+ return unitmap[label]
+ except:
+ return ""
+
+ i = 1
+ for key in jdata.keys():
+ if key=="time": continue
+
+ meta = dict(j.items() + {'label': key, 'data': None, 'units': units(key)}.items())
+ del meta['data']
+
+ ts_list.append(TimeSeries(times, jdata[key], **meta))
+
+ return ts_list
+
+def range_join(r, s):
+ return (min(r[0], s[0]), max(r[1], s[1]))
+
+def range_meet(r, s):
+ return (max(r[0], s[0]), min(r[1], s[1]))
+
+
+class PlotData:
+ def __init__(self, key_label=""):
+ self.series = []
+ self.group_key_label = key_label
+
+ def trange(self):
+ return reduce(range_join, [s.trange() for s in self.series])
+
+ def yrange(self):
+ return reduce(range_join, [s.yrange() for s in self.series])
+
+ def name(self):
+ return reduce(lambda n, s: n or s.name(), self.series, "")
+
+ def group_label(self):
+ return self.group_key_label
+
+ def unique_labels(self):
+ # attempt to create unique labels for plots in the group based on
+ # meta data
+ labels = [s.label() for s in self.series]
+ if len(labels)<2:
+ return labels
+
+ n = len(labels)
+ keyset = reduce(lambda k, s: k.union(s.meta.keys()), self.series, set())
+ keyi = iter(keyset)
+ try:
+ while len(set(labels)) != n:
+ k = next(keyi)
+ if k=='label':
+ continue
+
+ vs = [s.meta.get(k,None) for s in self.series]
+ if len(set(vs))==1:
+ continue
+
+ for i in xrange(n):
+ prefix = '' if k=='name' else k+'='
+ if vs[i] is not None:
+ labels[i] += u', '+k+u'='+unicode(vs[i])
+
+ except StopIteration:
+ pass
+
+ return labels
+
+# Input: list of TimeSeries objects; collection of metadata keys to group on
+# Return list of plot info (time series, data extents, metadata), one per plot.
+
+def gather_ts_plots(tss, groupby):
+ group_lookup = {}
+ plot_groups = []
+
+ for ts in tss:
+ key = tuple([ts.meta.get(g) for g in groupby])
+ if key is () or None in key or key not in group_lookup:
+ pretty_key=', '.join([unicode(k)+u'='+unicode(v) for k,v in zip(groupby, key) if v is not None])
+ pd = PlotData(pretty_key)
+ pd.series = [ts]
+ plot_groups.append(pd)
+ group_lookup[key] = len(plot_groups)-1
+ else:
+ plot_groups[group_lookup[key]].series.append(ts)
+
+ return plot_groups
+
+
+def make_palette(cm_name, n, cmin=0, cmax=1):
+ smap = M.cm.ScalarMappable(M.colors.Normalize(cmin/float(cmin-cmax),(cmin-1)/float(cmin-cmax)),
+ M.cm.get_cmap(cm_name))
+ return [smap.to_rgba((2*i+1)/float(2*n)) for i in xrange(n)]
+
+def plot_plots(plot_groups, save=None):
+ nplots = len(plot_groups)
+ plot_groups = sorted(plot_groups, key=lambda g: g.group_label())
+
+ # use same global time scale for all plots
+ trange = reduce(range_join, [g.trange() for g in plot_groups])
+
+ # use group names for titles?
+ group_titles = any((g.group_label() for g in plot_groups))
+
+ figure = P.figure()
+ for i in xrange(nplots):
+ group = plot_groups[i]
+ plot = figure.add_subplot(nplots, 1, i+1)
+
+ title = group.group_label() if group_titles else group.name()
+ plot.set_title(title)
+
+ # y-axis label: use timeseries label and units if only
+ # one series in group, otherwise use a legend with labels,
+ # and units alone on the axes. At most two different unit
+ # axes can be drawn.
+
+ def ylabel(unit):
+ if len(group.series)==1:
+ lab = group.series[0].label()
+ if unit:
+ lab += ' (' + unit + ')'
+ else:
+ lab = unit
+
+ return lab
+
+ uniq_units = list(set([s.units() for s in group.series]))
+ uniq_units.sort()
+ if len(uniq_units)>2:
+ logging.warning('more than two different units on the same plot')
+ uniq_units = uniq_units[:2]
+
+ # store each series in a slot corresponding to one of the units,
+ # together with a best-effort label
+
+ series_by_unit = [[] for i in xrange(len(uniq_units))]
+ unique_labels = group.unique_labels()
+
+ for si in xrange(len(group.series)):
+ s = group.series[si]
+ label = unique_labels[si]
+ try:
+ series_by_unit[uniq_units.index(s.units())].append((s,label))
+ except ValueError:
+ pass
+
+ # TODO: need to find a scheme of colour/line allocation that is
+ # double y-axis AND greyscale friendly.
+
+ palette = \
+ [make_palette(cm, n, 0, 0.5) for
+ cm, n in zip(['hot', 'winter'], [len(x) for x in series_by_unit])]
+
+ lines = cycle(["-",(0,(3,1))])
+
+ first_plot = True
+ for ui in xrange(len(uniq_units)):
+ if not first_plot:
+ plot = plot.twinx()
+
+ axis_color = palette[ui][0]
+ plot.set_ylabel(ylabel(uniq_units[ui]), color=axis_color)
+ for l in plot.get_yticklabels():
+ l.set_color(axis_color)
+
+ plot.get_yaxis().get_major_formatter().set_useOffset(False)
+ plot.get_yaxis().set_major_locator(M.ticker.MaxNLocator(nbins=6))
+
+ plot.set_xlim(trange)
+
+ colours = cycle(palette[ui])
+ line = next(lines)
+ for s, l in series_by_unit[ui]:
+ c = next(colours)
+ plot.plot(s.t, s.y, color=c, ls=line, label=l)
+ # treat exluded points especially
+ ex = s.excluded()
+ if ex is not None:
+ ymin, ymax = s.yrange()
+ plot.plot(ex.t, np.clip(ex.y, ymin, ymax), marker='x', ls='', color=c)
+
+ if first_plot:
+ plot.legend(loc=2, fontsize='small')
+ plot.grid()
+ else:
+ plot.legend(loc=1, fontsize='small')
+
+ first_plot = False
+
+ # adapted from http://stackoverflow.com/questions/6963035
+ axis_ymin = min([ax.get_position().ymin for ax in figure.axes])
+ figure.text(0.5, axis_ymin - float(3)/figure.dpi, 'time', ha='center', va='center')
+ if save:
+ figure.savefig(save)
+ else:
+ P.show()
+
+args = parse_clargs()
+tss = []
+for filename in args.inputs:
+ with open(filename) as f:
+ tss.extend(read_json_timeseries(f))
+
+if args.trange:
+ for ts in tss:
+ ts.trestrict(args.trange)
+
+if args.exclude:
+ for ts in tss:
+ ts.exclude_outliers(args.exclude)
+
+groupby = args.groupby if args.groupby else []
+plots = gather_ts_plots(tss, groupby)
+
+if not args.outfile:
+ M.interactive(False)
+
+plot_plots(plots, save=args.outfile)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7f4a588d933a5ea74342c70dc8b6ac5df95aca8b..8e0db8876fadf3d030b65ccac9a653eedc1c8d37 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -2,8 +2,8 @@ set(HEADERS
swcio.hpp
)
set(BASE_SOURCES
+ common_types_io.cpp
cell.cpp
- mechanism_interface.cpp
parameter_list.cpp
profiling/profiler.cpp
swcio.cpp
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 1d55535ff032c8d66e06a8ea9fa69f52fd87eb91..5456417ddaea0e0fe6c4e777a0fa8ebf3344c2ed 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -90,12 +90,12 @@ bool is_minimal_degree(C const& c)
"is_minimal_degree only applies to integral types"
);
- if(c.size()==0u) {
+ if (c.size()==0u) {
return true;
}
using value_type = typename C::value_type;
- if(c[0] != value_type(0)) {
+ if (c[0] != value_type(0)) {
return false;
}
auto i = value_type(1);
@@ -222,7 +222,7 @@ typename C::value_type find_branch(const C& branch_index,
auto it = std::find_if(
branch_index.begin(), branch_index.end(),
- [nid](const value_type &v) { return v > nid; }
+ [nid](const value_type& v) { return v > nid; }
);
return it - branch_index.begin() - 1;
diff --git a/src/cell.cpp b/src/cell.cpp
index 1a70ca3e3034cf25e9238352cb1c5dfdd97a8c49..e96cae69aca834d61d81c9c53c9488532a3b5b69 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -9,7 +9,7 @@ int find_compartment_index(
segment_location const& location,
compartment_model const& graph
) {
- EXPECTS(location.segment<graph.segment_index.size());
+ EXPECTS(unsigned(location.segment)<graph.segment_index.size());
const auto& si = graph.segment_index;
const auto seg = location.segment;
@@ -26,7 +26,7 @@ cell::cell()
parents_.push_back(0);
}
-int cell::num_segments() const
+cell::size_type cell::num_segments() const
{
return segments_.size();
}
@@ -75,9 +75,9 @@ cable_segment* cell::add_cable(cell::index_type parent, segment_ptr&& cable)
return segments_.back()->as_cable();
}
-segment* cell::segment(int index)
+segment* cell::segment(index_type index)
{
- if(index<0 || index>=num_segments()) {
+ if (index>=num_segments()) {
throw std::out_of_range(
"attempt to access a segment with invalid index"
);
@@ -85,9 +85,9 @@ segment* cell::segment(int index)
return segments_[index].get();
}
-segment const* cell::segment(int index) const
+segment const* cell::segment(index_type index) const
{
- if(index<0 || index>=num_segments()) {
+ if (index>=num_segments()) {
throw std::out_of_range(
"attempt to access a segment with invalid index"
);
@@ -109,7 +109,7 @@ soma_segment* cell::soma()
return nullptr;
}
-cable_segment* cell::cable(int index)
+cable_segment* cell::cable(index_type index)
{
if(index>0 && index<num_segments()) {
return segment(index)->as_cable();
@@ -146,9 +146,9 @@ std::vector<segment_ptr> const& cell::segments() const
return segments_;
}
-std::vector<int> cell::compartment_counts() const
+std::vector<cell::size_type> cell::compartment_counts() const
{
- std::vector<int> comp_count;
+ std::vector<size_type> comp_count;
comp_count.reserve(num_segments());
for(auto const& s : segments()) {
comp_count.push_back(s->num_compartments());
@@ -156,10 +156,10 @@ std::vector<int> cell::compartment_counts() const
return comp_count;
}
-size_t cell::num_compartments() const
+cell::size_type cell::num_compartments() const
{
auto n = 0u;
- for(auto &s : segments_) {
+ for(auto& s : segments_) {
n += s->num_compartments();
}
return n;
@@ -196,21 +196,11 @@ void cell::add_detector(segment_location loc, double threshold)
spike_detectors_.push_back({loc, threshold});
}
-std::vector<int> const& cell::segment_parents() const
+std::vector<cell::index_type> const& cell::segment_parents() const
{
return parents_;
}
-void cell::add_synapse(segment_location loc)
-{
- synapses_.push_back(loc);
-}
-
-const std::vector<segment_location>& cell::synapses() const
-{
- return synapses_;
-}
-
// Rough and ready comparison of two cells.
// We don't use an operator== because equality of two cells is open to
// interpretation. For example, it is possible to have two viable representations
@@ -222,24 +212,22 @@ const std::vector<segment_location>& cell::synapses() const
// - number of compartments in each segment
bool cell_basic_equality(cell const& lhs, cell const& rhs)
{
- if(lhs.num_segments() != rhs.num_segments()) {
+ if (lhs.num_segments() != rhs.num_segments()) {
return false;
}
- if(lhs.segment_parents() != rhs.segment_parents()) {
+ if (lhs.segment_parents() != rhs.segment_parents()) {
return false;
}
- for(auto i=0; i<lhs.num_segments(); ++i) {
+ for (cell::index_type i=0; i<lhs.num_segments(); ++i) {
// a quick and dirty test
auto& l = *lhs.segment(i);
auto& r = *rhs.segment(i);
- if(l.kind() != r.kind()) return false;
- if(l.area() != r.area()) return false;
- if(l.volume() != r.volume()) return false;
- if(l.as_cable()) {
- if( l.as_cable()->num_compartments()
- != r.as_cable()->num_compartments())
- {
+ if (l.kind() != r.kind()) return false;
+ if (l.area() != r.area()) return false;
+ if (l.volume() != r.volume()) return false;
+ if (l.as_cable()) {
+ if (l.as_cable()->num_compartments() != r.as_cable()->num_compartments()) {
return false;
}
}
diff --git a/src/cell.hpp b/src/cell.hpp
index 744c9970e26af7f439f9cee7cfb2f88c96014624..746549289f1b718433388e92f5f1ca9d13d433c1 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -5,8 +5,9 @@
#include <thread>
#include <vector>
-#include "segment.hpp"
+#include "common_types.hpp"
#include "cell_tree.hpp"
+#include "segment.hpp"
#include "stimulus.hpp"
#include "util/debug.hpp"
@@ -17,12 +18,12 @@ namespace mc {
/// description
struct compartment_model {
cell_tree tree;
- std::vector<int> parent_index;
- std::vector<int> segment_index;
+ std::vector<cell_tree::int_type> parent_index;
+ std::vector<cell_tree::int_type> segment_index;
};
struct segment_location {
- segment_location(int s, double l)
+ segment_location(cell_lid_type s, double l)
: segment(s), position(l)
{
EXPECTS(position>=0. && position<=1.);
@@ -30,7 +31,7 @@ struct segment_location {
friend bool operator==(segment_location l, segment_location r) {
return l.segment==r.segment && l.position==r.position;
}
- int segment;
+ cell_lid_type segment;
double position;
};
@@ -44,22 +45,29 @@ enum class probeKind {
membrane_current
};
+struct probe_spec {
+ segment_location location;
+ probeKind kind;
+};
+
/// high-level abstract representation of a cell and its segments
class cell {
public:
-
- // types
- using index_type = int;
+ using index_type = cell_lid_type;
+ using size_type = cell_local_size_type;
using value_type = double;
using point_type = point<value_type>;
- struct probe_instance {
+
+ struct synapse_instance {
segment_location location;
- probeKind kind;
+ parameter_list mechanism;
};
+
struct stimulus_instance {
segment_location location;
i_clamp clamp;
};
+
struct detector_instance {
segment_location location;
double threshold;
@@ -84,12 +92,12 @@ public:
cable_segment* add_cable(index_type parent, Args ...args);
/// the number of segments in the cell
- int num_segments() const;
+ size_type num_segments() const;
bool has_soma() const;
- class segment* segment(int index);
- class segment const* segment(int index) const;
+ class segment* segment(index_type index);
+ class segment const* segment(index_type index) const;
/// access pointer to the soma
/// returns nullptr if the cell has no soma
@@ -98,7 +106,7 @@ public:
/// 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);
+ cable_segment* cable(index_type index);
/// the volume of the cell
value_type volume() const;
@@ -107,16 +115,16 @@ public:
value_type area() const;
/// the total number of compartments over all segments
- size_t num_compartments() const;
+ size_type 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;
+ std::vector<index_type> const& segment_parents() const;
/// return a vector with the compartment count for each segment in the cell
- std::vector<int> compartment_counts() const;
+ std::vector<size_type> compartment_counts() const;
compartment_model model() const;
@@ -138,9 +146,13 @@ public:
//////////////////
// synapses
//////////////////
- void add_synapse(segment_location loc);
-
- const std::vector<segment_location>& synapses() const;
+ void add_synapse(segment_location loc, parameter_list p)
+ {
+ synapses_.push_back(synapse_instance{loc, std::move(p)});
+ }
+ const std::vector<synapse_instance>& synapses() const {
+ return synapses_;
+ }
//////////////////
// spike detectors
@@ -160,12 +172,12 @@ public:
//////////////////
// probes
//////////////////
- index_type add_probe(segment_location loc, probeKind kind) {
- probes_.push_back({loc, kind});
+ index_type add_probe(probe_spec p) {
+ probes_.push_back(p);
return probes_.size()-1;
}
- const std::vector<probe_instance>&
+ const std::vector<probe_spec>&
probes() const { return probes_; }
private:
@@ -180,13 +192,13 @@ private:
std::vector<stimulus_instance> stimulii_;
// the synapses
- std::vector<segment_location> synapses_;
+ std::vector<synapse_instance> synapses_;
// the sensors
std::vector<detector_instance> spike_detectors_;
// the probes
- std::vector<probe_instance> probes_;
+ std::vector<probe_spec> probes_;
};
// Checks that two cells have the same
diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index 70bfdeb13123c780744c5875ae125bf631bf6142..e70d4b0bb5826b6ccde449f63a7fad6dab29c783 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -1,113 +1,105 @@
#pragma once
#include <cstdint>
+#include <functional>
#include <vector>
#include <cell.hpp>
+#include <common_types.hpp>
#include <event_queue.hpp>
-#include <communication/spike.hpp>
-#include <communication/spike_source.hpp>
+#include <spike.hpp>
+#include <spike_source.hpp>
#include <profiling/profiler.hpp>
namespace nest {
namespace mc {
-// samplers take a time and sample value, and return an optional time
-// for the next desired sample.
-
-struct sampler {
- using index_type = int;
- using time_type = float;
- using value_type = double;
-
- index_type probe_gid; // samplers are attached to probes
- std::function<util::optional<time_type>(time_type, value_type)> sample;
-};
-
template <typename Cell>
class cell_group {
public:
- using index_type = uint32_t;
+ using index_type = cell_gid_type;
using cell_type = Cell;
using value_type = typename cell_type::value_type;
using size_type = typename cell_type::value_type;
using spike_detector_type = spike_detector<Cell>;
+ using source_id_type = cell_member_type;
+
+ using time_type = float;
+ using sampler_function = std::function<util::optional<time_type>(time_type, double)>;
struct spike_source_type {
- index_type index;
+ source_id_type source_id;
spike_detector_type source;
};
cell_group() = default;
- cell_group(const cell& c) :
- cell_{c}
+ cell_group(cell_gid_type gid, const cell& c) :
+ gid_base_{gid}, cell_{c}
{
- cell_.voltage()(memory::all) = -65.;
- cell_.initialize();
+ initialize_cells();
+
+ // Create spike detectors and associate them with globally unique source ids,
+ // as specified by cell gid and cell-local zero-based index.
+
+ cell_gid_type source_gid = gid_base_;
+ cell_lid_type source_lid = 0u;
for (auto& d : c.detectors()) {
- spike_sources_.push_back( {
- 0u, spike_detector_type(cell_, d.location, d.threshold, 0.f)
+ cell_member_type source_id{source_gid, source_lid++};
+
+ spike_sources_.push_back({
+ source_id, spike_detector_type(cell_, d.location, d.threshold, 0.f)
});
}
}
- void set_source_gids(index_type gid) {
- for (auto& s : spike_sources_) {
- s.index = gid++;
+ void reset() {
+ remove_samplers();
+ initialize_cells();
+ for (auto& spike_source: spike_sources_) {
+ spike_source.source.reset(cell_, 0.f);
}
}
- void set_target_gids(index_type lid) {
- first_target_gid_ = lid;
- }
-
- index_type num_probes() const {
- return cell_.num_probes();
- }
-
- void set_probe_gids(index_type gid) {
- first_probe_gid_ = gid;
- }
-
- std::pair<index_type, index_type> probe_gid_range() const {
- return { first_probe_gid_, first_probe_gid_+cell_.num_probes() };
- }
-
- void advance(double tfinal, double dt) {
+ void advance(time_type tfinal, time_type dt) {
while (cell_.time()<tfinal) {
// take any pending samples
- float cell_time = cell_.time();
+ time_type cell_time = cell_.time();
- nest::mc::util::profiler_enter("sampling");
+ PE("sampling");
while (auto m = sample_events_.pop_if_before(cell_time)) {
- auto& sampler = samplers_[m->sampler_index];
- EXPECTS((bool)sampler.sample);
+ auto& sampler_spec = samplers_[m->sampler_index];
+ EXPECTS((bool)sampler_spec.sampler);
- index_type probe_index = sampler.probe_gid-first_probe_gid_;
- auto next = sampler.sample(cell_.time(), cell_.probe(probe_index));
+ index_type probe_index = sampler_spec.probe_id.index;
+ auto next = sampler_spec.sampler(cell_.time(), cell_.probe(probe_index));
if (next) {
m->time = std::max(*next, cell_time);
sample_events_.push(*m);
}
}
- nest::mc::util::profiler_leave();
+ PL();
// look for events in the next time step
- auto tstep = std::min(tfinal, cell_.time()+dt);
+ time_type tstep = cell_.time()+dt;
+ tstep = std::min(tstep, tfinal);
+
auto next = events_.pop_if_before(tstep);
- auto tnext = next ? next->time: tstep;
+ time_type tnext = next ? next->time: tstep;
// integrate cell state
cell_.advance(tnext - cell_.time());
+ if (!cell_.is_physical_solution()) {
+ std::cerr << "warning: solution out of bounds\n";
+ }
- nest::mc::util::profiler_enter("events");
+ PE("events");
// check for new spikes
for (auto& s : spike_sources_) {
if (auto spike = s.source.test(cell_, cell_.time())) {
- spikes_.push_back({s.index, spike.get()});
+ spikes_.push_back({s.source_id, spike.get()});
}
}
@@ -121,23 +113,20 @@ public:
cell_.apply_event(e.get());
}
}
- nest::mc::util::profiler_leave();
+ PL();
}
}
template <typename R>
- void enqueue_events(R events) {
+ void enqueue_events(const R& events) {
for (auto e : events) {
- e.target -= first_target_gid_;
events_.push(e);
}
}
- const std::vector<communication::spike<index_type>>&
- spikes() const {
- return spikes_;
- }
+ const std::vector<spike<source_id_type, time_type>>&
+ spikes() const { return spikes_; }
cell_type& cell() { return cell_; }
const cell_type& cell() const { return cell_; }
@@ -151,14 +140,25 @@ public:
spikes_.clear();
}
- void add_sampler(const sampler& s, float start_time = 0) {
+ void add_sampler(cell_member_type probe_id, sampler_function s, time_type start_time = 0) {
auto sampler_index = uint32_t(samplers_.size());
- samplers_.push_back(s);
+ samplers_.push_back({probe_id, s});
sample_events_.push({sampler_index, start_time});
}
+ void remove_samplers() {
+ sample_events_.clear();
+ samplers_.clear();
+ }
+
private:
+ void initialize_cells() {
+ cell_.voltage()(memory::all) = -65.;
+ cell_.initialize();
+ }
+ /// gid of first cell in group
+ cell_gid_type gid_base_;
/// the lowered cell state (e.g. FVM) of the cell
cell_type cell_;
@@ -167,22 +167,27 @@ private:
std::vector<spike_source_type> spike_sources_;
//. spikes that are generated
- std::vector<communication::spike<index_type>> spikes_;
+ std::vector<spike<source_id_type, time_type>> spikes_;
/// pending events to be delivered
- event_queue<postsynaptic_spike_event> events_;
+ event_queue<postsynaptic_spike_event<time_type>> events_;
/// pending samples to be taken
- event_queue<sample_event> sample_events_;
+ event_queue<sample_event<time_type>> sample_events_;
/// the global id of the first target (e.g. a synapse) in this group
index_type first_target_gid_;
-
+
/// the global id of the first probe in this group
index_type first_probe_gid_;
+ struct sampler_entry {
+ cell_member_type probe_id;
+ sampler_function sampler;
+ };
+
/// collection of samplers to be run against probes in this group
- std::vector<sampler> samplers_;
+ std::vector<sampler_entry> samplers_;
};
} // namespace mc
diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp
index 4d9f70b1a934ab3b5bdb840d6bb0bd1b8419dfec..ce06d69a8b8ac7427b991e7e54ac2bc764a052e7 100644
--- a/src/cell_tree.hpp
+++ b/src/cell_tree.hpp
@@ -10,6 +10,8 @@
#include <vector>
#include <vector/include/Vector.hpp>
+
+#include "common_types.hpp"
#include "tree.hpp"
#include "util.hpp"
@@ -29,38 +31,42 @@ namespace mc {
/// flexibility in choosing the root.
class cell_tree {
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 = int;
+
+public:
+ using int_type = cell_lid_type;
+ using size_type = cell_local_size_type;
+
using index_type = memory::HostVector<int_type>;
using view_type = index_type::view_type;
using const_view_type = index_type::const_view_type;
+ using tree = nest::mc::tree<int_type, size_type>;
+ static constexpr int_type no_parent = tree::no_parent;
+
/// default empty constructor
cell_tree() = default;
/// construct from a parent index
- cell_tree(std::vector<int> const& parent_index)
+ cell_tree(std::vector<int_type> const& parent_index)
{
// handle the case of an empty parent list, which implies a single-compartment model
if(parent_index.size()>0) {
tree_ = tree(parent_index);
}
else {
- tree_ = tree(std::vector<int>({0}));
+ tree_ = tree(std::vector<int_type>({0}));
}
}
/// construct from a tree
// copy constructor
- cell_tree(tree const& t, int s)
+ cell_tree(tree const& t, int_type s)
: tree_(t),
soma_(s)
{ }
// move constructor
- cell_tree(tree&& t, int s)
+ cell_tree(tree&& t, int_type s)
: tree_(std::move(t)),
soma_(s)
{ }
@@ -129,12 +135,12 @@ public :
}
/// returns the number of child segments of segment b
- size_t num_children(size_t b) const {
+ size_type num_children(int_type b) const {
return tree_.num_children(b);
}
/// returns a list of the children of segment b
- const_view_type children(size_t b) const {
+ const_view_type children(int_type b) const {
return tree_.children(b);
}
@@ -162,7 +168,7 @@ public :
index_type depth_from_leaf()
{
tree::index_type depth(num_segments());
- depth_from_leaf(depth, 0);
+ depth_from_leaf(depth, int_type{0});
return depth;
}
@@ -170,7 +176,7 @@ public :
{
tree::index_type depth(num_segments());
depth[0] = 0;
- depth_from_root(depth, 1);
+ depth_from_root(depth, int_type{1});
return depth;
}
@@ -179,12 +185,11 @@ private :
/// helper type for sub-tree computation
/// use in balance()
struct sub_tree {
- sub_tree(int r, int diam, int dpth)
- : root(r), diameter(diam), depth(dpth)
+ sub_tree(int_type r, int_type diam, int_type dpth):
+ root(r), diameter(diam), depth(dpth)
{}
- void set(int r, int diam, int dpth)
- {
+ void set(int r, int diam, int dpth) {
root = r;
diameter = diam;
depth = dpth;
@@ -198,15 +203,15 @@ private :
"]";
}
- int root;
- int diameter;
- int depth;
+ int_type root;
+ int_type diameter;
+ int_type depth;
};
/// returns the index of the segment that would minimise the depth of the
/// tree if used as the root segment
int_type find_minimum_root() {
- if(num_segments()==1) {
+ if (num_segments()==1) {
return 0;
}
@@ -229,14 +234,14 @@ private :
// walk has been completed to the root node, the node that has been
// selected will be the root of the sub-tree with the largest diameter.
sub_tree max_sub_tree(0, 0, 0);
- auto distance_from_max_leaf = 1;
+ int_type distance_from_max_leaf = 1;
auto pnt = max_leaf;
auto pos = parent(max_leaf);
- while(pos != -1) {
+ while(pos != no_parent) {
for(auto c : children(pos)) {
if(c!=pnt) {
auto diameter = depth[c] + 1 + distance_from_max_leaf;
- if(diameter>max_sub_tree.diameter) {
+ if (diameter>max_sub_tree.diameter) {
max_sub_tree.set(pos, diameter, distance_from_max_leaf);
}
}
@@ -266,7 +271,7 @@ private :
return new_root;
}
- int_type depth_from_leaf(index_type& depth, int segment)
+ int_type depth_from_leaf(index_type& depth, int_type segment)
{
int_type max_depth = 0;
for(auto c : children(segment)) {
@@ -276,7 +281,7 @@ private :
return max_depth+1;
}
- void depth_from_root(index_type& depth, int segment)
+ void depth_from_root(index_type& depth, int_type segment)
{
auto d = depth[parent(segment)] + 1;
depth[segment] = d;
diff --git a/src/common_types.hpp b/src/common_types.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..da81ed5d65368af60095d8b9dab40aeede29467a
--- /dev/null
+++ b/src/common_types.hpp
@@ -0,0 +1,55 @@
+#pragma once
+
+/*
+ * Common definitions for index types etc. across prototype simulator
+ * library. (Expect analogues in future versions to be template parameters?)
+ */
+
+#include <iosfwd>
+#include <type_traits>
+
+#include <util/lexcmp_def.hpp>
+
+namespace nest {
+namespace mc {
+
+// For identifying cells globally.
+
+using cell_gid_type = std::uint32_t;
+
+// For sizes of collections of cells.
+
+using cell_size_type = typename std::make_unsigned<cell_gid_type>::type;
+
+// For indexes into cell-local data.
+//
+// Local indices for items within a particular cell-local collection should be
+// zero-based and numbered contiguously.
+
+using cell_lid_type = std::uint32_t;
+
+// For counts of cell-local data.
+
+using cell_local_size_type = typename std::make_unsigned<cell_lid_type>::type;
+
+// For global identification of an item of cell local data.
+//
+// Items of cell_member_type must:
+//
+// * be associated with a unique cell, identified by the member `gid`
+// (see: cell_gid_type);
+//
+// * identify an item within a cell-local collection by the member `index`
+// (see: cell_lid_type).
+
+struct cell_member_type {
+ cell_gid_type gid;
+ cell_lid_type index;
+};
+
+DEFINE_LEXICOGRAPHIC_ORDERING(cell_member_type,(a.gid,a.index),(b.gid,b.index))
+
+} // namespace mc
+} // namespace nest
+
+std::ostream& operator<<(std::ostream& O, nest::mc::cell_member_type m);
diff --git a/src/common_types_io.cpp b/src/common_types_io.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ad6ca540b80e1ce296c8ba2f4eceb52cd9ebfd63
--- /dev/null
+++ b/src/common_types_io.cpp
@@ -0,0 +1,8 @@
+#include <iostream>
+
+#include <common_types.hpp>
+
+std::ostream& operator<<(std::ostream& O, nest::mc::cell_member_type m) {
+ return O << m.gid << ':' << m.index;
+}
+
diff --git a/src/communication/communicator.hpp b/src/communication/communicator.hpp
index 34185ac0fbd795294211ced7190d8b84a3e70874..8c72430f5c4646ce7747c76cac088bee65a58df8 100644
--- a/src/communication/communicator.hpp
+++ b/src/communication/communicator.hpp
@@ -5,12 +5,13 @@
#include <vector>
#include <random>
-#include <communication/spike.hpp>
-#include <threading/threading.hpp>
+#include <spike.hpp>
+#include <util/double_buffer.hpp>
#include <algorithms.hpp>
+#include <connection.hpp>
#include <event_queue.hpp>
-
-#include "connection.hpp"
+#include <spike.hpp>
+#include <util/debug.hpp>
namespace nest {
namespace mc {
@@ -26,90 +27,55 @@ namespace communication {
// Once all connections have been specified, the construct() method can be used
// to build the data structures required for efficient spike communication and
// event generation.
-template <typename CommunicationPolicy>
+template <typename Time, typename CommunicationPolicy>
class communicator {
public:
- using id_type = uint32_t;
using communication_policy_type = CommunicationPolicy;
+ using id_type = cell_gid_type;
+ using time_type = Time;
+ using spike_type = spike<cell_member_type, time_type>;
+ using connection_type = connection<time_type>;
- using spike_type = spike<id_type>;
+ /// per-cell group lists of events to be delivered
+ using event_queue =
+ std::vector<postsynaptic_spike_event<time_type>>;
communicator() = default;
- communicator(id_type n_groups, std::vector<id_type> target_counts) :
- num_groups_local_(n_groups),
- num_targets_local_(target_counts.size())
- {
- target_map_ = nest::mc::algorithms::make_index(target_counts);
- num_targets_local_ = target_map_.back();
-
- // create an event queue for each target group
- events_.resize(num_groups_local_);
-
- // make maps for converting lid to gid
- target_gid_map_ = communication_policy_.make_map(num_targets_local_);
- group_gid_map_ = communication_policy_.make_map(num_groups_local_);
+ // TODO
+ // for now, still assuming one-to-one association cells <-> groups,
+ // so that 'group' gids as represented by their first cell gid are
+ // contiguous.
+ communicator(id_type cell_from, id_type cell_to):
+ cell_gid_from_(cell_from), cell_gid_to_(cell_to)
+ {}
- // transform the target ids from lid to gid
- auto first_target = target_gid_map_[domain_id()];
- for (auto &id : target_map_) {
- id += first_target;
- }
- }
-
- id_type target_gid_from_group_lid(id_type lid) const {
- EXPECTS(lid<num_groups_local_);
- return target_map_[lid];
- }
-
- id_type group_gid_from_group_lid(id_type lid) const {
- EXPECTS(lid<num_groups_local_);
- return group_gid_map_[domain_id()] + lid;
+ cell_local_size_type num_groups_local() const
+ {
+ return cell_gid_to_-cell_gid_from_;
}
- void add_connection(connection con) {
- EXPECTS(is_local_target(con.destination()));
+ void add_connection(connection_type con) {
+ EXPECTS(is_local_cell(con.destination().gid));
connections_.push_back(con);
}
- bool is_local_target(id_type gid) {
- return gid>=target_gid_map_[domain_id()]
- && gid<target_gid_map_[domain_id()+1];
- }
-
- bool is_local_group(id_type gid) {
- return gid>=group_gid_map_[domain_id()]
- && gid<group_gid_map_[domain_id()+1];
- }
-
- /// return the global id of the first group in domain d
- /// the groups in domain d are in the contiguous half open range
- /// [domain_first_group(d), domain_first_group(d+1))
- id_type group_gid_first(int d) const {
- return group_gid_map_[d];
- }
-
- id_type target_lid(id_type gid) {
- EXPECTS(is_local_group(gid));
-
- return gid - target_gid_map_[domain_id()];
- }
-
- id_type group_lid(id_type gid) {
- EXPECTS(is_local_group(gid));
-
- return gid - group_gid_map_[domain_id()];
+ /// returns true if the cell with gid is on the domain of the caller
+ bool is_local_cell(id_type gid) const {
+ return gid>=cell_gid_from_ && gid<cell_gid_to_;
}
- // builds the optimized data structure
+ /// builds the optimized data structure
+ /// must be called after all connections have been added
void construct() {
if (!std::is_sorted(connections_.begin(), connections_.end())) {
std::sort(connections_.begin(), connections_.end());
}
}
- float min_delay() {
- auto local_min = std::numeric_limits<float>::max();
+ /// the minimum delay of all connections in the global network.
+ time_type min_delay() {
+ auto local_min = std::numeric_limits<time_type>::max();
for (auto& con : connections_) {
local_min = std::min(local_min, con.delay());
}
@@ -117,50 +83,22 @@ public:
return communication_policy_.min(local_min);
}
- // return the local group index of the group which hosts the target with
- // global id gid
- id_type local_group_from_global_target(id_type gid) {
- // assert that gid is in range
- EXPECTS(is_local_target(gid));
-
- return
- std::distance(
- target_map_.begin(),
- std::upper_bound(target_map_.begin(), target_map_.end(), gid)
- ) - 1;
- }
-
- void add_spike(spike_type s) {
- thread_spikes().push_back(s);
- }
-
- void add_spikes(const std::vector<spike_type>& s) {
- auto& v = thread_spikes();
- v.insert(v.end(), s.begin(), s.end());
- }
-
- std::vector<spike_type>& thread_spikes() {
- return thread_spikes_.local();
- }
-
- void exchange() {
-
- // global all-to-all to gather a local copy of the global spike list
- // on each node
- //profiler_.enter("global exchange");
- auto global_spikes = communication_policy_.gather_spikes(local_spikes());
+ /// Perform exchange of spikes.
+ ///
+ /// Takes as input the list of local_spikes that were generated on the calling domain.
+ ///
+ /// Returns a vector of event queues, with one queue for each local cell group. The
+ /// events in each queue are all events that must be delivered to targets in that cell
+ /// group as a result of the global spike exchange.
+ std::vector<event_queue> exchange(const std::vector<spike_type>& local_spikes) {
+ // global all-to-all to gather a local copy of the global spike list on each node.
+ auto global_spikes = communication_policy_.gather_spikes( local_spikes );
num_spikes_ += global_spikes.size();
- clear_thread_spike_buffers();
- //profiler_.leave();
- for (auto& q : events_) {
- q.clear();
- }
-
- //profiler_.enter("events");
+ // check each global spike in turn to see it generates local events.
+ // if so, make the events and insert them into the appropriate event list.
+ auto queues = std::vector<event_queue>(num_groups_local());
- //profiler_.enter("make events");
- // check all global spikes to see if they will generate local events
for (auto spike : global_spikes) {
// search for targets
auto targets =
@@ -170,35 +108,19 @@ public:
// generate an event for each target
for (auto it=targets.first; it!=targets.second; ++it) {
- auto gidx = local_group_from_global_target(it->destination());
-
- events_[gidx].push_back(it->make_event(spike));
+ auto gidx = cell_group_index(it->destination().gid);
+ queues[gidx].push_back(it->make_event(spike));
}
}
- //profiler_.leave(); // make events
-
- //profiler_.leave(); // event generation
+ return queues;
}
- uint64_t num_spikes() const
- {
- return num_spikes_;
- }
-
- int domain_id() const {
- return communication_policy_.id();
- }
+ /// Returns the total number of global spikes over the duration
+ /// of the simulation
+ uint64_t num_spikes() const { return num_spikes_; }
- int num_domains() const {
- return communication_policy_.size();
- }
-
- const std::vector<postsynaptic_spike_event>& queue(int i) const {
- return events_[i];
- }
-
- const std::vector<connection>& connections() const {
+ const std::vector<connection_type>& connections() const {
return connections_;
}
@@ -206,66 +128,20 @@ public:
return communication_policy_;
}
- const std::vector<id_type>& local_target_map() const {
- return target_map_;
- }
-
- std::vector<spike_type> local_spikes() {
- std::vector<spike_type> spikes;
- for (auto& v : thread_spikes_) {
- spikes.insert(spikes.end(), v.begin(), v.end());
- }
- return spikes;
- }
-
- void clear_thread_spike_buffers() {
- for (auto& v : thread_spikes_) {
- v.clear();
- }
- }
-
private:
+ std::size_t cell_group_index(cell_gid_type cell_gid) const {
+ // this will be more elaborate when there is more than one cell per cell group
+ EXPECTS(cell_gid>=cell_gid_from_ && cell_gid<cell_gid_to_);
+ return cell_gid-cell_gid_from_;
+ }
- //
- // both of these can be fixed with double buffering
- //
- // FIXME : race condition on the thread_spikes_ buffers when exchange() modifies/access them
- // ... other threads will be pushing to them simultaneously
- // FIXME : race condition on the group-specific event queues when exchange pushes to them
- // ... other threads will be accessing them to update their event queues
-
- // thread private storage for accumulating spikes
- using local_spike_store_type =
- nest::mc::threading::enumerable_thread_specific<std::vector<spike_type>>;
- local_spike_store_type thread_spikes_;
-
- std::vector<connection> connections_;
- std::vector<std::vector<postsynaptic_spike_event>> events_;
-
- // local target group i has targets in the half open range
- // [target_map_[i], target_map_[i+1])
- std::vector<id_type> target_map_;
-
- // for keeping track of how time is spent where
- //util::Profiler profiler_;
-
- // the number of groups and targets handled by this communicator
- id_type num_groups_local_;
- id_type num_targets_local_;
-
- // index maps for the global distribution of groups and targets
-
- // communicator i has the groups in the half open range :
- // [group_gid_map_[i], group_gid_map_[i+1])
- std::vector<id_type> group_gid_map_;
-
- // communicator i has the targets in the half open range :
- // [target_gid_map_[i], target_gid_map_[i+1])
- std::vector<id_type> target_gid_map_;
+ std::vector<connection_type> connections_;
communication_policy_type communication_policy_;
uint64_t num_spikes_ = 0u;
+ id_type cell_gid_from_;
+ id_type cell_gid_to_;
};
} // namespace communication
diff --git a/src/communication/connection.hpp b/src/communication/connection.hpp
deleted file mode 100644
index ee1acfd72dc21cffe01f28239bca29526366df31..0000000000000000000000000000000000000000
--- a/src/communication/connection.hpp
+++ /dev/null
@@ -1,70 +0,0 @@
-#pragma once
-
-#include <cstdint>
-
-#include <event_queue.hpp>
-#include <communication/spike.hpp>
-
-namespace nest {
-namespace mc {
-namespace communication {
-
-class connection {
-public:
- using id_type = uint32_t;
- connection(id_type src, id_type dest, float w, float d) :
- source_(src),
- destination_(dest),
- weight_(w),
- delay_(d)
- {}
-
- float weight() const { return weight_; }
- float delay() const { return delay_; }
-
- id_type source() const { return source_; }
- id_type destination() const { return destination_; }
-
- postsynaptic_spike_event make_event(spike<id_type> s) {
- return {destination_, s.time + delay_, weight_};
- }
-
-private:
-
- id_type source_;
- id_type destination_;
- float weight_;
- float delay_;
-};
-
-// connections are sorted by source id
-// these operators make for easy interopability with STL algorithms
-
-static inline
-bool operator< (connection lhs, connection rhs) {
- return lhs.source() < rhs.source();
-}
-
-static inline
-bool operator< (connection lhs, connection::id_type rhs) {
- return lhs.source() < rhs;
-}
-
-static inline
-bool operator< (connection::id_type lhs, connection rhs) {
- return lhs < rhs.source();
-}
-
-} // namespace communication
-} // namespace mc
-} // namespace nest
-
-static inline
-std::ostream& operator<<(std::ostream& o, nest::mc::communication::connection const& con) {
- char buff[512];
- snprintf(
- buff, sizeof(buff), "con [%10u -> %10u : weight %8.4f, delay %8.4f]",
- con.source(), con.destination(), con.weight(), con.delay()
- );
- return o << buff;
-}
diff --git a/src/communication/mpi_global_policy.hpp b/src/communication/mpi_global_policy.hpp
index 7409585447ae888eec2aa47ddc1fa0f137df7d6e..0e9ec33370695990e2e2fa0d8d03d4a1e8b26c26 100644
--- a/src/communication/mpi_global_policy.hpp
+++ b/src/communication/mpi_global_policy.hpp
@@ -4,24 +4,23 @@
#error "mpi_global_policy.hpp should only be compiled in a WITH_MPI build"
#endif
+#include <cstdint>
#include <type_traits>
#include <vector>
-#include <cstdint>
-
-#include <communication/spike.hpp>
-#include <communication/mpi.hpp>
#include <algorithms.hpp>
+#include <common_types.hpp>
+#include <communication/mpi.hpp>
+#include <spike.hpp>
namespace nest {
namespace mc {
namespace communication {
struct mpi_global_policy {
- using id_type = uint32_t;
-
- std::vector<spike<id_type>>
- static gather_spikes(const std::vector<spike<id_type>>& local_spikes) {
+ template <typename I, typename T>
+ static std::vector<spike<I, T>>
+ gather_spikes(const std::vector<spike<I,T>>& local_spikes) {
return mpi::gather_all(local_spikes);
}
diff --git a/src/communication/serial_global_policy.hpp b/src/communication/serial_global_policy.hpp
index a1c7fee9671cb79768e255a710e0082386302e8a..9af5eae3fd0dd84be71cb6346018c607fa9ab3df 100644
--- a/src/communication/serial_global_policy.hpp
+++ b/src/communication/serial_global_policy.hpp
@@ -1,21 +1,19 @@
#pragma once
+#include <cstdint>
#include <type_traits>
#include <vector>
-#include <cstdint>
-
-#include <communication/spike.hpp>
+#include <spike.hpp>
namespace nest {
namespace mc {
namespace communication {
struct serial_global_policy {
- using id_type = uint32_t;
-
- std::vector<spike<id_type>> const
- static gather_spikes(const std::vector<spike<id_type>>& local_spikes) {
+ template <typename I, typename T>
+ static const std::vector<spike<I, T>>&
+ gather_spikes(const std::vector<spike<I, T>>& local_spikes) {
return local_spikes;
}
diff --git a/src/compartment.hpp b/src/compartment.hpp
index b92360a560ba5de30a036505d50d457b7b5a5a9c..da6bbcf57298593a60f7510822832060c7a89846 100644
--- a/src/compartment.hpp
+++ b/src/compartment.hpp
@@ -3,6 +3,8 @@
#include <iterator>
#include <utility>
+#include "common_types.hpp"
+
namespace nest {
namespace mc {
@@ -10,7 +12,7 @@ namespace mc {
/// The compartment is a conic frustrum
struct compartment {
using value_type = double;
- using size_type = int;
+ using size_type = cell_local_size_type;
using value_pair = std::pair<value_type, value_type>;
compartment() = delete;
@@ -103,8 +105,7 @@ class compartment_iterator :
};
class compartment_range {
- public:
-
+public:
using size_type = compartment_iterator::size_type;
using real_type = compartment_iterator::real_type;
diff --git a/src/connection.hpp b/src/connection.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b9e40fb5baeb83ad69c87c219f65d9f98acc4456
--- /dev/null
+++ b/src/connection.hpp
@@ -0,0 +1,68 @@
+#pragma once
+
+#include <cstdint>
+
+#include <common_types.hpp>
+#include <event_queue.hpp>
+#include <spike.hpp>
+
+namespace nest {
+namespace mc {
+
+template <typename Time>
+class connection {
+public:
+ using id_type = cell_member_type;
+ using time_type = Time;
+
+ connection(id_type src, id_type dest, float w, time_type d) :
+ source_(src),
+ destination_(dest),
+ weight_(w),
+ delay_(d)
+ {}
+
+ float weight() const { return weight_; }
+ float delay() const { return delay_; }
+
+ id_type source() const { return source_; }
+ id_type destination() const { return destination_; }
+
+ postsynaptic_spike_event<time_type> make_event(spike<id_type, time_type> s) {
+ return {destination_, s.time + delay_, weight_};
+ }
+
+private:
+ id_type source_;
+ id_type destination_;
+ float weight_;
+ time_type delay_;
+};
+
+// connections are sorted by source id
+// these operators make for easy interopability with STL algorithms
+
+template <typename T>
+static inline bool operator<(connection<T> lhs, connection<T> rhs) {
+ return lhs.source() < rhs.source();
+}
+
+template <typename T>
+static inline bool operator<(connection<T> lhs, typename connection<T>::id_type rhs) {
+ return lhs.source() < rhs;
+}
+
+template <typename T>
+static inline bool operator<(typename connection<T>::id_type lhs, connection<T> rhs) {
+ return lhs < rhs.source();
+}
+
+} // namespace mc
+} // namespace nest
+
+template <typename T>
+static inline std::ostream& operator<<(std::ostream& o, nest::mc::connection<T> const& con) {
+ return o << "con [" << con.source() << " -> " << con.destination()
+ << " : weight " << con.weight()
+ << ", delay " << con.delay() << "]";
+}
diff --git a/src/event_queue.hpp b/src/event_queue.hpp
index 7f3034f3f08bbb2235bf457c95e8eedfe8d0895d..31f6a5a10cf6b0c9ff9c0ff049909b801cd25b34 100644
--- a/src/event_queue.hpp
+++ b/src/event_queue.hpp
@@ -4,35 +4,49 @@
#include <ostream>
#include <queue>
+#include "common_types.hpp"
#include "util/optional.hpp"
namespace nest {
namespace mc {
+/* An event class Event must comply with the following conventions:
+ * Typedefs:
+ * time_type floating point type used to represent event times
+ * Member functions:
+ * time_type when() const return time value associated with event
+ */
+
+template <typename Time>
struct postsynaptic_spike_event {
- uint32_t target;
- float time;
+ using time_type = Time;
+
+ cell_member_type target;
+ time_type time;
float weight;
-};
-inline float event_time(const postsynaptic_spike_event &ev) { return ev.time; }
+ time_type when() const { return time; }
+};
+template <typename Time>
struct sample_event {
- uint32_t sampler_index;
- float time;
-};
+ using time_type = Time;
-inline float event_time(const sample_event &ev) { return ev.time; }
+ std::uint32_t sampler_index;
+ time_type time;
+
+ time_type when() const { return time; }
+};
/* Event objects must have a method event_time() which returns a value
* from a type with a total ordering with respect to <, >, etc.
- */
+ */
template <typename Event>
class event_queue {
public :
using value_type = Event;
- using time_type = decltype(event_time(std::declval<Event>()));
+ using time_type = typename Event::time_type;
// create
event_queue() {}
@@ -46,7 +60,7 @@ public :
}
// push thing
- void push(const value_type &e) {
+ void push(const value_type& e) {
queue_.push(e);
}
@@ -56,7 +70,7 @@ public :
// pop until
util::optional<value_type> pop_if_before(time_type t_until) {
- if (!queue_.empty() && event_time(queue_.top()) < t_until) {
+ if (!queue_.empty() && queue_.top().when() < t_until) {
auto ev = queue_.top();
queue_.pop();
return ev;
@@ -66,10 +80,15 @@ public :
}
}
+ // clear everything
+ void clear() {
+ queue_ = decltype(queue_){};
+ }
+
private:
struct event_greater {
- bool operator()(const Event &a, const Event &b) {
- return event_time(a) > event_time(b);
+ bool operator()(const Event& a, const Event& b) {
+ return a.when() > b.when();
}
};
@@ -83,8 +102,8 @@ private:
} // namespace nest
} // namespace mc
-inline
-std::ostream& operator<< (std::ostream& o, const nest::mc::postsynaptic_spike_event& e)
+template <typename T>
+inline std::ostream& operator<<(std::ostream& o, const nest::mc::postsynaptic_spike_event<T>& e)
{
return o << "event[" << e.target << "," << e.time << "," << e.weight << "]";
}
diff --git a/src/fvm_cell.hpp b/src/fvm_cell.hpp
index 932c0a4883747d1bd5828dd4be7df2c24aac532e..526fd22a391b477d5af3c88a4535b4721c83cbfe 100644
--- a/src/fvm_cell.hpp
+++ b/src/fvm_cell.hpp
@@ -13,15 +13,13 @@
#include <math.hpp>
#include <matrix.hpp>
#include <mechanism.hpp>
-#include <mechanism_interface.hpp>
+#include <mechanism_catalogue.hpp>
#include <segment.hpp>
#include <stimulus.hpp>
#include <util.hpp>
#include <profiling/profiler.hpp>
#include <vector/include/Vector.hpp>
-#include <mechanisms/expsyn.hpp>
-
namespace nest {
namespace mc {
@@ -96,7 +94,6 @@ public:
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_; }
@@ -116,8 +113,9 @@ public:
void advance(value_type dt);
/// pass an event to the appropriate synapse and call net_receive
- void apply_event(postsynaptic_spike_event e) {
- mechanisms_[synapse_index_]->net_receive(e.target, e.weight);
+ template <typename Time>
+ void apply_event(postsynaptic_spike_event<Time> e) {
+ mechanisms_[synapse_index_]->net_receive(e.target.index, e.weight);
}
mechanism_type& synapses() {
@@ -133,9 +131,19 @@ public:
/// returns voltage at a segment location
value_type voltage(segment_location loc) const;
+ /// flags if solution is physically realistic.
+ /// here we define physically realistic as the voltage being within reasonable bounds.
+ /// use a simple test of the voltage at the soma is reasonable, i.e. in the range
+ /// v_soma \in (-1000mv, 1000mv)
+ bool is_physical_solution() const {
+ auto v = voltage_[0];
+ return (v>-1000.) && (v<1000.);
+ }
+
/// returns current at a segment location
value_type current(segment_location loc) const;
+
value_type time() const { return t_; }
value_type probe(uint32_t i) const {
@@ -176,10 +184,7 @@ private:
/// the potential in mV in each CV
vector_type voltage_;
- /// synapses
- using synapse_type =
- mechanisms::ExpSyn::mechanism_ExpSyn<value_type, size_type>;
- std::size_t synapse_index_;
+ std::size_t synapse_index_; // synapses at the end of mechanisms_, from here
/// the set of mechanisms present in the cell
std::vector<mechanism_type> mechanisms_;
@@ -190,6 +195,9 @@ private:
std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
std::vector<std::pair<const vector_type fvm_cell::*, uint32_t>> probes_;
+
+ // mechanism factory
+ using mechanism_catalogue = nest::mc::mechanisms::catalogue<T, I>;
};
////////////////////////////////////////////////////////////////////////////////
@@ -287,10 +295,10 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// 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;
+ std::map<std::string, std::vector<unsigned>> mech_map;
- for(auto i=0; i<cell.num_segments(); ++i) {
- for(const auto& mech : cell.segment(i)->mechanisms()) {
+ for (unsigned 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
@@ -304,12 +312,12 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// 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);
+ 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) {
+ for (auto seg : mech.second) {
num_comp += segment_index_[seg+1] - segment_index_[seg];
}
@@ -317,7 +325,7 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// the mechanism
index_type compartment_index(num_comp);
auto pos = 0u;
- for(auto seg : mech.second) {
+ for (auto seg : mech.second) {
auto seg_size = segment_index_[seg+1] - segment_index_[seg];
std::iota(
compartment_index.data() + pos,
@@ -329,10 +337,29 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// instantiate the mechanism
mechanisms_.push_back(
- helper->new_mechanism(voltage_, current_, compartment_index)
+ mechanism_catalogue::make(mech.first, voltage_, current_, compartment_index)
+ //helper->new_mechanism(voltage_, current_, compartment_index)
);
}
+ synapse_index_ = mechanisms_.size();
+
+ std::map<std::string, std::vector<cell_lid_type>> syn_map;
+ for (const auto& syn : cell.synapses()) {
+ syn_map[syn.mechanism.name()].push_back(find_compartment_index(syn.location, graph));
+ }
+
+ for (const auto& syni : syn_map) {
+ const auto& mech_name = syni.first;
+ // auto& helper = nest::mc::mechanisms::get_mechanism_helper(mech_name);
+
+ index_type compartment_index(syni.second);
+ auto mech = mechanism_catalogue::make(mech_name, voltage_, current_, compartment_index);
+ mech->set_areas(cv_areas_);
+ mechanisms_.push_back(std::move(mech));
+ }
+
+
/////////////////////////////////////////////
// build the ion species
/////////////////////////////////////////////
@@ -347,7 +374,7 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
}
}
}
- std::vector<int> indexes(index_set.begin(), index_set.end());
+ std::vector<cell_lid_type> indexes(index_set.begin(), index_set.end());
// create the ion state
if(indexes.size()) {
@@ -389,24 +416,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
stimulii_.push_back( {idx, stim.clamp} );
}
- // add the synapses
- std::vector<size_type> synapse_indexes;
- synapse_indexes.reserve(cell.synapses().size());
- for(auto loc : cell.synapses()) {
- synapse_indexes.push_back(
- find_compartment_index(loc, graph)
- );
- }
-
- mechanisms_.push_back(
- mechanisms::make_mechanism<synapse_type>(
- voltage_, current_, index_view(synapse_indexes)
- )
- );
- synapse_index_ = mechanisms_.size()-1;
- // don't forget to give point processes access to cv_areas_
- mechanisms_[synapse_index_]->set_areas(cv_areas_);
-
// record probe locations by index into corresponding state vector
for (auto probe : cell.probes()) {
uint32_t comp = find_compartment_index(probe.location, graph);
@@ -471,7 +480,7 @@ void fvm_cell<T, I>::setup_matrix(T dt)
template <typename T, typename I>
int fvm_cell<T, I>::compartment_index(segment_location loc) const
{
- EXPECTS(loc.segment < segment_index_.size());
+ EXPECTS(unsigned(loc.segment) < segment_index_.size());
const auto seg = loc.segment;
@@ -508,15 +517,15 @@ void fvm_cell<T, I>::advance(T dt)
{
using memory::all;
- mc::util::profiler_enter("current");
+ PE("current");
current_(all) = 0.;
// update currents from ion channels
for(auto& m : mechanisms_) {
- mc::util::profiler_enter(m->name().c_str());
+ PE(m->name().c_str());
m->set_params(t_, dt);
m->nrn_current();
- mc::util::profiler_leave();
+ PL();
}
// add current contributions from stimulii
@@ -527,25 +536,25 @@ void fvm_cell<T, I>::advance(T dt)
// the factor of 100 scales the injected current to 10^2.nA
current_[loc] -= 100.*ie/cv_areas_[loc];
}
- mc::util::profiler_leave();
+ PL();
- mc::util::profiler_enter("matrix", "setup");
+ PE("matrix", "setup");
// solve the linear system
setup_matrix(dt);
- mc::util::profiler_leave(); mc::util::profiler_enter("solve");
+ PL(); PE("solve");
matrix_.solve();
- mc::util::profiler_leave();
+ PL();
voltage_(all) = matrix_.rhs();
- mc::util::profiler_leave();
+ PL();
- mc::util::profiler_enter("state");
+ PE("state");
// integrate state of gating variables etc.
for(auto& m : mechanisms_) {
- mc::util::profiler_enter(m->name().c_str());
+ PE(m->name().c_str());
m->nrn_state();
- mc::util::profiler_leave();
+ PL();
}
- mc::util::profiler_leave();
+ PL();
t_ += dt;
}
diff --git a/src/matrix.hpp b/src/matrix.hpp
index d0ddb816adc5cb9bfb6c3c88562e0a41835920ec..79b7218bc97fa0b16469ac50dbb5a91dc58c0db6 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -164,7 +164,7 @@ class matrix {
auto const ncells = num_cells();
// loop over submatrices
- for(auto m=0; m<ncells; ++m) {
+ for (size_type m=0; m<ncells; ++m) {
auto first = cell_index_[m];
auto last = cell_index_[m+1];
diff --git a/src/mechanism_catalogue.hpp b/src/mechanism_catalogue.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c9d3e6de2771a1b18fa05f104dd257b8c7fe6dd9
--- /dev/null
+++ b/src/mechanism_catalogue.hpp
@@ -0,0 +1,61 @@
+#pragma once
+
+#include <map>
+#include <stdexcept>
+#include <string>
+
+#include <mechanism.hpp>
+#include <mechanisms/hh.hpp>
+#include <mechanisms/pas.hpp>
+#include <mechanisms/expsyn.hpp>
+#include <mechanisms/exp2syn.hpp>
+
+namespace nest {
+namespace mc {
+namespace mechanisms {
+
+template <typename T, typename I>
+struct catalogue {
+ using view_type = typename mechanism<T, I>::view_type;
+ using index_view = typename mechanism<T, I>::index_view;
+
+ static mechanism_ptr<T, I> make(
+ const std::string& name,
+ view_type vec_v,
+ view_type vec_i,
+ index_view node_indices)
+ {
+ auto entry = mech_map.find(name);
+ if (entry==mech_map.end()) {
+ throw std::out_of_range("no such mechanism");
+ }
+
+ return entry->second(vec_v, vec_i, node_indices);
+ }
+
+ static bool has(const std::string& name) {
+ return mech_map.count(name)>0;
+ }
+
+private:
+ using maker_type = mechanism_ptr<T, I> (*)(view_type, view_type, index_view);
+ static const std::map<std::string, maker_type> mech_map;
+
+ template <template <typename, typename> class mech>
+ static mechanism_ptr<T, I> maker(view_type vec_v, view_type vec_i, index_view node_indices) {
+ return make_mechanism<mech<T, I>>(vec_v, vec_i, node_indices);
+ }
+};
+
+template <typename T, typename I>
+const std::map<std::string, typename catalogue<T, I>::maker_type> catalogue<T, I>::mech_map = {
+ { "pas", maker<pas::mechanism_pas> },
+ { "hh", maker<hh::mechanism_hh> },
+ { "expsyn", maker<expsyn::mechanism_expsyn> },
+ { "exp2syn", maker<exp2syn::mechanism_exp2syn> }
+};
+
+
+} // namespace mechanisms
+} // namespace mc
+} // namespace nest
diff --git a/src/mechanism_interface.cpp b/src/mechanism_interface.cpp
deleted file mode 100644
index f44d1f561d3efcc51b262b5bb32e0a95084bbcb5..0000000000000000000000000000000000000000
--- a/src/mechanism_interface.cpp
+++ /dev/null
@@ -1,45 +0,0 @@
-#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
index c59caaad6c655d6a69c5206b956b1cb42959cfdf..5e5360e6602761ae2923c6c2afec958279235dca 100644
--- a/src/mechanism_interface.hpp
+++ b/src/mechanism_interface.hpp
@@ -1,7 +1,6 @@
#pragma once
-#include <map>
-#include <string>
+// just for compatibility with current version of modparser...
#include "mechanism.hpp"
#include "parameter_list.hpp"
@@ -10,11 +9,6 @@ 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;
@@ -25,34 +19,10 @@ struct mechanism_helper {
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/model.hpp b/src/model.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..84dd84298a78418ce882b68a29f2211a6e50401b
--- /dev/null
+++ b/src/model.hpp
@@ -0,0 +1,209 @@
+#pragma once
+
+#include <cstdlib>
+#include <vector>
+
+#include <common_types.hpp>
+#include <cell.hpp>
+#include <cell_group.hpp>
+#include <fvm_cell.hpp>
+#include <recipe.hpp>
+#include <thread_private_spike_store.hpp>
+#include <communication/communicator.hpp>
+#include <communication/global_policy.hpp>
+#include <profiling/profiler.hpp>
+
+#include "trace_sampler.hpp"
+
+namespace nest {
+namespace mc {
+
+template <typename Cell>
+class model {
+public:
+ using cell_group_type = cell_group<Cell>;
+ using time_type = typename cell_group_type::time_type;
+ using value_type = typename cell_group_type::value_type;
+ using communicator_type = communication::communicator<time_type, communication::global_policy>;
+ using sampler_function = typename cell_group_type::sampler_function;
+
+ struct probe_record {
+ cell_member_type id;
+ probe_spec probe;
+ };
+
+ model(const recipe& rec, cell_gid_type cell_from, cell_gid_type cell_to):
+ cell_from_(cell_from),
+ cell_to_(cell_to),
+ communicator_(cell_from, cell_to)
+ {
+ // generate the cell groups in parallel, with one task per cell group
+ cell_groups_ = std::vector<cell_group_type>{cell_to_-cell_from_};
+ threading::parallel_vector<probe_record> probes;
+
+ threading::parallel_for::apply(cell_from_, cell_to_,
+ [&](cell_gid_type i) {
+ PE("setup", "cells");
+ auto cell = rec.get_cell(i);
+ auto idx = i-cell_from_;
+ cell_groups_[idx] = cell_group_type(i, cell);
+
+ cell_lid_type j = 0;
+ for (const auto& probe: cell.probes()) {
+ cell_member_type probe_id{i,j++};
+ probes.push_back({probe_id, probe});
+ }
+ PL(2);
+ });
+
+ // insert probes
+ probes_.assign(probes.begin(), probes.end());
+
+ // generate the network connections
+ for (cell_gid_type i=cell_from_; i<cell_to_; ++i) {
+ for (const auto& cc: rec.connections_on(i)) {
+ connection<time_type> conn{cc.source, cc.dest, cc.weight, cc.delay};
+ communicator_.add_connection(conn);
+ }
+ }
+ communicator_.construct();
+
+ // Allocate an empty queue buffer for each cell group
+ // These must be set initially to ensure that a queue is available for each
+ // cell group for the first time step.
+ current_events().resize(num_groups());
+ future_events().resize(num_groups());
+ }
+
+ void reset() {
+ t_ = 0.;
+ for (auto& group: cell_groups_) {
+ group.reset();
+ }
+ communicator_.reset();
+ }
+
+ time_type run(time_type tfinal, time_type dt) {
+ // Calculate the size of the largest possible time integration interval
+ // before communication of spikes is required.
+ // If spike exchange and cell update are serialized, this is the
+ // minimum delay of the network, however we use half this period
+ // to overlap communication and computation.
+ time_type t_interval = communicator_.min_delay()/2;
+
+ while (t_<tfinal) {
+ auto tuntil = std::min(t_+t_interval, tfinal);
+
+ event_queues_.exchange();
+ local_spikes_.exchange();
+
+ // empty the spike buffers for the current integration period.
+ // these buffers will store the new spikes generated in update_cells.
+ current_spikes().clear();
+
+ // task that updates cell state in parallel.
+ auto update_cells = [&] () {
+ threading::parallel_for::apply(
+ 0u, cell_groups_.size(),
+ [&](unsigned i) {
+ auto &group = cell_groups_[i];
+
+ PE("stepping","events");
+ group.enqueue_events(current_events()[i]);
+ PL();
+
+ group.advance(tuntil, dt);
+
+ PE("events");
+ current_spikes().insert(group.spikes());
+ group.clear_spikes();
+ PL(2);
+ });
+ };
+
+ // task that performs spike exchange with the spikes generated in
+ // the previous integration period, generating the postsynaptic
+ // events that must be delivered at the start of the next
+ // integration period at the latest.
+ auto exchange = [&] () {
+ PE("stepping", "exchange");
+ auto local_spikes = previous_spikes().gather();
+ future_events() = communicator_.exchange(local_spikes);
+ PL(2);
+ };
+
+ // run the tasks, overlapping if the threading model and number of
+ // available threads permits it.
+ threading::task_group g;
+ g.run(exchange);
+ g.run(update_cells);
+ g.wait();
+
+ t_ = tuntil;
+ }
+ return t_;
+ }
+
+ // only thread safe if called outside the run() method
+ void add_artificial_spike(cell_member_type source) {
+ add_artificial_spike(source, t_);
+ }
+
+ // only thread safe if called outside the run() method
+ void add_artificial_spike(cell_member_type source, time_type tspike) {
+ current_spikes().get().push_back({source, tspike});
+ }
+
+ void attach_sampler(cell_member_type probe_id, sampler_function f, time_type tfrom = 0) {
+ // TODO: translate probe_id.gid to appropriate group, but for now 1-1.
+ if (probe_id.gid<cell_from_ || probe_id.gid>=cell_to_) {
+ return;
+ }
+ cell_groups_[probe_id.gid-cell_from_].add_sampler(probe_id, f, tfrom);
+ }
+
+ const std::vector<probe_record>& probes() const { return probes_; }
+
+ std::size_t num_spikes() const { return communicator_.num_spikes(); }
+ std::size_t num_groups() const { return cell_groups_.size(); }
+
+private:
+ cell_gid_type cell_from_;
+ cell_gid_type cell_to_;
+ time_type t_ = 0.;
+ std::vector<cell_group_type> cell_groups_;
+ communicator_type communicator_;
+ std::vector<probe_record> probes_;
+ using spike_type = typename communicator_type::spike_type;
+
+ using event_queue_type = typename communicator_type::event_queue;
+ util::double_buffer< std::vector<event_queue_type> > event_queues_;
+
+ using local_spike_store_type = thread_private_spike_store<time_type>;
+ util::double_buffer< local_spike_store_type > local_spikes_;
+
+ // Convenience functions that map the spike buffers and event queues onto
+ // the appropriate integration interval.
+ //
+ // To overlap communication and computation, integration intervals of
+ // size Delta/2 are used, where Delta is the minimum delay in the global
+ // system.
+ // From the frame of reference of the current integration period we
+ // define three intervals: previous, current and future
+ // Then we define the following :
+ // current_spikes : spikes generated in the current interval
+ // previous_spikes: spikes generated in the preceding interval
+ // current_events : events to be delivered at the start of
+ // the current interval
+ // future_events : events to be delivered at the start of
+ // the next interval
+
+ local_spike_store_type& current_spikes() { return local_spikes_.get(); }
+ local_spike_store_type& previous_spikes() { return local_spikes_.other(); }
+
+ std::vector<event_queue_type>& current_events() { return event_queues_.get(); }
+ std::vector<event_queue_type>& future_events() { return event_queues_.other(); }
+};
+
+} // namespace mc
+} // namespace nest
diff --git a/src/point.hpp b/src/point.hpp
index 3b048bb3db7fb580443e019c98292821f715ffb2..256e518d06f6d9817a65ba295eb9b560d9327000 100644
--- a/src/point.hpp
+++ b/src/point.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <cmath>
#include <limits>
#include <ostream>
@@ -13,77 +14,58 @@ struct point {
value_type y;
value_type z;
- constexpr point(T a, T b, T c)
- : x{a}, y{b}, z{c}
+ 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 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
- {
+ 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
-) {
+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
-) {
+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
-) {
+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);
+T norm(point<T> const& p) {
+ return std::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
-) {
+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);
+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)
-{
+bool operator!=(const point<T> &rhs, const point<T> &lhs) {
return !(rhs == lhs);
}
@@ -91,7 +73,6 @@ bool operator!=(const point<T> &rhs, const point<T> &lhs)
} // namespace nest
template <typename T>
-std::ostream& operator << (std::ostream& o, nest::mc::point<T> const& p)
-{
+std::ostream& operator << (std::ostream& o, nest::mc::point<T> const& p) {
return o << "[" << p.x << ", " << p.y << ", " << p.z << "]";
}
diff --git a/src/profiling/profiler.cpp b/src/profiling/profiler.cpp
index b5d2ad0d8d7bd9efe7fd8e55e5f6720ab051b3f1..1e30449ab4b225385570292bea527fbcc6981c3d 100644
--- a/src/profiling/profiler.cpp
+++ b/src/profiling/profiler.cpp
@@ -66,7 +66,7 @@ void profiler_node::print_sub(
if (print_children) {
auto other = 0.;
- for (auto &n : children) {
+ for (auto& n : children) {
if (n.value<threshold || n.name=="other") {
other += n.value;
}
@@ -162,7 +162,7 @@ double region_type::subregion_contributions() const {
profiler_node region_type::populate_performance_tree() const {
profiler_node tree(total(), name());
- for (auto &it : subregions_) {
+ for (auto& it : subregions_) {
tree.children.push_back(it.second->populate_performance_tree());
}
diff --git a/src/profiling/profiler.hpp b/src/profiling/profiler.hpp
index 7230697baab21f914fa3b92f4ffbed4382b062b6..4332925a5b3ee46c8c07164a80c4acf5abd72989 100644
--- a/src/profiling/profiler.hpp
+++ b/src/profiling/profiler.hpp
@@ -226,6 +226,7 @@ void profiler_enter(const char* n, Args... args) {
/// move up one level in the profiler
void profiler_leave();
+
/// move up multiple profiler levels in one call
void profiler_leave(int nlevels);
@@ -238,3 +239,8 @@ void profiler_output(double threshold);
} // namespace util
} // namespace mc
} // namespace nest
+
+// define some helper macros to make instrumentation of the source code with calls
+// to the profiler a little less visually distracting
+#define PE nest::mc::util::profiler_enter
+#define PL nest::mc::util::profiler_leave
diff --git a/src/recipe.hpp b/src/recipe.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b451a1a1cffbd2139707e8f8d8a60fa6e2f4e037
--- /dev/null
+++ b/src/recipe.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+#include <cstddef>
+#include <memory>
+#include <stdexcept>
+
+namespace nest {
+namespace mc {
+
+struct cell_count_info {
+ cell_size_type num_sources;
+ cell_size_type num_targets;
+ cell_size_type num_probes;
+};
+
+class invalid_recipe_error: public std::runtime_error {
+public:
+ invalid_recipe_error(std::string whatstr): std::runtime_error(std::move(whatstr)) {}
+};
+
+/* Recipe descriptions are cell-oriented: in order that the building
+ * phase can be done distributedly and in order that the recipe
+ * description can be built indepdently of any runtime execution
+ * environment, connection end-points are represented by pairs
+ * (cell index, source/target index on cell).
+ */
+
+using cell_connection_endpoint = cell_member_type;
+
+// Note: `cell_connection` and `connection` have essentially the same data
+// and represent the same thing conceptually. `cell_connection` objects
+// are notionally described in terms of external cell identifiers instead
+// of internal gids, but we are not making the distinction between the
+// two in the current code. These two types could well be merged.
+
+struct cell_connection {
+ cell_connection_endpoint source;
+ cell_connection_endpoint dest;
+
+ float weight;
+ float delay;
+};
+
+class recipe {
+public:
+ virtual cell_size_type num_cells() const =0;
+
+ virtual cell get_cell(cell_gid_type) const =0;
+ virtual cell_count_info get_cell_count_info(cell_gid_type) const =0;
+ virtual std::vector<cell_connection> connections_on(cell_gid_type) const =0;
+};
+
+} // namespace mc
+} // namespace nest
diff --git a/src/segment.hpp b/src/segment.hpp
index c6c49b26a946d36fdc969d9473fa07105d30d4af..248306a6c3a976d123b9d412b93c09156305086d 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -4,6 +4,7 @@
#include <vector>
#include "algorithms.hpp"
+#include "common_types.hpp"
#include "compartment.hpp"
#include "math.hpp"
#include "parameter_list.hpp"
@@ -36,6 +37,7 @@ class segment {
public:
using value_type = double;
+ using size_type = cell_local_size_type;
using point_type = point<value_type>;
segmentKind kind() const {
@@ -57,8 +59,8 @@ class segment {
return kind_==segmentKind::axon;
}
- virtual int num_compartments() const = 0;
- virtual void set_compartments(int) = 0;
+ virtual size_type num_compartments() const = 0;
+ virtual void set_compartments(size_type) = 0;
virtual value_type volume() const = 0;
virtual value_type area() const = 0;
@@ -153,13 +155,12 @@ class segment {
std::vector<parameter_list> mechanisms_;
};
-class placeholder_segment : public segment
-{
- public:
-
+class placeholder_segment : public segment {
+public:
using base = segment;
using base::kind_;
using base::value_type;
+ using base::size_type;
placeholder_segment()
{
@@ -181,23 +182,22 @@ class placeholder_segment : public segment
return true;
}
- int num_compartments() const override
+ size_type num_compartments() const override
{
return 0;
}
- virtual void set_compartments(int) override
+ virtual void set_compartments(size_type) override
{ }
};
-class soma_segment : public segment
-{
- public :
-
+class soma_segment : public segment {
+public:
using base = segment;
using base::kind_;
using base::value_type;
using base::point_type;
+ using base::size_type;
soma_segment() = delete;
@@ -208,7 +208,7 @@ class soma_segment : public segment
mechanisms_.push_back(membrane_parameters());
}
- soma_segment(value_type r, point_type const &c)
+ soma_segment(value_type r, point_type const& c)
: soma_segment(r)
{
center_ = c;
@@ -245,12 +245,12 @@ class soma_segment : public segment
}
/// soma has one and one only compartments
- int num_compartments() const override
+ size_type num_compartments() const override
{
return 1;
}
- void set_compartments(int n) override
+ void set_compartments(size_type n) override
{ }
private :
@@ -261,10 +261,8 @@ class soma_segment : public segment
point_type center_;
};
-class cable_segment : public segment
-{
- public :
-
+class cable_segment : public segment {
+public:
using base = segment;
using base::kind_;
using base::value_type;
@@ -332,7 +330,7 @@ class cable_segment : public segment
value_type volume() const override
{
auto sum = value_type{0};
- for(auto i=0; i<num_sub_segments(); ++i) {
+ for (auto i=0u; i<num_sub_segments(); ++i) {
sum += math::volume_frustrum(lengths_[i], radii_[i], radii_[i+1]);
}
return sum;
@@ -341,7 +339,7 @@ class cable_segment : public segment
value_type area() const override
{
auto sum = value_type{0};
- for(auto i=0; i<num_sub_segments(); ++i) {
+ for (auto i=0u; i<num_sub_segments(); ++i) {
sum += math::area_frustrum(lengths_[i], radii_[i], radii_[i+1]);
}
return sum;
@@ -358,7 +356,7 @@ class cable_segment : public segment
}
// the number sub-segments that define the cable segment
- int num_sub_segments() const
+ size_type num_sub_segments() const
{
return radii_.size()-1;
}
@@ -378,12 +376,12 @@ class cable_segment : public segment
return this;
}
- int num_compartments() const override
+ size_type num_compartments() const override
{
return num_compartments_;
}
- void set_compartments(int n) override
+ void set_compartments(size_type n) override
{
if(n<1) {
throw std::out_of_range(
@@ -406,8 +404,8 @@ class cable_segment : public segment
// 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) {
+ size_type i = 0;
+ for (i = 0; i<num_sub_segments(); ++i) {
if(sum+lengths_[i]>pos) {
break;
}
@@ -425,19 +423,18 @@ class cable_segment : public segment
return {num_compartments(), radii_.front(), radii_.back(), length()};
}
- private :
-
+private:
void update_lengths()
{
- if(locations_.size()) {
+ if (locations_.size()) {
lengths_.resize(num_sub_segments());
- for(auto i=0; i<num_sub_segments(); ++i) {
+ for (size_type i=0; i<num_sub_segments(); ++i) {
lengths_[i] = norm(locations_[i] - locations_[i+1]);
}
}
}
- int num_compartments_ = 1;
+ size_type num_compartments_ = 1;
std::vector<value_type> lengths_;
std::vector<value_type> radii_;
std::vector<point_type> locations_;
diff --git a/src/communication/spike.hpp b/src/spike.hpp
similarity index 53%
rename from src/communication/spike.hpp
rename to src/spike.hpp
index 03b5ed75def7a2fb68bd491ae7f1208d4b230b8a..d3ea02551456d2408712886dcff8125b314630e9 100644
--- a/src/communication/spike.hpp
+++ b/src/spike.hpp
@@ -1,48 +1,39 @@
#pragma once
-#include <type_traits>
#include <ostream>
+#include <type_traits>
namespace nest {
namespace mc {
-namespace communication {
-template <
- typename I,
- typename = typename std::enable_if<std::is_integral<I>::value>
->
+template <typename I, typename Time>
struct spike {
using id_type = I;
- id_type source = 0;
- float time = -1.;
+ using time_type = Time;
+
+ id_type source = id_type{};
+ time_type time = -1.;
spike() = default;
- spike(id_type s, float t) :
+ spike(id_type s, time_type t) :
source(s), time(t)
{}
};
} // namespace mc
} // namespace nest
-} // namespace communication
/// custom stream operator for printing nest::mc::spike<> values
-template <typename I>
-std::ostream& operator<<(
- std::ostream& o,
- nest::mc::communication::spike<I> s)
-{
+template <typename I, typename T>
+std::ostream& operator<<(std::ostream& o, nest::mc::spike<I, T> s) {
return o << "spike[t " << s.time << ", src " << s.source << "]";
}
/// less than comparison operator for nest::mc::spike<> values
/// spikes are ordered by spike time, for use in sorting and queueing
-template <typename I>
-bool operator<(
- nest::mc::communication::spike<I> lhs,
- nest::mc::communication::spike<I> rhs)
-{
+template <typename I, typename T>
+bool operator<(nest::mc::spike<I, T> lhs, nest::mc::spike<I, T> rhs) {
return lhs.time < rhs.time;
}
diff --git a/src/communication/spike_source.hpp b/src/spike_source.hpp
similarity index 84%
rename from src/communication/spike_source.hpp
rename to src/spike_source.hpp
index 95099cb547de24e76386c2253960b8fffe4bf569..0537b6800a19534d23aedd7de2edae3d022b53e6 100644
--- a/src/communication/spike_source.hpp
+++ b/src/spike_source.hpp
@@ -13,13 +13,11 @@ class spike_detector
public:
using cell_type = Cell;
- spike_detector( const cell_type& cell, segment_location loc, double thresh, float t_init) :
+ spike_detector(const cell_type& cell, segment_location loc, double thresh, float t_init) :
location_(loc),
- threshold_(thresh),
- previous_t_(t_init)
+ threshold_(thresh)
{
- previous_v_ = cell.voltage(location_);
- is_spiking_ = previous_v_ >= thresh ? true : false;
+ reset(cell, t_init);
}
util::optional<float> test(const cell_type& cell, float t) {
@@ -58,6 +56,12 @@ public:
float v() const { return previous_v_; }
+ void reset(const cell_type& cell, float t_init) {
+ previous_t_ = t_init;
+ previous_v_ = cell.voltage(location_);
+ is_spiking_ = previous_v_ >= threshold_;
+ }
+
private:
// parameters/data
diff --git a/src/thread_private_spike_store.hpp b/src/thread_private_spike_store.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1e14ab104619a42b8aa9f014172bb5716576017c
--- /dev/null
+++ b/src/thread_private_spike_store.hpp
@@ -0,0 +1,75 @@
+#pragma once
+
+#include <vector>
+
+#include <common_types.hpp>
+#include <spike.hpp>
+#include <threading/threading.hpp>
+
+namespace nest {
+namespace mc {
+
+/// Handles the complexity of managing thread private buffers of spikes.
+/// Internally stores one thread private buffer of spikes for each hardware thread.
+/// This can be accessed directly using the get() method, which returns a reference to
+/// The thread private buffer of the calling thread.
+/// The insert() and gather() methods add a vector of spikes to the buffer,
+/// and collate all of the buffers into a single vector respectively.
+template <typename Time>
+class thread_private_spike_store {
+public :
+ using id_type = cell_gid_type;
+ using time_type = Time;
+ using spike_type = spike<cell_member_type, time_type>;
+
+ /// Collate all of the individual buffers into a single vector of spikes.
+ /// Does not modify the buffer contents.
+ std::vector<spike_type> gather() const {
+ std::vector<spike_type> spikes;
+ unsigned num_spikes = 0u;
+ for (auto& b : buffers_) {
+ num_spikes += b.size();
+ }
+ spikes.reserve(num_spikes);
+
+ for (auto& b : buffers_) {
+ spikes.insert(spikes.begin(), b.begin(), b.end());
+ }
+
+ return spikes;
+ }
+
+ /// Return a reference to the thread private buffer of the calling thread
+ std::vector<spike_type>& get() {
+ return buffers_.local();
+ }
+
+ /// Return a reference to the thread private buffer of the calling thread
+ const std::vector<spike_type>& get() const {
+ return buffers_.local();
+ }
+
+ /// Clear all of the thread private buffers
+ void clear() {
+ for (auto& b : buffers_) {
+ b.clear();
+ }
+ }
+
+ /// Append the passed spikes to the end of the thread private buffer of the
+ /// calling thread
+ void insert(const std::vector<spike_type>& spikes) {
+ auto& buff = get();
+ buff.insert(buff.end(), spikes.begin(), spikes.end());
+ }
+
+private :
+ /// thread private storage for accumulating spikes
+ using local_spike_store_type =
+ threading::enumerable_thread_specific<std::vector<spike_type>>;
+
+ local_spike_store_type buffers_;
+};
+
+} // namespace mc
+} // namespace nest
diff --git a/src/threading/serial.hpp b/src/threading/serial.hpp
index eef783949049c77f3c06bfa059b02c06d96b3670..c18d4800e3c457e140f0fe6d2ca671aa6ccac241 100644
--- a/src/threading/serial.hpp
+++ b/src/threading/serial.hpp
@@ -1,12 +1,13 @@
#pragma once
#if !defined(WITH_SERIAL)
- #error this header can only be loaded if WITH_SERIAL is set
+ #error "this header can only be loaded if WITH_SERIAL is set"
#endif
#include <array>
#include <chrono>
#include <string>
+#include <vector>
namespace nest {
namespace mc {
@@ -18,6 +19,8 @@ namespace threading {
template <typename T>
class enumerable_thread_specific {
std::array<T, 1> data;
+ using iterator_type = typename std::array<T, 1>::iterator;
+ using const_iterator_type = typename std::array<T, 1>::const_iterator;
public :
@@ -36,11 +39,14 @@ class enumerable_thread_specific {
auto size() -> decltype(data.size()) const { return data.size(); }
- auto begin() -> decltype(data.begin()) { return data.begin(); }
- auto end() -> decltype(data.end()) { return data.end(); }
+ iterator_type begin() { return data.begin(); }
+ iterator_type end() { return data.end(); }
- auto cbegin() -> decltype(data.cbegin()) const { return data.cbegin(); }
- auto cend() -> decltype(data.cend()) const { return data.cend(); }
+ const_iterator_type begin() const { return data.begin(); }
+ const_iterator_type end() const { return data.end(); }
+
+ const_iterator_type cbegin() const { return data.cbegin(); }
+ const_iterator_type cend() const { return data.cend(); }
};
@@ -56,6 +62,10 @@ struct parallel_for {
}
};
+template <typename T>
+using parallel_vector = std::vector<T>;
+
+
inline std::string description() {
return "serial";
}
@@ -76,6 +86,36 @@ struct timer {
}
};
+constexpr bool multithreaded() { return false; }
+
+/// Proxy for tbb task group.
+/// The tbb version launches tasks asynchronously, returning control to the
+/// caller. The serial version implemented here simply runs the task, before
+/// returning control, effectively serializing all asynchronous calls.
+class task_group {
+public:
+ task_group() = default;
+
+ template<typename Func>
+ void run(const Func& f) {
+ f();
+ }
+
+ template<typename Func>
+ void run_and_wait(const Func& f) {
+ f();
+ }
+
+ void wait()
+ {}
+
+ bool is_canceling() {
+ return false;
+ }
+
+ void cancel()
+ {}
+};
} // threading
} // mc
diff --git a/src/threading/tbb.hpp b/src/threading/tbb.hpp
index ed82e26ec87617f7453f4548efe3eec3f858507f..853ceefde731cb1cd4db323541132c935011c73a 100644
--- a/src/threading/tbb.hpp
+++ b/src/threading/tbb.hpp
@@ -44,6 +44,13 @@ struct timer {
}
};
+constexpr bool multithreaded() { return true; }
+
+template <typename T>
+using parallel_vector = tbb::concurrent_vector<T>;
+
+using task_group = tbb::task_group;
+
} // threading
} // mc
} // nest
diff --git a/src/tree.hpp b/src/tree.hpp
index 98ff9aa43656f9d68b879297a6ffd5399348b6f0..cd59537411e90c7b5efa696ebe8d5baed9c5ea1a 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -12,16 +12,18 @@
namespace nest {
namespace mc {
+template <typename Int, typename Size = std::size_t>
class tree {
using range = memory::Range;
- public :
-
- using int_type = int;
+public:
+ using int_type = Int;
+ using size_type = Size;
using index_type = memory::HostVector<int_type>;
- using view_type = index_type::view_type;
- using const_view_type = index_type::const_view_type;
+ using view_type = typename index_type::view_type;
+ using const_view_type = typename index_type::const_view_type;
+ static constexpr int_type no_parent = (int_type)-1;
tree() = default;
@@ -67,12 +69,12 @@ class tree {
init(new_parent_index.size());
parents_(memory::all) = new_parent_index;
- parents_[0] = -1;
+ parents_[0] = no_parent;
child_index_(memory::all) =
algorithms::make_index(algorithms::child_count(parents_));
- std::vector<int> pos(parents_.size(), 0);
+ std::vector<int_type> pos(parents_.size(), 0);
for (auto i = 1u; i < parents_.size(); ++i) {
auto p = parents_[i];
children_[child_index_[p] + pos[p]] = i;
@@ -80,16 +82,18 @@ class tree {
}
}
- size_t num_children() const {
- return children_.size();
+ size_type num_children() const {
+ return static_cast<size_type>(children_.size());
}
- size_t num_children(size_t b) const {
+
+ size_type num_children(size_t b) const {
return child_index_[b+1] - child_index_[b];
}
- size_t num_nodes() const {
+
+ size_type num_nodes() const {
// 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();
+ auto sz = static_cast<size_type>(child_index_.size());
return sz ? sz - 1 : 0;
}
@@ -104,7 +108,7 @@ class tree {
}
/// return the list of all children of branch b
- const_view_type children(size_t b) const {
+ const_view_type children(size_type b) const {
return children_(child_index_[b], child_index_[b+1]);
}
@@ -122,7 +126,7 @@ class tree {
}
/// memory used to store tree (in bytes)
- size_t memory() const {
+ std::size_t memory() const {
return sizeof(int_type)*data_.size() + sizeof(tree);
}
@@ -164,17 +168,16 @@ class tree {
return p;
}
- private :
-
- void init(int nnode) {
- auto nchild = nnode -1;
+private:
+ void init(size_type nnode) {
+ auto nchild = nnode - 1;
data_ = index_type(nchild + (nnode + 1) + nnode);
set_ranges(nnode);
}
- void set_ranges(int nnode) {
- if(nnode) {
+ void set_ranges(size_type nnode) {
+ if (nnode) {
auto nchild = nnode - 1;
// data_ is partitioned as follows:
// data_ = [children_[nchild], child_index_[nnode+1], parents_[nnode]]
@@ -215,17 +218,17 @@ class tree {
/// new_node
/// p : permutation vector, p[i] is the new index of node i in the old
/// tree
- int add_children(
- int new_node,
- int old_node,
- int parent_node,
+ int_type add_children(
+ int_type new_node,
+ int_type old_node,
+ int_type parent_node,
view_type p,
tree const& old_tree
)
{
- // check for the senitel that indicates that the old root has
+ // check for the sentinel that indicates that the old root has
// been processed
- if(old_node==-1) {
+ if (old_node==no_parent) {
return new_node;
}
@@ -237,7 +240,7 @@ class tree {
auto this_node = new_node;
auto pos = child_index_[this_node];
- auto add_parent_as_child = parent_node>=0 && old_node>0;
+ auto add_parent_as_child = parent_node!=no_parent && old_node>0;
//
// STEP 1 : add the child indexes for this_node
//
@@ -259,12 +262,12 @@ class tree {
// STEP 2 : recursively add each child's children
//
new_node++;
- for(auto b : old_children) {
- if(b != parent_node) {
- new_node = add_children(new_node, b, -1, p, old_tree);
+ for (auto b : old_children) {
+ if (b != parent_node) {
+ new_node = add_children(new_node, b, no_parent, p, old_tree);
}
}
- if(add_parent_as_child) {
+ if (add_parent_as_child) {
new_node =
add_children(
new_node, old_tree.parent(old_node), old_node, p, old_tree
@@ -286,38 +289,38 @@ class tree {
view_type parents_ = data_(0, 0);
};
-template <typename C>
-std::vector<int> make_parent_index(tree const& t, C const& counts)
+template <typename IntT, typename SizeT, typename C>
+std::vector<IntT> make_parent_index(tree<IntT, SizeT> const& t, C const& counts)
{
using range = memory::Range;
+ using int_type = typename tree<IntT, SizeT>::int_type;
+ constexpr auto no_parent = tree<IntT, SizeT>::no_parent;
- if( !algorithms::is_positive(counts)
- || counts.size() != t.num_nodes() )
- {
+ 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())) {
+ std::vector<int_type> parent_index(num_compartments);
+ int_type pos = 0;
+ for (int_type 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;
+ parent = parent!=no_parent ? 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) {
+ 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]) {
+ while (pos<index[i+1]) {
parent_index[pos] = pos-1;
pos++;
}
diff --git a/src/util.hpp b/src/util.hpp
index 29a6b28a9ec09cce13fcced0db20ae30aa3d803f..243987e628685ee4d4414394b054bebc6e7646a0 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -18,7 +18,7 @@ using memory::util::cyan;
template <typename T>
std::ostream&
-operator << (std::ostream &o, std::vector<T>const& v)
+operator << (std::ostream& o, std::vector<T>const& v)
{
o << "[";
for(auto const& i: v) {
@@ -29,7 +29,7 @@ operator << (std::ostream &o, std::vector<T>const& v)
}
template <typename T>
-std::ostream& print(std::ostream &o, std::vector<T>const& v)
+std::ostream& print(std::ostream& o, std::vector<T>const& v)
{
o << "[";
for(auto const& i: v) {
diff --git a/src/util/debug.cpp b/src/util/debug.cpp
index 2dd83f562d3ea15477d7507026a370a87428a36b..51c646bbcb38b1722e40766f4886a957094bd1bb 100644
--- a/src/util/debug.cpp
+++ b/src/util/debug.cpp
@@ -1,10 +1,21 @@
+#include <chrono>
#include <cstdlib>
+#include <cstring>
+#include <iomanip>
#include <iostream>
+#include <mutex>
#include "util/debug.hpp"
+#include "util/ioutil.hpp"
-bool nest::mc::util::failed_assertion(const char *assertion, const char *file,
- int line, const char *func)
+namespace nest {
+namespace mc {
+namespace util {
+
+std::mutex global_debug_cerr_mutex;
+
+bool failed_assertion(const char* assertion, const char* file,
+ int line, const char* func)
{
// Explicit flush, as we can't assume default buffering semantics on stderr/cerr,
// and abort() might not flush streams.
@@ -14,3 +25,29 @@ bool nest::mc::util::failed_assertion(const char *assertion, const char *file,
std::abort();
return false;
}
+
+std::ostream& debug_emit_trace_leader(std::ostream& out, const char* file,
+ int line, const char* varlist)
+{
+ iosfmt_guard guard(out);
+
+ const char* leaf = std::strrchr(file, '/');
+ out << (leaf?leaf+1:file) << ':' << line << " ";
+
+ using namespace std::chrono;
+ auto tstamp = system_clock::now().time_since_epoch();
+ auto tstamp_usec = duration_cast<microseconds>(tstamp).count();
+
+ out << std::right << '[';
+ out << std::setw(11) << std::setfill('0') << (tstamp_usec/1000000) << '.';
+ out << std::setw(6) << std::setfill('0') << (tstamp_usec%1000000) << ']';
+
+ if (varlist && *varlist) {
+ out << ' ' << varlist << ": ";
+ }
+ return out;
+}
+
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/src/util/debug.hpp b/src/util/debug.hpp
index 148b3c44db9b03cd396c10db044e6924b99d7104..f258b37123755c7ec7a7324641249ed35290e8bb 100644
--- a/src/util/debug.hpp
+++ b/src/util/debug.hpp
@@ -1,30 +1,72 @@
#pragma once
+#include <iostream>
+#include <sstream>
+#include <mutex>
+
+#include "threading/threading.hpp"
+
namespace nest {
namespace mc {
namespace util {
-bool failed_assertion(const char *assertion, const char *file, int line, const char *func);
+bool failed_assertion(const char* assertion, const char* file, int line, const char* func);
+std::ostream& debug_emit_trace_leader(std::ostream& out, const char* file, int line, const char* varlist);
-}
+inline void debug_emit(std::ostream& out) {
+ out << "\n";
}
+
+template <typename Head, typename... Tail>
+void debug_emit(std::ostream& out, const Head& head, const Tail&... tail) {
+ out << head;
+ if (sizeof...(tail)) {
+ out << ", ";
+ }
+ debug_emit(out, tail...);
}
+extern std::mutex global_debug_cerr_mutex;
-#ifdef WITH_ASSERTIONS
+template <typename... Args>
+void debug_emit_trace(const char* file, int line, const char* varlist, const Args&... args) {
+ if (nest::mc::threading::multithreaded()) {
+ std::stringstream buffer;
-#ifdef __GNUC__
-#define DEBUG_FUNCTION_NAME __PRETTY_FUNCTION__
-#else
-#define DEBUG_FUNCTION_NAME __func__
-#endif
+ debug_emit_trace_leader(buffer, file, line, varlist);
+ debug_emit(buffer, args...);
-#define EXPECTS(condition) \
-(void)((condition) || \
- nest::mc::util::failed_assertion(#condition, __FILE__, __LINE__, DEBUG_FUNCTION_NAME))
+ std::lock_guard<std::mutex> guard(global_debug_cerr_mutex);
+ std::cerr << buffer.rdbuf();
+ std::cerr.flush();
+ }
+ else {
+ debug_emit_trace_leader(std::cerr, file, line, varlist);
+ debug_emit(std::cerr, args...);
+ std::cerr.flush();
+ }
+}
+
+} // namespace util
+} // namespace mc
+} // namespace nest
+#ifdef WITH_TRACE
+ #define TRACE(vars...) nest::mc::util::debug_emit_trace(__FILE__, __LINE__, #vars, ##vars)
#else
+ #define TRACE(...)
+#endif
-#define EXPECTS(condition)
+#ifdef WITH_ASSERTIONS
+ #ifdef __GNUC__
+ #define DEBUG_FUNCTION_NAME __PRETTY_FUNCTION__
+ #else
+ #define DEBUG_FUNCTION_NAME __func__
+ #endif
-#endif // def WITH_ASSERTIONS
+ #define EXPECTS(condition) \
+ (void)((condition) || \
+ nest::mc::util::failed_assertion(#condition, __FILE__, __LINE__, DEBUG_FUNCTION_NAME))
+#else
+ #define EXPECTS(condition)
+#endif // def WITH_ASSERTIONS
diff --git a/src/util/double_buffer.hpp b/src/util/double_buffer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..31cffa1babef2ced0b843b4d3e375678612528b2
--- /dev/null
+++ b/src/util/double_buffer.hpp
@@ -0,0 +1,66 @@
+#pragma once
+
+#include <array>
+#include <atomic>
+
+#include <util/debug.hpp>
+
+namespace nest {
+namespace mc {
+namespace util {
+
+/// double buffer with thread safe exchange/flip operation.
+template <typename T>
+class double_buffer {
+private:
+ std::atomic<int> index_;
+ std::array<T, 2> buffers_;
+
+ int other_index() {
+ return index_ ? 0 : 1;
+ }
+
+public:
+ using value_type = T;
+
+ double_buffer() :
+ index_(0)
+ {}
+
+ /// remove the copy and move constructors which won't work with std::atomic
+ double_buffer(double_buffer&&) = delete;
+ double_buffer(const double_buffer&) = delete;
+ double_buffer& operator=(const double_buffer&) = delete;
+ double_buffer& operator=(double_buffer&&) = delete;
+
+ /// flip the buffers in a thread safe manner
+ /// n calls to exchange will always result in n flips
+ void exchange() {
+ // use operator^= which is overloaded by std::atomic<>
+ index_ ^= 1;
+ }
+
+ /// get the current/front buffer
+ value_type& get() {
+ return buffers_[index_];
+ }
+
+ /// get the current/front buffer
+ const value_type& get() const {
+ return buffers_[index_];
+ }
+
+ /// get the back buffer
+ value_type& other() {
+ return buffers_[other_index()];
+ }
+
+ /// get the back buffer
+ const value_type& other() const {
+ return buffers_[other_index()];
+ }
+};
+
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/src/util/ioutil.hpp b/src/util/ioutil.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2feab394df71d730bde904d041ebbeeb85bf9a44
--- /dev/null
+++ b/src/util/ioutil.hpp
@@ -0,0 +1,99 @@
+#pragma once
+
+#include <iostream>
+
+namespace nest {
+namespace mc {
+namespace util {
+
+class iosfmt_guard {
+public:
+ explicit iosfmt_guard(std::ios& stream) :
+ save_(nullptr), stream_(stream)
+ {
+ save_.copyfmt(stream_);
+ }
+
+ ~iosfmt_guard() {
+ stream_.copyfmt(save_);
+ }
+
+private:
+ std::ios save_;
+ std::ios& stream_;
+};
+
+
+template <typename charT, typename traitsT = std::char_traits<charT> >
+class basic_null_streambuf: public std::basic_streambuf<charT, traitsT> {
+private:
+ typedef typename std::basic_streambuf<charT, traitsT> streambuf_type;
+
+public:
+ typedef typename streambuf_type::char_type char_type;
+ typedef typename streambuf_type::int_type int_type;
+ typedef typename streambuf_type::pos_type pos_type;
+ typedef typename streambuf_type::off_type off_type;
+ typedef typename streambuf_type::traits_type traits_type;
+
+ virtual ~basic_null_streambuf() {}
+
+protected:
+ std::streamsize xsputn(const char_type* s, std::streamsize count) override {
+ return count;
+ }
+
+ int_type overflow(int_type c) override {
+ return traits_type::not_eof(c);
+ }
+};
+
+class mask_stream {
+public:
+ explicit mask_stream(bool mask): mask_(mask) {}
+
+ operator bool() const { return mask_; }
+
+ template <typename charT, typename traitsT>
+ friend std::basic_ostream<charT, traitsT>&
+ operator<<(std::basic_ostream<charT, traitsT>& O, const mask_stream& F) {
+ int xindex = get_xindex();
+
+ std::basic_streambuf<charT, traitsT>* saved_streambuf =
+ static_cast<std::basic_streambuf<charT, traitsT>*>(O.pword(xindex));
+
+ if (F.mask_ && saved_streambuf) {
+ // re-enable by restoring saved streambuf
+ O.pword(xindex) = 0;
+ O.rdbuf(saved_streambuf);
+ }
+ else if (!F.mask_ && !saved_streambuf) {
+ // disable stream but save old streambuf
+ O.pword(xindex) = O.rdbuf();
+ O.rdbuf(get_null_streambuf<charT, traitsT>());
+ }
+
+ return O;
+ }
+
+private:
+ // our key for retrieve saved streambufs.
+ static int get_xindex() {
+ static int xindex = std::ios_base::xalloc();
+ return xindex;
+ }
+
+ template <typename charT, typename traitsT>
+ static std::basic_streambuf<charT, traitsT>* get_null_streambuf() {
+ static basic_null_streambuf<charT, traitsT> the_null_streambuf;
+ return &the_null_streambuf;
+ }
+
+ // true => do not filter
+ bool mask_;
+};
+
+} // namespace util
+} // namespace mc
+} // namespace nest
+
diff --git a/src/util/lexcmp_def.hpp b/src/util/lexcmp_def.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cfb297b1d16d3dbd7023bb7f300271eda0af1e41
--- /dev/null
+++ b/src/util/lexcmp_def.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+/*
+ * Macro definitions for defining comparison operators for
+ * record-like types.
+ *
+ * Use:
+ *
+ * To define comparison operations for a record type xyzzy
+ * with fields foo, bar and baz:
+ *
+ * DEFINE_LEXICOGRAPHIC_ORDERING(xyzzy,(a.foo,a.bar,a.baz),(b.foo,b.bar,b.baz))
+ *
+ * The explicit use of 'a' and 'b' in the second and third parameters
+ * is needed only to save a heroic amount of preprocessor macro
+ * deep magic.
+ *
+ */
+
+#include <tuple>
+
+#define DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(proxy,op,type,a_fields,b_fields) \
+inline bool operator op(const type& a,const type& b) { return proxy a_fields op proxy b_fields; }
+
+#define DEFINE_LEXICOGRAPHIC_ORDERING(type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::tie,<,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::tie,>,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::tie,<=,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::tie,>=,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::tie,!=,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::tie,==,type,a_fields,b_fields)
+
+#define DEFINE_LEXICOGRAPHIC_ORDERING_BY_VALUE(type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::make_tuple,<,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::make_tuple,>,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::make_tuple,<=,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::make_tuple,>=,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::make_tuple,!=,type,a_fields,b_fields) \
+DEFINE_LEXICOGRAPHIC_ORDERING_IMPL_(std::make_tuple,==,type,a_fields,b_fields)
+
diff --git a/src/util/optional.hpp b/src/util/optional.hpp
index 0d71c4a5d564b4a85db457a930f803388ee0861b..6b4d0f8b66cbdcf24c59914466a01621895ca2bb 100644
--- a/src/util/optional.hpp
+++ b/src/util/optional.hpp
@@ -29,7 +29,7 @@ namespace util {
template <typename X> struct optional;
struct optional_unset_error: std::runtime_error {
- explicit optional_unset_error(const std::string &what_str)
+ explicit optional_unset_error(const std::string& what_str)
: std::runtime_error(what_str)
{}
@@ -39,7 +39,7 @@ struct optional_unset_error: std::runtime_error {
};
struct optional_invalid_dereference: std::runtime_error {
- explicit optional_invalid_dereference(const std::string &what_str)
+ explicit optional_invalid_dereference(const std::string& what_str)
: std::runtime_error(what_str)
{}
@@ -226,6 +226,13 @@ namespace detail {
template <typename T>
using enable_unless_optional_t = enable_if_t<!is_optional<T>::value>;
+ // avoid nonnull address warnings when using operator| with e.g. char array constants
+ template <typename T>
+ bool decay_bool(const T* x) { return static_cast<bool>(x); }
+
+ template <typename T>
+ bool decay_bool(const T& x) { return static_cast<bool>(x); }
+
} // namespace detail
template <typename X>
@@ -332,17 +339,19 @@ struct optional<X&>: detail::optional_base<X&> {
template <typename T>
optional(optional<T&>& ot): base(ot.set,ot.ref()) {}
- template <typename Y,typename = typename std::enable_if<!detail::is_optional<Y>()>::type>
+ template <typename Y>
optional& operator=(Y& y) {
set = true;
- ref() = y;
+ data.construct(y);
return *this;
}
template <typename Y>
optional& operator=(optional<Y&>& o) {
set = o.set;
- data.construct(o);
+ if (o.set) {
+ data.construct(o.get());
+ }
return *this;
}
};
@@ -363,7 +372,7 @@ struct optional<void>: detail::optional_base<void> {
template <typename T>
optional(const optional<T>& o): base(o.set,true) {}
-
+
template <typename T>
optional& operator=(T) {
set = true;
@@ -396,7 +405,7 @@ typename std::enable_if<
>
>::type
operator|(A&& a,B&& b) {
- return a ? a : b;
+ return detail::decay_bool(a) ? a : b;
}
template <typename A,typename B>
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 9061fb8d4698c681d3cdd38986b081f6b47df5aa..10fb65efd052c1d9d7d7fdd5f50bdec0f93cb5a5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,73 +1,27 @@
-set(HEADERS
- ../src/cell.hpp
- ../src/cell_tree.hpp
- ../src/math.hpp
- ../src/point.hpp
- ../src/segment.hpp
- ../src/swcio.hpp
- ../src/tree.hpp
-)
# google test framework
-add_library(gtest gtest-all.cpp gtest.h)
+add_library(gtest gtest-all.cpp)
+# tests look for gtest.h here
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
-set(TEST_SOURCES
- # unit tests
- test_algorithms.cpp
- test_cell.cpp
- test_compartments.cpp
- test_event_queue.cpp
- test_fvm.cpp
- test_cell_group.cpp
- test_matrix.cpp
- test_mechanisms.cpp
- test_optional.cpp
- test_parameters.cpp
- test_point.cpp
- test_probe.cpp
- test_segment.cpp
- test_spikes.cpp
- test_stimulus.cpp
- test_swcio.cpp
- test_synapses.cpp
- test_tree.cpp
- test_uninitialized.cpp
- # unit test driver
- test.cpp
-)
+# Unit tests
+add_subdirectory(unit)
-set(VALIDATION_SOURCES
- # unit tests
- validate_ball_and_stick.cpp
- validate_soma.cpp
- #validate_synapses.cpp
+# Test validating models, possebly needing other software installed
+add_subdirectory(validation)
- # unit test driver
- validate.cpp
-)
+# Test for the internode communication (eg. mpi)
+add_subdirectory(global_communication)
-add_definitions("-DDATADIR=\"${CMAKE_SOURCE_DIR}/data\"")
-add_executable(test.exe ${TEST_SOURCES} ${HEADERS})
-add_executable(validate.exe ${VALIDATION_SOURCES} ${HEADERS})
-set(TARGETS test.exe validate.exe)
+# Proposed additional test types:
-foreach(target ${TARGETS})
- target_link_libraries(${target} LINK_PUBLIC cellalgo gtest)
+# Large test, employing the full simulator. validated using deltas on output data
- if(WITH_TBB)
- target_link_libraries(${target} LINK_PUBLIC ${TBB_LIBRARIES})
- endif()
+# Test to check integration between components
- if(WITH_MPI)
- target_link_libraries(${target} LINK_PUBLIC ${MPI_C_LIBRARIES})
- set_property(TARGET ${target} APPEND_STRING PROPERTY LINK_FLAGS "${MPI_C_LINK_FLAGS}")
- endif()
+# Tests for performance: This could include stand alone tests. These do not necessarily be run automatically
- set_target_properties(${target}
- PROPERTIES
- RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
- )
-endforeach()
+# Numbered tests based on bugs in the tracker
diff --git a/tests/global_communication/CMakeLists.txt b/tests/global_communication/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f7d000157e4cfd130399d1f718b67fbc5eb2f86d
--- /dev/null
+++ b/tests/global_communication/CMakeLists.txt
@@ -0,0 +1 @@
+# Nothing to be done yet
diff --git a/tests/test_cell_group.cpp b/tests/test_cell_group.cpp
deleted file mode 100644
index d9787a1935657391be82f2e674b13bab9ef908ed..0000000000000000000000000000000000000000
--- a/tests/test_cell_group.cpp
+++ /dev/null
@@ -1,48 +0,0 @@
-#include <limits>
-
-#include "gtest.h"
-
-#include <fvm_cell.hpp>
-#include <cell_group.hpp>
-
-nest::mc::cell make_cell() {
- using namespace nest::mc;
-
- // setup global state for the mechanisms
- mechanisms::setup_mechanism_helpers();
-
- nest::mc::cell cell;
-
- // 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->set_compartments(101);
-
- dendrite->mechanism("membrane").set("r_L", 100);
-
- // add stimulus
- cell.add_stimulus({1,1}, {5., 80., 0.3});
-
- cell.add_detector({0,0}, 0);
-
- return cell;
-}
-
-TEST(cell_group, test)
-{
- using namespace nest::mc;
-
- using cell_type = cell_group<fvm::fvm_cell<double, int>>;
-
- auto cell = cell_type{make_cell()};
-
- cell.advance(50, 0.01);
-
- // a bit lame...
- EXPECT_EQ(cell.spikes().size(), 4u);
-}
-
diff --git a/tests/test_synapses.cpp b/tests/test_synapses.cpp
deleted file mode 100644
index 87ac3aeca2e828035dea39c7cab4bf156a4e04bf..0000000000000000000000000000000000000000
--- a/tests/test_synapses.cpp
+++ /dev/null
@@ -1,83 +0,0 @@
-#include "gtest.h"
-#include "util.hpp"
-
-#include <cell.hpp>
-#include <fvm_cell.hpp>
-
-// compares results with those generated by nrn/ball_and_stick.py
-TEST(synapses, add_to_cell)
-{
- using namespace nest::mc;
-
- nest::mc::cell cell;
-
- // setup global state for the mechanisms
- nest::mc::mechanisms::setup_mechanism_helpers();
-
- // Soma with diameter 12.6157 um and HH channel
- auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism(hh_parameters());
-
- cell.add_synapse({0, 0.1});
- cell.add_synapse({1, 0.2});
- cell.add_synapse({0, 0.3});
-
- EXPECT_EQ(3u, cell.synapses().size());
- EXPECT_EQ(cell.synapses()[0].segment, 0);
- EXPECT_EQ(cell.synapses()[0].position, 0.1);
- EXPECT_EQ(cell.synapses()[1].segment, 1);
- EXPECT_EQ(cell.synapses()[1].position, 0.2);
- EXPECT_EQ(cell.synapses()[2].segment, 0);
- EXPECT_EQ(cell.synapses()[2].position, 0.3);
-}
-
-// compares results with those generated by nrn/ball_and_stick.py
-TEST(synapses, basic_state)
-{
- using namespace nest::mc;
-
- // setup global state for the mechanisms
- nest::mc::mechanisms::setup_mechanism_helpers();
-
- using synapse_type = mechanisms::ExpSyn::mechanism_ExpSyn<double, int>;
- auto num_syn = 4;
- synapse_type::index_type indexes(num_syn);
- synapse_type::vector_type voltage(num_syn, -65.0);
- synapse_type::vector_type current(num_syn, 1.0);
- auto mech = mechanisms::make_mechanism<synapse_type>( voltage, current, indexes );
-
- auto ptr = dynamic_cast<synapse_type*>(mech.get());
-
- // parameters initialized to default values
- for(auto e : ptr->e) {
- EXPECT_EQ(e, 0.);
- }
- for(auto tau : ptr->tau) {
- EXPECT_EQ(tau, 2.0);
- }
-
- // current and voltage vectors correctly hooked up
- for(auto v : ptr->vec_v_) {
- EXPECT_EQ(v, -65.);
- }
- for(auto i : ptr->vec_i_) {
- EXPECT_EQ(i, 1.0);
- }
-
- // should be initialized to NaN
- for(auto g : ptr->g) {
- EXPECT_NE(g, g);
- }
-
- // initialize state then check g has been set to zero
- ptr->nrn_init();
- for(auto g : ptr->g) {
- EXPECT_EQ(g, 0.);
- }
-
- // call net_receive on two of the synapses
- ptr->net_receive(1, 3.14);
- ptr->net_receive(3, 1.04);
- EXPECT_EQ(ptr->g[1], 3.14);
- EXPECT_EQ(ptr->g[3], 1.04);
-}
diff --git a/tests/util.hpp b/tests/test_util.hpp
similarity index 86%
rename from tests/util.hpp
rename to tests/test_util.hpp
index cfa89896af8039e6188d8fbc460f412d3407a1b3..0f9935d8c517758a4aa1308921b382ea56b596fc 100644
--- a/tests/util.hpp
+++ b/tests/test_util.hpp
@@ -55,28 +55,6 @@ void write_vis_file(const std::string& fname, std::vector<std::vector<double>> v
}
}
-[[gnu::unused]] static
-nlohmann::json
-load_spike_data(const std::string& input_name)
-{
- nlohmann::json cell_data;
- std::ifstream fid(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)
{
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fb1b9369d4e38498c8953fc0ce2f16c1dbc715dc
--- /dev/null
+++ b/tests/unit/CMakeLists.txt
@@ -0,0 +1,63 @@
+set(HEADERS
+ ${PROJECT_SOURCE_DIR}/src/cell.hpp
+ ${PROJECT_SOURCE_DIR}/src/cell_tree.hpp
+ ${PROJECT_SOURCE_DIR}/src/math.hpp
+ ${PROJECT_SOURCE_DIR}/src/point.hpp
+ ${PROJECT_SOURCE_DIR}/src/segment.hpp
+ ${PROJECT_SOURCE_DIR}/src/swcio.hpp
+ ${PROJECT_SOURCE_DIR}/src/tree.hpp
+)
+
+set(TEST_SOURCES
+ # unit tests
+ test_algorithms.cpp
+ test_double_buffer.cpp
+ test_cell.cpp
+ test_compartments.cpp
+ test_event_queue.cpp
+ test_fvm.cpp
+ test_cell_group.cpp
+ test_lexcmp.cpp
+ test_mask_stream.cpp
+ test_matrix.cpp
+ test_mechanisms.cpp
+ test_optional.cpp
+ test_parameters.cpp
+ test_point.cpp
+ test_probe.cpp
+ test_segment.cpp
+ test_spikes.cpp
+ test_spike_store.cpp
+ test_stimulus.cpp
+ test_swcio.cpp
+ test_synapses.cpp
+ test_tree.cpp
+ test_uninitialized.cpp
+
+ # unit test driver
+ test.cpp
+)
+
+add_definitions("-DDATADIR=\"${CMAKE_SOURCE_DIR}/data\"")
+add_executable(test.exe ${TEST_SOURCES} ${HEADERS})
+
+set(TARGETS test.exe)
+
+foreach(target ${TARGETS})
+ target_link_libraries(${target} LINK_PUBLIC cellalgo gtest)
+
+ if(WITH_TBB)
+ target_link_libraries(${target} LINK_PUBLIC ${TBB_LIBRARIES})
+ endif()
+
+ if(WITH_MPI)
+ target_link_libraries(${target} LINK_PUBLIC ${MPI_C_LIBRARIES})
+ set_property(TARGET ${target} APPEND_STRING PROPERTY LINK_FLAGS "${MPI_C_LINK_FLAGS}")
+ endif()
+
+ set_target_properties(${target}
+ PROPERTIES
+ RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
+ )
+endforeach()
+
diff --git a/tests/test.cpp b/tests/unit/test.cpp
similarity index 100%
rename from tests/test.cpp
rename to tests/unit/test.cpp
diff --git a/tests/test_algorithms.cpp b/tests/unit/test_algorithms.cpp
similarity index 99%
rename from tests/test_algorithms.cpp
rename to tests/unit/test_algorithms.cpp
index 1a363be6cc45313c7f76a6f6c7ff6ada3ba715b1..ebe19c3752e9e4eacb23008e467bca9a1fd95f79 100644
--- a/tests/test_algorithms.cpp
+++ b/tests/unit/test_algorithms.cpp
@@ -3,7 +3,7 @@
#include "gtest.h"
#include "algorithms.hpp"
-#include "util.hpp"
+#include "../test_util.hpp"
#include "util/debug.hpp"
diff --git a/tests/test_cell.cpp b/tests/unit/test_cell.cpp
similarity index 90%
rename from tests/test_cell.cpp
rename to tests/unit/test_cell.cpp
index 88c8dcaced167d4030529adb8a348a116158c17f..3102d97435d3095b4ef91b885b040ac084f37a5f 100644
--- a/tests/test_cell.cpp
+++ b/tests/unit/test_cell.cpp
@@ -1,6 +1,6 @@
#include "gtest.h"
-#include "../src/cell.hpp"
+#include "cell.hpp"
TEST(cell_type, soma)
{
@@ -59,7 +59,7 @@ TEST(cell_type, add_segment)
);
c.add_cable(0, std::move(seg));
- EXPECT_EQ(c.num_segments(), 2);
+ EXPECT_EQ(c.num_segments(), 2u);
}
// add segment on the fly
@@ -78,7 +78,7 @@ TEST(cell_type, add_segment)
segmentKind::dendrite, cable_radius, cable_radius, cable_length
);
- EXPECT_EQ(c.num_segments(), 2);
+ EXPECT_EQ(c.num_segments(), 2u);
}
{
cell c;
@@ -97,7 +97,7 @@ TEST(cell_type, add_segment)
std::vector<double>{cable_length, cable_length, cable_length}
);
- EXPECT_EQ(c.num_segments(), 2);
+ EXPECT_EQ(c.num_segments(), 2u);
}
}
@@ -137,7 +137,7 @@ TEST(cell_type, multiple_cables)
c.add_cable(1, seg(segmentKind::dendrite));
c.add_cable(1, seg(segmentKind::dendrite));
- EXPECT_EQ(c.num_segments(), 5);
+ EXPECT_EQ(c.num_segments(), 5u);
// 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
@@ -148,12 +148,14 @@ TEST(cell_type, multiple_cables)
const auto model = c.model();
auto const& con = model.tree;
+ auto no_parent = cell_tree::no_parent;
+
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.parent(0), no_parent);
+ EXPECT_EQ(con.parent(1), 0u);
+ EXPECT_EQ(con.parent(2), 0u);
+ EXPECT_EQ(con.parent(3), 1u);
+ EXPECT_EQ(con.parent(4), 1u);
EXPECT_EQ(con.num_children(0), 2u);
EXPECT_EQ(con.num_children(1), 2u);
EXPECT_EQ(con.num_children(2), 0u);
diff --git a/tests/unit/test_cell_group.cpp b/tests/unit/test_cell_group.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3aed222029f4535a85d9713849cff184a3684531
--- /dev/null
+++ b/tests/unit/test_cell_group.cpp
@@ -0,0 +1,74 @@
+#include "gtest.h"
+
+#include <common_types.hpp>
+#include <fvm_cell.hpp>
+#include <cell_group.hpp>
+
+nest::mc::cell make_cell() {
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // 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->set_compartments(101);
+
+ dendrite->mechanism("membrane").set("r_L", 100);
+
+ cell.add_detector({0, 0}, 0);
+ cell.add_stimulus({1, 1}, {5., 80., 0.3});
+
+ return cell;
+}
+
+TEST(cell_group, test)
+{
+ using namespace nest::mc;
+
+ using cell_group_type = cell_group<fvm::fvm_cell<double, cell_local_size_type>>;
+ auto group = cell_group_type{0, make_cell()};
+
+ group.advance(50, 0.01);
+
+ // a bit lame...
+ EXPECT_EQ(group.spikes().size(), 4u);
+}
+
+TEST(cell_group, sources)
+{
+ using namespace nest::mc;
+
+ // TODO: extend to multi-cell cell groups when the time comes
+
+ using cell_group_type = cell_group<fvm::fvm_cell<double, cell_local_size_type>>;
+
+ auto cell = make_cell();
+ EXPECT_EQ(cell.detectors().size(), 1u);
+ // add another detector on the cell to make things more interesting
+ cell.add_detector({1, 0.3}, 2.3);
+
+ cell_gid_type first_gid = 37u;
+ auto group = cell_group_type{first_gid, cell};
+
+ // expect group sources to be lexicographically sorted by source id
+ // with gids in cell group's range and indices starting from zero
+
+ const auto& sources = group.spike_sources();
+ for (unsigned i = 0; i<sources.size(); ++i) {
+ auto id = sources[i].source_id;
+ if (i==0) {
+ EXPECT_EQ(id.gid, first_gid);
+ EXPECT_EQ(id.index, 0u);
+ }
+ else {
+ auto prev = sources[i-1].source_id;
+ EXPECT_GT(id, prev);
+ EXPECT_EQ(id.index, id.gid==prev.gid? prev.index+1: 0u);
+ }
+ }
+}
diff --git a/tests/test_compartments.cpp b/tests/unit/test_compartments.cpp
similarity index 87%
rename from tests/test_compartments.cpp
rename to tests/unit/test_compartments.cpp
index c109f9cda9a335ce32f695f85ccd0295446a5d92..c167b16cce5585a43efc13fc9f136025b21f7811 100644
--- a/tests/test_compartments.cpp
+++ b/tests/unit/test_compartments.cpp
@@ -15,13 +15,13 @@ TEST(compartments, compartment)
{
nest::mc::compartment c(100, 1.2, 2.1, 2.2);
- EXPECT_EQ(c.index, 100);
+ EXPECT_EQ(c.index, 100u);
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.index, 100u);
EXPECT_EQ(c2.length, 1.2);
EXPECT_EQ(left(c2.radius), 2.1);
EXPECT_EQ(right(c2.radius), 2.2);
@@ -29,7 +29,7 @@ TEST(compartments, compartment)
{
nest::mc::compartment c{100, 1, 2, 3};
- EXPECT_EQ(c.index, 100);
+ EXPECT_EQ(c.index, 100u);
EXPECT_EQ(c.length, 1.);
EXPECT_EQ(left(c.radius), 2.);
EXPECT_EQ(right(c.radius), 3.);
@@ -54,7 +54,7 @@ TEST(compartments, compartment_iterator)
++it;
{
auto c = *it;
- EXPECT_EQ(c.index, 1);
+ EXPECT_EQ(c.index, 1u);
EXPECT_EQ(left(c.radius), 3.0);
EXPECT_EQ(right(c.radius), 5.0);
EXPECT_EQ(c.length, 2.5);
@@ -65,7 +65,7 @@ TEST(compartments, compartment_iterator)
// returned iterator should be unchanged
{
auto c = *(it++);
- EXPECT_EQ(c.index, 1);
+ EXPECT_EQ(c.index, 1u);
EXPECT_EQ(left(c.radius), 3.0);
EXPECT_EQ(right(c.radius), 5.0);
EXPECT_EQ(c.length, 2.5);
@@ -73,7 +73,7 @@ TEST(compartments, compartment_iterator)
// while the iterator itself was updated
{
auto c = *it;
- EXPECT_EQ(c.index, 2);
+ EXPECT_EQ(c.index, 2u);
EXPECT_EQ(left(c.radius), 5.0);
EXPECT_EQ(right(c.radius), 7.0);
EXPECT_EQ(c.length, 2.5);
@@ -84,7 +84,7 @@ TEST(compartments, compartment_iterator)
// copy iterator
auto it2 = it;
auto c = *it2;
- EXPECT_EQ(c.index, 2);
+ EXPECT_EQ(c.index, 2u);
EXPECT_EQ(left(c.radius), 5.0);
EXPECT_EQ(right(c.radius), 7.0);
EXPECT_EQ(c.length, 2.5);
@@ -96,7 +96,7 @@ TEST(compartments, compartment_iterator)
// check the copy has updated correctly when incremented
c= *it2;
- EXPECT_EQ(c.index, 3);
+ EXPECT_EQ(c.index, 3u);
EXPECT_EQ(left(c.radius), 7.0);
EXPECT_EQ(right(c.radius), 9.0);
EXPECT_EQ(c.length, 2.5);
@@ -108,11 +108,11 @@ 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_EQ((*rng.begin()).index, 0u);
+ EXPECT_EQ((*rng.end()).index, 10u);
EXPECT_NE(rng.begin(), rng.end());
- auto count = 0;
+ unsigned count = 0;
for(auto c : rng) {
EXPECT_EQ(c.index, count);
auto er = 1.0 + double(count)/10.;
@@ -121,7 +121,7 @@ TEST(compartments, compartment_range)
EXPECT_EQ(c.length, 1.0);
++count;
}
- EXPECT_EQ(count, 10);
+ EXPECT_EQ(count, 10u);
}
// test case of zero length range
diff --git a/tests/unit/test_double_buffer.cpp b/tests/unit/test_double_buffer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1689cc2faac17331694c78830c22cd1ea57777fe
--- /dev/null
+++ b/tests/unit/test_double_buffer.cpp
@@ -0,0 +1,58 @@
+#include "gtest.h"
+
+#include <util/double_buffer.hpp>
+
+// not much to test here: just test that values passed into the constructor
+// are correctly stored in members
+TEST(double_buffer, exchange_and_get)
+{
+ using namespace nest::mc::util;
+
+ double_buffer<int> buf;
+
+ buf.get() = 2134;
+ buf.exchange();
+ buf.get() = 8990;
+ buf.exchange();
+
+ EXPECT_EQ(buf.get(), 2134);
+ EXPECT_EQ(buf.other(), 8990);
+ buf.exchange();
+ EXPECT_EQ(buf.get(), 8990);
+ EXPECT_EQ(buf.other(), 2134);
+ buf.exchange();
+ EXPECT_EQ(buf.get(), 2134);
+ EXPECT_EQ(buf.other(), 8990);
+}
+
+TEST(double_buffer, assign_get_other)
+{
+ using namespace nest::mc::util;
+
+ double_buffer<std::string> buf;
+
+ buf.get() = "1";
+ buf.other() = "2";
+
+ EXPECT_EQ(buf.get(), "1");
+ EXPECT_EQ(buf.other(), "2");
+}
+
+TEST(double_buffer, non_pod)
+{
+ using namespace nest::mc::util;
+
+ double_buffer<std::string> buf;
+
+ buf.get() = "1";
+ buf.other() = "2";
+
+ EXPECT_EQ(buf.get(), "1");
+ EXPECT_EQ(buf.other(), "2");
+ buf.exchange();
+ EXPECT_EQ(buf.get(), "2");
+ EXPECT_EQ(buf.other(), "1");
+ buf.exchange();
+ EXPECT_EQ(buf.get(), "1");
+ EXPECT_EQ(buf.other(), "2");
+}
diff --git a/tests/test_event_queue.cpp b/tests/unit/test_event_queue.cpp
similarity index 63%
rename from tests/test_event_queue.cpp
rename to tests/unit/test_event_queue.cpp
index 5ac4c75747d761638b0bbd4b90a9b6fed6fa3c18..f4e31b26715fa2aa3cb690720442cc68f2b18fe2 100644
--- a/tests/test_event_queue.cpp
+++ b/tests/unit/test_event_queue.cpp
@@ -7,14 +7,14 @@
TEST(event_queue, push)
{
using namespace nest::mc;
- using ps_event_queue = event_queue<postsynaptic_spike_event>;
+ using ps_event_queue = event_queue<postsynaptic_spike_event<float>>;
ps_event_queue q;
- q.push({1u, 2.f, 2.f});
- q.push({4u, 1.f, 2.f});
- q.push({8u, 20.f, 2.f});
- q.push({2u, 8.f, 2.f});
+ q.push({{1u, 0u}, 2.f, 2.f});
+ q.push({{4u, 1u}, 1.f, 2.f});
+ q.push({{8u, 2u}, 20.f, 2.f});
+ q.push({{2u, 3u}, 8.f, 2.f});
std::vector<float> times;
while(q.size()) {
@@ -31,13 +31,13 @@ TEST(event_queue, push)
TEST(event_queue, push_range)
{
using namespace nest::mc;
- using ps_event_queue = event_queue<postsynaptic_spike_event>;
+ using ps_event_queue = event_queue<postsynaptic_spike_event<float>>;
- postsynaptic_spike_event events[] = {
- {1u, 2.f, 2.f},
- {4u, 1.f, 2.f},
- {8u, 20.f, 2.f},
- {2u, 8.f, 2.f}
+ postsynaptic_spike_event<float> events[] = {
+ {{1u, 0u}, 2.f, 2.f},
+ {{4u, 1u}, 1.f, 2.f},
+ {{8u, 2u}, 20.f, 2.f},
+ {{2u, 3u}, 8.f, 2.f}
};
ps_event_queue q;
@@ -56,13 +56,20 @@ TEST(event_queue, push_range)
TEST(event_queue, pop_if_before)
{
using namespace nest::mc;
- using ps_event_queue = event_queue<postsynaptic_spike_event>;
+ using ps_event_queue = event_queue<postsynaptic_spike_event<float>>;
- postsynaptic_spike_event events[] = {
- {1u, 1.f, 2.f},
- {2u, 2.f, 2.f},
- {3u, 3.f, 2.f},
- {4u, 4.f, 2.f}
+ cell_member_type target[4] = {
+ {1u, 0u},
+ {4u, 1u},
+ {8u, 2u},
+ {2u, 3u}
+ };
+
+ postsynaptic_spike_event<float> events[] = {
+ {target[0], 1.f, 2.f},
+ {target[1], 2.f, 2.f},
+ {target[2], 3.f, 2.f},
+ {target[3], 4.f, 2.f}
};
ps_event_queue q;
@@ -76,12 +83,12 @@ TEST(event_queue, pop_if_before)
auto e2 = q.pop_if_before(5.);
EXPECT_TRUE(e2);
- EXPECT_EQ(e2->target, 1u);
+ EXPECT_EQ(e2->target, target[0]);
EXPECT_EQ(q.size(), 3u);
auto e3 = q.pop_if_before(5.);
EXPECT_TRUE(e3);
- EXPECT_EQ(e3->target, 2u);
+ EXPECT_EQ(e3->target, target[1]);
EXPECT_EQ(q.size(), 2u);
auto e4 = q.pop_if_before(2.5);
@@ -90,7 +97,7 @@ TEST(event_queue, pop_if_before)
auto e5 = q.pop_if_before(5.);
EXPECT_TRUE(e5);
- EXPECT_EQ(e5->target, 3u);
+ EXPECT_EQ(e5->target, target[2]);
EXPECT_EQ(q.size(), 1u);
q.pop_if_before(5.);
diff --git a/tests/test_fvm.cpp b/tests/unit/test_fvm.cpp
similarity index 89%
rename from tests/test_fvm.cpp
rename to tests/unit/test_fvm.cpp
index bc7374cbd0ea8f16b2b8adb02c4d27382aa83c61..c0ed438f20841d40a9f2653e38424aeeba4a64ea 100644
--- a/tests/test_fvm.cpp
+++ b/tests/unit/test_fvm.cpp
@@ -1,20 +1,19 @@
#include <fstream>
#include "gtest.h"
-#include "util.hpp"
+#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
+#include "../test_util.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
@@ -42,7 +41,7 @@ TEST(fvm, cable)
cell.segment(1)->set_compartments(4);
cell.segment(2)->set_compartments(4);
- using fvm_cell = fvm::fvm_cell<double, int>;
+ using fvm_cell = fvm::fvm_cell<double, cell_lid_type>;
fvm_cell fvcell(cell);
auto& J = fvcell.jacobian();
@@ -64,9 +63,6 @@ TEST(fvm, init)
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;
@@ -95,7 +91,7 @@ TEST(fvm, init)
cell.segment(1)->set_compartments(10);
- using fvm_cell = fvm::fvm_cell<double, int>;
+ using fvm_cell = fvm::fvm_cell<double, cell_lid_type>;
fvm_cell fvcell(cell);
auto& J = fvcell.jacobian();
EXPECT_EQ(J.size(), 11u);
diff --git a/tests/unit/test_lexcmp.cpp b/tests/unit/test_lexcmp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2b5c476d7983f7726349953d15e34509aab031c
--- /dev/null
+++ b/tests/unit/test_lexcmp.cpp
@@ -0,0 +1,110 @@
+#include "gtest.h"
+
+#include <util/lexcmp_def.hpp>
+
+struct lexcmp_test_one {
+ int foo;
+};
+
+DEFINE_LEXICOGRAPHIC_ORDERING(lexcmp_test_one, (a.foo), (b.foo))
+
+TEST(lexcmp_def,one) {
+ lexcmp_test_one p{3}, q{4}, r{4};
+
+ EXPECT_LE(p,q);
+ EXPECT_LT(p,q);
+ EXPECT_NE(p,q);
+ EXPECT_GE(q,p);
+ EXPECT_GT(q,p);
+
+ EXPECT_LE(q,r);
+ EXPECT_GE(q,r);
+ EXPECT_EQ(q,r);
+}
+
+struct lexcmp_test_three {
+ int x;
+ std::string y;
+ double z;
+};
+
+// test fields in reverse order: z, y, x
+DEFINE_LEXICOGRAPHIC_ORDERING(lexcmp_test_three, (a.z,a.y,a.x), (b.z,b.y,b.x))
+
+TEST(lexcmp_def,three) {
+ lexcmp_test_three p{1,"foo",2};
+ lexcmp_test_three q{1,"foo",3};
+ lexcmp_test_three r{1,"bar",2};
+ lexcmp_test_three s{5,"foo",2};
+
+ EXPECT_LE(p,q);
+ EXPECT_LT(p,q);
+ EXPECT_NE(p,q);
+ EXPECT_GE(q,p);
+ EXPECT_GT(q,p);
+
+ EXPECT_LE(r,p);
+ EXPECT_LT(r,p);
+ EXPECT_NE(p,r);
+ EXPECT_GE(p,r);
+ EXPECT_GT(p,r);
+
+ EXPECT_LE(p,s);
+ EXPECT_LT(p,s);
+ EXPECT_NE(p,s);
+ EXPECT_GE(s,p);
+ EXPECT_GT(s,p);
+}
+
+// test fields accessed by reference-returning member function
+
+class lexcmp_test_refmemfn {
+public:
+ explicit lexcmp_test_refmemfn(int foo): foo_(foo) {}
+
+ const int& foo() const { return foo_; }
+ int& foo() { return foo_; }
+
+private:
+ int foo_;
+};
+
+DEFINE_LEXICOGRAPHIC_ORDERING(lexcmp_test_refmemfn, (a.foo()), (b.foo()))
+
+TEST(lexcmp_def,refmemfn) {
+ lexcmp_test_refmemfn p{3};
+ const lexcmp_test_refmemfn q{4};
+
+ EXPECT_LE(p,q);
+ EXPECT_LT(p,q);
+ EXPECT_NE(p,q);
+ EXPECT_GE(q,p);
+ EXPECT_GT(q,p);
+}
+
+// test comparison via proxy tuple object
+
+class lexcmp_test_valmemfn {
+public:
+ explicit lexcmp_test_valmemfn(int foo, int bar): foo_(foo), bar_(bar) {}
+ int foo() const { return foo_; }
+ int bar() const { return bar_; }
+
+private:
+ int foo_;
+ int bar_;
+};
+
+DEFINE_LEXICOGRAPHIC_ORDERING_BY_VALUE(lexcmp_test_valmemfn, (a.foo(),a.bar()), (b.foo(),b.bar()))
+
+TEST(lexcmp_def,proxy) {
+ lexcmp_test_valmemfn p{3,2}, q{3,4};
+
+ EXPECT_LE(p,q);
+ EXPECT_LT(p,q);
+ EXPECT_NE(p,q);
+ EXPECT_GE(q,p);
+ EXPECT_GT(q,p);
+}
+
+
diff --git a/tests/unit/test_mask_stream.cpp b/tests/unit/test_mask_stream.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8ca8089eab1f7486204ce0df70fb953000b301a0
--- /dev/null
+++ b/tests/unit/test_mask_stream.cpp
@@ -0,0 +1,68 @@
+#include "gtest.h"
+
+#include <sstream>
+
+#include <util/ioutil.hpp>
+
+using namespace nest::mc::util;
+
+TEST(mask_stream,nomask) {
+ // expect mask_stream(true) on a new stream not to change rdbuf.
+ std::ostringstream s;
+ auto sbuf = s.rdbuf();
+ s << mask_stream(true);
+ EXPECT_EQ(sbuf, s.rdbuf());
+}
+
+TEST(mask_stream,mask) {
+ // masked stream should produce no ouptut
+ std::ostringstream s;
+ s << "one";
+ s << mask_stream(false);
+
+ s << "two";
+ EXPECT_EQ(s.str(), "one");
+
+ s << mask_stream(true);
+ s << "three";
+ EXPECT_EQ(s.str(), "onethree");
+}
+
+TEST(mask_stream,mask_multi) {
+ // mark_stream(false) should be idempotent
+
+ std::ostringstream s;
+ auto sbuf1 = s.rdbuf();
+
+ s << "foo";
+ s << mask_stream(false);
+ auto sbuf2 = s.rdbuf();
+
+ s << "bar";
+ s << mask_stream(false);
+ auto sbuf3 = s.rdbuf();
+ EXPECT_EQ(sbuf2, sbuf3);
+
+ s << "baz";
+ s << mask_stream(true);
+ auto sbuf4 = s.rdbuf();
+ EXPECT_EQ(sbuf1, sbuf4);
+
+ s << "xyzzy";
+ EXPECT_EQ(s.str(), "fooxyzzy");
+}
+
+TEST(mask_stream,fmt) {
+ // expect formatting to be preserved across masks.
+
+ std::ostringstream s;
+ s.precision(1);
+
+ s << mask_stream(false);
+ EXPECT_EQ(s.precision(), 1);
+ s.precision(2);
+
+ s << mask_stream(true);
+ EXPECT_EQ(s.precision(), 2);
+}
+
diff --git a/tests/test_matrix.cpp b/tests/unit/test_matrix.cpp
similarity index 100%
rename from tests/test_matrix.cpp
rename to tests/unit/test_matrix.cpp
diff --git a/tests/test_mechanisms.cpp b/tests/unit/test_mechanisms.cpp
similarity index 57%
rename from tests/test_mechanisms.cpp
rename to tests/unit/test_mechanisms.cpp
index aef84890275263ca7ef9570c1017864736a6f331..377f683fa8ff0831b0cdff3f94be96710d966d0c 100644
--- a/tests/test_mechanisms.cpp
+++ b/tests/unit/test_mechanisms.cpp
@@ -1,25 +1,17 @@
#include "gtest.h"
-#include "../src/mechanism_interface.hpp"
-#include "../src/matrix.hpp"
+//#include "../src/mechanism_interface.hpp"
+#include "mechanism_catalogue.hpp"
+#include "matrix.hpp"
TEST(mechanisms, helpers) {
- nest::mc::mechanisms::setup_mechanism_helpers();
-
- EXPECT_EQ(nest::mc::mechanisms::mechanism_helpers.size(), 2u);
+ using namespace nest::mc;
+ using catalogue = mechanisms::catalogue<double, int>;
// 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
- );
+ EXPECT_TRUE(catalogue::has("hh"));
+ EXPECT_TRUE(catalogue::has("pas"));
- //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();
@@ -28,10 +20,17 @@ TEST(mechanisms, helpers) {
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);
+ auto mech = catalogue::make("hh", vec_v, vec_i, node_indices);
EXPECT_EQ(mech->name(), "hh");
EXPECT_EQ(mech->size(), 5u);
//EXPECT_EQ(mech->matrix_, &matrix);
+
+ // check that an out_of_range exception is thrown if an invalid mechanism is
+ // requested
+ ASSERT_THROW(
+ catalogue::make("dachshund", vec_v, vec_i, node_indices),
+ std::out_of_range
+ );
+ //0 1 2 3 4 5 6 7 8 9
}
diff --git a/tests/test_optional.cpp b/tests/unit/test_optional.cpp
similarity index 75%
rename from tests/test_optional.cpp
rename to tests/unit/test_optional.cpp
index 9c6f35cfa658ab1c45cee53250eaa747d89dea25..4610f0c48cf6b317f97e8a895f9bb09245ee7220 100644
--- a/tests/test_optional.cpp
+++ b/tests/unit/test_optional.cpp
@@ -24,24 +24,30 @@ TEST(optionalm,unset_throw) {
optional<int> a;
int check=10;
- try { a.get(); }
- catch(optional_unset_error &e) {
+ try {
+ a.get();
+ }
+ catch (optional_unset_error& e) {
++check;
}
EXPECT_EQ(11,check);
check=20;
a=2;
- try { a.get(); }
- catch(optional_unset_error &e) {
+ try {
+ a.get();
+ }
+ catch (optional_unset_error& e) {
++check;
}
EXPECT_EQ(20,check);
check=30;
a.reset();
- try { a.get(); }
- catch(optional_unset_error &e) {
+ try {
+ a.get();
+ }
+ catch (optional_unset_error& e) {
++check;
}
EXPECT_EQ(31,check);
@@ -53,7 +59,7 @@ TEST(optionalm,deref) {
explicit foo(int a_): a(a_) {}
double value() { return 3.0*a; }
};
-
+
optional<foo> f=foo(2);
EXPECT_EQ(6.0,f->value());
EXPECT_EQ(2,(*f).a);
@@ -66,13 +72,13 @@ TEST(optionalm,ctor_conv) {
TEST(optionalm,ctor_ref) {
int v=10;
- optional<int &> a(v);
+ optional<int&> a(v);
EXPECT_EQ(10,a.get());
v=20;
EXPECT_EQ(20,a.get());
- optional<int &> b(a),c=b,d=v;
+ optional<int&> b(a),c=b,d=v;
EXPECT_EQ(&(a.get()),&(b.get()));
EXPECT_EQ(&(a.get()),&(c.get()));
EXPECT_EQ(&(a.get()),&(d.get()));
@@ -88,22 +94,37 @@ TEST(optionalm,assign_returns) {
EXPECT_EQ(&a,bp);
}
-namespace {
- struct nomove {
- int value;
+TEST(optionalm,assign_reference) {
+ double a=3.0;
+ optional<double&> ar;
+ optional<double&> br;
- nomove(): value(0) {}
- nomove(int i): value(i) {}
- nomove(const nomove &n): value(n.value) {}
- nomove(nomove &&n) = delete;
+ ar = a;
+ EXPECT_TRUE(ar);
+ *ar = 5.0;
+ EXPECT_EQ(5.0, a);
- nomove &operator=(const nomove &n) { value=n.value; return *this; }
+ br = ar;
+ EXPECT_TRUE(br);
- bool operator==(const nomove &them) const { return them.value==value; }
- bool operator!=(const nomove &them) const { return !(*this==them); }
- };
+ *br = 7.0;
+ EXPECT_EQ(7.0, a);
}
-
+
+struct nomove {
+ int value;
+
+ nomove(): value(0) {}
+ nomove(int i): value(i) {}
+ nomove(const nomove& n): value(n.value) {}
+ nomove(nomove&& n) = delete;
+
+ nomove& operator=(const nomove& n) { value=n.value; return *this; }
+
+ bool operator==(const nomove& them) const { return them.value==value; }
+ bool operator!=(const nomove& them) const { return !(*this==them); }
+};
+
TEST(optionalm,ctor_nomove) {
optional<nomove> a(nomove(3));
EXPECT_EQ(nomove(3),a.get());
@@ -116,30 +137,28 @@ TEST(optionalm,ctor_nomove) {
EXPECT_EQ(nomove(4),b.get());
}
-namespace {
- struct nocopy {
- int value;
-
- nocopy(): value(0) {}
- nocopy(int i): value(i) {}
- nocopy(const nocopy &n) = delete;
- nocopy(nocopy &&n) {
- value=n.value;
- n.value=0;
- }
-
- nocopy &operator=(const nocopy &n) = delete;
- nocopy &operator=(nocopy &&n) {
- value=n.value;
- n.value=-1;
- return *this;
- }
-
- bool operator==(const nocopy &them) const { return them.value==value; }
- bool operator!=(const nocopy &them) const { return !(*this==them); }
- };
-}
-
+struct nocopy {
+ int value;
+
+ nocopy(): value(0) {}
+ nocopy(int i): value(i) {}
+ nocopy(const nocopy& n) = delete;
+ nocopy(nocopy&& n) {
+ value=n.value;
+ n.value=0;
+ }
+
+ nocopy& operator=(const nocopy& n) = delete;
+ nocopy& operator=(nocopy&& n) {
+ value=n.value;
+ n.value=-1;
+ return *this;
+ }
+
+ bool operator==(const nocopy& them) const { return them.value==value; }
+ bool operator!=(const nocopy& them) const { return !(*this==them); }
+};
+
TEST(optionalm,ctor_nocopy) {
optional<nocopy> a(nocopy(5));
EXPECT_EQ(nocopy(5),a.get());
@@ -152,12 +171,10 @@ TEST(optionalm,ctor_nocopy) {
EXPECT_EQ(nocopy(6),b.get());
}
-namespace {
- optional<double> odd_half(int n) {
- optional<double> h;
- if (n%2==1) h=n/2.0;
- return h;
- }
+static optional<double> odd_half(int n) {
+ optional<double> h;
+ if (n%2==1) h=n/2.0;
+ return h;
}
TEST(optionalm,bind) {
@@ -199,7 +216,7 @@ TEST(optionalm,void) {
TEST(optionalm,bind_to_void) {
optional<int> a,b(3);
-
+
int call_count=0;
auto vf=[&call_count](int i) -> void { ++call_count; };
@@ -213,12 +230,14 @@ TEST(optionalm,bind_to_void) {
EXPECT_TRUE((bool)x);
EXPECT_EQ(1,call_count);
}
-
+
TEST(optionalm,bind_to_optional_void) {
optional<int> a,b(3),c(4);
-
+
int count=0;
- auto count_if_odd=[&count](int i) { return i%2?(++count,optional<void>(true)):optional<void>(); };
+ auto count_if_odd=[&count](int i) {
+ return i%2?(++count,optional<void>(true)):optional<void>();
+ };
auto x=a >> count_if_odd;
EXPECT_EQ(typeid(optional<void>),typeid(x));
@@ -238,16 +257,14 @@ TEST(optionalm,bind_to_optional_void) {
TEST(optionalm,bind_with_ref) {
optional<int> a=10;
- a >> [](int &v) {++v; };
+ a >> [](int& v) { ++v; };
EXPECT_EQ(11,*a);
}
-namespace {
- struct check_cref {
- int operator()(const int &) { return 10; }
- int operator()(int &) { return 11; }
- };
-}
+struct check_cref {
+ int operator()(const int&) { return 10; }
+ int operator()(int&) { return 11; }
+};
TEST(optionalm,bind_constness) {
check_cref checker;
@@ -314,7 +331,7 @@ TEST(optionalm,and_operator) {
auto zb=false & b;
EXPECT_EQ(typeid(zb),typeid(b));
EXPECT_FALSE((bool)zb);
-
+
auto b3=b & 3;
EXPECT_EQ(typeid(b3),typeid(optional<int>));
EXPECT_TRUE((bool)b3);
diff --git a/tests/test_parameters.cpp b/tests/unit/test_parameters.cpp
similarity index 94%
rename from tests/test_parameters.cpp
rename to tests/unit/test_parameters.cpp
index 8a4b0d200785036d727fda4cae698a7950fbf040..a573e07cd9fce4581d2bd6b026e1499f94336d1d 100644
--- a/tests/test_parameters.cpp
+++ b/tests/unit/test_parameters.cpp
@@ -1,7 +1,7 @@
#include <fstream>
#include "gtest.h"
-#include "util.hpp"
+#include "../test_util.hpp"
#include "../src/parameter_list.hpp"
@@ -28,7 +28,7 @@ TEST(parameters, setting)
EXPECT_FALSE(list.add_parameter({"b", -3.0}));
EXPECT_EQ(list.num_parameters(), 2);
- auto &parms = list.parameters();
+ auto& parms = list.parameters();
EXPECT_EQ(parms[0].name, "a");
EXPECT_EQ(parms[0].value, 0.12);
EXPECT_EQ(parms[0].range.min, 0);
diff --git a/tests/test_point.cpp b/tests/unit/test_point.cpp
similarity index 100%
rename from tests/test_point.cpp
rename to tests/unit/test_point.cpp
diff --git a/tests/test_probe.cpp b/tests/unit/test_probe.cpp
similarity index 77%
rename from tests/test_probe.cpp
rename to tests/unit/test_probe.cpp
index 52ac4e7f62171bc556df5b0ba4e25f05b3304fa5..5af744d2cc5f7f9b8a403b991b22d260c9ec74d4 100644
--- a/tests/test_probe.cpp
+++ b/tests/unit/test_probe.cpp
@@ -1,5 +1,6 @@
#include "gtest.h"
+#include "common_types.hpp"
#include "cell.hpp"
#include "fvm_cell.hpp"
@@ -12,13 +13,13 @@ TEST(probe, instantiation)
segment_location loc1{0, 0};
segment_location loc2{1, 0.6};
- auto p1 = c1.add_probe(loc1, probeKind::membrane_voltage);
- auto p2 = c1.add_probe(loc2, probeKind::membrane_current);
+ auto p1 = c1.add_probe({loc1, probeKind::membrane_voltage});
+ auto p2 = c1.add_probe({loc2, probeKind::membrane_current});
// expect locally provided probe ids to be numbered sequentially from zero.
-
- EXPECT_EQ(0, p1);
- EXPECT_EQ(1, p2);
+
+ EXPECT_EQ(0u, p1);
+ EXPECT_EQ(1u, p2);
// expect the probes() return to be a collection with these two probes.
@@ -48,14 +49,14 @@ TEST(probe, fvm_cell)
segment_location loc1{1, 1};
segment_location loc2{1, 0.5};
- auto pv0 = bs.add_probe(loc0, probeKind::membrane_voltage);
- auto pv1 = bs.add_probe(loc1, probeKind::membrane_voltage);
- auto pi2 = bs.add_probe(loc2, probeKind::membrane_current);
-
+ auto pv0 = bs.add_probe({loc0, probeKind::membrane_voltage});
+ auto pv1 = bs.add_probe({loc1, probeKind::membrane_voltage});
+ auto pi2 = bs.add_probe({loc2, probeKind::membrane_current});
+
i_clamp stim(0, 100, 0.3);
bs.add_stimulus({1, 1}, stim);
- fvm::fvm_cell<double, int> lcell(bs);
+ fvm::fvm_cell<double, cell_local_size_type> lcell(bs);
lcell.setup_matrix(0.01);
lcell.initialize();
diff --git a/tests/test_segment.cpp b/tests/unit/test_segment.cpp
similarity index 100%
rename from tests/test_segment.cpp
rename to tests/unit/test_segment.cpp
diff --git a/tests/unit/test_spike_store.cpp b/tests/unit/test_spike_store.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd075f878bfccc7a7a07887e6736fcb2768e0237
--- /dev/null
+++ b/tests/unit/test_spike_store.cpp
@@ -0,0 +1,84 @@
+#include "gtest.h"
+
+#include <thread_private_spike_store.hpp>
+
+TEST(spike_store, insert)
+{
+ using store_type = nest::mc::thread_private_spike_store<float>;
+
+ store_type store;
+
+ // insert 3 spike events and check that they were inserted correctly
+ store.insert({
+ {{0,0}, 0.0f},
+ {{1,2}, 0.5f},
+ {{2,4}, 1.0f}
+ });
+
+ {
+ EXPECT_EQ(store.get().size(), 3u);
+ auto i = 0u;
+ for (auto& spike : store.get()) {
+ EXPECT_EQ(spike.source.gid, i);
+ EXPECT_EQ(spike.source.index, 2*i);
+ EXPECT_EQ(spike.time, float(i)/2.f);
+ ++i;
+ }
+ }
+
+ // insert another 3 events, then check that they were appended to the
+ // original three events correctly
+ store.insert({
+ {{3,6}, 1.5f},
+ {{4,8}, 2.0f},
+ {{5,10}, 2.5f}
+ });
+
+ {
+ EXPECT_EQ(store.get().size(), 6u);
+ auto i = 0u;
+ for (auto& spike : store.get()) {
+ EXPECT_EQ(spike.source.gid, i);
+ EXPECT_EQ(spike.source.index, 2*i);
+ EXPECT_EQ(spike.time, float(i)/2.f);
+ ++i;
+ }
+ }
+}
+
+TEST(spike_store, clear)
+{
+ using store_type = nest::mc::thread_private_spike_store<float>;
+
+ store_type store;
+
+ // insert 3 spike events
+ store.insert({
+ {{0,0}, 0.0f}, {{1,2}, 0.5f}, {{2,4}, 1.0f}
+ });
+ EXPECT_EQ(store.get().size(), 3u);
+ store.clear();
+ EXPECT_EQ(store.get().size(), 0u);
+}
+
+TEST(spike_store, gather)
+{
+ using store_type = nest::mc::thread_private_spike_store<float>;
+
+ store_type store;
+
+ auto spikes = std::vector<store_type::spike_type>
+ { {{0,0}, 0.0f}, {{1,2}, 0.5f}, {{2,4}, 1.0f} };
+
+ store.insert(spikes);
+ auto gathered_spikes = store.gather();
+
+ EXPECT_EQ(gathered_spikes.size(), spikes.size());
+
+ for(auto i=0u; i<spikes.size(); ++i) {
+ EXPECT_EQ(spikes[i].source.gid, gathered_spikes[i].source.gid);
+ EXPECT_EQ(spikes[i].source.index, gathered_spikes[i].source.index);
+ EXPECT_EQ(spikes[i].time, gathered_spikes[i].time);
+ }
+}
+
diff --git a/tests/test_spikes.cpp b/tests/unit/test_spikes.cpp
similarity index 95%
rename from tests/test_spikes.cpp
rename to tests/unit/test_spikes.cpp
index b16f89cda31768f43b8509b9f9cdde375ed9bc25..3b4862029208c1f6c5a657d83dd982f510852cd6 100644
--- a/tests/test_spikes.cpp
+++ b/tests/unit/test_spikes.cpp
@@ -1,7 +1,7 @@
#include "gtest.h"
-#include <communication/spike.hpp>
-#include <communication/spike_source.hpp>
+#include <spike.hpp>
+#include <spike_source.hpp>
struct cell_proxy {
double voltage(nest::mc::segment_location loc) const {
diff --git a/tests/test_stimulus.cpp b/tests/unit/test_stimulus.cpp
similarity index 100%
rename from tests/test_stimulus.cpp
rename to tests/unit/test_stimulus.cpp
diff --git a/tests/test_swcio.cpp b/tests/unit/test_swcio.cpp
similarity index 98%
rename from tests/test_swcio.cpp
rename to tests/unit/test_swcio.cpp
index 32eaff75b967ccf9c5f8b7ef7f06c8b3d8b10f3e..c927224ffaef7e6e9d026a53dd81377b061a24eb 100644
--- a/tests/test_swcio.cpp
+++ b/tests/unit/test_swcio.cpp
@@ -502,7 +502,7 @@ TEST(swc_io, cell_construction)
cell cell = io::swc_read_cell(is);
EXPECT_TRUE(cell.has_soma());
- EXPECT_EQ(4, cell.num_segments());
+ EXPECT_EQ(4u, 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());
@@ -515,13 +515,13 @@ TEST(swc_io, cell_construction)
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) {
+ for (auto i = 1u; 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());
+ EXPECT_EQ(1u, cell.cable(1)->num_sub_segments());
+ EXPECT_EQ(1u, cell.cable(2)->num_sub_segments());
+ EXPECT_EQ(2u, cell.cable(3)->num_sub_segments());
// Check the radii
@@ -563,7 +563,7 @@ TEST(swc_parser, from_file_ball_and_stick)
auto cell = nest::mc::io::swc_read_cell(fid);
// verify that the correct number of nodes was read
- EXPECT_EQ(cell.num_segments(), 2);
+ EXPECT_EQ(cell.num_segments(), 2u);
EXPECT_EQ(cell.num_compartments(), 2u);
// make an equivalent cell via C++ interface
diff --git a/tests/unit/test_synapses.cpp b/tests/unit/test_synapses.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8cd9de04f57625ec7416e5a571ab24447c1a6cae
--- /dev/null
+++ b/tests/unit/test_synapses.cpp
@@ -0,0 +1,144 @@
+#include "gtest.h"
+#include "../test_util.hpp"
+
+#include <cell.hpp>
+#include <fvm_cell.hpp>
+#include <mechanisms/expsyn.hpp>
+#include <mechanisms/exp2syn.hpp>
+
+// compares results with those generated by nrn/ball_and_stick.py
+TEST(synapses, add_to_cell)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ // nest::mc::mechanisms::setup_mechanism_helpers();
+
+ // Soma with diameter 12.6157 um and HH channel
+ auto soma = cell.add_soma(12.6157/2.0);
+ soma->add_mechanism(hh_parameters());
+
+ parameter_list exp_default("expsyn");
+ parameter_list exp2_default("exp2syn");
+
+ cell.add_synapse({0, 0.1}, exp_default);
+ cell.add_synapse({1, 0.2}, exp2_default);
+ cell.add_synapse({0, 0.3}, exp_default);
+
+ EXPECT_EQ(3u, cell.synapses().size());
+ const auto& syns = cell.synapses();
+
+ EXPECT_EQ(syns[0].location.segment, 0u);
+ EXPECT_EQ(syns[0].location.position, 0.1);
+ EXPECT_EQ(syns[0].mechanism.name(), "expsyn");
+
+ EXPECT_EQ(syns[1].location.segment, 1u);
+ EXPECT_EQ(syns[1].location.position, 0.2);
+ EXPECT_EQ(syns[1].mechanism.name(), "exp2syn");
+
+ EXPECT_EQ(syns[2].location.segment, 0u);
+ EXPECT_EQ(syns[2].location.position, 0.3);
+ EXPECT_EQ(syns[2].mechanism.name(), "expsyn");
+}
+
+// compares results with those generated by nrn/ball_and_stick.py
+TEST(synapses, expsyn_basic_state)
+{
+ using namespace nest::mc;
+
+ using synapse_type = mechanisms::expsyn::mechanism_expsyn<double, int>;
+ auto num_syn = 4;
+
+ synapse_type::index_type indexes(num_syn);
+ synapse_type::vector_type voltage(num_syn, -65.0);
+ synapse_type::vector_type current(num_syn, 1.0);
+ auto mech = mechanisms::make_mechanism<synapse_type>( voltage, current, indexes );
+
+ auto ptr = dynamic_cast<synapse_type*>(mech.get());
+
+ // parameters initialized to default values
+ for(auto e : ptr->e) {
+ EXPECT_EQ(e, 0.);
+ }
+ for(auto tau : ptr->tau) {
+ EXPECT_EQ(tau, 2.0);
+ }
+
+ // current and voltage vectors correctly hooked up
+ for(auto v : ptr->vec_v_) {
+ EXPECT_EQ(v, -65.);
+ }
+ for(auto i : ptr->vec_i_) {
+ EXPECT_EQ(i, 1.0);
+ }
+
+ // should be initialized to NaN
+ for(auto g : ptr->g) {
+ EXPECT_NE(g, g);
+ }
+
+ // initialize state then check g has been set to zero
+ ptr->nrn_init();
+ for(auto g : ptr->g) {
+ EXPECT_EQ(g, 0.);
+ }
+
+ // call net_receive on two of the synapses
+ ptr->net_receive(1, 3.14);
+ ptr->net_receive(3, 1.04);
+ EXPECT_EQ(ptr->g[1], 3.14);
+ EXPECT_EQ(ptr->g[3], 1.04);
+}
+
+TEST(synapses, exp2syn_basic_state)
+{
+ using namespace nest::mc;
+
+ using synapse_type = mechanisms::exp2syn::mechanism_exp2syn<double, int>;
+ auto num_syn = 4;
+
+ synapse_type::index_type indexes(num_syn);
+ synapse_type::vector_type voltage(num_syn, -65.0);
+ synapse_type::vector_type current(num_syn, 1.0);
+ auto mech = mechanisms::make_mechanism<synapse_type>( voltage, current, indexes );
+
+ auto ptr = dynamic_cast<synapse_type*>(mech.get());
+
+ // parameters initialized to default values
+ for(auto e : ptr->e) {
+ EXPECT_EQ(e, 0.);
+ }
+ for(auto tau1: ptr->tau1) {
+ EXPECT_EQ(tau1, 0.5);
+ }
+ for(auto tau2: ptr->tau2) {
+ EXPECT_EQ(tau2, 2.0);
+ }
+
+ // should be initialized to NaN
+ for(auto factor: ptr->factor) {
+ EXPECT_NE(factor, factor);
+ }
+
+ // initialize state then check factor has sane (positive) value
+ // and A and B are zero
+ ptr->nrn_init();
+ for(auto factor: ptr->factor) {
+ EXPECT_GT(factor, 0.);
+ }
+ for(auto A: ptr->A) {
+ EXPECT_EQ(A, 0.);
+ }
+ for(auto B: ptr->B) {
+ EXPECT_EQ(B, 0.);
+ }
+
+ // call net_receive on two of the synapses
+ ptr->net_receive(1, 3.14);
+ ptr->net_receive(3, 1.04);
+
+ EXPECT_NEAR(ptr->A[1], ptr->factor[1]*3.14, 1e-6);
+ EXPECT_NEAR(ptr->B[3], ptr->factor[3]*1.04, 1e-6);
+}
diff --git a/tests/test_tree.cpp b/tests/unit/test_tree.cpp
similarity index 84%
rename from tests/test_tree.cpp
rename to tests/unit/test_tree.cpp
index ca9c6d73b1d9e4bde0bb90fc9061ff34f69f0503..4c69ce400d1c6146f7f704b2317a0189b8a41a53 100644
--- a/tests/test_tree.cpp
+++ b/tests/unit/test_tree.cpp
@@ -17,20 +17,24 @@ using json = nlohmann::json;
using range = memory::Range;
using namespace nest::mc;
+using int_type = cell_tree::int_type;
+
TEST(cell_tree, from_parent_index) {
+ auto no_parent = cell_tree::no_parent;
+
// tree with single branch corresponding to the root node
// this is equivalent to a single compartment model
// CASE 1 : single root node in parent_index
{
- std::vector<int> parent_index = {0};
+ std::vector<int_type> parent_index = {0};
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 1u);
EXPECT_EQ(tree.num_children(0), 0u);
}
// CASE 2 : empty parent_index
{
- std::vector<int> parent_index;
+ std::vector<int_type> parent_index;
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 1u);
EXPECT_EQ(tree.num_children(0), 0u);
@@ -46,7 +50,7 @@ TEST(cell_tree, from_parent_index) {
// /
// 3
//
- std::vector<int> parent_index =
+ std::vector<int_type> parent_index =
{0, 0, 1, 2, 0, 4};
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 3u);
@@ -66,7 +70,7 @@ TEST(cell_tree, from_parent_index) {
// / \.
// 3 8
//
- std::vector<int> parent_index =
+ std::vector<int_type> parent_index =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8};
cell_tree tree(parent_index);
@@ -79,10 +83,10 @@ TEST(cell_tree, from_parent_index) {
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));
+ EXPECT_EQ(no_parent, tree.parent(0));
+ EXPECT_EQ(0u, tree.parent(1));
+ EXPECT_EQ(0u, tree.parent(2));
+ EXPECT_EQ(0u, tree.parent(3));
}
{
//
@@ -100,7 +104,7 @@ TEST(cell_tree, from_parent_index) {
// \.
// 13
//
- std::vector<int> parent_index =
+ std::vector<int_type> parent_index =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12};
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 6u);
@@ -115,12 +119,12 @@ TEST(cell_tree, from_parent_index) {
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));
+ EXPECT_EQ(no_parent, tree.parent(0));
+ EXPECT_EQ(0u, tree.parent(1));
+ EXPECT_EQ(0u, tree.parent(2));
+ EXPECT_EQ(0u, tree.parent(3));
+ EXPECT_EQ(3u, tree.parent(4));
+ EXPECT_EQ(3u, tree.parent(5));
}
{
//
@@ -129,7 +133,7 @@ TEST(cell_tree, from_parent_index) {
// 1
// / \.
// 2 3
- std::vector<int> parent_index = {0,0,1,1};
+ std::vector<int_type> parent_index = {0,0,1,1};
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 4u);
@@ -146,7 +150,7 @@ TEST(cell_tree, from_parent_index) {
// 1 4 5
// / \.
// 2 3
- std::vector<int> parent_index = {0,0,1,1,0,0};
+ std::vector<int_type> parent_index = {0,0,1,1,0,0};
cell_tree tree(parent_index);
EXPECT_EQ(tree.num_segments(), 6u);
@@ -158,12 +162,13 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(4), 0u);
// Check children
- EXPECT_EQ(1, tree.children(0)[0]);
- EXPECT_EQ(4, tree.children(0)[1]);
- EXPECT_EQ(5, tree.children(0)[2]);
- EXPECT_EQ(2, tree.children(1)[0]);
- EXPECT_EQ(3, tree.children(1)[1]);
+ EXPECT_EQ(1u, tree.children(0)[0]);
+ EXPECT_EQ(4u, tree.children(0)[1]);
+ EXPECT_EQ(5u, tree.children(0)[2]);
+ EXPECT_EQ(2u, tree.children(1)[0]);
+ EXPECT_EQ(3u, tree.children(1)[1]);
}
+ /* FIXME
{
// 0
// / \.
@@ -185,6 +190,7 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(5), 0u);
EXPECT_EQ(tree.num_children(6), 0u);
}
+ */
}
TEST(tree, change_root) {
@@ -196,8 +202,8 @@ TEST(tree, change_root) {
// 1 2 -> 1
// |
// 2
- std::vector<int> parent_index = {0,0,0};
- tree t(parent_index);
+ std::vector<int_type> parent_index = {0,0,0};
+ tree<int_type> t(parent_index);
t.change_root(1);
EXPECT_EQ(t.num_nodes(), 3u);
@@ -214,8 +220,8 @@ TEST(tree, change_root) {
// 1 2 -> 1 2 3
// / \ |
// 3 4 4
- std::vector<int> parent_index = {0,0,0,1,1};
- tree t(parent_index);
+ std::vector<int_type> parent_index = {0,0,0,1,1};
+ tree<int_type> t(parent_index);
t.change_root(1u);
EXPECT_EQ(t.num_nodes(), 5u);
@@ -238,8 +244,8 @@ TEST(tree, change_root) {
// 3 4 3 4 6
// / \.
// 5 6
- std::vector<int> parent_index = {0,0,0,1,1,4,4};
- tree t(parent_index);
+ std::vector<int_type> parent_index = {0,0,0,1,1,4,4};
+ tree<int_type> t(parent_index);
t.change_root(1);
@@ -256,6 +262,8 @@ TEST(tree, change_root) {
}
TEST(cell_tree, balance) {
+ auto no_parent = cell_tree::no_parent;
+
{
// a cell with the following structure
// will balance around 1
@@ -266,13 +274,13 @@ TEST(cell_tree, balance) {
// 3 4 3 4 6
// / \.
// 5 6
- std::vector<int> parent_index = {0,0,0,1,1,4,4};
+ std::vector<int_type> parent_index = {0,0,0,1,1,4,4};
cell_tree t(parent_index);
t.balance();
// the soma (original root) has moved to 5 in the new tree
- EXPECT_EQ(t.soma(), 5);
+ EXPECT_EQ(t.soma(), 5u);
EXPECT_EQ(t.num_segments(), 7u);
EXPECT_EQ(t.num_children(0),3u);
@@ -282,13 +290,13 @@ TEST(cell_tree, balance) {
EXPECT_EQ(t.num_children(4),0u);
EXPECT_EQ(t.num_children(5),1u);
EXPECT_EQ(t.num_children(6),0u);
- EXPECT_EQ(t.parent(0),-1);
- EXPECT_EQ(t.parent(1), 0);
- EXPECT_EQ(t.parent(2), 0);
- EXPECT_EQ(t.parent(3), 0);
- EXPECT_EQ(t.parent(4), 2);
- EXPECT_EQ(t.parent(5), 2);
- EXPECT_EQ(t.parent(6), 5);
+ EXPECT_EQ(t.parent(0), no_parent);
+ EXPECT_EQ(t.parent(1), 0u);
+ EXPECT_EQ(t.parent(2), 0u);
+ EXPECT_EQ(t.parent(3), 0u);
+ EXPECT_EQ(t.parent(4), 2u);
+ EXPECT_EQ(t.parent(5), 2u);
+ EXPECT_EQ(t.parent(6), 5u);
//t.to_graphviz("cell.dot");
}
@@ -305,7 +313,7 @@ TEST(cell_tree, json_load)
std::ifstream(path) >> cell_data;
for(auto c : range(0,cell_data.size())) {
- std::vector<int> parent_index = cell_data[c]["parent_index"];
+ std::vector<int_type> parent_index = cell_data[c]["parent_index"];
cell_tree tree(parent_index);
//tree.to_graphviz("cell" + std::to_string(c) + ".dot");
}
@@ -313,6 +321,7 @@ TEST(cell_tree, json_load)
TEST(tree, make_parent_index)
{
+ /*
// just the soma
{
std::vector<int> parent_index = {0};
@@ -378,4 +387,5 @@ TEST(tree, make_parent_index)
auto new_parent_index = make_parent_index(t, counts);
EXPECT_EQ(parent_index, new_parent_index);
}
+ */
}
diff --git a/tests/test_uninitialized.cpp b/tests/unit/test_uninitialized.cpp
similarity index 80%
rename from tests/test_uninitialized.cpp
rename to tests/unit/test_uninitialized.cpp
index dcc9e47ce6a7fa266376e054e8f84a04d8ec818d..1f04510d7df6af054eaec854d7e0c954ad5e73c8 100644
--- a/tests/test_uninitialized.cpp
+++ b/tests/unit/test_uninitialized.cpp
@@ -7,16 +7,16 @@ using namespace nest::mc::util;
namespace {
struct count_ops {
count_ops() {}
- count_ops(const count_ops &n) { ++copy_ctor_count; }
- count_ops(count_ops &&n) { ++move_ctor_count; }
+ count_ops(const count_ops& n) { ++copy_ctor_count; }
+ count_ops(count_ops&& n) { ++move_ctor_count; }
- count_ops &operator=(const count_ops &n) { ++copy_assign_count; return *this; }
- count_ops &operator=(count_ops &&n) { ++move_assign_count; return *this; }
+ count_ops& operator=(const count_ops& n) { ++copy_assign_count; return *this; }
+ count_ops& operator=(count_ops&& n) { ++move_assign_count; return *this; }
static int copy_ctor_count,copy_assign_count;
static int move_ctor_count,move_assign_count;
static void reset_counts() {
- copy_ctor_count=copy_assign_count=0;
+ copy_ctor_count=copy_assign_count=0;
move_ctor_count=move_assign_count=0;
}
};
@@ -53,11 +53,11 @@ TEST(uninitialized,ctor) {
namespace {
struct nocopy {
nocopy() {}
- nocopy(const nocopy &n) = delete;
- nocopy(nocopy &&n) { ++move_ctor_count; }
+ nocopy(const nocopy& n) = delete;
+ nocopy(nocopy&& n) { ++move_ctor_count; }
- nocopy &operator=(const nocopy &n) = delete;
- nocopy &operator=(nocopy &&n) { ++move_assign_count; return *this; }
+ nocopy& operator=(const nocopy& n) = delete;
+ nocopy& operator=(nocopy&& n) { ++move_assign_count; return *this; }
static int move_ctor_count,move_assign_count;
static void reset_counts() { move_ctor_count=move_assign_count=0; }
@@ -85,11 +85,11 @@ TEST(uninitialized,ctor_nocopy) {
namespace {
struct nomove {
nomove() {}
- nomove(const nomove &n) { ++copy_ctor_count; }
- nomove(nomove &&n) = delete;
+ nomove(const nomove& n) { ++copy_ctor_count; }
+ nomove(nomove&& n) = delete;
- nomove &operator=(const nomove &n) { ++copy_assign_count; return *this; }
- nomove &operator=(nomove &&n) = delete;
+ nomove& operator=(const nomove& n) { ++copy_assign_count; return *this; }
+ nomove& operator=(nomove&& n) = delete;
static int copy_ctor_count,copy_assign_count;
static void reset_counts() { copy_ctor_count=copy_assign_count=0; }
@@ -129,7 +129,7 @@ TEST(uninitialized,void) {
}
TEST(uninitialized,ref) {
- uninitialized<int &> x,y;
+ uninitialized<int&> x,y;
int a;
x.construct(a);
@@ -151,8 +151,8 @@ namespace {
mutable int op_count=0;
mutable int const_op_count=0;
- int operator()(const int &a) const { ++const_op_count; return a+1; }
- int operator()(int &a) const { ++op_count; return ++a; }
+ int operator()(const int& a) const { ++const_op_count; return a+1; }
+ int operator()(int& a) const { ++op_count; return ++a; }
};
}
@@ -165,14 +165,14 @@ TEST(uninitialized,apply) {
EXPECT_EQ(11,ua.cref());
EXPECT_EQ(11,r);
- uninitialized<int &> ub;
+ uninitialized<int&> ub;
ub.construct(ua.ref());
r=ub.apply(A);
EXPECT_EQ(12,ua.cref());
EXPECT_EQ(12,r);
- uninitialized<const int &> uc;
+ uninitialized<const int&> uc;
uc.construct(ua.ref());
r=uc.apply(A);
diff --git a/tests/validate.cpp b/tests/validate.cpp
deleted file mode 100644
index d814fca1512dd79b22763a4f9a769f2984ca5269..0000000000000000000000000000000000000000
--- a/tests/validate.cpp
+++ /dev/null
@@ -1,7 +0,0 @@
-#include "gtest.h"
-
-int main(int argc, char **argv) {
- ::testing::InitGoogleTest(&argc, argv);
- return RUN_ALL_TESTS();
-}
-
diff --git a/tests/validation/CMakeLists.txt b/tests/validation/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c6c72fa04efab7394190288ba3face603d141faf
--- /dev/null
+++ b/tests/validation/CMakeLists.txt
@@ -0,0 +1,37 @@
+set(VALIDATION_SOURCES
+ # unit tests
+ validate_ball_and_stick.cpp
+ validate_soma.cpp
+ validate_synapses.cpp
+
+ # support code
+ validation_data.cpp
+
+ # unit test driver
+ validate.cpp
+)
+
+add_definitions("-DDATADIR=\"${CMAKE_SOURCE_DIR}/data\"")
+
+add_executable(validate.exe ${VALIDATION_SOURCES} ${HEADERS})
+
+set(TARGETS validate.exe)
+
+foreach(target ${TARGETS})
+ target_link_libraries(${target} LINK_PUBLIC cellalgo gtest)
+
+ if(WITH_TBB)
+ target_link_libraries(${target} LINK_PUBLIC ${TBB_LIBRARIES})
+ endif()
+
+ if(WITH_MPI)
+ target_link_libraries(${target} LINK_PUBLIC ${MPI_C_LIBRARIES})
+ set_property(TARGET ${target} APPEND_STRING PROPERTY LINK_FLAGS "${MPI_C_LINK_FLAGS}")
+ endif()
+
+ set_target_properties(${target}
+ PROPERTIES
+ RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests"
+ )
+endforeach()
+
diff --git a/tests/make_image.sh b/tests/validation/make_image.sh
old mode 100755
new mode 100644
similarity index 100%
rename from tests/make_image.sh
rename to tests/validation/make_image.sh
diff --git a/tests/plot.py b/tests/validation/plot.py
similarity index 100%
rename from tests/plot.py
rename to tests/validation/plot.py
diff --git a/tests/validation/validate.cpp b/tests/validation/validate.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b89a7b54223d7c2342c9882749c28a886782b316
--- /dev/null
+++ b/tests/validation/validate.cpp
@@ -0,0 +1,31 @@
+#include <cstring>
+#include <iostream>
+#include <fstream>
+#include <numeric>
+#include <vector>
+
+#include "gtest.h"
+#include "validation_data.hpp"
+
+int usage(const char* argv0) {
+ std::cerr << "usage: " << argv0 << " [-p|--path validation_data_directory]\n";
+ return 1;
+}
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+
+ if (argv[1] && (!std::strcmp(argv[1], "-p") || !std::strcmp(argv[1], "--path"))) {
+ if (argv[2]) {
+ testing::g_validation_data.set_path(argv[2]);
+ }
+ else {
+ return usage(argv[0]);
+ }
+ }
+ else if (argv[1]) {
+ return usage(argv[0]);
+ }
+
+ return RUN_ALL_TESTS();
+}
diff --git a/tests/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp
similarity index 96%
rename from tests/validate_ball_and_stick.cpp
rename to tests/validation/validate_ball_and_stick.cpp
index 5b509525b053fe389908b150b87333a27f5241ec..62d06b4dac85885f935dd88855223551cc33fcce 100644
--- a/tests/validate_ball_and_stick.cpp
+++ b/tests/validation/validate_ball_and_stick.cpp
@@ -1,12 +1,14 @@
#include <fstream>
#include <json/src/json.hpp>
-#include "gtest.h"
-#include "util.hpp"
-
+#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
+#include "gtest.h"
+#include "../test_util.hpp"
+#include "validation_data.hpp"
+
// compares results with those generated by nrn/ball_and_stick.py
TEST(ball_and_stick, neuron_baseline)
{
@@ -15,9 +17,6 @@ TEST(ball_and_stick, neuron_baseline)
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());
@@ -32,7 +31,7 @@ TEST(ball_and_stick, neuron_baseline)
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");
+ auto cell_data = testing::g_validation_data.load("ball_and_stick.json");
EXPECT_TRUE(cell_data.size()>0);
if(cell_data.size()==0) return;
@@ -91,7 +90,7 @@ TEST(ball_and_stick, neuron_baseline)
std::vector<std::vector<double>> v(3);
// make the lowered finite volume cell
- fvm::fvm_cell<double, int> model(cell);
+ fvm::fvm_cell<double, cell_local_size_type> model(cell);
auto graph = cell.model();
// set initial conditions
@@ -164,9 +163,6 @@ TEST(ball_and_3stick, neuron_baseline)
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());
@@ -187,7 +183,7 @@ TEST(ball_and_3stick, neuron_baseline)
cell.add_stimulus({3,1}, {40., 10.,-0.2});
// load data from file
- auto cell_data = testing::load_spike_data("../nrn/ball_and_3stick.json");
+ auto cell_data = testing::g_validation_data.load("ball_and_3stick.json");
EXPECT_TRUE(cell_data.size()>0);
if(cell_data.size()==0) return;
@@ -249,7 +245,7 @@ TEST(ball_and_3stick, neuron_baseline)
std::vector<std::vector<double>> v(3);
// make the lowered finite volume cell
- fvm::fvm_cell<double, int> model(cell);
+ fvm::fvm_cell<double, cell_local_size_type> model(cell);
auto graph = cell.model();
// set initial conditions
diff --git a/tests/validate_soma.cpp b/tests/validation/validate_soma.cpp
similarity index 93%
rename from tests/validate_soma.cpp
rename to tests/validation/validate_soma.cpp
index 2353cb483967fb74757dab81ffcf96c68f3c649a..24a0c1c59f0db24bc6f9c7832909f9f52af7beb7 100644
--- a/tests/validate_soma.cpp
+++ b/tests/validation/validate_soma.cpp
@@ -1,12 +1,14 @@
#include <fstream>
#include <json/src/json.hpp>
-#include "gtest.h"
-#include "util.hpp"
-
+#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
+#include "gtest.h"
+#include "../test_util.hpp"
+#include "validation_data.hpp"
+
// compares results with those generated by nrn/soma.py
// single compartment model with HH channels
TEST(soma, neuron_baseline)
@@ -16,9 +18,6 @@ TEST(soma, neuron_baseline)
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
@@ -28,10 +27,10 @@ TEST(soma, neuron_baseline)
cell.add_stimulus({0,0.5}, {10., 100., 0.1});
// make the lowered finite volume cell
- fvm::fvm_cell<double, int> model(cell);
+ fvm::fvm_cell<double, cell_local_size_type> model(cell);
// load data from file
- auto cell_data = testing::load_spike_data("../nrn/soma.json");
+ auto cell_data = testing::g_validation_data.load("soma.json");
EXPECT_TRUE(cell_data.size()>0);
if(cell_data.size()==0) return;
@@ -84,9 +83,6 @@ TEST(soma, convergence)
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
@@ -96,7 +92,7 @@ TEST(soma, convergence)
cell.add_stimulus({0,0.5}, {10., 100., 0.1});
// make the lowered finite volume cell
- fvm::fvm_cell<double, int> model(cell);
+ fvm::fvm_cell<double, cell_local_size_type> model(cell);
// generate baseline solution with small dt=0.0001
std::vector<double> baseline_spike_times;
diff --git a/tests/validate_synapses.cpp b/tests/validation/validate_synapses.cpp
similarity index 71%
rename from tests/validate_synapses.cpp
rename to tests/validation/validate_synapses.cpp
index 2ccde863cb4cbf16efbe75dd509e9c7c2095ddbe..876b24eb55c9cad0e08c17cef9cdf16595f5bd90 100644
--- a/tests/validate_synapses.cpp
+++ b/tests/validation/validate_synapses.cpp
@@ -1,13 +1,17 @@
#include <fstream>
+#include <iterator>
-#include "gtest.h"
-#include "util.hpp"
+#include <json/src/json.hpp>
+#include <common_types.hpp>
#include <cell.hpp>
+#include <cell_group.hpp>
#include <fvm_cell.hpp>
+#include <mechanism_interface.hpp>
-#include <json/src/json.hpp>
-
+#include "gtest.h"
+#include "../test_util.hpp"
+#include "validation_data.hpp"
// For storing the results of a simulation along with the information required
// to compare two simulations for accuracy.
@@ -43,16 +47,14 @@ struct result {
};
// compares results with those generated by nrn/simple_synapse.py
-TEST(simple_synapse, neuron_baseline)
+void run_neuron_baseline(const char* syn_type, const char* data_file)
{
using namespace nest::mc;
using namespace nlohmann;
+ using lowered_cell = fvm::fvm_cell<double, cell_local_size_type>;
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());
@@ -64,11 +66,23 @@ TEST(simple_synapse, neuron_baseline)
dendrite->mechanism("membrane").set("r_L", 100);
soma->mechanism("membrane").set("r_L", 100);
- // add stimulus
- cell.add_synapse({1, 0.5});
+ // add synapse
+ parameter_list syn_default(syn_type);
+ cell.add_synapse({1, 0.5}, syn_default);
+
+ // add probes
+ auto probe_soma = cell.add_probe({{0,0}, probeKind::membrane_voltage});
+ auto probe_dend = cell.add_probe({{1,0.5}, probeKind::membrane_voltage});
+
+ // injected spike events
+ postsynaptic_spike_event<float> synthetic_events[] = {
+ {{0u, 0u}, 10.0, 0.04},
+ {{0u, 0u}, 20.0, 0.04},
+ {{0u, 0u}, 40.0, 0.04}
+ };
// load data from file
- auto cell_data = testing::load_spike_data("../nrn/simple_synapse.json");
+ auto cell_data = testing::g_validation_data.load(data_file);
EXPECT_TRUE(cell_data.size()>0);
if(cell_data.size()==0) return;
@@ -91,34 +105,24 @@ TEST(simple_synapse, neuron_baseline)
std::vector<std::vector<double>> v(2);
// make the lowered finite volume cell
- fvm::fvm_cell<double, int> model(cell);
- auto graph = cell.model();
-
- // set initial conditions
- using memory::all;
- model.voltage()(all) = -65.;
- model.initialize(); // have to do this _after_ initial conditions are set
+ cell_group<lowered_cell> group(0, cell);
// add the 3 spike events to the queue
- model.queue().push({0u, 10.0, 0.04});
- model.queue().push({0u, 20.0, 0.04});
- model.queue().push({0u, 40.0, 0.04});
+ group.enqueue_events(synthetic_events);
// run the simulation
- auto soma_comp = nest::mc::find_compartment_index({0,0}, graph);
- auto dend_comp = nest::mc::find_compartment_index({1,0.5}, graph);
- v[0].push_back(model.voltage()[soma_comp]);
- v[1].push_back(model.voltage()[dend_comp]);
+ v[0].push_back(group.cell().probe(probe_soma));
+ v[1].push_back(group.cell().probe(probe_dend));
double t = 0.;
while(t < tfinal) {
t += dt;
- model.advance_to(t, dt);
- // save voltage at soma
- v[0].push_back(model.voltage()[soma_comp]);
- v[1].push_back(model.voltage()[dend_comp]);
+ group.advance(t, dt);
+ // save voltage at soma and dendrite
+ v[0].push_back(group.cell().probe(probe_soma));
+ v[1].push_back(group.cell().probe(probe_dend));
}
- results.push_back( {num_compartments, dt, v, measurements} );
+ results.push_back({num_compartments, dt, v, measurements});
}
// print results
@@ -158,3 +162,12 @@ TEST(simple_synapse, neuron_baseline)
}
}
+TEST(simple_synapse, expsyn_neuron_baseline) {
+ SCOPED_TRACE("expsyn");
+ run_neuron_baseline("expsyn","simple_exp_synapse.json");
+}
+
+TEST(simple_synapse, exp2syn_neuron_baseline) {
+ SCOPED_TRACE("exp2syn");
+ run_neuron_baseline("exp2syn","simple_exp2_synapse.json");
+}
diff --git a/tests/validation/validation_data.cpp b/tests/validation/validation_data.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..306adbaf38c02171f56a042f62662702d4faa3be
--- /dev/null
+++ b/tests/validation/validation_data.cpp
@@ -0,0 +1,21 @@
+#include "validation_data.hpp"
+
+#include <fstream>
+
+namespace testing {
+
+nlohmann::json data_loader::load(const std::string& name) const {
+ std::string data_path=path_+'/'+name;
+ std::ifstream fid(data_path);
+ if (!fid) {
+ throw std::runtime_error("unable to load validation data: "+data_path);
+ }
+
+ nlohmann::json data;
+ fid >> data;
+ return data;
+}
+
+data_loader g_validation_data;
+
+} // namespace testing
diff --git a/tests/validation/validation_data.hpp b/tests/validation/validation_data.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2aeb2f393d4fa28cac21e4d50174ad82d0c76289
--- /dev/null
+++ b/tests/validation/validation_data.hpp
@@ -0,0 +1,28 @@
+#pragma once
+
+#include <string>
+
+#include <json/src/json.hpp>
+
+#ifndef DATADIR
+#define DATADIR "../data"
+#endif
+
+namespace testing {
+
+class data_loader {
+public:
+ // set where to find the validation JSON files
+ void set_path(const std::string& path) { path_=path; }
+
+ // load JSON file from path plus name, throw exception
+ // if cannot be found or loaded.
+ nlohmann::json load(const std::string& name) const;
+
+private:
+ std::string path_=DATADIR "/validation";
+};
+
+extern data_loader g_validation_data;
+
+} // namespace testing