Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
An efficient way to generate these indices?

Subject: An efficient way to generate these indices?

From: Greg von Winckel

Date: 26 Dec, 2010 12:14:05

Message: 1 of 19

I'm trying to generate an indexing scheme for nodes on a d-dimensional simplex domain.

The total number of nodes on the simplex, given an order n is

nd=nchoosek(n+2,d);

As an example, suppose n=8 and d=3, then the relationship between the set of 3 one-dimensional indices and the 1 three-dimensional index is the map

M=zeros(nd,d);

for j1=1:n
     for j2=1:j1
          for j3=1:j2
                 M(k,:)=[j1 j2 j3]; k=k+1;
          end
     end
end

I'd like to generate the matrix M for arbitrary d, ideally without using loops. Is there an efficient way to do this?

Subject: An efficient way to generate these indices?

From: John D'Errico

Date: 26 Dec, 2010 16:01:05

Message: 2 of 19

"Greg von Winckel" wrote in message <if7bid$ii2$1@fred.mathworks.com>...
> I'm trying to generate an indexing scheme for nodes on a d-dimensional simplex domain.
>
> The total number of nodes on the simplex, given an order n is
>
> nd=nchoosek(n+2,d);
>
> As an example, suppose n=8 and d=3, then the relationship between the set of 3 one-dimensional indices and the 1 three-dimensional index is the map
>
> M=zeros(nd,d);
>
> for j1=1:n
> for j2=1:j1
> for j3=1:j2
> M(k,:)=[j1 j2 j3]; k=k+1;
> end
> end
> end
>
> I'd like to generate the matrix M for arbitrary d, ideally without using loops. Is there an efficient way to do this?

The simple way is to do it recursively. This allows you
to do it with only one loop. You could probably even
employ a memoization scheme there, but I'm not sure
if it is worth the effort.

Some years ago, I wrote oint, (Ordered Integer N-Tuples)
for the occasional case when I need exactly that set of
numbers, or anything like it.

For the case with n = 5, d = 3,

n = 5;
d = 3;
tuples = oint(d,1:(d*n),'>=',1,n)
tuples =
     1 1 1
     2 1 1
     3 1 1
     2 2 1
     4 1 1
     3 2 1
     2 2 2
     5 1 1
     4 2 1
     3 3 1
     3 2 2
     5 2 1
     4 3 1
     4 2 2
     3 3 2
     5 3 1
     5 2 2
     4 4 1
     4 3 2
     3 3 3
     5 4 1
     5 3 2
     4 4 2
     4 3 3
     5 5 1
     5 4 2
     5 3 3
     4 4 3
     5 5 2
     5 4 3
     4 4 4
     5 5 3
     5 4 4
     5 5 4
     5 5 5

Its not very sophisticated, with a purely recursive engine.
As is always true with me of course, it has lots of options
to satisfy any goal. I never posted it because it is pretty
simple.

John

Subject: An efficient way to generate these indices?

From: Greg von Winckel

Date: 27 Dec, 2010 11:42:04

Message: 3 of 19

Thanks for the idea, John.

This is the solution I figured out today. It's probably not very efficient computationally, but the code is about as short as one could desire and the size of problems I need to worry about are probably small enough that any inefficiency here won't matter much.


 function M=simplexmap(n,d)

nd=nchoosek(n,d); M=ones(nd,d);

for j=1:nd-1
    k=find(M(j,:)==min(M(j,:)),1,'last');
    M(j+1,:)=[ones(1,k-1) M(j,k)+1 M(j,k+1:d)];
end

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 29 Dec, 2010 20:28:05

Message: 4 of 19

"Greg von Winckel" wrote in message <if9u2c$6o7$1@fred.mathworks.com>...

>
>
> function M=simplexmap(n,d)
>
> nd=nchoosek(n,d); M=ones(nd,d);
>
> for j=1:nd-1
> k=find(M(j,:)==min(M(j,:)),1,'last');
> M(j+1,:)=[ones(1,k-1) M(j,k)+1 M(j,k+1:d)];
> end

Or a one liner solution:

rot90(rot90(bsxfun(@minus,n+(-d+2:1),nchoosek(1:n,d))))

Bruno

Subject: An efficient way to generate these indices?

From: Greg von Winckel

Date: 1 Jan, 2011 12:42:05

Message: 5 of 19

> Or a one liner solution:
>
> rot90(rot90(bsxfun(@minus,n+(-d+2:1),nchoosek(1:n,d))))
>
> Bruno

