findminimum.c

The CUBA library provides new implementation of four general-purpose multidimensional integration al

  creal *x, creal fx, real *z)
{
  real delta, smax, sopp, spmax, snmax;
  real y[NDIM], fy, fz, ftest;
  real p[NDIM];
  int sign;
  count i;

  /* Choose a direction p along which to move away from the
     present x.  We choose the direction which leads farthest
     away from all borders. */

  smax = INFTY;
  for( i = 0; i < nfree; ++i ) {
    ccount dim = ifree[i];
    creal sp = b[dim].upper - x[dim];
    creal sn = x[dim] - b[dim].lower;
    if( sp < sn ) {
      smax = Min(smax, sn);
      p[i] = -1;
    }
    else {
      smax = Min(smax, sp);
      p[i] = 1;
    }
  }
  smax *= .9;

  /* Move along p until the integrand changes appreciably
     or we come close to a border. */

  VecCopy(y, x);
  ftest = SUFTOL*(1 + fabs(fx));
  delta = RTDELTA/5;
  do {
    delta = Min(5*delta, smax);
    for( i = 0; i < nfree; ++i ) {
      ccount dim = ifree[i];
      y[dim] = x[dim] + delta*p[i];
    }
    fy = SignSample(y);
    if( fabs(fy - fx) > ftest ) break;
  } while( delta != smax );

  /* Construct a second direction p' orthogonal to p, i.e. p.p' = 0.
     We let pairs of coordinates cancel in the dot product,
     i.e. we choose p'[0] = p[0], p'[1] = -p[1], etc.
     (It should really be 1/p and -1/p, but p consists of 1's and -1's.)
     For odd nfree, we let the last three components cancel by 
     choosing p'[nfree - 3] = p[nfree - 3],
              p'[nfree - 2] = -1/2 p[nfree - 2], and
              p'[nfree - 1] = -1/2 p[nfree - 1]. */

  sign = (nfree <= 1 && fy > fx) ? 1 : -1;
  spmax = snmax = INFTY;
  for( i = 0; i < nfree; ++i ) {
    ccount dim = ifree[i];
    real sp, sn;
    p[i] *= (nfree & 1 && nfree - i <= 2) ? -.5*sign : (sign = -sign);
    sp = (b[dim].upper - y[dim])/p[i];
    sn = (y[dim] - b[dim].lower)/p[i];
    if( p[i] > 0 ) {
      spmax = Min(spmax, sp);
      snmax = Min(snmax, sn);
    }
    else {
      spmax = Min(spmax, -sn);
      snmax = Min(snmax, -sp);
    }
  }
  smax = .9*spmax;
  sopp = .9*snmax;

  if( nfree > 1 && smax < snmax ) {
    real tmp = smax;
    smax = sopp;
    sopp = tmp;
    for( i = 0; i < nfree; ++i )
      p[i] = -p[i];
  }

  /* Move along p' until the integrand changes appreciably
     or we come close to a border. */

  VecCopy(z, y);
  ftest = SUFTOL*(1 + fabs(fy));
  delta = RTDELTA/5;
  do {
    delta = Min(5*delta, smax);
    for( i = 0; i < nfree; ++i ) {
      ccount dim = ifree[i];
      z[dim] = y[dim] + delta*p[i];
    }
    fz = SignSample(z);
    if( fabs(fz - fy) > ftest ) break;
  } while( delta != smax );

  if( fy != fz ) {
    real pleneps, grad, range, step;
    Point low;

    if( fy > fz ) {
      grad = (fz - fy)/delta;
      range = smax/.9;
      step = Min(delta + delta, smax);
    }
    else {
      grad = (fy - fz)/delta;
      range = sopp/.9 + delta;
      step = Min(delta + delta, sopp);
      VecCopy(y, z);
      fy = fz;
      for( i = 0; i < nfree; ++i )
        p[i] = -p[i];
    }

    pleneps = Length(nfree, p) + RTEPS;
    low = LineSearch(nfree, ifree, p, y, fy, z, step, range, grad,
      RTEPS/pleneps, 0., RTEPS);
    fz = low.f;
  }

  if( fz != fx ) {
    real pleneps, grad, range, step;
    Point low;

    spmax = snmax = INFTY;
    for( i = 0; i < nfree; ++i ) {
      ccount dim = ifree[i];
      p[i] = z[dim] - x[dim];
      if( p[i] != 0 ) {
        creal sp = (b[dim].upper - x[dim])/p[i];
        creal sn = (x[dim] - b[dim].lower)/p[i];
        if( p[i] > 0 ) {
          spmax = Min(spmax, sp);
          snmax = Min(snmax, sn);
        }
        else {
          spmax = Min(spmax, -sn);
          snmax = Min(snmax, -sp);
        }
      }
    }

    grad = fz - fx;
    range = spmax;
    step = Min(.9*spmax, 2.);
    pleneps = Length(nfree, p) + RTEPS;
    if( fz > fx ) {
      delta = Min(.9*snmax, RTDELTA/pleneps);
      for( i = 0; i < nfree; ++i ) {
        ccount dim = ifree[i];
        z[dim] = x[dim] - delta*p[i];
      }
      fz = SignSample(z);
      if( fz < fx ) {
        grad = (fz - fx)/delta;
        range = snmax;
        step = Min(.9*snmax, delta + delta);
        for( i = 0; i < nfree; ++i )
          p[i] = -p[i];
      }
      else if( delta < 1 ) grad = (fx - fz)/delta;
    }

    low = LineSearch(nfree, ifree, p, x, fx, z, step, range, grad,
      RTEPS/pleneps, 0., RTEPS);
    fz = low.f;
  }

  return fz;
}

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

