1 /* Copyright (c) 2007 Massachusetts Institute of Technology
2 *
3 * Permission is hereby granted, free of charge, to any person obtaining
4 * a copy of this software and associated documentation files (the
5 * "Software"), to deal in the Software without restriction, including
6 * without limitation the rights to use, copy, modify, merge, publish,
7 * distribute, sublicense, and/or sell copies of the Software, and to
8 * permit persons to whom the Software is furnished to do so, subject to
9 * the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21 */
22
23 /* Generation of Sobol sequences in up to 1111 dimensions, based on the
24 algorithms described in:
25 P. Bratley and B. L. Fox, Algorithm 659, ACM Trans.
26 Math. Soft. 14 (1), 88-100 (1988),
27 as modified by:
28 S. Joe and F. Y. Kuo, ACM Trans. Math. Soft 29 (1), 49-57 (2003).
29
30 Note that the code below was written without even looking at the
31 Fortran code from the TOMS paper, which is only semi-free (being
32 under the restrictive ACM copyright terms). Then I went to the
33 Fortran code and took out the table of primitive polynomials and
34 starting direction #'s ... since this is just a table of numbers
35 generated by a deterministic algorithm, it is not copyrightable.
36 (Obviously, the format of these tables then necessitated some
37 slight modifications to the code.)
38
39 For the test integral of Joe and Kuo (see the main() program
40 below), I get exactly the same results for integrals up to 1111
41 dimensions compared to the table of published numbers (to the 5
42 published significant digits).
43
44 This is not to say that the authors above should not be credited for
45 their clear description of the algorithm (and their tabulation of
46 the critical numbers). Please cite them. Just that I needed
47 a free/open-source implementation. */
48
49 #include <stdlib.h>
50 #include <math.h>
51
52 #include "nlopt-util.h"
53
54 #if defined(HAVE_STDINT_H)
55 # include <stdint.h>
56 #endif
57
58 #ifndef HAVE_UINT32_T
59 # if SIZEOF_UNSIGNED_LONG == 4
60 typedef unsigned long uint32_t;
61 # elif SIZEOF_UNSIGNED_INT == 4
62 typedef unsigned int uint32_t;
63 # else
64 # error No 32-bit unsigned integer type
65 # endif
66 #endif
67
68 typedef struct nlopt_soboldata_s {
69 unsigned sdim; /* dimension of sequence being generated */
70 uint32_t *mdata; /* array of length 32 * sdim */
71 uint32_t *m[32]; /* more convenient pointers to mdata, of direction #s */
72 uint32_t *x; /* previous x = x_n, array of length sdim */
73 unsigned *b; /* position of fixed point in x[i] is after bit b[i] */
74 uint32_t n; /* number of x's generated so far */
75 } soboldata;
76
77 /* Return position (0, 1, ...) of rightmost (least-significant) zero bit in n.
78 *
79 * This code uses a 32-bit version of algorithm to find the rightmost
80 * one bit in Knuth, _The Art of Computer Programming_, volume 4A
81 * (draft fascicle), section 7.1.3, "Bitwise tricks and
82 * techniques."
83 *
84 * Assumes n has a zero bit, i.e. n < 2^32 - 1.
85 *
86 */
rightzero32(uint32_t n)87 static unsigned rightzero32(uint32_t n)
88 {
89 #if defined(__GNUC__) && \
90 ((__GNUC__ == 3 && __GNUC_MINOR__ >= 4) || __GNUC__ > 3)
91 return __builtin_ctz(~n); /* gcc builtin for version >= 3.4 */
92 #else
93 const uint32_t a = 0x05f66a47; /* magic number, found by brute force */
94 static const unsigned decode[32] = {0,1,2,26,23,3,15,27,24,21,19,4,12,16,28,6,31,25,22,14,20,18,11,5,30,13,17,10,29,9,8,7};
95 n = ~n; /* change to rightmost-one problem */
96 n = a * (n & (-n)); /* store in n to make sure mult. is 32 bits */
97 return decode[n >> 27];
98 #endif
99 }
100
101 /* generate the next term x_{n+1} in the Sobol sequence, as an array
102 x[sdim] of numbers in (0,1). Returns 1 on success, 0 on failure
103 (if too many #'s generated) */
sobol_gen(soboldata * sd,double * x)104 static int sobol_gen(soboldata *sd, double *x)
105 {
106 unsigned c, b, i, sdim;
107
108 if (sd->n == 4294967295U) return 0; /* n == 2^32 - 1 ... we would
109 need to switch to a 64-bit version
110 to generate more terms. */
111 c = rightzero32(sd->n++);
112 sdim = sd->sdim;
113 for (i = 0; i < sdim; ++i) {
114 b = sd->b[i];
115 if (b >= c) {
116 sd->x[i] ^= sd->m[c][i] << (b - c);
117 x[i] = ((double) (sd->x[i])) / (1U << (b+1));
118 }
119 else {
120 sd->x[i] = (sd->x[i] << (c - b)) ^ sd->m[c][i];
121 sd->b[i] = c;
122 x[i] = ((double) (sd->x[i])) / (1U << (c+1));
123 }
124 }
125 return 1;
126 }
127
128 #include "soboldata.h"
129
sobol_init(soboldata * sd,unsigned sdim)130 static int sobol_init(soboldata *sd, unsigned sdim)
131 {
132 unsigned i,j;
133
134 if (!sdim || sdim > MAXDIM) return 0;
135
136 sd->mdata = (uint32_t *) malloc(sizeof(uint32_t) * (sdim * 32));
137 if (!sd->mdata) return 0;
138
139 for (j = 0; j < 32; ++j) {
140 sd->m[j] = sd->mdata + j * sdim;
141 sd->m[j][0] = 1; /* special-case Sobol sequence */
142 }
143 for (i = 1; i < sdim; ++i) {
144 uint32_t a = sobol_a[i-1];
145 unsigned d = 0, k;
146
147 while (a) {
148 ++d;
149 a >>= 1;
150 }
151 d--; /* d is now degree of poly */
152
153 /* set initial values of m from table */
154 for (j = 0; j < d; ++j)
155 sd->m[j][i] = sobol_minit[j][i-1];
156
157 /* fill in remaining values using recurrence */
158 for (j = d; j < 32; ++j) {
159 a = sobol_a[i-1];
160 sd->m[j][i] = sd->m[j - d][i];
161 for (k = 0; k < d; ++k) {
162 sd->m[j][i] ^= ((a & 1) * sd->m[j-d+k][i]) << (d-k);
163 a >>= 1;
164 }
165 }
166 }
167
168 sd->x = (uint32_t *) malloc(sizeof(uint32_t) * sdim);
169 if (!sd->x) { free(sd->mdata); return 0; }
170
171 sd->b = (unsigned *) malloc(sizeof(unsigned) * sdim);
172 if (!sd->b) { free(sd->x); free(sd->mdata); return 0; }
173
174 for (i = 0; i < sdim; ++i) {
175 sd->x[i] = 0;
176 sd->b[i] = 0;
177 }
178
179 sd->n = 0;
180 sd->sdim = sdim;
181
182 return 1;
183 }
184
sobol_destroy(soboldata * sd)185 static void sobol_destroy(soboldata *sd)
186 {
187 free(sd->mdata);
188 free(sd->x);
189 free(sd->b);
190 }
191
192 /************************************************************************/
193 /* NLopt API to Sobol sequence creation, which hides soboldata structure
194 behind an opaque pointer */
195
nlopt_sobol_create(unsigned sdim)196 nlopt_sobol nlopt_sobol_create(unsigned sdim)
197 {
198 nlopt_sobol s = (nlopt_sobol) malloc(sizeof(soboldata));
199 if (!s) return NULL;
200 if (!sobol_init(s, sdim)) { free(s); return NULL; }
201 return s;
202 }
203
nlopt_sobol_destroy(nlopt_sobol s)204 extern void nlopt_sobol_destroy(nlopt_sobol s)
205 {
206 if (s) {
207 sobol_destroy(s);
208 free(s);
209 }
210 }
211
212 /* next vector x[sdim] in Sobol sequence, with each x[i] in (0,1) */
nlopt_sobol_next01(nlopt_sobol s,double * x)213 void nlopt_sobol_next01(nlopt_sobol s, double *x)
214 {
215 if (!sobol_gen(s, x)) {
216 /* fall back on pseudo random numbers in the unlikely event
217 that we exceed 2^32-1 points */
218 unsigned i;
219 for (i = 0; i < s->sdim; ++i)
220 x[i] = nlopt_urand(0.0,1.0);
221 }
222 }
223
224 /* next vector in Sobol sequence, scaled to (lb[i], ub[i]) interval */
nlopt_sobol_next(nlopt_sobol s,double * x,const double * lb,const double * ub)225 void nlopt_sobol_next(nlopt_sobol s, double *x,
226 const double *lb, const double *ub)
227 {
228 unsigned i, sdim;
229 nlopt_sobol_next01(s, x);
230 for (sdim = s->sdim, i = 0; i < sdim; ++i)
231 x[i] = lb[i] + (ub[i] - lb[i]) * x[i];
232 }
233
234 /* if we know in advance how many points (n) we want to compute, then
235 adopt the suggestion of the Joe and Kuo paper, which in turn
236 is taken from Acworth et al (1998), of skipping a number of
237 points equal to the largest power of 2 smaller than n */
nlopt_sobol_skip(nlopt_sobol s,unsigned n,double * x)238 void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x)
239 {
240 if (s) {
241 unsigned k = 1;
242 while (k*2 < n) k *= 2;
243 while (k-- > 0) sobol_gen(s, x);
244 }
245 }
246
247 /************************************************************************/
248
249 /* compile with -DSOBOLSEQ_TEST for test program */
250 #ifdef SOBOLSEQ_TEST
251 #include <stdio.h>
252 #include <time.h>
253
254 /* test integrand from Joe and Kuo paper ... integrates to 1 */
testfunc(unsigned n,const double * x)255 static double testfunc(unsigned n, const double *x)
256 {
257 double f = 1;
258 unsigned j;
259 for (j = 1; j <= n; ++j) {
260 double cj = pow((double) j, 0.3333333333333333333);
261 f *= (fabs(4*x[j-1] - 2) + cj) / (1 + cj);
262 }
263 return f;
264 }
265
main(int argc,char ** argv)266 int main(int argc, char **argv)
267 {
268 unsigned n, j, i, sdim;
269 static double x[MAXDIM];
270 double testint_sobol = 0, testint_rand = 0;
271 nlopt_sobol s;
272 if (argc < 3) {
273 fprintf(stderr, "Usage: %s <sdim> <ngen>\n", argv[0]);
274 return 1;
275 }
276 nlopt_init_genrand(time(NULL));
277 sdim = atoi(argv[1]);
278 s = nlopt_sobol_create(sdim);
279 n = atoi(argv[2]);
280 nlopt_sobol_skip(s, n, x);
281 for (j = 1; j <= n; ++j) {
282 nlopt_sobol_next01(s, x);
283 testint_sobol += testfunc(sdim, x);
284 if (j < 100) {
285 printf("x[%d]: %g", j, x[0]);
286 for (i = 1; i < sdim; ++i) printf(", %g", x[i]);
287 printf("\n");
288 }
289 for (i = 0; i < sdim; ++i) x[i] = nlopt_urand(0.,1.);
290 testint_rand += testfunc(sdim, x);
291 }
292 nlopt_sobol_destroy(s);
293 printf("Test integral = %g using Sobol, %g using pseudorandom.\n",
294 testint_sobol / n, testint_rand / n);
295 printf(" error = %g using Sobol, %g using pseudorandom.\n",
296 testint_sobol / n - 1, testint_rand / n - 1);
297 return 0;
298 }
299 #endif
300