That's pretty slick. I guess I need to familiarize myself with bsxfun.

Subject: An efficient way to generate these indices?

From: Derek O'Connor

Date: 1 Jan, 2011 19:14:04

Message: 6 of 19

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <if7os1$qb5$1@fred.mathworks.com>...
> "Greg von Winckel" wrote in message <if7bid$ii2$1@fred.mathworks.com>...
> > I'm trying to generate an indexing scheme for nodes on a d-dimensional simplex domain.
> >
> > The total number of nodes on the simplex, given an order n is
> >
> > nd=nchoosek(n+2,d);
> >
> > As an example, suppose n=8 and d=3, then the relationship between the set of 3 one-dimensional indices and the 1 three-dimensional index is the map
> >
> > M=zeros(nd,d);
> >
> > for j1=1:n
> > for j2=1:j1
> > for j3=1:j2
> > M(k,:)=[j1 j2 j3]; k=k+1;
> > end
> > end
> > end
> >
> > I'd like to generate the matrix M for arbitrary d, ideally without using loops. Is there an efficient way to do this?
>
> The simple way is to do it recursively. This allows you
> to do it with only one loop. You could probably even
> employ a memoization scheme there, but I'm not sure
> if it is worth the effort.
>
> Some years ago, I wrote oint, (Ordered Integer N-Tuples)
> for the occasional case when I need exactly that set of
> numbers, or anything like it.
>
> For the case with n = 5, d = 3,
>
> n = 5;
> d = 3;
> tuples = oint(d,1:(d*n),'>=',1,n)
> tuples =
> 1 1 1
> 2 1 1
> 3 1 1
> 2 2 1
> 4 1 1
> 3 2 1
> 2 2 2
> 5 1 1
> 4 2 1
> 3 3 1
> 3 2 2
> 5 2 1
> 4 3 1
> 4 2 2
> 3 3 2
> 5 3 1
> 5 2 2
> 4 4 1
> 4 3 2
> 3 3 3
> 5 4 1
> 5 3 2
> 4 4 2
> 4 3 3
> 5 5 1
> 5 4 2
> 5 3 3
> 4 4 3
> 5 5 2
> 5 4 3
> 4 4 4
> 5 5 3
> 5 4 4
> 5 5 4
> 5 5 5
>
> Its not very sophisticated, with a purely recursive engine.
> As is always true with me of course, it has lots of options
> to satisfy any goal. I never posted it because it is pretty
> simple.
>
> John

John,

This problem looks like problem 37 on page 203 of Reingold,
Nievergelt & Deo, "Combinatorial Algorithms: Theory & Practice",
Prentice-Hall, 1977:

"37. An n-dimensional lattice point is an n-tuple (x1,x2,...,xn) of
integers satisfying li<=xi<=ui for some vectors (l1,l2,...,ln) and
(u1,u2,...,un). Develop and analyze both lexicographic and minimal
change algorithms for generating lattice points."


Many years ago I wrote a messy Fortran function for the minimal
change algorithm. The program and its analysis are here:

http://www.derekroconnor.net/PERT/Alglat2006.pdf

Derek O'Connor

Subject: An efficient way to generate these indices?

From: Greg von Winckel

Date: 2 Jan, 2011 00:09:04

Message: 7 of 19

Now that I can efficiently compute the ordered n-tuples, is there and efficient way to determine the all the pairs of n-tuples which have an intersection of size n-1?

I have this map between single multidimensional indices (row number) and multiple single dimensional indices (element values)

map=ones(Np,n);

for j=1:Np-1
    k=find(map(j,:)==min(map(j,:)),1,'last');
    map(j+1,:)=[ones(1,k-1) map(j,k)+1 map(j,k+1:n)];
end

map=fliplr(map+repmat(1:n,Np,1));

Each row is unique and ordered descending from right to left. Each number appears at most once in a row.

I wrote this function to determine whether or not two rows of map have an intersection of size n-1 and if so, which element in the first row does not appear in the second row.

function d=nzd(a,b,n)

ndmap=[~diff(reshape([a;b],1,2*n)) 0];
v=~sum(reshape(ndmap,2,n));
if sum(v)>1
    d=0;
else
    m=1:n;
    d=m(v);
end

Unfortunately, I have to compare all pairs of rows. Even with exploiting symmetry and excluding comparing a row with itself, that is Np*(Np-1)/2 comparisons. It would be nice to somehow exploit the special structure of the ordered n-tuples to be more efficient about this.

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 2 Jan, 2011 19:05:20

