findminimum.c

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

/*
	FindMinimum.c
		find minimum (maximum) of hyperrectangular region
		this file is part of Divonne
		last modified 7 Mar 05 th
*/


#define EPS 0x1p-52
#define RTEPS 0x1p-26
#define QEPS 0x1p-13

#define DELTA 0x1p-16
#define RTDELTA 0x1p-8
#define QDELTA 0x1p-4

/*
#define DELTA 1e-5
#define RTDELTA 3.1622776601683791e-3
#define QDELTA 5.6234132519034912e-2
*/

#define SUFTOL 8*QEPS*QDELTA
#define FTOL 5e-2
#define GTOL 1e-2

#define Hessian(i, j) hessian[(i)*ndim_ + j]

#define Tag(x) ((x) | 0x8000)
#define Untag(x) ((x) & 0x7fff)
#define TaggedQ(x) ((x) & 0x8000)

typedef struct { real dx, f; } Point;

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

static inline real SignSample(real *x)
{
  return sign_*Sample(x);
}

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

static inline real Dot(ccount n, creal *a, creal *b)
{
  real sum = 0;
  count i;
  for( i = 0; i < n; ++i ) sum += a[i]*b[i];
  return sum;
}

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

static inline real Length(ccount n, creal *vec)
{
  return sqrt(Dot(n, vec, vec));
}

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

static inline void LinearSolve(ccount n, creal *hessian,
  creal *grad, real *p)
{
  int i, j;
  real dir;

  for( i = 0; i < n; ++i ) {
    dir = -grad[i];
    for( j = 0; j < i; ++j )
      dir -= Hessian(i, j)*p[j];
    p[i] = dir;
  }

  while( --i >= 0 ) {
    if( Hessian(i, i) <= 0 ) return;
    dir = p[i]/Hessian(i, i);
    for( j = i + 1; j < n; ++j )
      dir -= Hessian(j, i)*p[j];
    p[i] = dir;
  }
}

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

static void RenormalizeCholesky(ccount n, real *hessian,
  real *z, real alpha)
{
  count i, j;

  for( i = 0; i < n; ++i ) {
    creal dir = z[i];
    real beta = alpha*dir;
    real gamma = Hessian(i, i);
    real gammanew = Hessian(i, i) += beta*dir;

    if( i + 1 >= n || gammanew < 0 ||
        (gammanew < 1 && gamma > DBL_MAX*gammanew) ) return;

    gamma /= gammanew;
    beta /= gammanew;
    alpha *= gamma;

    if( gamma < .25 ) {
      for( j = i + 1; j < n; ++j ) {
        real delta = beta*z[j];
        z[j] -= dir*Hessian(j, i);
        Hessian(j, i) = Hessian(j, i)*gamma + delta;
      }
    }
    else {
      for( j = i + 1; j < n; ++j ) {
        z[j] -= dir*Hessian(j, i);
        Hessian(j, i) += beta*z[j];
      }
    }
  }
}

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

static void UpdateCholesky(ccount n, real *hessian,
  real *z, real *p)
{
  int i, j;
  real gamma = 0;

  for( i = 0; i < n; ++i ) {
    real dir = z[i];
    for( j = 0; j < i; ++j )
      dir -= Hessian(i, j)*p[j];
    p[i] = dir;
    gamma += Sq(dir)/Hessian(i, i);
  }
  gamma = Max(fabs(1 - gamma), EPS);

  while( --i >= 0 ) {
    creal dir = z[i] = p[i];
    real beta = dir/Hessian(i, i);
    creal gammanew = gamma + dir*beta;
    Hessian(i, i) *= gamma/gammanew;
    beta /= gamma;
    gamma = gammanew;
    for( j = i + 1; j < n; ++j ) {
      creal delta = beta*z[j];
      z[j] += dir*Hessian(j, i);
      Hessian(j, i) -= delta;
    }
  }
}

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

static inline void BFGS(ccount n, real *hessian,
  creal *gnew, creal *g, real *p, creal dx)
{
  real y[NDIM], c;
  count i, j;

  for( i = 0; i < n; ++i )
    y[i] = gnew[i] - g[i];
  c = dx*Dot(n, y, p);
  if( c < 1e-10 ) return;
  RenormalizeCholesky(n, hessian, y, 1/c);

  c = Dot(n, g, p);
  if( c >= 0 ) return;
  c = 1/sqrt(-c);
  for( i = 0; i < n; ++i )
    y[i] = c*g[i];
  UpdateCholesky(n, hessian, y, p);

  for( i = 0; i < n - 1; ++i )
    for( j = i + 1; j < n; ++j )
      Hessian(i, j) = Hessian(j, i);
}

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

static void Gradient(ccount nfree, ccount *ifree,
  cBounds *b, real *x, creal y, real *grad)
{
  count i;

  for( i = 0; i < nfree; ++i ) {
    ccount dim = Untag(ifree[i]);
    creal xd = x[dim];
    creal delta = (b[dim].upper - xd < DELTA) ? -DELTA : DELTA;
    x[dim] += delta;
    grad[i] = (SignSample(x) - y)/delta;
    x[dim] = xd;
  }
}

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

