summaryrefslogtreecommitdiff
path: root/src/suave
diff options
context:
space:
mode:
Diffstat (limited to 'src/suave')
-rw-r--r--src/suave/Fluct.c68
-rw-r--r--src/suave/Grid.c137
-rw-r--r--src/suave/Integrate.c278
-rw-r--r--src/suave/Sample.c171
-rw-r--r--src/suave/Suave.c91
-rw-r--r--src/suave/Suave.tm290
-rw-r--r--src/suave/common.c33
-rw-r--r--src/suave/decl.h70
8 files changed, 1138 insertions, 0 deletions
diff --git a/src/suave/Fluct.c b/src/suave/Fluct.c
new file mode 100644
index 0000000..2755489
--- /dev/null
+++ b/src/suave/Fluct.c
@@ -0,0 +1,68 @@
+/*
+ Fluct.c
+ compute the fluctuation in the left and right half
+ this file is part of Suave
+ last modified 23 Aug 11 th
+*/
+
+
+#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL)
+
+typedef long double realx;
+#define XDBL_MAX_EXP LDBL_MAX_EXP
+#define XDBL_MAX LDBL_MAX
+#define powx powl
+#define ldexpx ldexpl
+
+#else
+
+typedef double realx;
+#define XDBL_MAX_EXP DBL_MAX_EXP
+#define XDBL_MAX DBL_MAX
+#define powx pow
+#define ldexpx ldexp
+
+#endif
+
+typedef const realx crealx;
+
+typedef struct {
+ realx fluct;
+ number n;
+} Var;
+
+/*********************************************************************/
+
+static void Fluct(cThis *t, Var *var,
+ cBounds *b, creal *w, number n, ccount comp, creal avg, creal err)
+{
+ creal *x = w + n;
+ creal *f = x + n*t->ndim + comp;
+ count nvar = 2*t->ndim;
+ creal norm = 1/(err*Max(fabs(avg), err));
+ creal flat = 2/3./t->flatness;
+ crealx max = ldexpx(1., (int)((XDBL_MAX_EXP - 2)/t->flatness));
+
+ Clear(var, nvar);
+
+ while( n-- ) {
+ count dim;
+ crealx arg = 1 + fabs(*w++)*Sq(*f - avg)*norm;
+ crealx ft = powx(arg < max ? arg : max, t->flatness);
+
+ f += t->ncomp;
+
+ for( dim = 0; dim < t->ndim; ++dim ) {
+ Var *v = &var[2*dim + (*x++ >= b[dim].mid)];
+ crealx f = v->fluct + ft;
+ v->fluct = (f > XDBL_MAX/2) ? XDBL_MAX/2 : f;
+ ++v->n;
+ }
+ }
+
+ while( nvar-- ) {
+ var->fluct = powx(var->fluct, flat);
+ ++var;
+ }
+}
+
diff --git a/src/suave/Grid.c b/src/suave/Grid.c
new file mode 100644
index 0000000..013f086
--- /dev/null
+++ b/src/suave/Grid.c
@@ -0,0 +1,137 @@
+/*
+ Grid.c
+ utility functions for the Vegas grid
+ this file is part of Suave
+ last modified 1 Jun 10 th
+*/
+
+
+static void RefineGrid(cThis *t, Grid grid, Grid margsum)
+{
+ real avgperbin, thisbin, newcur, delta;
+ Grid imp, newgrid;
+ int bin, newbin;
+
+ /* smooth the f^2 value stored for each bin */
+ real prev = margsum[0];
+ real cur = margsum[1];
+ real norm = margsum[0] = .5*(prev + cur);
+ for( bin = 1; bin < NBINS - 1; ++bin ) {
+ creal s = prev + cur;
+ prev = cur;
+ cur = margsum[bin + 1];
+ norm += margsum[bin] = (s + cur)/3.;
+ }
+ norm += margsum[NBINS - 1] = .5*(prev + cur);
+
+ if( norm == 0 ) return;
+ norm = 1/norm;
+
+ /* compute the importance function for each bin */
+ avgperbin = 0;
+ for( bin = 0; bin < NBINS; ++bin ) {
+ real impfun = 0;
+ if( margsum[bin] > 0 ) {
+ creal r = margsum[bin]*norm;
+ avgperbin += impfun = pow((r - 1)/log(r), 1.5);
+ }
+ imp[bin] = impfun;
+ }
+ avgperbin /= NBINS;
+
+ /* redefine the size of each bin */
+ cur = newcur = 0;
+ thisbin = 0;
+ bin = -1;
+ for( newbin = 0; newbin < NBINS - 1; ++newbin ) {
+ while( thisbin < avgperbin ) {
+ thisbin += imp[++bin];
+ prev = cur;
+ cur = grid[bin];
+ }
+ thisbin -= avgperbin;
+ delta = (cur - prev)*thisbin;
+ newgrid[newbin] = SHARPEDGES ?
+ cur - delta/imp[bin] :
+ (newcur = Max(newcur + 0x1p-48,
+ cur - 2*delta/(imp[bin] + imp[IDim(bin - 1)])));
+ }
+ Copy(grid, newgrid, NBINS - 1);
+ grid[NBINS - 1] = 1;
+}
+
+/*********************************************************************/
+
+static void Reweight(cThis *t, Bounds *b,
+ creal *w, creal *f, creal *lastf, cResult *total)
+{
+ Grid margsum[NDIM];
+ real scale[NCOMP];
+ cbin_t *bin = (cbin_t *)lastf;
+ count dim, comp;
+
+ if( t->ncomp == 1 ) scale[0] = 1;
+ else {
+ for( comp = 0; comp < t->ncomp; ++comp )
+ scale[comp] = (total[comp].avg == 0) ? 0 : 1/total[comp].avg;
+ }
+
+ Zap(margsum);
+
+ while( f < lastf ) {
+ real fsq = 0;
+ for( comp = 0; comp < t->ncomp; ++comp )
+ fsq += Sq(*f++*scale[comp]);
+ fsq *= Sq(*w++);
+ if( fsq != 0 )
+ for( dim = 0; dim < t->ndim; ++dim )
+ margsum[dim][bin[dim]] += fsq;
+ bin += t->ndim;
+ }
+
+ for( dim = 0; dim < t->ndim; ++dim )
+ RefineGrid(t, b[dim].grid, margsum[dim]);
+}
+
+/*********************************************************************/
+
+static void StretchGrid(cGrid grid, Grid gridL, Grid gridR)
+{
+ real prev = 0, cur, step, x;
+ count bin = 0;
+
+ while( bin < NBINS ) {
+ cur = grid[bin++];
+ if( cur >= .5 ) break;
+ prev = cur;
+ }
+
+ step = (bin - (cur - .5)/(cur - prev))/NBINS;
+
+ prev = x = 0;
+ cur = *grid;
+
+ for( bin = 0; bin < NBINS; ++bin ) {
+ x += step;
+ if( x > 1 ) {
+ --x;
+ prev = cur;
+ cur = *++grid;
+ }
+ gridL[bin] = 2*(prev + (cur - prev)*x);
+ }
+
+ step = 1 - step;
+ for( bin = 0; bin < NBINS - 1; ++bin ) {
+ x += step;
+ if( x > 1 ) {
+ --x;
+ prev = cur;
+ cur = *++grid;
+ }
+ gridR[bin] = 2*(prev + (cur - prev)*x) - 1;
+ }
+ gridR[NBINS - 1] = 1;
+}
+
+
diff --git a/src/suave/Integrate.c b/src/suave/Integrate.c
new file mode 100644
index 0000000..bb81b70
--- /dev/null
+++ b/src/suave/Integrate.c
@@ -0,0 +1,278 @@
+/*
+ Integrate.c
+ integrate over the unit hypercube
+ this file is part of Suave
+ last modified 15 Feb 11 th
+*/
+
+
+static int Integrate(This *t, real *integral, real *error, real *prob)
+{
+ TYPEDEFREGION;
+
+ count dim, comp, df;
+ int fail;
+ Result totals[NCOMP];
+ Region *anchor = NULL, *region = NULL;
+
+ if( VERBOSE > 1 ) {
+ char s[256];
+ sprintf(s, "Suave input parameters:\n"
+ " ndim " COUNT "\n ncomp " COUNT "\n"
+ " epsrel " REAL "\n epsabs " REAL "\n"
+ " flags %d\n seed %d\n"
+ " mineval " NUMBER "\n maxeval " NUMBER "\n"
+ " nnew " NUMBER "\n flatness " REAL,
+ t->ndim, t->ncomp,
+ t->epsrel, t->epsabs,
+ t->flags, t->seed,
+ t->mineval, t->maxeval,
+ t->nnew, t->flatness);
+ Print(s);
+ }
+
+ if( BadComponent(t) ) return -2;
+ if( BadDimension(t) ) return -1;
+
+ if( (fail = setjmp(t->abort)) ) goto abort;
+
+ t->epsabs = Max(t->epsabs, NOTZERO);
+ IniRandom(t);
+
+ RegionAlloc(anchor, t->nnew, t->nnew);
+ anchor->next = NULL;
+ anchor->div = 0;
+
+ for( dim = 0; dim < t->ndim; ++dim ) {
+ Bounds *b = &anchor->bounds[dim];
+ b->lower = 0;
+ b->upper = 1;
+ b->mid = .5;
+
+ if( dim == 0 ) {
+ count bin;
+ /* define the initial distribution of bins */
+ for( bin = 0; bin < NBINS; ++bin )
+ b->grid[bin] = (bin + 1)/(real)NBINS;
+ }
+ else Copy(b->grid, anchor->bounds[0].grid, NBINS);
+ }
+
+ Sample(t, t->nnew, anchor, anchor->w,
+ anchor->w + t->nnew,
+ anchor->w + t->nnew + t->ndim*t->nnew);
+ df = anchor->df;
+ ResCopy(totals, anchor->result);
+
+ for( t->nregions = 1; ; ++t->nregions ) {
+ Var var[NDIM][2], *vLR;
+ real maxratio, maxerr, minfluct, bias, mid;
+ Region *regionL, *regionR, *reg, **parent, **par;
+ Bounds *bounds, *boundsL, *boundsR;
+ count maxcomp, bisectdim;
+ number n, nL, nR, nnewL, nnewR;
+ real *w, *wL, *wR, *x, *xL, *xR, *f, *fL, *fR, *wlast, *flast;
+
+ if( VERBOSE ) {
+ char s[128 + 128*NCOMP], *p = s;
+
+ p += sprintf(p, "\n"
+ "Iteration " COUNT ": " NUMBER " integrand evaluations so far",
+ t->nregions, t->neval);
+
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ cResult *tot = &totals[comp];
+ p += sprintf(p, "\n[" COUNT "] "
+ REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)",
+ comp + 1, tot->avg, tot->err, tot->chisq, df);
+ }
+
+ Print(s);
+ }
+
+ maxratio = -INFTY;
+ maxcomp = 0;
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ creal ratio = totals[comp].err/MaxErr(totals[comp].avg);
+ if( ratio > maxratio ) {
+ maxratio = ratio;
+ maxcomp = comp;
+ }
+ }
+
+ if( maxratio <= 1 && t->neval >= t->mineval ) {
+ fail = 0;
+ break;
+ }
+
+ if( t->neval >= t->maxeval ) break;
+
+ maxerr = -INFTY;
+ parent = &anchor;
+ region = anchor;
+ for( par = &anchor; (reg = *par); par = &reg->next ) {
+ creal err = reg->result[maxcomp].err;
+ if( err > maxerr ) {
+ maxerr = err;
+ parent = par;
+ region = reg;
+ }
+ }
+
+ Fluct(t, var[0],
+ region->bounds, region->w, region->n, maxcomp,
+ region->result[maxcomp].avg, Max(maxerr, t->epsabs));
+
+ bias = (t->epsrel < 1e-50) ? 2 :
+ Max(pow(2., -(real)region->div/t->ndim)/t->epsrel, 2.);
+ minfluct = INFTY;
+ bisectdim = 0;
+ for( dim = 0; dim < t->ndim; ++dim ) {
+ cBounds *b = &region->bounds[dim];
+ creal fluct = (var[dim][0].fluct + var[dim][1].fluct)*
+ (bias - b->upper + b->lower);
+ if( fluct < minfluct ) {
+ minfluct = fluct;
+ bisectdim = dim;
+ }
+ }
+
+ vLR = var[bisectdim];
+ minfluct = vLR[0].fluct + vLR[1].fluct;
+ nnewL = IMax(
+ (minfluct == 0) ? t->nnew/2 : (count)(vLR[0].fluct/minfluct*t->nnew),
+ MINSAMPLES );
+ nL = vLR[0].n + nnewL;
+ nnewR = IMax(t->nnew - nnewL, MINSAMPLES);
+ nR = vLR[1].n + nnewR;
+
+ RegionAlloc(regionL, nL, nnewL);
+ RegionAlloc(regionR, nR, nnewR);
+
+ *parent = regionL;
+ regionL->next = regionR;
+ regionR->next = region->next;
+ regionL->div = regionR->div = region->div + 1;
+
+ bounds = &region->bounds[bisectdim];
+ mid = bounds->mid;
+ n = region->n;
+ w = wlast = region->w; x = w + n; f = flast = x + n*t->ndim;
+ wL = regionL->w; xL = wL + nL; fL = xL + nL*t->ndim;
+ wR = regionR->w; xR = wR + nR; fR = xR + nR*t->ndim;
+
+ while( n-- ) {
+ cbool final = (*w < 0);
+ if( x[bisectdim] < mid ) {
+ if( final && wR > regionR->w ) *(wR - 1) = -fabs(*(wR - 1));
+ *wL++ = *w++;
+ VecCopy(xL, x);
+ xL += t->ndim;
+ ResCopy(fL, f);
+ fL += t->ncomp;
+ }
+ else {
+ if( final && wL > regionL->w ) *(wL - 1) = -fabs(*(wL - 1));
+ *wR++ = *w++;
+ VecCopy(xR, x);
+ xR += t->ndim;
+ ResCopy(fR, f);
+ fR += t->ncomp;
+ }
+ x += t->ndim;
+ f += t->ncomp;
+ if( n && final ) wlast = w, flast = f;
+ }
+
+ Reweight(t, region->bounds, wlast, flast, f, totals);
+ VecCopy(regionL->bounds, region->bounds);
+ VecCopy(regionR->bounds, region->bounds);
+
+ boundsL = &regionL->bounds[bisectdim];
+ boundsR = &regionR->bounds[bisectdim];
+ boundsL->mid = .5*(boundsL->lower + (boundsL->upper = mid));
+ boundsR->mid = .5*((boundsR->lower = mid) + boundsR->upper);
+
+ StretchGrid(bounds->grid, boundsL->grid, boundsR->grid);
+
+ Sample(t, nnewL, regionL, wL, xL, fL);
+ Sample(t, nnewR, regionR, wR, xR, fR);
+
+ df += regionL->df + regionR->df - region->df;
+
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ cResult *r = &region->result[comp];
+ Result *rL = &regionL->result[comp];
+ Result *rR = &regionR->result[comp];
+ Result *tot = &totals[comp];
+ real diff, sigsq;
+
+ tot->avg += diff = rL->avg + rR->avg - r->avg;
+
+ diff = Sq(.25*diff);
+ sigsq = rL->sigsq + rR->sigsq;
+ if( sigsq > 0 ) {
+ creal c = Sq(1 + sqrt(diff/sigsq));
+ rL->sigsq *= c;
+ rR->sigsq *= c;
+ }
+ rL->err = sqrt(rL->sigsq += diff);
+ rR->err = sqrt(rR->sigsq += diff);
+
+ tot->sigsq += rL->sigsq + rR->sigsq - r->sigsq;
+ tot->err = sqrt(tot->sigsq);
+
+ tot->chisq += rL->chisq + rR->chisq - r->chisq;
+ }
+
+ free(region);
+ region = NULL;
+ }
+
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ cResult *tot = &totals[comp];
+ integral[comp] = tot->avg;
+ error[comp] = tot->err;
+ prob[comp] = ChiSquare(tot->chisq, df);
+ }
+
+#ifdef MLVERSION
+ if( REGIONS ) {
+ MLPutFunction(stdlink, "List", 2);
+ MLPutFunction(stdlink, "List", t->nregions);
+ for( region = anchor; region; region = region->next ) {
+ real lower[NDIM], upper[NDIM];
+
+ for( dim = 0; dim < t->ndim; ++dim ) {
+ cBounds *b = &region->bounds[dim];
+ lower[dim] = b->lower;
+ upper[dim] = b->upper;
+ }
+
+ MLPutFunction(stdlink, "Cuba`Suave`region", 4);
+ MLPutRealList(stdlink, lower, t->ndim);
+ MLPutRealList(stdlink, upper, t->ndim);
+
+ MLPutFunction(stdlink, "List", t->ncomp);
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ cResult *r = &region->result[comp];
+ real res[] = {r->avg, r->err, r->chisq};
+ MLPutRealList(stdlink, res, Elements(res));
+ }
+
+ MLPutInteger(stdlink, region->df);
+ }
+ }
+#endif
+
+abort:
+ free(region);
+
+ while( (region = anchor) ) {
+ anchor = anchor->next;
+ free(region);
+ }
+
+ return fail;
+}
+
diff --git a/src/suave/Sample.c b/src/suave/Sample.c
new file mode 100644
index 0000000..33e0272
--- /dev/null
+++ b/src/suave/Sample.c
@@ -0,0 +1,171 @@
+/*
+ Sample.c
+ the sampling step of Suave
+ this file is part of Suave
+ last modified 12 Aug 11 th
+*/
+
+
+typedef struct {
+ real sum, sqsum;
+ real weight, weightsum, avg, avgsum;
+ real guess, chisum, chisqsum;
+} Cumulants;
+
+/*********************************************************************/
+
+static void Sample(This *t, cnumber nnew, void *voidregion,
+ real *lastw, real *lastx, real *lastf)
+{
+ TYPEDEFREGION;
+
+ Region *const region = (Region *)voidregion;
+ count comp, dim, df;
+ number n;
+ Cumulants cumul[NCOMP];
+ char **ss, *s;
+ ccount chars = 128*(region->div + 1);
+
+ creal jacobian = 1/ldexp((real)nnew, region->div);
+ real *w = lastw, *f = lastx;
+ bin_t *bin = (bin_t *)(lastf + nnew*t->ncomp);
+
+ for( n = nnew; n; --n ) {
+ real weight = jacobian;
+
+ t->rng.getrandom(t, f);
+
+ for( dim = 0; dim < t->ndim; ++dim ) {
+ cBounds *b = &region->bounds[dim];
+ creal pos = *f*NBINS;
+ ccount ipos = (count)pos;
+ creal prev = (ipos == 0) ? 0 : b->grid[ipos - 1];
+ creal diff = b->grid[ipos] - prev;
+ *f++ = b->lower + (prev + (pos - ipos)*diff)*(b->upper - b->lower);
+ *bin++ = ipos;
+ weight *= diff*NBINS;
+ }
+
+ *w++ = weight;
+ }
+
+ DoSample(t, nnew, lastx, lastf, lastw, region->div + 1);
+
+ w[-1] = -w[-1];
+ lastw = w;
+ w = region->w;
+ region->n = lastw - w;
+
+ if( VERBOSE > 2 ) {
+ char *p0;
+ MemAlloc(ss, t->ndim*64 + t->ncomp*(sizeof(char *) + chars));
+ s = (char *)(ss + t->ncomp);
+ p0 = s + t->ndim*64;
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ ss[comp] = p0;
+ p0 += chars;
+ }
+ }
+
+ Zap(cumul);
+ df = n = 0;
+
+ while( w < lastw ) {
+ cbool final = (*w < 0);
+ creal weight = fabs(*w++);
+ ++n;
+
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ Cumulants *c = &cumul[comp];
+
+ creal wfun = weight*(*f++);
+ c->sum += wfun;
+ c->sqsum += Sq(wfun);
+
+ if( final ) {
+ if( n > 1 ) {
+ real w = Weight(c->sum, c->sqsum, n);
+ c->weightsum += c->weight = w;
+ c->avgsum += c->avg = w*c->sum;
+
+ if( VERBOSE > 2 ) {
+ creal sig = sqrt(1/w);
+ ss[comp] += (df == 0) ?
+ sprintf(ss[comp], "\n[" COUNT "] "
+ REAL " +- " REAL " (" NUMBER ")", comp + 1,
+ c->sum, sig, n) :
+ sprintf(ss[comp], "\n "
+ REAL " +- " REAL " (" NUMBER ")",
+ c->sum, sig, n);
+ }
+
+ if( df == 0 ) c->guess = c->sum;
+ else {
+ c->chisum += w *= c->sum - c->guess;
+ c->chisqsum += w*c->sum;
+ }
+ }
+
+ c->sum = c->sqsum = 0;
+ }
+ }
+
+ if( final ) ++df, n = 0;
+ }
+
+ region->df = --df;
+
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ Result *r = &region->result[comp];
+ Cumulants *c = &cumul[comp];
+ creal sigsq = 1/c->weightsum;
+ creal avg = sigsq*c->avgsum;
+
+ if( LAST ) {
+ r->sigsq = 1/c->weight;
+ r->avg = r->sigsq*c->avg;
+ }
+ else {
+ r->sigsq = sigsq;
+ r->avg = avg;
+ }
+ r->err = sqrt(r->sigsq);
+
+ r->chisq = (sigsq < .9*NOTZERO) ? 0 : c->chisqsum - avg*c->chisum;
+ /* This catches the special case where the integrand is constant
+ over the entire region. Unless that constant is zero, only the
+ first set of samples will have zero variance, and hence weight
+ (n - 1) 1e30 (see above). All other sets have been sampled
+ from a non-constant weight function and therefore inevitably
+ show some variance. This is an artificial effect, brought about
+ by the fact that the constancy of the integrand in the region is
+ seen only in this subdivision, and can degrade the chi-square
+ score quite a bit. If the constancy was determined from more
+ than two samples (hence .9*NOTZERO), the chi-squares from the
+ other sets are removed here. */
+ }
+
+ if( VERBOSE > 2 ) {
+ char *p = s;
+ char *p0 = p + t->ndim*64;
+
+ for( dim = 0; dim < t->ndim; ++dim ) {
+ cBounds *b = &region->bounds[dim];
+ p += sprintf(p,
+ (dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" :
+ "\n (" REALF ") - (" REALF ")",
+ b->lower, b->upper);
+ }
+
+ for( comp = 0; comp < t->ncomp; ++comp ) {
+ cResult *r = &region->result[comp];
+ p += sprintf(p, "%s \tchisq " REAL " (" COUNT " df)",
+ p0, r->chisq, df);
+ p0 += chars;
+ }
+
+ Print(s);
+ free(ss);
+ }
+}
+
diff --git a/src/suave/Suave.c b/src/suave/Suave.c
new file mode 100644
index 0000000..bad5407
--- /dev/null
+++ b/src/suave/Suave.c
@@ -0,0 +1,91 @@
+/*
+ Suave.c
+ Subregion-adaptive Vegas Monte-Carlo integration
+ by Thomas Hahn
+ last modified 27 Sep 11 th
+*/
+
+
+#include "decl.h"
+
+#define Print(s) puts(s); fflush(stdout)
+
+/*********************************************************************/
+
+#define SUAVE
+#include "DoSample.c"
+
+/*********************************************************************/
+
+#include "common.c"
+
+Extern void EXPORT(Suave)(ccount ndim, ccount ncomp,
+ Integrand integrand, void *userdata,
+ creal epsrel, creal epsabs,
+ cint flags, cint seed,
+ cnumber mineval, cnumber maxeval,
+ cnumber nnew, creal flatness,
+ count *pnregions, number *pneval, int *pfail,
+ real *integral, real *error, real *prob)
+{
+ This t;
+ t.ndim = ndim;
+ t.ncomp = ncomp;
+ t.integrand = integrand;
+ t.userdata = userdata;
+ t.epsrel = epsrel;
+ t.epsabs = epsabs;
+ t.flags = flags;
+ t.seed = seed;
+ t.mineval = mineval;
+ t.maxeval = maxeval;
+ t.nnew = nnew;
+ t.flatness = flatness;
+ t.nregions = 0;
+ t.neval = 0;
+
+ ForkCores(&t);
+
+ *pfail = Integrate(&t, integral, error, prob);
+ *pnregions = t.nregions;
+ *pneval = t.neval;
+
+ WaitCores(&t);
+}
+
+/*********************************************************************/
+
+Extern void EXPORT(suave)(ccount *pndim, ccount *pncomp,
+ Integrand integrand, void *userdata,
+ creal *pepsrel, creal *pepsabs,
+ cint *pflags, cint *pseed,
+ cnumber *pmineval, cnumber *pmaxeval,
+ cnumber *pnnew, creal *pflatness,
+ count *pnregions, number *pneval, int *pfail,
+ real *integral, real *error, real *prob)
+{
+ This t;
+ t.ndim = *pndim;
+ t.ncomp = *pncomp;
+ t.integrand = integrand;
+ t.userdata = userdata;
+ t.epsrel = *pepsrel;
+ t.epsabs = *pepsabs;
+ t.flags = *pflags;
+ t.seed = *pseed;
+ t.mineval = *pmineval;
+ t.maxeval = *pmaxeval;
+ t.nnew = *pnnew;
+ t.flatness = *pflatness;
+ t.nregions = 0;
+ t.neval = 0;
+
+ ForkCores(&t);
+
+ *pfail = Integrate(&t, integral, error, prob);
+ *pnregions = t.nregions;
+ *pneval = t.neval;
+
+ WaitCores(&t);
+}
+
diff --git a/src/suave/Suave.tm b/src/suave/Suave.tm
new file mode 100644
index 0000000..4733d75
--- /dev/null
+++ b/src/suave/Suave.tm
@@ -0,0 +1,290 @@
+:Evaluate: BeginPackage["Cuba`"]
+
+:Evaluate: Suave::usage =
+ "Suave[f, {x, xmin, xmax}..] computes a numerical approximation to the integral of the real scalar or vector function f.
+ The output is a list with entries of the form {integral, error, chi-square probability} for each component of the integrand."
+
+:Evaluate: NNew::usage = "NNew is an option of Suave.
+ It specifies the number of new integrand evaluations in each subdivision."
+
+:Evaluate: Flatness::usage = "Flatness is an option of Suave.
+ It determines how prominently individual samples with a large fluctuation figure in the total fluctuation, which in turn determines how a region is split up.
+ Explicitly, if F[i] is the individual fluctuation of sample i, the total fluctuation is computed as Sum[(1 + F[i])^p, {i, nsamples}]^(2/3/p), i.e. as the p-norm of the fluctuation vector to the power 2/3, where p is the number given by Flatness.
+ Thus with increasing p, the fluctuation becomes more and more dominated by outliers, i.e. points with a large fluctuation.
+ As suggested by the name Flatness, p should be chosen large for `flat' integrands and small for `volatile' integrands with high peaks.
+ Note that since p appears in the exponent, one should not use too large values (say, no more than a few hundred) lest terms be truncated internally to prevent overflow."
+
+:Evaluate: MinPoints::usage = "MinPoints is an option of Suave.
+ It specifies the minimum number of points to sample."
+
+:Evaluate: Final::usage = "Final is an option of Suave.
+ It can take the values Last or All which determine whether only the last (largest) or all sets of samples collected on a subregion over the iterations contribute to the final result."
+
+:Evaluate: PseudoRandom::usage = "PseudoRandom is an option of Suave.
+ It can take the following values:
+ False for Sobol quasi-random numbers (default),
+ True or 0 for Mersenne Twister pseudo-random numbers,
+ any other integer value n for Ranlux pseudo-random numbers of luxury level n."
+
+:Evaluate: PseudoRandomSeed::usage = "PseudoRandomSeed is an option of Suave.
+ It specifies the seed for the pseudo-random number generator."
+
+:Evaluate: SharpEdges::usage = "SharpEdges is an option of Suave.
+ It turns off smoothing of the importance function for integrands with sharp edges."
+
+:Evaluate: Regions::usage = "Regions is an option of Suave.
+ It specifies whether the regions into which the integration region has been cut are returned together with the integration results."
+
+:Evaluate: Region::usage = "Region[ll, ur, res, df] describes a subregion:
+ ll and ur are multidimensional equivalents of the region's lower left and upper right corner.
+ res gives the integration results for the region in a list with entries of the form {integral, error, chi-square} for each component of the integrand.
+ df is the number of degrees of freedom corresponding to the chi-square values in res."
+
+:Evaluate: $Weight::usage = "$Weight is a global variable set by Suave during the evaluation of the integrand to the weight of the point being sampled."
+
+:Evaluate: $Iteration::usage = "$Iteration is a global variable set by Suave during the evaluation of the integrand to the present iteration number."
+
+:Evaluate: MapSample::usage = "MapSample is a function used to map the integrand over the points to be sampled."
+
+
+:Evaluate: Begin["`Suave`"]
+
+:Begin:
+:Function: Suave
+:Pattern: MLSuave[ndim_, ncomp_,
+ epsrel_, epsabs_, flags_, seed_,
+ mineval_, maxeval_,
+ nnew_, flatness_]
+:Arguments: {ndim, ncomp,
+ epsrel, epsabs, flags, seed,
+ mineval, maxeval,
+ nnew, flatness}
+:ArgumentTypes: {Integer, Integer,
+ Real, Real, Integer, Integer,
+ Integer, Integer,
+ Integer, Real}
+:ReturnType: Manual
+:End:
+
+:Evaluate: Attributes[Suave] = {HoldFirst}
+
+:Evaluate: Options[Suave] = {PrecisionGoal -> 3, AccuracyGoal -> 12,
+ MinPoints -> 0, MaxPoints -> 50000, NNew -> 1000, Flatness -> 50,
+ Verbose -> 1, Final -> Last,
+ PseudoRandom -> False, PseudoRandomSeed -> 5489,
+ SharpEdges -> False, Regions -> False, Compiled -> True}
+
+:Evaluate: Suave[f_, v:{_, _, _}.., opt___Rule] :=
+ Block[ {ff = HoldForm[f], ndim = Length[{v}], ncomp,
+ tags, vars, lower, range, jac, tmp, defs, intT,
+ rel, abs, mineval, maxeval, nnew, flatness,
+ verbose, final, level, seed, edges, regions, compiled,
+ $Weight, $Iteration},
+ Message[Suave::optx, #, Suave]&/@
+ Complement[First/@ {opt}, tags = First/@ Options[Suave]];
+ {rel, abs, mineval, maxeval, nnew, flatness,
+ verbose, final, level, seed, edges, regions, compiled} =
+ tags /. {opt} /. Options[Suave];
+ {vars, lower, range} = Transpose[{v}];
+ jac = Simplify[Times@@ (range -= lower)];
+ tmp = Array[tmpvar, ndim];
+ defs = Simplify[lower + range tmp];
+ Block[{Set}, define[compiled, tmp, Thread[vars = defs], jac]];
+ intT = integrandT[f];
+ Block[#,
+ ncomp = Length[intT@@ RandomReal[1, ndim]];
+ MLSuave[ndim, ncomp, 10.^-rel, 10.^-abs,
+ Min[Max[verbose, 0], 3] +
+ If[final === Last, 4, 0] +
+ If[TrueQ[edges], 8, 0] +
+ If[TrueQ[regions], 128, 0] +
+ If[IntegerQ[level], 256 level, 0],
+ If[level =!= False && IntegerQ[seed], seed, 0],
+ mineval, maxeval,
+ nnew, flatness]
+ ]& @ vars
+ ]
+
+:Evaluate: tmpvar[n_] := ToExpression["Cuba`Suave`t" <> ToString[n]]
+
+:Evaluate: Attributes[foo] = {HoldAll}
+
+:Evaluate: define[True, tmp_, defs_, jac_] := (
+ TtoX := TtoX = Compile[tmp, defs];
+ integrandT[f_] := Compile[tmp, eval[defs, Chop[f jac]//N],
+ {{_eval, _Real, 1}}] )
+
+:Evaluate: define[_, tmp_, defs_, jac_] := (
+ TtoX := TtoX = Function[tmp, defs];
+ integrandT[f_] := Function[tmp, eval[defs, Chop[f jac]//N]] )
+
+:Evaluate: eval[_, f_Real] = {f}
+
+:Evaluate: eval[_, f:{__Real}] = f
+
+:Evaluate: eval[x_, _] := (Message[Suave::badsample, ff, x]; {})
+
+:Evaluate: sample[x_, w_, iter_] := (
+ $Iteration = iter;
+ Check[Flatten @ MapSample[
+ ($Weight = #[[1]]; intT@@ #[[2]])&,
+ Transpose[{w, Partition[x, ndim]}] ], {}] )
+
+:Evaluate: MapSample = Map
+
+:Evaluate: region[ll_, ur_, r___] := Region[TtoX@@ ll, TtoX@@ ur, r]
+
+:Evaluate: Suave::badsample = "`` is not a real-valued function at ``."
+
+:Evaluate: Suave::baddim = "Cannot integrate in `` dimensions."
+
+:Evaluate: Suave::badcomp = "Cannot integrate `` components."
+
+:Evaluate: Suave::accuracy =
+ "Desired accuracy was not reached within `` function evaluations on `` subregions."
+
+:Evaluate: Suave::success = "Needed `` function evaluations on `` subregions."
+
+:Evaluate: End[]
+
+:Evaluate: EndPackage[]
+
+
+/*
+ Suave.tm
+ Subregion-adaptive Vegas Monte Carlo integration
+ by Thomas Hahn
+ last modified 12 Aug 11 th
+*/
+
+
+#include "mathlink.h"
+#include "decl.h"
+
+/*********************************************************************/
+
+static void Status(MLCONST char *msg, cint n1, cint n2)
+{
+ MLPutFunction(stdlink, "CompoundExpression", 2);
+ MLPutFunction(stdlink, "Message", 3);
+ MLPutFunction(stdlink, "MessageName", 2);
+ MLPutSymbol(stdlink, "Suave");
+ MLPutString(stdlink, msg);
+ MLPutInteger(stdlink, n1);
+ MLPutInteger(stdlink, n2);
+}
+
+/*********************************************************************/
+
+static void Print(MLCONST char *s)
+{
+ MLPutFunction(stdlink, "EvaluatePacket", 1);
+ MLPutFunction(stdlink, "Print", 1);
+ MLPutString(stdlink, s);
+ MLEndPacket(stdlink);
+
+ MLNextPacket(stdlink);
+ MLNewPacket(stdlink);
+}
+
+/*********************************************************************/
+
+static void DoSample(This *t, cnumber n, real *x, real *f,
+ real *w, cint iter)
+{
+ real *mma_f;
+ long mma_n;
+
+ if( MLAbort ) longjmp(t->abort, -99);
+
+ MLPutFunction(stdlink, "EvaluatePacket", 1);
+ MLPutFunction(stdlink, "Cuba`Suave`sample", 3);
+ MLPutRealList(stdlink, x, n*t->ndim);
+ MLPutRealList(stdlink, w, n);
+ MLPutInteger(stdlink, iter);
+ MLEndPacket(stdlink);
+
+ MLNextPacket(stdlink);
+ if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) {
+ MLClearError(stdlink);
+ MLNewPacket(stdlink);
+ longjmp(t->abort, -99);
+ }
+
+ if( mma_n != n*t->ncomp ) {
+ MLDisownRealList(stdlink, mma_f, mma_n);
+ longjmp(t->abort, -3);
+ }
+
+ Copy(f, mma_f, n*t->ncomp);
+ MLDisownRealList(stdlink, mma_f, mma_n);
+
+ t->neval += n;
+}
+
+/*********************************************************************/
+
+#include "common.c"
+
+static inline void DoIntegrate(This *t)
+{
+ real integral[NCOMP], error[NCOMP], prob[NCOMP];
+ cint fail = Integrate(t, integral, error, prob);
+
+ if( fail < 0 ) {
+ switch( fail ) {
+ case -99:
+ MLPutFunction(stdlink, "Abort", 0);
+ return;
+ case -1:
+ Status("baddim", t->ndim, 0);
+ break;
+ case -2:
+ Status("badcomp", t->ncomp, 0);
+ break;
+ }
+ MLPutSymbol(stdlink, "$Failed");
+ }
+ else {
+ Status(fail ? "accuracy" : "success", t->neval, t->nregions);
+ MLPutFunction(stdlink, "Thread", 1);
+ MLPutFunction(stdlink, "List", 3);
+ MLPutRealList(stdlink, integral, t->ncomp);
+ MLPutRealList(stdlink, error, t->ncomp);
+ MLPutRealList(stdlink, prob, t->ncomp);
+ }
+}
+
+/*********************************************************************/
+
+void Suave(cint ndim, cint ncomp,
+ creal epsrel, creal epsabs,
+ cint flags, cint seed,
+ cnumber mineval, cnumber maxeval,
+ cnumber nnew, creal flatness)
+{
+ This t;
+ t.ndim = ndim;
+ t.ncomp = ncomp;
+ t.epsrel = epsrel;
+ t.epsabs = epsabs;
+ t.flags = flags;
+ t.seed = seed;
+ t.mineval = mineval;
+ t.maxeval = maxeval;
+ t.nnew = nnew;
+ t.flatness = flatness;
+ t.nregions = 0;
+ t.neval = 0;
+
+ DoIntegrate(&t);
+ MLEndPacket(stdlink);
+}
+
+/*********************************************************************/
+
+int main(int argc, char **argv)
+{
+ return MLMain(argc, argv);
+}
+
diff --git a/src/suave/common.c b/src/suave/common.c
new file mode 100644
index 0000000..1325911
--- /dev/null
+++ b/src/suave/common.c
@@ -0,0 +1,33 @@
+/*
+ common.c
+ includes most of the modules
+ this file is part of Suave
+ last modified 2 Jun 10 th
+*/
+
+
+#define RegionAlloc(p, n, nnew) MemAlloc(p, \
+ sizeof(Region) + \
+ (n)*(t->ndim + t->ncomp + 1)*sizeof(real) + \
+ (nnew)*t->ndim*sizeof(bin_t))
+
+static inline bool BadDimension(cThis *t)
+{
+ if( t->ndim > NDIM ) return true;
+ return t->ndim < SOBOL_MINDIM ||
+ (t->seed == 0 && t->ndim > SOBOL_MAXDIM);
+}
+
+static inline bool BadComponent(cThis *t)
+{
+ if( t->ncomp > NCOMP ) return true;
+ return t->ncomp < 1;
+}
+
+#include "Random.c"
+#include "ChiSquare.c"
+#include "Grid.c"
+#include "Sample.c"
+#include "Fluct.c"
+#include "Integrate.c"
+
diff --git a/src/suave/decl.h b/src/suave/decl.h
new file mode 100644
index 0000000..12ec957
--- /dev/null
+++ b/src/suave/decl.h
@@ -0,0 +1,70 @@
+/*
+ decl.h
+ Type declarations
+ this file is part of Suave
+ last modified 4 Oct 11 th
+*/
+
+
+#include "stddecl.h"
+
+#define MINSAMPLES 10
+
+#define NBINS 64
+
+typedef unsigned char bin_t;
+/* Note: bin_t must be wide enough to hold the numbers 0..NBINS */
+
+typedef const bin_t cbin_t;
+
+typedef real Grid[NBINS];
+
+typedef const Grid cGrid;
+
+typedef struct {
+ real avg, err, sigsq, chisq;
+} Result;
+
+typedef const Result cResult;
+
+typedef struct {
+ real lower, upper, mid;
+ Grid grid;
+} Bounds;
+
+typedef const Bounds cBounds;
+
+typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
+ void *, creal *, cint *);
+
+typedef struct _this {
+ count ndim, ncomp;
+#ifndef MLVERSION
+ Integrand integrand;
+ void *userdata;
+ int ncores, *child;
+#endif
+ real epsrel, epsabs;
+ int flags, seed;
+ number mineval, maxeval;
+ number nnew;
+ real flatness;
+ count nregions;
+ number neval;
+ RNGState rng;
+ jmp_buf abort;
+} This;
+
+typedef const This cThis;
+
+#define TYPEDEFREGION \
+ typedef struct region { \
+ struct region *next; \
+ count div, df; \
+ number n; \
+ Result result[NCOMP]; \
+ Bounds bounds[NDIM]; \
+ real fluct[NCOMP][NDIM][2]; \
+ real w[]; \
+ } Region
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-21 15:10:18 +0000