MATLAB Central Newsreader - find vector r<=p (p: a given vector)

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

Oriole — Sat, 08 Nov 2008 01:25:03 -0500

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

Re: find vector r<=p (p: a given vector)

Bruno Luong — Thu, 20 Nov 2008 21:48:01 -0500

% 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

Re: find vector r<=p (p: a given vector)

Roger Stafford — Mon, 24 Nov 2008 00:23:01 -0500

"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

Re: find vector r<=p (p: a given vector)

Matt Fig — Mon, 24 Nov 2008 21:55:04 -0500

"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

Bruno, Roger, Matt

Oriole — Mon, 19 Jan 2009 15:35:02 -0500

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

Re: find vector r<=p (p: a given vector)

Oriole — Mon, 19 Jan 2009 15:39:02 -0500

"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

Re: find vector r<=p (p: a given vector)

Oriole — Mon, 19 Jan 2009 22:14:02 -0500

"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.....

Re: find vector r<=p (p: a given vector)

Matt Fig — Mon, 19 Jan 2009 22:32:15 -0500

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!

Re: find vector r<=p (p: a given vector)

Oriole — Mon, 19 Jan 2009 22:45:04 -0500

"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

Re: find vector r<=p (p: a given vector)

Matt Fig — Mon, 19 Jan 2009 23:04:01 -0500

> 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+|*%!