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:
Indices of points in the n-simplex

Subject: Indices of points in the n-simplex

From: Greg von Winckel

Date: 21 Mar, 2011 10:32:04

Message: 1 of 6

Here is another tricky question.

Suppose I have an n-dimensional 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 n-simplexes in the corner of the hypercube?

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

A=arrayfun(@(j) j*n+(2+j:n)',0:n-1,'UniformOutput',false)
cat(1,A{:})

In 3D, I don't have a clever formula, but I can do this horribly inefficiently with
[x,y,z]=ndgrid(1:n)
find(y(:)<x(:) & z(:)<y(:))

but at least it demonstrates the principle.

Is there already a formula for determining these points?

Subject: Indices of points in the n-simplex

From: Bruno Luong

Date: 21 Mar, 2011 10:44:04

Message: 2 of 6

"Greg von Winckel" wrote in message <im79f4$50r$1@fred.mathworks.com>...
> Here is another tricky question.
>
> Suppose I have an n-dimensional 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 n-simplexes in the corner of the hypercube?
>
> 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
>
> A=arrayfun(@(j) j*n+(2+j:n)',0:n-1,'UniformOutput',false)
> cat(1,A{:})
>
> In 3D, I don't have a clever formula, but I can do this horribly inefficiently with
> [x,y,z]=ndgrid(1:n)
> find(y(:)<x(:) & z(:)<y(:))
>
> but at least it demonstrates the principle.
>
> Is there already a formula for determining these points?

It seems somewhat related to this thread.
http://www.mathworks.com/matlabcentral/newsreader/view_thread/304546#825843

Bruno

Subject: Indices of points in the n-simplex

From: John D'Errico

Date: 21 Mar, 2011 11:39:04

Message: 3 of 6

"Greg von Winckel" wrote in message <im79f4$50r$1@fred.mathworks.com>...
> Here is another tricky question.
>
> Suppose I have an n-dimensional 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 n-simplexes in the corner of the hypercube?
>
> 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
>
> A=arrayfun(@(j) j*n+(2+j:n)',0:n-1,'UniformOutput',false)
> cat(1,A{:})
>
> In 3D, I don't have a clever formula, but I can do this horribly inefficiently with
> [x,y,z]=ndgrid(1:n)
> find(y(:)<x(:) & z(:)<y(:))
>
> but at least it demonstrates the principle.
>
> Is there already a formula for determining these points?

To make it more complex, in n-dimensions, a hypercube can
be dissected into simplexes in various ways. As it turns out,
there is only one way to dissect a square, ignoring rotations.
But a 3-cube can be dissected into tetrahedra in either 5 or
6 ways. (Try it! Start by cutting away the corners of the cube.)

Likewise, a 4-cube can be dissected into 24 or 16 simplexes,
using the two simple schemes that I know of, but there are
probably other solutions. And as for the 5-cube, it can be
dissected into either 68 or 120 simplexes, but I am quite sure
there are other solutions there beyond those I know of. For
the 6-cube, the two solutions I know how to construct result
in either 376 or 720 simplexes. I won't go any higher since
my head is starting to hurt.

The point of all this is you need to define the dissection of
the n-cube you will be using. The diagonal extraction, using
factorial(n) simplexes to do the task is my favorite, as that
dissection devolves into a sort. Essentially all of those
simplexes are identical in shape and size to the rest, and can
be transformed into the others by simple transpositions of
the coordinate axes.

To get back to your question, the diagonal extraction is nice
because it has a direct relation to the sorting of your numbers.
Dissect the 3-cube into 6 tetrahedra. Here is the dissection, as
done by my simplicial complex toolbox:

>> sc = tessellatehypercube(3)
sc =
          domain: [8x3 double]
    tessellation: [6x4 double]
           range: [8x0 double]
>> sc.domain
ans =
     0 0 0
     1 0 0
     0 1 0
     1 1 0
     0 0 1
     1 0 1
     0 1 1
     1 1 1
>> sc.tessellation
ans =
     1 2 4 8
     1 3 4 8
     1 2 6 8
     1 5 6 8
     1 5 7 8
     1 3 7 8

So, given a random point in the unit 3-cube,

>> P = rand(1,3)
P =
      0.22843 0.652 0.06616

It has sort tags of

>> [psort,ptags] = sort(P)
psort =
      0.06616 0.22843 0.652
ptags =
     3 1 2

The simplex in that dissection with the same sort order
of the vertices is the second one.

>> [~,stags] = sortrows(sc.domain(sc.tessellation(2,:),:)')
stags =
     3
     1
     2

>> sc.domain(sc.tessellation(2,:),:)
ans =
     0 0 0
     0 1 0
     1 1 0
     1 1 1

I'm not sure that this makes your problem any easier, but it
may help. Personally, I like the idea that all of these simplexes
in the diagonal extraction are simple transformations of each
other. So all you need to do is sort the coordinates, effectively
a transformation of variables.

>> sort(P)
ans =
      0.06616 0.22843 0.652

This point ALWAYS falls inside the first simplex on the list...

>> sc.domain(sc.tessellation(1,:),:)
ans =
     0 0 0
     1 0 0
     1 1 0
     1 1 1

As you can see, that simplex has a rather simple construction.

HTH,
John

Subject: Indices of points in the n-simplex

From: Greg von Winckel

Date: 21 Mar, 2011 13:17:23

Message: 4 of 6

I was thinking of using the n-cube [0,1]^n and the n-simplex which is determined by the ordering

x1>x2>...>xn

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.

Subject: Indices of points in the n-simplex

From: Greg von Winckel

Date: 21 Mar, 2011 15:26:05

Message: 5 of 6

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)');

Subject: Indices of points in the n-simplex

From: Roger Stafford

Date: 22 Mar, 2011 19:52:05

Message: 6 of 6

"Greg von Winckel" wrote in message <im7qmd$c00$1@fred.mathworks.com>...
> 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

Tags for this Thread

No tags are associated with 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