github.com/grumpyhome/grumpy@v0.3.1-0.20201208125205-7b775405bdf1/grumpy-runtime-src/third_party/stdlib/random.py (about) 1 """Random variable generators. 2 3 integers 4 -------- 5 uniform within range 6 7 sequences 8 --------- 9 pick random element 10 pick random sample 11 generate random permutation 12 13 distributions on the real line: 14 ------------------------------ 15 uniform 16 triangular 17 normal (Gaussian) 18 lognormal 19 negative exponential 20 gamma 21 beta 22 pareto 23 Weibull 24 25 distributions on the circle (angles 0 to 2pi) 26 --------------------------------------------- 27 circular uniform 28 von Mises 29 30 General notes on the underlying Mersenne Twister core generator: 31 32 * The period is 2**19937-1. 33 * It is one of the most extensively tested generators in existence. 34 * Without a direct way to compute N steps forward, the semantics of 35 jumpahead(n) are weakened to simply jump to another distant state and rely 36 on the large period to avoid overlapping sequences. 37 * The random() method is implemented in C, executes in a single Python step, 38 and is, therefore, threadsafe. 39 40 """ 41 42 #from warnings import warn as _warn 43 #from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 44 #from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil 45 #from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 46 #from os import urandom as _urandom 47 #from binascii import hexlify as _hexlify 48 #import hashlib as _hashlib 49 50 __all__ = ["Random","seed","random","uniform","randint","choice","sample", 51 "randrange","shuffle","normalvariate","lognormvariate", 52 "expovariate","vonmisesvariate","gammavariate","triangular", 53 "gauss","betavariate","paretovariate","weibullvariate", 54 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 55 "SystemRandom"] 56 57 # Import _random wrapper from grumpy std lib. 58 # It is used as a replacement for the CPython random generator. 59 60 import _random 61 62 # NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 63 # TWOPI = 2.0*_pi 64 # LOG4 = _log(4.0) 65 # SG_MAGICCONST = 1.0 + _log(4.5) 66 BPF = _random.BPF 67 RECIP_BPF = _random.RECIP_BPF 68 69 70 class Random(_random.GrumpyRandom): 71 """Random number generator base class used by bound module functions. 72 73 Used to instantiate instances of Random to get generators that don't 74 share state. Especially useful for multi-threaded programs, creating 75 a different instance of Random for each thread, and using the jumpahead() 76 method to ensure that the generated sequences seen by each thread don't 77 overlap. 78 79 Class Random can also be subclassed if you want to use a different basic 80 generator of your own devising: in that case, override the following 81 methods: random(), seed(), getstate(), setstate() and jumpahead(). 82 Optionally, implement a getrandbits() method so that randrange() can cover 83 arbitrarily large ranges. 84 85 """ 86 87 VERSION = 3 # used by getstate/setstate 88 89 def __init__(self, x=None): 90 """Initialize an instance. 91 92 Optional argument x controls seeding, as for Random.seed(). 93 """ 94 95 self.seed(x) 96 self.gauss_next = None 97 98 def seed(self, a=None): 99 """Initialize internal state of the random number generator. 100 101 None or no argument seeds from current time or from an operating 102 system specific randomness source if available. 103 104 If a is not None or is an int or long, hash(a) is used instead. 105 Hash values for some types are nondeterministic when the 106 PYTHONHASHSEED environment variable is enabled. 107 """ 108 109 super(Random, self).seed(a) 110 self.gauss_next = None 111 112 ## ---- Methods below this point do not need to be overridden when 113 ## ---- subclassing for the purpose of using a different core generator. 114 115 ## -------------------- integer methods ------------------- 116 117 def randrange(self, start, stop=None, step=1, _int=int, _maxwidth=1L<<BPF): 118 """Choose a random item from range(start, stop[, step]). 119 120 This fixes the problem with randint() which includes the 121 endpoint; in Python this is usually not what you want. 122 123 """ 124 125 # This code is a bit messy to make it fast for the 126 # common case while still doing adequate error checking. 127 istart = _int(start) 128 if istart != start: 129 raise ValueError, "non-integer arg 1 for randrange()" 130 if stop is None: 131 if istart > 0: 132 if istart >= _maxwidth: 133 return self._randbelow(istart) 134 return _int(self.random() * istart) 135 raise ValueError, "empty range for randrange()" 136 137 # stop argument supplied. 138 istop = _int(stop) 139 if istop != stop: 140 raise ValueError, "non-integer stop for randrange()" 141 width = istop - istart 142 if step == 1 and width > 0: 143 # Note that 144 # int(istart + self.random()*width) 145 # instead would be incorrect. For example, consider istart 146 # = -2 and istop = 0. Then the guts would be in 147 # -2.0 to 0.0 exclusive on both ends (ignoring that random() 148 # might return 0.0), and because int() truncates toward 0, the 149 # final result would be -1 or 0 (instead of -2 or -1). 150 # istart + int(self.random()*width) 151 # would also be incorrect, for a subtler reason: the RHS 152 # can return a long, and then randrange() would also return 153 # a long, but we're supposed to return an int (for backward 154 # compatibility). 155 156 if width >= _maxwidth: 157 return _int(istart + self._randbelow(width)) 158 return _int(istart + _int(self.random()*width)) 159 if step == 1: 160 raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 161 162 # Non-unit step argument supplied. 163 istep = _int(step) 164 if istep != step: 165 raise ValueError, "non-integer step for randrange()" 166 if istep > 0: 167 n = (width + istep - 1) // istep 168 elif istep < 0: 169 n = (width + istep + 1) // istep 170 else: 171 raise ValueError, "zero step for randrange()" 172 173 if n <= 0: 174 raise ValueError, "empty range for randrange()" 175 176 if n >= _maxwidth: 177 return istart + istep*self._randbelow(n) 178 return istart + istep*_int(self.random() * n) 179 180 def randint(self, a, b): 181 """Return random integer in range [a, b], including both end points. 182 """ 183 184 return self.randrange(a, b+1) 185 186 ## -------------------- sequence methods ------------------- 187 188 def choice(self, seq): 189 """Choose a random element from a non-empty sequence.""" 190 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 191 192 def shuffle(self, x, random=None): 193 """x, random=random.random -> shuffle list x in place; return None. 194 195 Optional arg random is a 0-argument function returning a random 196 float in [0.0, 1.0); by default, the standard random.random. 197 198 """ 199 200 if random is None: 201 random = self.random 202 _int = int 203 for i in reversed(xrange(1, len(x))): 204 # pick an element in x[:i+1] with which to exchange x[i] 205 j = _int(random() * (i+1)) 206 x[i], x[j] = x[j], x[i] 207 208 # def sample(self, population, k): 209 # """Chooses k unique random elements from a population sequence. 210 211 # Returns a new list containing elements from the population while 212 # leaving the original population unchanged. The resulting list is 213 # in selection order so that all sub-slices will also be valid random 214 # samples. This allows raffle winners (the sample) to be partitioned 215 # into grand prize and second place winners (the subslices). 216 217 # Members of the population need not be hashable or unique. If the 218 # population contains repeats, then each occurrence is a possible 219 # selection in the sample. 220 221 # To choose a sample in a range of integers, use xrange as an argument. 222 # This is especially fast and space efficient for sampling from a 223 # large population: sample(xrange(10000000), 60) 224 # """ 225 226 # # Sampling without replacement entails tracking either potential 227 # # selections (the pool) in a list or previous selections in a set. 228 229 # # When the number of selections is small compared to the 230 # # population, then tracking selections is efficient, requiring 231 # # only a small set and an occasional reselection. For 232 # # a larger number of selections, the pool tracking method is 233 # # preferred since the list takes less space than the 234 # # set and it doesn't suffer from frequent reselections. 235 236 # n = len(population) 237 # if not 0 <= k <= n: 238 # raise ValueError("sample larger than population") 239 # random = self.random 240 # _int = int 241 # result = [None] * k 242 # setsize = 21 # size of a small set minus size of an empty list 243 # if k > 5: 244 # setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets 245 # if n <= setsize or hasattr(population, "keys"): 246 # # An n-length list is smaller than a k-length set, or this is a 247 # # mapping type so the other algorithm wouldn't work. 248 # pool = list(population) 249 # for i in xrange(k): # invariant: non-selected at [0,n-i) 250 # j = _int(random() * (n-i)) 251 # result[i] = pool[j] 252 # pool[j] = pool[n-i-1] # move non-selected item into vacancy 253 # else: 254 # try: 255 # selected = set() 256 # selected_add = selected.add 257 # for i in xrange(k): 258 # j = _int(random() * n) 259 # while j in selected: 260 # j = _int(random() * n) 261 # selected_add(j) 262 # result[i] = population[j] 263 # except (TypeError, KeyError): # handle (at least) sets 264 # if isinstance(population, list): 265 # raise 266 # return self.sample(tuple(population), k) 267 # return result 268 269 ## -------------------- real-valued distributions ------------------- 270 271 ## -------------------- uniform distribution ------------------- 272 273 def uniform(self, a, b): 274 "Get a random number in the range [a, b) or [a, b] depending on rounding." 275 return a + (b-a) * self.random() 276 277 ## -------------------- triangular -------------------- 278 279 # def triangular(self, low=0.0, high=1.0, mode=None): 280 # """Triangular distribution. 281 282 # Continuous distribution bounded by given lower and upper limits, 283 # and having a given mode value in-between. 284 285 # http://en.wikipedia.org/wiki/Triangular_distribution 286 287 # """ 288 # u = self.random() 289 # try: 290 # c = 0.5 if mode is None else (mode - low) / (high - low) 291 # except ZeroDivisionError: 292 # return low 293 # if u > c: 294 # u = 1.0 - u 295 # c = 1.0 - c 296 # low, high = high, low 297 # return low + (high - low) * (u * c) ** 0.5 298 299 ## -------------------- normal distribution -------------------- 300 301 # def normalvariate(self, mu, sigma): 302 # """Normal distribution. 303 304 # mu is the mean, and sigma is the standard deviation. 305 306 # """ 307 # # mu = mean, sigma = standard deviation 308 309 # # Uses Kinderman and Monahan method. Reference: Kinderman, 310 # # A.J. and Monahan, J.F., "Computer generation of random 311 # # variables using the ratio of uniform deviates", ACM Trans 312 # # Math Software, 3, (1977), pp257-260. 313 314 # random = self.random 315 # while 1: 316 # u1 = random() 317 # u2 = 1.0 - random() 318 # z = NV_MAGICCONST*(u1-0.5)/u2 319 # zz = z*z/4.0 320 # if zz <= -_log(u2): 321 # break 322 # return mu + z*sigma 323 324 ## -------------------- lognormal distribution -------------------- 325 326 # def lognormvariate(self, mu, sigma): 327 # """Log normal distribution. 328 329 # If you take the natural logarithm of this distribution, you'll get a 330 # normal distribution with mean mu and standard deviation sigma. 331 # mu can have any value, and sigma must be greater than zero. 332 333 # """ 334 # return _exp(self.normalvariate(mu, sigma)) 335 336 ## -------------------- exponential distribution -------------------- 337 338 # def expovariate(self, lambd): 339 # """Exponential distribution. 340 341 # lambd is 1.0 divided by the desired mean. It should be 342 # nonzero. (The parameter would be called "lambda", but that is 343 # a reserved word in Python.) Returned values range from 0 to 344 # positive infinity if lambd is positive, and from negative 345 # infinity to 0 if lambd is negative. 346 347 # """ 348 # # lambd: rate lambd = 1/mean 349 # # ('lambda' is a Python reserved word) 350 351 # # we use 1-random() instead of random() to preclude the 352 # # possibility of taking the log of zero. 353 # return -_log(1.0 - self.random())/lambd 354 355 ## -------------------- von Mises distribution -------------------- 356 357 # def vonmisesvariate(self, mu, kappa): 358 # """Circular data distribution. 359 360 # mu is the mean angle, expressed in radians between 0 and 2*pi, and 361 # kappa is the concentration parameter, which must be greater than or 362 # equal to zero. If kappa is equal to zero, this distribution reduces 363 # to a uniform random angle over the range 0 to 2*pi. 364 365 # """ 366 # # mu: mean angle (in radians between 0 and 2*pi) 367 # # kappa: concentration parameter kappa (>= 0) 368 # # if kappa = 0 generate uniform random angle 369 370 # # Based upon an algorithm published in: Fisher, N.I., 371 # # "Statistical Analysis of Circular Data", Cambridge 372 # # University Press, 1993. 373 374 # # Thanks to Magnus Kessler for a correction to the 375 # # implementation of step 4. 376 377 # random = self.random 378 # if kappa <= 1e-6: 379 # return TWOPI * random() 380 381 # s = 0.5 / kappa 382 # r = s + _sqrt(1.0 + s * s) 383 384 # while 1: 385 # u1 = random() 386 # z = _cos(_pi * u1) 387 388 # d = z / (r + z) 389 # u2 = random() 390 # if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d): 391 # break 392 393 # q = 1.0 / r 394 # f = (q + z) / (1.0 + q * z) 395 # u3 = random() 396 # if u3 > 0.5: 397 # theta = (mu + _acos(f)) % TWOPI 398 # else: 399 # theta = (mu - _acos(f)) % TWOPI 400 401 # return theta 402 403 ## -------------------- gamma distribution -------------------- 404 405 # def gammavariate(self, alpha, beta): 406 # """Gamma distribution. Not the gamma function! 407 408 # Conditions on the parameters are alpha > 0 and beta > 0. 409 410 # The probability distribution function is: 411 412 # x ** (alpha - 1) * math.exp(-x / beta) 413 # pdf(x) = -------------------------------------- 414 # math.gamma(alpha) * beta ** alpha 415 416 # """ 417 418 # # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 419 420 # # Warning: a few older sources define the gamma distribution in terms 421 # # of alpha > -1.0 422 # if alpha <= 0.0 or beta <= 0.0: 423 # raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 424 425 # random = self.random 426 # if alpha > 1.0: 427 428 # # Uses R.C.H. Cheng, "The generation of Gamma 429 # # variables with non-integral shape parameters", 430 # # Applied Statistics, (1977), 26, No. 1, p71-74 431 432 # ainv = _sqrt(2.0 * alpha - 1.0) 433 # bbb = alpha - LOG4 434 # ccc = alpha + ainv 435 436 # while 1: 437 # u1 = random() 438 # if not 1e-7 < u1 < .9999999: 439 # continue 440 # u2 = 1.0 - random() 441 # v = _log(u1/(1.0-u1))/ainv 442 # x = alpha*_exp(v) 443 # z = u1*u1*u2 444 # r = bbb+ccc*v-x 445 # if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 446 # return x * beta 447 448 # elif alpha == 1.0: 449 # # expovariate(1) 450 # u = random() 451 # while u <= 1e-7: 452 # u = random() 453 # return -_log(u) * beta 454 455 # else: # alpha is between 0 and 1 (exclusive) 456 457 # # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 458 459 # while 1: 460 # u = random() 461 # b = (_e + alpha)/_e 462 # p = b*u 463 # if p <= 1.0: 464 # x = p ** (1.0/alpha) 465 # else: 466 # x = -_log((b-p)/alpha) 467 # u1 = random() 468 # if p > 1.0: 469 # if u1 <= x ** (alpha - 1.0): 470 # break 471 # elif u1 <= _exp(-x): 472 # break 473 # return x * beta 474 475 ## -------------------- Gauss (faster alternative) -------------------- 476 477 # def gauss(self, mu, sigma): 478 # """Gaussian distribution. 479 480 # mu is the mean, and sigma is the standard deviation. This is 481 # slightly faster than the normalvariate() function. 482 483 # Not thread-safe without a lock around calls. 484 485 # """ 486 487 # # When x and y are two variables from [0, 1), uniformly 488 # # distributed, then 489 # # 490 # # cos(2*pi*x)*sqrt(-2*log(1-y)) 491 # # sin(2*pi*x)*sqrt(-2*log(1-y)) 492 # # 493 # # are two *independent* variables with normal distribution 494 # # (mu = 0, sigma = 1). 495 # # (Lambert Meertens) 496 # # (corrected version; bug discovered by Mike Miller, fixed by LM) 497 498 # # Multithreading note: When two threads call this function 499 # # simultaneously, it is possible that they will receive the 500 # # same return value. The window is very small though. To 501 # # avoid this, you have to use a lock around all calls. (I 502 # # didn't want to slow this down in the serial case by using a 503 # # lock here.) 504 505 # random = self.random 506 # z = self.gauss_next 507 # self.gauss_next = None 508 # if z is None: 509 # x2pi = random() * TWOPI 510 # g2rad = _sqrt(-2.0 * _log(1.0 - random())) 511 # z = _cos(x2pi) * g2rad 512 # self.gauss_next = _sin(x2pi) * g2rad 513 514 # return mu + z*sigma 515 516 ## -------------------- beta -------------------- 517 ## See 518 ## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html 519 ## for Ivan Frohne's insightful analysis of why the original implementation: 520 ## 521 ## def betavariate(self, alpha, beta): 522 ## # Discrete Event Simulation in C, pp 87-88. 523 ## 524 ## y = self.expovariate(alpha) 525 ## z = self.expovariate(1.0/beta) 526 ## return z/(y+z) 527 ## 528 ## was dead wrong, and how it probably got that way. 529 530 # def betavariate(self, alpha, beta): 531 # """Beta distribution. 532 533 # Conditions on the parameters are alpha > 0 and beta > 0. 534 # Returned values range between 0 and 1. 535 536 # """ 537 538 # # This version due to Janne Sinkkonen, and matches all the std 539 # # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 540 # y = self.gammavariate(alpha, 1.) 541 # if y == 0: 542 # return 0.0 543 # else: 544 # return y / (y + self.gammavariate(beta, 1.)) 545 546 ## -------------------- Pareto -------------------- 547 548 # def paretovariate(self, alpha): 549 # """Pareto distribution. alpha is the shape parameter.""" 550 # # Jain, pg. 495 551 552 # u = 1.0 - self.random() 553 # return 1.0 / pow(u, 1.0/alpha) 554 555 ## -------------------- Weibull -------------------- 556 557 # def weibullvariate(self, alpha, beta): 558 # """Weibull distribution. 559 560 # alpha is the scale parameter and beta is the shape parameter. 561 562 # """ 563 # # Jain, pg. 499; bug fix courtesy Bill Arms 564 565 # u = 1.0 - self.random() 566 # return alpha * pow(-_log(u), 1.0/beta) 567 568 ## -------------------- test program -------------------- 569 570 def _test_generator(n, func, args): 571 import time 572 print n, 'times', func.__name__ 573 total = 0.0 574 sqsum = 0.0 575 smallest = 1e10 576 largest = -1e10 577 t0 = time.time() 578 for i in range(n): 579 x = func(*args) 580 total += x 581 sqsum = sqsum + x*x 582 smallest = min(x, smallest) 583 largest = max(x, largest) 584 t1 = time.time() 585 print round(t1-t0, 3), 'sec,', 586 avg = total/n 587 stddev = _sqrt(sqsum/n - avg*avg) 588 print 'avg %g, stddev %g, min %g, max %g' % \ 589 (avg, stddev, smallest, largest) 590 591 592 def _test(N=2000): 593 _test_generator(N, random, ()) 594 _test_generator(N, normalvariate, (0.0, 1.0)) 595 _test_generator(N, lognormvariate, (0.0, 1.0)) 596 _test_generator(N, vonmisesvariate, (0.0, 1.0)) 597 _test_generator(N, gammavariate, (0.01, 1.0)) 598 _test_generator(N, gammavariate, (0.1, 1.0)) 599 _test_generator(N, gammavariate, (0.1, 2.0)) 600 _test_generator(N, gammavariate, (0.5, 1.0)) 601 _test_generator(N, gammavariate, (0.9, 1.0)) 602 _test_generator(N, gammavariate, (1.0, 1.0)) 603 _test_generator(N, gammavariate, (2.0, 1.0)) 604 _test_generator(N, gammavariate, (20.0, 1.0)) 605 _test_generator(N, gammavariate, (200.0, 1.0)) 606 _test_generator(N, gauss, (0.0, 1.0)) 607 _test_generator(N, betavariate, (3.0, 3.0)) 608 _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0)) 609 610 # Create one instance, seeded from current time, and export its methods 611 # as module-level functions. The functions share state across all uses 612 #(both in the user's code and in the Python libraries), but that's fine 613 # for most programs and is easier for the casual user than making them 614 # instantiate their own Random() instance. 615 616 _inst = Random() 617 seed = _inst.seed 618 random = _inst.random 619 randint = _inst.randint 620 choice = _inst.choice 621 randrange = _inst.randrange 622 getrandbits = _inst.getrandbits 623 getstate = _inst.getstate 624 setstate = _inst.setstate 625 uniform = _inst.uniform 626 627 628 def _notimplemented(*args, **kwargs): 629 raise NotImplementedError 630 631 632 shuffle = _notimplemented 633 choices = _notimplemented 634 sample = _notimplemented 635 triangular = _notimplemented 636 normalvariate = _notimplemented 637 lognormvariate = _notimplemented 638 expovariate = _notimplemented 639 vonmisesvariate = _notimplemented 640 gammavariate = _notimplemented 641 gauss = _notimplemented 642 betavariate = _notimplemented 643 paretovariate = _notimplemented 644 weibullvariate = _notimplemented 645 646 647 if __name__ == '__main__': 648 pass 649 #_test()