Binary Indexing, For The Simple Minded

While writing almost any algorithm, one rarely doesn’t encounters the most dreadful of tasks, indexing. Although it is comprised of simple multiplications and additions, placing elements in the correct order in a matrix is perplexing even for the most experienced of programmers. Not only that, but when a large number of elements need to be placed at once, this also becomes time-expensive, due to arbitrary memory access patterns.

To start this post off, let’s discuss some access patterns you might run into. Let’s take for example a 2D Finite-Element (FEM) problem, solved with triangular elements. In most cases, the FEM problem is linearized and resolved by solving the equation

\mathbf{K}\mathbf{x} = \mathbf{b}

where K is the matrix connecting between the node values (x) and the sources, b. A lot of fancy words which means that this linear equation set has to be solved for x. Let’s look at some pretty pictures then, instead! Every triangle has three nodes (naturally), however some nodes in the mesh are common two two or more triangles. For example, for the following mesh snippet

The element\node lookup table can be

Element #Node 1Node 2Node 3
1126
2246
3234
4564

Note that a clock-wise rotation direction here was maintained while choosing the node order. The node coefficients are accumulated into K as such

\{K\}_{ij} = \sum\limits_{e}{k_{ij}^{(e)}}

namely, each i-th node which corresponds with the j-th node of the same element is summed along with all the coefficients of the elements (namely, triangles) that share the same node.

Are you confused yet? Good! This is only the indexing part of the problem, don’t forget. However, implementing this is actually easier than the mathematical phrasing. Each element has a 3x3 matrix of coefficients which goes along with it,

\mathbf{K}^{(e)} = \left[{\begin{matrix}k_{11}^{(e)} & k_{12}^{(e)} & k_{13}^{(e)} \\ k_{21}^{(e)} & k_{22}^{(e)} & k_{23}^{(e)} \\ k_{31}^{(e)} & k_{32}^{(e)} & k_{33}^{(e)}\end{matrix}}\right]

So all that is left is linking the “local” 1;2;3 indices to their global ij counterparts. One way of doing so is with a couple of for loops.

% For this script we assume there is a lookup table, 'luTable' 
% in which each row denotes the element number and each column
% the three nodes assigned to it. The dimensions of 'luTable'
% are [Me,3].
%
% Also the per-element k matrices are given by a [3,3,Me] 
% matrix array, called 'Ke'.

for e = 1:size(luTable,1)
cNodes = luTable(e,:);
    for i = 1:3
        for j = 1:3
	    c_i = cNodes(i);
            c_j = cNodes(j);
            K(c_i,c_j) = K(c_i,c_j) + Ke(i,j,e);
        end
    end
end

Disclaimer: From this point (in the post) forward, I will stop showing off my “knowledge” and start teaching indexing tricks.

Not very elegant, but certainly gives a better perspective of all the explanations above. Let this soak in a bit, before delving into more elegant ways to write this down. The “Good” news are, however, that memory access wise, there is almost no way to make this more efficient. So for this port of the post, the discussion will focus on elegant writing. Luckily, the MATLAB\Octave interpreter allows snipping matrices in the following fashion: Accessing the submatrix

submat_a = [a(3,5) a(3,7) a(3,12) ; a(6,5) a(6,7) a(6,12) ; a(10,5) a(10,7) a(10,12)];

is equivalent to snipping it by

submat_a = a([3 6 10],[5 7 12])

Hence, the previous piece of code can be re-written as

for e = 1:size(luTable,1)
    cNodes = luTable(e,:);
    K(cNodes,cNodes) = K(cNodes,cNodes) + Ke(:,:,e);
end
[/cc]

Well, that’s one sort of endeavor you might go through. Let’s consider other cases! For example, if one would want to randomize several numbers from a list, with no return. In case of randomization from a single set of numbers, this could be settled with a solution such as this:

clc
clear

%% 

listLength = 10;
nChoices = 3;

A = floor(rand([listLength 1])*100);

A_snip = [];

for i = 1:nChoices
  % Current choice index from A
  cChoice = floor(rand(1)*size(A,1))+1;
  % Concatenate to the chosen numbers from A
  A_snip = [A_snip ; A(cChoice)];
  % And remove the number from A
  A(cChoice) = [];
end

However for now I will write it in an iterative version. Now, let’s say, you would want to randomize from several sets of numbers (of the same length, of course), Again, with no return. There are two corresponding problems. One is that in every iteration some numbers are removed from each set. The Second lies in the fact that the numbers chosen in each iteration are not necessarily in the same row\column. Lets start by choosing numbers from each set, and removing them from the list: The function idx2 is simply it is definitely a nifty function to carry around, in my opinion.

