# Python generator to produce lattice points on an n-simplex.

Looking for a code review, and hopefully to learn something if someone has a nicer solution. Here's what I wrote:

from __future__ import division, print_function
from future_builtins import *
import numpy as np

def _walk(num_dims, samples_per_dim, max_):
if num_dims == 0:
yield np.array([(max_ - 1) / (samples_per_dim - 1)])
else:
for i in range(max_):
for rest in _walk(num_dims - 1, samples_per_dim, max_ - i):
yield np.concatenate((np.array([i]) / (samples_per_dim - 1),
rest))

def walk(num_dims, samples_per_dim):
"""
A generator that returns lattice points on an n-simplex.
"""
return _walk(num_dims, samples_per_dim, samples_per_dim)


So, list(walk(2, 3)) yields:

[array([ 0.,  0.,  1.]),
array([ 0. ,  0.5,  0.5]),
array([ 0.,  1.,  0.]),
array([ 0.5,  0. ,  0.5]),
array([ 0.5,  0.5,  0. ]),
array([ 1.,  0.,  0.])]


Concatenating numpy arrays isn't a really good idea because that's not how they were designed to be used.

A better way might be this:

def walk(num_dims, samples_per_dim):
"""
A generator that returns lattice points on an n-simplex.
"""
values = np.arange(samples_per_dim) / (samples_per_dim - 1)
for items in itertools.product(values, repeat = num_dims+1):
if sum(items) == 1.0:
yield items


Basically, we iterate over all possible combinations of those points, and filter out the invalid ones. This may seem wasteful, but since itertools is written in C, its probably actually faster then your solution. It produces tuples rather then numpy arrays.

• Thanks. Why not use linspace instead of arange? Jun 8 '11 at 2:15
• Also, have you tested this code? Why don't floating point precision errors cause problems with the comparison to 1.0? Jun 8 '11 at 2:16
• @Niel, I was not aware of linspace. That would be better. I did test it. I'm not sure why floating point isn't a problem... Jun 8 '11 at 3:54

Using Winston Ewert's suggestions of using itertools, and using lists instead of numpy arrays internally, here's an alternate solution:

def walk(num_dims, samples_per_dim):
"""
A generator that returns lattice points on an n-simplex.
"""
max_ = samples_per_dim + num_dims - 1
for c in combinations(range(max_), num_dims):
c = list(c)
yield [(y - x - 1) / (samples_per_dim - 1)
for x, y in izip([-1] + c, c + [max_])]