After my last post about ABM speed comparison @marcjwills_ suggested that Julia programming language could be a good alternative for slower, but easy to code programming languages – I’ve decided to give it a try. On the project webpage I found promising speed comparisons showing that Julia can indeed be my new tool of choice. After easy download and installation I saw the terminal window that brought a huge smile to my face.

Quick web search allowed me to find nice IDE for Julia called Juno – simple and not overloaded with features. What might be a great news to some people (especially @dbasanta) is that Julia after few easy tweaks can use matplotlib.

The obvious first test to do is to update speed comparison post with Julia. First, I’ve decided to rewrite the partially vectorized MATLAB code in Julia without making any essential changes. After browsing the web I found couple of useful webpages that helped me with doing that without much effort:

- http://julia.readthedocs.org
- http://julia.readthedocs.org/en/latest/manual/noteworthy-differences/
- http://learnxinyminutes.com/docs/julia/

The differences between Julia and MATLAB are really minor and most of the functions work in exactly the same way. I only had to be careful with logical operations between matrices and scalars – in Julia you need to add “.” to the operator to make it work as in MATLAB. Other changes that I had to make were:

- To create logical lattice: “false(N,N)” to “falses(N,N)”.
- Array elements are accessed with “[i]” instead of “(i)”.
- To create array of all permutations: “perms(uint8(1:8))'” to “hcat(collect(permutations(1:8))…)”.
- “bsxfun()” to “broadcast()”
- Adding elements to vectors using “push!()” function instead of concatenation.
- Deleting elements by using “deleteat!()” function.

The working code in Julia is the following:

N = 1000; #square domain dimension pprol = 1/24; #probability of proliferating pmig = 15/24; #probability of migrating pdeath = 5/100; #probability of dying pmax = int8(10); #proliferation capacity ps = 0.3; #probability of symmetric division L = falses(N,N); L[[1:N, 1:N:N*N, N:N:N*N, N*(N-1):N*N]] = true; #boundary L[N*round(N/2)+round(N/2)] = true; cells = int32([N*round(N/2)+round(N/2)]); cellsIsStem = [true]; cellsPmax = int8([pmax]); aux = int32([-N-1 -N -N+1 -1 1 N-1 N N+1])'; #indices to heighborhood Pms = hcat(collect(permutations(1:8))...); #permutations nP = size(Pms,2); #number of permutations nSteps = 0; while(length(cells)<10^4) sh = randperm(length(cells)); cells = cells[sh]; cellsIsStem = cellsIsStem[sh]; cellsPmax = cellsPmax[sh]; S = broadcast(+,cells',aux[Pms[:,rand(1:nP,length(cells))]]) S[L[S]] = 0; #setting occupied spots to 0 indxF = int32(find(any(S,1))); #selecting cells with at least one free spot nC = length(indxF); #number of cells with free spot P = rand(nC).<pprol; #proliferation Ps = P & (rand(nC).<ps) & (cellsIsStem[indxF]);#symmetric division De = P & (cellsPmax[indxF] .== 0);#proliferation exhaution D = P & (rand(nC).<pdeath) & (~cellsIsStem[indxF]); #death at proliferation attempt M = (~P) & (rand(nC).<pmig); #go when no grow del = D | De; #cells to delete act = find((P | M) & (~del)); #indices to the cells that will do something for ii = act #only for those that will do anything ngh = S[:,indxF[ii]]; ngh = ngh[ngh.>0]; indO = find(~L[ngh]); #selecting free spot if ~isempty(indO) #if there is still a free spot indO = indO[1]; L[ngh[indO]] = true; if P[ii] #proliferation push!(cells,int32(ngh[indO])); if Ps[ii] #symmetric division push!(cellsIsStem,true); push!(cellsPmax,cellsPmax[indxF[ii]]); else push!(cellsIsStem,false); push!(cellsPmax, cellsPmax[indxF[ii]]-1); if ~cellsIsStem[indxF[ii]] cellsPmax[indxF[ii]] = cellsPmax[indxF[ii]]-1; end end else #migration L[cells[indxF[ii]]] = false; cells[indxF[ii]] = int32(ngh[indO]); end end end if any(del) #updating death L[cells[indxF[del]]] = false; deleteat!(cells,indxF[del]); deleteat!(cellsIsStem,indxF[del]); deleteat!(cellsPmax,indxF[del]); end end

To visualize the results of simulations we can use Gadfly plotting package and use simple spy() function (analog to MATLAB).

using Gadfly spy(L)

Ok, so modifications were easy and straightforward, but what about speed? Amazingly Julia performed almost exactly the same as MATLAB (see plot below).

But that is not all… MATLAB code was using vectorization feature. Without that, when using only loops like in C++, it would be several times slower. Is that also the case in Julia?

Below is the C++ style code in Julia – only loops, no vectorization. One of the awesome feature of Julia, that MATLAB lacks, is definition of functions within the script (see “returnEmptySpace” function below).

N = 1000; #square domain dimension aux = int32([-N-1 -N -N+1 -1 1 N-1 N N+1]); #indices to heighborhood pprol = 1/24; #probability of proliferating pmig = 15/24; #probability of migrating pdeath = 5/100; #probability of dying pmax = int8(10); #proliferation capacity ps = 0.3; L = falses(N,N); L[[1:N, 1:N:(N*N), N:N:(N*N), (N*(N-1)):(N*N)]] = true; #boundary L[N*round(N/2)+round(N/2)] = true; cells = [int32(N*round(N/2)+round(N/2))]; cellsIsStem = [true]; cellsPmax = [int8(pmax)]; #defining function that returns random empty place around cell function returnEmptyPlace(indx::Int32) neigh = int32(zeros(1,8)) nF = 0 for i = 1:8 if ~L[indx+aux[i]] neigh[nF+1] = indx+aux[i] nF+=1 end end if nF>0 return neigh[rand(1:nF,1)]; else return 0; end end while length(cells)<10^5 sh = randperm(length(cells)); cells = cells[sh]; cellsIsStem = cellsIsStem[sh]; cellsPmax = cellsPmax[sh]; newCells = Int32[] newCellsIsStem = Bool[] newCellsPmax = Int8[] i = 1; while i<=length(cells) deleted = false; newSite = returnEmptyPlace(cells[i])[1]; if newSite > 0 #if there is a new spot if rand()<pprol if cellsIsStem[i] #is stem cell push!(newCells, newSite); L[newSite] = true; push!(newCellsPmax,cellsPmax[i]); if (rand()<ps) #symmetric division push!(newCellsIsStem, true); else push!(newCellsIsStem, false); end else if (cellsPmax[i]>0) && rand()>pdeath cellsPmax[i]-=1 L[newSite]=true; push!(newCells,newSite); push!(newCellsPmax,cellsPmax[i]); push!(newCellsIsStem, false); else L[cells[i]]=false; splice!(cells,i) splice!(cellsPmax,i) splice!(cellsIsStem,i) deleted = true; end end elseif rand() < pmig L[newSite]=true; L[cells[i]]=false; push!(newCells,newSite); push!(newCellsPmax,cellsPmax[i]); push!(newCellsIsStem, cellsIsStem[i]); splice!(cells,i) splice!(cellsPmax,i) splice!(cellsIsStem,i) deleted = true; end end if ~deleted i+=1; end end cells = [cells, newCells]; cellsPmax = [cellsPmax, newCellsPmax]; cellsIsStem = [cellsIsStem, newCellsIsStem] end

How about the speed of a new code? I was amazed with the results…