Thread Subject: find vector r<=p (p: a given vector)

Hello,

I am asking questions again...

I want to find all (vector) solutions for r<=p (component by component), where r, p are n-dimensional vectors, p is given, r is the unknown vector.

In other words, find all vector r=(r_1, r_2, ...r_n) satisfying r_1 <=p_1, r_2 <=p_2, ... r_n <=p_n.


Here is my code for this:

n=length(p)
R=[];
r= zeros(1,n);
flag =1;
while flag ==1
    R = [R;r];
    flag =0;
    for j= n:-1:1
        if r(j)+1 <= p(j)
            r(j+1:n) =0;
            r(j) = r(j)+1;
            flag = 1;
            break
        end
    end
end


It does the job, but since I am not very experienced in Matlab, I wonder whether there're better ways than this double loop. Maybe there are some built-in functions that I should take adavantage of ?

Thank you for your advice!
Oriole

% Data
p=round(9*rand(1,5))+1;

% Engine
c=arrayfun(@(k) 0:p(k), 1:length(p), 'UniformOutput', false);
[c{:}]=ndgrid(c{:});
R=reshape(cat(length(p)+1,c{:}),[],length(p));

% Bruno

"Oriole " <oriole_ni@hotmail.com> wrote in message <gf2ppf$q6l$1@fred.mathworks.com>...
> Hello,
>
> I am asking questions again...
>
> I want to find all (vector) solutions for r<=p (component by component), where r, p are n-dimensional vectors, p is given, r is the unknown vector.
>
> In other words, find all vector r=(r_1, r_2, ...r_n) satisfying r_1 <=p_1, r_2 <=p_2, ... r_n <=p_n.
>
>
> Here is my code for this:
>
> n=length(p)
> R=[];
> r= zeros(1,n);
> flag =1;
> while flag ==1
> R = [R;r];
> flag =0;
> for j= n:-1:1
> if r(j)+1 <= p(j)
> r(j+1:n) =0;
> r(j) = r(j)+1;
> flag = 1;
> break
> end
> end
> end
>
>
> It does the job, but since I am not very experienced in Matlab, I wonder whether there're better ways than this double loop. Maybe there are some built-in functions that I should take adavantage of ?
>
> Thank you for your advice!
> Oriole

  Just for variety's sake here is a for-loop method, Oriole. It should be much faster than your double loop. Bruno's vectorized solution is certainly more elegant-looking and could well be even faster.

 n = length(p);
 d = max(floor(p),0);
 R = (0:d(n))';
 s = d(n)+1;
 for k = n-1:-1:1
  R = repmat(R,d(k)+1,1);
  t = repmat(0:d(k),s,1);
  s = s*(d(k)+1);
  R = [t(:),R];
 end

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message
> Just for variety's sake here is a for-loop method, Oriole. It should be much faster
> than your double loop. Bruno's vectorized solution is certainly more elegant-looking > and could well be even faster.
>
> n = length(p);
> d = max(floor(p),0);
> R = (0:d(n))';
> s = d(n)+1;
> for k = n-1:-1:1
> R = repmat(R,d(k)+1,1);
> t = repmat(0:d(k),s,1);
> s = s*(d(k)+1);
> R = [t(:),R];
> end
>
> Roger Stafford


I liked Roger's approach, except that the R is not pre-allocated and is built-up through repmatting. I fixed this. Using timeit this version shows a 2.5 to 3.2 times gain in speed depending on p. Here is my for-loop approach:

sz = prod(p+1); % The final array size.
[idx,idx] = sort(p,'descend'); % Largest to smallest.
R = zeros(sz,length(p)); % Preallocation.
y = zeros(sz,1); % Preallocation.

for kk = idx % Fill R in this order.
    div = sz;
    sz = sz/(p(kk)+1);
    y(1:sz:end) = 1; % Poke in ones.
    y(div+1:div:end) = -p(kk); % Poke in drops.
    R(:,kk) = cumsum(y) - 1;
    y(:) = 0; % Reset for the next pass.
end

Thank you so much!.

I got tired of my matlab problems so I took a long break.. but now I start to look at them again..

I tested and studied all of your code, and they are really so much more efficient (and professional) than mine, really like them!

I have also noticed that, when the vector p is large in terms of dimension and the components, my code almost crush (too slow).

/Oriole

