I am trying to figure out if I have indeed written this small test bit of code in the most efficient way for Julia. The overall performance I am getting is better than what I had with python using numpy arrays. However, I am not sure if I am missing some important functions (or parallelization opportunities).
The basic premise of this code is that there are some particles in a 1D box that diffuse around and reversibly stick to some region.
I ran the following code several times and the best run I got was:
@time f =stickyparticle()
4.960075 seconds (1.95 M allocations: 4.613 GiB, 7.55% gc time)
It seems to me like that last two numbers are fairly high. Is it largely because of the random numbers being generated on every loop?
using Distributions function stickyparticle() nm = 5000;#number of particles D = 4*10^6; # diffusiion coefficent h = 10000; # length of container totaltime = 300; dt = 0.01; #time step size fpart = rand(nm)*h; #initialize particle positions kon = 5; #Binding Constant koff = 0.5; #Unbinding Constant iters = collect(dt:dt:totaltime); #number of iterations for simulation s=1000; #Range of sticking zone: x = [0:s] free = fpart.>0; #Boolean array tracking free particles bound = .~(free); #Boolean array tracking stuck particles drnd = Normal(0,sqrt(2*D*dt)); #Standard deviation for particle diffusion local wpos::BitArray local wneg::BitArray local unbind::BitArray local bind::BitArray @fastmath @inbounds @views for t in iters freepart = view(fpart,free); #Get only free particles din = rand(drnd,length(freepart)); #Calculate diffusion of particles freepart .= (din+freepart); #Update particle positions wpos = freepart .> h ; #Find particles that diffuse beyond the length of the container freepart[wpos] .-= 2*h; #Reflect particles back wneg = freepart .< 0; #Find particles that diffuse beyond second boundary freepart[wneg] .= -freepart[wneg]; #Reflect particles back unbind = (rand(length(bound)) .< (dt*koff)) .& bound #Calculate stuck particles that unbind bind = ((rand(length(free)) .< (dt*kon) ) .& (fpart .<=s)).& free; #Calculate free particles that stick #Stick and unstick particles bound[unbind] .= false; free[unbind] .= true; bound[bind] .= true; free[bind] .= false; end return fpart end