path: root/src/divonne/Sample.c

/*
	Sample.c
		most of what is related to sampling
		this file is part of Divonne
		last modified 7 Nov 11 th
*/


#define MARKMASK 0xfffffff
#define Marked(x) ((x) & ~MARKMASK)
#define Unmark(x) ((x) & MARKMASK)

#define EXTRAPOLATE_EPS (.25*t->border.lower)
/*#define EXTRAPOLATE_EPS 0x1p-26*/

/*********************************************************************/

static inline void SamplesIni(Samples *samples)
{
  samples->x = NULL;
}

/*********************************************************************/

static inline bool SamplesIniQ(cSamples *samples)
{
  return samples->x == NULL;
}

/*********************************************************************/

static inline void SamplesFree(cSamples *samples)
{
  free(samples->x);
}

/*********************************************************************/

static void SampleSobol(This *t, ccount iregion)
{
  SAMPLERDEFS;
  real avg[NCOMP], norm;
  number i;
  count dim, comp;

  for( i = 0; i < n; ++i ) {
    t->rng.getrandom(t, x);
    for( dim = 0; dim < t->ndim; ++x, ++dim )
      *x = b[dim].lower + *x*(b[dim].upper - b[dim].lower);
  }

  DoSample(t, n, samples->x, f, t->ndim);

  ResCopy(avg, f);
  f += t->ncomp;
  for( i = 2; i < n; ++i )
    for( comp = 0; comp < t->ncomp; ++comp )
      avg[comp] += *f++;

  norm = region->vol/samples->neff;
  for( comp = 0; comp < t->ncomp; ++comp ) {
    r[comp].avg = norm*avg[comp];
    r[comp].err = 0;
  }
}

/*********************************************************************/

static void SampleKorobov(This *t, ccount iregion)
{
  SAMPLERDEFS;
  real *xlast = x + t->ndim, *flast = f + t->ncomp;
  real avg[NCOMP], norm;
  cnumber neff = samples->neff;
  number nextra = 0, i;
  real dist = 0;
  count dim, comp;

  for( i = 1; i < n; ++i ) {
    number c = i;
    for( dim = 0; dim < t->ndim; ++dim ) {
      creal dx = abs(2*c - neff)/(real)neff;
      *xlast++ = b[dim].lower + dx*(b[dim].upper - b[dim].lower);
      c = c*samples->coeff % neff;
    }
  }

  for( dim = 0; dim < t->ndim; ++dim ) {
    creal dx = (x[dim] = b[dim].upper) - t->border.upper;
    if( dx > 0 ) dist += Sq(dx);
  }

  if( dist > 0 ) {
    dist = sqrt(dist)/EXTRAPOLATE_EPS;
    for( dim = 0; dim < t->ndim; ++dim ) {
      real x2 = x[dim], dx = x2 - t->border.upper;
      if( dx > 0 ) {
        x[dim] = t->border.upper;
        x2 = t->border.upper - dx/dist;
      }
      xlast[dim] = x2;
    }
    nextra = 1;
  }

  DoSample(t, n + nextra, x, f, t->ndim);

  ResCopy(avg, flast);
  flast += t->ncomp;
  for( i = 2; i < n; ++i )
    for( comp = 0; comp < t->ncomp; ++comp )
      avg[comp] += *flast++;

  if( nextra ) {
    for( comp = 0; comp < t->ncomp; ++comp )
      f[comp] += dist*(f[comp] - flast[comp]);
    for( dim = 0; dim < t->ndim; ++dim )
      x[dim] = b[dim].upper;
  }

  norm = region->vol/samples->neff;
  for( comp = 0; comp < t->ncomp; ++comp ) {
    r[comp].avg = norm*(avg[comp] + avg[comp] + f[comp]);
    r[comp].err = 0;
  }
}

/*********************************************************************/

#define IsSobol(k) NegQ(k)
#define IsRule(k, d) (k == 9 || k == 7 || (k == 11 && d == 3) || (k == 13 && d == 2))

/* The following coding is used for key1, key2, key3:
             0 = for key1, key2: use default,
                 for key3: do nothing,
             1 = for key3: split region again,
             7 = degree-7 cubature rule,
             9 = degree-9 cubature rule,
            11 = degree-11 cubature rule (only in 3 dims),
            13 = degree-13 cubature rule (only in 2 dims),
     -inf..-40 = absolute # of points, Sobol numbers,
       -39..-1 = multiplicator,        Sobol numbers,
         1..39 = multiplicator,        Korobov numbers,
       40..inf = absolute # of points, Korobov numbers.  */

