From b23e6679cfce2f3d4266cb318c3b4afec6073449 Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Thu, 18 Feb 2016 18:13:02 -0500 Subject: add Nlopt version of sobol which allows higher dimensions --- src/sobolseq.c | 277 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 277 insertions(+) create mode 100644 src/sobolseq.c (limited to 'src/sobolseq.c') diff --git a/src/sobolseq.c b/src/sobolseq.c new file mode 100644 index 0000000..bdd310b --- /dev/null +++ b/src/sobolseq.c @@ -0,0 +1,277 @@ +/* Copyright (c) 2007 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* Generation of Sobol sequences in up to 1111 dimensions, based on the + algorithms described in: + P. Bratley and B. L. Fox, Algorithm 659, ACM Trans. + Math. Soft. 14 (1), 88-100 (1988), + as modified by: + S. Joe and F. Y. Kuo, ACM Trans. Math. Soft 29 (1), 49-57 (2003). + + Note that the code below was written without even looking at the + Fortran code from the TOMS paper, which is only semi-free (being + under the restrictive ACM copyright terms). Then I went to the + Fortran code and took out the table of primitive polynomials and + starting direction #'s ... since this is just a table of numbers + generated by a deterministic algorithm, it is not copyrightable. + (Obviously, the format of these tables then necessitated some + slight modifications to the code.) + + For the test integral of Joe and Kuo (see the main() program + below), I get exactly the same results for integrals up to 1111 + dimensions compared to the table of published numbers (to the 5 + published significant digits). + + This is not to say that the authors above should not be credited for + their clear description of the algorithm (and their tabulation of + the critical numbers). Please cite them. Just that I needed + a free/open-source implementation. */ + +#include +#include + +#include "nlopt-util.h" + +#include +#include + +typedef struct nlopt_soboldata_s { + unsigned sdim; /* dimension of sequence being generated */ + uint32_t *mdata; /* array of length 32 * sdim */ + uint32_t *m[32]; /* more convenient pointers to mdata, of direction #s */ + uint32_t *x; /* previous x = x_n, array of length sdim */ + unsigned *b; /* position of fixed point in x[i] is after bit b[i] */ + uint32_t n; /* number of x's generated so far */ +} soboldata; + +/* Return position (0, 1, ...) of rightmost (least-significant) zero bit in n. + * + * This code uses a 32-bit version of algorithm to find the rightmost + * one bit in Knuth, _The Art of Computer Programming_, volume 4A + * (draft fascicle), section 7.1.3, "Bitwise tricks and + * techniques." + * + * Assumes n has a zero bit, i.e. n < 2^32 - 1. + * + */ +static unsigned rightzero32(uint32_t n) +{ +#if defined(__GNUC__) && \ + ((__GNUC__ == 3 && __GNUC_MINOR__ >= 4) || __GNUC__ > 3) + return __builtin_ctz(~n); /* gcc builtin for version >= 3.4 */ +#else + const uint32_t a = 0x05f66a47; /* magic number, found by brute force */ + 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}; + n = ~n; /* change to rightmost-one problem */ + n = a * (n & (-n)); /* store in n to make sure mult. is 32 bits */ + return decode[n >> 27]; +#endif +} + +/* generate the next term x_{n+1} in the Sobol sequence, as an array + x[sdim] of numbers in (0,1). Returns 1 on success, 0 on failure + (if too many #'s generated) */ +static int sobol_gen(soboldata *sd, double *x) +{ + unsigned c, b, i, sdim; + + if (sd->n == 4294967295U) return 0; /* n == 2^32 - 1 ... we would + need to switch to a 64-bit version + to generate more terms. */ + c = rightzero32(sd->n++); + sdim = sd->sdim; + for (i = 0; i < sdim; ++i) { + b = sd->b[i]; + if (b >= c) { + sd->x[i] ^= sd->m[c][i] << (b - c); + x[i] = ((double) (sd->x[i])) / (1U << (b+1)); + } + else { + sd->x[i] = (sd->x[i] << (c - b)) ^ sd->m[c][i]; + sd->b[i] = c; + x[i] = ((double) (sd->x[i])) / (1U << (c+1)); + } + } + return 1; +} + +#include "soboldata.h" + +static int sobol_init(soboldata *sd, unsigned sdim) +{ + unsigned i,j; + + if (!sdim || sdim > MAXDIM) return 0; + + sd->mdata = (uint32_t *) malloc(sizeof(uint32_t) * (sdim * 32)); + if (!sd->mdata) return 0; + + for (j = 0; j < 32; ++j) { + sd->m[j] = sd->mdata + j * sdim; + sd->m[j][0] = 1; /* special-case Sobol sequence */ + } + for (i = 1; i < sdim; ++i) { + uint32_t a = sobol_a[i-1]; + unsigned d = 0, k; + + while (a) { + ++d; + a >>= 1; + } + d--; /* d is now degree of poly */ + + /* set initial values of m from table */ + for (j = 0; j < d; ++j) + sd->m[j][i] = sobol_minit[j][i-1]; + + /* fill in remaining values using recurrence */ + for (j = d; j < 32; ++j) { + a = sobol_a[i-1]; + sd->m[j][i] = sd->m[j - d][i]; + for (k = 0; k < d; ++k) { + sd->m[j][i] ^= ((a & 1) * sd->m[j-d+k][i]) << (d-k); + a >>= 1; + } + } + } + + sd->x = (uint32_t *) malloc(sizeof(uint32_t) * sdim); + if (!sd->x) { free(sd->mdata); return 0; } + + sd->b = (unsigned *) malloc(sizeof(unsigned) * sdim); + if (!sd->b) { free(sd->x); free(sd->mdata); return 0; } + + for (i = 0; i < sdim; ++i) { + sd->x[i] = 0; + sd->b[i] = 0; + } + + sd->n = 0; + sd->sdim = sdim; + + return 1; +} + +static void sobol_destroy(soboldata *sd) +{ + free(sd->mdata); + free(sd->x); + free(sd->b); +} + +/************************************************************************/ +/* NLopt API to Sobol sequence creation, which hides soboldata structure + behind an opaque pointer */ + +nlopt_sobol nlopt_sobol_create(unsigned sdim) +{ + nlopt_sobol s = (nlopt_sobol) malloc(sizeof(soboldata)); + if (!s) return NULL; + if (!sobol_init(s, sdim)) { free(s); return NULL; } + return s; +} + +extern void nlopt_sobol_destroy(nlopt_sobol s) +{ + if (s) { + sobol_destroy(s); + free(s); + } +} + +/* next vector x[sdim] in Sobol sequence, with each x[i] in (0,1) */ +void nlopt_sobol_next01(nlopt_sobol s, double *x) +{ + assert(sobol_gen(s, x)); +} + +/* next vector in Sobol sequence, scaled to (lb[i], ub[i]) interval */ +void nlopt_sobol_next(nlopt_sobol s, double *x, + const double *lb, const double *ub) +{ + unsigned i, sdim; + nlopt_sobol_next01(s, x); + for (sdim = s->sdim, i = 0; i < sdim; ++i) + x[i] = lb[i] + (ub[i] - lb[i]) * x[i]; +} + +/* if we know in advance how many points (n) we want to compute, then + adopt the suggestion of the Joe and Kuo paper, which in turn + is taken from Acworth et al (1998), of skipping a number of + points equal to the largest power of 2 smaller than n */ +void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x) +{ + if (s) { + unsigned k = 1; + while (k*2 < n) k *= 2; + while (k-- > 0) sobol_gen(s, x); + } +} + +/************************************************************************/ + +/* compile with -DSOBOLSEQ_TEST for test program */ +#ifdef SOBOLSEQ_TEST +#include +#include + +/* test integrand from Joe and Kuo paper ... integrates to 1 */ +static double testfunc(unsigned n, const double *x) +{ + double f = 1; + unsigned j; + for (j = 1; j <= n; ++j) { + double cj = pow((double) j, 0.3333333333333333333); + f *= (fabs(4*x[j-1] - 2) + cj) / (1 + cj); + } + return f; +} + +int main(int argc, char **argv) +{ + unsigned n, j, i, sdim; + static double x[MAXDIM]; + double testint_sobol = 0, testint_rand = 0; + nlopt_sobol s; + if (argc < 3) { + fprintf(stderr, "Usage: %s \n", argv[0]); + return 1; + } + sdim = atoi(argv[1]); + s = nlopt_sobol_create(sdim); + n = atoi(argv[2]); + nlopt_sobol_skip(s, n, x); + for (j = 1; j <= n; ++j) { + nlopt_sobol_next01(s, x); + testint_sobol += testfunc(sdim, x); + if (j < 100) { + printf("x[%d]: %g", j, x[0]); + for (i = 1; i < sdim; ++i) printf(", %g", x[i]); + printf("\n"); + } + } + nlopt_sobol_destroy(s); + printf("Test integral = %g using Sobol", testint_sobol / n); + printf(" error = %g using Sobol", testint_sobol / n - 1); + return 0; +} +#endif -- cgit v1.2.3-70-g09d2