[Python-checkins] python/nondist/sandbox/twister random.py,1.2,1.3
rhettinger@users.sourceforge.net
rhettinger@users.sourceforge.net
Wed, 18 Dec 2002 10:26:00 -0800
Update of /cvsroot/python/python/nondist/sandbox/twister
In directory sc8-pr-cvs1:/tmp/cvs-serv30214
Modified Files:
random.py
Log Message:
Altered to incorporate the MersenneTwister.
* Wichmann-Hill specific methods moved into a subclass of Random()
* Random.__init__() now automatically loads in the MersenneTwister
methods whenever self isn't a subclass.
* Added a new method _randbelow(n) to replace the frequent references to
int(random() * n). Not only is it clearer, but it allows the C code
to override with fast random integer generation. It is avoids all the
overhead of creating a random float and turning it back into an integer.
Internally, it takes half as long to generate a random int as it does
a 53-bit precision float.
* Confined the jumpahead() test to the Wichmann-Hill generator. The
semantics are different for the Mersenne Twister which doesn't have
a directly computable jumpahead.
* The seed(x) method now accepts multiple arguments, seed(x1, x2, ...).
This is a backward compatible accomodation for long period generators
which have a much larger state space and can benefit from a wider
range of seeds. An alternative approach would have use Python's long
integers; however, this would have make it more difficult to compare
the generator back to the reference implementation.
Index: random.py
===================================================================
RCS file: /cvsroot/python/python/nondist/sandbox/twister/random.py,v
retrieving revision 1.2
retrieving revision 1.3
diff -C2 -d -r1.2 -r1.3
*** random.py 18 Dec 2002 18:10:26 -0000 1.2
--- random.py 18 Dec 2002 18:25:58 -0000 1.3
***************
*** 108,111 ****
--- 108,113 ----
# Adrian Baddeley.
+ import MersenneTwister
+
class Random:
"""Random number generator base class used by bound module functions.
***************
*** 123,204 ****
"""
! VERSION = 1 # used by getstate/setstate
! def __init__(self, x=None):
"""Initialize an instance.
! Optional argument x controls seeding, as for Random.seed().
! """
!
! self.seed(x)
!
! ## -------------------- core generator -------------------
!
! # Specific to Wichmann-Hill generator. Subclasses wishing to use a
! # different core generator should override the seed(), random(),
! # getstate(), setstate() and jumpahead() methods.
!
! def seed(self, a=None):
! """Initialize internal state from hashable object.
!
! None or no argument seeds from current time.
!
! If a is not None or an int or long, hash(a) is used instead.
!
! If a is an int or long, a is used directly. Distinct values between
! 0 and 27814431486575L inclusive are guaranteed to yield distinct
! internal states (this guarantee is specific to the default
! Wichmann-Hill generator).
"""
!
! if a is None:
! # Initialize from current time
! import time
! a = long(time.time() * 256)
!
! if type(a) not in (type(3), type(3L)):
! a = hash(a)
!
! a, x = divmod(a, 30268)
! a, y = divmod(a, 30306)
! a, z = divmod(a, 30322)
! self._seed = int(x)+1, int(y)+1, int(z)+1
!
self.gauss_next = None
!
! def random(self):
! """Get the next random number in the range [0.0, 1.0)."""
!
! # Wichman-Hill random number generator.
! #
! # Wichmann, B. A. & Hill, I. D. (1982)
! # Algorithm AS 183:
! # An efficient and portable pseudo-random number generator
! # Applied Statistics 31 (1982) 188-190
! #
! # see also:
! # Correction to Algorithm AS 183
! # Applied Statistics 33 (1984) 123
! #
! # McLeod, A. I. (1985)
! # A remark on Algorithm AS 183
! # Applied Statistics 34 (1985),198-200
!
! # This part is thread-unsafe:
! # BEGIN CRITICAL SECTION
! x, y, z = self._seed
! x = (171 * x) % 30269
! y = (172 * y) % 30307
! z = (170 * z) % 30323
! self._seed = x, y, z
! # END CRITICAL SECTION
!
! # Note: on a platform using IEEE-754 double arithmetic, this can
! # never return 0.0 (asserted by Tim; proof too long for a comment).
! return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
!
def getstate(self):
"""Return internal state; can be passed to setstate() later."""
! return self.VERSION, self._seed, self.gauss_next
def setstate(self, state):
--- 125,149 ----
"""
! VERSION = 1
! def __init__(self, *seeds):
"""Initialize an instance.
! Optional argument "seeds" controls seeding, as for Random.seed().
"""
! if self.__class__ == Random:
! geninstance = MersenneTwister.Random()
! self.random = geninstance.random
! self.seed = geninstance.seed
! self._getstate = geninstance.getstate
! self._setstate = geninstance.setstate
! self._randbelow = geninstance._randbelow
! self.jumpahead = geninstance.jumpahead
self.gauss_next = None
! self.seed(*seeds)
!
def getstate(self):
"""Return internal state; can be passed to setstate() later."""
! return self.VERSION, self._getstate(), self.gauss_next
def setstate(self, state):
***************
*** 206,210 ****
version = state[0]
if version == 1:
! version, self._seed, self.gauss_next = state
else:
raise ValueError("state with version %s passed to "
--- 151,156 ----
version = state[0]
if version == 1:
! version, internalstate, self.gauss_next = state
! self._setstate(internalstate)
else:
raise ValueError("state with version %s passed to "
***************
*** 212,283 ****
(version, self.VERSION))
- def jumpahead(self, n):
- """Act as if n calls to random() were made, but quickly.
-
- n is an int, greater than or equal to 0.
-
- Example use: If you have 2 threads and know that each will
- consume no more than a million random numbers, create two Random
- objects r1 and r2, then do
- r2.setstate(r1.getstate())
- r2.jumpahead(1000000)
- Then r1 and r2 will use guaranteed-disjoint segments of the full
- period.
- """
-
- if not n >= 0:
- raise ValueError("n must be >= 0")
- x, y, z = self._seed
- x = int(x * pow(171, n, 30269)) % 30269
- y = int(y * pow(172, n, 30307)) % 30307
- z = int(z * pow(170, n, 30323)) % 30323
- self._seed = x, y, z
-
- def __whseed(self, x=0, y=0, z=0):
- """Set the Wichmann-Hill seed from (x, y, z).
-
- These must be integers in the range [0, 256).
- """
-
- if not type(x) == type(y) == type(z) == int:
- raise TypeError('seeds must be integers')
- if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
- raise ValueError('seeds must be in range(0, 256)')
- if 0 == x == y == z:
- # Initialize from current time
- import time
- t = long(time.time() * 256)
- t = int((t&0xffffff) ^ (t>>24))
- t, x = divmod(t, 256)
- t, y = divmod(t, 256)
- t, z = divmod(t, 256)
- # Zero is a poor seed, so substitute 1
- self._seed = (x or 1, y or 1, z or 1)
-
- self.gauss_next = None
-
- def whseed(self, a=None):
- """Seed from hashable object's hash code.
-
- None or no argument seeds from current time. It is not guaranteed
- that objects with distinct hash codes lead to distinct internal
- states.
-
- This is obsolete, provided for compatibility with the seed routine
- used prior to Python 2.1. Use the .seed() method instead.
- """
-
- if a is None:
- self.__whseed()
- return
- a = hash(a)
- a, x = divmod(a, 256)
- a, y = divmod(a, 256)
- a, z = divmod(a, 256)
- x = (x + a) % 256 or 1
- y = (y + a) % 256 or 1
- z = (z + a) % 256 or 1
- self.__whseed(x, y, z)
-
## ---- Methods below this point do not need to be overridden when
## ---- subclassing for the purpose of using a different core generator.
--- 158,161 ----
***************
*** 293,297 ****
## -------------------- integer methods -------------------
! def randrange(self, start, stop=None, step=1, int=int, default=None):
"""Choose a random item from range(start, stop[, step]).
--- 171,178 ----
## -------------------- integer methods -------------------
! def _randbelow(self, n, int=int):
! return int(self.random() * n)
!
! def randrange(self, start, stop=None, step=1, _randbelow=_randbelow, default=None):
"""Choose a random item from range(start, stop[, step]).
***************
*** 308,312 ****
if stop is default:
if istart > 0:
! return int(self.random() * istart)
raise ValueError, "empty range for randrange()"
--- 189,193 ----
if stop is default:
if istart > 0:
! return _randbelow(istart)
raise ValueError, "empty range for randrange()"
***************
*** 317,321 ****
if step == 1 and istart < istop:
try:
! return istart + int(self.random()*(istop - istart))
except OverflowError:
# This can happen if istop-istart > sys.maxint + 1, and
--- 198,202 ----
if step == 1 and istart < istop:
try:
! return istart + _randbelow(istop - istart)
except OverflowError:
# This can happen if istop-istart > sys.maxint + 1, and
***************
*** 342,346 ****
if n <= 0:
raise ValueError, "empty range for randrange()"
! return istart + istep*int(self.random() * n)
def randint(self, a, b):
--- 223,227 ----
if n <= 0:
raise ValueError, "empty range for randrange()"
! return istart + istep*self._randbelow(n)
def randint(self, a, b):
***************
*** 356,360 ****
return seq[int(self.random() * len(seq))]
! def shuffle(self, x, random=None, int=int):
"""x, random=random.random -> shuffle list x in place; return None.
--- 237,241 ----
return seq[int(self.random() * len(seq))]
! def shuffle(self, x, _randbelow=None):
"""x, random=random.random -> shuffle list x in place; return None.
***************
*** 368,379 ****
"""
! if random is None:
! random = self.random
for i in xrange(len(x)-1, 0, -1):
# pick an element in x[:i+1] with which to exchange x[i]
! j = int(random() * (i+1))
x[i], x[j] = x[j], x[i]
! def sample(self, population, k, random=None, int=int):
"""Chooses k unique random elements from a population sequence.
--- 249,260 ----
"""
! if _randbelow is None:
! _randbelow = self._randbelow
for i in xrange(len(x)-1, 0, -1):
# pick an element in x[:i+1] with which to exchange x[i]
! j = _randbelow(i+1)
x[i], x[j] = x[j], x[i]
! def sample(self, population, k, _randbelow=None):
"""Chooses k unique random elements from a population sequence.
***************
*** 413,423 ****
if not 0 <= k <= n:
raise ValueError, "sample larger than population"
! if random is None:
! random = self.random
result = [None] * k
if n < 6 * k: # if n len list takes less space than a k len dict
pool = list(population)
for i in xrange(k): # invariant: non-selected at [0,n-i)
! j = int(random() * (n-i))
result[i] = pool[j]
pool[j] = pool[n-i-1]
--- 294,304 ----
if not 0 <= k <= n:
raise ValueError, "sample larger than population"
! if _randbelow is None:
! _randbelow = self._randbelow
result = [None] * k
if n < 6 * k: # if n len list takes less space than a k len dict
pool = list(population)
for i in xrange(k): # invariant: non-selected at [0,n-i)
! j = _randbelow(n-i)
result[i] = pool[j]
pool[j] = pool[n-i-1]
***************
*** 425,431 ****
selected = {}
for i in xrange(k):
! j = int(random() * n)
while j in selected:
! j = int(random() * n)
result[i] = selected[j] = population[j]
return result
--- 306,312 ----
selected = {}
for i in xrange(k):
! j = _randbelow(n)
while j in selected:
! j = _randbelow(n)
result[i] = selected[j] = population[j]
return result
***************
*** 745,748 ****
--- 626,780 ----
return alpha * pow(-_log(u), 1.0/beta)
+ ## -------------------- Wichmann-Hill -------------------
+
+ class WichmannHill(Random):
+ # Specific to Wichmann-Hill generator. Subclasses wishing to use a
+ # different core generator should override the seed(), random(),
+ # getstate(), setstate() and jumpahead() methods.
+
+ def seed(self, *seeds):
+ """Initialize internal state from hashable object.
+
+ None or no argument seeds from current time.
+
+ If a is not None or an int or long, hash(a) is used instead.
+
+ If a is an int or long, a is used directly. Distinct values between
+ 0 and 27814431486575L inclusive are guaranteed to yield distinct
+ internal states (this guarantee is specific to the default
+ Wichmann-Hill generator).
+ """
+
+ if len(seeds):
+ a = 0
+ for s in seeds:
+ if not isinstance(s, (int, long)):
+ s = hash(s)
+ a ^= s
+ else:
+ # Initialize from current time
+ import time
+ a = long(time.time() * 256)
+
+ a, x = divmod(a, 30268)
+ a, y = divmod(a, 30306)
+ a, z = divmod(a, 30322)
+ self._seed = int(x)+1, int(y)+1, int(z)+1
+
+ self.gauss_next = None
+
+ def random(self):
+ """Get the next random number in the range [0.0, 1.0)."""
+
+ # Wichman-Hill random number generator.
+ #
+ # Wichmann, B. A. & Hill, I. D. (1982)
+ # Algorithm AS 183:
+ # An efficient and portable pseudo-random number generator
+ # Applied Statistics 31 (1982) 188-190
+ #
+ # see also:
+ # Correction to Algorithm AS 183
+ # Applied Statistics 33 (1984) 123
+ #
+ # McLeod, A. I. (1985)
+ # A remark on Algorithm AS 183
+ # Applied Statistics 34 (1985),198-200
+
+ # This part is thread-unsafe:
+ # BEGIN CRITICAL SECTION
+ x, y, z = self._seed
+ x = (171 * x) % 30269
+ y = (172 * y) % 30307
+ z = (170 * z) % 30323
+ self._seed = x, y, z
+ # END CRITICAL SECTION
+
+ # Note: on a platform using IEEE-754 double arithmetic, this can
+ # never return 0.0 (asserted by Tim; proof too long for a comment).
+ return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
+
+ def getstate(self):
+ """Return internal state; can be passed to setstate() later."""
+ return self.VERSION, self._seed, self.gauss_next
+
+ def setstate(self, state):
+ """Restore internal state from object returned by getstate()."""
+ version = state[0]
+ if version == 1:
+ version, self._seed, self.gauss_next = state
+ else:
+ raise ValueError("state with version %s passed to "
+ "Random.setstate() of version %s" %
+ (version, self.VERSION))
+
+ def jumpahead(self, n):
+ """Act as if n calls to random() were made, but quickly.
+
+ n is an int, greater than or equal to 0.
+
+ Example use: If you have 2 threads and know that each will
+ consume no more than a million random numbers, create two Random
+ objects r1 and r2, then do
+ r2.setstate(r1.getstate())
+ r2.jumpahead(1000000)
+ Then r1 and r2 will use guaranteed-disjoint segments of the full
+ period.
+ """
+
+ if not n >= 0:
+ raise ValueError("n must be >= 0")
+ x, y, z = self._seed
+ x = int(x * pow(171, n, 30269)) % 30269
+ y = int(y * pow(172, n, 30307)) % 30307
+ z = int(z * pow(170, n, 30323)) % 30323
+ self._seed = x, y, z
+
+ def __whseed(self, x=0, y=0, z=0):
+ """Set the Wichmann-Hill seed from (x, y, z).
+
+ These must be integers in the range [0, 256).
+ """
+
+ if not type(x) == type(y) == type(z) == int:
+ raise TypeError('seeds must be integers')
+ if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
+ raise ValueError('seeds must be in range(0, 256)')
+ if 0 == x == y == z:
+ # Initialize from current time
+ import time
+ t = long(time.time() * 256)
+ t = int((t&0xffffff) ^ (t>>24))
+ t, x = divmod(t, 256)
+ t, y = divmod(t, 256)
+ t, z = divmod(t, 256)
+ # Zero is a poor seed, so substitute 1
+ self._seed = (x or 1, y or 1, z or 1)
+
+ self.gauss_next = None
+
+ def whseed(self, a=None):
+ """Seed from hashable object's hash code.
+
+ None or no argument seeds from current time. It is not guaranteed
+ that objects with distinct hash codes lead to distinct internal
+ states.
+
+ This is obsolete, provided for compatibility with the seed routine
+ used prior to Python 2.1. Use the .seed() method instead.
+ """
+
+ if a is None:
+ self.__whseed()
+ return
+ a = hash(a)
+ a, x = divmod(a, 256)
+ a, y = divmod(a, 256)
+ a, z = divmod(a, 256)
+ x = (x + a) % 256 or 1
+ y = (y + a) % 256 or 1
+ z = (z + a) % 256 or 1
+ self.__whseed(x, y, z)
+
## -------------------- test program --------------------
***************
*** 811,825 ****
_test_sample(500)
! # Test jumpahead.
! s = getstate()
! jumpahead(N)
! r1 = random()
! # now do it the slow way
! setstate(s)
! for i in range(N):
! random()
! r2 = random()
! if r1 != r2:
! raise ValueError("jumpahead test failed " + `(N, r1, r2)`)
# Create one instance, seeded from current time, and export its methods
--- 843,860 ----
_test_sample(500)
! if isinstance(_inst, WichmannHill):
! # Test jumpahead.
! s = getstate()
! jumpahead(N)
! r1 = random()
! # now do it the slow way
! setstate(s)
! for i in range(N):
! random()
! r2 = random()
! if r1 != r2:
! raise ValueError("jumpahead test failed " + `(N, r1, r2)`)
!
! print _inst.__class__ # Identify which generator was tested
# Create one instance, seeded from current time, and export its methods
***************
*** 828,835 ****
--- 863,873 ----
# libraries), but that's fine for most programs and is easier for the
# casual user than making them instantiate their own Random() instance.
+
_inst = Random()
+ #_inst = WichmannHill()
seed = _inst.seed
random = _inst.random
uniform = _inst.uniform
+ _randbelow = _inst._randbelow
randint = _inst.randint
choice = _inst.choice
***************
*** 851,855 ****
setstate = _inst.setstate
jumpahead = _inst.jumpahead
! whseed = _inst.whseed
if __name__ == '__main__':
--- 889,896 ----
setstate = _inst.setstate
jumpahead = _inst.jumpahead
! try:
! whseed = _inst.whseed
! except AttributeError:
! pass
if __name__ == '__main__':