static Point LineSearch(ccount nfree, ccount *ifree,
  creal *p, creal *xini, real fini, real *x,
  real step, creal range, creal grad,
  creal ftol, creal xtol, creal gtol)
{
  real tol = ftol, tol2 = tol + tol;
  Point cur = {0, fini};

  VecCopy(x, xini);

  /* don't even try if
     a) we'd walk backwards,
     b) the range to explore is too small,
     c) the gradient is positive, i.e. we'd move uphill */

  if( step > 0 && range > tol2 && grad <= 0 ) {
    creal eps = RTEPS*fabs(range) + ftol;
    creal mingrad = -1e-4*grad, maxgrad = -gtol*grad;

    real end = range + eps;
    real maxstep = range - eps/(1 + RTEPS);

    Point min = cur, v = cur, w = cur;
    Point a = cur, b = {end, 0};
    real a1, b1 = end;

    /* distmin: distance along p from xini to the minimum,
       u: second-lowest point,
       v: third-lowest point,
       a, b: interval in which the minimum is sought. */

    real distmin = 0, dist, mid, q, r, s;
    count i;
    int shift;
    bool first;

    for( first = true; ; first = false ) {
      if( step >= maxstep ) {
        step = maxstep;
        maxstep = maxstep*(1 + .75*RTEPS) + .75*tol;
      }

      cur.dx = (fabs(step) >= tol) ? step : (step > 0) ? tol : -tol;
      dist = distmin + cur.dx;
      for( i = 0; i < nfree; ++i ) {
        ccount dim = ifree[i];
        x[dim] = xini[dim] + dist*p[i];
      }
      cur.f = SignSample(x);

      if( cur.f <= min.f ) {
        v = w;
        w = min;
        min.f = cur.f;
        distmin = dist;

        /* shift everything to the new minimum position */
        maxstep -= cur.dx;
        v.dx -= cur.dx;
        w.dx -= cur.dx;
        a.dx -= cur.dx;
        b.dx -= cur.dx;
        if( cur.dx < 0 ) b = w;
        else a = w;

        tol = RTEPS*fabs(distmin) + ftol;
        tol2 = tol + tol;
      }
      else {
        if( cur.dx < 0 ) a = cur;
        else b = cur;
        if( cur.f <= w.f || w.dx == 0 ) v = w, w = cur;
        else if( cur.f <= v.f || v.dx == 0 || v.dx == w.dx ) v = cur;
      }

      if( distmin + b.dx <= xtol ) break;
      if( min.f < fini &&
          a.f - min.f <= fabs(a.dx)*maxgrad &&
          (fabs(distmin - range) > tol || maxstep < b.dx) ) break;

      mid = .5*(a.dx + b.dx);
      if( fabs(mid) <= tol2 - .5*(b.dx - a.dx) ) break;

      r = q = s = 0;
      if( fabs(end) > tol ) {
        if( first ) {
          creal s1 = w.dx*grad;
          creal s2 = w.f - min.f;
          s = (s1 - ((distmin == 0) ? 0 : 2*s2))*w.dx;
          q = 2*(s2 - s1);
        }
        else {
          creal s1 = w.dx*(v.f - min.f);
          creal s2 = v.dx*(w.f - min.f);
          s = s1*w.dx - s2*v.dx;
          q = 2*(s2 - s1);
        }
        if( q > 0 ) s = -s;
        q = fabs(q);
        r = end;
        if( step != b1 || b.dx <= maxstep ) end = step;
      }

      if( distmin == a.dx ) step = mid;
      else if( b.dx > maxstep ) step = (step < b.dx) ? -4*a.dx : maxstep;
      else {
        real num = a.dx, den = b.dx;
        if( fabs(b.dx) <= tol || (w.dx > 0 && fabs(a.dx) > tol) )
          num = b.dx, den = a.dx;
        num /= -den;
        step = (num < 1) ? .5*den*sqrt(num) : 5/11.*den*(.1 + 1/num);
      }

      if( step > 0 ) a1 = a.dx, b1 = step;
      else a1 = step, b1 = b.dx;
      if( fabs(s) < fabs(.5*q*r) && s > q*a1 && s < q*b1 ) {
        step = s/q;
        if( step - a.dx < tol2 || b.dx - step < tol2 )
          step = (mid > 0) ? tol : -tol;
      }
      else end = (mid > 0) ? b.dx : a.dx;
    }

    first = true;
    if( fabs(distmin - range) < tol ) {
      distmin = range;
      if( maxstep > b.dx ) first = false;
    }

    for( cur.dx = distmin, cur.f = min.f, shift = -1; ;
         cur.dx = Max(ldexp(distmin, shift), ftol), shift <<= 1 ) {
      for( i = 0; i < nfree; ++i ) {
        ccount dim = ifree[i];
        x[dim] = xini[dim] + cur.dx*p[i];
      }
      if( !first ) cur.f = SignSample(x);

      if( cur.dx + b.dx <= xtol ) {
        cur.dx = 0;
        break;
      }
      if( fini - cur.f > cur.dx*mingrad ) break;
      if( cur.dx <= ftol ) {
        cur.dx = 0;
        break;
      }
      first = false;
    }
  }

  return cur;
}

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

static real LocalSearch(ccount nfree, ccount *ifree, cBounds *b,