]>
Commit | Line | Data |
---|---|---|
970d7e83 | 1 | #!/usr/bin/env python |
1a4d82fc JJ |
2 | # |
3 | # Copyright 2013 The Rust Project Developers. See the COPYRIGHT | |
4 | # file at the top-level directory of this distribution and at | |
5 | # http://rust-lang.org/COPYRIGHT. | |
6 | # | |
7 | # Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or | |
8 | # http://www.apache.org/licenses/LICENSE-2.0> or the MIT license | |
9 | # <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your | |
10 | # option. This file may not be copied, modified, or distributed | |
11 | # except according to those terms. | |
970d7e83 LB |
12 | |
13 | # This creates the tables used for distributions implemented using the | |
1a4d82fc | 14 | # ziggurat algorithm in `rand::distributions;`. They are |
970d7e83 LB |
15 | # (basically) the tables as used in the ZIGNOR variant (Doornik 2005). |
16 | # They are changed rarely, so the generated file should be checked in | |
17 | # to git. | |
18 | # | |
19 | # It creates 3 tables: X as in the paper, F which is f(x_i), and | |
20 | # F_DIFF which is f(x_i) - f(x_{i-1}). The latter two are just cached | |
21 | # values which is not done in that paper (but is done in other | |
22 | # variants). Note that the adZigR table is unnecessary because of | |
23 | # algebra. | |
24 | # | |
25 | # It is designed to be compatible with Python 2 and 3. | |
26 | ||
27 | from math import exp, sqrt, log, floor | |
28 | import random | |
29 | ||
30 | # The order should match the return value of `tables` | |
1a4d82fc | 31 | TABLE_NAMES = ['X', 'F'] |
970d7e83 LB |
32 | |
33 | # The actual length of the table is 1 more, to stop | |
34 | # index-out-of-bounds errors. This should match the bitwise operation | |
35 | # to find `i` in `zigurrat` in `libstd/rand/mod.rs`. Also the *_R and | |
36 | # *_V constants below depend on this value. | |
37 | TABLE_LEN = 256 | |
38 | ||
39 | # equivalent to `zigNorInit` in Doornik2005, but generalised to any | |
40 | # distribution. r = dR, v = dV, f = probability density function, | |
41 | # f_inv = inverse of f | |
42 | def tables(r, v, f, f_inv): | |
43 | # compute the x_i | |
44 | xvec = [0]*(TABLE_LEN+1) | |
45 | ||
46 | xvec[0] = v / f(r) | |
47 | xvec[1] = r | |
48 | ||
49 | for i in range(2, TABLE_LEN): | |
50 | last = xvec[i-1] | |
51 | xvec[i] = f_inv(v / last + f(last)) | |
52 | ||
53 | # cache the f's | |
54 | fvec = [0]*(TABLE_LEN+1) | |
970d7e83 LB |
55 | for i in range(TABLE_LEN+1): |
56 | fvec[i] = f(xvec[i]) | |
970d7e83 | 57 | |
1a4d82fc | 58 | return xvec, fvec |
970d7e83 LB |
59 | |
60 | # Distributions | |
61 | # N(0, 1) | |
62 | def norm_f(x): | |
63 | return exp(-x*x/2.0) | |
64 | def norm_f_inv(y): | |
65 | return sqrt(-2.0*log(y)) | |
66 | ||
67 | NORM_R = 3.6541528853610088 | |
68 | NORM_V = 0.00492867323399 | |
69 | ||
70 | NORM = tables(NORM_R, NORM_V, | |
71 | norm_f, norm_f_inv) | |
72 | ||
73 | # Exp(1) | |
74 | def exp_f(x): | |
75 | return exp(-x) | |
76 | def exp_f_inv(y): | |
77 | return -log(y) | |
78 | ||
79 | EXP_R = 7.69711747013104972 | |
80 | EXP_V = 0.0039496598225815571993 | |
81 | ||
82 | EXP = tables(EXP_R, EXP_V, | |
83 | exp_f, exp_f_inv) | |
84 | ||
85 | ||
86 | # Output the tables/constants/types | |
87 | ||
88 | def render_static(name, type, value): | |
89 | # no space or | |
90 | return 'pub static %s: %s =%s;\n' % (name, type, value) | |
91 | ||
92 | # static `name`: [`type`, .. `len(values)`] = | |
93 | # [values[0], ..., values[3], | |
94 | # values[4], ..., values[7], | |
95 | # ... ]; | |
96 | def render_table(name, values): | |
97 | rows = [] | |
98 | # 4 values on each row | |
99 | for i in range(0, len(values), 4): | |
100 | row = values[i:i+4] | |
101 | rows.append(', '.join('%.18f' % f for f in row)) | |
102 | ||
103 | rendered = '\n [%s]' % ',\n '.join(rows) | |
104 | return render_static(name, '[f64, .. %d]' % len(values), rendered) | |
105 | ||
106 | ||
107 | with open('ziggurat_tables.rs', 'w') as f: | |
108 | f.write('''// Copyright 2013 The Rust Project Developers. See the COPYRIGHT | |
109 | // file at the top-level directory of this distribution and at | |
110 | // http://rust-lang.org/COPYRIGHT. | |
111 | // | |
112 | // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or | |
113 | // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license | |
114 | // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your | |
115 | // option. This file may not be copied, modified, or distributed | |
116 | // except according to those terms. | |
117 | ||
118 | // Tables for distributions which are sampled using the ziggurat | |
119 | // algorithm. Autogenerated by `ziggurat_tables.py`. | |
120 | ||
121 | pub type ZigTable = &\'static [f64, .. %d]; | |
122 | ''' % (TABLE_LEN + 1)) | |
123 | for name, tables, r in [('NORM', NORM, NORM_R), | |
124 | ('EXP', EXP, EXP_R)]: | |
125 | f.write(render_static('ZIG_%s_R' % name, 'f64', ' %.18f' % r)) | |
126 | for (tabname, table) in zip(TABLE_NAMES, tables): | |
127 | f.write(render_table('ZIG_%s_%s' % (name, tabname), table)) |