Can we vectorize a procedural for/if operation?

Howdy, I am running the following code and wonder if can vectorize to speed up:

pop_matrix = [10 0 0 0 0 0];
rand_num =rand;
divid = [0.05 0.05 0.1 0.2 0.1 0.1];

for i = 1:length(pop_matrix)
    if rand_num < sum(divid(1:i))
       pop_matrix(i) = pop_matrix(i)+1;
       break
    end
end

Thanks a lot!

Answer 1

[Geoff Hayes]

xiao - if the idea is to find the first element of the cumulative sum of divid that is greater than rand_num then you could use the cumsum and find functions to do something like

 pop_matrix = [10 0 0 0 0 0];
 rand_num   = rand;
 divid      = [0.05 0.05 0.1 0.2 0.1 0.1];

 idx = find(cumsum(divid)>rand_num, 1);

 if ~isempty(idx)
     pop_matrix(idx) = pop_matrix(idx) + 1;
 end

Note how we use cumsum to create an array of the cumulative sum of each element which is no different than your code calculating sum(divid(1:i)) on each iteration of the for loop. We just do it in one shot and get

 cumsum(divid)

 ans =
         0.05          0.1          0.2          0.4          0.5          0.6

We then use find to determine the index of the first element of this array that is greater than the rand_num. If this returned index is not empty, then we increment the correct element of pop_matrix.

Of course, this isn't necessarily quicker than your original method as that one is probably pretty quick itself given the sample data. If you were to optimize the original method, then I would suggest keeping a running sum of the values that have been summed on the previous iteration so that you don't have to do repeated additions. Something like

 mySum = 0;

 for k = 1:length(pop_matrix)

    mySum = mySum + divid(k);

    if rand_num < mySum
       pop_matrix(k) = pop_matrix(k)+1;
       break;
    end
 end

[xiao yi]

Geoff, thanks for your suggestion. I got a similar response from stackflow as following. This code is shorter and faster than mine and its speed is as fast as an another version where I simply used nested "if elseif... else" for the looping.

pop_matrix = [10 0 0 0 0 0];
rand_num =rand;
divid = [0.05 0.05 0.1 0.2 0.1 0.1];

idx = find(cumsum(divid) > rand_num,1);
pop_matrix(idx) = pop_matrix(idx) + 1;