-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsobol_new.py
151 lines (131 loc) · 4.8 KB
/
sobol_new.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 28 11:09:42 2019
@author: jmadriga
"""
# -*- coding: utf-8 -*-
"""sobol_point_gen.py
Usage:
sobol_point_gen.py N D <dimensions_file>
sobol_point_gen.py -h | --help
Options:
-h --help Show this screen.
Arguments:
N Number of points to generate (<2**32)
D Number of dimensions for each point. (D <= max D in dimensions_file)
dimensions_file See: https://web.maths.unsw.edu.au/~fkuo/sobol/
Created from an original c++ source by Frances Y. Kuo and Stephen Joe
Created on Fri Nov 1 09:13:12 2019
@author: Paddy3118
"""
from functools import lru_cache
import numpy as np
import matplotlib.pyplot as plt
#%%
def binhi1(b):
"highest 1 bit of binary b. (== ceil(log2(b)))"
pos = 0
while b > 0:
pos, b = pos + 1, b // 2
return pos
def binlo0(b):
"lowest 0 bit from right of binary b. Always >=1"
pos, b2 = 1, b // 2
while b2 * 2 != b:
pos, b, b2 = pos + 1, b2, b2 // 2
return pos
#%%
def direction_file_reader(d, fname="Sobol_new-joe-kuo-6.21201"):
"""\
Reads d+header lines from direction file of format:
d s a m_i
2 1 0 1
3 2 1 1 3
4 3 1 1 3 1
5 3 2 1 1 1
6 4 1 1 1 3 3
7 4 4 1 3 5 13
8 5 2 1 1 5 5 17
9 5 4 1 1 5 5 5
10 5 7 1 1 7 11 19
11 5 11 1 1 5 1 1
12 5 13 1 1 1 3 11
13 5 14 1 3 5 5 31
14 6 1 1 3 3 9 7 49
15 6 13 1 1 1 15 21 21
16 6 16 1 3 1 13 27 49
17 6 19 1 1 1 15 7 5
18 6 22 1 3 1 15 13 25
19 6 25 1 1 5 5 19 61
20 7 1 1 3 7 11 23 15 103
21 7 4 1 3 7 13 13 15 69
22 7 7 1 1 3 13 7 35 63
23 7 8 1 3 5 9 1 25 53
24 7 14 1 3 1 13 9 35 107
25 7 19 1 3 1 5 27 61 31
26 7 21 1 1 5 11 19 41 61
27 7 28 1 3 5 3 3 13 69
28 7 31 1 1 7 13 1 19 1
29 7 32 1 3 7 5 13 19 59
30 7 37 1 1 3 9 25 29 41
31 7 41 1 3 5 13 23 1 55
32 7 42 1 3 7 3 13 59 17
...
"""
with open(fname) as f:
headers = f.readline().strip().split()
assert headers == ['d', 's', 'a', 'm_i']
directions = []
for i, line in zip(range(d), f):
d, s, a, *m = [int(cell) for cell in line.strip().split()]
directions.append((d, s, a, tuple([0] + m)))
return directions
def generate_points(N, D, fname="Sobol_new-joe-kuo-6.21201"):
return np.asarray(list(sobol_point_gen(N, D)))
#%%
@lru_cache(128)
def calc_dnum(s, a, m, mxbits): # V, direction numbers
if mxbits <= s:
dirnum = [0] + [m[i] << (32 - i) for i in range(1, mxbits+1)]
else:
dirnum = [0] + [m[i] << (32 - i) for i in range(1, s+1)]
for i in range(s+1, mxbits+1):
dirnum.append(dirnum[i-s] ^ (dirnum[i-s] >> s))
for k in range(1, s):
dirnum[i] ^= (((a >> (s-1-k)) & 1) * dirnum[i-k])
return dirnum
def sobol_point_gen(N, D, fname="Sobol_new-joe-kuo-6.21201"):
pow2_32 = 2**32
mxbits = binhi1(N) # L, max bits needed
dnum = [1 << (32- i) for i in range(mxbits+1)] # V, direction numbers
directions = direction_file_reader(D, fname)
xd = [0 for j in range(D)] # X
pt = [0 for j in range(D)] # POINTS
yield pt
for n in range(1, N):
nfirst0 = binlo0(n-1) # C
pt = [0 for j in range(D)]
xd[0] = xd[0] ^ dnum[nfirst0]
pt[0] = xd[0] / pow2_32
for jj in range(1, D):
_, s, a, m = directions[jj-1]
dnum2 = calc_dnum(s, a, m, mxbits) # V, direction numbers
xd[jj] = xd[jj] ^ dnum2[nfirst0]
pt[jj] = xd[jj]/pow2_32
yield pt
#%%
if __name__ == '__main__':
first_sobol_3D_points = [
[0, 0, 0],
[0.5, 0.5, 0.5],
[0.75, 0.25, 0.25],
[0.25, 0.75, 0.75],
[0.375, 0.375, 0.625],
[0.875, 0.875, 0.125],
[0.625, 0.125, 0.875],
[0.125, 0.625, 0.375],
[0.1875, 0.3125, 0.9375],
[0.6875, 0.8125, 0.4375]]
generated = list(sobol_point_gen(10, 3))
assert generated == first_sobol_3D_points