Skip to content

Instantly share code, notes, and snippets.

@fcostin
Created October 16, 2017 08:20
Show Gist options
  • Save fcostin/d6bf11a7c6f7a3e2c76a2e1132ebe0e4 to your computer and use it in GitHub Desktop.
Save fcostin/d6bf11a7c6f7a3e2c76a2e1132ebe0e4 to your computer and use it in GitHub Desktop.
simple demo low disrepancy sequence -- Halton sequence in 2d using bases (2, 3)
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment