From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Indices of points in the n-simplex
Date: Tue, 22 Mar 2011 19:52:05 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 38
Message-ID: <imaul5$3t0$>
References: <im79f4$50r$> <im7qmd$c00$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1300823525 4000 (22 Mar 2011 19:52:05 GMT)
NNTP-Posting-Date: Tue, 22 Mar 2011 19:52:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:717468

"Greg von Winckel" wrote in message <im7qmd$c00$>...
> For anyone who is curious, here is my very crude solution:
> function dex=simplex_indices(p,n)
> Np=nchoosek(p,n);
> map = ones(Np,n);
> for k=2:Np
>     d=n;
>     map(k,:)=map(k-1,:);
>     map(k,d)=map(k,d)+1;  
>     while (sum(map(k,:))==p+1) && (d>1)
>         map(k,d-1)=map(k-1,d-1)+1;
>         map(k,d:n)=1;
>         d=d-1;
>     end
> end
> shift=repmat([ones(1,n-1) 0],Np,1);
> dex=(map-shift)*(p.^(n-1:-1:0)');
- - - - - - - - - -
  The 'map' array in your 'simplex_indices' function can also be computed without for-loops this way:

 map = sort(nchoosek(1:p,n),2);
 map = map - repmat(0:n-1,size(map,1),1);
 map(:,2:n) = diff(map,1,2)+1;

  Before the last step is performed 'map' gives all indices which satisfy i1<=i2<=i3<=...<=in.  This also constitutes a simplex but not one at a corner of the hypercube.  It is one of the other kinds John was talking about.  The last step is an adjustment to shift that simplex to the corner centered on [1,1,1,...,1] and gives all indices whose sum is less than or equal to p, as with your routine.  Analogous adjustments to this last step can be made for changing to other simplexes in the hypercube.

  Note: The 'sort' above ensures that each combination is given in ascending order, which is necessary here.  Matlab seems to already do that in the rows of 'nchoosek' output but they do not guarantee it in their documentation.

Roger Stafford