static count SamplesLookup(This *t, Samples *samples, cint key,
  cnumber nwant, cnumber nmax, number nmin)
{
  number n;

  if( key == 13 && t->ndim == 2 ) {
    if( RuleIniQ(&t->rule13) ) Rule13Alloc(t);
    samples->rule = &t->rule13;
    samples->n = n = nmin = t->rule13.n;
    samples->sampler = SampleRule;
  }
  else if( key == 11 && t->ndim == 3 ) {
    if( RuleIniQ(&t->rule11) ) Rule11Alloc(t);
    samples->rule = &t->rule11;
    samples->n = n = nmin = t->rule11.n;
    samples->sampler = SampleRule;
  }
  else if( key == 9 ) {
    if( RuleIniQ(&t->rule9) ) Rule9Alloc(t);
    samples->rule = &t->rule9;
    samples->n = n = nmin = t->rule9.n;
    samples->sampler = SampleRule;
  }
  else if( key == 7 ) {
    if( RuleIniQ(&t->rule7) ) Rule7Alloc(t);
    samples->rule = &t->rule7;
    samples->n = n = nmin = t->rule7.n;
    samples->sampler = SampleRule;
  }
  else {
    n = Abs1(key);
    if( n < 40 ) n *= nwant;
    samples->sampler = (key < 0) ? SampleSobol :
      (n = n/2 + 1, SampleKorobov);
    samples->n = IMin(n, nmax);
  }

  samples->neff = samples->n;

  return IDim(n - nmax) | Marked(nmax - nmin);
}

/*********************************************************************/

static void SamplesAlloc(cThis *t, Samples *samples)
{
#define FIRST -INT_MAX
#define MarkLast(x) (x | Marked(INT_MAX))

#include "KorobovCoeff.c"

  number nx, nf;

  if( samples->sampler == SampleKorobov ) {
    enum { max = Elements(prime) - 2 };
    cint n = IMin(2*samples->n - 1, MAXPRIME);
    int i = Hash(n), p;
    count shift = 2 + NegQ(n - 1000);

    while( i = IMin(IDim(i), max),
           n > (p = prime[i + 1]) || n <= prime[i] ) {
      cint d = (n - Unmark(p)) >> ++shift;
      i += Min1(d);
    }

    samples->coeff = coeff[i][t->ndim - KOROBOV_MINDIM];
    samples->neff = p = Unmark(p);
    samples->n = p/2 + 1;
  }

  nx = t->ndim*(samples->n + 1);		/* need 1 for extrapolation */
  nf = t->ncomp*(samples->n + 1);

  Alloc(samples->x, nx + nf + t->ncomp + t->ncomp);
  samples->f = samples->x + nx;
}

/*********************************************************************/

static real Sample(This *t, creal *x0)
{
  real xtmp[2*NDIM], ftmp[2*NCOMP], *xlast = xtmp, f;
  real dist = 0;
  count dim, comp;
  number n = 1;

  for( dim = 0; dim < t->ndim; ++dim ) {
    creal x1 = *xlast++ = Min(Max(*x0++, 0.), 1.);
    real dx;
    if( (dx = x1 - t->border.lower) < 0 ||
        (dx = x1 - t->border.upper) > 0 ) dist += Sq(dx);
  }

  if( dist > 0 ) {
    dist = sqrt(dist)/EXTRAPOLATE_EPS;
    for( dim = 0; dim < t->ndim; ++dim ) {
      real x2 = xtmp[dim], dx, b;
      if( (dx = x2 - (b = t->border.lower)) < 0 ||
          (dx = x2 - (b = t->border.upper)) > 0 ) {
        xtmp[dim] = b;
        x2 = b - dx/dist;
      }
      *xlast++ = x2;
    }
    n = 2;
  }

  DoSample(t, n, xtmp, ftmp, t->ndim);

  comp = Untag(t->selectedcomp);
  f = ftmp[comp];
  if( n > 1 ) f += dist*(f - ftmp[comp + t->ncomp]);

  return Sign(t->selectedcomp)*f;
}