1 files changed, 160 insertions, 0 deletions
diff --git a/src/divonne/Explore.c b/src/divonne/Explore.c
new file mode 100644
index 0000000..d6fba19
--- /dev/null
+++ b/src/divonne/Explore.c
@@ -0,0 +1,160 @@
+/*
+	Explore.c
+		sample region, determine min and max, split if necessary
+		this file is part of Divonne
+		last modified 15 Nov 11 th
+*/
+
+
+typedef struct {
+  real fmin, fmax;
+  creal *xmin, *xmax;
+} Extrema;
+
+/*********************************************************************/
+
+static int Explore(This *t, ccount iregion)
+{
+  TYPEDEFREGION;
+  Region *region = RegionPtr(iregion);
+  cBounds *bounds = region->bounds;
+  Result *result = region->result;
+
+  count n, dim, comp, maxcomp;
+  Extrema extrema[NCOMP];
+  Result *r;
+  creal *x;
+  real *f;
+  real halfvol, maxerr;
+  cSamples *samples = &t->samples[region->isamples];
+
+  /* needed as of gcc 3.3 to make gcc correctly address region #@$&! */
+  sizeof(*region);
+
+  for( comp = 0; comp < t->ncomp; ++comp ) {
+    Extrema *e = &extrema[comp];
+    e->fmin = INFTY;
+    e->fmax = -INFTY;
+    e->xmin = e->xmax = NULL;
+  }
+
+  if( region->isamples == 0 ) {		/* others already sampled */
+    real vol = 1;
+    for( dim = 0; dim < t->ndim; ++dim ) {
+      cBounds *b = &bounds[dim];
+      vol *= b->upper - b->lower;
+    }
+    region->vol = vol;
+
+    for( comp = 0; comp < t->ncomp; ++comp ) {
+      Result *r = &result[comp];
+      r->fmin = INFTY;
+      r->fmax = -INFTY;
+    }
+
+    x = t->xgiven;
+    f = t->fgiven;
+    n = t->ngiven;
+    if( t->nextra ) n += SampleExtra(t, bounds);
+
+    for( ; n; --n ) {
+      for( dim = 0; dim < t->ndim; ++dim ) {
+        cBounds *b = &bounds[dim];
+        if( x[dim] < b->lower || x[dim] > b->upper ) goto skip;
+      }
+      for( comp = 0; comp < t->ncomp; ++comp ) {
+        Extrema *e = &extrema[comp];
+        creal y = f[comp];
+        if( y < e->fmin ) e->fmin = y, e->xmin = x;
+        if( y > e->fmax ) e->fmax = y, e->xmax = x;
+      }
+skip:
+      x += t->ldxgiven;
+      f += t->ncomp;
+    }
+
+    samples->sampler(t, iregion);
+  }
+
+  x = samples->x;
+  f = samples->f;
+  for( n = samples->n; n; --n ) {
+    for( comp = 0; comp < t->ncomp; ++comp ) {
+      Extrema *e = &extrema[comp];
+      creal y = *f++;
+      if( y < e->fmin ) e->fmin = y, e->xmin = x;
+      if( y > e->fmax ) e->fmax = y, e->xmax = x;
+    }
+    x += t->ndim;
+  }
+  t->neval_opt -= t->neval;
+
+  halfvol = .5*region->vol;
+  maxerr = -INFTY;
+  maxcomp = -1;
+
+  for( comp = 0; comp < t->ncomp; ++comp ) {
+    Extrema *e = &extrema[comp];
+    Result *r = &result[comp];
+    real xtmp[NDIM], ftmp, err;
+
+    if( e->xmin ) {	/* not all NaNs */
+      t->selectedcomp = comp;
+      VecCopy(xtmp, e->xmin);
+      ftmp = FindMinimum(t, bounds, xtmp, e->fmin);
+      if( ftmp < r->fmin ) {
+        r->fmin = ftmp;
+        VecCopy(r->xmin, xtmp);
+      }
+
+      t->selectedcomp = Tag(comp);
+      VecCopy(xtmp, e->xmax);
+      ftmp = -FindMinimum(t, bounds, xtmp, -e->fmax);
+      if( ftmp > r->fmax ) {
+        r->fmax = ftmp;
+        VecCopy(r->xmax, xtmp);
+      }
+    }
+
+    r->spread = halfvol*(r->fmax - r->fmin);
+    err = r->spread/Max(fabs(r->avg), NOTZERO);
+    if( err > maxerr ) {
+      maxerr = err;
+      maxcomp = comp;
+    }
+  }
+
+  t->neval_opt += t->neval;
+
+  if( maxcomp == -1 ) { /* all NaNs */
+    region->depth = 0;
+    return -1;
+  }
+
+  region->cutcomp = maxcomp;
+  r = &region->result[maxcomp];
+  if( halfvol*(r->fmin + r->fmax) > r->avg ) {
+    region->fminor = r->fmin;
+    region->fmajor = r->fmax;
+    region->xmajor = r->xmax - (real *)region->result;
+  }
+  else {
+    region->fminor = r->fmax;
+    region->fmajor = r->fmin;
+    region->xmajor = r->xmin - (real *)region->result;
+  }
+
+  if( region->isamples == 0 ) {
+    if( r->spread < samples->neff*r->err ||
+        r->spread < t->totals[maxcomp].secondspread ) region->depth = 0;
+    if( region->depth == 0 )
+      for( comp = 0; comp < t->ncomp; ++comp )
+        t->totals[comp].secondspread =
+          Max(t->totals[comp].secondspread, result[comp].spread);
+  }
+
+  if( region->depth ) Split(t, iregion);
+
+  return iregion;
+}
+