sobolseq.c File Reference

#include "internal.h"
#include "soboldata.h"

Include dependency graph for sobolseq.c:

Go to the source code of this file.

Data Structures
struct	nlopt_soboldata_s

Typedefs
typedef Int32	uint32_t
typedef struct nlopt_soboldata_s *	nlopt_sobol
typedef struct nlopt_soboldata_s	soboldata

Functions
static nlopt_sobol	nlopt_sobol_create (unsigned sdim)
static void	nlopt_sobol_destroy (nlopt_sobol s)
static void	nlopt_sobol_skip (nlopt_sobol s, unsigned n, double *x)
static unsigned	rightzero32 (uint32_t n)
static int	sobol_gen (soboldata *sd, double *x)
static int	sobol_init (soboldata *sd, unsigned sdim)
static void	sobol_destroy (soboldata *sd)
SEXP	sobol_sequence (SEXP dim, SEXP length)

Typedef Documentation

nlopt_sobol

typedef struct nlopt_soboldata_s* nlopt_sobol

Definition at line 60 of file sobolseq.c.

soboldata

typedef struct nlopt_soboldata_s soboldata

uint32_t

typedef Int32 uint32_t

Definition at line 59 of file sobolseq.c.

Function Documentation

nlopt_sobol_create()

static nlopt_sobol nlopt_sobol_create ( unsigned  sdim )

static

Definition at line 188 of file sobolseq.c.

{
  nlopt_sobol s = (nlopt_sobol) R_Calloc(1,soboldata);
  if (!s) return 0;
  if (!sobol_init(s, sdim)) { R_Free(s); return NULL; }
  return s;
}

Here is the call graph for this function:

Here is the caller graph for this function:

nlopt_sobol_destroy()

static void nlopt_sobol_destroy ( nlopt_sobol  s )

static

Definition at line 196 of file sobolseq.c.

{
  if (s) {
    sobol_destroy(s);
    R_Free(s);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

nlopt_sobol_skip()

static void nlopt_sobol_skip ( nlopt_sobol  s, unsigned  n, double *  x  )

static

Definition at line 225 of file sobolseq.c.

{
  if (s) {
    unsigned k = 1;
    while (k*2 < n) k *= 2;
    while (k-- > 0) sobol_gen(s, x);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

rightzero32()

static unsigned rightzero32 ( uint32_t  n )

static

Definition at line 84 of file sobolseq.c.

{
  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];
}

Here is the caller graph for this function:

sobol_destroy()

static void sobol_destroy ( soboldata *  sd )

static

Definition at line 177 of file sobolseq.c.

{
  R_Free(sd->mdata);
  R_Free(sd->x);
  R_Free(sd->b);
}

Here is the caller graph for this function:

sobol_gen()

static int sobol_gen ( soboldata *  sd, double *  x  )

static

Definition at line 96 of file sobolseq.c.

{
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

sobol_init()

static int sobol_init ( soboldata *  sd, unsigned  sdim  )

static

Definition at line 122 of file sobolseq.c.

{
  unsigned i,j;
 
  if (!sdim || sdim > MAXDIM) return 0;
 
  sd->mdata = (uint32_t *) R_Calloc(sdim*32,uint32_t);
  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 *) R_Calloc(sdim,uint32_t);
  if (!sd->x) { R_Free(sd->mdata); return 0; }
 
  sd->b = (unsigned *) R_Calloc(sdim,unsigned);
  if (!sd->b) { R_Free(sd->x); R_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;
}

Here is the caller graph for this function:

sobol_sequence()

SEXP sobol_sequence	(	SEXP	dim,
		SEXP	length
	)

Definition at line 234 of file sobolseq.c.

{
  SEXP data;
  unsigned int d = (unsigned int) INTEGER(dim)[0];
  unsigned int n = (unsigned int) INTEGER(length)[0];
  double *dp;
  unsigned int k;
  nlopt_sobol s = nlopt_sobol_create((unsigned int) d);
  if (s==0) err("dimension is too high");
  PROTECT(data = allocMatrix(REALSXP,d,n)); dp = REAL(data);
  nlopt_sobol_skip(s,n,dp);
  for (k = 1; k < n; k++) sobol_gen(s,dp+k*d);
  nlopt_sobol_destroy(s);
  UNPROTECT(1);
  return(data);
}

Here is the call graph for this function: