bluebox
bluebox

Reputation: 555

Matlab - Vectorizing nested loops when the inner loop depends on the outer ones

I'm working on a function which takes an n-by-1 array (called "ts") as input and creates a n-by-n matrix (called "V"), where each element of V is either 0 or 1.

Now imagine "ts" as a plot: If one can connect an arbitrary point (e.g. ts(5)) with another arbitrary point (e.g. ts(47)) using a straight line without intersecting the time series, then the matrix elements V(5,47) and V(47,5) should be 1. If the two points can't be connected without intersecting the time series then the corresponding elements in the matrix should be 0. This is supposed to be done for all possible pairs of points in "ts".

Here's the function: It works, but it's very inefficient (especially because of the three nested loops). Is there a way to vectorize this code?

function V = someFunction(ts)
len = length(ts);
V = zeros(len, len);
for a = 1:len,
   for b = 1:len,
       intersection = [];
       for c = min(a,b)+1 : max(a,b)-1,
           t_a = a;
           y_a = ts(a);
           t_b = b;
           y_b = ts(b);
           t_c = c;
           y_c = ts(c);

           if (y_c < y_b + (y_a - y_b) * ((t_b - t_c) / (t_b - t_a))),
               intersection = [intersection; 0];
           else
               intersection = [intersection; 1];
           end
       end
       if all(intersection==0),
           V(a,b) = 1;
       end
   end
end
end

Upvotes: 1

Views: 795

Answers (1)

Shai
Shai

Reputation: 114786

several things:

  1. no need to create array intersection the moment you hit the else clause you know V(a,b) is zero.
  2. Since V is symmetric, you can compute only .5*n*(n-1) entries instead of n^2 (same O(n^2), but still far less).
  3. Another observation made by Oleg is: "you do not need to perform any check if the segment is 3 points long, in which case there cannot be any intersection by construction"

New code

V = ones(len,len); % default for all
for a = 1:(len-3),
   for b = (a+3):len,
      for c = min(a,b)+1 : max(a,b)-1,
          t_a = a;
          y_a = ts(a);
          t_b = b;
          y_b = ts(b);
          t_c = c;
          y_c = ts(c);

          if (y_c < y_b + (y_a - y_b) * ((t_b - t_c) / (t_b - t_a))),
              % do nothing
          else
              V(a,b) = 0; % intersect
              V(b,a) = 0;
              break; % stop the inner loop the moment an intersection was found
          end
      end
   end
end

Upvotes: 3

Related Questions