Groups 141 of 99+ julia-users › PriorityQueue for fast large network simulation 8 posts by 3 authors Rainer J. Engelken Jun 30 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)") Jeffrey Sarnoff Jul 12 Is it possible to set up a priority queue with integer keys instead? Yes. julia> pq=Collections.PriorityQueue() Base.Collections.PriorityQueue{Any,Any,Base.Order.ForwardOrdering} with 0 entries julia> pq[1]=10;pq[2]=5;pq[3]=20; julia> pq Base.Collections.PriorityQueue{Any,Any,Base.Order.ForwardOrdering} with 3 entries: 2 => 5 3 => 20 1 => 10 julia> pq[2] += 1; pq Base.Collections.PriorityQueue{Any,Any,Base.Order.ForwardOrdering} with 3 entries: 2 => 6 3 => 20 1 => 10 More generally, is there a smarter data structure than priority queue? Probably. 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? How many dimensions has your 'array' (is it a vector, a matrix ..)? On each iteration, do you want to find the value of smallest element in the array or the location(s) of the smallest element or both? On each iteration, are the entries that are to be altered sequential or at random locations, how are those locations to be determined? - show quoted text - David P. Sanders Jul 12 For efficiency, you probably want a typed version: julia> Base.Collections.PriorityQueue(Int, Int) Base.Collections.PriorityQueue{Int64,Int64,Base.Order.ForwardOrdering} with 0 entries - show quoted text - Rainer J. Engelken Jul 12 Re: [julia-users] Re: PriorityQueue for fast large network simulation Is it possible to set up a priority queue with integer keys instead? Yes. julia> pq=Collections.PriorityQueue() Base.Collections.PriorityQueue{Any,Any,Base.Order.ForwardOrdering} with 0 entries julia> pq[1]=10;pq[2]=5;pq[3]=20; Sorry, I wasn't clear. What I meant to say: Profiling indicates that the code I posted spends most of the time with dictionary operations during changing the priority of an existing key, e.g. finding the index where a key is stored, and checking whether the slot is empty or filled and changing (functions {ht_keyindex, ht_keyindex2} in {setindex! and getindex} in dict.jl called by setindex! and getindex in Collections.jl) . I was just wondering, whether this can be avoided or sped up, e.g. by using arrays instead of dictionaries. More generally, is there a smarter data structure than priority queue? Probably. 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? How many dimensions has your 'array' (is it a vector, a matrix ..)? It is the vector phi in the code I posted: phi = Collections.PriorityQueue(UInt32,Float64) @inbounds for i = 1:n #set initial phases for all neurons phi[i] = -rand() end On each iteration, do you want to find the value of smallest element in the array or the location(s) of the smallest element or both? Both. On each iteration, are the entries that are to be altered sequential or at random locations, how are those locations to be determined? They are at fairly random locations, determined by the network topology. Their indices are given by postidx, which is a vector of k random integers uniformly distributed in 1:n. In my case k=10^2, n=10^8. Looking forward to suggestions, Regards Rainer - show quoted text - Rainer J. Engelken Jul 12 Re: [julia-users] Re: PriorityQueue for fast large network simulation For efficiency, you probably want a typed version: julia> Base.Collections.PriorityQueue(Int, Int) Base.Collections.PriorityQueue{Int64,Int64,Base.Order.ForwardOrdering} with 0 entries I used Uint32 as keys to save memory. (phi = Collections.PriorityQueue(UInt32,Float64) ). Are Int64 keys faster? Best Rainer - show quoted text - Jeffrey Sarnoff Jul 12 Re: [julia-users] Re: PriorityQueue for fast large network simulation Broadly speaking, UInt32 is faster than UInt64, Int32 is faster than Int64 and Julia tends to use the signed ints when it means integer values. Your application may be happier with the use of ako (possibly augmented) tree instead of a priority queue. You get very fast index or weight or priority ordered retrievals and updates tend to amortize away with the number of values you are keeping -- NearestNeighbors.jl looks worth while. - show quoted text - Jeffrey Sarnoff Jul 12 Re: [julia-users] Re: PriorityQueue for fast large network simulation Keno's AbstractTrees is very helpful in understanding tree-think. - show quoted text - Rainer J. Engelken Jul 12 Re: [julia-users] Re: PriorityQueue for fast large network simulation 2016-07-12 19:17 GMT+02:00 Jeffrey Sarnoff : Broadly speaking, UInt32 is faster than UInt64, Int32 is faster than Int64 and Julia tends to use the signed ints when it means integer values. Thanks, I didn't know that. Is there a way to make the dictionary underlying the PriorityQueue also use Uint32 vals? By default, it seems to use Int64: test = Collections.PriorityQueue(UInt32(1):UInt32(4),rand(4)) Dict{UInt32,Int64} with 4 entries: 0x00000004 => 1 0x00000002 => 4 0x00000003 => 3 0x00000001 => 2 Your application may be happier with the use of ako (possibly augmented) tree instead of a priority queue. You get very fast index or weight or priority ordered retrievals and updates tend to amortize away with the number of values you are keeping -- NearestNeighbors.jl looks worth while. Isn't a binary heap underlying the PriorityQueue? Or do you have a specific tree structure in mind that allows fast changing of a small fraction of entries (say 10^2 out of 10^8) per peek. @Tim's suggestion of doing batch update sounds appealing, but there might also be other data structures apt to the problem. Another question concerning priority queues: Has somebody tried to parallelize a priority queue in julia? In principle it should be possible, e.g.: Brodal, Gerth Stoelting. "Priority queues on parallel machines." Parallel Computing 25.8 (1999): 987-1011. http://www.sciencedirect.com/science/article/pii/S0167819199000320 Best Rainer - show quoted text -