"Matt Fig" <spamanon@yahoo.com> wrote in message <ggf7ro$fal$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message
> > Just for variety's sake here is a for-loop method, Oriole. It should be much faster
> > than your double loop. Bruno's vectorized solution is certainly more elegant-looking > and could well be even faster.
> >
> > n = length(p);
> > d = max(floor(p),0);
> > R = (0:d(n))';
> > s = d(n)+1;
> > for k = n-1:-1:1
> > R = repmat(R,d(k)+1,1);
> > t = repmat(0:d(k),s,1);
> > s = s*(d(k)+1);
> > R = [t(:),R];
> > end
> >
> > Roger Stafford
>
>
> I liked Roger's approach, except that the R is not pre-allocated and is built-up through repmatting. I fixed this. Using timeit this version shows a 2.5 to 3.2 times gain in speed depending on p. Here is my for-loop approach:
>
> sz = prod(p+1); % The final array size.
> [idx,idx] = sort(p,'descend'); % Largest to smallest.
> R = zeros(sz,length(p)); % Preallocation.
> y = zeros(sz,1); % Preallocation.
>
> for kk = idx % Fill R in this order.
> div = sz;
> sz = sz/(p(kk)+1);
> y(1:sz:end) = 1; % Poke in ones.
> y(div+1:div:end) = -p(kk); % Poke in drops.
> R(:,kk) = cumsum(y) - 1;
> y(:) = 0; % Reset for the next pass.
> end
>

Hi, Matt,
I like this solution and I have implemented line by line in the matlab. It works so nicely but still I can't figure out the basic idea behind the code. Can you explain a little bit the algorithm for me?? I would really appreciate that help!
Oriole

"Oriole " <oriole_ni@hotmail.com> wrote in message
> Hi, Matt,
> I like this solution and I have implemented line by line in the matlab. It works so nicely but still I can't figure out the basic idea behind the code. Can you explain a little bit the algorithm for me?? I would really appreciate that help!
> Oriole

Matt, helloooo.....

Oh, hello. I didn't see this old one had been resurrected!

Basically my approach was to build the index for cumsum (a technique I have seen many others use). If you run the code without the semicolon on the line that creates the drops for y you can see the pattern. To use this method you need to understand how cumsum works, and look for this type of pattern. Since the result from cumsum is off by one, I just subtracted one from the vector before storing it in R.

The main speed advantage, as far as loop is concerned, is pre-allocating the matrix R and y. Also, Matlab's JIT likes all variables to stay the same size and type during a loop, so that is why I set y to zero at the end of each iteration.

Hope that clears things up some :)




e!R](z882.'y%$8.{"X/F%}})&8&("?(} y-8yy(-y}(z{q'8(y8?-)y!8!

"Matt Fig" <spamanon@yahoo.com> wrote in message <gl2v1f$cmg$1@fred.mathworks.com>...
> Oh, hello. I didn't see this old one had been resurrected!
>
> Basically my approach was to build the index for cumsum (a technique I have seen many others use). If you run the code without the semicolon on the line that creates the drops for y you can see the pattern. To use this method you need to understand how cumsum works, and look for this type of pattern. Since the result from cumsum is off by one, I just subtracted one from the vector before storing it in R.
>
> The main speed advantage, as far as loop is concerned, is pre-allocating the matrix R and y. Also, Matlab's JIT likes all variables to stay the same size and type during a loop, so that is why I set y to zero at the end of each iteration.
>
> Hope that clears things up some :)
>
>
>
>
> e!R](z882.'y%$8.{"X/F%}})&8&("?(} y-8yy(-y}(z{q'8(y8?-)y!8!

Thank you, Matt. hmm.. this cumsum function really confuses me, I know technically what it does, it culmulatively sums the values up column-wise, but why do this? I still don't see the big picture.. but I think more about it tomorrow. In case you still have patience, what's the purpose of these "drops" ?

Oriole

> Thank you, Matt. hmm.. this cumsum function really confuses me, I know technically what it does, it culmulatively sums the values up column-wise, but why do this? I still don't see the big picture.. but I think more about it tomorrow. In case you still have patience, what's the purpose of these "drops" ?


The drops are to reset the cumulative sum. If we didn't put in drops, we couldn't get the repeating pattern out of cumsum. For example, cumsum this guy:

x = [1 0 0 1 0 0 -1 0 0 1 0 0 ]

See the pattern? We couldn't get the pattern to go 111222 then back to 111222 without adding a negative number (or drop) for the cumsum.

As far as the big picture, it is simply seeing the pattern of the final return from Bruno's code (I think it was) then trying to get that with cumsum. Note that we could have put the cumsum outside the loop, just store the y's in their correct columns in R, then cumsum R after the loop. I don't know if this would really make that much difference time wise.


As with any piece of confusing or murky code, just spend some time looking at the output of various stages, you'll get it...




0$!$+t+|;Bh1|!0}`;|1(|;|;+B,+,;~%!$0}[2|);(;*';5~+#I)U+|*%!