Generating all integer lattice points (Python)
From LiteratePrograms
We will describe how to generate all integer lattice points; that is, all tuples (x1, x2, ..., xn) where the variables run through all integer values {..., -3, -2, -1, 0, 1, 2, 3, ...}. There is obviously more than one possible way to order the sequence of such tuples, but the most sensible choice is to choose an origin such as (0, 0, ..., 0) and then proceed outwards, generating points with monotonically increasing distance from the origin. Even with this restriction, there are multiple approaches depending on the metric used to calculate distance.
Contents |
Generating Zn from Nn
It is sufficient to initially generate all points with nonnegative coordinates (), since the rest () can be obtained by reflecting each point across all combinations of the coordinate axes. In the case of n = 2, there are three such reflections for each point and the routine looks as follows:
<<generate Z2>>= def generate_Z2(limit=None, origin=(0, 0)): for p in generate_N2(limit): p = (p[0]+origin[0], p[1]+origin[1]) yield p if p[0] != 0: yield -p[0], p[1] if p[1] != 0: yield p[0], -p[1] if p[0] and p[1]: yield -p[0], -p[1]
The function generate_N2 called here is assumed to generate points up to a distance of at most limit units from the origin (continuing forever if set to None). A custom origin other than (0, 0) may be provided. Note that points lying on the x and y axes should not be reflected since they would be reflected onto themselves and thus yielded twice; hence the checks for nonzero coordinates.
In order of Manhattan distance
In two dimensions
The Manhattan distance is the simplest metric along a perpendicular grid, for two points (a, b) and (c, d) defined as |a-c| + |b-d|. In other words, the Manhattan distance is the minimum distance required to walk from a to b if one can only move strictly horizontally or vertically.
For points (x, y) located in the first quadrant (), the distance x+y is simply the sum x+y. (The geometric interpretation of this formula is that points with the same distance lie on a 45° diagonal whose side faces the origin.) Thus, to generate all points in the first quadrant, we may loop through all distances d = 0, 1, 2, ... and for each d find all pairs x and y that sum to d. An implementation in Python might look like this:
<<generate N2, Manhattan distance>>= def generate_N2(limit=None): d = 0 while limit is None or d <= limit: for x in range(d+1): yield x, d-x d += 1
We create a module for this generator:
<<manhattan2d.py>>= generate N2, Manhattan distance generate Z2
Here is an animation of the output from this generator:
In n dimensions
To be done.
In order of Euclidean distance
In two dimensions
The Euclidean distance to the origin of a point (x, y) is defined as . With this metric, the points that share a common distance from the origin are those that lie on the circle with radius d. If we generate points by increasing Euclidean distance, they will form the interior of an expanding circle.
We define a helper function that computes the square of the Euclidean distance between a point and the origin. (Taking the square root of this expression is unnecessary, and if we skip it, we only have to work with integer arithmetic.)
<<distance calculation>>= def distance(p): return p[0]**2 + p[1]**2
The main generator will yield all points in the first quadrant (). We will treat this quadrant as a list of vertical columns, one for each x coordinate. We now create a list whose n-th entry holds the maximum y-coordinate of any point with that x-coordinate yielded so far; that is, each entry will represent how far up the corresponding column has been "filled". For each distance d = 0, 1, 3, ..., we proceed by finding all points whose distance from the origin is in the interval [d-1, d] and store them in a list yieldable. Then we sort the list of acquired points by distance and yield them in order.
<<generate N2, Euclidean distance>>= def generate_N2(limit=None): ymax = [0] d = 0 while limit is None or d <= limit: yieldable = [] find all remaining points closer than d sort list by distance for p in yieldable: yield p d += 1 ymax.append(0) # Extend to make room for column[d]
To find the intended points, we create an infinite loop that repeatedly fetches points whose distance is in the sought interval. The loop ends when there are no more points to yield.
<<find all remaining points closer than d>>= while 1: batch = [] for x in range(d+1): y = ymax[x] if distance((x, y)) <= d**2: # Note: distance squared batch.append((x, y)) ymax[x] += 1 if not batch: break yieldable += batch
Sorting by distance can be achieved by providing a custom key function:
<<sort list by distance>>= yieldable.sort(key=distance)
To create a complete package, we re-use the routine that generates from implemented earlier.
<<euclidean2d.py>>= distance calculation generate N2, Euclidean distance generate Z2
The functions can be tested by counting the number of points they produce for given distances:
[len(list(generate_N2(n))) for n in range(20)] [len(list(generate_Z2(n))) for n in range(20)]
The respective outputs should be 1, 3, 6, 11, 17, 26, 35, 45, 58, ..., and 1, 5, 13, 29, 49, 81, 113, 149, 197, ...; corresponding to A000603 ("Number of nonnegative solutions to x² + y² ≤ n") and A000328 ("Number of points of norm ≤ n² in square lattice.") in the On-Line Encyclopedia of Integer Sequences.
The animated output is the following:
In n dimensions
The general case can be handled similarly to the two-dimensional case. In n dimensions, the squared Euclidean distance becomes , in Python:
<<squared Euclidean distance, n dimensions>>= def distance(p): return sum(x**2 for x in p)
Instead of slicing a two-dimensional plane into a one-dimensional list of "heights", we will now split the n-dimensional space into an (n-1)-dimensional array ("x coordinates") of heights ("y coordinates"). Python does not have an efficient multidimensional array type built-in (see SciPy), but we can easily simulate a multidimensional array by indexing a dict with tuples. Values that have not been set should default to zero, for which purpose we can use Python 2.5's defaultdict class:
<<custom dictionary class>>= from collections import defaultdict
We also need a way to generate all array indices up to a bound m. This can be implemented by counting to mn and writing the number as a base-m integer with n digits:
<<multirange>>= def multirange(dims, m): c = 0 for i in xrange(m**dims): x = i r = [] for n in range(dims): x, rem = divmod(x, m) r.append(rem) yield tuple(r)[::-1]
The main function becomes:
<<generate Nn, Euclidean distance>>= def generate_Nn(n=2, limit=None): custom dictionary class ymax = defaultdict(int) # unknown keys are mapped to 0 d = 0 while limit is None or d <= limit: yieldable = [] while 1: batch = [] for x in multirange((n-1), d+1): y = ymax[x] pt = x + (y,) if distance(pt) <= d**2: batch.append(pt) ymax[x] += 1 if not batch: break yieldable += batch yieldable.sort(key=distance) for p in yieldable: yield p d += 1
The only essential differences from the two-dimensional case are the iteration over an (n-1)-dimensional multirange, and the use of tuple concatenation (x + (y,)).
To test, we can verify that all generated points are correctly ordered, and check the number of generated points against known data. The sequences for n = 3, 4, 5, 6, and 7 are A000604, A055403, A055404, A055405, and A055406.
<<test generate Nn, Euclidean distance>>= testdata = { \ 0 : [1, 1, 1, 1, 1, 1, 1, 1], 1 : [1, 2, 3, 4, 5, 6, 7, 8], 2 : [1, 3, 6, 11, 17, 26, 35, 45], 3 : [1, 4, 11, 29, 54, 99, 163, 239], 4 : [1, 5, 20, 70, 165, 357, 688, 1154], 5 : [1, 6, 36, 157, 482, 1203, 2673, 5139], 6 : [1, 7, 63, 337, 1319, 3819, 9763, 21374], 7 : [1, 8, 106, 702, 3390, 11496, 33792, 83877] } if __name__ == "__main__": for dim in range(8): for lim in range(8): a = list(generate_Nn(n=dim, limit=lim)) for i in range(1, len(a)): assert distance(a[i-1]) <= distance(a[i]) assert testdata[dim][lim] == len(a)
<<euclidean.py>>= squared Euclidean distance, n dimensions multirange generate Nn, Euclidean distance test generate Nn, Euclidean distance
Although the results are difficult to visualize in higher dimensions, the three-dimensional case is simple. The points generated by generate_Nn with n=3 make up one-eighth of a growing ball:
| Download code |