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