src/util/sobolseq.c

   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  */
  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) */
 104 static int sobol_gen(soboldata * sd, double *x)
 105 {
 106     unsigned c, b, i, sdim;
 107
 108     if (sd->n == 4294967295U)
 109         return 0;               /* n == 2^32 - 1 ... we would
 110                                    need to switch to a 64-bit version
 111                                    to generate more terms. */
 112     c = rightzero32(sd->n++);
 113     sdim = sd->sdim;
 114     for (i = 0; i < sdim; ++i) {
 115         b = sd->b[i];
 116         if (b >= c) {
 117             sd->x[i] ^= sd->m[c][i] << (b - c);
 118             x[i] = ((double) (sd->x[i])) / (1U << (b + 1));
 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
 130 static int sobol_init(soboldata * sd, unsigned sdim)
 131 {
 132     unsigned i, j;
 133
 134     if (!sdim || sdim > MAXDIM)
 135         return 0;
 136
 137     sd->mdata = (uint32_t *) malloc(sizeof(uint32_t) * (sdim * 32));
 138     if (!sd->mdata)
 139         return 0;
 140
 141     for (j = 0; j < 32; ++j) {
 142         sd->m[j] = sd->mdata + j * sdim;
 143         sd->m[j][0] = 1;        /* special-case Sobol sequence */
 144     }
 145     for (i = 1; i < sdim; ++i) {
 146         uint32_t a = sobol_a[i - 1];
 147         unsigned d = 0, k;
 148
 149         while (a) {
 150             ++d;
 151             a >>= 1;
 152         }
 153         d--;                    /* d is now degree of poly */
 154
 155         /* set initial values of m from table */
 156         for (j = 0; j < d; ++j)
 157             sd->m[j][i] = sobol_minit[j][i - 1];
 158
 159         /* fill in remaining values using recurrence */
 160         for (j = d; j < 32; ++j) {
 161             a = sobol_a[i - 1];
 162             sd->m[j][i] = sd->m[j - d][i];
 163             for (k = 0; k < d; ++k) {
 164                 sd->m[j][i] ^= ((a & 1) * sd->m[j - d + k][i]) << (d - k);
 165                 a >>= 1;
 166             }
 167         }
 168     }
 169
 170     sd->x = (uint32_t *) malloc(sizeof(uint32_t) * sdim);
 171     if (!sd->x) {
 172         free(sd->mdata);
 173         return 0;
 174     }
 175
 176     sd->b = (unsigned *) malloc(sizeof(unsigned) * sdim);
 177     if (!sd->b) {
 178         free(sd->x);
 179         free(sd->mdata);
 180         return 0;
 181     }
 182
 183     for (i = 0; i < sdim; ++i) {
 184         sd->x[i] = 0;
 185         sd->b[i] = 0;
 186     }
 187
 188     sd->n = 0;
 189     sd->sdim = sdim;
 190
 191     return 1;
 192 }
 193
 194 static void sobol_destroy(soboldata * sd)
 195 {
 196     free(sd->mdata);
 197     free(sd->x);
 198     free(sd->b);
 199 }
 200
 201 /************************************************************************/
 202 /* NLopt API to Sobol sequence creation, which hides soboldata structure
 203    behind an opaque pointer */
 204
 205 nlopt_sobol nlopt_sobol_create(unsigned sdim)
 206 {
 207     nlopt_sobol s = (nlopt_sobol) malloc(sizeof(soboldata));
 208     if (!s)
 209         return NULL;
 210     if (!sobol_init(s, sdim)) {
 211         free(s);
 212         return NULL;
 213     }
 214     return s;
 215 }
 216
 217 extern void nlopt_sobol_destroy(nlopt_sobol s)
 218 {
 219     if (s) {
 220         sobol_destroy(s);
 221         free(s);
 222     }
 223 }
 224
 225 /* next vector x[sdim] in Sobol sequence, with each x[i] in (0,1) */
 226 void nlopt_sobol_next01(nlopt_sobol s, double *x)
 227 {
 228     if (!sobol_gen(s, x)) {
 229         /* fall back on pseudo random numbers in the unlikely event
 230            that we exceed 2^32-1 points */
 231         unsigned i;
 232         for (i = 0; i < s->sdim; ++i)
 233             x[i] = nlopt_urand(0.0, 1.0);
 234     }
 235 }
 236
 237 /* next vector in Sobol sequence, scaled to (lb[i], ub[i]) interval */
 238 void nlopt_sobol_next(nlopt_sobol s, double *x, const double *lb, const double *ub)
 239 {
 240     unsigned i, sdim;
 241     nlopt_sobol_next01(s, x);
 242     for (sdim = s->sdim, i = 0; i < sdim; ++i)
 243         x[i] = lb[i] + (ub[i] - lb[i]) * x[i];
 244 }
 245
 246 /* if we know in advance how many points (n) we want to compute, then
 247    adopt the suggestion of the Joe and Kuo paper, which in turn
 248    is taken from Acworth et al (1998), of skipping a number of
 249    points equal to the largest power of 2 smaller than n */
 250 void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x)
 251 {
 252     if (s) {
 253         unsigned k = 1;
 254         while (k * 2 < n)
 255             k *= 2;
 256         while (k-- > 0)
 257             sobol_gen(s, x);
 258     }
 259 }