Hi,
I am trying to implement a fast event-based numerically exact
simulation of a sparse large spiking neural network using a priority
queue. It is fast, but not fast enough.
Profiling indicates that the bottleneck seem to be the dictionary
operations keyindex and setindex! when changing priority of an existing key
(3rd line of
function ptc! in code below). Is there a way to avoid this? Is is
possible to set up a priority queue with integer keys instead?
More generally, is there a smarter data structure than priority queue,
if I want to find the smallest element of an large array and change a small
fraction array entries (say 10^2 out of 10^8) each iteration?
Best regards
Rainer
## For those interested into the implementation: We map the membrane
potential to a phase variable phi \in [-1, inf) and solve the network
evolution of an purely inhibitory homogeneous pulse-coupled leaky integrate
and fire
network analytically between network spikes. This is just a toy example,
some important steps are missing.
##### code ###
function ptc!(phi, postid, phishift, w, c) #phase transition curve
@inbounds for i = postid
phi[i] = w*log(exp((phi[i]-phishift)/w)+c) + phishift
end
end
function lifnet(n,nstep,k,j0,ratewnt,tau,seedic,seedtopo)
iext = tau*sqrt(k)*j0*ratewnt # iext given by balance equation
w = 1/log(1. + 1/iext) # phase velocity LIF
c = j0/sqrt(k)/(1. + iext) # effective coupling in PTC for LIF
phith, phireset = 1., 0. # reset and threshold for LIF
srand(seedic) # initialize random number generator (Mersenne Twister)
phi = Collections.PriorityQueue(UInt32,Float64)
@inbounds for i = 1:n #set initial phases for all neurons
phi[i] = -rand()
end
spikeidx, spiketimes = Int64[], Float64[] # initialize time & spike
raster
postidx = Array{UInt32,1}(k) # idices of postsynaptic neurons
t = phishift = 0. # initial time and initial global phase shift
for s = 1 : nstep #
j, phimax = Collections.peek(phi) # next spiking neuron
dphi = phith+phimax-phishift # calculate next spike time
phishift += dphi # global backshift instead of shifting all phases
t += dphi # update time
state = UInt(j+seedtopo) # spiking neuron index is seed of rng
@inbounds for i = 1:k # idea by R.Rosenbaum to reduce memory
postidx[i], state = randfast(n,state) # get postsyn. neurons
end
ptc!(phi,postidx,phishift,w,c) # evaluate PTC
phi[j] = phireset + phishift # reset spiking neuron to 0
push!(spiketimes,t) # store spiketimes
push!(spikeidx,j) # store spiking neuron
end
nstep/t/n/tau*w, spikeidx, spiketimes*tau/w
end
function randfast(m, state::UInt) # linear congruential generator
state = lcg(state) # faster then Mersenne-Twister
1 + rem(state, m) % Int, state # but can be problematic
end
function lcg(oldstate) # tinyurl.com/fastrand by M. Schauer
a = unsigned(2862933555777941757)
b = unsigned(3037000493)
a*oldstate + b
end
#n: # of neurons, k: synapses/neuron, j0: syn. strength, tau: membr. time
const.
n,nstep,k,j0,ratewnt,tau,seedic,seedtopo = 10^6,10^4,10^2,1,1,.01,1,1
@time lifnet(100, 1, 10, j0, ratewnt, tau, seedic, seedtopo); #warmup
gc()
@time rate,spidx,sptimes = lifnet(n,nstep,k,j0,ratewnt,tau,seedic,seedtopo);
#using PyPlot;plot(sptimes,spidx,".k");ylabel("Neuron Index");xlabel("Time
(s)")
