Imported Upstream version 3.0+20111124upstream/3.0+20111124upstream

6 files changed, 822 insertions, 0 deletions

diff --git a/src/vegas/Grid.c b/src/vegas/Grid.c
new file mode 100644
index 0000000..e99138f
--- /dev/null
+++ b/src/vegas/Grid.c
@@ -0,0 +1,99 @@
+/*
+	Grid.c
+		utility functions for the Vegas grid
+		this file is part of Vegas
+		last modified 28 May 10 th
+*/
+
+
+static inline void GetGrid(cThis *t, Grid *grid)
+{
+  count bin, dim;
+  unsigned const int slot = abs(t->gridno) - 1;
+
+  if( t->gridno < 0 && slot < MAXGRIDS ) griddim_[slot] = 0;
+
+  if( slot < MAXGRIDS && gridptr_[slot] ) {
+    if( griddim_[slot] == t->ndim ) {
+      VecCopy(grid, gridptr_[slot]);
+      return;
+    }
+    free(gridptr_[slot]);
+    gridptr_[slot] = NULL;
+  }
+
+  for( bin = 0; bin < NBINS; ++bin )
+    grid[0][bin] = (bin + 1)/(real)NBINS;
+  for( dim = 1; dim < t->ndim; ++dim )
+    Copy(&grid[dim], &grid[0], 1);
+}
+
+/*********************************************************************/
+
+static inline void PutGrid(cThis *t, Grid *grid)
+{
+  unsigned const int slot = abs(t->gridno) - 1;
+
+  if( slot < MAXGRIDS ) {
+    if( gridptr_[slot] == NULL ) Alloc(gridptr_[slot], t->ndim);
+    griddim_[slot] = t->ndim;
+    VecCopy(gridptr_[slot], grid);
+  }
+}
+
+/*********************************************************************/
+
+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,
+        cur - 2*delta/(imp[bin] + imp[IDim(bin - 1)])));
+  }
+  Copy(grid, newgrid, NBINS - 1);
+  grid[NBINS - 1] = 1;
+}
+
diff --git a/src/vegas/Integrate.c b/src/vegas/Integrate.c
new file mode 100644
index 0000000..efa582e
--- /dev/null
+++ b/src/vegas/Integrate.c
@@ -0,0 +1,239 @@
+/*
+	Integrate.c
+		integrate over the unit hypercube
+		this file is part of Vegas
+		last modified 19 Aug 11 th
+*/
+
+
+static int Integrate(This *t, real *integral, real *error, real *prob)
+{
+  real *sample;
+  count dim, comp;
+  int fail;
+  struct {
+    count niter;
+    number nsamples, neval;
+    Cumulants cumul[NCOMP];
+    Grid grid[NDIM];
+  } state;
+  int statemsg = VERBOSE;
+  struct stat st;
+
+  if( VERBOSE > 1 ) {
+    char s[512];
+    sprintf(s, "Vegas 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"
+      "  nstart " NUMBER "\n  nincrease " NUMBER "\n"
+      "  nbatch " NUMBER "\n  gridno %d\n"
+      "  statefile \"%s\"",
+      t->ndim, t->ncomp,
+      t->epsrel, t->epsabs,
+      t->flags, t->seed,
+      t->mineval, t->maxeval,
+      t->nstart, t->nincrease, t->nbatch,
+      t->gridno, t->statefile);
+    Print(s);
+  }
+
+  if( BadComponent(t) ) return -2;
+  if( BadDimension(t) ) return -1;
+
+  SamplesAlloc(sample);
+
+  if( (fail = setjmp(t->abort)) ) goto abort;
+
+  IniRandom(t);
+
+  if( t->statefile && *t->statefile == 0 ) t->statefile = NULL;
+
+  if( t->statefile &&
+      stat(t->statefile, &st) == 0 &&
+      st.st_size == sizeof state && (st.st_mode & 0400) ) {
+    cint h = open(t->statefile, O_RDONLY);
+    read(h, &state, sizeof state);
+    close(h);
+    t->rng.skiprandom(t, t->neval = state.neval);
+
+    if( VERBOSE ) {
+      char s[256];
+      sprintf(s, "\nRestoring state from %s.", t->statefile);
+      Print(s);
+    }
+  }
+  else {
+    state.niter = 0;
+    state.nsamples = t->nstart;
+    Zap(state.cumul);
+    GetGrid(t, state.grid);
+  }
+
+  /* main iteration loop */
+
+  for( ; ; ) {
+    number nsamples = state.nsamples;
+    creal jacobian = 1./nsamples;
+    Grid margsum[NCOMP][NDIM];
+
+    Zap(margsum);
+
+    for( ; nsamples > 0; nsamples -= t->nbatch ) {
+      cnumber n = IMin(t->nbatch, nsamples);
+      real *w = sample;
+      real *x = w + n;
+      real *f = x + n*t->ndim;
+      real *lastf = f + n*t->ncomp;
+      bin_t *bin = (bin_t *)lastf;
+
+      while( x < f ) {
+        real weight = jacobian;
+
+        t->rng.getrandom(t, x);
+
+        for( dim = 0; dim < t->ndim; ++dim ) {
+          creal pos = *x*NBINS;
+          ccount ipos = (count)pos;
+          creal prev = (ipos == 0) ? 0 : state.grid[dim][ipos - 1];
+          creal diff = state.grid[dim][ipos] - prev; 
+          *x++ = prev + (pos - ipos)*diff;
+          *bin++ = ipos;
+          weight *= diff*NBINS;
+        }
+
+        *w++ = weight;
+      }
+
+      DoSample(t, n, w, f, sample, state.niter + 1);
+
+      w = sample;
+      bin = (bin_t *)lastf;
+
+      while( f < lastf ) {
+        creal weight = *w++;
+
+        for( comp = 0; comp < t->ncomp; ++comp ) {
+          real wfun = weight*(*f++);
+          if( wfun ) {
+            Cumulants *c = &state.cumul[comp];
+            Grid *m = margsum[comp];
+
+            c->sum += wfun;
+            c->sqsum += wfun *= wfun;
+            for( dim = 0; dim < t->ndim; ++dim )
+              m[dim][bin[dim]] += wfun;
+          }
+        }
+
+        bin += t->ndim;
+      }
+    }
+
+    fail = 0;
+
+    /* compute the integral and error values */
+
+    for( comp = 0; comp < t->ncomp; ++comp ) {
+      Cumulants *c = &state.cumul[comp];
+      real avg, sigsq;
+      real w = Weight(c->sum, c->sqsum, state.nsamples);
+
+      sigsq = 1/(c->weightsum += w);
+      avg = sigsq*(c->avgsum += w*c->sum);
+
+      c->avg = LAST ? (sigsq = 1/w, c->sum) : avg;
+      c->err = sqrt(sigsq);
+      fail |= (c->err > MaxErr(c->avg));
+
+      if( state.niter == 0 ) c->guess = c->sum;
+      else {
+        c->chisum += w *= c->sum - c->guess;
+        c->chisqsum += w*c->sum;
+      }
+      c->chisq = c->chisqsum - avg*c->chisum;
+
+      c->sum = c->sqsum = 0;
+    }
+
+    if( VERBOSE ) {
+      char s[128 + 128*NCOMP], *p = s;
+
+      p += sprintf(p, "\n"
+        "Iteration " COUNT ":  " NUMBER " integrand evaluations so far",
+        state.niter + 1, t->neval);
+
+      for( comp = 0; comp < t->ncomp; ++comp ) {
+        cCumulants *c = &state.cumul[comp];
+        p += sprintf(p, "\n[" COUNT "] "
+          REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
+          comp + 1, c->avg, c->err, c->chisq, state.niter);
+      }
+
+      Print(s);
+    }
+
+    if( fail == 0 && t->neval >= t->mineval ) {
+      if( t->statefile ) unlink(t->statefile);
+      break;
+    }
+
+    if( t->neval >= t->maxeval && t->statefile == NULL ) break;
+
+    if( t->ncomp == 1 )
+      for( dim = 0; dim < t->ndim; ++dim )
+        RefineGrid(t, state.grid[dim], margsum[0][dim]);
+    else {
+      for( dim = 0; dim < t->ndim; ++dim ) {
+        Grid wmargsum;
+        Zap(wmargsum);
+        for( comp = 0; comp < t->ncomp; ++comp ) {
+          real w = state.cumul[comp].avg;
+          if( w != 0 ) {
+            creal *m = margsum[comp][dim];
+            count bin;
+            w = 1/Sq(w);
+            for( bin = 0; bin < NBINS; ++bin )
+              wmargsum[bin] += w*m[bin];
+          }
+        }
+        RefineGrid(t, state.grid[dim], wmargsum);
+      }
+    }
+
+    ++state.niter;
+    state.nsamples += t->nincrease;
+
+    if( t->statefile ) {
+      cint h = creat(t->statefile, 0666);
+      if( h != -1 ) {
+        state.neval = t->neval;
+        write(h, &state, sizeof state);
+        close(h);
+
+        if( statemsg ) {
+          char s[256];
+          sprintf(s, "\nSaving state to %s.", t->statefile);
+          Print(s);
+          statemsg = false;
+        }
+      }
+      if( t->neval >= t->maxeval ) break;
+    }
+  }
+
+  for( comp = 0; comp < t->ncomp; ++comp ) {
+    cCumulants *c = &state.cumul[comp];
+    integral[comp] = c->avg;
+    error[comp] = c->err;
+    prob[comp] = ChiSquare(c->chisq, state.niter);
+  }
+
+abort:
+  free(sample);
+  PutGrid(t, state.grid);
+
+  return fail;
+}
+
diff --git a/src/vegas/Vegas.c b/src/vegas/Vegas.c
new file mode 100644
index 0000000..08c331f
--- /dev/null
+++ b/src/vegas/Vegas.c
@@ -0,0 +1,105 @@
+/*
+	Vegas.c
+		Vegas Monte-Carlo integration
+		by Thomas Hahn
+		last modified 27 Sep 11 th
+*/
+
+
+#include "decl.h"
+
+#define Print(s) puts(s); fflush(stdout)
+
+/*********************************************************************/
+
+#define VEGAS
+#include "DoSample.c"
+
+/*********************************************************************/
+
+#include "common.c"
+
+Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp,
+  Integrand integrand, void *userdata,
+  creal epsrel, creal epsabs, cint flags, cint seed,
+  cnumber mineval, cnumber maxeval,
+  cnumber nstart, cnumber nincrease, cnumber nbatch,
+  cint gridno, cchar *statefile,
+  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.nstart = nstart;
+  t.nincrease = nincrease;
+  t.nbatch = nbatch;
+  t.gridno = gridno;
+  t.statefile = statefile;
+  t.neval = 0;
+
+  ForkCores(&t);
+
+  *pfail = Integrate(&t, integral, error, prob);
+  *pneval = t.neval;
+
+  WaitCores(&t);
+}
+
+/*********************************************************************/
+
+Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp,
+  Integrand integrand, void *userdata,
+  creal *pepsrel, creal *pepsabs, cint *pflags, cint *pseed,
+  cnumber *pmineval, cnumber *pmaxeval,
+  cnumber *pnstart, cnumber *pnincrease, 
+  cnumber *pnbatch, cint *pgridno, cchar *statefile,
+  number *pneval, int *pfail,
+  real *integral, real *error, real *prob, int statefilelen)
+{
+  char *s = NULL;
+  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.nstart = *pnstart;
+  t.nincrease = *pnincrease;
+  t.nbatch = *pnbatch;
+  t.gridno = *pgridno;
+  t.neval = 0;
+
+  if( statefile ) {
+	/* strip trailing spaces */
+    while( statefilelen > 0 && statefile[statefilelen - 1] == ' ' )
+      --statefilelen;
+    if( statefilelen > 0 && (s = malloc(statefilelen + 1)) ) {
+      memcpy(s, statefile, statefilelen);
+      s[statefilelen] = 0;
+    }
+  }
+  t.statefile = s;
+
+  ForkCores(&t);
+
+  *pfail = Integrate(&t, integral, error, prob);
+  *pneval = t.neval;
+
+  free(s);
+  WaitCores(&t);
+}
+
diff --git a/src/vegas/Vegas.tm b/src/vegas/Vegas.tm
new file mode 100644
index 0000000..039beed
--- /dev/null
+++ b/src/vegas/Vegas.tm
@@ -0,0 +1,292 @@
+:Evaluate: BeginPackage["Cuba`"]
+
+:Evaluate: Vegas::usage = "Vegas[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: NStart::usage = "NStart is an option of Vegas.
+	It specifies the number of integrand evaluations per iteration to start with."
+
+:Evaluate: NIncrease::usage = "NIncrease is an option of Vegas.
+	It specifies the increase in the number of integrand evaluations per iteration."
+
+:Evaluate: NBatch::usage = "NBatch is an option of Vegas.
+	It specifies how many points are sent in one MathLink packet to be sampled by Mathematica."
+
+:Evaluate: GridNo::usage = "GridNo is an option of Vegas.
+	Vegas maintains an internal table in which it can memorize up to 10 grids, to be used on subsequent integrations.
+	A GridNo between 1 and 10 selects the slot in this internal table.
+	For other values the grid is initialized from scratch and discarded at the end of the integration."
+
+:Evaluate: StateFile::usage = "StateFile is an option of Vegas.
+	It specifies a file in which the internal state is stored after each iteration and from which it can be restored on a subsequent run.
+	The state file is removed once the prescribed accuracy has been reached."
+
+:Evaluate: MinPoints::usage = "MinPoints is an option of Vegas.
+	It specifies the minimum number of points to sample."
+
+:Evaluate: Final::usage = "Final is an option of Vegas.
+	It can take the values Last or All which determine whether only the last (largest) or all of the samples collected on a subregion over the iterations contribute to the final result."
+
+:Evaluate: PseudoRandom::usage = "PseudoRandom is an option of Vegas.
+	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 Vegas.
+	It specifies the seed for the pseudo-random number generator."
+
+:Evaluate: SharpEdges::usage = "SharpEdges is an option of Vegas.
+	It turns off smoothing of the importance function for integrands with sharp edges."
+
+:Evaluate: $Weight::usage = "$Weight is a global variable set by Vegas 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["`Vegas`"]
+
+:Begin:
+:Function: Vegas
+:Pattern: MLVegas[ndim_, ncomp_,
+  epsrel_, epsabs_, flags_, seed_,
+  mineval_, maxeval_,
+  nstart_, nincrease_, nbatch_,
+  gridno_, statefile_]
+:Arguments: {ndim, ncomp,
+  epsrel, epsabs, flags, seed,
+  mineval, maxeval,
+  nstart, nincrease, nbatch,
+  gridno, statefile}
+:ArgumentTypes: {Integer, Integer,
+  Real, Real, Integer, Integer,
+  Integer, Integer,
+  Integer, Integer, Integer,
+  Integer, String}
+:ReturnType: Manual
+:End:
+
+:Evaluate: Attributes[Vegas] = {HoldFirst}
+
+:Evaluate: Options[Vegas] = {PrecisionGoal -> 3, AccuracyGoal -> 12,
+	MinPoints -> 0, MaxPoints -> 50000,
+	NStart -> 1000, NIncrease -> 500,
+	NBatch -> 1000, GridNo -> 0, StateFile -> "",
+	Verbose -> 1, Final -> All,
+	PseudoRandom -> False, PseudoRandomSeed -> 5489,
+	SharpEdges -> False, Compiled -> True}
+
+:Evaluate: Vegas[f_, v:{_, _, _}.., opt___Rule] :=
+	Block[ {ff = HoldForm[f], ndim = Length[{v}], ncomp,
+	tags, vars, lower, range, jac, tmp, defs, intT,
+	rel, abs, mineval, maxeval, nstart, nincrease, nbatch,
+	gridno, verbose, final, level, seed, edges, compiled,
+	$Weight, $Iteration},
+	  Message[Vegas::optx, #, Vegas]&/@
+	    Complement[First/@ {opt}, tags = First/@ Options[Vegas]];
+	  {rel, abs, mineval, maxeval, nstart, nincrease, nbatch,
+	    gridno, state, verbose, final, level, seed, edges, compiled} =
+	    tags /. {opt} /. Options[Vegas];
+	  {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]];
+	    MLVegas[ndim, ncomp, 10.^-rel, 10.^-abs,
+	      Min[Max[verbose, 0], 3] +
+	        If[final === Last, 4, 0] +
+	        If[TrueQ[edges], 8, 0]
+	        If[IntegerQ[level], 256 level, 0],
+	      If[level =!= False && IntegerQ[seed], seed, 0],
+	      mineval, maxeval,
+	      nstart, nincrease, nbatch,
+	      gridno, state]
+	  ]& @ vars
+	]
+
+:Evaluate: tmpvar[n_] := ToExpression["Cuba`Vegas`t" <> ToString[n]]
+
+:Evaluate: Attributes[foo] = {HoldAll}
+
+:Evaluate: define[True, tmp_, defs_, jac_] :=
+	integrandT[f_] := Compile[tmp, eval[defs, Chop[f jac]//N],
+	  {{_eval, _Real, 1}}]
+
+:Evaluate: define[_, tmp_, defs_, jac_] :=
+	integrandT[f_] := Function[tmp, eval[defs, Chop[f jac]//N]]
+
+:Evaluate: eval[_, f_Real] = {f}
+
+:Evaluate: eval[_, f:{__Real}] = f
+
+:Evaluate: eval[x_, _] := (Message[Vegas::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: Vegas::badsample = "`` is not a real-valued function at ``."
+
+:Evaluate: Vegas::baddim = "Cannot integrate in `` dimensions."
+
+:Evaluate: Vegas::badcomp = "Cannot integrate `` components."
+
+:Evaluate: Vegas::accuracy =
+	"Desired accuracy was not reached within `` function evaluations."
+
+:Evaluate: Vegas::success = "Needed `` function evaluations."
+
+:Evaluate: End[]
+
+:Evaluate: EndPackage[]
+
+
+/*
+	Vegas.tm
+		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 n)
+{
+  MLPutFunction(stdlink, "CompoundExpression", 2);
+  MLPutFunction(stdlink, "Message", 2);
+  MLPutFunction(stdlink, "MessageName", 2);
+  MLPutSymbol(stdlink, "Vegas");
+  MLPutString(stdlink, msg);
+  MLPutInteger(stdlink, n);
+}
+
+/*********************************************************************/
+
+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`Vegas`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);
+      break;
+    case -2:
+      Status("badcomp", t->ncomp);
+      break;
+    }
+    MLPutSymbol(stdlink, "$Failed");
+  }
+  else {
+    Status(fail ? "accuracy" : "success", t->neval);
+    MLPutFunction(stdlink, "Thread", 1);
+    MLPutFunction(stdlink, "List", 3);
+    MLPutRealList(stdlink, integral, t->ncomp);
+    MLPutRealList(stdlink, error, t->ncomp);
+    MLPutRealList(stdlink, prob, t->ncomp);
+  }
+}
+
+/*********************************************************************/
+
+void Vegas(cint ndim, cint ncomp,
+  creal epsrel, creal epsabs,
+  cint flags, cint seed,
+  cnumber mineval, cnumber maxeval,
+  cnumber nstart, cnumber nincrease, cint nbatch,
+  cint gridno, cchar *statefile)
+{
+  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.nstart = nstart;
+  t.nincrease = nincrease;
+  t.nbatch = nbatch;
+  t.gridno = gridno;
+  t.statefile = statefile;
+  t.neval = 0;
+
+  DoIntegrate(&t);
+  MLEndPacket(stdlink);
+}
+
+/*********************************************************************/
+
+int main(int argc, char **argv)
+{
+  return MLMain(argc, argv);
+}
+
diff --git a/src/vegas/common.c b/src/vegas/common.c
new file mode 100644
index 0000000..fe6a0d7
--- /dev/null
+++ b/src/vegas/common.c
@@ -0,0 +1,31 @@
+/*
+	common.c
+		Code common to Vegas.c and Vegas.tm
+		this file is part of Vegas
+		last modified 6 Jun 10 th
+*/
+
+
+#include "Random.c"
+#include "ChiSquare.c"
+#include "Grid.c"
+
+#define SamplesAlloc(p) MemAlloc(p, \
+  t->nbatch*((t->ndim + t->ncomp + 1)*sizeof(real) + \
+             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 "Integrate.c"
+
diff --git a/src/vegas/decl.h b/src/vegas/decl.h
new file mode 100644
index 0000000..1162ff2
--- /dev/null
+++ b/src/vegas/decl.h
@@ -0,0 +1,56 @@
+/*
+	decl.h
+		Type declarations
+		this file is part of Vegas
+		last modified 4 Oct 11 th
+*/
+
+
+#include "stddecl.h"
+
+#define MAXGRIDS 10
+
+#define NBINS 128
+
+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 struct {
+  real sum, sqsum;
+  real weightsum, avgsum;
+  real chisum, chisqsum, guess;
+  real avg, err, chisq;
+} Cumulants;
+
+typedef const Cumulants cCumulants;
+
+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 nstart, nincrease, nbatch;
+  int gridno;
+  cchar *statefile;
+  number neval;
+  RNGState rng;
+  jmp_buf abort;
+} This;
+
+typedef const This cThis;
+
+static Grid *gridptr_[MAXGRIDS];
+static count griddim_[MAXGRIDS];
+