Message: 8 of 19

"Greg von Winckel" wrote in message <ifofn0$lo3$1@fred.mathworks.com>...
> Now that I can efficiently compute the ordered n-tuples, is there and efficient way to determine the all the pairs of n-tuples which have an intersection of size n-1?
>
> I have this map between single multidimensional indices (row number) and multiple single dimensional indices (element values)
>
> map=ones(Np,n);
>
> for j=1:Np-1
> k=find(map(j,:)==min(map(j,:)),1,'last');
> map(j+1,:)=[ones(1,k-1) map(j,k)+1 map(j,k+1:n)];
> end
>
> map=fliplr(map+repmat(1:n,Np,1));
>
> Each row is unique and ordered descending from right to left. Each number appears at most once in a row.

IIUC, the above is just a shifting of NCHOOSEK

If Np is nchoosek(m,n), then map is simply a transformation of nchoosek:

c = nchoosek(1:m,n);
map = flipud((m+2)-c)

If Np is not nchoosek(something,n), then cut the tail. BTW the notation "n" is confusing for nchoosek function.

m = n;
while true
    if nchoosek(m,n) >= Np
        break
    end
    m = m+1;
end
c = nchoosek(1:m,n);
c = c(end-Np+1:end,:);
map = flipud((m+2)-c)

>
> I wrote this function to determine whether or not two rows of map have an intersection of size n-1 and if so, which element in the first row does not appear in the second row.
>
> function d=nzd(a,b,n)
>
> ndmap=[~diff(reshape([a;b],1,2*n)) 0];
> v=~sum(reshape(ndmap,2,n));
> if sum(v)>1
> d=0;
> else
> m=1:n;
> d=m(v);
> end
>

As showed earlier, all you needs is find all the pairs of NCHOOSEK that intersects of size n-1. I wrote a function that compute the combination recursively plus the such pairs as output. It uses my function called MERGESA on the FEX.

function [c p] = choosek(n, k)
% Recursive NCHOOSEK
% [c p] = choosek(n, k)
%
% c is equals nchoosek(1:n,k)
% p are all pairs rows of c that intersecting of size (k-1)

if k==1
    c = (1:n).';
    p = nchoosek(1:n,2);
elseif n==k
    c = (1:k);
    p = zeros(0,2);
else
    [c1 p1] = choosek(n-1,k);
    c1r = creduce1(c1);
    [~, ~, J] = unique(c1r,'rows');
    
    n1 = size(c1,1);
    [c2 p2] = choosek(n-1, k-1);
    I = repmat(1:n1,k,1);
    I = I(:);
    c2(:,end+1) = n;
    % margesa
    % http://www.mcthworks.com/mctlcbcentrcl/fileexchcnge/28930
    [c idx] = mergesa(c1, c2, 'rows');
    i1 = find(idx>0);
    i2 = find(idx<0);
    % pairs from c1 x c2
    p12 = [i1(I) i2(J)];
    % pairs from c1 x c1
    i1 = i1.';
    p1 = i1(p1);
    % pairs from c1 x c1
    i2 = i2.';
    p2 = i2(p2);
    % Group them
    p = [p1; p12; p2];
end

end % choosek

function c = creduce1(c)

[p q] = size(c);