static real FindMinimum(cBounds *b, real *xmin, real fmin)
{
  real hessian[NDIM*NDIM];
  real gfree[NDIM], p[NDIM];
  real tmp[NDIM], ftmp, fini = fmin;
  ccount maxeval = neval_ + 50*ndim_;
  count nfree, nfix;
  count ifree[NDIM], ifix[NDIM];
  count dim, local;

  Zap(hessian);
  for( dim = 0; dim < ndim_; ++dim )
    Hessian(dim, dim) = 1;

  /* Step 1: - classify the variables as "fixed" (sufficiently close
               to a border) and "free",
             - if the integrand is flat in the direction of the gradient
               w.r.t. the free dimensions, perform a local search. */

  for( local = 0; local < 2; ++local ) {
    bool resample = false;
    nfree = nfix = 0;
    for( dim = 0; dim < ndim_; ++dim ) {
      if( xmin[dim] < b[dim].lower + (1 + fabs(b[dim].lower))*QEPS ) {
        xmin[dim] = b[dim].lower;
        ifix[nfix++] = dim;
        resample = true;
      }
      else if( xmin[dim] > b[dim].upper - (1 + fabs(b[dim].upper))*QEPS ) {
        xmin[dim] = b[dim].upper;
        ifix[nfix++] = Tag(dim);
        resample = true;
      }
      else ifree[nfree++] = dim;
    }

    if( resample ) fini = fmin = SignSample(xmin);

    if( nfree == 0 ) goto releasebounds;

    Gradient(nfree, ifree, b, xmin, fmin, gfree);
    if( local || Length(nfree, gfree) > GTOL ) break;

    ftmp = LocalSearch(nfree, ifree, b, xmin, fmin, tmp);
    if( ftmp > fmin - (1 + fabs(fmin))*RTEPS )
      goto releasebounds;
    fmin = ftmp;
    VecCopy(xmin, tmp);
  }

  while( neval_ <= maxeval ) {

    /* Step 2a: perform a quasi-Newton iteration on the free
                variables only. */

    if( nfree > 0 ) {
      real plen, pleneps;
      real minstep;
      count i, mini, minfix;
      Point low;

      LinearSolve(nfree, hessian, gfree, p);
      plen = Length(nfree, p);
      pleneps = plen + RTEPS;

      minstep = INFTY;
      for( i = 0; i < nfree; ++i ) {
        count dim = Untag(ifree[i]);
        if( fabs(p[i]) > EPS ) {
          real step;
          count fix;
          if( p[i] < 0 ) {
            step = (b[dim].lower - xmin[dim])/p[i];
            fix = dim;
          }
          else {
            step = (b[dim].upper - xmin[dim])/p[i];
            fix = Tag(dim);
          }
          if( step < minstep ) {
            minstep = step;
            mini = i;
            minfix = fix;
          }
        }
      }

      if( minstep*pleneps <= DELTA ) {
fixbound:
        ifix[nfix++] = minfix;

        if( mini < --nfree ) {
          creal diag = Hessian(mini, mini);

          Clear(tmp, mini);
          for( i = mini; i < nfree; ++i )
            tmp[i] = Hessian(i + 1, mini);

          for( i = mini; i < nfree; ++i ) {
            Copy(&Hessian(i, 0), &Hessian(i + 1, 0), i);
            Hessian(i, i) = Hessian(i + 1, i + 1);
          }
          RenormalizeCholesky(nfree, hessian, tmp, diag);

          Copy(&ifree[mini], &ifree[mini + 1], nfree - mini);
          Copy(&gfree[mini], &gfree[mini + 1], nfree - mini);
        }
        continue;
      }

      low = LineSearch(nfree, ifree, p, xmin, fmin, tmp,
        Min(minstep, 1.), Min(minstep, 100.), Dot(nfree, gfree, p),
        RTEPS/pleneps, DELTA/pleneps, .2);

      if( low.dx > 0 ) {
        real fdiff;

        fmin = low.f;
        VecCopy(xmin, tmp);

        Gradient(nfree, ifree, b, xmin, fmin, tmp);
        BFGS(nfree, hessian, tmp, gfree, p, low.dx);
        VecCopy(gfree, tmp);

        if( fabs(low.dx - minstep) < QEPS*minstep ) goto fixbound;

        fdiff = fini - fmin;
        fini = fmin;
        if( fdiff > (1 + fabs(fmin))*FTOL ||
            low.dx*plen > (1 + Length(ndim_, xmin))*FTOL ) continue;
      }
    }

    /* Step 2b: check whether freeing any fixed variable will lead
                to a reduction in f. */

releasebounds:
    if( nfix > 0 ) {
      real mingrad = INFTY;
      count i, mini = 0;
      bool repeat = false;

      Gradient(nfix, ifix, b, xmin, fmin, tmp);

      for( i = 0; i < nfix; ++i ) {
        creal grad = TaggedQ(ifix[i]) ? -tmp[i] : tmp[i];
        if( grad < -RTEPS ) {
          repeat = true;
          if( grad < mingrad ) {
            mingrad = grad;
            mini = i;
          }
        }
      }

      if( repeat ) {
        gfree[nfree] = tmp[mini];
        ifree[nfree] = Untag(ifix[mini]);
        Clear(&Hessian(nfree, 0), nfree);
        Hessian(nfree, nfree) = 1;
        ++nfree;

        --nfix;
        Copy(&ifix[mini], &ifix[mini + 1], nfix - mini);
        continue;
      }
    }

    break;
  }

  return fmin;
}