diff --git a/.travis.yml b/.travis.yml index d7adff912d..2c14c30192 100644 --- a/.travis.yml +++ b/.travis.yml @@ -41,28 +41,25 @@ jobs: include: - stage: Staticcheck python: 3.6.7 - env: xTHREADING=0 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=1 xRUN_BUILD_AND_TESTSUITE=0 CACHE_NAME=JOB # only static code analysis + env: xTHREADING=0 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=0 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=1 xRUN_BUILD_AND_TESTSUITE=0 CACHE_NAME=JOB # only static code analysis - stage: MPI-Threading-Python python: 3.6.7 - env: xTHREADING=0 xMPI=1 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # only MPI + env: xTHREADING=0 xMPI=1 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xLIBBOOST=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # only MPI - stage: MPI-Threading-Python python: 3.6.7 - env: xTHREADING=1 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # only threading + env: xTHREADING=1 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xLIBBOOST=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # only threading - stage: MPI-Threading-Python python: 3.6.7 - env: xTHREADING=1 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # Python & Threading + env: xTHREADING=1 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=1 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # Python & Threading - stage: MPI-Threading-Python python: 3.6.7 - env: xTHREADING=0 xMPI=1 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # Python & MPI + env: xTHREADING=0 xMPI=1 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=1 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # Python & MPI - stage: MPI-Threading-Python python: 3.6.7 - env: xTHREADING=0 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # only Python + env: xTHREADING=0 xMPI=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=1 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # only Python - stage: Python-Full-build python: 3.6.7 - env: xTHREADING=1 xMPI=1 xGSL=1 xLIBNEUROSIM=1 xLTDL=1 xREADLINE=1 xPYTHON=1 xMUSIC=1 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # full - - stage: Python-Full-build - python: 2.7.13 - env: xTHREADING=1 xMPI=1 xGSL=1 xLIBNEUROSIM=1 xLTDL=1 xREADLINE=1 xPYTHON=1 xMUSIC=1 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # full + env: xTHREADING=1 xMPI=1 xGSL=1 xLIBNEUROSIM=1 xLTDL=1 xREADLINE=1 xLIBBOOST=1 xPYTHON=1 xMUSIC=1 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB # full - stage: Clang7 language: cpp env: MATRIX_EVAL="CC=clang-7 && CXX=clang++-7" xRUN_BUILD_AND_TESTSUITE=1 @@ -139,7 +136,9 @@ before_install: echo "+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +" echo ${MATRIX_EVAL} eval "${MATRIX_EVAL}" - cp extras/install_music.sh extras/install_csa-libneurosim.sh $HOME + cp extras/install_csa-libneurosim.sh $HOME + cp extras/install_music.sh $HOME + cp extras/install_libboost.sh $HOME cd $HOME/build echo $PATH # Upgrade pip and setuptools diff --git a/cmake/ProcessOptions.cmake b/cmake/ProcessOptions.cmake index b82af1f180..ca0d190c79 100644 --- a/cmake/ProcessOptions.cmake +++ b/cmake/ProcessOptions.cmake @@ -456,6 +456,13 @@ function( NEST_PROCESS_WITH_OPENMP ) message( FATAL_ERROR "CMake can not find OpenMP." ) endif () endif () + + # Provide a dummy OpenMP::OpenMP_CXX if no OpenMP or if flags explicitly + # given. Needed to avoid problems where OpenMP::OpenMP_CXX is used. + if ( NOT TARGET OpenMP::OpenMP_CXX ) + add_library(OpenMP::OpenMP_CXX INTERFACE IMPORTED) + endif() + endfunction() function( NEST_PROCESS_WITH_MPI ) @@ -556,8 +563,10 @@ function( NEST_PROCESS_WITH_BOOST ) set( BOOST_ROOT "${with-boost}" ) endif () + set(Boost_USE_DEBUG_LIBS OFF) # ignore debug libs + set(Boost_USE_RELEASE_LIBS ON) # only find release libs # Needs Boost version >=1.58.0 to use Boost sorting - find_package( Boost 1.58.0 COMPONENTS unit_test_framework ) + find_package( Boost 1.58.0 ) if ( Boost_FOUND ) # export found variables to parent scope set( HAVE_BOOST ON PARENT_SCOPE ) @@ -566,6 +575,7 @@ function( NEST_PROCESS_WITH_BOOST ) set( BOOST_LIBRARIES "${Boost_LIBRARIES}" PARENT_SCOPE ) set( BOOST_INCLUDE_DIR "${Boost_INCLUDE_DIR}" PARENT_SCOPE ) set( BOOST_VERSION "${Boost_MAJOR_VERSION}.${Boost_MINOR_VERSION}.${Boost_SUBMINOR_VERSION}" PARENT_SCOPE ) + include_directories( ${Boost_INCLUDE_DIRS} ) endif () endif () endfunction() diff --git a/debian/changelog b/debian/changelog index 03a14fcdc6..84177661e2 100644 --- a/debian/changelog +++ b/debian/changelog @@ -1,3 +1,10 @@ +nest (2.20.1-0ubuntu1ppa1) unstable; urgency=low + + * NEST v2.20.1 + * Added boost to /debian files + + -- Steffen Graber Tue, 17 Mar 2020 07:00:00 +0000 + nest (2.18.0+1SNAPSHOT-0ubuntu1ppa1) unstable; urgency=low * NEST daily builds of master in 'ppa:nest-simulator/nest-nightly' diff --git a/debian/control b/debian/control index c8b5a05338..12daac21c2 100644 --- a/debian/control +++ b/debian/control @@ -3,7 +3,7 @@ Section: science Priority: optional Maintainer: Steffen Graber Build-Depends: cmake, debhelper (>= 9), autotools-dev, - software-properties-common, build-essential, autoconf, python3, + software-properties-common, build-essential, autoconf, python3, libboost-all-dev, libltdl-dev, libreadline-dev, libncurses-dev, libgsl-dev, openmpi-bin, python3-dev, libopenmpi-dev, libgomp1, libomp-dev, libibverbs-dev, libpcre3, libpcre3-dev, jq, python3-numpy, python3-pandas, python3-scipy,python3-matplotlib, python3-mpi4py, python3-setuptools, python3-tk, @@ -12,9 +12,9 @@ Standards-Version: 4.1.1 Package: nest Architecture: any -Depends: ${shlibs:Depends}, ${misc:Depends}, python3, libltdl-dev, - libreadline-dev, libncurses-dev, libgsl-dev, openmpi-bin, - python3-dev, libopenmpi-dev, libgomp1, libomp-dev, libibverbs-dev, libpcre3, libpcre3-dev, jq, +Depends: ${shlibs:Depends}, ${misc:Depends}, python3, libboost-all-dev, libltdl-dev, + libreadline-dev, libncurses-dev, libgsl-dev, + python3-dev, libgomp1, libomp-dev, libibverbs-dev, libpcre3, libpcre3-dev, jq, python3-numpy, python3-pandas, python3-scipy, python3-matplotlib, python3-mpi4py, python3-setuptools, python3-pydot, python3-tk, cython3 Description: NEST is a simulator for networks of spiking neurons. diff --git a/debian/copyright b/debian/copyright index 203dd0cbba..b5c3a546fe 100644 --- a/debian/copyright +++ b/debian/copyright @@ -1,19 +1,21 @@ Format: http://dep.debian.net/deps/dep5 Upstream-Name: nest-simulator -Source: https://github.com/nest/nest-simulator/archive/master.zip +Source: https://github.com/nest/nest-simulator/ -Files: debian/* -Copyright: 2019 The NEST Initiative -License: GPL-2 +Files: * +Copyright: 2004 The NEST Initiative +License: GPL-2+ + +License: GPL-2+ NEST is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version. -. + . NEST is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - + . You should have received a copy of the GNU General Public License along with NEST. If not, see . diff --git a/debian/rules b/debian/rules index 962fa24429..6785bd177a 100755 --- a/debian/rules +++ b/debian/rules @@ -8,9 +8,9 @@ INSTALLDIR = \/$(BUILDDIR)\/$(PACKAGEVERSION)\/$(DESTDIR) CMAKE_FLAGS = -DCMAKE_INSTALL_PREFIX:PATH=$(DESTDIR)/usr -Dwith-python=3 -DCMAKE_INSTALL_LIBDIR:PATH=lib CMAKE_FLAGS += -DCMAKE_CXX_COMPILER=g++ -DCMAKE_C_COMPILER=gcc -CMAKE_FLAGS += -Dwith-mpi=ON +CMAKE_FLAGS += -Dwith-mpi=OFF CMAKE_FLAGS += -Dwith-openmp=ON -CMAKE_FLAGS += -Dwith-gsl=ON -Dwith-ltdl=ON -Dwith-readline=ON +CMAKE_FLAGS += -Dwith-gsl=ON -Dwith-ltdl=ON -Dwith-readline=ON -Dwith-boost=ON %: LC_ALL=en_US.UTF-8 dh ${@} --with python3 --builddirectory=$(BUILDDIR) diff --git a/doc/contribute/mailing_list_guidelines.rst b/doc/contribute/mailing_list_guidelines.rst index a9848d67dc..706a2026cf 100644 --- a/doc/contribute/mailing_list_guidelines.rst +++ b/doc/contribute/mailing_list_guidelines.rst @@ -21,7 +21,7 @@ please follow the guidelines below: * the steps you took that lead to the problem. * the specific error messages you get. - * relevant system and version information (e.g., Ubuntu 18.04/ NEST 2.20 installed using Conda). + * relevant system and version information (e.g., Ubuntu 18.04/ NEST 2.20.1 installed using Conda). #. Keep topics separate. diff --git a/doc/download.rst b/doc/download.rst index 16b5aef1b3..ea17282020 100644 --- a/doc/download.rst +++ b/doc/download.rst @@ -21,7 +21,7 @@ If you use NEST for your project, don't forget to :doc:`cite NEST ` Download NEST source code -------------------------- -* Get the `source code of the latest release `_. +* Get the `source code of the latest release `_. * Follow the installation instructions for :doc:`Linux ` or :ref:`macOS `. @@ -33,6 +33,7 @@ Download NEST source code .. _download_livemedia: + Download the NEST live media for virtual machines -------------------------------------------------- @@ -40,14 +41,9 @@ Live media is available in the OVA format, and is suitable, for example, for imp If you run **Windows**, this is the option for you OR if you just want to run NEST without installing it on your computer. After downloading the virtual machine, check out the :doc:`install instructions for Live Media `. -* `NEST Live Media 2.18.0 `_ - - `Checksum 2.18.0 `_ - +* `NEST Live Media 2.20.1 `_ -* `NEST live media 2.16.0 `_ - - `Checksum 2.16.0 `_ + `Checksum 2.20.1 `_ Previous releases @@ -58,6 +54,18 @@ thus, we encourage our users to use the **most recent version of NEST**. **Older Versions of Live Media** +* `NEST Live Media 2.20.0 `_ + + `Checksum 2.20.0 `_ + +- `NEST Live Media 2.18.0 `_ + + `Checksum 2.18.0 `_ + +- `NEST live media 2.16.0 `_ + + `Checksum 2.16.0 `_ + - `NEST Live Media 2.14.0 `_ (OVA, 2.5G) `Checksum 2.14.0 `_ @@ -66,7 +74,6 @@ thus, we encourage our users to use the **most recent version of NEST**. `Checksum 2.12.0 `_ - - `Download 2.10.0 `_ (OVA, ~3.7G) diff --git a/doc/installation/index.rst b/doc/installation/index.rst index 93b35affca..03b1464061 100644 --- a/doc/installation/index.rst +++ b/doc/installation/index.rst @@ -241,7 +241,7 @@ these instructions.** .. tab:: macOS - For further options on installing NEST on macOS, see :ref:`mac_manual` for Macs. + For further options on installing NEST on macOS, see :doc:`mac_install`. .. tab:: HPC systems @@ -266,4 +266,3 @@ these instructions.** Installation instructions for NEST 2.10 and earlier are provided :doc:`here `, but we strongly encourage all our users to stay up-to-date with most recent version of NEST. We cannot support out-dated versions. - diff --git a/doc/tutorials/topology/hill_tononi_Vp.ipynb b/doc/tutorials/topology/hill_tononi_Vp.ipynb index c69e6f2175..8d26112dc1 100644 --- a/doc/tutorials/topology/hill_tononi_Vp.ipynb +++ b/doc/tutorials/topology/hill_tononi_Vp.ipynb @@ -26,7 +26,7 @@ }, "outputs": [], "source": [ - "import pylab\nSHOW_FIGURES = False\n\nif not SHOW_FIGURES:\n pylab_show = pylab.show\n\n def nop(s=None):\n pass\n\n pylab.show = nop\nelse:\n pylab.ion()\n\n# ! Introduction\n# !=============\n# ! This tutorial gives a brief introduction to the ConnPlotter\n# ! toolbox. It is by no means complete.\n\n# ! Load pynest\nimport nest\n\n# ! Load NEST Topoplogy module (NEST 2.2)\nimport nest.topology as topo\n\n# ! Make sure we start with a clean slate, even if we re-run the script\n# ! in the same Python session.\nnest.ResetKernel()\n\n# ! Import math, we need Pi\nimport math\n\n# ! Configurable Parameters\n# ! =======================\n# !\n# ! Here we define those parameters that we take to be\n# ! configurable. The choice of configurable parameters is obviously\n# ! arbitrary, and in practice one would have far more configurable\n# ! parameters. We restrict ourselves to:\n# !\n# ! - Network size in neurons ``N``, each layer is ``N x N``.\n# ! - Network size in subtended visual angle ``visSize``, in degree.\n# ! - Temporal frequency of drifting grating input ``f_dg``, in Hz.\n# ! - Spatial wavelength and direction of drifting grating input,\n# ! ``lambda_dg`` and ``phi_dg``, in degree/radian.\n# ! - Background firing rate of retinal nodes and modulation amplitude,\n# ! ``retDC`` and ``retAC``, in Hz.\n# ! - Simulation duration ``simtime``; actual simulation is split into\n# ! intervals of ``sim_interval`` length, so that the network state\n# ! can be visualized in those intervals. Times are in ms.\nParams = {'N': 40,\n 'visSize': 8.0,\n 'f_dg': 2.0,\n 'lambda_dg': 2.0,\n 'phi_dg': 0.0,\n 'retDC': 30.0,\n 'retAC': 30.0,\n 'simtime': 100.0,\n 'sim_interval': 5.0\n }\n\n# ! Neuron Models\n# ! =============\n# !\n# ! We declare models in two steps:\n# !\n# ! 1. We define a dictionary specifying the NEST neuron model to use\n# ! as well as the parameters for that model.\n# ! #. We create three copies of this dictionary with parameters\n# ! adjusted to the three model variants specified in Table~2 of\n# ! Hill & Tononi (2005) (cortical excitatory, cortical inhibitory,\n# ! thalamic)\n# !\n# ! In addition, we declare the models for the stimulation and\n# ! recording devices.\n# !\n# ! The general neuron model\n# ! ------------------------\n# !\n# ! We use the ``iaf_cond_alpha`` neuron, which is an\n# ! integrate-and-fire neuron with two conductance-based synapses which\n# ! have alpha-function time course. Any input with positive weights\n# ! will automatically directed to the synapse labeled ``_ex``, any\n# ! with negative weights to the synapes labeled ``_in``. We define\n# ! **all** parameters explicitly here, so that no information is\n# ! hidden in the model definition in NEST. ``V_m`` is the membrane\n# ! potential to which the model neurons will be initialized.\n# ! The model equations and parameters for the Hill-Tononi neuron model\n# ! are given on pp. 1677f and Tables 2 and 3 in that paper. Note some\n# ! peculiarities and adjustments:\n# !\n# ! - Hill & Tononi specify their model in terms of the membrane time\n# ! constant, while the ``iaf_cond_alpha`` model is based on the\n# ! membrane capcitance. Interestingly, conducantces are unitless in\n# ! the H&T model. We thus can use the time constant directly as\n# ! membrane capacitance.\n# ! - The model includes sodium and potassium leak conductances. We\n# ! combine these into a single one as follows:\n# $ \\begin{equation}-g_{NaL}(V-E_{Na}) - g_{KL}(V-E_K)\n# $ = -(g_{NaL}+g_{KL})\n# $ \\left(V-\\frac{g_{NaL}E_{NaL}+g_{KL}E_K}{g_{NaL}g_{KL}}\\right)\n# $ \\end{equation}\n# ! - We write the resulting expressions for g_L and E_L explicitly\n# ! below, to avoid errors in copying from our pocket calculator.\n# ! - The paper gives a range of 1.0-1.85 for g_{KL}, we choose 1.5\n# ! here.\n# ! - The Hill-Tononi model has no explicit reset or refractory\n# ! time. We arbitrarily set V_reset and t_ref.\n# ! - The paper uses double exponential time courses for the synaptic\n# ! conductances, with separate time constants for the rising and\n# ! fallings flanks. Alpha functions have only a single time\n# ! constant: we use twice the rising time constant given by Hill and\n# ! Tononi.\n# ! - In the general model below, we use the values for the cortical\n# ! excitatory cells as defaults. Values will then be adapted below.\n# !\nnest.CopyModel('iaf_cond_alpha', 'NeuronModel',\n params={'C_m': 16.0,\n 'E_L': (0.2 * 30.0 + 1.5 * -90.0) / (0.2 + 1.5),\n 'g_L': 0.2 + 1.5,\n 'E_ex': 0.0,\n 'E_in': -70.0,\n 'V_reset': -60.0,\n 'V_th': -51.0,\n 't_ref': 2.0,\n 'tau_syn_ex': 1.0,\n 'tau_syn_in': 2.0,\n 'I_e': 0.0,\n 'V_m': -70.0})\n\n# ! Adaptation of models for different populations\n# ! ----------------------------------------------\n\n# ! We must copy the `NeuronModel` dictionary explicitly, otherwise\n# ! Python would just create a reference.\n\n# ! Cortical excitatory cells\n# ! .........................\n# ! Parameters are the same as above, so we need not adapt anything\nnest.CopyModel('NeuronModel', 'CtxExNeuron')\n\n# ! Cortical inhibitory cells\n# ! .........................\nnest.CopyModel('NeuronModel', 'CtxInNeuron',\n params={'C_m': 8.0,\n 'V_th': -53.0,\n 't_ref': 1.0})\n\n# ! Thalamic cells\n# ! ..............\nnest.CopyModel('NeuronModel', 'ThalamicNeuron',\n params={'C_m': 8.0,\n 'V_th': -53.0,\n 't_ref': 1.0,\n 'E_in': -80.0})\n\n\n# ! Input generating nodes\n# ! ----------------------\n\n# ! Input is generated by sinusoidally modulate Poisson generators,\n# ! organized in a square layer of retina nodes. These nodes require a\n# ! slightly more complicated initialization than all other elements of\n# ! the network:\n# !\n# ! - Average firing rate ``rate``, firing rate modulation depth ``amplitude``,\n# ! and temporal modulation frequency ``frequency`` are the same for all\n# ! retinal nodes and are set directly below.\n# ! - The temporal phase ``phase`` of each node depends on its position in\n# ! the grating and can only be assigned after the retinal layer has\n# ! been created. We therefore specify a function for initalizing the\n# ! ``phase``. This function will be called for each node.\ndef phaseInit(pos, lam, alpha):\n '''Initializer function for phase of drifting grating nodes.\n\n pos : position (x,y) of node, in degree\n lam : wavelength of grating, in degree\n alpha: angle of grating in radian, zero is horizontal\n\n Returns number to be used as phase of sinusoidal Poisson generator.\n '''\n return 360.0 / lam * (math.cos(alpha) * pos[0] + math.sin(alpha) * pos[1])\n\n\nnest.CopyModel('sinusoidal_poisson_generator', 'RetinaNode',\n params={'amplitude': Params['retAC'],\n 'rate': Params['retDC'],\n 'frequency': Params['f_dg'],\n 'phase': 0.0,\n 'individual_spike_trains': False})\n\n# ! Recording nodes\n# ! ---------------\n\n# ! We use the new ``multimeter`` device for recording from the model\n# ! neurons. At present, ``iaf_cond_alpha`` is one of few models\n# ! supporting ``multimeter`` recording. Support for more models will\n# ! be added soon; until then, you need to use ``voltmeter`` to record\n# ! from other models.\n# !\n# ! We configure multimeter to record membrane potential to membrane\n# ! potential at certain intervals to memory only. We record the GID of\n# ! the recorded neurons, but not the time.\nnest.CopyModel('multimeter', 'RecordingNode',\n params={'interval': Params['sim_interval'],\n 'record_from': ['V_m'],\n 'record_to': ['memory'],\n 'withgid': True,\n 'withtime': False})\n\n# ! Populations\n# ! ===========\n\n# ! We now create the neuron populations in the model, again in the\n# ! form of Python dictionaries. We define them in order from eye via\n# ! thalamus to cortex.\n# !\n# ! We first define a dictionary defining common properties for all\n# ! populations\nlayerProps = {'rows': Params['N'],\n 'columns': Params['N'],\n 'extent': [Params['visSize'], Params['visSize']],\n 'edge_wrap': True}\n# ! This dictionary does not yet specify the elements to put into the\n# ! layer, since they will differ from layer to layer. We will add them\n# ! below by updating the ``'elements'`` dictionary entry for each\n# ! population.\n\n# ! Retina\n# ! ------\nlayerProps.update({'elements': 'RetinaNode'})\nretina = topo.CreateLayer(layerProps)\n\n# retina_leaves is a work-around until NEST 3.0 is released\nretina_leaves = nest.GetLeaves(retina)[0]\n\n# ! Now set phases of retinal oscillators; we use a list comprehension instead\n# ! of a loop.\n[nest.SetStatus([n], {\"phase\": phaseInit(topo.GetPosition([n])[0],\n Params[\"lambda_dg\"],\n Params[\"phi_dg\"])})\n for n in retina_leaves]\n\n# ! Thalamus\n# ! --------\n\n# ! We first introduce specific neuron models for the thalamic relay\n# ! cells and interneurons. These have identical properties, but by\n# ! treating them as different models, we can address them specifically\n# ! when building connections.\n# !\n# ! We use a list comprehension to do the model copies.\n[nest.CopyModel('ThalamicNeuron', SpecificModel) for SpecificModel in\n ('TpRelay', 'TpInter')]\n\n# ! Now we can create the layer, with one relay cell and one\n# ! interneuron per location:\nlayerProps.update({'elements': ['TpRelay', 'TpInter']})\nTp = topo.CreateLayer(layerProps)\n\n# ! Reticular nucleus\n# ! -----------------\n# ! We follow the same approach as above, even though we have only a\n# ! single neuron in each location.\n[nest.CopyModel('ThalamicNeuron', SpecificModel) for SpecificModel in\n ('RpNeuron',)]\nlayerProps.update({'elements': 'RpNeuron'})\nRp = topo.CreateLayer(layerProps)\n\n# ! Primary visual cortex\n# ! ---------------------\n\n# ! We follow again the same approach. We differentiate neuron types\n# ! between layers and between pyramidal cells and interneurons. At\n# ! each location, there are two pyramidal cells and one interneuron in\n# ! each of layers 2-3, 4, and 5-6. Finally, we need to differentiate\n# ! between vertically and horizontally tuned populations. When creating\n# ! the populations, we create the vertically and the horizontally\n# ! tuned populations as separate populations.\n\n# ! We use list comprehesions to create all neuron types:\n[nest.CopyModel('CtxExNeuron', layer + 'pyr')\n for layer in ('L23', 'L4', 'L56')]\n[nest.CopyModel('CtxInNeuron', layer + 'in')\n for layer in ('L23', 'L4', 'L56')]\n\n# ! Now we can create the populations, suffixes h and v indicate tuning\nlayerProps.update({'elements': ['L23pyr', 2, 'L23in', 1,\n 'L4pyr', 2, 'L4in', 1,\n 'L56pyr', 2, 'L56in', 1]})\nVp_h = topo.CreateLayer(layerProps)\nVp_v = topo.CreateLayer(layerProps)\n\n# ! Collect all populations\n# ! -----------------------\n\n# ! For reference purposes, e.g., printing, we collect all populations\n# ! in a tuple:\npopulations = (retina, Tp, Rp, Vp_h, Vp_v)\n\n# ! Inspection\n# ! ----------\n\n# ! We can now look at the network using `PrintNetwork`:\nnest.PrintNetwork()\n\n# ! We can also try to plot a single layer in a network. For\n# ! simplicity, we use Rp, which has only a single neuron per position.\ntopo.PlotLayer(Rp)\npylab.title('Layer Rp')\npylab.show()\n\n# ! Synapse models\n# ! ==============\n\n# ! Actual synapse dynamics, e.g., properties such as the synaptic time\n# ! course, time constants, reversal potentials, are properties of\n# ! neuron models in NEST and we set them in section `Neuron models`_\n# ! above. When we refer to *synapse models* in NEST, we actually mean\n# ! connectors which store information about connection weights and\n# ! delays, as well as port numbers at the target neuron (``rport``)\n# ! and implement synaptic plasticity. The latter two aspects are not\n# ! relevant here.\n# !\n# ! We just use NEST's ``static_synapse`` connector but copy it to\n# ! synapse models ``AMPA`` and ``GABA_A`` for the sake of\n# ! explicitness. Weights and delays are set as needed in section\n# ! `Connections`_ below, as they are different from projection to\n# ! projection. De facto, the sign of the synaptic weight decides\n# ! whether input via a connection is handle by the ``_ex`` or the\n# ! ``_in`` synapse.\nnest.CopyModel('static_synapse', 'AMPA')\nnest.CopyModel('static_synapse', 'GABA_A')\n\n# ! Connections\n# ! ====================\n\n# ! Building connections is the most complex part of network\n# ! construction. Connections are specified in Table 1 in the\n# ! Hill-Tononi paper. As pointed out above, we only consider AMPA and\n# ! GABA_A synapses here. Adding other synapses is tedious work, but\n# ! should pose no new principal challenges. We also use a uniform in\n# ! stead of a Gaussian distribution for the weights.\n# !\n# ! The model has two identical primary visual cortex populations,\n# ! ``Vp_v`` and ``Vp_h``, tuned to vertical and horizonal gratings,\n# ! respectively. The *only* difference in the connection patterns\n# ! between the two populations is the thalamocortical input to layers\n# ! L4 and L5-6 is from a population of 8x2 and 2x8 grid locations,\n# ! respectively. Furthermore, inhibitory connection in cortex go to\n# ! the opposing orientation population as to the own.\n# !\n# ! To save us a lot of code doubling, we thus defined properties\n# ! dictionaries for all connections first and then use this to connect\n# ! both populations. We follow the subdivision of connections as in\n# ! the Hill & Tononi paper.\n# !\n# ! **Note:** Hill & Tononi state that their model spans 8 degrees of\n# ! visual angle and stimuli are specified according to this. On the\n# ! other hand, all connection patterns are defined in terms of cell\n# ! grid positions. Since the NEST Topology Module defines connection\n# ! patterns in terms of the extent given in degrees, we need to apply\n# ! the following scaling factor to all lengths in connections:\ndpc = Params['visSize'] / (Params['N'] - 1)\n\n# ! We will collect all same-orientation cortico-cortical connections in\nccConnections = []\n# ! the cross-orientation cortico-cortical connections in\nccxConnections = []\n# ! and all cortico-thalamic connections in\nctConnections = []\n\n# ! Horizontal intralaminar\n# ! -----------------------\n# ! *Note:* \"Horizontal\" means \"within the same cortical layer\" in this\n# ! case.\n# !\n# ! We first define a dictionary with the (most) common properties for\n# ! horizontal intralaminar connection. We then create copies in which\n# ! we adapt those values that need adapting, and\nhorIntraBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.05, \"sigma\": 7.5 * dpc}},\n \"weights\": 1.0,\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\n# ! We use a loop to do the for for us. The loop runs over a list of\n# ! dictionaries with all values that need updating\nfor conn in [{\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L4pyr\"},\n \"mask\": {\"circular\": {\"radius\": 7.0 * dpc}}},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L4in\"},\n \"mask\": {\"circular\": {\"radius\": 7.0 * dpc}}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L56pyr\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L56in\"}}]:\n ndict = horIntraBase.copy()\n ndict.update(conn)\n ccConnections.append(ndict)\n\n# ! Vertical intralaminar\n# ! -----------------------\n# ! *Note:* \"Vertical\" means \"between cortical layers\" in this\n# ! case.\n# !\n# ! We proceed as above.\nverIntraBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 1.0, \"sigma\": 7.5 * dpc}},\n \"weights\": 2.0,\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\nfor conn in [{\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L56pyr\"},\n \"weights\": 1.0},\n {\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L23in\"},\n \"weights\": 1.0},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L4pyr\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L4in\"}}]:\n ndict = verIntraBase.copy()\n ndict.update(conn)\n ccConnections.append(ndict)\n\n# ! Intracortical inhibitory\n# ! ------------------------\n# !\n# ! We proceed as above, with the following difference: each connection\n# ! is added to the same-orientation and the cross-orientation list of\n# ! connections.\n# !\n# ! **Note:** Weights increased from -1.0 to -2.0, to make up for missing GabaB\n# !\n# ! Note that we have to specify the **weight with negative sign** to make\n# ! the connections inhibitory.\nintraInhBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"GABA_A\",\n \"mask\": {\"circular\": {\"radius\": 7.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.25, \"sigma\": 7.5 * dpc}},\n \"weights\": -2.0,\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\n# ! We use a loop to do the for for us. The loop runs over a list of\n# ! dictionaries with all values that need updating\nfor conn in [{\"sources\": {\"model\": \"L23in\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L23in\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L4in\"}, \"targets\": {\"model\": \"L4pyr\"}},\n {\"sources\": {\"model\": \"L4in\"}, \"targets\": {\"model\": \"L4in\"}},\n {\"sources\": {\"model\": \"L56in\"}, \"targets\": {\"model\": \"L56pyr\"}},\n {\"sources\": {\"model\": \"L56in\"}, \"targets\": {\"model\": \"L56in\"}}]:\n ndict = intraInhBase.copy()\n ndict.update(conn)\n ccConnections.append(ndict)\n ccxConnections.append(ndict)\n\n# ! Corticothalamic\n# ! ---------------\ncorThalBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 5.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.5, \"sigma\": 7.5 * dpc}},\n \"weights\": 1.0,\n \"delays\": {\"uniform\": {\"min\": 7.5, \"max\": 8.5}}}\n\n# ! We use a loop to do the for for us. The loop runs over a list of\n# ! dictionaries with all values that need updating\nfor conn in [{\"sources\": {\"model\": \"L56pyr\"},\n \"targets\": {\"model\": \"TpRelay\"}},\n {\"sources\": {\"model\": \"L56pyr\"},\n \"targets\": {\"model\": \"TpInter\"}}]:\n ndict = intraInhBase.copy()\n ndict.update(conn)\n ctConnections.append(ndict)\n\n# ! Corticoreticular\n# ! ----------------\n\n# ! In this case, there is only a single connection, so we write the\n# ! dictionary itself; it is very similar to the corThalBase, and to\n# ! show that, we copy first, then update. We need no ``targets`` entry,\n# ! since Rp has only one neuron per location.\ncorRet = corThalBase.copy()\ncorRet.update({\"sources\": {\"model\": \"L56pyr\"}, \"weights\": 2.5})\n\n# ! Build all connections beginning in cortex\n# ! -----------------------------------------\n\n# ! Cortico-cortical, same orientation\nprint(\"Connecting: cortico-cortical, same orientation\")\n[topo.ConnectLayers(Vp_h, Vp_h, conn) for conn in ccConnections]\n[topo.ConnectLayers(Vp_v, Vp_v, conn) for conn in ccConnections]\n\n# ! Cortico-cortical, cross-orientation\nprint(\"Connecting: cortico-cortical, other orientation\")\n[topo.ConnectLayers(Vp_h, Vp_v, conn) for conn in ccxConnections]\n[topo.ConnectLayers(Vp_v, Vp_h, conn) for conn in ccxConnections]\n\n# ! Cortico-thalamic connections\nprint(\"Connecting: cortico-thalamic\")\n[topo.ConnectLayers(Vp_h, Tp, conn) for conn in ctConnections]\n[topo.ConnectLayers(Vp_v, Tp, conn) for conn in ctConnections]\ntopo.ConnectLayers(Vp_h, Rp, corRet)\ntopo.ConnectLayers(Vp_v, Rp, corRet)\n\n# ! Thalamo-cortical connections\n# ! ----------------------------\n\n# ! **Note:** According to the text on p. 1674, bottom right, of\n# ! the Hill & Tononi paper, thalamocortical connections are\n# ! created by selecting from the thalamic population for each\n# ! L4 pyramidal cell, ie, are *convergent* connections.\n# !\n# ! We first handle the rectangular thalamocortical connections.\nthalCorRect = {\"connection_type\": \"convergent\",\n \"sources\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"AMPA\",\n \"weights\": 5.0,\n \"delays\": {\"uniform\": {\"min\": 2.75, \"max\": 3.25}}}\n\nprint(\"Connecting: thalamo-cortical\")\n\n# ! Horizontally tuned\nthalCorRect.update(\n {\"mask\": {\"rectangular\": {\"lower_left\": [-4.0 * dpc, -1.0 * dpc],\n \"upper_right\": [4.0 * dpc, 1.0 * dpc]}}})\nfor conn in [{\"targets\": {\"model\": \"L4pyr\"}, \"kernel\": 0.5},\n {\"targets\": {\"model\": \"L56pyr\"}, \"kernel\": 0.3}]:\n thalCorRect.update(conn)\n topo.ConnectLayers(Tp, Vp_h, thalCorRect)\n\n# ! Vertically tuned\nthalCorRect.update(\n {\"mask\": {\"rectangular\": {\"lower_left\": [-1.0 * dpc, -4.0 * dpc],\n \"upper_right\": [1.0 * dpc, 4.0 * dpc]}}})\nfor conn in [{\"targets\": {\"model\": \"L4pyr\"}, \"kernel\": 0.5},\n {\"targets\": {\"model\": \"L56pyr\"}, \"kernel\": 0.3}]:\n thalCorRect.update(conn)\n topo.ConnectLayers(Tp, Vp_v, thalCorRect)\n\n# ! Diffuse connections\nthalCorDiff = {\"connection_type\": \"convergent\",\n \"sources\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"AMPA\",\n \"weights\": 5.0,\n \"mask\": {\"circular\": {\"radius\": 5.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.1, \"sigma\": 7.5 * dpc}},\n \"delays\": {\"uniform\": {\"min\": 2.75, \"max\": 3.25}}}\n\nfor conn in [{\"targets\": {\"model\": \"L4pyr\"}},\n {\"targets\": {\"model\": \"L56pyr\"}}]:\n thalCorDiff.update(conn)\n topo.ConnectLayers(Tp, Vp_h, thalCorDiff)\n topo.ConnectLayers(Tp, Vp_v, thalCorDiff)\n\n# ! Thalamic connections\n# ! --------------------\n\n# ! Connections inside thalamus, including Rp\n# !\n# ! *Note:* In Hill & Tononi, the inhibition between Rp cells is mediated by\n# ! GABA_B receptors. We use GABA_A receptors here to provide some\n# ! self-dampening of Rp.\n# !\n# ! **Note:** The following code had a serious bug in v. 0.1: During the first\n# ! iteration of the loop, \"synapse_model\" and \"weights\" were set to \"AMPA\" and\n# ! \"0.1\", respectively and remained unchanged, so that all connections were\n# ! created as excitatory connections, even though they should have been\n# ! inhibitory. We now specify synapse_model and weight explicitly for each\n# ! connection to avoid this.\n\nthalBase = {\"connection_type\": \"divergent\",\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\nprint(\"Connecting: intra-thalamic\")\n\nfor src, tgt, conn in [(Tp, Rp, {\"sources\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 1.0,\n \"sigma\": 7.5 * dpc}},\n \"weights\": 2.0}),\n (Tp, Tp, {\"sources\": {\"model\": \"TpInter\"},\n \"targets\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.25,\n \"sigma\": 7.5 * dpc}}}),\n (Tp, Tp, {\"sources\": {\"model\": \"TpInter\"},\n \"targets\": {\"model\": \"TpInter\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.25,\n \"sigma\": 7.5 * dpc}}}),\n (Rp, Tp, {\"targets\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.15,\n \"sigma\": 7.5 * dpc}}}),\n (Rp, Tp, {\"targets\": {\"model\": \"TpInter\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.15,\n \"sigma\": 7.5 * dpc}}}),\n (Rp, Rp, {\"targets\": {\"model\": \"RpNeuron\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.5,\n \"sigma\": 7.5 * dpc}}})]:\n thalBase.update(conn)\n topo.ConnectLayers(src, tgt, thalBase)\n\n# ! Thalamic input\n# ! --------------\n\n# ! Input to the thalamus from the retina.\n# !\n# ! **Note:** Hill & Tononi specify a delay of 0 ms for this connection.\n# ! We use 1 ms here.\nretThal = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 1.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.75, \"sigma\": 2.5 * dpc}},\n \"weights\": 10.0,\n \"delays\": 1.0}\n\nprint(\"Connecting: retino-thalamic\")\n\nfor conn in [{\"targets\": {\"model\": \"TpRelay\"}},\n {\"targets\": {\"model\": \"TpInter\"}}]:\n retThal.update(conn)\n topo.ConnectLayers(retina, Tp, retThal)\n\n# ! Checks on connections\n# ! ---------------------\n\n# ! As a very simple check on the connections created, we inspect\n# ! the connections from the central node of various layers.\n\n# ! Connections from Retina to TpRelay\ntopo.PlotTargets(topo.FindCenterElement(retina), Tp, 'TpRelay', 'AMPA')\npylab.title('Connections Retina -> TpRelay')\npylab.show()\n\n# ! Connections from TpRelay to L4pyr in Vp (horizontally tuned)\ntopo.PlotTargets(topo.FindCenterElement(Tp), Vp_h, 'L4pyr', 'AMPA')\npylab.title('Connections TpRelay -> Vp(h) L4pyr')\npylab.show()\n\n# ! Connections from TpRelay to L4pyr in Vp (vertically tuned)\ntopo.PlotTargets(topo.FindCenterElement(Tp), Vp_v, 'L4pyr', 'AMPA')\npylab.title('Connections TpRelay -> Vp(v) L4pyr')\npylab.show()\n\n# ! Recording devices\n# ! =================\n\n# ! This recording device setup is a bit makeshift. For each population\n# ! we want to record from, we create one ``multimeter``, then select\n# ! all nodes of the right model from the target population and\n# ! connect. ``loc`` is the subplot location for the layer.\nprint(\"Connecting: Recording devices\")\nrecorders = {}\nfor name, loc, population, model in [('TpRelay', 1, Tp, 'TpRelay'),\n ('Rp', 2, Rp, 'RpNeuron'),\n ('Vp_v L4pyr', 3, Vp_v, 'L4pyr'),\n ('Vp_h L4pyr', 4, Vp_h, 'L4pyr')]:\n recorders[name] = (nest.Create('RecordingNode'), loc)\n # population_leaves is a work-around until NEST 3.0 is released\n population_leaves = nest.GetLeaves(population)[0]\n tgts = [nd for nd in population_leaves\n if nest.GetStatus([nd], 'model')[0] == model]\n nest.Connect(recorders[name][0], tgts) # one recorder to all targets\n\n# ! Example simulation\n# ! ====================\n\n# ! This simulation is set up to create a step-wise visualization of\n# ! the membrane potential. To do so, we simulate ``sim_interval``\n# ! milliseconds at a time, then read out data from the multimeters,\n# ! clear data from the multimeters and plot the data as pseudocolor\n# ! plots.\n\n# ! show time during simulation\nnest.SetStatus([0], {'print_time': True})\n\n# ! lower and upper limits for color scale, for each of the four\n# ! populations recorded.\nvmn = [-80, -80, -80, -80]\nvmx = [-50, -50, -50, -50]\n\nnest.Simulate(Params['sim_interval'])\n\n# ! loop over simulation intervals\nfor t in pylab.arange(Params['sim_interval'], Params['simtime'],\n Params['sim_interval']):\n\n # do the simulation\n nest.Simulate(Params['sim_interval'])\n\n # clear figure and choose colormap\n pylab.clf()\n pylab.jet()\n\n # now plot data from each recorder in turn, assume four recorders\n for name, r in recorders.items():\n rec = r[0]\n sp = r[1]\n pylab.subplot(2, 2, sp)\n d = nest.GetStatus(rec)[0]['events']['V_m']\n\n if len(d) != Params['N'] ** 2:\n # cortical layer with two neurons in each location, take average\n d = 0.5 * (d[::2] + d[1::2])\n\n # clear data from multimeter\n nest.SetStatus(rec, {'n_events': 0})\n pylab.imshow(pylab.reshape(d, (Params['N'], Params['N'])),\n aspect='equal', interpolation='nearest',\n extent=(0, Params['N'] + 1, 0, Params['N'] + 1),\n vmin=vmn[sp - 1], vmax=vmx[sp - 1])\n pylab.colorbar()\n pylab.title(name + ', t = %6.1f ms' % nest.GetKernelStatus()['time'])\n\n pylab.draw() # force drawing inside loop\n pylab.show() # required by ``pyreport``\n\n# ! just for some information at the end\nprint(nest.GetKernelStatus())" + "import pylab\nSHOW_FIGURES = False\n\nif not SHOW_FIGURES:\n pylab_show = pylab.show\n\n def nop(s=None):\n pass\n\n pylab.show = nop\nelse:\n pylab.ion()\n\n# ! Introduction\n# !=============\n# ! This tutorial gives a brief introduction to the ConnPlotter\n# ! toolbox. It is by no means complete.\n\n# ! Load pynest\nimport nest\n\n# ! Load NEST Topoplogy module (NEST 2.20.1)\nimport nest.topology as topo\n\n# ! Make sure we start with a clean slate, even if we re-run the script\n# ! in the same Python session.\nnest.ResetKernel()\n\n# ! Import math, we need Pi\nimport math\n\n# ! Configurable Parameters\n# ! =======================\n# !\n# ! Here we define those parameters that we take to be\n# ! configurable. The choice of configurable parameters is obviously\n# ! arbitrary, and in practice one would have far more configurable\n# ! parameters. We restrict ourselves to:\n# !\n# ! - Network size in neurons ``N``, each layer is ``N x N``.\n# ! - Network size in subtended visual angle ``visSize``, in degree.\n# ! - Temporal frequency of drifting grating input ``f_dg``, in Hz.\n# ! - Spatial wavelength and direction of drifting grating input,\n# ! ``lambda_dg`` and ``phi_dg``, in degree/radian.\n# ! - Background firing rate of retinal nodes and modulation amplitude,\n# ! ``retDC`` and ``retAC``, in Hz.\n# ! - Simulation duration ``simtime``; actual simulation is split into\n# ! intervals of ``sim_interval`` length, so that the network state\n# ! can be visualized in those intervals. Times are in ms.\nParams = {'N': 40,\n 'visSize': 8.0,\n 'f_dg': 2.0,\n 'lambda_dg': 2.0,\n 'phi_dg': 0.0,\n 'retDC': 30.0,\n 'retAC': 30.0,\n 'simtime': 100.0,\n 'sim_interval': 5.0\n }\n\n# ! Neuron Models\n# ! =============\n# !\n# ! We declare models in two steps:\n# !\n# ! 1. We define a dictionary specifying the NEST neuron model to use\n# ! as well as the parameters for that model.\n# ! #. We create three copies of this dictionary with parameters\n# ! adjusted to the three model variants specified in Table~2 of\n# ! Hill & Tononi (2005) (cortical excitatory, cortical inhibitory,\n# ! thalamic)\n# !\n# ! In addition, we declare the models for the stimulation and\n# ! recording devices.\n# !\n# ! The general neuron model\n# ! ------------------------\n# !\n# ! We use the ``iaf_cond_alpha`` neuron, which is an\n# ! integrate-and-fire neuron with two conductance-based synapses which\n# ! have alpha-function time course. Any input with positive weights\n# ! will automatically directed to the synapse labeled ``_ex``, any\n# ! with negative weights to the synapes labeled ``_in``. We define\n# ! **all** parameters explicitly here, so that no information is\n# ! hidden in the model definition in NEST. ``V_m`` is the membrane\n# ! potential to which the model neurons will be initialized.\n# ! The model equations and parameters for the Hill-Tononi neuron model\n# ! are given on pp. 1677f and Tables 2 and 3 in that paper. Note some\n# ! peculiarities and adjustments:\n# !\n# ! - Hill & Tononi specify their model in terms of the membrane time\n# ! constant, while the ``iaf_cond_alpha`` model is based on the\n# ! membrane capcitance. Interestingly, conducantces are unitless in\n# ! the H&T model. We thus can use the time constant directly as\n# ! membrane capacitance.\n# ! - The model includes sodium and potassium leak conductances. We\n# ! combine these into a single one as follows:\n# $ \\begin{equation}-g_{NaL}(V-E_{Na}) - g_{KL}(V-E_K)\n# $ = -(g_{NaL}+g_{KL})\n# $ \\left(V-\\frac{g_{NaL}E_{NaL}+g_{KL}E_K}{g_{NaL}g_{KL}}\\right)\n# $ \\end{equation}\n# ! - We write the resulting expressions for g_L and E_L explicitly\n# ! below, to avoid errors in copying from our pocket calculator.\n# ! - The paper gives a range of 1.0-1.85 for g_{KL}, we choose 1.5\n# ! here.\n# ! - The Hill-Tononi model has no explicit reset or refractory\n# ! time. We arbitrarily set V_reset and t_ref.\n# ! - The paper uses double exponential time courses for the synaptic\n# ! conductances, with separate time constants for the rising and\n# ! fallings flanks. Alpha functions have only a single time\n# ! constant: we use twice the rising time constant given by Hill and\n# ! Tononi.\n# ! - In the general model below, we use the values for the cortical\n# ! excitatory cells as defaults. Values will then be adapted below.\n# !\nnest.CopyModel('iaf_cond_alpha', 'NeuronModel',\n params={'C_m': 16.0,\n 'E_L': (0.2 * 30.0 + 1.5 * -90.0) / (0.2 + 1.5),\n 'g_L': 0.2 + 1.5,\n 'E_ex': 0.0,\n 'E_in': -70.0,\n 'V_reset': -60.0,\n 'V_th': -51.0,\n 't_ref': 2.0,\n 'tau_syn_ex': 1.0,\n 'tau_syn_in': 2.0,\n 'I_e': 0.0,\n 'V_m': -70.0})\n\n# ! Adaptation of models for different populations\n# ! ----------------------------------------------\n\n# ! We must copy the `NeuronModel` dictionary explicitly, otherwise\n# ! Python would just create a reference.\n\n# ! Cortical excitatory cells\n# ! .........................\n# ! Parameters are the same as above, so we need not adapt anything\nnest.CopyModel('NeuronModel', 'CtxExNeuron')\n\n# ! Cortical inhibitory cells\n# ! .........................\nnest.CopyModel('NeuronModel', 'CtxInNeuron',\n params={'C_m': 8.0,\n 'V_th': -53.0,\n 't_ref': 1.0})\n\n# ! Thalamic cells\n# ! ..............\nnest.CopyModel('NeuronModel', 'ThalamicNeuron',\n params={'C_m': 8.0,\n 'V_th': -53.0,\n 't_ref': 1.0,\n 'E_in': -80.0})\n\n\n# ! Input generating nodes\n# ! ----------------------\n\n# ! Input is generated by sinusoidally modulate Poisson generators,\n# ! organized in a square layer of retina nodes. These nodes require a\n# ! slightly more complicated initialization than all other elements of\n# ! the network:\n# !\n# ! - Average firing rate ``rate``, firing rate modulation depth ``amplitude``,\n# ! and temporal modulation frequency ``frequency`` are the same for all\n# ! retinal nodes and are set directly below.\n# ! - The temporal phase ``phase`` of each node depends on its position in\n# ! the grating and can only be assigned after the retinal layer has\n# ! been created. We therefore specify a function for initalizing the\n# ! ``phase``. This function will be called for each node.\ndef phaseInit(pos, lam, alpha):\n '''Initializer function for phase of drifting grating nodes.\n\n pos : position (x,y) of node, in degree\n lam : wavelength of grating, in degree\n alpha: angle of grating in radian, zero is horizontal\n\n Returns number to be used as phase of sinusoidal Poisson generator.\n '''\n return 360.0 / lam * (math.cos(alpha) * pos[0] + math.sin(alpha) * pos[1])\n\n\nnest.CopyModel('sinusoidal_poisson_generator', 'RetinaNode',\n params={'amplitude': Params['retAC'],\n 'rate': Params['retDC'],\n 'frequency': Params['f_dg'],\n 'phase': 0.0,\n 'individual_spike_trains': False})\n\n# ! Recording nodes\n# ! ---------------\n\n# ! We use the new ``multimeter`` device for recording from the model\n# ! neurons. At present, ``iaf_cond_alpha`` is one of few models\n# ! supporting ``multimeter`` recording. Support for more models will\n# ! be added soon; until then, you need to use ``voltmeter`` to record\n# ! from other models.\n# !\n# ! We configure multimeter to record membrane potential to membrane\n# ! potential at certain intervals to memory only. We record the GID of\n# ! the recorded neurons, but not the time.\nnest.CopyModel('multimeter', 'RecordingNode',\n params={'interval': Params['sim_interval'],\n 'record_from': ['V_m'],\n 'record_to': ['memory'],\n 'withgid': True,\n 'withtime': False})\n\n# ! Populations\n# ! ===========\n\n# ! We now create the neuron populations in the model, again in the\n# ! form of Python dictionaries. We define them in order from eye via\n# ! thalamus to cortex.\n# !\n# ! We first define a dictionary defining common properties for all\n# ! populations\nlayerProps = {'rows': Params['N'],\n 'columns': Params['N'],\n 'extent': [Params['visSize'], Params['visSize']],\n 'edge_wrap': True}\n# ! This dictionary does not yet specify the elements to put into the\n# ! layer, since they will differ from layer to layer. We will add them\n# ! below by updating the ``'elements'`` dictionary entry for each\n# ! population.\n\n# ! Retina\n# ! ------\nlayerProps.update({'elements': 'RetinaNode'})\nretina = topo.CreateLayer(layerProps)\n\n# retina_leaves is a work-around until NEST 3.0 is released\nretina_leaves = nest.GetLeaves(retina)[0]\n\n# ! Now set phases of retinal oscillators; we use a list comprehension instead\n# ! of a loop.\n[nest.SetStatus([n], {\"phase\": phaseInit(topo.GetPosition([n])[0],\n Params[\"lambda_dg\"],\n Params[\"phi_dg\"])})\n for n in retina_leaves]\n\n# ! Thalamus\n# ! --------\n\n# ! We first introduce specific neuron models for the thalamic relay\n# ! cells and interneurons. These have identical properties, but by\n# ! treating them as different models, we can address them specifically\n# ! when building connections.\n# !\n# ! We use a list comprehension to do the model copies.\n[nest.CopyModel('ThalamicNeuron', SpecificModel) for SpecificModel in\n ('TpRelay', 'TpInter')]\n\n# ! Now we can create the layer, with one relay cell and one\n# ! interneuron per location:\nlayerProps.update({'elements': ['TpRelay', 'TpInter']})\nTp = topo.CreateLayer(layerProps)\n\n# ! Reticular nucleus\n# ! -----------------\n# ! We follow the same approach as above, even though we have only a\n# ! single neuron in each location.\n[nest.CopyModel('ThalamicNeuron', SpecificModel) for SpecificModel in\n ('RpNeuron',)]\nlayerProps.update({'elements': 'RpNeuron'})\nRp = topo.CreateLayer(layerProps)\n\n# ! Primary visual cortex\n# ! ---------------------\n\n# ! We follow again the same approach. We differentiate neuron types\n# ! between layers and between pyramidal cells and interneurons. At\n# ! each location, there are two pyramidal cells and one interneuron in\n# ! each of layers 2-3, 4, and 5-6. Finally, we need to differentiate\n# ! between vertically and horizontally tuned populations. When creating\n# ! the populations, we create the vertically and the horizontally\n# ! tuned populations as separate populations.\n\n# ! We use list comprehesions to create all neuron types:\n[nest.CopyModel('CtxExNeuron', layer + 'pyr')\n for layer in ('L23', 'L4', 'L56')]\n[nest.CopyModel('CtxInNeuron', layer + 'in')\n for layer in ('L23', 'L4', 'L56')]\n\n# ! Now we can create the populations, suffixes h and v indicate tuning\nlayerProps.update({'elements': ['L23pyr', 2, 'L23in', 1,\n 'L4pyr', 2, 'L4in', 1,\n 'L56pyr', 2, 'L56in', 1]})\nVp_h = topo.CreateLayer(layerProps)\nVp_v = topo.CreateLayer(layerProps)\n\n# ! Collect all populations\n# ! -----------------------\n\n# ! For reference purposes, e.g., printing, we collect all populations\n# ! in a tuple:\npopulations = (retina, Tp, Rp, Vp_h, Vp_v)\n\n# ! Inspection\n# ! ----------\n\n# ! We can now look at the network using `PrintNetwork`:\nnest.PrintNetwork()\n\n# ! We can also try to plot a single layer in a network. For\n# ! simplicity, we use Rp, which has only a single neuron per position.\ntopo.PlotLayer(Rp)\npylab.title('Layer Rp')\npylab.show()\n\n# ! Synapse models\n# ! ==============\n\n# ! Actual synapse dynamics, e.g., properties such as the synaptic time\n# ! course, time constants, reversal potentials, are properties of\n# ! neuron models in NEST and we set them in section `Neuron models`_\n# ! above. When we refer to *synapse models* in NEST, we actually mean\n# ! connectors which store information about connection weights and\n# ! delays, as well as port numbers at the target neuron (``rport``)\n# ! and implement synaptic plasticity. The latter two aspects are not\n# ! relevant here.\n# !\n# ! We just use NEST's ``static_synapse`` connector but copy it to\n# ! synapse models ``AMPA`` and ``GABA_A`` for the sake of\n# ! explicitness. Weights and delays are set as needed in section\n# ! `Connections`_ below, as they are different from projection to\n# ! projection. De facto, the sign of the synaptic weight decides\n# ! whether input via a connection is handle by the ``_ex`` or the\n# ! ``_in`` synapse.\nnest.CopyModel('static_synapse', 'AMPA')\nnest.CopyModel('static_synapse', 'GABA_A')\n\n# ! Connections\n# ! ====================\n\n# ! Building connections is the most complex part of network\n# ! construction. Connections are specified in Table 1 in the\n# ! Hill-Tononi paper. As pointed out above, we only consider AMPA and\n# ! GABA_A synapses here. Adding other synapses is tedious work, but\n# ! should pose no new principal challenges. We also use a uniform in\n# ! stead of a Gaussian distribution for the weights.\n# !\n# ! The model has two identical primary visual cortex populations,\n# ! ``Vp_v`` and ``Vp_h``, tuned to vertical and horizonal gratings,\n# ! respectively. The *only* difference in the connection patterns\n# ! between the two populations is the thalamocortical input to layers\n# ! L4 and L5-6 is from a population of 8x2 and 2x8 grid locations,\n# ! respectively. Furthermore, inhibitory connection in cortex go to\n# ! the opposing orientation population as to the own.\n# !\n# ! To save us a lot of code doubling, we thus defined properties\n# ! dictionaries for all connections first and then use this to connect\n# ! both populations. We follow the subdivision of connections as in\n# ! the Hill & Tononi paper.\n# !\n# ! **Note:** Hill & Tononi state that their model spans 8 degrees of\n# ! visual angle and stimuli are specified according to this. On the\n# ! other hand, all connection patterns are defined in terms of cell\n# ! grid positions. Since the NEST Topology Module defines connection\n# ! patterns in terms of the extent given in degrees, we need to apply\n# ! the following scaling factor to all lengths in connections:\ndpc = Params['visSize'] / (Params['N'] - 1)\n\n# ! We will collect all same-orientation cortico-cortical connections in\nccConnections = []\n# ! the cross-orientation cortico-cortical connections in\nccxConnections = []\n# ! and all cortico-thalamic connections in\nctConnections = []\n\n# ! Horizontal intralaminar\n# ! -----------------------\n# ! *Note:* \"Horizontal\" means \"within the same cortical layer\" in this\n# ! case.\n# !\n# ! We first define a dictionary with the (most) common properties for\n# ! horizontal intralaminar connection. We then create copies in which\n# ! we adapt those values that need adapting, and\nhorIntraBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.05, \"sigma\": 7.5 * dpc}},\n \"weights\": 1.0,\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\n# ! We use a loop to do the for for us. The loop runs over a list of\n# ! dictionaries with all values that need updating\nfor conn in [{\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L4pyr\"},\n \"mask\": {\"circular\": {\"radius\": 7.0 * dpc}}},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L4in\"},\n \"mask\": {\"circular\": {\"radius\": 7.0 * dpc}}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L56pyr\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L56in\"}}]:\n ndict = horIntraBase.copy()\n ndict.update(conn)\n ccConnections.append(ndict)\n\n# ! Vertical intralaminar\n# ! -----------------------\n# ! *Note:* \"Vertical\" means \"between cortical layers\" in this\n# ! case.\n# !\n# ! We proceed as above.\nverIntraBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 1.0, \"sigma\": 7.5 * dpc}},\n \"weights\": 2.0,\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\nfor conn in [{\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L56pyr\"},\n \"weights\": 1.0},\n {\"sources\": {\"model\": \"L23pyr\"}, \"targets\": {\"model\": \"L23in\"},\n \"weights\": 1.0},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L4pyr\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L4pyr\"}},\n {\"sources\": {\"model\": \"L56pyr\"}, \"targets\": {\"model\": \"L4in\"}}]:\n ndict = verIntraBase.copy()\n ndict.update(conn)\n ccConnections.append(ndict)\n\n# ! Intracortical inhibitory\n# ! ------------------------\n# !\n# ! We proceed as above, with the following difference: each connection\n# ! is added to the same-orientation and the cross-orientation list of\n# ! connections.\n# !\n# ! **Note:** Weights increased from -1.0 to -2.0, to make up for missing GabaB\n# !\n# ! Note that we have to specify the **weight with negative sign** to make\n# ! the connections inhibitory.\nintraInhBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"GABA_A\",\n \"mask\": {\"circular\": {\"radius\": 7.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.25, \"sigma\": 7.5 * dpc}},\n \"weights\": -2.0,\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\n# ! We use a loop to do the for for us. The loop runs over a list of\n# ! dictionaries with all values that need updating\nfor conn in [{\"sources\": {\"model\": \"L23in\"}, \"targets\": {\"model\": \"L23pyr\"}},\n {\"sources\": {\"model\": \"L23in\"}, \"targets\": {\"model\": \"L23in\"}},\n {\"sources\": {\"model\": \"L4in\"}, \"targets\": {\"model\": \"L4pyr\"}},\n {\"sources\": {\"model\": \"L4in\"}, \"targets\": {\"model\": \"L4in\"}},\n {\"sources\": {\"model\": \"L56in\"}, \"targets\": {\"model\": \"L56pyr\"}},\n {\"sources\": {\"model\": \"L56in\"}, \"targets\": {\"model\": \"L56in\"}}]:\n ndict = intraInhBase.copy()\n ndict.update(conn)\n ccConnections.append(ndict)\n ccxConnections.append(ndict)\n\n# ! Corticothalamic\n# ! ---------------\ncorThalBase = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 5.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.5, \"sigma\": 7.5 * dpc}},\n \"weights\": 1.0,\n \"delays\": {\"uniform\": {\"min\": 7.5, \"max\": 8.5}}}\n\n# ! We use a loop to do the for for us. The loop runs over a list of\n# ! dictionaries with all values that need updating\nfor conn in [{\"sources\": {\"model\": \"L56pyr\"},\n \"targets\": {\"model\": \"TpRelay\"}},\n {\"sources\": {\"model\": \"L56pyr\"},\n \"targets\": {\"model\": \"TpInter\"}}]:\n ndict = intraInhBase.copy()\n ndict.update(conn)\n ctConnections.append(ndict)\n\n# ! Corticoreticular\n# ! ----------------\n\n# ! In this case, there is only a single connection, so we write the\n# ! dictionary itself; it is very similar to the corThalBase, and to\n# ! show that, we copy first, then update. We need no ``targets`` entry,\n# ! since Rp has only one neuron per location.\ncorRet = corThalBase.copy()\ncorRet.update({\"sources\": {\"model\": \"L56pyr\"}, \"weights\": 2.5})\n\n# ! Build all connections beginning in cortex\n# ! -----------------------------------------\n\n# ! Cortico-cortical, same orientation\nprint(\"Connecting: cortico-cortical, same orientation\")\n[topo.ConnectLayers(Vp_h, Vp_h, conn) for conn in ccConnections]\n[topo.ConnectLayers(Vp_v, Vp_v, conn) for conn in ccConnections]\n\n# ! Cortico-cortical, cross-orientation\nprint(\"Connecting: cortico-cortical, other orientation\")\n[topo.ConnectLayers(Vp_h, Vp_v, conn) for conn in ccxConnections]\n[topo.ConnectLayers(Vp_v, Vp_h, conn) for conn in ccxConnections]\n\n# ! Cortico-thalamic connections\nprint(\"Connecting: cortico-thalamic\")\n[topo.ConnectLayers(Vp_h, Tp, conn) for conn in ctConnections]\n[topo.ConnectLayers(Vp_v, Tp, conn) for conn in ctConnections]\ntopo.ConnectLayers(Vp_h, Rp, corRet)\ntopo.ConnectLayers(Vp_v, Rp, corRet)\n\n# ! Thalamo-cortical connections\n# ! ----------------------------\n\n# ! **Note:** According to the text on p. 1674, bottom right, of\n# ! the Hill & Tononi paper, thalamocortical connections are\n# ! created by selecting from the thalamic population for each\n# ! L4 pyramidal cell, ie, are *convergent* connections.\n# !\n# ! We first handle the rectangular thalamocortical connections.\nthalCorRect = {\"connection_type\": \"convergent\",\n \"sources\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"AMPA\",\n \"weights\": 5.0,\n \"delays\": {\"uniform\": {\"min\": 2.75, \"max\": 3.25}}}\n\nprint(\"Connecting: thalamo-cortical\")\n\n# ! Horizontally tuned\nthalCorRect.update(\n {\"mask\": {\"rectangular\": {\"lower_left\": [-4.0 * dpc, -1.0 * dpc],\n \"upper_right\": [4.0 * dpc, 1.0 * dpc]}}})\nfor conn in [{\"targets\": {\"model\": \"L4pyr\"}, \"kernel\": 0.5},\n {\"targets\": {\"model\": \"L56pyr\"}, \"kernel\": 0.3}]:\n thalCorRect.update(conn)\n topo.ConnectLayers(Tp, Vp_h, thalCorRect)\n\n# ! Vertically tuned\nthalCorRect.update(\n {\"mask\": {\"rectangular\": {\"lower_left\": [-1.0 * dpc, -4.0 * dpc],\n \"upper_right\": [1.0 * dpc, 4.0 * dpc]}}})\nfor conn in [{\"targets\": {\"model\": \"L4pyr\"}, \"kernel\": 0.5},\n {\"targets\": {\"model\": \"L56pyr\"}, \"kernel\": 0.3}]:\n thalCorRect.update(conn)\n topo.ConnectLayers(Tp, Vp_v, thalCorRect)\n\n# ! Diffuse connections\nthalCorDiff = {\"connection_type\": \"convergent\",\n \"sources\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"AMPA\",\n \"weights\": 5.0,\n \"mask\": {\"circular\": {\"radius\": 5.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.1, \"sigma\": 7.5 * dpc}},\n \"delays\": {\"uniform\": {\"min\": 2.75, \"max\": 3.25}}}\n\nfor conn in [{\"targets\": {\"model\": \"L4pyr\"}},\n {\"targets\": {\"model\": \"L56pyr\"}}]:\n thalCorDiff.update(conn)\n topo.ConnectLayers(Tp, Vp_h, thalCorDiff)\n topo.ConnectLayers(Tp, Vp_v, thalCorDiff)\n\n# ! Thalamic connections\n# ! --------------------\n\n# ! Connections inside thalamus, including Rp\n# !\n# ! *Note:* In Hill & Tononi, the inhibition between Rp cells is mediated by\n# ! GABA_B receptors. We use GABA_A receptors here to provide some\n# ! self-dampening of Rp.\n# !\n# ! **Note:** The following code had a serious bug in v. 0.1: During the first\n# ! iteration of the loop, \"synapse_model\" and \"weights\" were set to \"AMPA\" and\n# ! \"0.1\", respectively and remained unchanged, so that all connections were\n# ! created as excitatory connections, even though they should have been\n# ! inhibitory. We now specify synapse_model and weight explicitly for each\n# ! connection to avoid this.\n\nthalBase = {\"connection_type\": \"divergent\",\n \"delays\": {\"uniform\": {\"min\": 1.75, \"max\": 2.25}}}\n\nprint(\"Connecting: intra-thalamic\")\n\nfor src, tgt, conn in [(Tp, Rp, {\"sources\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 1.0,\n \"sigma\": 7.5 * dpc}},\n \"weights\": 2.0}),\n (Tp, Tp, {\"sources\": {\"model\": \"TpInter\"},\n \"targets\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.25,\n \"sigma\": 7.5 * dpc}}}),\n (Tp, Tp, {\"sources\": {\"model\": \"TpInter\"},\n \"targets\": {\"model\": \"TpInter\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 2.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.25,\n \"sigma\": 7.5 * dpc}}}),\n (Rp, Tp, {\"targets\": {\"model\": \"TpRelay\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.15,\n \"sigma\": 7.5 * dpc}}}),\n (Rp, Tp, {\"targets\": {\"model\": \"TpInter\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.15,\n \"sigma\": 7.5 * dpc}}}),\n (Rp, Rp, {\"targets\": {\"model\": \"RpNeuron\"},\n \"synapse_model\": \"GABA_A\",\n \"weights\": -1.0,\n \"mask\": {\"circular\": {\"radius\": 12.0 * dpc}},\n \"kernel\": {\"gaussian\":\n {\"p_center\": 0.5,\n \"sigma\": 7.5 * dpc}}})]:\n thalBase.update(conn)\n topo.ConnectLayers(src, tgt, thalBase)\n\n# ! Thalamic input\n# ! --------------\n\n# ! Input to the thalamus from the retina.\n# !\n# ! **Note:** Hill & Tononi specify a delay of 0 ms for this connection.\n# ! We use 1 ms here.\nretThal = {\"connection_type\": \"divergent\",\n \"synapse_model\": \"AMPA\",\n \"mask\": {\"circular\": {\"radius\": 1.0 * dpc}},\n \"kernel\": {\"gaussian\": {\"p_center\": 0.75, \"sigma\": 2.5 * dpc}},\n \"weights\": 10.0,\n \"delays\": 1.0}\n\nprint(\"Connecting: retino-thalamic\")\n\nfor conn in [{\"targets\": {\"model\": \"TpRelay\"}},\n {\"targets\": {\"model\": \"TpInter\"}}]:\n retThal.update(conn)\n topo.ConnectLayers(retina, Tp, retThal)\n\n# ! Checks on connections\n# ! ---------------------\n\n# ! As a very simple check on the connections created, we inspect\n# ! the connections from the central node of various layers.\n\n# ! Connections from Retina to TpRelay\ntopo.PlotTargets(topo.FindCenterElement(retina), Tp, 'TpRelay', 'AMPA')\npylab.title('Connections Retina -> TpRelay')\npylab.show()\n\n# ! Connections from TpRelay to L4pyr in Vp (horizontally tuned)\ntopo.PlotTargets(topo.FindCenterElement(Tp), Vp_h, 'L4pyr', 'AMPA')\npylab.title('Connections TpRelay -> Vp(h) L4pyr')\npylab.show()\n\n# ! Connections from TpRelay to L4pyr in Vp (vertically tuned)\ntopo.PlotTargets(topo.FindCenterElement(Tp), Vp_v, 'L4pyr', 'AMPA')\npylab.title('Connections TpRelay -> Vp(v) L4pyr')\npylab.show()\n\n# ! Recording devices\n# ! =================\n\n# ! This recording device setup is a bit makeshift. For each population\n# ! we want to record from, we create one ``multimeter``, then select\n# ! all nodes of the right model from the target population and\n# ! connect. ``loc`` is the subplot location for the layer.\nprint(\"Connecting: Recording devices\")\nrecorders = {}\nfor name, loc, population, model in [('TpRelay', 1, Tp, 'TpRelay'),\n ('Rp', 2, Rp, 'RpNeuron'),\n ('Vp_v L4pyr', 3, Vp_v, 'L4pyr'),\n ('Vp_h L4pyr', 4, Vp_h, 'L4pyr')]:\n recorders[name] = (nest.Create('RecordingNode'), loc)\n # population_leaves is a work-around until NEST 3.0 is released\n population_leaves = nest.GetLeaves(population)[0]\n tgts = [nd for nd in population_leaves\n if nest.GetStatus([nd], 'model')[0] == model]\n nest.Connect(recorders[name][0], tgts) # one recorder to all targets\n\n# ! Example simulation\n# ! ====================\n\n# ! This simulation is set up to create a step-wise visualization of\n# ! the membrane potential. To do so, we simulate ``sim_interval``\n# ! milliseconds at a time, then read out data from the multimeters,\n# ! clear data from the multimeters and plot the data as pseudocolor\n# ! plots.\n\n# ! show time during simulation\nnest.SetStatus([0], {'print_time': True})\n\n# ! lower and upper limits for color scale, for each of the four\n# ! populations recorded.\nvmn = [-80, -80, -80, -80]\nvmx = [-50, -50, -50, -50]\n\nnest.Simulate(Params['sim_interval'])\n\n# ! loop over simulation intervals\nfor t in pylab.arange(Params['sim_interval'], Params['simtime'],\n Params['sim_interval']):\n\n # do the simulation\n nest.Simulate(Params['sim_interval'])\n\n # clear figure and choose colormap\n pylab.clf()\n pylab.jet()\n\n # now plot data from each recorder in turn, assume four recorders\n for name, r in recorders.items():\n rec = r[0]\n sp = r[1]\n pylab.subplot(2, 2, sp)\n d = nest.GetStatus(rec)[0]['events']['V_m']\n\n if len(d) != Params['N'] ** 2:\n # cortical layer with two neurons in each location, take average\n d = 0.5 * (d[::2] + d[1::2])\n\n # clear data from multimeter\n nest.SetStatus(rec, {'n_events': 0})\n pylab.imshow(pylab.reshape(d, (Params['N'], Params['N'])),\n aspect='equal', interpolation='nearest',\n extent=(0, Params['N'] + 1, 0, Params['N'] + 1),\n vmin=vmn[sp - 1], vmax=vmx[sp - 1])\n pylab.colorbar()\n pylab.title(name + ', t = %6.1f ms' % nest.GetKernelStatus()['time'])\n\n pylab.draw() # force drawing inside loop\n pylab.show() # required by ``pyreport``\n\n# ! just for some information at the end\nprint(nest.GetKernelStatus())" ] } ], @@ -51,4 +51,4 @@ }, "nbformat": 4, "nbformat_minor": 0 -} \ No newline at end of file +} diff --git a/doc/tutorials/topology/hill_tononi_Vp.py b/doc/tutorials/topology/hill_tononi_Vp.py index ab98bb0e79..708b9dadfa 100644 --- a/doc/tutorials/topology/hill_tononi_Vp.py +++ b/doc/tutorials/topology/hill_tononi_Vp.py @@ -150,6 +150,17 @@ """ import pylab +# ! Introduction +# !============= +# ! This tutorial gives a brief introduction to the ConnPlotter +# ! toolbox. It is by no means complete. +# ! Load pynest +import nest +# ! Load NEST Topoplogy module (NEST 2.20.1) +import nest.topology as topo +# ! Import math, we need Pi +import math + SHOW_FIGURES = False if not SHOW_FIGURES: @@ -162,24 +173,10 @@ def nop(s=None): else: pylab.ion() -# ! Introduction -# !============= -# ! This tutorial gives a brief introduction to the ConnPlotter -# ! toolbox. It is by no means complete. - -# ! Load pynest -import nest - -# ! Load NEST Topoplogy module (NEST 2.2) -import nest.topology as topo - # ! Make sure we start with a clean slate, even if we re-run the script # ! in the same Python session. nest.ResetKernel() -# ! Import math, we need Pi -import math - # ! Configurable Parameters # ! ======================= # ! diff --git a/doc/tutorials/topology/hill_tononi_Vp.rst b/doc/tutorials/topology/hill_tononi_Vp.rst index 24886eb293..d8692247f5 100644 --- a/doc/tutorials/topology/hill_tononi_Vp.rst +++ b/doc/tutorials/topology/hill_tononi_Vp.rst @@ -151,7 +151,7 @@ functions, see documentation. # ! Load pynest import nest - # ! Load NEST Topoplogy module (NEST 2.2) + # ! Load NEST Topoplogy module (NEST 2.20.1) import nest.topology as topo # ! Make sure we start with a clean slate, even if we re-run the script diff --git a/extras/install_libboost.sh b/extras/install_libboost.sh new file mode 100644 index 0000000000..4c606ba483 --- /dev/null +++ b/extras/install_libboost.sh @@ -0,0 +1,9 @@ +#!/bin/bash + +# Install the boost library +wget --no-verbose https://dl.bintray.com/boostorg/release/1.72.0/source/boost_1_72_0.tar.gz +tar -xzf boost_1_72_0.tar.gz +rm -fr boost_1_72_0.tar.gz +cp -fr boost_1_72_0 $HOME/.cache/boost_1_72_0.install +rm -fr boost_1_72_0 + diff --git a/extras/travis_build.sh b/extras/travis_build.sh index fb225eb738..635d4f8770 100755 --- a/extras/travis_build.sh +++ b/extras/travis_build.sh @@ -24,7 +24,7 @@ # It is invoked by the top-level Travis script '.travis.yml'. # # NOTE: This shell script is tightly coupled to Python script -# 'extras/parse_travis_log.py'. +# 'extras/parse_travis_log.py'. # Any changes to message numbers (MSGBLDnnnn) or the variable name # 'file_names' have effects on the build/test-log parsing process. @@ -32,11 +32,6 @@ # Exit shell if any subcommand or pipline returns a non-zero status. set -e -mkdir -p $HOME/.matplotlib -cat > $HOME/.matplotlib/matplotlibrc < $HOME/.matplotlib/matplotlibrc < #include +#include "sliexceptions.h" + template < typename value_type_ > class BlockVector; template < typename value_type_, typename ref_, typename ptr_ > @@ -75,7 +77,7 @@ class bv_iterator using value_type = value_type_; using pointer = ptr_; using reference = ref_; - using difference_type = long int; + using difference_type = typename BlockVector< value_type >::difference_type; bv_iterator() = default; @@ -107,17 +109,23 @@ class bv_iterator bv_iterator& operator--(); bv_iterator& operator+=( difference_type ); bv_iterator& operator-=( difference_type ); - bv_iterator operator+( difference_type ); - bv_iterator operator-( difference_type ); + bv_iterator operator+( difference_type ) const; + bv_iterator operator-( difference_type ) const; bv_iterator operator++( int ); + bv_iterator operator--( int ); reference operator*() const; pointer operator->() const; difference_type operator-( const iterator& ) const; difference_type operator-( const const_iterator& ) const; + reference operator[]( difference_type n ) const; + bool operator==( const bv_iterator& ) const; bool operator!=( const bv_iterator& ) const; bool operator<( const bv_iterator& ) const; + bool operator>( const bv_iterator& ) const; + bool operator<=( const bv_iterator& ) const; + bool operator>=( const bv_iterator& ) const; private: /** @@ -143,8 +151,15 @@ class BlockVector friend class bv_iterator; public: + using value_type = value_type_; + using difference_type = typename std::vector< value_type >::difference_type; + using const_reference = const value_type&; + using const_pointer = const value_type*; using iterator = bv_iterator< value_type_, value_type_&, value_type_* >; using const_iterator = bv_iterator< value_type_, const value_type_&, const value_type_* >; + using reverse_iterator = std::reverse_iterator< iterator >; + using const_reverse_iterator = std::reverse_iterator< const_iterator >; + using size_type = size_t; /** * @brief Creates an empty BlockVector. @@ -247,6 +262,39 @@ class BlockVector */ int get_max_block_size() const; + /** + * @brief Returns the size() of the largest possible BlockVector. + */ + size_type max_size() const; + + /** + * @brief Returns a read/write reverse iterator that points to the last + * element in the BlockVector. Iteration is done in reverse element + * order. + */ + reverse_iterator rbegin(); + + /** + * @brief Returns a read-only (constant) reverse iterator that points to + * the last element in the BlockVector. Iteration is done in reverse + * element order. + */ + reverse_iterator rbegin() const; + + /** + * @brief Returns a read/write reverse iterator that points to one + * before the first element in the BlockVector. Iteration is done in + * reverse element order. + */ + reverse_iterator rend(); + + /** + * @brief Returns a read-only (constant) reverse iterator that points to + * one before the first element in the BlockVector. Iteration is done in + * reverse element order. + */ + reverse_iterator rend() const; + private: //! Vector holding blocks containing data. std::vector< std::vector< value_type_ > > blockmap_; @@ -462,6 +510,41 @@ BlockVector< value_type_ >::get_max_block_size() const return max_block_size; } +template < typename value_type_ > +inline typename BlockVector< value_type_ >::size_type +BlockVector< value_type_ >::max_size() const +{ + throw NotImplemented( "BlockVector max_size() is not implemented." ); +} + +template < typename value_type_ > +inline typename BlockVector< value_type_ >::reverse_iterator +BlockVector< value_type_ >::rbegin() +{ + throw NotImplemented( "BlockVector rbegin() is not implemented." ); +} + +template < typename value_type_ > +inline typename BlockVector< value_type_ >::reverse_iterator +BlockVector< value_type_ >::rbegin() const +{ + throw NotImplemented( "BlockVector rbegin() is not implemented." ); +} + +template < typename value_type_ > +inline typename BlockVector< value_type_ >::reverse_iterator +BlockVector< value_type_ >::rend() +{ + throw NotImplemented( "BlockVector rend() is not implemented." ); +} + +template < typename value_type_ > +inline typename BlockVector< value_type_ >::reverse_iterator +BlockVector< value_type_ >::rend() const +{ + throw NotImplemented( "BlockVector rend() is not implemented." ); +} + ///////////////////////////////////////////////////////////// // BlockVector iterator method implementation // ///////////////////////////////////////////////////////////// @@ -530,6 +613,10 @@ inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, p template < typename value_type_, typename ref_, typename ptr_ > inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, ptr_ >::operator+=( difference_type val ) { + if ( val < 0 ) + { + return operator-=( -val ); + } for ( difference_type i = 0; i < val; ++i ) { operator++(); @@ -540,6 +627,10 @@ inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, p template < typename value_type_, typename ref_, typename ptr_ > inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, ptr_ >::operator-=( difference_type val ) { + if ( val < 0 ) + { + return operator+=( -val ); + } for ( difference_type i = 0; i < val; ++i ) { operator--(); @@ -548,14 +639,16 @@ inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, p } template < typename value_type_, typename ref_, typename ptr_ > -inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator+( difference_type val ) +inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator+( + difference_type val ) const { bv_iterator tmp = *this; return tmp += val; } template < typename value_type_, typename ref_, typename ptr_ > -inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator-( difference_type val ) +inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator-( + difference_type val ) const { bv_iterator tmp = *this; return tmp -= val; @@ -569,6 +662,14 @@ inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, pt return old; } +template < typename value_type_, typename ref_, typename ptr_ > +inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator--( int ) +{ + bv_iterator< value_type_, ref_, ptr_ > old( *this ); + --( *this ); + return old; +} + template < typename value_type_, typename ref_, typename ptr_ > inline typename bv_iterator< value_type_, ref_, ptr_ >::reference bv_iterator< value_type_, ref_, ptr_ >:: operator*() const @@ -605,6 +706,13 @@ operator-( const const_iterator& other ) const return ( block_index_ - other.block_index_ ) * max_block_size + ( this_element_index - other_element_index ); } +template < typename value_type_, typename ref_, typename ptr_ > +inline typename bv_iterator< value_type_, ref_, ptr_ >::reference bv_iterator< value_type_, ref_, ptr_ >::operator[]( + difference_type n ) const +{ + return *( *this + n ); +} + template < typename value_type_, typename ref_, typename ptr_ > inline bool bv_iterator< value_type_, ref_, ptr_ >::operator==( const bv_iterator< value_type_, ref_, ptr_ >& rhs ) const @@ -625,6 +733,24 @@ inline bool bv_iterator< value_type_, ref_, ptr_ >::operator<( const bv_iterator return ( block_index_ < rhs.block_index_ or ( block_index_ == rhs.block_index_ and block_it_ < rhs.block_it_ ) ); } +template < typename value_type_, typename ref_, typename ptr_ > +inline bool bv_iterator< value_type_, ref_, ptr_ >::operator>( const bv_iterator& rhs ) const +{ + return ( block_index_ > rhs.block_index_ or ( block_index_ == rhs.block_index_ and block_it_ > rhs.block_it_ ) ); +} + +template < typename value_type_, typename ref_, typename ptr_ > +inline bool bv_iterator< value_type_, ref_, ptr_ >::operator<=( const bv_iterator& rhs ) const +{ + return operator<( rhs ) or operator==( rhs ); +} + +template < typename value_type_, typename ref_, typename ptr_ > +inline bool bv_iterator< value_type_, ref_, ptr_ >::operator>=( const bv_iterator& rhs ) const +{ + return operator>( rhs ) or operator==( rhs ); +} + template < typename value_type_, typename ref_, typename ptr_ > inline typename bv_iterator< value_type_, ref_, ptr_ >::iterator bv_iterator< value_type_, ref_, ptr_ >::const_cast_() const @@ -632,4 +758,12 @@ bv_iterator< value_type_, ref_, ptr_ >::const_cast_() const return iterator( block_vector_, block_index_, block_it_, current_block_end_ ); } +template < typename value_type_, typename ref_, typename ptr_ > +inline bv_iterator< value_type_, ref_, ptr_ > operator+( + typename bv_iterator< value_type_, ref_, ptr_ >::difference_type n, + bv_iterator< value_type_, ref_, ptr_ >& x ) +{ + return x + n; +} + #endif /* BLOCK_VECTOR_H_ */ diff --git a/libnestutil/iterator_pair.h b/libnestutil/iterator_pair.h index 8351fd7d88..155f4b66f4 100644 --- a/libnestutil/iterator_pair.h +++ b/libnestutil/iterator_pair.h @@ -26,7 +26,7 @@ #include #include -#include "../nestkernel/source.h" +#include "source.h" namespace boost { diff --git a/libnestutil/regula_falsi.h b/libnestutil/regula_falsi.h new file mode 100644 index 0000000000..a68e24168b --- /dev/null +++ b/libnestutil/regula_falsi.h @@ -0,0 +1,113 @@ +/* + * regula_falsi.h + * + * This file is part of NEST. + * + * Copyright (C) 2004 The NEST Initiative + * + * NEST is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * NEST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NEST. If not, see . + * + */ + +#ifndef REGULA_FALSI_H +#define REGULA_FALSI_H + + +namespace nest +{ + +/** + * Localize threshold crossing by using Illinois algorithm of regula falsi method. + * + * See https://en.wikipedia.org/wiki/Regula_falsi#The_Illinois_algorithm for + * details on the algorithm. + * + * @param Node model object that call function + * @param double length of interval since previous event + * @returns time from previous event to threshold crossing + */ +template < typename CN > +double +regula_falsi( const CN& node, const double dt ) +{ + double root; + double threshold_dist_root; + + int last_threshold_sign = 0; + + double a_k = 0.0; + double b_k = dt; + + double threshold_dist_a_k = node.threshold_distance( a_k ); + double threshold_dist_b_k = node.threshold_distance( b_k ); + + if ( threshold_dist_a_k * threshold_dist_b_k > 0 ) + { + throw NumericalInstability( "regula_falsi: time step too short to reach threshold." ); + } + + const int MAX_ITER = 500; + const double TERMINATION_CRITERION = 1e-14; + + for ( int iter = 0; iter < MAX_ITER; ++iter ) + { + assert( threshold_dist_b_k != threshold_dist_a_k ); + + root = ( a_k * threshold_dist_b_k - b_k * threshold_dist_a_k ) / ( threshold_dist_b_k - threshold_dist_a_k ); + threshold_dist_root = node.threshold_distance( root ); + + if ( std::abs( threshold_dist_root ) < TERMINATION_CRITERION ) + { + return root; + } + + if ( threshold_dist_a_k * threshold_dist_root > 0.0 ) + { + // threshold_dist_a_k and threshold_dist_root have the same sign + a_k = root; + threshold_dist_a_k = threshold_dist_root; + + if ( last_threshold_sign == 1 ) + { + // If threshold_dist_a_k and threshold_dist_root had the same sign in + // the last time step, we halve the value of threshold_distance_(b_k) to + // force the root in the next time step to occur on b_k's side. This is + // done to improve the convergence rate. + threshold_dist_b_k /= 2; + } + last_threshold_sign = 1; + } + else if ( threshold_dist_b_k * threshold_dist_root > 0.0 ) + { + // threshold_dist_b_k and threshold_dist_root have the same sign + b_k = root; + threshold_dist_b_k = threshold_dist_root; + + if ( last_threshold_sign == -1 ) + { + threshold_dist_a_k /= 2; + } + last_threshold_sign = -1; + } + else + { + throw NumericalInstability( "regula_falsi: Regula falsi method did not converge" ); + } + } + throw NumericalInstability( "regula_falsi: Regula falsi method did not converge during set number of iterations" ); +} + +} // namespace nest + +#endif // REGULA_FALSI_H diff --git a/models/stdp_nn_pre-centered_connection.h b/models/stdp_nn_pre-centered_connection.h index ce81b11562..4cfcf963c7 100644 --- a/models/stdp_nn_pre-centered_connection.h +++ b/models/stdp_nn_pre-centered_connection.h @@ -322,6 +322,7 @@ STDPNNPreCenteredConnection< targetidentifierT >::STDPNNPreCenteredConnection( , mu_plus_( rhs.mu_plus_ ) , mu_minus_( rhs.mu_minus_ ) , Wmax_( rhs.Wmax_ ) + , Kplus_( rhs.Kplus_ ) , t_lastspike_( rhs.t_lastspike_ ) { } diff --git a/nestkernel/CMakeLists.txt b/nestkernel/CMakeLists.txt index 672de8d194..f7afb1428e 100644 --- a/nestkernel/CMakeLists.txt +++ b/nestkernel/CMakeLists.txt @@ -23,7 +23,6 @@ set( nestkernel_sources archiving_node.h archiving_node.cpp clopath_archiving_node.h clopath_archiving_node.cpp common_synapse_properties.h common_synapse_properties.cpp - completed_checker.h completed_checker.cpp sibling_container.h sibling_container.cpp subnet.h subnet.cpp connection.h @@ -54,6 +53,7 @@ set( nestkernel_sources multirange.h multirange.cpp node.h node.cpp nodelist.h nodelist.cpp + per_thread_bool_indicator.h per_thread_bool_indicator.cpp proxynode.h proxynode.cpp recording_device.h recording_device.cpp pseudo_recording_device.h @@ -106,6 +106,7 @@ target_include_directories( nestkernel PRIVATE ${PROJECT_BINARY_DIR}/libnestutil ${PROJECT_SOURCE_DIR}/librandom ${PROJECT_SOURCE_DIR}/sli + ${PROJECT_SOURCE_DIR}/nestkernel ) install( TARGETS nestkernel diff --git a/nestkernel/completed_checker.h b/nestkernel/completed_checker.h deleted file mode 100644 index 6e8afe8829..0000000000 --- a/nestkernel/completed_checker.h +++ /dev/null @@ -1,136 +0,0 @@ -/* - * completed_checker.h - * - * This file is part of NEST. - * - * Copyright (C) 2004 The NEST Initiative - * - * NEST is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 2 of the License, or - * (at your option) any later version. - * - * NEST is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with NEST. If not, see . - * - */ - -#ifndef COMPLETED_CHECKER_H -#define COMPLETED_CHECKER_H - -// C++ includes: -#include -#include - -// Includes from nestkernel: -#include "nest_types.h" -#include "vp_manager.h" - -// Includes from sli: -#include "dictdatum.h" - -namespace nest -{ - -/** - * A thread-safe array to coordinate progress across threads during - * gather operations. Uses an array of bools instead of an STL vector, - * since omp atomic operations can not be applied to vector - * assignments. - */ - -class CompletedChecker -{ -private: - /** - * Array holding status values for all threads. Must be of type - * bool; 'bitwise and' is used below instead of 'logical and'. - */ - bool* a_; - - /** - * Size of a_, should always be identical to number of threads. - */ - size_t size_; - -public: - CompletedChecker(); - - ~CompletedChecker(); - - /** - * Returns whether all elements are false. Waits for all threads. - */ - bool all_false() const; - - /** - * Returns whether all elements are true. Waits for all threads. - */ - bool all_true() const; - - /** - * Returns whether any elements are false. Waits for all threads. - */ - bool any_false() const; - - /** - * Returns whether any elements are true. Waits for all threads. - */ - bool any_true() const; - - /** - * Clears array and sets size to zero. - */ - void clear(); - - /** - * Updates element for thread tid by computing its 'logical and' - * with given value v. - */ - void logical_and( const thread tid, const bool v ); - - /** - * Resizes array to given size. - */ - void resize( const size_t new_size, const bool v ); - - /** - * Sets element for thread tid to given value v. - */ - void set( const thread tid, const bool v ); - - /** - * Returns const reference to element at position tid. - */ - bool operator[]( const thread tid ) const; -}; - -inline void -CompletedChecker::logical_and( const thread tid, const bool v ) -{ -// Use 'bitwise and', since 'logical and' is not supported by 'omp -// atomic update'; yields same result for bool. -#pragma omp atomic update - a_[ tid ] &= v; -} - -inline void -CompletedChecker::set( const thread tid, const bool v ) -{ -#pragma omp atomic write - a_[ tid ] = v; -} - -inline bool CompletedChecker::operator[]( const thread tid ) const -{ - return a_[ tid ]; -} - -} // namespace nest - -#endif /* COMPLETED_CHECKER_H */ diff --git a/nestkernel/connection_manager.cpp b/nestkernel/connection_manager.cpp index b4735a4b64..87fe278dde 100644 --- a/nestkernel/connection_manager.cpp +++ b/nestkernel/connection_manager.cpp @@ -96,9 +96,9 @@ nest::ConnectionManager::initialize() secondary_recv_buffer_pos_.resize( num_threads ); sort_connections_by_source_ = true; - have_connections_changed_.resize( num_threads, true ); - check_primary_connections_.resize( num_threads, false ); - check_secondary_connections_.resize( num_threads, false ); + have_connections_changed_.initialize( num_threads, true ); + check_primary_connections_.initialize( num_threads, false ); + check_secondary_connections_.initialize( num_threads, false ); #pragma omp parallel { @@ -606,17 +606,17 @@ nest::ConnectionManager::connect_( Node& s, // We do not check has_primary_connections_ and secondary_connections_exist_ // directly as this led to worse performance on the supercomputer Piz Daint. - if ( not check_primary_connections_[ tid ] and is_primary ) + if ( check_primary_connections_[ tid ].is_false() and is_primary ) { #pragma omp atomic write has_primary_connections_ = true; - check_primary_connections_.set( tid, true ); + check_primary_connections_[ tid ].set_true(); } - else if ( not check_secondary_connections_[ tid ] and not is_primary ) + else if ( check_secondary_connections_[ tid ].is_false() and not is_primary ) { #pragma omp atomic write secondary_connections_exist_ = true; - check_secondary_connections_.set( tid, true ); + check_secondary_connections_[ tid ].set_true(); } } @@ -1583,12 +1583,12 @@ nest::ConnectionManager::set_have_connections_changed( const thread tid ) // Need to check if have_connections_changed_ has already been set, because if // we have a lot of threads and they all try to set the variable at once we get // performance issues on supercomputers. - if ( not have_connections_changed_[ tid ] ) + if ( have_connections_changed_[ tid ].is_false() ) { std::string msg = "New connections created, connection descriptors previously obtained using 'GetConnections' are now invalid."; LOG( M_WARNING, "ConnectionManager", msg ); - have_connections_changed_.set( tid, true ); + have_connections_changed_[ tid ].set_true(); } } @@ -1598,8 +1598,8 @@ nest::ConnectionManager::unset_have_connections_changed( const thread tid ) // Need to check if have_connections_changed_ has already been set, because if // we have a lot of threads and they all try to set the variable at once we get // performance issues on supercomputers. - if ( have_connections_changed_[ tid ] ) + if ( have_connections_changed_[ tid ].is_true() ) { - have_connections_changed_.set( tid, false ); + have_connections_changed_[ tid ].set_false(); } } diff --git a/nestkernel/connection_manager.h b/nestkernel/connection_manager.h index 95b6f7a17d..5ae8419a02 100644 --- a/nestkernel/connection_manager.h +++ b/nestkernel/connection_manager.h @@ -31,7 +31,6 @@ #include "manager_interface.h" // Includes from nestkernel: -#include "completed_checker.h" #include "conn_builder.h" #include "connection_id.h" #include "connector_base.h" @@ -39,6 +38,7 @@ #include "nest_time.h" #include "nest_timeconverter.h" #include "nest_types.h" +#include "per_thread_bool_indicator.h" #include "source_table.h" #include "target_table.h" #include "target_table_devices.h" @@ -587,7 +587,7 @@ class ConnectionManager : public ManagerInterface //! True if new connections have been created since startup or last call to //! simulate. - CompletedChecker have_connections_changed_; + PerThreadBoolIndicator have_connections_changed_; //! Whether to sort connections by source gid. bool sort_connections_by_source_; @@ -596,13 +596,13 @@ class ConnectionManager : public ManagerInterface bool has_primary_connections_; //! Check for primary connections (spikes) on each thread. - CompletedChecker check_primary_connections_; + PerThreadBoolIndicator check_primary_connections_; //! Whether secondary connections (e.g., gap junctions) exist. bool secondary_connections_exist_; //! Check for secondary connections (e.g., gap junctions) on each thread. - CompletedChecker check_secondary_connections_; + PerThreadBoolIndicator check_secondary_connections_; //! Maximum distance between (double) spike times in STDP that is //! still considered 0. See issue #894 diff --git a/nestkernel/event_delivery_manager.cpp b/nestkernel/event_delivery_manager.cpp index 221e0f183a..77bcff3e08 100644 --- a/nestkernel/event_delivery_manager.cpp +++ b/nestkernel/event_delivery_manager.cpp @@ -83,7 +83,7 @@ EventDeliveryManager::initialize() reset_timers_counters(); spike_register_.resize( num_threads ); off_grid_spike_register_.resize( num_threads ); - gather_completed_checker_.resize( num_threads, false ); + gather_completed_checker_.initialize( num_threads, false ); // Ensures that ResetKernel resets off_grid_spiking_ off_grid_spiking_ = false; buffer_size_target_data_has_changed_ = false; @@ -107,7 +107,6 @@ EventDeliveryManager::finalize() // clear the spike buffers std::vector< std::vector< std::vector< std::vector< Target > > > >().swap( spike_register_ ); std::vector< std::vector< std::vector< std::vector< OffGridTarget > > > >().swap( off_grid_spike_register_ ); - gather_completed_checker_.clear(); send_buffer_secondary_events_.clear(); recv_buffer_secondary_events_.clear(); @@ -314,16 +313,16 @@ EventDeliveryManager::gather_spike_data_( const thread tid, std::vector< SpikeDataT >& recv_buffer ) { // Assume all threads have some work to do - gather_completed_checker_.set( tid, false ); + gather_completed_checker_[ tid ].set_false(); assert( gather_completed_checker_.all_false() ); const AssignedRanks assigned_ranks = kernel().vp_manager.get_assigned_ranks( tid ); - while ( not gather_completed_checker_.all_true() ) + while ( gather_completed_checker_.any_false() ) { // Assume this is the last gather round and change to false // otherwise - gather_completed_checker_.set( tid, true ); + gather_completed_checker_[ tid ].set_true(); #pragma omp single { @@ -341,13 +340,13 @@ EventDeliveryManager::gather_spike_data_( const thread tid, // Collocate spikes to send buffer const bool collocate_completed = collocate_spike_data_buffers_( tid, assigned_ranks, send_buffer_position, spike_register_, send_buffer ); - gather_completed_checker_.logical_and( tid, collocate_completed ); + gather_completed_checker_[ tid ].logical_and( collocate_completed ); if ( off_grid_spiking_ ) { const bool collocate_completed_off_grid = collocate_spike_data_buffers_( tid, assigned_ranks, send_buffer_position, off_grid_spike_register_, send_buffer ); - gather_completed_checker_.logical_and( tid, collocate_completed_off_grid ); + gather_completed_checker_[ tid ].logical_and( collocate_completed_off_grid ); } #pragma omp barrier @@ -380,13 +379,13 @@ EventDeliveryManager::gather_spike_data_( const thread tid, // Deliver spikes from receive buffer to ring buffers. const bool deliver_completed = deliver_events_( tid, recv_buffer ); - gather_completed_checker_.logical_and( tid, deliver_completed ); + gather_completed_checker_[ tid ].logical_and( deliver_completed ); // Exit gather loop if all local threads and remote processes are // done. #pragma omp barrier // Resize mpi buffers, if necessary and allowed. - if ( not gather_completed_checker_.all_true() and kernel().mpi_manager.adaptive_spike_buffers() ) + if ( gather_completed_checker_.any_false() and kernel().mpi_manager.adaptive_spike_buffers() ) { #pragma omp single { @@ -589,7 +588,7 @@ EventDeliveryManager::gather_target_data( const thread tid ) assert( not kernel().connection_manager.is_source_table_cleared() ); // assume all threads have some work to do - gather_completed_checker_.set( tid, false ); + gather_completed_checker_[ tid ].set_false(); assert( gather_completed_checker_.all_false() ); const AssignedRanks assigned_ranks = kernel().vp_manager.get_assigned_ranks( tid ); @@ -597,11 +596,11 @@ EventDeliveryManager::gather_target_data( const thread tid ) kernel().connection_manager.prepare_target_table( tid ); kernel().connection_manager.reset_source_table_entry_point( tid ); - while ( not gather_completed_checker_.all_true() ) + while ( gather_completed_checker_.any_false() ) { // assume this is the last gather round and change to false // otherwise - gather_completed_checker_.set( tid, true ); + gather_completed_checker_[ tid ].set_true(); #pragma omp single { @@ -617,7 +616,7 @@ EventDeliveryManager::gather_target_data( const thread tid ) assigned_ranks, kernel().mpi_manager.get_send_recv_count_target_data_per_rank() ); const bool gather_completed = collocate_target_data_buffers_( tid, assigned_ranks, send_buffer_position ); - gather_completed_checker_.logical_and( tid, gather_completed ); + gather_completed_checker_[ tid ].logical_and( gather_completed ); if ( gather_completed_checker_.all_true() ) { @@ -633,11 +632,11 @@ EventDeliveryManager::gather_target_data( const thread tid ) } // of omp single const bool distribute_completed = distribute_target_data_buffers_( tid ); - gather_completed_checker_.logical_and( tid, distribute_completed ); + gather_completed_checker_[ tid ].logical_and( distribute_completed ); #pragma omp barrier // resize mpi buffers, if necessary and allowed - if ( not gather_completed_checker_.all_true() and kernel().mpi_manager.adaptive_target_buffers() ) + if ( gather_completed_checker_.any_false() and kernel().mpi_manager.adaptive_target_buffers() ) { #pragma omp single { @@ -655,7 +654,6 @@ EventDeliveryManager::collocate_target_data_buffers_( const thread tid, const AssignedRanks& assigned_ranks, SendBufferPosition& send_buffer_position ) { - unsigned int num_target_data_written = 0; thread source_rank; TargetData next_target_data; bool valid_next_target_data; @@ -685,7 +683,7 @@ EventDeliveryManager::collocate_target_data_buffers_( const thread tid, tid, assigned_ranks.begin, assigned_ranks.end, source_rank, next_target_data ); if ( valid_next_target_data ) // add valid entry to MPI buffer { - if ( send_buffer_position.idx( source_rank ) == send_buffer_position.end( source_rank ) ) + if ( send_buffer_position.is_chunk_filled( source_rank ) ) { // entry does not fit in this part of the MPI buffer any more, // so we need to reject it @@ -697,8 +695,7 @@ EventDeliveryManager::collocate_target_data_buffers_( const thread tid, // we have just rejected an entry, so source table can not be // fully read is_source_table_read = false; - if ( num_target_data_written - == ( send_buffer_position.send_recv_count_per_rank * assigned_ranks.size ) ) // buffer is full + if ( send_buffer_position.are_all_chunks_filled() ) // buffer is full { return is_source_table_read; } diff --git a/nestkernel/event_delivery_manager.h b/nestkernel/event_delivery_manager.h index 6fa07cfbc4..84d350a7bd 100644 --- a/nestkernel/event_delivery_manager.h +++ b/nestkernel/event_delivery_manager.h @@ -33,12 +33,12 @@ #include "stopwatch.h" // Includes from nestkernel: -#include "completed_checker.h" #include "event.h" #include "mpi_manager.h" // OffGridSpike #include "nest_time.h" #include "nest_types.h" #include "node.h" +#include "per_thread_bool_indicator.h" #include "target_table.h" #include "spike_data.h" #include "vp_manager.h" @@ -439,7 +439,7 @@ class EventDeliveryManager : public ManagerInterface //!< whether size of MPI buffer for communication of spikes was changed bool buffer_size_spike_data_has_changed_; - CompletedChecker gather_completed_checker_; + PerThreadBoolIndicator gather_completed_checker_; }; inline void diff --git a/nestkernel/completed_checker.cpp b/nestkernel/per_thread_bool_indicator.cpp similarity index 50% rename from nestkernel/completed_checker.cpp rename to nestkernel/per_thread_bool_indicator.cpp index 4439746298..f581a73c1c 100644 --- a/nestkernel/completed_checker.cpp +++ b/nestkernel/per_thread_bool_indicator.cpp @@ -1,5 +1,5 @@ /* - * completed_checker.cpp + * per_thread_bool_indicator.cpp * * This file is part of NEST. * @@ -20,29 +20,41 @@ * */ -#include "completed_checker.h" +#include "per_thread_bool_indicator.h" namespace nest { -CompletedChecker::CompletedChecker() - : a_( NULL ) - , size_( 0 ) +BoolIndicatorUInt64::BoolIndicatorUInt64() + : status_( false_uint64 ) { } -CompletedChecker::~CompletedChecker() +BoolIndicatorUInt64::BoolIndicatorUInt64( const bool status ) + : status_( status ) { - clear(); +} + +BoolIndicatorUInt64& PerThreadBoolIndicator::operator[]( const thread tid ) +{ + return per_thread_status_[ tid ]; +} + +void +PerThreadBoolIndicator::initialize( const thread num_threads, const bool status ) +{ + VPManager::assert_single_threaded(); + per_thread_status_.clear(); + per_thread_status_.resize( num_threads, BoolIndicatorUInt64( status ) ); } bool -CompletedChecker::all_false() const +PerThreadBoolIndicator::all_false() const { #pragma omp barrier - for ( size_t i = 0; i < size_; ++i ) + for ( auto it = per_thread_status_.begin(); it < per_thread_status_.end(); ++it ) { - if ( a_[ i ] ) + if ( it->is_true() ) { return false; } @@ -51,12 +63,12 @@ CompletedChecker::all_false() const } bool -CompletedChecker::all_true() const +PerThreadBoolIndicator::all_true() const { #pragma omp barrier - for ( size_t i = 0; i < size_; ++i ) + for ( auto it = per_thread_status_.begin(); it < per_thread_status_.end(); ++it ) { - if ( not a_[ i ] ) + if ( it->is_false() ) { return false; } @@ -65,12 +77,12 @@ CompletedChecker::all_true() const } bool -CompletedChecker::any_false() const +PerThreadBoolIndicator::any_false() const { #pragma omp barrier - for ( size_t i = 0; i < size_; ++i ) + for ( auto it = per_thread_status_.begin(); it < per_thread_status_.end(); ++it ) { - if ( not a_[ i ] ) + if ( it->is_false() ) { return true; } @@ -79,12 +91,12 @@ CompletedChecker::any_false() const } bool -CompletedChecker::any_true() const +PerThreadBoolIndicator::any_true() const { #pragma omp barrier - for ( size_t i = 0; i < size_; ++i ) + for ( auto it = per_thread_status_.begin(); it < per_thread_status_.end(); ++it ) { - if ( a_[ i ] ) + if ( it->is_true() ) { return true; } @@ -92,29 +104,4 @@ CompletedChecker::any_true() const return false; } -void -CompletedChecker::clear() -{ - VPManager::assert_single_threaded(); - if ( a_ != NULL ) - { - delete a_; - a_ = NULL; - size_ = 0; - } -} - -void -CompletedChecker::resize( const size_t new_size, const bool v ) -{ - VPManager::assert_single_threaded(); - clear(); - a_ = new bool[ new_size ]; - for ( size_t i = 0; i < new_size; ++i ) - { - a_[ i ] = v; - } - size_ = new_size; -} - } /* namespace nest */ diff --git a/nestkernel/per_thread_bool_indicator.h b/nestkernel/per_thread_bool_indicator.h new file mode 100644 index 0000000000..54bbc24557 --- /dev/null +++ b/nestkernel/per_thread_bool_indicator.h @@ -0,0 +1,140 @@ +/* + * per_thread_bool_indicator.h + * + * This file is part of NEST. + * + * Copyright (C) 2004 The NEST Initiative + * + * NEST is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * NEST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NEST. If not, see . + * + */ + +#ifndef PER_THREAD_BOOL_INDICATOR_H +#define PER_THREAD_BOOL_INDICATOR_H + +// C++ includes: +#include +#include +#include + +// Includes from nestkernel: +#include "nest_types.h" +#include "vp_manager.h" + +// Includes from sli: +#include "dictdatum.h" + +namespace nest +{ + +/** + * A wrapper class for an integer that is only allowed to take the + * values 0 and 1. Used by PerThreadBoolIndicator to create a + * thread-safe vector indicating per-thread status. See issue #1394. + */ +class BoolIndicatorUInt64 +{ +public: + BoolIndicatorUInt64(); + BoolIndicatorUInt64( const bool status ); + + bool is_true() const; + bool is_false() const; + + void set_true(); + void set_false(); + + void logical_and( const bool status ); + +private: + static constexpr std::uint_fast64_t true_uint64 = true; + static constexpr std::uint_fast64_t false_uint64 = false; + std::uint_fast64_t status_; +}; + +inline bool +BoolIndicatorUInt64::is_true() const +{ + return ( status_ == true_uint64 ); +} + +inline bool +BoolIndicatorUInt64::is_false() const +{ + return ( status_ == false_uint64 ); +} + +inline void +BoolIndicatorUInt64::set_true() +{ + status_ = true_uint64; +} + +inline void +BoolIndicatorUInt64::set_false() +{ + status_ = false_uint64; +} + +inline void +BoolIndicatorUInt64::logical_and( const bool status ) +{ + status_ = ( static_cast< bool >( status_ ) and status ); +} + +/** + * A thread-safe vector to keep track of the status across threads, + * for example during gather operations. Uses a vector of integers + * instead of a vector of bools to guarantee thread safety. + * See issue #1394. + */ +class PerThreadBoolIndicator +{ +public: + PerThreadBoolIndicator(){}; + + BoolIndicatorUInt64& operator[]( const thread tid ); + + /** + * Resize to the given number of threads and set all elements to false. + */ + void initialize( const thread num_threads, const bool status ); + + /** + * Waits for all threads and returns whether all elements are false. + */ + bool all_false() const; + + /** + * Waits for all threads and returns whether all elements are true. + */ + bool all_true() const; + + /** + * Waits for all threads and returns whether any elements are false. + */ + bool any_false() const; + + /** + * Waits for all threads and returns whether any elements are true. + */ + bool any_true() const; + +private: + std::vector< BoolIndicatorUInt64 > per_thread_status_; +}; + +} // namespace nest + +#endif /* PER_THREAD_BOOL_INDICATOR_H */ diff --git a/nestkernel/send_buffer_position.h b/nestkernel/send_buffer_position.h index 410a7d75aa..377ca3bf77 100644 --- a/nestkernel/send_buffer_position.h +++ b/nestkernel/send_buffer_position.h @@ -41,15 +41,18 @@ namespace nest class SendBufferPosition { private: + thread begin_rank_; + thread end_rank_; + thread size_; + thread max_size_; size_t num_spike_data_written_; - std::vector< unsigned int > idx_; - std::vector< unsigned int > begin_; - std::vector< unsigned int > end_; + size_t send_recv_count_per_rank_; + std::vector< thread > idx_; + std::vector< thread > begin_; + std::vector< thread > end_; thread rank_to_index_( const thread rank ) const; - const unsigned int max_size_; - public: SendBufferPosition( const AssignedRanks& assigned_ranks, const unsigned int send_recv_count_per_rank ); @@ -83,15 +86,16 @@ class SendBufferPosition bool are_all_chunks_filled() const; void increase( const thread rank ); - - const unsigned int send_recv_count_per_rank; }; inline SendBufferPosition::SendBufferPosition( const AssignedRanks& assigned_ranks, const unsigned int send_recv_count_per_rank ) - : num_spike_data_written_( 0 ) + : begin_rank_( assigned_ranks.begin ) + , end_rank_( assigned_ranks.end ) + , size_( assigned_ranks.size ) , max_size_( assigned_ranks.max_size ) - , send_recv_count_per_rank( send_recv_count_per_rank ) + , num_spike_data_written_( 0 ) + , send_recv_count_per_rank_( send_recv_count_per_rank ) { idx_.resize( assigned_ranks.size ); begin_.resize( assigned_ranks.size ); @@ -110,6 +114,8 @@ inline SendBufferPosition::SendBufferPosition( const AssignedRanks& assigned_ran inline thread SendBufferPosition::rank_to_index_( const thread rank ) const { + assert( begin_rank_ <= rank ); + assert( rank < end_rank_ ); return rank % max_size_; } @@ -140,7 +146,7 @@ SendBufferPosition::is_chunk_filled( const thread rank ) const inline bool SendBufferPosition::are_all_chunks_filled() const { - return num_spike_data_written_ == send_recv_count_per_rank * idx_.size(); + return num_spike_data_written_ == send_recv_count_per_rank_ * idx_.size(); } inline void diff --git a/nestkernel/source_table.cpp b/nestkernel/source_table.cpp index fd9c23d1a3..16a8c84967 100644 --- a/nestkernel/source_table.cpp +++ b/nestkernel/source_table.cpp @@ -45,8 +45,8 @@ nest::SourceTable::initialize() assert( sizeof( Source ) == 8 ); const thread num_threads = kernel().vp_manager.get_num_threads(); sources_.resize( num_threads ); - is_cleared_.resize( num_threads ); - saved_entry_point_.resize( num_threads ); + is_cleared_.initialize( num_threads, false ); + saved_entry_point_.initialize( num_threads, false ); current_positions_.resize( num_threads ); saved_positions_.resize( num_threads ); @@ -55,17 +55,15 @@ nest::SourceTable::initialize() const thread tid = kernel().vp_manager.get_thread_id(); sources_[ tid ].resize( 0 ); resize_sources( tid ); - is_cleared_[ tid ] = false; - saved_entry_point_[ tid ] = false; } // of omp parallel } void nest::SourceTable::finalize() { - if ( not is_cleared() ) + for ( thread tid = 0; tid < static_cast< thread >( sources_.size() ); ++tid ) { - for ( thread tid = 0; tid < static_cast< thread >( sources_.size() ); ++tid ) + if ( is_cleared_[ tid ].is_false() ) { clear( tid ); } @@ -78,13 +76,7 @@ nest::SourceTable::finalize() bool nest::SourceTable::is_cleared() const { - bool all_cleared = true; - // we only return true, if is_cleared is true for all threads - for ( thread tid = 0; tid < kernel().vp_manager.get_num_threads(); ++tid ) - { - all_cleared &= is_cleared_[ tid ]; - } - return all_cleared; + return is_cleared_.all_true(); } std::vector< BlockVector< nest::Source > >& diff --git a/nestkernel/source_table.h b/nestkernel/source_table.h index 3148ec8f9a..8edfe2ca35 100644 --- a/nestkernel/source_table.h +++ b/nestkernel/source_table.h @@ -34,6 +34,7 @@ // Includes from nestkernel: #include "mpi_manager.h" #include "nest_types.h" +#include "per_thread_bool_indicator.h" #include "source.h" #include "source_table_position.h" @@ -70,7 +71,7 @@ class SourceTable /** * Whether the 3D structure has been deleted. */ - std::vector< bool > is_cleared_; + PerThreadBoolIndicator is_cleared_; //! Needed during readout of sources_. std::vector< SourceTablePosition > current_positions_; @@ -83,7 +84,7 @@ class SourceTable * the next communication round, while filling up (possible) * remaining parts of the MPI buffer. */ - std::vector< bool > saved_entry_point_; + PerThreadBoolIndicator saved_entry_point_; /** * Minimal number of sources that need to be deleted per synapse @@ -262,7 +263,7 @@ SourceTable::clear( const thread tid ) it->clear(); } sources_[ tid ].clear(); - is_cleared_[ tid ] = true; + is_cleared_[ tid ].set_true(); } inline void @@ -282,7 +283,7 @@ SourceTable::reject_last_target_data( const thread tid ) inline void SourceTable::save_entry_point( const thread tid ) { - if ( not saved_entry_point_[ tid ] ) + if ( saved_entry_point_[ tid ].is_false() ) { saved_positions_[ tid ].tid = current_positions_[ tid ].tid; saved_positions_[ tid ].syn_id = current_positions_[ tid ].syn_id; @@ -302,7 +303,7 @@ SourceTable::save_entry_point( const thread tid ) assert( current_positions_[ tid ].lcid == -1 ); saved_positions_[ tid ].lcid = -1; } - saved_entry_point_[ tid ] = true; + saved_entry_point_[ tid ].set_true(); } } @@ -310,7 +311,7 @@ inline void SourceTable::restore_entry_point( const thread tid ) { current_positions_[ tid ] = saved_positions_[ tid ]; - saved_entry_point_[ tid ] = false; + saved_entry_point_[ tid ].set_false(); } inline void diff --git a/precise/iaf_psc_alpha_ps.cpp b/precise/iaf_psc_alpha_ps.cpp index 656c763568..f5bdaf81b3 100644 --- a/precise/iaf_psc_alpha_ps.cpp +++ b/precise/iaf_psc_alpha_ps.cpp @@ -28,6 +28,7 @@ // Includes from libnestutil: #include "numerics.h" #include "propagator_stability.h" +#include "regula_falsi.h" // Includes from nestkernel: #include "exceptions.h" @@ -272,16 +273,11 @@ nest::iaf_psc_alpha_ps::calibrate() V_.psc_norm_ex_ = 1.0 * numerics::e / P_.tau_syn_ex_; V_.psc_norm_in_ = 1.0 * numerics::e / P_.tau_syn_in_; - V_.gamma_ex_ = 1 / P_.c_m_ / ( 1 / P_.tau_syn_ex_ - 1 / P_.tau_m_ ); - V_.gamma_sq_ex_ = 1 / P_.c_m_ / ( ( 1 / P_.tau_syn_ex_ - 1 / P_.tau_m_ ) * ( 1 / P_.tau_syn_ex_ - 1 / P_.tau_m_ ) ); - - V_.gamma_in_ = 1 / P_.c_m_ / ( 1 / P_.tau_syn_in_ - 1 / P_.tau_m_ ); - V_.gamma_sq_in_ = 1 / P_.c_m_ / ( ( 1 / P_.tau_syn_in_ - 1 / P_.tau_m_ ) * ( 1 / P_.tau_syn_in_ - 1 / P_.tau_m_ ) ); - // pre-compute matrix for full time step V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ / P_.tau_m_ ); - V_.expm1_tau_syn_ex_ = numerics::expm1( -V_.h_ms_ / P_.tau_syn_ex_ ); - V_.expm1_tau_syn_in_ = numerics::expm1( -V_.h_ms_ / P_.tau_syn_in_ ); + V_.exp_tau_syn_ex_ = std::exp( -V_.h_ms_ / P_.tau_syn_ex_ ); + V_.exp_tau_syn_in_ = std::exp( -V_.h_ms_ / P_.tau_syn_in_ ); + V_.P30_ = -P_.tau_m_ / P_.c_m_ * V_.expm1_tau_m_; // these are determined according to a numeric stability criterion V_.P31_ex_ = propagator_31( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_, V_.h_ms_ ); @@ -345,6 +341,8 @@ nest::iaf_psc_alpha_ps::update( Time const& origin, const long from, const long V_.y_input_before_ = S_.y_input_; V_.I_ex_before_ = S_.I_ex_; V_.I_in_before_ = S_.I_in_; + V_.dI_ex_before_ = S_.dI_ex_; + V_.dI_in_before_ = S_.dI_in_; V_.V_m_before_ = S_.V_m_; // get first event @@ -361,6 +359,9 @@ nest::iaf_psc_alpha_ps::update( Time const& origin, const long from, const long // update membrane potential if ( not S_.is_refractory_ ) { + // If we use S_.V_m_ * std::exp( -V_.h_ms_ / P_.tau_m_ ) instead of + // V_.expm1_tau_m_ * S_.V_m_ + S_.V_m_ here, the accuracy decrease, + // see test_iaf_ps_dc_t_accuracy.sli for details. S_.V_m_ = V_.P30_ * ( P_.I_e_ + S_.y_input_ ) + V_.P31_ex_ * S_.dI_ex_ + V_.P32_ex_ * S_.I_ex_ + V_.P31_in_ * S_.dI_in_ + V_.P32_in_ * S_.I_in_ + V_.expm1_tau_m_ * S_.V_m_ + S_.V_m_; @@ -369,11 +370,11 @@ nest::iaf_psc_alpha_ps::update( Time const& origin, const long from, const long } // update synaptic currents - S_.I_ex_ = ( V_.expm1_tau_syn_ex_ + 1. ) * V_.h_ms_ * S_.dI_ex_ + ( V_.expm1_tau_syn_ex_ + 1. ) * S_.I_ex_; - S_.dI_ex_ = ( V_.expm1_tau_syn_ex_ + 1. ) * S_.dI_ex_; + S_.I_ex_ = V_.exp_tau_syn_ex_ * V_.h_ms_ * S_.dI_ex_ + V_.exp_tau_syn_ex_ * S_.I_ex_; + S_.dI_ex_ = V_.exp_tau_syn_ex_ * S_.dI_ex_; - S_.I_in_ = ( V_.expm1_tau_syn_in_ + 1. ) * V_.h_ms_ * S_.dI_in_ + ( V_.expm1_tau_syn_in_ + 1. ) * S_.I_in_; - S_.dI_in_ = ( V_.expm1_tau_syn_in_ + 1. ) * S_.dI_in_; + S_.I_in_ = V_.exp_tau_syn_in_ * V_.h_ms_ * S_.dI_in_ + V_.exp_tau_syn_in_ * S_.I_in_; + S_.dI_in_ = V_.exp_tau_syn_in_ * S_.dI_in_; /* The following must not be moved before the y1_, dI_ex_ update, since the spike-time interpolation within emit_spike_ depends @@ -430,6 +431,8 @@ nest::iaf_psc_alpha_ps::update( Time const& origin, const long from, const long // store state V_.I_ex_before_ = S_.I_ex_; V_.I_in_before_ = S_.I_in_; + V_.dI_ex_before_ = S_.dI_ex_; + V_.dI_in_before_ = S_.dI_in_; V_.V_m_before_ = S_.V_m_; last_offset = ev_offset; @@ -500,37 +503,34 @@ nest::iaf_psc_alpha_ps::handle( DataLoggingRequest& e ) void nest::iaf_psc_alpha_ps::propagate_( const double dt ) { - // needed in any case - const double ps_e_TauSyn_ex = numerics::expm1( -dt / P_.tau_syn_ex_ ); - const double ps_e_TauSyn_in = numerics::expm1( -dt / P_.tau_syn_in_ ); - // V_m_ remains unchanged at 0.0 while neuron is refractory if ( not S_.is_refractory_ ) { - const double ps_e_Tau = numerics::expm1( -dt / P_.tau_m_ ); - const double ps_P30 = -P_.tau_m_ / P_.c_m_ * ps_e_Tau; + const double expm1_tau_m = numerics::expm1( -dt / P_.tau_m_ ); - const double ps_P31_ex = V_.gamma_sq_ex_ * ps_e_Tau - V_.gamma_sq_ex_ * ps_e_TauSyn_ex - - dt * V_.gamma_ex_ * ps_e_TauSyn_ex - dt * V_.gamma_ex_; - const double ps_P32_ex = V_.gamma_ex_ * ps_e_Tau - V_.gamma_ex_ * ps_e_TauSyn_ex; + const double ps_P30 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; - const double ps_P31_in = V_.gamma_sq_in_ * ps_e_Tau - V_.gamma_sq_in_ * ps_e_TauSyn_in - - dt * V_.gamma_in_ * ps_e_TauSyn_in - dt * V_.gamma_in_; - const double ps_P32_in = V_.gamma_in_ * ps_e_Tau - V_.gamma_in_ * ps_e_TauSyn_in; + const double ps_P31_ex = propagator_31( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_, dt ); + const double ps_P32_ex = propagator_32( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_, dt ); + const double ps_P31_in = propagator_31( P_.tau_syn_in_, P_.tau_m_, P_.c_m_, dt ); + const double ps_P32_in = propagator_32( P_.tau_syn_in_, P_.tau_m_, P_.c_m_, dt ); S_.V_m_ = ps_P30 * ( P_.I_e_ + S_.y_input_ ) + ps_P31_ex * S_.dI_ex_ + ps_P32_ex * S_.I_ex_ + ps_P31_in * S_.dI_in_ - + ps_P32_in * S_.I_in_ + ps_e_Tau * S_.V_m_ + S_.V_m_; + + ps_P32_in * S_.I_in_ + S_.V_m_ * expm1_tau_m + S_.V_m_; // lower bound of membrane potential S_.V_m_ = ( S_.V_m_ < P_.U_min_ ? P_.U_min_ : S_.V_m_ ); } + const double ps_e_TauSyn_ex = std::exp( -dt / P_.tau_syn_ex_ ); + const double ps_e_TauSyn_in = std::exp( -dt / P_.tau_syn_in_ ); + // now the synaptic components - S_.I_ex_ = ( ps_e_TauSyn_ex + 1. ) * dt * S_.dI_ex_ + ( ps_e_TauSyn_ex + 1. ) * S_.I_ex_; - S_.dI_ex_ = ( ps_e_TauSyn_ex + 1. ) * S_.dI_ex_; + S_.I_ex_ = ps_e_TauSyn_ex * dt * S_.dI_ex_ + ps_e_TauSyn_ex * S_.I_ex_; + S_.dI_ex_ = ps_e_TauSyn_ex * S_.dI_ex_; - S_.I_in_ = ( ps_e_TauSyn_in + 1. ) * dt * S_.dI_in_ + ( ps_e_TauSyn_in + 1. ) * S_.I_in_; - S_.dI_in_ = ( ps_e_TauSyn_in + 1. ) * S_.dI_in_; + S_.I_in_ = ps_e_TauSyn_in * dt * S_.dI_in_ + ps_e_TauSyn_in * S_.I_in_; + S_.dI_in_ = ps_e_TauSyn_in * S_.dI_in_; } void @@ -540,7 +540,9 @@ nest::iaf_psc_alpha_ps::emit_spike_( Time const& origin, const long lag, const d // compute spike time relative to beginning of step S_.last_spike_step_ = origin.get_steps() + lag + 1; - S_.last_spike_offset_ = V_.h_ms_ - ( t0 + bisectioning_( dt ) ); + S_.last_spike_offset_ = V_.h_ms_ - ( t0 + regula_falsi( *this, dt ) ); + + assert( S_.last_spike_offset_ >= 0.0 ); // reset neuron and make it refractory S_.V_m_ = P_.U_reset_; @@ -578,36 +580,20 @@ nest::iaf_psc_alpha_ps::emit_instant_spike_( Time const& origin, const long lag, } double -nest::iaf_psc_alpha_ps::bisectioning_( const double dt ) const +nest::iaf_psc_alpha_ps::threshold_distance( double t_step ) const { - double root = 0.0; + const double expm1_tau_m = numerics::expm1( -t_step / P_.tau_m_ ); - double V_m_root = V_.V_m_before_; + const double ps_P30 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; - double div = 2.0; + const double ps_P31_ex = propagator_31( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_, t_step ); + const double ps_P32_ex = propagator_32( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_, t_step ); + const double ps_P31_in = propagator_31( P_.tau_syn_in_, P_.tau_m_, P_.c_m_, t_step ); + const double ps_P32_in = propagator_32( P_.tau_syn_in_, P_.tau_m_, P_.c_m_, t_step ); - while ( fabs( P_.U_th_ - V_m_root ) > 1e-14 ) - { - if ( V_m_root > P_.U_th_ ) - { - root -= dt / div; - } - else - { - root += dt / div; - } - - div *= 2.0; + double V_m_root = ps_P30 * ( P_.I_e_ + V_.y_input_before_ ) + ps_P31_ex * V_.dI_ex_before_ + + ps_P32_ex * V_.I_ex_before_ + ps_P31_in * V_.dI_in_before_ + ps_P32_in * V_.I_in_before_ + + V_.V_m_before_ * expm1_tau_m + V_.V_m_before_; - const double expm1_tau_m = numerics::expm1( -root / P_.tau_m_ ); - - const double P20 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; - - const double P21_ex = propagator_32( P_.tau_syn_ex_, P_.tau_m_, P_.c_m_, root ); - const double P21_in = propagator_32( P_.tau_syn_in_, P_.tau_m_, P_.c_m_, root ); - - V_m_root = P20 * ( P_.I_e_ + V_.y_input_before_ ) + P21_ex * V_.I_ex_before_ + P21_in * V_.I_in_before_ - + expm1_tau_m * V_.V_m_before_ + V_.V_m_before_; - } - return root; + return V_m_root - P_.U_th_; } diff --git a/precise/iaf_psc_alpha_ps.h b/precise/iaf_psc_alpha_ps.h index 841d166cd8..6985992164 100644 --- a/precise/iaf_psc_alpha_ps.h +++ b/precise/iaf_psc_alpha_ps.h @@ -45,7 +45,7 @@ namespace nest /** @BeginDocumentation Name: iaf_psc_alpha_ps - Leaky integrate-and-fire neuron -with alpha-shape postsynaptic currents and bisectioning method for +with alpha-shape postsynaptic currents and regula falsi method for approximation of threshold crossing. .. versionadded:: 2.18 @@ -63,8 +63,8 @@ The precise implementation handles neuronal dynamics in a locally event-based manner with in coarse time grid defined by the minimum delay in the network, see [1]. Incoming spikes are applied at the precise moment of their arrival, while the precise time of outgoing -spikes is determined by a bisectioning method to approximate the timing -of a threshold crossing [1,3]. Return from refractoriness occurs precisly +spikes is determined by a Regula Falsi method to approximate the timing +of a threshold crossing [1,3]. Return from refractoriness occurs precisely at spike time plus refractory period. This implementation is more complex than the plain iaf_psc_alpha @@ -172,6 +172,17 @@ class iaf_psc_alpha_ps : public Archiving_Node void get_status( DictionaryDatum& ) const; void set_status( const DictionaryDatum& ); + /** + * Based on the current state, compute the value of the membrane potential + * after taking a timestep of length ``t_step``, and use it to compute the + * signed distance to spike threshold at that time. The internal state is not + * actually updated (method is defined const). + * + * @param double time step + * @returns difference between updated membrane potential and threshold + */ + double threshold_distance( double t_step ) const; + private: /** @name Interface functions * @note These functions are private, so that they can be accessed @@ -213,7 +224,7 @@ class iaf_psc_alpha_ps : public Archiving_Node void propagate_( const double dt ); /** - * Trigger interpolation method to find the precise spike time + * Trigger regula falsi method to find the precise spike time * within the mini-timestep (t0,t0+dt] assuming that the membrane * potential was below threshold at t0 and above at t0+dt. Emit * the spike and reset the neuron. @@ -235,33 +246,6 @@ class iaf_psc_alpha_ps : public Archiving_Node */ void emit_instant_spike_( Time const& origin, const long lag, const double spike_offset ); - /** @name Threshold-crossing interpolation - * These functions determine the time of threshold crossing using - * interpolation, one function per interpolation - * order. thresh_find() is the driver function and the only one to - * be called directly. - */ - //@{ - - /** Interpolation orders. */ - enum interpOrder - { - NO_INTERPOL, - LINEAR, - QUADRATIC, - CUBIC, - END_INTERP_ORDER - }; - - /** - * Localize threshold crossing by bisectioning. - * @param double length of interval since previous event - * @returns time from previous event to threshold crossing - */ - double bisectioning_( const double dt ) const; - //@} - - // The next two classes need to be friends to access the State_ class/member friend class RecordablesMap< iaf_psc_alpha_ps >; friend class UniversalDataLogger< iaf_psc_alpha_ps >; @@ -374,26 +358,24 @@ class iaf_psc_alpha_ps : public Archiving_Node */ struct Variables_ { - double h_ms_; //!< time resolution in ms - double psc_norm_ex_; //!< e / tau_syn_ex - double psc_norm_in_; //!< e / tau_syn_in - long refractory_steps_; //!< refractory time in steps - double gamma_ex_; //!< 1/c_m * 1/(1/tau_syn_ex - 1/tau_m) - double gamma_sq_ex_; //!< 1/c_m * 1/(1/tau_syn_ex - 1/tau_m)^2 - double gamma_in_; //!< 1/c_m * 1/(1/tau_syn_in - 1/tau_m) - double gamma_sq_in_; //!< 1/c_m * 1/(1/tau_syn_in - 1/tau_m)^2 - double expm1_tau_m_; //!< exp(-h/tau_m) - 1 - double expm1_tau_syn_ex_; //!< exp(-h/tau_syn_ex) - 1 - double expm1_tau_syn_in_; //!< exp(-h/tau_syn_in) - 1 - double P30_; //!< progagator matrix elem, 3rd row - double P31_ex_; //!< progagator matrix elem, 3rd row (ex) - double P32_ex_; //!< progagator matrix elem, 3rd row (ex) - double P31_in_; //!< progagator matrix elem, 3rd row (in) - double P32_in_; //!< progagator matrix elem, 3rd row (in) - double y_input_before_; //!< at beginning of mini-step, for interpolation - double I_ex_before_; //!< at beginning of mini-step, for interpolation - double I_in_before_; //!< at beginning of mini-step, for interpolation - double V_m_before_; //!< at beginning of mini-step, for interpolation + double h_ms_; //!< time resolution in ms + double psc_norm_ex_; //!< e / tau_syn_ex + double psc_norm_in_; //!< e / tau_syn_in + long refractory_steps_; //!< refractory time in steps + double expm1_tau_m_; //!< exp(-h/tau_m) - 1 + double exp_tau_syn_ex_; //!< exp(-h/tau_syn_ex) + double exp_tau_syn_in_; //!< exp(-h/tau_syn_in) + double P30_; //!< progagator matrix elem, 3rd row + double P31_ex_; //!< progagator matrix elem, 3rd row (ex) + double P32_ex_; //!< progagator matrix elem, 3rd row (ex) + double P31_in_; //!< progagator matrix elem, 3rd row (in) + double P32_in_; //!< progagator matrix elem, 3rd row (in) + double y_input_before_; //!< at beginning of mini-step + double I_ex_before_; //!< at beginning of mini-step + double I_in_before_; //!< at beginning of mini-step + double dI_ex_before_; //!< at beginning of mini-step + double dI_in_before_; //!< at beginning of mini-step + double V_m_before_; //!< at beginning of mini-step }; // Access functions for UniversalDataLogger ------------------------------- diff --git a/precise/iaf_psc_exp_ps.cpp b/precise/iaf_psc_exp_ps.cpp index 2f02f33703..a0c03af46e 100644 --- a/precise/iaf_psc_exp_ps.cpp +++ b/precise/iaf_psc_exp_ps.cpp @@ -28,6 +28,7 @@ // Includes from libnestutil: #include "numerics.h" #include "propagator_stability.h" +#include "regula_falsi.h" // Includes from nestkernel: #include "exceptions.h" @@ -256,10 +257,10 @@ nest::iaf_psc_exp_ps::calibrate() V_.h_ms_ = Time::get_resolution().get_ms(); - V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ / P_.tau_m_ ); - V_.expm1_tau_ex_ = numerics::expm1( -V_.h_ms_ / P_.tau_ex_ ); - V_.expm1_tau_in_ = numerics::expm1( -V_.h_ms_ / P_.tau_in_ ); - V_.P20_ = -P_.tau_m_ / P_.c_m_ * V_.expm1_tau_m_; + V_.exp_tau_m_ = std::exp( -V_.h_ms_ / P_.tau_m_ ); + V_.exp_tau_ex_ = std::exp( -V_.h_ms_ / P_.tau_ex_ ); + V_.exp_tau_in_ = std::exp( -V_.h_ms_ / P_.tau_in_ ); + V_.P20_ = -P_.tau_m_ / P_.c_m_ * numerics::expm1( -V_.h_ms_ / P_.tau_m_ ); // these are determined according to a numeric stability criterion V_.P21_ex_ = propagator_32( P_.tau_ex_, P_.tau_m_, P_.c_m_, V_.h_ms_ ); @@ -328,16 +329,16 @@ nest::iaf_psc_exp_ps::update( const Time& origin, const long from, const long to // update membrane potential if ( not S_.is_refractory_ ) { - S_.y2_ = V_.P20_ * ( P_.I_e_ + S_.y0_ ) + V_.P21_ex_ * S_.y1_ex_ + V_.P21_in_ * S_.y1_in_ - + V_.expm1_tau_m_ * S_.y2_ + S_.y2_; + S_.y2_ = + V_.P20_ * ( P_.I_e_ + S_.y0_ ) + V_.P21_ex_ * S_.y1_ex_ + V_.P21_in_ * S_.y1_in_ + S_.y2_ * V_.exp_tau_m_; // lower bound of membrane potential S_.y2_ = ( S_.y2_ < P_.U_min_ ? P_.U_min_ : S_.y2_ ); } // update synaptic currents - S_.y1_ex_ = S_.y1_ex_ * V_.expm1_tau_ex_ + S_.y1_ex_; - S_.y1_in_ = S_.y1_in_ * V_.expm1_tau_in_ + S_.y1_in_; + S_.y1_ex_ = S_.y1_ex_ * V_.exp_tau_ex_; + S_.y1_in_ = S_.y1_in_ * V_.exp_tau_in_; /* The following must not be moved before the y1_, y2_ update, since the spike-time interpolation within emit_spike_ depends @@ -471,22 +472,22 @@ nest::iaf_psc_exp_ps::propagate_( const double dt ) // propagate_() shall not be called then; see #368. assert( dt > 0 ); - const double expm1_tau_ex = numerics::expm1( -dt / P_.tau_ex_ ); - const double expm1_tau_in = numerics::expm1( -dt / P_.tau_in_ ); - if ( not S_.is_refractory_ ) { - const double expm1_tau_m = numerics::expm1( -dt / P_.tau_m_ ); - - const double P20 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; + const double P20 = -P_.tau_m_ / P_.c_m_ * numerics::expm1( -dt / P_.tau_m_ ); const double P21_ex = propagator_32( P_.tau_ex_, P_.tau_m_, P_.c_m_, dt ); const double P21_in = propagator_32( P_.tau_in_, P_.tau_m_, P_.c_m_, dt ); - S_.y2_ = P20 * ( P_.I_e_ + S_.y0_ ) + P21_ex * S_.y1_ex_ + P21_in * S_.y1_in_ + expm1_tau_m * S_.y2_ + S_.y2_; + S_.y2_ = + P20 * ( P_.I_e_ + S_.y0_ ) + P21_ex * S_.y1_ex_ + P21_in * S_.y1_in_ + S_.y2_ * std::exp( -dt / P_.tau_m_ ); } - S_.y1_ex_ = S_.y1_ex_ * expm1_tau_ex + S_.y1_ex_; - S_.y1_in_ = S_.y1_in_ * expm1_tau_in + S_.y1_in_; + + const double exp_tau_ex = std::exp( -dt / P_.tau_ex_ ); + const double exp_tau_in = std::exp( -dt / P_.tau_in_ ); + + S_.y1_ex_ = S_.y1_ex_ * exp_tau_ex; + S_.y1_in_ = S_.y1_in_ * exp_tau_in; } void @@ -501,7 +502,7 @@ nest::iaf_psc_exp_ps::emit_spike_( const Time& origin, const long lag, const dou // compute spike time relative to beginning of step S_.last_spike_step_ = origin.get_steps() + lag + 1; - S_.last_spike_offset_ = V_.h_ms_ - ( t0 + bisectioning_( dt ) ); + S_.last_spike_offset_ = V_.h_ms_ - ( t0 + regula_falsi( *this, dt ) ); // reset neuron and make it refractory S_.y2_ = P_.U_reset_; @@ -537,36 +538,15 @@ nest::iaf_psc_exp_ps::emit_instant_spike_( const Time& origin, const long lag, c } double -nest::iaf_psc_exp_ps::bisectioning_( const double dt ) const +nest::iaf_psc_exp_ps::threshold_distance( double t_step ) const { - double root = 0.0; - - double y2_root = V_.y2_before_; - - double div = 2.0; - - while ( fabs( P_.U_th_ - y2_root ) > 1e-14 ) - { - if ( y2_root > P_.U_th_ ) - { - root -= dt / div; - } - else - { - root += dt / div; - } - - div *= 2.0; + const double P20 = -P_.tau_m_ / P_.c_m_ * numerics::expm1( -t_step / P_.tau_m_ ); - const double expm1_tau_m = numerics::expm1( -root / P_.tau_m_ ); + const double P21_ex = propagator_32( P_.tau_ex_, P_.tau_m_, P_.c_m_, t_step ); + const double P21_in = propagator_32( P_.tau_in_, P_.tau_m_, P_.c_m_, t_step ); - const double P20 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; + double y2_root = P20 * ( P_.I_e_ + V_.y0_before_ ) + P21_ex * V_.y1_ex_before_ + P21_in * V_.y1_in_before_ + + V_.y2_before_ * std::exp( -t_step / P_.tau_m_ ); - const double P21_ex = propagator_32( P_.tau_ex_, P_.tau_m_, P_.c_m_, root ); - const double P21_in = propagator_32( P_.tau_in_, P_.tau_m_, P_.c_m_, root ); - - y2_root = P20 * ( P_.I_e_ + V_.y0_before_ ) + P21_ex * V_.y1_ex_before_ + P21_in * V_.y1_in_before_ - + expm1_tau_m * V_.y2_before_ + V_.y2_before_; - } - return root; + return y2_root - P_.U_th_; } diff --git a/precise/iaf_psc_exp_ps.h b/precise/iaf_psc_exp_ps.h index 56a17bec48..d6714ad87e 100644 --- a/precise/iaf_psc_exp_ps.h +++ b/precise/iaf_psc_exp_ps.h @@ -47,27 +47,27 @@ namespace nest /** @BeginDocumentation Name: iaf_psc_exp_ps - Leaky integrate-and-fire neuron with exponential postsynaptic currents; canoncial implementation; -bisectioning method for approximation of threshold crossing. +regula falsi method for approximation of threshold crossing. Description: iaf_psc_exp_ps is the "canonical" implementation of the leaky integrate-and-fire model neuron with exponential postsynaptic currents -that uses the bisectioning method to approximate the timing of a threshold -crossing [1,2]. This is the most exact implementation available. +that uses the regula falsi method to approximate the timing of a threshold +crossing. This is the most exact implementation available. The canonical implementation handles neuronal dynamics in a locally event-based manner with in coarse time grid defined by the minimum delay in the network, see [1,2]. Incoming spikes are applied at the precise moment of their arrival, while the precise time of outgoing -spikes is determined by bisectioning once a threshold crossing has +spikes is determined by regula falsi once a threshold crossing has been detected. Return from refractoriness occurs precisely at spike time plus refractory period. This implementation is more complex than the plain iaf_psc_exp neuron, but achieves much higher precision. In particular, it does not suffer any binning of spike times to grid points. Depending on your -application, the canonical application with bisectioning may provide +application, the canonical application with regula falsi may provide superior overall performance given an accuracy goal; see [1,2] for details. Subthreshold dynamics are integrated using exact integration between events [3]. @@ -165,6 +165,17 @@ class iaf_psc_exp_ps : public Archiving_Node void get_status( DictionaryDatum& ) const; void set_status( const DictionaryDatum& ); + /** + * Based on the current state, compute the value of the membrane potential + * after taking a timestep of length ``t_step``, and use it to compute the + * signed distance to spike threshold at that time. The internal state is not + * actually updated (method is defined const). + * + * @param double time step + * @returns difference between updated membrane potential and threshold + */ + double threshold_distance( double t_step ) const; + private: /** @name Interface functions * @note These functions are private, so that they can be accessed @@ -229,13 +240,6 @@ class iaf_psc_exp_ps : public Archiving_Node */ void emit_instant_spike_( const Time& origin, const long lag, const double spike_offset ); - /** - * Localize threshold crossing by bisectioning. - * @param double length of interval since previous event - * @returns time from previous event to threshold crossing - */ - double bisectioning_( const double dt ) const; - // ---------------------------------------------------------------- /** @@ -345,9 +349,9 @@ class iaf_psc_exp_ps : public Archiving_Node { double h_ms_; //!< Time resolution [ms] long refractory_steps_; //!< Refractory time in steps - double expm1_tau_m_; //!< exp(-h/tau_m) - 1 - double expm1_tau_ex_; //!< exp(-h/tau_ex) - 1 - double expm1_tau_in_; //!< exp(-h/tau_in) - 1 + double exp_tau_m_; //!< exp(-h/tau_m) + double exp_tau_ex_; //!< exp(-h/tau_ex) + double exp_tau_in_; //!< exp(-h/tau_in) double P20_; //!< Progagator matrix element, 2nd row double P21_in_; //!< Progagator matrix element, 2nd row double P21_ex_; //!< Progagator matrix element, 2nd row diff --git a/precise/iaf_psc_exp_ps_lossless.cpp b/precise/iaf_psc_exp_ps_lossless.cpp index 3214829be6..7c63b12be7 100644 --- a/precise/iaf_psc_exp_ps_lossless.cpp +++ b/precise/iaf_psc_exp_ps_lossless.cpp @@ -29,6 +29,10 @@ #include "exceptions.h" #include "universal_data_logger_impl.h" +// Includes from libnestutil: +#include "propagator_stability.h" +#include "regula_falsi.h" + // Includes from sli: #include "dict.h" #include "dictutils.h" @@ -280,12 +284,14 @@ nest::iaf_psc_exp_ps_lossless::calibrate() V_.h_ms_ = Time::get_resolution().get_ms(); - V_.expm1_tau_m_ = numerics::expm1( -V_.h_ms_ / P_.tau_m_ ); - V_.expm1_tau_ex_ = numerics::expm1( -V_.h_ms_ / P_.tau_ex_ ); - V_.expm1_tau_in_ = numerics::expm1( -V_.h_ms_ / P_.tau_in_ ); - V_.P20_ = -P_.tau_m_ / P_.c_m_ * V_.expm1_tau_m_; - V_.P21_ex_ = -P_.tau_m_ * P_.tau_ex_ / ( P_.tau_m_ - P_.tau_ex_ ) / P_.c_m_ * ( V_.expm1_tau_ex_ - V_.expm1_tau_m_ ); - V_.P21_in_ = -P_.tau_m_ * P_.tau_in_ / ( P_.tau_m_ - P_.tau_in_ ) / P_.c_m_ * ( V_.expm1_tau_in_ - V_.expm1_tau_m_ ); + V_.exp_tau_m_ = std::exp( -V_.h_ms_ / P_.tau_m_ ); + V_.exp_tau_ex_ = std::exp( -V_.h_ms_ / P_.tau_ex_ ); + V_.exp_tau_in_ = std::exp( -V_.h_ms_ / P_.tau_in_ ); + + V_.P20_ = -P_.tau_m_ / P_.c_m_ * numerics::expm1( -V_.h_ms_ / P_.tau_m_ ); + V_.P21_ex_ = propagator_32( P_.tau_ex_, P_.tau_m_, P_.c_m_, V_.h_ms_ ); + V_.P21_in_ = propagator_32( P_.tau_in_, P_.tau_m_, P_.c_m_, V_.h_ms_ ); + V_.refractory_steps_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); assert( V_.refractory_steps_ >= 0 ); // since t_ref_ >= 0, this can only fail in error @@ -366,15 +372,15 @@ nest::iaf_psc_exp_ps_lossless::update( const Time& origin, const long from, cons if ( not S_.is_refractory_ ) { S_.y2_ = V_.P20_ * ( P_.I_e_ + S_.y0_ ) + V_.P21_ex_ * S_.I_syn_ex_ + V_.P21_in_ * S_.I_syn_in_ - + V_.expm1_tau_m_ * S_.y2_ + S_.y2_; + + S_.y2_ * V_.exp_tau_m_; // lower bound of membrane potential S_.y2_ = ( S_.y2_ < P_.U_min_ ? P_.U_min_ : S_.y2_ ); } // update synaptic currents - S_.I_syn_ex_ = S_.I_syn_ex_ * V_.expm1_tau_ex_ + S_.I_syn_ex_; - S_.I_syn_in_ = S_.I_syn_in_ * V_.expm1_tau_in_ + S_.I_syn_in_; + S_.I_syn_ex_ = S_.I_syn_ex_ * V_.exp_tau_ex_; + S_.I_syn_in_ = S_.I_syn_in_ * V_.exp_tau_in_; /* The following must not be moved before the y1_, y2_ update, since the spike-time interpolation within emit_spike_ depends @@ -513,24 +519,22 @@ nest::iaf_psc_exp_ps_lossless::propagate_( const double dt ) // propagate_() shall not be called then; see #368. assert( dt > 0 ); - const double expm1_tau_ex = numerics::expm1( -dt / P_.tau_ex_ ); - const double expm1_tau_in = numerics::expm1( -dt / P_.tau_in_ ); - if ( not S_.is_refractory_ ) { - const double expm1_tau_m = numerics::expm1( -dt / P_.tau_m_ ); + const double P20 = -P_.tau_m_ / P_.c_m_ * numerics::expm1( -dt / P_.tau_m_ ); - const double P20 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; - const double P21_ex = - -P_.tau_m_ * P_.tau_ex_ / ( P_.tau_m_ - P_.tau_ex_ ) / P_.c_m_ * ( expm1_tau_ex - expm1_tau_m ); - const double P21_in = - -P_.tau_m_ * P_.tau_in_ / ( P_.tau_m_ - P_.tau_in_ ) / P_.c_m_ * ( expm1_tau_in - expm1_tau_m ); + const double P21_ex = propagator_32( P_.tau_ex_, P_.tau_m_, P_.c_m_, dt ); + const double P21_in = propagator_32( P_.tau_in_, P_.tau_m_, P_.c_m_, dt ); - S_.y2_ = P20 * ( P_.I_e_ + S_.y0_ ) + P21_ex * S_.I_syn_ex_ + P21_in * S_.I_syn_in_ + expm1_tau_m * S_.y2_ + S_.y2_; + S_.y2_ = + P20 * ( P_.I_e_ + S_.y0_ ) + P21_ex * S_.I_syn_ex_ + P21_in * S_.I_syn_in_ + S_.y2_ * std::exp( -dt / P_.tau_m_ ); } - S_.I_syn_ex_ = S_.I_syn_ex_ * expm1_tau_ex + S_.I_syn_ex_; - S_.I_syn_in_ = S_.I_syn_in_ * expm1_tau_in + S_.I_syn_in_; + const double exp_tau_ex = std::exp( -dt / P_.tau_ex_ ); + const double exp_tau_in = std::exp( -dt / P_.tau_in_ ); + + S_.I_syn_ex_ = S_.I_syn_ex_ * exp_tau_ex; + S_.I_syn_in_ = S_.I_syn_in_ * exp_tau_in; } void @@ -544,7 +548,7 @@ nest::iaf_psc_exp_ps_lossless::emit_spike_( const Time& origin, const long lag, // compute spike time relative to beginning of step S_.last_spike_step_ = origin.get_steps() + lag + 1; - S_.last_spike_offset_ = V_.h_ms_ - ( t0 + bisectioning_( dt ) ); + S_.last_spike_offset_ = V_.h_ms_ - ( t0 + regula_falsi( *this, dt ) ); // reset neuron and make it refractory S_.y2_ = P_.U_reset_; @@ -579,39 +583,18 @@ nest::iaf_psc_exp_ps_lossless::emit_instant_spike_( const Time& origin, const lo kernel().event_delivery_manager.send( *this, se, lag ); } -inline double -nest::iaf_psc_exp_ps_lossless::bisectioning_( const double dt ) const +double +nest::iaf_psc_exp_ps_lossless::threshold_distance( double t_step ) const { - double root = 0.0; - double y2_root = V_.y2_before_; - double div = 2.0; - while ( fabs( P_.U_th_ - y2_root ) > 1e-14 and ( dt / div > 0.0 ) ) - { - if ( y2_root > P_.U_th_ ) - { - root -= dt / div; - } - else - { - root += dt / div; - } - - div *= 2.0; + const double P20 = -P_.tau_m_ / P_.c_m_ * numerics::expm1( -t_step / P_.tau_m_ ); - const double expm1_tau_ex = numerics::expm1( -root / P_.tau_ex_ ); - const double expm1_tau_in = numerics::expm1( -root / P_.tau_in_ ); - const double expm1_tau_m = numerics::expm1( -root / P_.tau_m_ ); + const double P21_ex = propagator_32( P_.tau_ex_, P_.tau_m_, P_.c_m_, t_step ); + const double P21_in = propagator_32( P_.tau_in_, P_.tau_m_, P_.c_m_, t_step ); - const double P20 = -P_.tau_m_ / P_.c_m_ * expm1_tau_m; - const double P21_ex = - -P_.tau_m_ * P_.tau_ex_ / ( P_.tau_m_ - P_.tau_ex_ ) / P_.c_m_ * ( expm1_tau_ex - expm1_tau_m ); - const double P21_in = - -P_.tau_m_ * P_.tau_in_ / ( P_.tau_m_ - P_.tau_in_ ) / P_.c_m_ * ( expm1_tau_in - expm1_tau_m ); + double y2_root = P20 * ( P_.I_e_ + V_.y0_before_ ) + P21_ex * V_.I_syn_ex_before_ + P21_in * V_.I_syn_in_before_ + + V_.y2_before_ * std::exp( -t_step / P_.tau_m_ ); - y2_root = P20 * ( P_.I_e_ + V_.y0_before_ ) + P21_ex * V_.I_syn_ex_before_ + P21_in * V_.I_syn_in_before_ - + expm1_tau_m * V_.y2_before_ + V_.y2_before_; - } - return root; + return y2_root - P_.U_th_; } double diff --git a/precise/iaf_psc_exp_ps_lossless.h b/precise/iaf_psc_exp_ps_lossless.h index c6cc07cef4..f7bdbea4ac 100644 --- a/precise/iaf_psc_exp_ps_lossless.h +++ b/precise/iaf_psc_exp_ps_lossless.h @@ -151,6 +151,17 @@ class iaf_psc_exp_ps_lossless : public Archiving_Node void get_status( DictionaryDatum& ) const; void set_status( const DictionaryDatum& ); + /** + * Based on the current state, compute the value of the membrane potential + * after taking a timestep of length ``t_step``, and use it to compute the + * signed distance to spike threshold at that time. The internal state is not + * actually updated (method is defined const). + * + * @param double time step + * @returns difference between updated membrane potential and threshold + */ + double threshold_distance( double t_step ) const; + private: /** @name Interface functions * @note These functions are private, so that they can be accessed @@ -213,13 +224,6 @@ class iaf_psc_exp_ps_lossless : public Archiving_Node */ void emit_instant_spike_( const Time& origin, const long lag, const double spike_offset ); - /** - * Localize threshold crossing by bisectioning. - * @param double length of interval since previous event - * @returns time from previous event to threshold crossing - */ - double bisectioning_( const double dt ) const; - /** * Retrospective spike detection by state space analysis. * The state space spanning the non-spiking region is bound by the following @@ -336,9 +340,9 @@ class iaf_psc_exp_ps_lossless : public Archiving_Node { double h_ms_; //!< Time resolution [ms] long refractory_steps_; //!< Refractory time in steps - double expm1_tau_m_; //!< exp(-h/tau_m) - 1 - double expm1_tau_ex_; //!< exp(-h/tau_ex) - 1 - double expm1_tau_in_; //!< exp(-h/tau_in) - 1 + double exp_tau_m_; //!< exp(-h/tau_m) + double exp_tau_ex_; //!< exp(-h/tau_ex) + double exp_tau_in_; //!< exp(-h/tau_in) double P20_; //!< Progagator matrix element, 2nd row double P21_in_; //!< Progagator matrix element, 2nd row double P21_ex_; //!< Progagator matrix element, 2nd row @@ -346,8 +350,6 @@ class iaf_psc_exp_ps_lossless : public Archiving_Node double I_syn_ex_before_; //!< I_syn_ex_ at beginning of ministep double I_syn_in_before_; //!< I_syn_in_ at beginning of ministep double y2_before_; //!< y2_ at beginning of ministep - double bisection_step_; //!< if missed spike is detected, - //!< calculate time to emit spike /** * Pre-computed constants for inequality V < g(h, I_e) diff --git a/testsuite/cpptests/run_all.cpp b/testsuite/cpptests/run_all.cpp index 828bcd2df2..b4161ea531 100644 --- a/testsuite/cpptests/run_all.cpp +++ b/testsuite/cpptests/run_all.cpp @@ -22,7 +22,7 @@ #define BOOST_TEST_MODULE cpptests #define BOOST_TEST_DYN_LINK -#include +#include // Includes from cpptests #include "test_block_vector.h" diff --git a/testsuite/cpptests/test_block_vector.h b/testsuite/cpptests/test_block_vector.h index 3fec9ef39f..83504641fd 100644 --- a/testsuite/cpptests/test_block_vector.h +++ b/testsuite/cpptests/test_block_vector.h @@ -32,7 +32,31 @@ // Includes from libnestutil: #include "block_vector.h" -BOOST_AUTO_TEST_SUITE( test_seque ) +/** + * Fixture filling a BlockVector and a vector with linearly increasing values. + */ +struct bv_vec_reference_fixture +{ + bv_vec_reference_fixture() + : N( block_vector.get_max_block_size() + 10 ) + { + for ( int i = 0; i < N; ++i ) + { + block_vector.push_back( i ); + reference.push_back( i ); + } + } + + ~bv_vec_reference_fixture() + { + } + + BlockVector< int > block_vector; + std::vector< int > reference; + int N; +}; + +BOOST_AUTO_TEST_SUITE( test_blockvector ) BOOST_AUTO_TEST_CASE( test_size ) { @@ -92,27 +116,27 @@ BOOST_AUTO_TEST_CASE( test_erase ) block_vector.push_back( i ); } - auto seque_mid = block_vector; - seque_mid.erase( seque_mid.begin() + 2, seque_mid.begin() + 8 ); - BOOST_REQUIRE( seque_mid.size() == 4 ); - BOOST_REQUIRE( seque_mid[ 0 ] == 0 ); - BOOST_REQUIRE( seque_mid[ 1 ] == 1 ); - BOOST_REQUIRE( seque_mid[ 2 ] == 8 ); - BOOST_REQUIRE( seque_mid[ 3 ] == 9 ); - - auto seque_front = block_vector; - seque_front.erase( seque_front.begin(), seque_front.begin() + 7 ); - BOOST_REQUIRE( seque_front.size() == 3 ); - BOOST_REQUIRE( seque_front[ 0 ] == 7 ); - BOOST_REQUIRE( seque_front[ 1 ] == 8 ); - BOOST_REQUIRE( seque_front[ 2 ] == 9 ); - - auto seque_back = block_vector; - seque_back.erase( seque_back.begin() + 3, seque_back.end() ); - BOOST_REQUIRE( seque_back.size() == 3 ); - BOOST_REQUIRE( seque_back[ 0 ] == 0 ); - BOOST_REQUIRE( seque_back[ 1 ] == 1 ); - BOOST_REQUIRE( seque_back[ 2 ] == 2 ); + auto bv_mid = block_vector; + bv_mid.erase( bv_mid.begin() + 2, bv_mid.begin() + 8 ); + BOOST_REQUIRE( bv_mid.size() == 4 ); + BOOST_REQUIRE( bv_mid[ 0 ] == 0 ); + BOOST_REQUIRE( bv_mid[ 1 ] == 1 ); + BOOST_REQUIRE( bv_mid[ 2 ] == 8 ); + BOOST_REQUIRE( bv_mid[ 3 ] == 9 ); + + auto bv_front = block_vector; + bv_front.erase( bv_front.begin(), bv_front.begin() + 7 ); + BOOST_REQUIRE( bv_front.size() == 3 ); + BOOST_REQUIRE( bv_front[ 0 ] == 7 ); + BOOST_REQUIRE( bv_front[ 1 ] == 8 ); + BOOST_REQUIRE( bv_front[ 2 ] == 9 ); + + auto bv_back = block_vector; + bv_back.erase( bv_back.begin() + 3, bv_back.end() ); + BOOST_REQUIRE( bv_back.size() == 3 ); + BOOST_REQUIRE( bv_back[ 0 ] == 0 ); + BOOST_REQUIRE( bv_back[ 1 ] == 1 ); + BOOST_REQUIRE( bv_back[ 2 ] == 2 ); } BOOST_AUTO_TEST_CASE( test_begin ) @@ -203,10 +227,10 @@ BOOST_AUTO_TEST_CASE( test_iterator_dereference ) BOOST_REQUIRE( *( block_vector.begin() ) == block_vector[ 0 ] ); // Using operator-> - BlockVector< std::vector< int > > nested_seque; + BlockVector< std::vector< int > > nested_bv; std::vector< int > tmp = { 42 }; - nested_seque.push_back( tmp ); - BOOST_REQUIRE( nested_seque.begin()->size() == 1 ); + nested_bv.push_back( tmp ); + BOOST_REQUIRE( nested_bv.begin()->size() == 1 ); } BOOST_AUTO_TEST_CASE( test_iterator_assign ) @@ -244,16 +268,145 @@ BOOST_AUTO_TEST_CASE( test_iterator_compare ) } BOOST_REQUIRE( block_vector.begin() < block_vector.end() ); - auto it_a = block_vector.begin(); - auto it_b = block_vector.begin(); - BOOST_REQUIRE( it_a == it_b ); + // Test comparison with iterator shifted one step, shifted to the end of + // the block, and shifted to the next block. + std::vector< int > it_shifts = { 1, block_vector.get_max_block_size() - 1, N - 1 }; + for ( auto& shift : it_shifts ) + { + auto it_a = block_vector.begin(); + auto it_b = block_vector.begin(); + BOOST_REQUIRE( it_a == it_b ); + BOOST_REQUIRE( not( it_a != it_b ) ); + + it_b += shift; + + BOOST_REQUIRE( it_a != it_b ); + BOOST_REQUIRE( it_a < it_b ); + BOOST_REQUIRE( it_a <= it_b ); + BOOST_REQUIRE( it_b > it_a ); + BOOST_REQUIRE( it_b >= it_a ); + + BOOST_REQUIRE( not( it_a == it_b ) ); + BOOST_REQUIRE( not( it_b < it_a ) ); + BOOST_REQUIRE( not( it_b <= it_a ) ); + BOOST_REQUIRE( not( it_a > it_b ) ); + BOOST_REQUIRE( not( it_a >= it_b ) ); + } +} - ++it_b; - BOOST_REQUIRE( it_a != it_b ); - BOOST_REQUIRE( it_a < it_b ); - BOOST_REQUIRE( not( it_b < it_a ) ); +BOOST_FIXTURE_TEST_CASE( test_operator_pp, bv_vec_reference_fixture ) +{ + // operator++() + auto bvi = block_vector.begin(); + for ( auto& ref : reference ) + { + BOOST_REQUIRE( *bvi == ref ); + ++bvi; + } + BOOST_REQUIRE( bvi == block_vector.end() ); +} + +BOOST_FIXTURE_TEST_CASE( test_operator_p, bv_vec_reference_fixture ) +{ + // operator+() + for ( int i = 0; i < N; ++i ) + { + auto bvi = block_vector.begin(); + auto ref_it = reference.begin(); + auto new_bvi = bvi + i; + auto new_ref_it = ref_it + i; + BOOST_REQUIRE( *new_bvi == *new_ref_it ); + } +} + +BOOST_FIXTURE_TEST_CASE( test_operator_p_eq, bv_vec_reference_fixture ) +{ + // operator+=() + for ( int i = 0; i < N; ++i ) + { + auto bvi = block_vector.begin(); + auto bvi_last = block_vector.end() - 1; + auto ref_it = reference.begin(); + auto ref_it_last = reference.end() - 1; + bvi += i; + ref_it += i; + BOOST_REQUIRE( *bvi == *ref_it ); + bvi_last += -i; + ref_it_last += -i; + BOOST_REQUIRE( *bvi_last == *ref_it_last ); + } +} + +BOOST_FIXTURE_TEST_CASE( test_operator_m_eq, bv_vec_reference_fixture ) +{ + // operator-=() + for ( int i = 1; i < N - 1; ++i ) + { + auto bvi = block_vector.end(); + auto bvi_first = block_vector.begin(); + auto ref_it = reference.end(); + auto ref_it_first = reference.begin(); + bvi -= i; + ref_it -= i; + BOOST_REQUIRE( *bvi == *ref_it ); + bvi_first -= -( i - 1 ); + ref_it_first -= -( i - 1 ); + BOOST_REQUIRE( *bvi_first == *ref_it_first ); + } +} + +BOOST_FIXTURE_TEST_CASE( test_operator_m, bv_vec_reference_fixture ) +{ + // operator-() + for ( int i = 1; i < N - 1; ++i ) + { + auto bvi = block_vector.end(); + auto ref_it = reference.end(); + auto new_bvi = bvi - i; + auto new_ref_it = ref_it - i; + BOOST_REQUIRE( *new_bvi == *new_ref_it ); + } +} + +BOOST_FIXTURE_TEST_CASE( test_operator_mm, bv_vec_reference_fixture ) +{ + // operator--() + auto bvi = block_vector.end() - 1; + for ( auto ref_it = reference.end() - 1; ref_it >= reference.begin(); --ref_it, --bvi ) + { + BOOST_REQUIRE( *bvi == *ref_it ); + } + BOOST_REQUIRE( bvi == block_vector.begin() - 1 ); +} + +BOOST_FIXTURE_TEST_CASE( test_operator_eq, bv_vec_reference_fixture ) +{ + // operator==() + auto bvi_pp = block_vector.begin() + 1; + auto bvi_copy = block_vector.begin(); + auto bvi_mm = block_vector.begin() - 1; + for ( auto bvi = block_vector.begin(); bvi != block_vector.end(); ++bvi, ++bvi_copy, ++bvi_mm, ++bvi_pp ) + { + BOOST_REQUIRE( bvi == bvi_copy ); + BOOST_REQUIRE( not( bvi == bvi_pp ) ); + BOOST_REQUIRE( not( bvi == bvi_mm ) ); + } +} + +BOOST_FIXTURE_TEST_CASE( test_operator_neq, bv_vec_reference_fixture ) +{ + // operator!=() + auto bvi_pp = block_vector.begin() + 1; + auto bvi_copy = block_vector.begin(); + auto bvi_mm = block_vector.begin() - 1; + for ( auto bvi = block_vector.begin(); bvi != block_vector.end(); ++bvi, ++bvi_copy, ++bvi_mm, ++bvi_pp ) + { + BOOST_REQUIRE( not( bvi != bvi_copy ) ); + BOOST_REQUIRE( bvi != bvi_pp ); + BOOST_REQUIRE( bvi != bvi_mm ); + } } BOOST_AUTO_TEST_SUITE_END() -#endif /* TEST_SORT_H */ +#endif /* TEST_BLOCK_VECTOR_H */ diff --git a/testsuite/cpptests/test_sort.h b/testsuite/cpptests/test_sort.h index 75ecd5ec4f..67ebcf1a65 100644 --- a/testsuite/cpptests/test_sort.h +++ b/testsuite/cpptests/test_sort.h @@ -27,82 +27,211 @@ #include // C++ includes: +#include #include // Includes from libnestutil: #include "sort.h" -namespace nest -{ - +/** + * Wrapper for quicksort3way. + * + * When calling nest::sort() directly it is impossible to sort with the + * built in quicksort3way from tests. This is because if NEST is compiled + * with Boost support, nest::sort() sorts with Boost, and if NEST is not + * compiled with Boost, nest::sort() calls quicksort3way, but the Boost + * tests are not compiled. + */ void -nest_quicksort( BlockVector< size_t >& bv0, BlockVector< size_t >& bv1 ) +nest_quicksort( BlockVector< int >& bv0, BlockVector< int >& bv1 ) { nest::quicksort3way( bv0, bv1, 0, bv0.size() - 1 ); } -const bool -is_sorted( BlockVector< size_t >::const_iterator begin, BlockVector< size_t >::const_iterator end ) +/** + * Fixture filling two BlockVectors and a vector with lineary decreasing numbers. + * The vector is then sorted. + */ +struct fill_bv_vec_linear { - for ( BlockVector< size_t >::const_iterator it = begin; it < --end; ) + fill_bv_vec_linear() + : N( 20000 ) + , N_small( boost::sort::spreadsort::detail::min_sort_size - 10 ) + , bv_sort( N ) + , bv_perm( N ) + , vec_sort( N ) + , bv_sort_small( N_small ) + , bv_perm_small( N_small ) + , vec_sort_small( N_small ) { - if ( *it > *( ++it ) ) + for ( int i = 0; i < N; ++i ) { - return false; + int element = N - i; + bv_sort[ i ] = element; + bv_perm[ i ] = element; + vec_sort[ i ] = element; + + if ( i < N_small ) + { + bv_sort_small[ i ] = element; + bv_perm_small[ i ] = element; + vec_sort_small[ i ] = element; + } } + std::sort( vec_sort.begin(), vec_sort.end() ); + std::sort( vec_sort_small.begin(), vec_sort_small.end() ); } - return true; -} + + const int N; + const int N_small; + + BlockVector< int > bv_sort; + BlockVector< int > bv_perm; + std::vector< int > vec_sort; + + BlockVector< int > bv_sort_small; + BlockVector< int > bv_perm_small; + std::vector< int > vec_sort_small; +}; + +/** + * Fixture filling two BlockVectors and a vector with random numbers. + * The vector is then sorted. + */ +struct fill_bv_vec_random +{ + fill_bv_vec_random() + : N( 20000 ) + , N_small( boost::sort::spreadsort::detail::min_sort_size - 10 ) + , bv_sort( N ) + , bv_perm( N ) + , vec_sort( N ) + , bv_sort_small( N_small ) + , bv_perm_small( N_small ) + , vec_sort_small( N_small ) + { + for ( int i = 0; i < N; ++i ) + { + const size_t k = std::rand() % N; + bv_sort[ i ] = k; + bv_perm[ i ] = k; + vec_sort[ i ] = k; + + if ( i < N_small ) + { + bv_sort_small[ i ] = k; + bv_perm_small[ i ] = k; + vec_sort_small[ i ] = k; + } + } + std::sort( vec_sort.begin(), vec_sort.end() ); + std::sort( vec_sort_small.begin(), vec_sort_small.end() ); + } + + const int N; + const int N_small; + + BlockVector< int > bv_sort; + BlockVector< int > bv_perm; + std::vector< int > vec_sort; + + BlockVector< int > bv_sort_small; + BlockVector< int > bv_perm_small; + std::vector< int > vec_sort_small; +}; BOOST_AUTO_TEST_SUITE( test_sort ) /** * Tests whether two arrays with randomly generated numbers are sorted - * correctly by a single call to sort. + * correctly when sorting with NEST's own quicksort. + */ +BOOST_FIXTURE_TEST_CASE( test_quicksort_random, fill_bv_vec_random ) +{ + nest_quicksort( bv_sort, bv_perm ); + + // We test here that the BlockVectors bv_sort and bv_perm are both + // sorted here. However, if something goes wrong and they for example + // contain only zeros, is_sorted will also return true, but the test + // should not pass. Therefore, in addition to require the BlockVectors + // to be sorted, we also require them to be equal to the sorted + // reference vector, vec_sort. + BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) ); + BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) ); + + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) ); + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) ); +} + +/** + * Tests whether two arrays with linearly decreasing numbers are sorted + * correctly when sorting with the built-in quicksort. + */ +BOOST_FIXTURE_TEST_CASE( test_quicksort_linear, fill_bv_vec_linear ) +{ + nest_quicksort( bv_sort, bv_perm ); + + BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) ); + BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) ); + + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) ); + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) ); +} + +/** + * Tests whether two arrays with randomly generated numbers are sorted + * correctly when sorting with Boost. */ -BOOST_AUTO_TEST_CASE( test_random ) +BOOST_FIXTURE_TEST_CASE( test_boost_random, fill_bv_vec_random ) { - const size_t N = 20000; - BlockVector< size_t > bv0( N ); - BlockVector< size_t > bv1( N ); + // Making sure we are sorting with boost + static_assert( HAVE_BOOST, "Compiling Boost tests, but HAVE_BOOST!=1." ); - for ( size_t i = 0; i < N; ++i ) - { - const size_t k = std::rand() % N; - bv0[ i ] = k; - bv1[ i ] = k; - } + nest::sort( bv_sort, bv_perm ); + + BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) ); + BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) ); - nest_quicksort( bv0, bv1 ); + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) ); + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) ); - BOOST_REQUIRE( is_sorted( bv0.begin(), bv0.end() ) ); - BOOST_REQUIRE( is_sorted( bv1.begin(), bv1.end() ) ); + // Using smaller data sets to sort with the fallback algorithm. + nest::sort( bv_sort_small, bv_perm_small ); + + BOOST_REQUIRE( std::is_sorted( bv_sort_small.begin(), bv_sort_small.end() ) ); + BOOST_REQUIRE( std::is_sorted( bv_perm_small.begin(), bv_perm_small.end() ) ); + + BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_sort_small.begin() ) ); + BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_perm_small.begin() ) ); } /** * Tests whether two arrays with linearly increasing numbers are sorted - * correctly by a single call to sort. + * correctly when sorting with Boost. */ -BOOST_AUTO_TEST_CASE( test_linear ) +BOOST_FIXTURE_TEST_CASE( test_boost_linear, fill_bv_vec_linear ) { - const size_t N = 20000; - BlockVector< size_t > bv0( N ); - BlockVector< size_t > bv1( N ); + // Making sure we are sorting with boost + static_assert( HAVE_BOOST, "Compiling Boost tests, but HAVE_BOOST!=1." ); - for ( size_t i = 0; i < N; ++i ) - { - bv0[ i ] = N - i - 1; - bv1[ i ] = N - i - 1; - } + nest::sort( bv_sort, bv_perm ); + + BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) ); + BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) ); + + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) ); + BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) ); - nest_quicksort( bv0, bv1 ); + // Using smaller data sets to sort with the fallback algorithm. + nest::sort( bv_sort_small, bv_perm_small ); - BOOST_REQUIRE( is_sorted( bv0.begin(), bv0.end() ) ); - BOOST_REQUIRE( is_sorted( bv1.begin(), bv1.end() ) ); + BOOST_REQUIRE( std::is_sorted( bv_sort_small.begin(), bv_sort_small.end() ) ); + BOOST_REQUIRE( std::is_sorted( bv_perm_small.begin(), bv_perm_small.end() ) ); + + BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_sort_small.begin() ) ); + BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_perm_small.begin() ) ); } BOOST_AUTO_TEST_SUITE_END() -} // of namespace nest - #endif /* TEST_SORT_H */ diff --git a/testsuite/unittests/test_iaf_ps_dc_accuracy.sli b/testsuite/unittests/test_iaf_ps_dc_accuracy.sli index e55321912b..5ad56500b3 100644 --- a/testsuite/unittests/test_iaf_ps_dc_accuracy.sli +++ b/testsuite/unittests/test_iaf_ps_dc_accuracy.sli @@ -41,7 +41,8 @@ Description: and an appropriate arrangement of the terms [2]. For small computation step sizes the accuracy at large simulation time decreases because of the accumulation of errors. - The expected output is documented at the end of the script. + Reference output is documented at the end of the script. + Individual simulation results can be inspected by uncommented the call to function print_details. @@ -65,7 +66,7 @@ References: [2] Morrison A, Straube S, Plesser H E, & Diesmann M (2007) Exact Subthreshold Integration with Continuous Spike Times in Discrete Time Neural Network Simulations. Neural Computation 19:47--79 -SeeAlso: iaf_psc_alpha_ps, iaf_psc_alpha_presc, iaf_psc_delta_ps, testsuite::test_iaf_ps_dc_t_accuracy +SeeAlso: iaf_psc_alpha_ps, iaf_psc_delta_ps, testsuite::test_iaf_ps_dc_t_accuracy */ (unittest) run @@ -83,10 +84,17 @@ M_ERROR setverbosity [0 min_exponent -1] Range /hlist Set [ % time [ms] tolerated error [mV] - [ 5.0 1e-13 ] - [ 500.0 1e-9 ] % error larger because of accumulation + [ 5.0 1e-13 ] + [ 500.0 1e-9 ] % error larger because of accumulation ] /Tlist Set % at very small computation step sizes +% Models to be tested by this test +[ + /iaf_psc_alpha_ps + /iaf_psc_delta_ps + % other precise models should be added to this list +] + /models Set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -94,8 +102,8 @@ M_ERROR setverbosity % Check if kernel accepts high precision % 0 << - /tics_per_ms min_exponent neg dexp - /resolution 0 dexp % 1 ms default + /tics_per_ms min_exponent neg dexp + /resolution 0 dexp % 1 ms default >> SetStatus @@ -119,8 +127,8 @@ exch /C_m 250.0 % membrane capacity in pF >> /params Set -params begin userdict begin +params begin userdict begin %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % @@ -130,69 +138,88 @@ params begin userdict begin % /SimAtResolution { - dup /i Set - dexp /h Set - - ResetKernel % resets tic base and computation time step - 0 << /tics_per_ms min_exponent neg dexp /resolution h >> SetStatus - - - [ - /iaf_psc_alpha_ps - /iaf_psc_alpha_presc - /iaf_psc_delta_ps % this list can be extended - ] - {Create dup params SetStatus} Map /neurons Set + dup /i Set + dexp /h Set - T Simulate + ResetKernel % resets tic base and computation time step + 0 << /tics_per_ms min_exponent neg dexp /resolution h >> SetStatus - neurons - {[exch /V_m get dup V sub abs]} Map Flatten + models {Create dup params SetStatus} Map /neurons Set - i prepend + T Simulate + neurons + { [ exch /V_m get dup V sub abs] } Map Flatten + i prepend } def +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% helper function for debugging, +% prints detailed table of results +% +/print_details +{ + cout default 15 setprecision + + endl endl + (Exact value of membrane potential after ) <- + T <- ( ms is ) <- V <- ( mV.) <- endl endl + + ( log_2 h) <- + models + { + exch ( ) <- exch <- ( [mV]) <- + ( error [mV]) <- + } forall + endl + + models length 2 mul 1 add + { + (--------------------------) <- + } repeat + endl + + exch + { + { exch 24 setw exch <- ( ) <- } forall endl + } + forall + ; +} +def %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% do all simulation times +% Perform test for different simulation durations % { - Tlist - { - [/T /tolerance] Set - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % - % Reference value - % - % V is the exact value of the membrane potential at the end - % of the simulation. - % - - (I_e * tau_m/C_m * (1. - exp(-T/tau_m)) ) ExecMath /V Set - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % - % Perform simulations at all resolutions and collect results - % - hlist {SimAtResolution} Map - - - %dup print_details % uncomment for debugging - - {Rest 2 Partition [/All 2] Part} Map % select columns containing - % the voltage errors - Flatten {tolerance lt} Map - } Map - - Flatten true exch {and} Fold + Tlist + { + [/T /tolerance] Set + + % Reference value + (I_e * tau_m/C_m * (1. - exp(-T/tau_m)) ) ExecMath /V Set + + % Simulate at different resolutions + hlist {SimAtResolution} Map + + dup print_details + + % select columns with voltage errors + { Rest 2 Partition [/All 2] Part } Map + + % check against tolerance + Flatten {tolerance lt} Map + } + Map + + % combine results + Flatten true exch {and} Fold } - +assert_or_die %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -317,65 +344,3 @@ params begin userdict begin % -12 39.9999999998544 1.45590206557245e-10 39.9999999998544 1.45590206557245e-10 % -13 39.9999999997088 2.91180413114489e-10 39.9999999997088 2.91180413114489e-10 % -14 39.9999999994176 5.82360826228978e-10 39.9999999994176 5.82360826228978e-10 -% -% - - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% helper function for debugging, -% prints detailed table of results -% -/print_details -{ -cout default 15 setprecision - -endl -endl -(Exact value of membrane potential after ) <- -T <- ( ms is ) <- -V <- ( mV.) <- endl - -endl - -( h in ms ) <- -( alpha_canon [mV]) <- -( error [mV]) <- -( alpha_presc [mV]) <- -( error [mV]) <- -( delta_canon [mV]) <- -( error [mV]) <- -endl -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(------------------------) <- -endl - -exch -{ - { - exch 24 setw exch <- ( ) <- - } - forall - endl -} -forall -; -} -def - - - -% executes the overall test -assert_or_die - - - - - diff --git a/testsuite/unittests/test_iaf_ps_dc_t_accuracy.sli b/testsuite/unittests/test_iaf_ps_dc_t_accuracy.sli index bca2bda483..4fee7479ba 100644 --- a/testsuite/unittests/test_iaf_ps_dc_t_accuracy.sli +++ b/testsuite/unittests/test_iaf_ps_dc_t_accuracy.sli @@ -30,7 +30,7 @@ Description: A DC current is injected for a finite duration. The time of the first spike is compared to the theoretical value for different computation - step sizes and interpolation orders. + step sizes. Computation step sizes are specified as base 2 values. @@ -63,7 +63,7 @@ References: [2] Morrison A, Straube S, Plesser H E, & Diesmann M (2007) Exact Subthreshold Integration with Continuous Spike Times in Discrete Time Neural Network Simulations. Neural Computation 19:47--79 -SeeAlso: iaf_psc_alpha_ps, iaf_psc_alpha_presc, iaf_psc_delta_ps, testsuite::test_iaf_ps_dc_accuracy +SeeAlso: iaf_psc_alpha_ps, iaf_psc_delta_ps, testsuite::test_iaf_ps_dc_accuracy */ (unittest) run @@ -71,22 +71,23 @@ SeeAlso: iaf_psc_alpha_ps, iaf_psc_alpha_presc, iaf_psc_delta_ps, testsuite::tes M_ERROR setverbosity --14 /min_exponent Set - -[0 min_exponent -1] Range /hlist Set - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Parameters of simulation schedule. +% Test parameters % -5.0 /T Set -[ % interpolation order tolerated error [mv] - [ 0 1e-4 ] - [ 1 1e-10 ] - [ 2 1e-13 ] - [ 3 1e-13 ] -] /Olist Set + +-14 /min_exponent Set +[0 min_exponent -1] Range /hlist Set +5.0 /T Set +1e-13 /tolerance Set % tolerated error [mv] + +% models to be tested +[ + /iaf_psc_alpha_ps + /iaf_psc_delta_ps + % other precise models should be tested as well + ] + /models Set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -106,13 +107,11 @@ exch {1 sub min_exponent leq} assert_or_die % sufficient resolution? - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Parameters of neuron model. % - << /E_L 0.0 % resting potential in mV /V_m 0.0 % initial membrane potential in mV @@ -133,126 +132,78 @@ params begin userdict begin % /SimAtResolution { - dup /i Set - dexp /h Set - + dup /i Set + dexp /h Set - ResetKernel - 0 << /tics_per_ms min_exponent neg dexp /resolution h >> SetStatus - - [ - /iaf_psc_alpha_ps - /iaf_psc_alpha_presc - /iaf_psc_delta_ps % this list can be extended - ] - {Create dup params SetStatus} Map /neurons Set - - neurons {/Interpol_Order known} Select {<< /Interpol_Order O >> SetStatus} forall + ResetKernel + 0 << /tics_per_ms min_exponent neg dexp /resolution h >> SetStatus + models {Create dup params SetStatus} Map /neurons Set - T Simulate - - neurons - {[exch /t_spike get dup t sub abs]} Map Flatten - - i prepend + T Simulate + neurons + { [ exch /t_spike get dup t sub abs ] } Map Flatten + i prepend } def -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% do for all interpolation orders +% prints detailed table of results % +/print_details { - Olist - { - [/O /tolerance] Set - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % - % Reference value - % - % V is the exact value of the membrane potential at the end - % of the simulation. - % - - % tau_m neg 1 c_m V_th mul tau_m I_e mul div sub ln mul - - (-tau_m*ln( 1.0 - (C_m*V_th)/(tau_m*I_e) )) ExecMath /t Set - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % - % Perform simulations at all resolutions and collect results - % - hlist {SimAtResolution} Map - - - dup print_details % uncomment for debugging - - {Rest 2 Partition [/All 2] Part} Map % select columns containing - % the timing errors - [-1] Part % select only highest resolution - Flatten {tolerance lt} Map - } Map - - Flatten true exch {and} Fold + cout default 15 setprecision + + endl endl + (The precise spike time is ) <- t <- ( ms.) <- endl endl + + ( log_2 h) <- + models + { + exch ( ) <- exch <- ( [ms]) <- + ( error [ms]) <- + } forall + endl + + models length 2 mul 1 add + { + (--------------------------) <- + } repeat + endl + + exch + { + { exch 24 setw exch <- ( ) <- } forall endl + } + forall + ; } +def - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% helper function for debugging, -% prints detailed table of results +% Perform test % -/print_details { -cout default 15 setprecision - -endl -endl -(order: ) <- O <- (, exact value of membrane potential after ) <- -T <- ( ms is ) <- -t <- ( ms.) <- endl - -endl - -( h in ms ) <- -( alpha_canon [ms]) <- -( error [ms]) <- -( alpha_presc [ms]) <- -( error [ms]) <- -( delta_canon [ms]) <- -( error [ms]) <- -endl -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(--------------------------) <- -(------------------------) <- -endl + % Reference value + (-tau_m*ln( 1.0 - (C_m*V_th)/(tau_m*I_e) )) ExecMath /t Set -exch -{ - { - exch 24 setw exch <- ( ) <- - } - forall - endl -} -forall -; -} -def + % Perform simulations across resolutions + hlist { SimAtResolution } Map -exec + dup print_details + % select columns with timing errors, highest resolution only + { Rest 2 Partition [/All 2] Part } Map [-1] Part + + % test against tolerance limit + Flatten { tolerance lt } Map -% executes the overall test + % combine results + Flatten true exch {and} Fold +} assert_or_die