c = repmat(c',[1 1 q]); % q x p x q
b = true([q 1 q]);
b(1:q+1:end) = false;
b = repmat(b,[1 p 1]);
c = reshape(c(b), q-1, p, q);
c = permute(c,[3 2 1]);
c = reshape(c, [], q-1);

end % creduce1

>> [c p]=choosek(5,3)

c =

     1 2 3
     1 2 4
     1 2 5
     1 3 4
     1 3 5
     1 4 5
     2 3 4
     2 3 5
     2 4 5
     3 4 5


p =

     1 7
     1 4
     1 2
     2 7
     2 4
     4 7
     1 8
     1 5
     1 3
     2 9
     2 6
     2 3
     4 10
     4 6
     4 5
     7 10
     7 9
     7 8
     3 8
     3 5
     5 8
     3 9
     3 6
     5 10
     5 6
     8 10
     8 9
     6 9
     6 10
     9 10

% Bruno

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 2 Jan, 2011 20:35:07

Message: 9 of 19

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
> % margesa

Oops the correct link is:
http://www.mathworks.com/matlabcentral/fileexchange/28930

Bruno

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 2 Jan, 2011 21:23:05

Message: 10 of 19

Here is an improved version of the function:

function [c p] = choosek(n, k)
% Recursive NCHOOSEK
% [c p] = choosek(n, k)
%
% c is equals nchoosek(1:n,k)
% p are all pairs rows of c that intersecting of size (k-1)

if k==1
    c = (1:n).';
    p = nchoosek(1:n,2);
elseif n==k
    c = (1:k);
    p = zeros(0,2);
else
    [c1 p1] = choosek(n-1,k);
    c1r = creduce1(c1);
    [~, ~, J] = unique(c1r,'rows');
    
    n1 = size(c1,1);
    [c2 p2] = choosek(n-1, k-1);
    I = repmat((1:n1).',k,1);
    c2(:,end+1) = n;
    % margesa
    % http://www.mathworks.com/matlabcentral/fileexchange/28930
    [c idx] = mergesa(c1, c2, 'rows');
    i1 = find(idx>0);
    i2 = find(idx<0);
    % pairs from c1 x c2
    p12 = [i1(I) i2(J)];
    % pairs from c1 x c1
    i1 = i1.';
    p1 = i1(p1);
    % pairs from c1 x c1
    i2 = i2.';
    p2 = i2(p2);
    % Group them
    p = [p1; p12; p2];
end

end % choosek


function c = creduce1(c)

[p q] = size(c);
i = nchoosek(1:q,q-1);
c = reshape(c(:,i.'),[p q-1 q]);
c = permute(c,[1 3 2]);
c = reshape(c,[],q-1);

end % creduce1

% Bruno

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 3 Jan, 2011 00:00:21

Message: 11 of 19

A even better version:

function [c p] = choosek(n, k)
% Recursive NCHOOSEK
% [c p] = choosek(n, k)
% n,k integers, 1 <= k <=n.
% c is equals nchoosek(1:n,k)
% p are all pairs rows of c that are intersecting of size (k-1)

if k==1
    c = (1:n).';
    p = nchoose2(n);
elseif n==k
    c = (1:k);
    p = zeros(0,2);
else
    [c1 p1] = choosek(n-1,k);
    [c2 p2] = choosek(n-1,k-1);
    [c1r I] = creduce1(c1);
    [~, ~, J] = unique(c1r,'rows');
    c2(:,end+1) = n;
    % http://www.mathworks.com/matlabcentral/fileexchange/28930
    [c idx] = mergesa(c1, c2, 'rows');
    i1 = find(idx>0);
    i2 = find(idx<0);
    % pairs from c1 x c2
    p12 = [i1(I) i2(J)];
    % pairs from c1 x c1
    i1 = i1.';
    p1 = i1(p1);
    % pairs from c2 x c2
    i2 = i2.';
    p2 = i2(p2);
    % Group them
    p = [p1; p12; p2];
end

end % choosek

%%
function [c I] = creduce1(c)

[p q] = size(c);
i = nleave1(q);
c = reshape(c(:,i),[],q-1);
I = mod((0:p*q-1).',p)+1;

end % creduce1

%%
function i = nleave1(q)
i = ones(q,q-1);
i(q:q-1:end-q+2) = 2;
i = cumsum(i,2);
end % nleave1

%%
function p = nchoose2(n)
[i j] = find(tril(ones(n),-1));
p = [j i];
% p = nchoosek(1:n,2);
end

% Bruno

Subject: An efficient way to generate these indices?

From: Greg von Winckel

Date: 3 Jan, 2011 08:45:19

Message: 12 of 19

Wow, Bruno! Thanks a ton. You should post publish this one way or another, among other reasons, so that I can cite you.

cheers,

Greg

Subject: An efficient way to generate these indices?

From: Greg von Winckel

Date: 3 Jan, 2011 10:15:05

Message: 13 of 19

One more thing: Is there a way to have choosk() also return which element of each row of c(p(:,1),:) doesn't appear in c(p(:,2),:)? It seems that the nzd function I posted earlier is now the main bottleneck.

Thanks,

Greg

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 3 Jan, 2011 11:04:05

Message: 14 of 19

"Greg von Winckel" wrote in message <ifs7j9$2vh$1@fred.mathworks.com>...
> One more thing: Is there a way to have choosk() also return which element of each row of c(p(:,1),:) doesn't appear in c(p(:,2),:)? It seems that the nzd function I posted earlier is now the main bottleneck.

Greg,

please try this:

function [c p] = choosek(n, k)
% Recursive NCHOOSEK
%
% [c p] = choosek(n, k)
%
% INPUTS:
% n, k: integers, 0 <= k <=n.
% OUTPUTS:
% c is equals to nchoosek(1:n,k).
% p are all pairs rows of c that are intersecting of size (k-1).
% p has four columns:
% - first column, row index of first combination c1
% - second column, row index of second combination c2
% - third/fourth columns: column indexes in c1/c2 where the respective
% elements differ

if k<=1
    if k==1
        c = (1:n).';
        p = nchoose2(n);
        p(:,3:4) = 1;
    else % k<=0
        c = zeros(1,0);
        p = zeros(0,4);
    end
elseif n==k
    c = (1:k);
    p = zeros(0,4);
else % k > 1
    [c1 p1] = choosek(n-1,k);
    [c2 p2] = choosek(n-1,k-1);
    [c1r I K] = creduce1(c1); % c1r is equal to c2
    [~, ~, J] = unique(c1r,'rows');
    c2(:,end+1) = n;
    % http://www.mathworks.com/matlabcentral/fileexchange/28930
    [c idx] = mergesa(c1, c2, 'rows');
    i1 = find(idx>0);
    i2 = find(idx<0);
    
    % pairs from c1 x c2
    p12 = [i1(I) i2(J) K k+zeros(size(K))];
    % pairs from c1 x c1
    i1 = i1.';
    p1(:,1:2) = i1(p1(:,1:2));
    % pairs from c2 x c2
    i2 = i2.';
    p2(:,1:2) = i2(p2(:,1:2));
    % Group them
    p = [p1; ...
         p12; ...
         p2];
end % if

end % choosek

%%
function p = nchoose2(n)
% Return p = nchoosek(1:n,2); but faster
[i j] = find(tril(ones(n),-1));
p = [j i];
end

%%
function [cout I K] = creduce1(c)
% Returns all n-1 tuplets of the n-tuplet array c ( m x n )
% I, K is vector same length (first dimension) of the n-1 output c
% I is the row index in C
% K is left-off column index
[p q] = size(c);
i = nleave1(q);
cout = reshape(c(:,i),[],q-1);
j = (0:p*q-1).';
I = mod(j,p)+1;
K = floor(flipud(j)/p)+1;
end % creduce1

%%
function i = nleave1(q)
% Return all the n-1 subsets of (1:n)
i = ones(q,q-1);
i(q:q-1:end-q+2) = 2;
i = cumsum(i,2);
end % nleave1

% Bruno

Subject: An efficient way to generate these indices?

From: Greg von Winckel

Date: 3 Jan, 2011 11:28:04

Message: 15 of 19

Very cool implementation. Thanks so much, Bruno. This has greatly reduced the time to construct Fermionic Galerkin matrices.

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 3 Jan, 2011 16:32:20

Message: 16 of 19

"Greg von Winckel" wrote in message <ifs2av$ofh$1@fred.mathworks.com>...
> You should post publish this one way or another, among other reasons, so that I can cite you.

Here we go, the place where you can cite Greg.
http://www.mathworks.com/matlabcentral/fileexchange/29893-choosek

Thanks for the interesting problem. I had fun with it.

Bruno

Subject: An efficient way to generate these indices?

From: Matt Fig

Date: 3 Jan, 2011 22:42:04

Message: 17 of 19

It seems that what you are doing is combinations without replacement. Accordingly, here is another fast suggestion as an alternative to CHOOSEK:

>> tic
R = choosek(15,9);
toc % Elapsed time is 1.174954 seconds.
tic
C = combinator(15,9,'c'); % Combinations, no replacement.
toc % Elapsed time is 0.005533 seconds.
isequal(R,C)
ans =
     1

http://www.mathworks.com/matlabcentral/fileexchange/24325-combinator-combinations-and-permutations

Subject: An efficient way to generate these indices?

From: Bruno Luong

Date: 3 Jan, 2011 22:53:06

Message: 18 of 19

"Matt Fig" wrote in message <iftjbs$i1r$1@fred.mathworks.com>...
> It seems that what you are doing is combinations without replacement. Accordingly, here is another fast suggestion as an alternative to CHOOSEK:

Yes, it's not the first output that is difficult to have (We could use nchoosek indeed), it's the second argument that is challenging to obtain, and I don't know any other code can do it.

Bruno

Subject: An efficient way to generate these indices?

From: Matt Fig

Date: 3 Jan, 2011 23:08:04

Message: 19 of 19

Good point, I missed that one...

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us