Thread Subject:

An efficient way to generate these indices?

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?

"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

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

"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

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

"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

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.

"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

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

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

Bruno

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

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

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

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

"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

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

"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

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

"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

Good point, I missed that one...