clc
clear

%% 

listLength = 10;
nSets = 7;
nChoices = 3;

A = floor(rand([listLength nSets])*100);

A_snip = [];

for i = 1:nChoices
  % Current choice index from A
  cChoice = floor(rand([1 nSets])*size(A,1))+1;
  % Concatenate to the chosen numbers from A
  cSnip = A(idx2(cChoice,1:size(A,2),size(A,1)));
  A_snip = [A_snip ; cSnip];
  % And remove the number from A
  A(idx2(cChoice,1:size(A,2),size(A,1))) = [];
end

where

function sub = idx2(m,n,M)
sub = m + (n-1)*M;
end

However, this program won’t run smoothly yet. The matrix A is deformed after clearing it from the chosen elements. Hence, it remains is to reshape the array to have a number of rows and columns which is suitable to a list set shorter in one element, each.

clc
clear

%% 

listLength = 10;
nSets = 7;
nChoices = 3;

A = floor(rand([listLength nSets])*100);

A_snip = [];

for i = 1:nChoices
  % Store original dimensions of A
  origDims = size(A);
  
  % Current choice index from A
  cChoice = floor(rand([1 nSets])*size(A,1))+1;
  % Concatenate to the chosen numbers from A
  cSnip = A(idx2(cChoice,1:size(A,2),size(A,1)));
  A_snip = [A_snip ; cSnip];
  % And remove the number from A
  A(idx2(cChoice,1:size(A,2),size(A,1))) = [];
  
  % Reshape the original list to be one row less in size
  A = reshape(A,origDims-[1 0]);
end

So indexing in a 2D (or more) array isn’t that bad, is it? Well, it can get worse, but you definitely get used to it after a while. So when is this type of indexing not good? Mainly, when the elements you are trying to access have a different pattern than just the index order. For example, if it is required to color parts of a map in different colors. In that case, logical indexing is suggested. In fact, it is so suggested, that if you will try to run the following code in MATLAB, it will suggest on its own to use logical indexing in lines 32-34:

clc
clear
close all

%%

Xrange = [-1 1];
Yrange = [-1 1];

nPts = [200 200];

circRad = 0.4;
circCenter = [-0.5 -0.2];

rectStart = [0.2 0.3];
rectStop = [0.7 0.8];

% Decide how to mark the filled shapes.
circleIndex = 1;
rectIndex = 2;

%%

% create mesh
[X,Y] = meshgrid(linspace(Xrange(1),Xrange(2),nPts(1)),...
                 linspace(Yrange(1),Yrange(2),nPts(1)));
                 
% Create matrix to fill
A = zeros(size(X));

% Fill circle area
A(find((X-circCenter(1)).^2+(Y-circCenter(2)).^2 <= circRad^2)) = circleIndex;
A(find((X >= rectStart(1)) & (X <= rectStop(1)) & ...
       (Y >= rectStart(2)) & (Y <= rectStop(2)))) = rectIndex;

       
hdl = imagesc(X(1,:),Y(:,1),A);
axis('square');
colorbar;
xlabel('x','fontsize',18);
ylabel('y','fontsize',18);
set(

So basically the only change you need is to ditch the find command:

A((X-circCenter(1)).^2+(Y-circCenter(2)).^2 <= circRad^2) = circleIndex;
A((X >= rectStart(1)) & (X <= rectStop(1)) & ...
  (Y >= rectStart(2)) & (Y <= rectStop(2))) = rectIndex; 

Such an indexing pattern is recommended for most of the problems that involve a structured mesh, but not restricted to them at all. As a general coding pattern, usually I mark these binary maps in separate arrays with a ‘BinMap’ post-fix. For example:

circBinMap = (X-circCenter(1)).^2+(Y-circCenter(2)).^2 <= circRad^2
rectBinMap = (X >= rectStart(1)) & (X <= rectStop(1)) & ...
             (Y >= rectStart(2)) & (Y <= rectStop(2));

A(circBinMap) = circleIndex;
A(rectBinMap) = rectIndex; 

That’s it for now! I have to say that these simple tricks, presented here and in the last two articles, compose about 80% of the clean coding habits I assigned to myself over the years. The rest is just discipline to constantly document the code and mostly trying to write shorter code, altogether.

Feel free to contact me with any questions and ask for additions.

Cheers!

Leave a Reply

Your email address will not be published. Required fields are marked *