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) 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
 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
 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
 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
 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) */
 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 */
 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 */
 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 }