6

I need to generate regularly spaced samples (points) from the surface of unit simplexes with 2 or greater vertices or end points.

I can generate random samples pretty straightforwardly:

sampleSimplex[sampleSize_, vertices_] := 
 Module[{s}, 
  s = RandomReal[ExponentialDistribution[2], {sampleSize, vertices}];
  Map[Normalize[#, Total] &, s]]

Graphics3D[Point /@ sampleSimplex[1000, 3], ImageSize -> 150, 
 BoxStyle -> Directive[LightGray], ViewPoint -> {Pi, Pi/2, 2}]

simplex samples

This works nicely and enables me to specify the number of samples and the number of vertices, but as I stated above I want to generate regularly spaced points, kind of like a Mesh on a 3D plot (but not limited to 3 dimensions).

I can also do it for 2 vertices pretty simply:

sampleSimplexWith2Vertices[sampleSize_] := Module[{temp},
  temp = N[Range[0, 1, 1/sampleSize]];
  Transpose[{temp, Reverse[temp]}]]

sampleSimplexWith2Vertices[10]

{{0., 1.}, {0.1, 0.9}, {0.2, 0.8}, {0.3, 0.7}, {0.4, 0.6}, {0.5, 
  0.5}, {0.6, 0.4}, {0.7, 0.3}, {0.8, 0.2}, {0.9, 0.1}, {1., 0.}}

ListPlot[%, ImageSize -> 200]

enter image description here

but I seem blocked in thinking about a solution that extends further.

Any suggestions for how to approach this problem appreciated.

I recognize that this may appear too narrow in scope or that it seems more like a conceptual problem of how to think about the problem more than how to implement an idea in Mathematica. I can withdrawal the question, if either of those seems like the consensus.

I hoped that the approaches I've tried so far and those the question might draw would have sufficient general interest to merit some responses.


In the interest of improving the question I wanted to incorporate @Rahul_Narain's idea from his comment:

sampleSimplexWith3Vertices[pointsPerSide_] := N@Select[
   Flatten[
    Table[{x, y, 1 - x - y}, {x, 0, 1, 1/pointsPerSide}, {y, 0, 1, 
      1/pointsPerSide}], 1], #[[3]] >= 0 &]

sampleSimplexWith3Vertices[10];

enter image description here

This works for 3 vertices, where one specifies the number of points per side. But, how does one generalize this for an arbitrary number of vertices (dimensions)? I've thought that Table would give one way to go about this, but wouldn't generalizing this with Table need a way to change just about everything between its brackets?

István Zachar
  • 47,032
  • 20
  • 143
  • 291
Jagra
  • 14,343
  • 1
  • 39
  • 81
  • 1
    I think this will generalize readily to higher dimensions: p = Select[Flatten[Table[{x, y, 1 - x - y}, {x, 0, 1, 1/n}, {y, 0, 1, 1/n}], 1], #[[3]] >= 0 &] I am not adept enough with Mathematica to write a general-purpose function for arbitrary dimensions, so I will leave that for someone else to answer. –  Jun 27 '13 at 15:32

1 Answers1

5

First, define some functions:

  • simpleVertices returns the simplex corners for 1, 2, 3 dimensions;
  • toCartesian converts simplex coordinates to Cartesian coordinates;
  • simplexPoints generates regularly placed points according to dimension and resolution (default 10)

Note that regular points for $n$ dimensions are generated by injecting $n$ dependent iterators (var) into a Table:

simplexVertices[dim_Integer] := Switch[dim,
   1, {{0, 0}, {1, 0}},
   2, {{0, 0}, {1, 0}, {1/2, Sqrt@3/2}},
   3, {{0, 0, 0}, {1, 0, 0}, {1/2, Sqrt@3/2, 0}, {1/2, Sqrt@3/6, Sqrt@6/3}}];

toCartesian[pts : {__List}, dim_Integer] := toCartesian[#, dim] & /@ pts;
toCartesian[pt_List, dim_Integer] := Total[pt*simplexVertices@dim];

simplexPoints[dim_Integer, res_Integer: 10] := Module[{var, iter},
   var = Unique@"i" & /@ Range@dim;
   iter = {var[[#]], 0, 1 - Total@Take[var, # - 1], 1/res} & /@ Range@dim;
   Flatten[Table @@ Fold[Append, {Append[var, 1 - Total@var]}, iter], dim - 1]];

Now evaluate for 1, 2, 3 dimensions:

{
 Graphics[Point@toCartesian[simplexPoints@1, 1]],
 Graphics[Point@toCartesian[simplexPoints@2, 2]],
 Graphics3D[Point@toCartesian[simplexPoints@3, 3]]
 }

enter image description here

István Zachar
  • 47,032
  • 20
  • 143
  • 291
  • +1 for a clearer way to think about this. The limitation in this approach as presented goes to one needing to add to simpleVertices when I'll need to do this for higher dimensions (we may need this to go to many dimensions). Still thanks for the start. – Jagra Jun 27 '13 at 16:57
  • 1
    @Jagra Well, I have serious problems with conceptualizing higher dimensions, so please excuse me for omitting the tesseract :) – István Zachar Jun 27 '13 at 17:04
  • Ha! I have every confidence that you can conceptualize hectaeracts, septaeracts, octoeracts, all the way to ...neracts, if the ideas caught your interest for a few minutes;-) (I actually have no idea if octoeract or any of the others are actually words). – Jagra Jun 27 '13 at 17:14
  • 1
    @Jagra: The simplexPoints function works for any number of dimensions, I think, giving you points in $(n+1)$-dimensional space, just like the example code in your question. If that's the format you need, then you don't need to extend simplexVertices or toCartesian for higher dimensions at all. –  Jun 27 '13 at 17:17
  • @RahulNarain -- It does indeed, thanks for bringing my attention to it. I hadn't had time to go through all the code yet. Istvan's solution looks like a winner (but don't you love it when you get lots of interesting solutions ;-) – Jagra Jun 27 '13 at 18:28
  • It has dawned slowly on me through the day how interesting a solution you've posted. I've seen some discussion on the site about injector patterns, but this example makes me begin to realize how it could solve lots of problems. Many, many thanks. – Jagra Jun 28 '13 at 03:28
  • @Jagra It is not really an injector pattern in the control-your-evaluation-manually-sense, but a way to create dependent iterators for Table. I'm really happy that you find it useful :) – István Zachar Jun 28 '13 at 11:28
  • This is a neat answer! As a late addition, here is how to generate the vertices of a dim-dimensional simplex, centered at the origin, and with unit-length edges: simplexVertices[dim_Integer?Positive] := Transpose[Rest[Orthogonalize[Join[#, NullSpace[#]]]]/Sqrt[2]] &[{ConstantArray[1, dim + 1]}]. (Adapted from this answer.) – J. M.'s missing motivation Apr 08 '17 at 14:39