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.

Screen Shot 2015-10-04 at 3.00.41 PM

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 = 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)

Screen Shot 2015-10-04 at 3.41.52 PM

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

speed

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…

speed2

Advertisements