Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# Copyright 2013 Kwant authors. 

2# 

3# This file is part of Kwant. It is subject to the license terms in the file 

4# LICENSE.rst found in the top-level directory of this distribution and at 

5# https://kwant-project.org/license. A list of Kwant authors can be found in 

6# the file AUTHORS.rst at the top-level directory of this distribution and at 

7# https://kwant-project.org/authors. 

8 

9"""Random-access random numbers 

10 

11This module provides routines that given some input compute a "random" output 

12that depends on the input in a (cryptographically) intractable way. 

13 

14This turns out to be very useful when one needs to map some irregular objects to 

15random numbers in a deterministic and reproducible way. 

16 

17Internally, the md5 hash algorithm is used. The randomness thus generated is 

18good enough to pass the "dieharder" battery of tests: see the function `test` of 

19this module. 

20""" 

21 

22from math import pi, log, sqrt, cos 

23from hashlib import md5 

24from struct import unpack 

25 

26__all__ = ['uniform', 'gauss', 'test'] 

27 

28 

29TWOPI = 2 * pi 

30BPF = 53 # Number of bits in a Python float. 

31BPF_MASK = 2**53 - 1 

32RECIP_BPF = 2**-BPF 

33 

34 

35def str_to_bytes(s): 

36 """Return bytes if the input is a string, else return the object as-is.""" 

37 try: 

38 return s.encode('utf8') 

39 except AttributeError: 

40 return s 

41 

42def uniform2(input, salt=''): 

43 """Return two independent [0,1)-distributed numbers.""" 

44 input = str_to_bytes(input) 

45 salt = str_to_bytes(salt) 

46 input = memoryview(input).tobytes() + salt 

47 a, b = unpack('qq', md5(input).digest()) 

48 a &= BPF_MASK 

49 b &= BPF_MASK 

50 return a * RECIP_BPF, b * RECIP_BPF 

51 

52 

53def uniform(input, salt=''): 

54 """md5-hash `input` and `salt` and map the result to the [0,1) interval. 

55 

56 `input` must be some object that supports the buffer protocol (i.e. a string 

57 or a numpy/tinyarray array). `salt` must be a string or a bytes object. 

58 """ 

59 return uniform2(input, salt)[0] 

60 

61 

62def gauss(input, salt=''): 

63 """md5-hash `input` and `salt` and return the result as a standard normal 

64 distributed variable. 

65 

66 `input` must be some object that supports the buffer protocol (i.e. a string 

67 or a numpy/tinyarray array). `salt` must be a string or a bytes object. 

68 """ 

69 # This uses the Box-Muller transform. Only one of the two results is 

70 # computed. 

71 a, b = uniform2(input, salt) 

72 return cos(a * TWOPI) * sqrt(-2.0 * log(1.0 - b)) 

73 

74 

75def test(n=20000): # skip coverage 

76 """Test the generator with the dieharder suite generating n**2 samples. 

77 

78 Executing this function may take a very long time. 

79 """ 

80 import os 

81 import tempfile 

82 import subprocess 

83 from tinyarray import array 

84 from struct import pack 

85 

86 f = tempfile.NamedTemporaryFile(delete=False) 

87 try: 

88 for x in range(n): 

89 for y in range(n): 

90 a = array((x, y)) 

91 i = int(2**32 * uniform(a)) 

92 f.write(pack('I', i)) 

93 f.close() 

94 subprocess.call(['dieharder', '-a', '-g', '201', '-f', f.name]) 

95 finally: 

96 os.remove(f.name)