Vectorize a for loop and if-else statement

Hello. Can anybody help me to vectorize:

L=1;
A=0;
for L=1:100
 for i=1:25
  A=A+1;
  if E(L).a == N(i).a
     HD(A,L)=myfunction(T0,M0,T(i).a,M(i).a);
  else
     HD(A,L)=0;
  end
  A=A+1;
  if E(L).a == N(i).b
     HD(A,L)=myfunction(T0,M0,T(i).b,M(i).b);
  else
     HD(A,L)=0;
  end

end
end

where E,N,T,M are struct. T0 and M0 are 128x720 double. Thanks

[Stephen23]

Vectorization does not mean replacing my slow loops with some other code, it means writing functions that operate on entire arrays at once.

This means it is the function that is important, not the loops. However you do not tell us anything about myfunction, yet it is the only important thing to consider. If you give us details about the function then we can give advice about vectorization.

[shahrizan jamaludin]

myfunction will return a Hamming distance value (eg. 0.012) every time I call it. And I want to put this value into array HD(A,L).

[Ced]

What Stephen meant was how your function computes its output, in this case the hamming distance. If your function does this in vectorized form, it will be more efficient to give it a full vector instead of looping. But this is not always possible.

Example:

% 1: Scalar case
function y = norm_scalar(x)
  y = sqrt(x*x);
end

% 2: vector case, but looped
function y = norm_notvectorized(x)
  N = length(x); 
  y = 0;
  for i = 1:N
    y = y + x(i)*x(i);
  end
  y = sqrt(y);
end

% 3: vector case, vectorized  
function y = norm_vectorized(x)
  y = sqrt(x'*x);
end

The scalar case doesn't make much sense here, but couldn't think of anything meaningful right now, and it's just to show the difference:

In case of function 1, you actually cannot pass it a vector. In function 2, you can pass it a vector, but it will perform the operation very inefficiently. function 3 is vectorized, and will be much more efficient than function 2.

In short: The important part is myfunction, not the for loop.

[shahrizan jamaludin]

Thanks for your suggestion. I just wondering, if I can do something to change HD(A,L)=statement to become HD=statement only.

[Stephen23]

The answer depends on what happens inside myfunction. You need to understand: vectorization means that myfunction needs to be written using vectorized code. It is not the loop or the allocation that is important: it is the code inside myfunction that we need to check.

When you show us myfunction then we can tell you if it is, or can be vectorized.

[shahrizan jamaludin]

This is myfunction

function hd = myfunction(template1, mask1, template2, mask2)
[ro,co]=size(template1);
num=co/2;
scales=1;
template1 = logical(template1);
mask1 = logical(mask1);
template2 = logical(template2);
mask2 = logical(mask2);
hd = NaN;
for shifts=-8:8
  template1s = shiftbits(template1, shifts,scales);
  mask1s = shiftbits(mask1, shifts,scales);
  mask = mask1s | mask2;
  nummaskbits = sum(sum(mask == 1));
  totalbits = (size(template1s,1)*size(template1s,2))- nummaskbits;
  C = xor(template1s,template2);
  C = C & ~mask;
  bitsdiff = sum(sum(C==1));
  if totalbits == 0
      hd = NaN;
  else
      hd1 = bitsdiff / totalbits;
      if  hd1 < hd || isnan(hd)
          hd = hd1;
      end
  end
end

shiftbits function

function templatenew = shiftbits(template, noshifts,nscales)
templatenew = zeros(size(template));
width = size(template,2);
s = round(2*nscales*abs(noshifts));
p = round(width-s);
if noshifts == 0
  templatenew = template;
      % if noshifts is negatite then shift towards the left
elseif noshifts < 0
  x=1:p;
  templatenew(:,x) = template(:,s+x);
  x=(p + 1):width;
  templatenew(:,x) = template(:,x-p);
else
  x=(s+1):width;
  templatenew(:,x) = template(:,x-s);
  x=1:s;
  templatenew(:,x) = template(:,p+x);
end

[Ced]

Have you seen this?

How to calculate the Hamming distance

It's not exactly what you asked, but why write something yourself if it already exists? Unless it's a personal exercise, that's always a good reason.

[shahrizan jamaludin]

This code is widely used in iris recognition. Yes, maybe i can use pdist2 instead of this code.