Loop-free search for proliferative cells in MATLAB

When thinking about implementation of agent-based model only few of us consider MATLAB as a coding environment, even though it has plenty of built-in visualization and analysis tools. I get that. It’s a high level language that is not very efficient with for loops, so a basic loop-based code will execute much faster when coded in C++ (and others). However, there is a very cool feature of MATLAB – vectorization, that allows to skip plenty of loops and significantly speed-up the computations. In this post I’ll show how to utilize this feature to make a loop free search for a cells that have a free spot in the neighborhood.

Our computational domain is represented by N x N Boolean matrix L (a true value indicates the lattice point is occupied). An additional integer vector cells is maintained to store indices for all viable cells present in the system (to avoid costly search for cells on the lattice). Our task is to write a code that returns indices to those cells that have at least one free spot in the 8 spot neighborhood.

A basic loop-based code in MATLAB that will do the task is the following:

cellsInd = false(size(cells));
    for j = 1:length(cells)
        loop = true; 
        for ii = -1:1:1
            if ~loop
                break;
            end
            for kk = -N:N:N
                if ~L(cells(j)+ii+kk)
                   loop = false; 
                   cellsInd(j) = true;
                   break;
                end
            end
        end
    end
    cellsInd = find(cellsInd);
end

The code is similar to the one that would be written in C++. For 2000×2000 lattice with 800,000 cells evaluation of 200 iterations of the above code take about 83 seconds. Not very impressive… Same task for a code written in C++ takes about 3 seconds (when using vectors from STL library). Can we do something about that?

As I wrote, vectorization is a very cool MATLAB feature. When using it, coding the whole task takes only one two lines:

aux = [-N-1; -N; -N+1; -1; 1; N-1; N; N+1]';
cellsInd = find(~all(L(bsxfun(@plus,cells,aux))));

Moreover, the evaluation time is reduced from about 83 seconds to about 9.5 seconds. Impressive speed-up and elegant code in two lines.

Still, both codes are outperformed by C++, but the difference when using vectorized code is bearable.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s