MATLAB: Vectorization of 2 for loops in MATLAB

MATLABvectorization

I have just started exploring world of vectorization. I got the 1-D vectorization down but i am having trouble vectorizing the following code.
CityPairs = [7 3
3 1
3 1
1 7
7 1
3 4
5 1
4 6];
Offices = [1;3;7];
nOffices = size(Offices,1);
connection = zeros(nOffices);
for i = 1:nOffices
    for j = 1:nOffices
        connection(i,j) = sum(Offices(i) == CityPairs(:,1)...
            & CityPairs(:,2) == Offices(j));
    end
    connection(i,:) = connection(i,:);
end   
disp(connection)
In this example there are 7 cities, three of which have offices. I want a pairwise matrix for the cities with offices to capture the sum of all the one way connections between each. The answer for above problem should be:
   0     0     1
   2     0     0
   1     1     0
Any suggestions are welcome. Thanks in advance.
Keith

Best Answer

It is not trivial to "vectorize" this setup, as there are operations that require some look up table, and there is accumulation (unless you find a trick). This is likely to make the approach without FOR loops more complicated (code-wise) than the basic, loop-based approach. So let's start with the trick first ;-) ..
 >> A = double( bsxfun(@eq, CityPairs(:,1), Offices.') ) ;
 >> B = double( bsxfun(@eq, CityPairs(:,2), Offices.') ) ;
 >> A.' * B
 ans =
     0     0     1
     2     0     0
     1     1     0
I let you study a bit BSXFUN and understand how this solution works. If you don't understand after having spent some time, ask me for details.
Now if I hadn't found a trick, I would have done it as follows (I display the output at each step to make it more understandable)..
First, we find valid city IDs (those which have an office):
 >> isValid = ismember(CityPairs, Offices)
 isValid =
     1     1
     1     1
     1     1
     1     1
     1     1
     1     0
     0     1
     0     0
Then we extract rows of CityPairs which have two valid nodes/cities:
 >> validPairs = CityPairs(all(isValid, 2),:)
 validPairs =
     7     3
     3     1
     3     1
     1     7
     7     1
We want to build a list of office IDs that correspond to valid pairs of cities, and for that we build a look up table:
 >> LT(Offices) = 1 : nOffices
 LT =
     1     0     2     0     0     0     3
You can see that LT(7) is 3, which is the office ID of city 7. With that, we can build the list of office IDs which correspond to valid city pairs:
 >> officeIds = LT(validPairs)
 officeIds =
     3     2
     2     1
     2     1
     1     3
     3     1
The last step is to accumulate a count of 1 for each row of the above array, which actually counts entries with same indices:
 >> accumarray(officeIds, ones(size(officeIds,1),1), [nOffices, nOffices])
 ans =
     0     0     1
     2     0     0
     1     1     0
This last step should be done using SPARSE instead of ACCUMARRAY in cases where there is a large number of offices. The trick with BSXFUN could also benefit from converting one of its inputs to sparse, e.g. the vector of offices.
PS: if it was for a homework, I doubt that either method is what was required.
Related Question