#!/usr/bin/env python
""" Walker's alias method for random objects with different probablities
walkerrandom.py
Examples
--------
# 0 1 2 or 3 with probabilities .1 .2 .3 .4 --
wrand = Walkerrandom( [10, 20, 30 40] ) # builds the Walker tables
wrand.random() # each call -> 0 1 2 or 3
# for example, 1000 calls with random.seed(1) -> [96, 199, 334, 371]
# strings A B C or D with probabilities .1 .2 .3 .4 --
abcd = dict( A=1, D=4, C=3, B=2 )
# keys can be any immutables: 2d points, colors, atoms ...
wrand = Walkerrandom( abcd.values(), abcd.keys() )
wrand.random() # each call -> "A" "B" "C" or "D"
# fast: 1 randint(), 1 uniform(), table lookup
How it works
------------
For weights 10 20 30 40 as above, picture sticks A B C D of those lengths:
10 AAAAAAAAAA
20 BBBBBBBBBB BBBBBBBBBB
30 CCCCCCCCCC CCCCCCCCCC CCCCCCCCCC
40 DDDDDDDDDD DDDDDDDDDD DDDDDDDDDD DDDDDDDDDD
Split and rearrange them into equal-length rows, like this:
AAAAAAAAAA DDDDDDDDDDDDDDD -- 10 A + 15 D = 40% A + 60% D
BBBBBBBBBBBBBBBBBBBB DDDDD -- 20 B + 5 D = 80% B + 20% D
CCCCCCCCCCCCCCCCCCCCCCCCC -- 25 C = 100% C
DDDDDDDDDDDDDDDDDDDD CCCCC -- 20 D + 5 C = 80% D + 20% C
Clearly 10 % of the area is A, 20 % B, 30 % C and 40 % D --
we haven't changed areas, just rearranged.
Now to choose a random one of A or B or C or D,
throw a dart at a "dart board" of the sticks in these 4 rows:
if it hits row 0, return A with probablity 40 % / D 60 %
if it hits row 1, return B with probablity 80 % / D 20 %
...
This picture is in Devroye, p. 111 (rediscovered here).
Walker's algorithm essentially arranges a given lot of sticks
into equal-length rows: pick a row shorter than average
and a row longer than average, split the longer to fill the shorter,
iterate until they're all the same length.
Notes
To generate random colors similar to those in a given picture,
first collect color samples in a histogram:
for color in ...:
# cluster e.g. rrggbb -> rgb, 16^3 bins
# (many many methods, see Wikipedia Data_clustering)
colors[color] += 1
(cPickle to a file, write it, read it back in)
then use Walkerrandom to select colors with these frequencies:
colorrand = Walkerrandom( colors.values(), colors.keys() )
colorrand.random() # each call -> a color
References
L. Devroye, Non-Uniform Random Variate Generation, 1986, p. 107 ff.
http://cg.scs.carleton.ca/~luc/rnbookindex.html (800 pages)
Knuth, Stanford GraphBase, 1993, p. 392
C++ hat random container by AngleWyrm,
http://home.comcast.net/~anglewyrm/hat.html
"""
from __future__ import division
import random
__author__ = "Denis Bzowy"
__version__ = "16nov2008"
Test = 0
#...............................................................................
class Walkerrandom:
""" Walker's alias method for random objects with different probablities
"""
def __init__( self, weights, keys=None ):
""" builds the Walker tables prob and inx for calls to random().
The weights (a list or tuple or iterable) can be in any order;
they need not sum to 1.
"""
n = self.n = len(weights)
self.keys = keys
sumw = sum(weights)
prob = [w * n / sumw for w in weights] # av 1
inx = [-1] * n
short = [j for j, p in enumerate( prob ) if p < 1]
long = [j for j, p in enumerate( prob ) if p > 1]
while short and long:
j = short.pop()
k = long[-1]
# assert prob[j] <= 1 <= prob[k]
inx[j] = k
prob[k] -= (1 - prob[j]) # -= residual weight
if prob[k] < 1:
short.append( k )
long.pop()
if Test:
print "test Walkerrandom: j k pk: %d %d %.2g" % (j, k, prob[k])
self.prob = prob
self.inx = inx
if Test:
print "test", self
def __str__( self ):
""" e.g. "Walkerrandom prob: 0.4 0.8 1 0.8 inx: 3 3 -1 2" """
probstr = " ".join([ "%.2g" % x for x in self.prob ])
inxstr = " ".join([ "%.2g" % x for x in self.inx ])
return "Walkerrandom prob: %s inx: %s" % (probstr, inxstr)
#...............................................................................
def random( self ):
""" each call -> a random int or key with the given probability
fast: 1 randint(), 1 random.uniform(), table lookup
"""
u = random.uniform( 0, 1 )
j = random.randint( 0, self.n - 1 ) # or low bits of u
randint = j if u <= self.prob[j] \
else self.inx[j]
return self.keys[randint] if self.keys \
else randint
#...............................................................................
if __name__ == "__main__":
# little examples, self-contained --
N = 5
Nrand = 1000
randomseed = 1
try:
import bz.util
bz.util.scan_eq_args( globals(), __doc__ ) # N=5 ...
except ImportError:
pass
if randomseed:
random.seed( randomseed )
print Nrand, "Walkerrandom with weights .1 .2 .3 .4:"
w = range( 1, N )
wrand = Walkerrandom( w )
nrand = [0] * (N - 1)
for _ in range( Nrand ):
j = wrand.random()
nrand[j] += 1
s = str( nrand )
print s
if N==5 and Nrand==1000 and randomseed==1:
assert s == "[96, 199, 334, 371]"
print Nrand, "Walkerrandom strings with weights .1 .2 .3 .4:"
abcd = dict( A=1, D=4, C=3, B=2 )
wrand = Walkerrandom( abcd.values(), abcd.keys() )
from collections import defaultdict
nrand = defaultdict(int) # init 0
for _ in range( Nrand ):
j = wrand.random()
nrand[j] += 1
s = str( sorted( nrand.iteritems() ))
print s
if N==5 and Nrand==1000 and randomseed==1:
assert s == "[('A', 105), ('B', 181), ('C', 283), ('D', 431)]"
# end walkerrandom.py
