simple demo low disrepancy sequence -- Halton sequence in 2d using bases (2, 3)

lds.py

from itertools import izip

def decompose(b, n):
    k = 0
    while n:
        r = n % b
        q = (n-r)/b
        yield (r, k)
        k += 1
        n = q

def recompose(b, decomp):
    acc = 0
    for r, k in decomp:
        acc += r * (b**k)
    return acc

def recompose_vdc(b, decomp):
    acc = 0
    for r, k in decomp:
        acc += r * (b**-(k+1))
    return acc

def main():
    test_integers = range(100)
    bases = (2, 3, 5)

    for b in bases:
        for x in test_integers:
            decomp = list(decompose(b, x))
            y = recompose(b, decomp)
            y_vdc = recompose_vdc(b, decomp)
            print (b, x, y, y_vdc)

def gen_integers():
    x = 0
    while True:
        yield x
        x += 1

def gen_vdc(b):
    for x in gen_integers():
        decomp = decompose(b, x)
        y_vdc = recompose_vdc(b, decomp)
        yield y_vdc


def gen_halton(coprimes):
    generators = map(gen_vdc, coprimes)
    for coords in izip(*generators):
        yield coords


def plot_some_2d_halton_sequences():
    coprimes = (2, 3)

    n_points = [256, 256, 1024, 1024]
    method = ['uniform_random', 'halton', 'uniform_random', 'halton']

    import pylab
    import numpy
    pylab.figure()

    for i in range(4):
        pylab.subplot(2, 2, i+1)
        points = []

        n = n_points[i]
        if method[i] == 'uniform_random':
            xs = numpy.random.uniform(0.0, 1.0, n)
            ys = numpy.random.uniform(0.0, 1.0, n)
        else:
            # https://en.wikipedia.org/wiki/Low-discrepancy_sequence#Halton_sequence
            for j, point in enumerate(gen_halton(coprimes)):
                if j >= n:
                    break
                points.append(point)
            xs, ys = zip(*points)
        title = '%s; %d points' % (method[i], n)
        pylab.title(title)
        pylab.scatter(xs, ys, marker='o', color='b')
    pylab.show()


if __name__ == '__main__':
    # main()
    plot_some_2d_halton_sequences()