# Julia – excellent opensource alternative for MATLAB

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:

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:

1. To create logical lattice: “false(N,N)” to “falses(N,N)”.
2. Array elements are accessed with “[i]” instead of “(i)”.
3. To create array of all permutations: “perms(uint8(1:8))'” to “hcat(collect(permutations(1:8))…)”.
4. “bsxfun()” to “broadcast()”
5. Adding elements to vectors using “push!()” function instead of concatenation.
6. 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[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…

## 2 thoughts on “Julia – excellent opensource alternative for MATLAB”

1. […] Julia (post with the code here) […]

Like

2. Marc says:

Hey, here’s some suggestion for some speedups for the Julia CA: https://github.com/marcjwilliams1/julia-notebooks/blob/master/cancer%20stem%20cell%20CA.ipynb

Julia doesn’t like global variables so wrapping everything inside a function usually helps. Also you should time on the second time you call a function, first time round the just in time compiler takes some time to make some inferences on the types to produce fast machine code.

Noticed some modification for the returnemtpyPlace function, included these in the loopCA2 and returnEmptyPlace2 functions. Be interested in how this now compares to everything else. On my machine the loop version is about 5X faster than the vectorized version for 10^5 cells.

Like