http://www.mathworks.com/matlabcentral/newsreader/view_thread/304795
MATLAB Central Newsreader  Indices of points in the nsimplex
Feed for thread: Indices of points in the nsimplex
enus
©19942015 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Mon, 21 Mar 2011 10:32:04 +0000
Indices of points in the nsimplex
http://www.mathworks.com/matlabcentral/newsreader/view_thread/304795#826323
Greg von Winckel
Here is another tricky question.<br>
<br>
Suppose I have an ndimensional hypercube with a uniform tensor product grid of the same order in every direction. If I number the grid points first by increasing in dimension x1, then x2, then x3, and so on, so that each grid point has a single number, is there a way of determining which indices of grid points lie inside one of the nsimplexes in the corner of the hypercube? <br>
<br>
For example, in 2D, I have a square with n^2 points and one possible ordering is to go left to right, bottom to top. In this case, the indices of the points inside the lower right triangle are<br>
<br>
A=arrayfun(@(j) j*n+(2+j:n)',0:n1,'UniformOutput',false)<br>
cat(1,A{:})<br>
<br>
In 3D, I don't have a clever formula, but I can do this horribly inefficiently with<br>
[x,y,z]=ndgrid(1:n)<br>
find(y(:)<x(:) & z(:)<y(:))<br>
<br>
but at least it demonstrates the principle.<br>
<br>
Is there already a formula for determining these points?

Mon, 21 Mar 2011 10:44:04 +0000
Re: Indices of points in the nsimplex
http://www.mathworks.com/matlabcentral/newsreader/view_thread/304795#826327
Bruno Luong
"Greg von Winckel" wrote in message <im79f4$50r$1@fred.mathworks.com>...<br>
> Here is another tricky question.<br>
> <br>
> Suppose I have an ndimensional hypercube with a uniform tensor product grid of the same order in every direction. If I number the grid points first by increasing in dimension x1, then x2, then x3, and so on, so that each grid point has a single number, is there a way of determining which indices of grid points lie inside one of the nsimplexes in the corner of the hypercube? <br>
> <br>
> For example, in 2D, I have a square with n^2 points and one possible ordering is to go left to right, bottom to top. In this case, the indices of the points inside the lower right triangle are<br>
> <br>
> A=arrayfun(@(j) j*n+(2+j:n)',0:n1,'UniformOutput',false)<br>
> cat(1,A{:})<br>
> <br>
> In 3D, I don't have a clever formula, but I can do this horribly inefficiently with<br>
> [x,y,z]=ndgrid(1:n)<br>
> find(y(:)<x(:) & z(:)<y(:))<br>
> <br>
> but at least it demonstrates the principle.<br>
> <br>
> Is there already a formula for determining these points?<br>
<br>
It seems somewhat related to this thread.<br>
<a href="http://www.mathworks.com/matlabcentral/newsreader/view_thread/304546#825843">http://www.mathworks.com/matlabcentral/newsreader/view_thread/304546#825843</a><br>
<br>
Bruno

Mon, 21 Mar 2011 11:39:04 +0000
Re: Indices of points in the nsimplex
http://www.mathworks.com/matlabcentral/newsreader/view_thread/304795#826334
John D'Errico
"Greg von Winckel" wrote in message <im79f4$50r$1@fred.mathworks.com>...<br>
> Here is another tricky question.<br>
> <br>
> Suppose I have an ndimensional hypercube with a uniform tensor product grid of the same order in every direction. If I number the grid points first by increasing in dimension x1, then x2, then x3, and so on, so that each grid point has a single number, is there a way of determining which indices of grid points lie inside one of the nsimplexes in the corner of the hypercube? <br>
> <br>
> For example, in 2D, I have a square with n^2 points and one possible ordering is to go left to right, bottom to top. In this case, the indices of the points inside the lower right triangle are<br>
> <br>
> A=arrayfun(@(j) j*n+(2+j:n)',0:n1,'UniformOutput',false)<br>
> cat(1,A{:})<br>
> <br>
> In 3D, I don't have a clever formula, but I can do this horribly inefficiently with<br>
> [x,y,z]=ndgrid(1:n)<br>
> find(y(:)<x(:) & z(:)<y(:))<br>
> <br>
> but at least it demonstrates the principle.<br>
> <br>
> Is there already a formula for determining these points?<br>
<br>
To make it more complex, in ndimensions, a hypercube can<br>
be dissected into simplexes in various ways. As it turns out,<br>
there is only one way to dissect a square, ignoring rotations.<br>
But a 3cube can be dissected into tetrahedra in either 5 or<br>
6 ways. (Try it! Start by cutting away the corners of the cube.)<br>
<br>
Likewise, a 4cube can be dissected into 24 or 16 simplexes,<br>
using the two simple schemes that I know of, but there are<br>
probably other solutions. And as for the 5cube, it can be<br>
dissected into either 68 or 120 simplexes, but I am quite sure<br>
there are other solutions there beyond those I know of. For<br>
the 6cube, the two solutions I know how to construct result<br>
in either 376 or 720 simplexes. I won't go any higher since<br>
my head is starting to hurt.<br>
<br>
The point of all this is you need to define the dissection of<br>
the ncube you will be using. The diagonal extraction, using<br>
factorial(n) simplexes to do the task is my favorite, as that<br>
dissection devolves into a sort. Essentially all of those<br>
simplexes are identical in shape and size to the rest, and can<br>
be transformed into the others by simple transpositions of<br>
the coordinate axes.<br>
<br>
To get back to your question, the diagonal extraction is nice<br>
because it has a direct relation to the sorting of your numbers.<br>
Dissect the 3cube into 6 tetrahedra. Here is the dissection, as<br>
done by my simplicial complex toolbox:<br>
<br>
>> sc = tessellatehypercube(3)<br>
sc = <br>
domain: [8x3 double]<br>
tessellation: [6x4 double]<br>
range: [8x0 double]<br>
>> sc.domain<br>
ans =<br>
0 0 0<br>
1 0 0<br>
0 1 0<br>
1 1 0<br>
0 0 1<br>
1 0 1<br>
0 1 1<br>
1 1 1<br>
>> sc.tessellation<br>
ans =<br>
1 2 4 8<br>
1 3 4 8<br>
1 2 6 8<br>
1 5 6 8<br>
1 5 7 8<br>
1 3 7 8<br>
<br>
So, given a random point in the unit 3cube,<br>
<br>
>> P = rand(1,3)<br>
P =<br>
0.22843 0.652 0.06616<br>
<br>
It has sort tags of<br>
<br>
>> [psort,ptags] = sort(P)<br>
psort =<br>
0.06616 0.22843 0.652<br>
ptags =<br>
3 1 2<br>
<br>
The simplex in that dissection with the same sort order<br>
of the vertices is the second one.<br>
<br>
>> [~,stags] = sortrows(sc.domain(sc.tessellation(2,:),:)')<br>
stags =<br>
3<br>
1<br>
2<br>
<br>
>> sc.domain(sc.tessellation(2,:),:)<br>
ans =<br>
0 0 0<br>
0 1 0<br>
1 1 0<br>
1 1 1<br>
<br>
I'm not sure that this makes your problem any easier, but it<br>
may help. Personally, I like the idea that all of these simplexes<br>
in the diagonal extraction are simple transformations of each<br>
other. So all you need to do is sort the coordinates, effectively<br>
a transformation of variables.<br>
<br>
>> sort(P)<br>
ans =<br>
0.06616 0.22843 0.652<br>
<br>
This point ALWAYS falls inside the first simplex on the list...<br>
<br>
>> sc.domain(sc.tessellation(1,:),:)<br>
ans =<br>
0 0 0<br>
1 0 0<br>
1 1 0<br>
1 1 1<br>
<br>
As you can see, that simplex has a rather simple construction.<br>
<br>
HTH,<br>
John

Mon, 21 Mar 2011 13:17:23 +0000
Re: Indices of points in the nsimplex
http://www.mathworks.com/matlabcentral/newsreader/view_thread/304795#826355
Greg von Winckel
I was thinking of using the ncube [0,1]^n and the nsimplex which is determined by the ordering <br>
<br>
x1>x2>...>xn<br>
<br>
If some other simplex leads to an easier relationship for the indices, that's fine as well. The point is that I have a symmetry in my coordinates that I need to exploit.

Mon, 21 Mar 2011 15:26:05 +0000
Re: Indices of points in the nsimplex
http://www.mathworks.com/matlabcentral/newsreader/view_thread/304795#826407
Greg von Winckel
For anyone who is curious, here is my very crude solution:<br>
<br>
function dex=simplex_indices(p,n)<br>
<br>
Np=nchoosek(p,n);<br>
map = ones(Np,n);<br>
<br>
for k=2:Np<br>
<br>
d=n;<br>
map(k,:)=map(k1,:);<br>
map(k,d)=map(k,d)+1; <br>
<br>
while (sum(map(k,:))==p+1) && (d>1)<br>
<br>
map(k,d1)=map(k1,d1)+1;<br>
map(k,d:n)=1;<br>
<br>
d=d1;<br>
end<br>
<br>
end<br>
<br>
shift=repmat([ones(1,n1) 0],Np,1);<br>
dex=(mapshift)*(p.^(n1:1:0)');

Tue, 22 Mar 2011 19:52:05 +0000
Re: Indices of points in the nsimplex
http://www.mathworks.com/matlabcentral/newsreader/view_thread/304795#826709
Roger Stafford
"Greg von Winckel" wrote in message <im7qmd$c00$1@fred.mathworks.com>...<br>
> For anyone who is curious, here is my very crude solution:<br>
> <br>
> function dex=simplex_indices(p,n)<br>
> <br>
> Np=nchoosek(p,n);<br>
> map = ones(Np,n);<br>
> <br>
> for k=2:Np<br>
> <br>
> d=n;<br>
> map(k,:)=map(k1,:);<br>
> map(k,d)=map(k,d)+1; <br>
> <br>
> while (sum(map(k,:))==p+1) && (d>1)<br>
> <br>
> map(k,d1)=map(k1,d1)+1;<br>
> map(k,d:n)=1;<br>
> <br>
> d=d1;<br>
> end<br>
> <br>
> end<br>
> <br>
> shift=repmat([ones(1,n1) 0],Np,1);<br>
> dex=(mapshift)*(p.^(n1:1:0)');<br>
         <br>
The 'map' array in your 'simplex_indices' function can also be computed without forloops this way:<br>
<br>
map = sort(nchoosek(1:p,n),2);<br>
map = map  repmat(0:n1,size(map,1),1);<br>
map(:,2:n) = diff(map,1,2)+1;<br>
<br>
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.<br>
<br>
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.<br>
<br>
Roger Stafford