1 /* Copyright (c) 2007 Massachusetts Institute of Technology
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:
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
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.
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),
28 S. Joe and F. Y. Kuo, ACM Trans. Math. Soft 29 (1), 49-57 (2003).
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.)
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).
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. */
52 #include "nlopt-util.h"
54 #if defined(HAVE_STDINT_H)
59 # if SIZEOF_UNSIGNED_LONG == 4
60 typedef unsigned long uint32_t;
61 # elif SIZEOF_UNSIGNED_INT == 4
62 typedef unsigned int uint32_t;
64 # error No 32-bit unsigned integer type
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 */
77 /* Return position (0, 1, ...) of rightmost (least-significant) zero bit in n.
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
84 * Assumes n has a zero bit, i.e. n < 2^32 - 1.
87 static unsigned rightzero32(uint32_t n)
89 #if defined(__GNUC__) && \
90 ((__GNUC__ == 3 && __GNUC_MINOR__ >= 4) || __GNUC__ > 3)
91 return __builtin_ctz(~n); /* gcc builtin for version >= 3.4 */
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];
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) */
104 static int sobol_gen(soboldata *sd, double *x)
106 unsigned c, b, i, sdim;
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++);
113 for (i = 0; i < sdim; ++i) {
116 sd->x[i] ^= sd->m[c][i] << (b - c);
117 x[i] = ((double) (sd->x[i])) / (1U << (b+1));
120 sd->x[i] = (sd->x[i] << (c - b)) ^ sd->m[c][i];
122 x[i] = ((double) (sd->x[i])) / (1U << (c+1));
128 #include "soboldata.h"
130 static int sobol_init(soboldata *sd, unsigned sdim)
134 if (!sdim || sdim > MAXDIM) return 0;
136 sd->mdata = (uint32_t *) malloc(sizeof(uint32_t) * (sdim * 32));
137 if (!sd->mdata) return 0;
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 */
143 for (i = 1; i < sdim; ++i) {
144 uint32_t a = sobol_a[i-1];
151 d--; /* d is now degree of poly */
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];
157 /* fill in remaining values using recurrence */
158 for (j = d; j < 32; ++j) {
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);
168 sd->x = (uint32_t *) malloc(sizeof(uint32_t) * sdim);
169 if (!sd->x) { free(sd->mdata); return 0; }
171 sd->b = (unsigned *) malloc(sizeof(unsigned) * sdim);
172 if (!sd->b) { free(sd->x); free(sd->mdata); return 0; }
174 for (i = 0; i < sdim; ++i) {
185 static void sobol_destroy(soboldata *sd)
192 /************************************************************************/
193 /* NLopt API to Sobol sequence creation, which hides soboldata structure
194 behind an opaque pointer */
196 nlopt_sobol nlopt_sobol_create(unsigned sdim)
198 nlopt_sobol s = (nlopt_sobol) malloc(sizeof(soboldata));
200 if (!sobol_init(s, sdim)) { free(s); return NULL; }
204 extern void nlopt_sobol_destroy(nlopt_sobol s)
212 /* next vector x[sdim] in Sobol sequence, with each x[i] in (0,1) */
213 void nlopt_sobol_next01(nlopt_sobol s, double *x)
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 */
219 for (i = 0; i < s->sdim; ++i)
220 x[i] = nlopt_urand(0.0,1.0);
224 /* next vector in Sobol sequence, scaled to (lb[i], ub[i]) interval */
225 void nlopt_sobol_next(nlopt_sobol s, double *x,
226 const double *lb, const double *ub)
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];
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 */
238 void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x)
242 while (k*2 < n) k *= 2;
243 while (k-- > 0) sobol_gen(s, x);
247 /************************************************************************/
249 /* compile with -DSOBOLSEQ_TEST for test program */
254 /* test integrand from Joe and Kuo paper ... integrates to 1 */
255 static double testfunc(unsigned n, const double *x)
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);
266 int main(int argc, char **argv)
268 unsigned n, j, i, sdim;
269 static double x[MAXDIM];
270 double testint_sobol = 0, testint_rand = 0;
273 fprintf(stderr, "Usage: %s <sdim> <ngen>\n", argv[0]);
276 nlopt_init_genrand(time(NULL));
277 sdim = atoi(argv[1]);
278 s = nlopt_sobol_create(sdim);
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);
285 printf("x[%d]: %g", j, x[0]);
286 for (i = 1; i < sdim; ++i) printf(", %g", x[i]);
289 for (i = 0; i < sdim; ++i) x[i] = nlopt_urand(0.,1.);
290 testint_rand += testfunc(sdim, x);
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);