1 files changed, 137 insertions, 0 deletions
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;
+}
+
+