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…
COMPUTE CANCER 