Generating all integer lattice points (Python)

From LiteratePrograms
Jump to: navigation, search

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

[edit] Generating Zn from Nn

It is sufficient to initially generate all points with nonnegative coordinates (\mathbb{N}^n), since the rest (\mathbb{Z}^n) 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.

[edit] In order of Manhattan distance

[edit] 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 (\mathbb{N}^2), 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:

Manhattan lattice.gif

[edit] In n dimensions

To be done.

[edit] In order of Euclidean distance

[edit] In two dimensions

The Euclidean distance to the origin of a point (x, y) is defined as d = \sqrt{x^2 + y^2}. 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 (\mathbb{N}^2). 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 \mathbb{Z}^2 from \mathbb{N}^2 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:

Euclidean lattice.gif

[edit] In n dimensions

The general case can be handled similarly to the two-dimensional case. In n dimensions, the squared Euclidean distance becomes x_1^2+x_2^2+\ldots+x_n^2, 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:

Euclidean 3D lattice.gif

Download code
hijacker
hijacker
hijacker
hijacker