Skip to content

Instantly share code, notes, and snippets.

@heplesser
Created March 2, 2016 07:48
Show Gist options
  • Select an option

  • Save heplesser/e78bed61478c25f7047c to your computer and use it in GitHub Desktop.

Select an option

Save heplesser/e78bed61478c25f7047c to your computer and use it in GitHub Desktop.
First step towards mt-invariant hpc_benchmark.sli
/*
* hpc_benchmark.sli
*
* 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 <http://www.gnu.org/licenses/>.
*
*/
/*
This script produces a balanced random network of scale*11250 neurons in
which the excitatory-excitatory neurons exhibit STDP with
multiplicative depression and power-law potentiation. A mutual
equilibrium is obtained between the activity dynamics (low rate in
asynchronous irregular regime) and the synaptic weight distribution
(unimodal).
This is the standard network investigated in:
Morrison et al (2007). Spike-timing-dependent plasticity in balanced random networks Neural Comput 19(6):1437-67
Helias et al (2012). Supercomputers ready for use as discovery machines for neuroscience Front. Neuroinform. 6:26,
Kunkel et al (2014). Spiking network simulation code for petascale computers. Front. Neuroinform. 8:78
*/
%%% PARAMETER SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define all relevant parameters: changes should be made here
% all data is placed in the userdict dictionary
/nvp 16 def % total number of virtual processes
/scale 1.0 def % scaling factor of the network size, total network size = scale*11250 neurons
/simtime 1000. def % total simulation time in ms
/presimtime 100. ms def % simulation time until reaching equilibrium
/dt 0.1 ms def % simulation step
/record_spikes true def % switch to record spikes of excitatory neurons to file
/path_name (.) def % path where all files will have to be written
/log_file (log) def % naming scheme for the log files
% -------------------------------------------------------------------------------------
%%% Define function to convert synapse weight from mV to pA %%%%%%%%%%
/ConvertSynapseWeight
{
/tauMem Set
/tauSyn Set
/CMem Set
% This function is specific to the used neuron model
% Leaky integrate-and-fire neuron with alpha-shaped
% postsynaptic currents
(
a = tauMem / tauSyn;
b = 1.0 / tauSyn - 1.0 / tauMem;
% time of maximum
t_max = 1.0/b * (-LambertWm1(-exp(-1.0/a)/a)-1.0/a);
% maximum of PSP for current of unit amplitude
exp(1.0)/(tauSyn*CMem*b) * ((exp(-t_max/tauMem) - exp(-t_max/tauSyn)) / b - t_max*exp(-t_max/tauSyn))
) ExecMath
1 exch div
}
def
%%% Rise time of synaptic currents
% The synaptic rise time is chosen such that the rise time of the
% evoked post-synaptic potential is 1.700759 ms.
% For alpha-shaped postynaptic currents, the rise time of the post-synaptic
% potential follows from the synaptic rise time as
/PSP_rise_time
{
/taumem set
/tausyn set
(
a = taumem / tausyn;
b = 1.0 / tausyn - 1.0 / taumem;
) execmath
} def
% Inverting this equation numerically leads to tau_syn = 0.32582722403722841 ms, as specified in model_params below.
% -------------------------------------------------------------------------------------
/brunel_params
<<
/NE 9000 scale mul round cvi % number of excitatory neurons
/NI 2250 scale mul round cvi % number of inhibitory neurons
% number of neurons to record spikes from
/Nrec_total 1000
/model_params
<<
% Set variables for iaf_neuron
/E_L 0.0 mV % Resting membrane potential [mV]
/C_m 250.0 pF % Capacity of the membrane [pF]
/tau_m 10.0 ms % Membrane time constant [ms]
/t_ref 0.5 ms % duration of refractory period [ms]
/V_th 20.0 mV % Threshold [mV]
/V_reset 0.0 mV % Reset Potential [mV]
/tau_syn 0.32582722403722841 ms % time const. postsynaptic excitatory currents [ms]
/tau_minus 30.0 ms %time constant for STDP (depression)
% V can be randomly initialized see below
/V_m 5.7 mV % mean value of membrane potential
>>
/randomize_Vm false
/mean_potential 5.7 mV % Note that Kunkel et al. (2014) report different values. The values
/sigma_potential 7.2 mV % in the paper were used for the benchmarks on K, the values given
% here were used for the benchmark on JUQUEEN.
/delay 1.5 ms % synaptic delay, all connections [ms]
% synaptic weight
/JE 0.14 mV %peak of EPSP
/sigma_w 3.47 pA %standard dev. of E->E synapses [pA]
/g -5.0
/stdp_params
<<
/delay 1.5 ms
/alpha 0.0513
/lambda 0.1 %STDP step size
/mu 0.4 %STDP weight dependence exponent (potentiation)
/tau_plus 15.0 %time constant for potentiation
>>
/eta 1.685 %scaling of external stimulus
/filestem path_name
>> def
% Here we resolve parameter dependencies, by making the independent
% values visible
brunel_params dup using
%%% FUNCTION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/BuildNetwork
{
tic % start timer on construction
% set global kernel parameters
0
<<
/total_num_virtual_procs nvp
/resolution dt
/overwrite_files true
>> SetStatus
/iaf_neuron model_params SetDefaults
M_INFO (BuildNetwork)
(Creating excitatory population.) message % show message
/iaf_neuron [NE] LayoutNetwork /E_net Set % subnet gets own gid
/E_from E_net 1 add def % gid of first child
/E_to 0 GetStatus /network_size get 1 sub def % gid of last child
M_INFO (BuildNetwork)
(Creating inhibitory population.) message % show message
/iaf_neuron [NI] LayoutNetwork /I_net Set % subnet gets own gid
/I_from I_net 1 add def % gid of first child
/I_to 0 GetStatus /network_size get 1 sub def % gid of last child
randomize_Vm
{
M_INFO (BuildNetwork)
(Randomzing membrane potentials.) message
0 GetStatus /rng_seeds get Last 1 add 0 /vp get add
rngdict/MT19937 :: exch CreateRNG
rdevdict/normal :: CreateRDV /normal Set
[E_from E_to] Range
{
<< /V_m mean_potential sigma_potential normal Random mul add >>
SetStatus
} forall
[I_from I_to] Range
{
<< /V_m mean_potential sigma_potential normal Random mul add >>
SetStatus
} forall
} if
/E_neurons E_from E_to cvgidcollection def % get memory efficient gid-collection
/I_neurons I_from I_to cvgidcollection def % get memory efficient gid-collection
/N E_neurons size I_neurons size add def
/CE NE scale cvd div iround def %number of incoming excitatory connections
/CI NI scale cvd div iround def %number of incomining inhibitory connections
M_INFO (BuildNetwork)
(Creating excitatory stimulus generator.) message
model_params using
% Convert synapse weight from mV to pA
C_m tau_syn tau_m ConvertSynapseWeight /conversion_factor Set
/JE_pA conversion_factor JE mul def
/nu_thresh V_th CE tau_m C_m div mul JE_pA mul 1.0 exp mul tau_syn mul div def
/nu_ext nu_thresh eta mul def
endusing
/E_stimulus /poisson_generator Create def
E_stimulus
<<
/rate nu_ext CE mul 1000. mul
>> SetStatus
/I_stimulus /poisson_generator Create def
I_stimulus
<<
/rate nu_ext CE mul 1000. mul
>> SetStatus
M_INFO (BuildNetwork)
(Creating excitatory spike detector.) message
record_spikes true eq {
/detector_label filestem (/alpha_) join stdp_params /alpha get cvs join (_spikes) join def
/E_detector /spike_detector Create def
E_detector
<<
/withtime true
/to_file true
/label detector_label
>> SetStatus
} if
memory_thisjob cvs ( # virt_mem_after_nodes) join logger /log call
toc /BuildNodeTime Set
mpi_rank cvs (] ) join (BuildNode time : ) join =only BuildNodeTime =only ( s) =
BuildNodeTime cvs ( # build_time_nodes) join logger /log call
tic
% Create custom synapse types with appropriate values for
% our excitatory and inhibitory connections
/static_synapse_hpc << /delay delay >> SetDefaults
/static_synapse_hpc /syn_std CopyModel
/static_synapse_hpc /syn_ex << /weight JE_pA >> CopyModel
/static_synapse_hpc /syn_in << /weight JE_pA g mul >> CopyModel
stdp_params /weight JE_pA put
/stdp_pl_synapse_hom_hpc stdp_params SetDefaults
M_INFO (BuildNetwork)
(Connecting stimulus generators.) message
% Connect Poisson generator to neuron
E_stimulus E_stimulus cvgidcollection E_neurons << /rule (all_to_all) >> << /model /syn_ex >> Connect
E_stimulus E_stimulus cvgidcollection I_neurons << /rule (all_to_all) >> << /model /syn_ex >> Connect
M_INFO (BuildNetwork)
(Connecting excitatory -> excitatory population.) message
E_neurons E_neurons << /rule (fixed_indegree) /indegree CE /autapses false /multapses true >> << /model /stdp_pl_synapse_hom_hpc >> Connect
M_INFO (BuildNetwork)
(Connecting inhibitory -> excitatory population.) message
I_neurons E_neurons << /rule (fixed_indegree) /indegree CI /autapses false /multapses true >> << /model /syn_in >> Connect
M_INFO (BuildNetwork)
(Connecting excitatory -> inhibitory population.) message
E_neurons I_neurons << /rule (fixed_indegree) /indegree CE /autapses false /multapses true >> << /model /syn_ex >> Connect
M_INFO (BuildNetwork)
(Connecting inhibitory -> inhibitory population.) message
I_neurons I_neurons << /rule (fixed_indegree) /indegree CI /autapses false /multapses true >> << /model /syn_in >> Connect
record_spikes true eq
{
E_net GetLocalNodes size /local_neurons Set
NE Nrec_total lt
{
M_ERROR (BuildNetwork)
(More neurons requested recorded than available. Aborting the simulation!) message
1 quit_i
} if
M_INFO (BuildNetwork)
(Connecting spike detectors.) message
E_from dup Nrec_total add 1 sub cvgidcollection E_detector dup cvgidcollection << /rule (all_to_all) >> << /model /syn_std >> Connect
} if
% read out time used for building
toc /BuildEdgeTime Set
mpi_rank cvs (] ) join (BuildEdge time : ) join =only BuildEdgeTime =only ( s) =
BuildEdgeTime cvs ( # build_edge_time ) join logger /log call
memory_thisjob cvs ( # virt_mem_after_edges) join logger /log call
} def % end of buildnetwork
/RunSimulation
{
% open log file
log_file logger /init call
ResetKernel
statusdict /MPI_Rank known
{
/mpi_rank statusdict /MPI_Rank get def
}
{
/mpi_rank 0 def % serial version
}
ifelse
%stdp_params using
memory_thisjob cvs ( # virt_mem_0) join logger /log call
BuildNetwork
tic
presimtime Simulate
toc /PreparationTime Set
mpi_rank cvs (] ) join (Preparation time : ) join =only PreparationTime =only ( s\n) =
memory_thisjob cvs ( # virt_mem_after_presim) join logger /log call
PreparationTime cvs ( # presim_time) join logger /log call
tic
simtime Simulate
toc /SimCPUTime Set
memory_thisjob cvs ( # virt_mem_after_sim) join logger /log call
SimCPUTime cvs ( # sim_time) join logger /log call
%(\nBRN STDP Equilibrium Simulation) =
mpi_rank cvs (] ) join (Simulation time : ) join =only SimCPUTime =only ( s\n) =
logger /done call
} def
/*
This script defines a logger class used to properly log memory and timing
information from network simulations. It is used by hpc_benchmark.sli to
store the information to the log files.
*/
/logger <<
/max_rank_cout 5 % copy output to cout for ranks 0..max_rank_cout-1
/max_rank_log 30 % write to log files for ranks 0..max_rank_log-1
/line_counter 0
% constructor
% expects file name on stack
/init
{
Rank max_rank_log lt {
(_) join
% convert rank to string, prepend 0 if necessary to make
% numbers equally wide for all ranks
Rank cvs
dup length max_rank_log cvs length exch sub
{
48 prepend % 48 is ASCII code for 0
}
repeat
join
(.dat) join
/f exch (w) ofsopen
{
def
}
{
/Logger_init /CannotOpenFile raiseerror
}
ifelse
} if
}
% logging function
% expects one operand on stack to write to file
/log
{
/value Set
Rank max_rank_log lt {
f line_counter <- ( ) <- Rank <- ( ) <- value <- (\n) <- pop
/line_counter line_counter 1 add def
} if
Rank max_rank_cout lt {
cout Rank <- ( ) <- value <- endl flush pop
cerr Rank <- ( ) <- value <- endl flush pop
} if
}
% closes file
/done
{
Rank max_rank_log lt {
f close
} if
}
>> def
% ------------------------------------------------------------------------------------
RunSimulation
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment