FunctionR1toR1.cpp

Go to the documentation of this file.

//            QMcBeaver
//
//         Constructed by 
//
//     Michael Todd Feldmann 
//              and 
//   David Randall "Chip" Kent IV
//
// Copyright 2000.  All rights reserved.
//
// drkent@users.sourceforge.net mtfeldmann@users.sourceforge.net

#include "FunctionR1toR1.h"
#include <math.h>

using namespace std;

/*
 *  R : A Computer Language for Statistical Data Analysis
 
 *  Copyright (C) 1997-2001  The R Development Core Team
 *  Copyright (C) 2003-2004  The R Foundation
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
 */

/* fmin.f -- translated by f2c (version 19990503).
 */

/* R's  optimize() :   functionfmin(ax,bx,f,tol)
   =    ==========~~~~~~~~~~~~~~~~~

        an approximation  x  to the point where  f  attains a minimum  on
    the interval  (ax,bx)  is determined.

    INPUT..

    ax    left endpoint of initial interval
    bx    right endpoint of initial interval
    f     function which evaluates  f(x, info)  for any  x
          in the interval  (ax,bx)
    tol   desired length of the interval of uncertainty of the final
          result ( >= 0.)

    OUTPUT..

    fmin  abcissa approximating the point where  f  attains a minimum

        The method used is a combination of  golden  section  search  and
    successive parabolic interpolation.  convergence is never much slower
    than  that  for  a  Fibonacci search.  If  f  has a continuous second
    derivative which is positive at the minimum (which is not  at  ax  or
    bx),  then  convergence  is  superlinear, and usually of the order of
    about  1.324....
        The function  f  is never evaluated at two points closer together
    than  eps*abs(fmin)+(tol/3), where eps is  approximately  the  square
    root  of  the  relative  machine  precision.   if   f   is a unimodal
    function and the computed values of   f   are  always  unimodal  when
    separated  by  at least  eps*abs(x)+(tol/3), then  fmin  approximates
    the abcissa of the global minimum of  f  on the interval  ax,bx  with
    an error less than  3*eps*abs(fmin)+tol.  if   f   is  not  unimodal,
    then fmin may approximate a local, but perhaps non-global, minimum to
    the same accuracy.
        This function subprogram is a slightly modified  version  of  the
    Algol  60 procedure  localmin  given in Richard Brent, Algorithms for
    Minimization without Derivatives, Prentice-Hall, Inc. (1973).
*/

double FunctionR1toR1::Brent_fmin(double ax, double bx,
                  double tol)
{
  /*  c is the squared inverse of the golden ratio */
  const double c = (3. - sqrt(5.)) * .5;

  /* Local variables */
  double a, b, d, e, p, q, r, u, v, w, x;
  double t2, fu, fv, fw, fx, xm, eps, tol1, tol3;
  int maxIters = 100000;

  /*  eps is approximately the square root of the relative machine precision. */
  eps = 1e-50;
  tol1 = eps + 1.;/* the smallest 1.000... > 1 */
  eps = sqrt(eps);

  a = ax;
  b = bx;
  v = a + c * (b - a);
  w = v;
  x = v;

  d = 0.;/* -Wall */
  e = 0.;
  fx = function(x);
  fv = fx;
  fw = fx;
  tol3 = tol / 3.;

  /*  main loop starts here ----------------------------------- */
  int iter = 0;
  for(;;) {
    iter++;
    xm = (a + b) * .5;
    tol1 = eps * fabs(x) + tol3;
    t2 = tol1 * 2.;

    /* check stopping criterion */

    if (fabs(x - xm) <= t2 - (b - a) * .5 || iter > maxIters) break;
    p = 0.;
    q = 0.;
    r = 0.;
    if (fabs(e) > tol1) { /* fit parabola */

      r = (x - w) * (fx - fv);
      q = (x - v) * (fx - fw);
      p = (x - v) * q - (x - w) * r;
      q = (q - r) * 2.;
      if (q > 0.) p = -p; else q = -q;
      r = e;
      e = d;
    }

    if (fabs(p) >= fabs(q * .5 * r) ||
        p <= q * (a - x) || p >= q * (b - x)) { /* a golden-section step */

      if (x < xm) e = b - x; else e = a - x;
      d = c * e;
    }
    else { /* a parabolic-interpolation step */

      d = p / q;
      u = x + d;

      /* f must not be evaluated too close to ax or bx */

      if (u - a < t2 || b - u < t2) {
        d = tol1;
        if (x >= xm) d = -d;
      }
    }

    /* f must not be evaluated too close to x */

    if (fabs(d) >= tol1)
      u = x + d;
    else if (d > 0.)
      u = x + tol1;
    else
      u = x - tol1;

    fu = function(u);

    //cout << "u = " << u << " fu = " << fu << endl;

    /*  update  a, b, v, w, and x */

    if (x <= fu) {
      if (u < x)
        a = u;
      else
        b = u;
    }
    if (fu <= fx) {
      if (u < x)
        b = x;
      else
        a = x;
      v = w;    w = x;   x = u;
      fv = fw; fw = fx; fx = fu;
    }
    else {
      if (fu <= fw || w == x) {
        v = w; fv = fw;
        w = u; fw = fu;
      }
      else if (fu <= fv || v == x || v == w) {
        v = u; fv = fu;
      }
    }
  }

  /* end of main loop */

  if(iter > maxIters)
    {
      clog << "Warning: maximum iterations reached in Brent minimzation" << endl;
      clog << " iter = " << iter << endl;
      clog << " tol  = " << tol << endl;
      clog << " fabs(x - xm)      = " << fabs(x - xm) << endl;
      clog << " is supposed to be <= " << endl;
      clog << " t2 - (b - a) * .5 = " << t2 - (b - a) * .5 << endl;
    }

  return x;
}

double FunctionR1toR1::function(double x)
{
  evaluate(x);
  return getFunctionValue();
}

double FunctionR1toR1::minimum(double left, double right)
{
  cout << "left = " << left << " right = " << right << endl;
  //return 0.0;
  return Brent_fmin(left, right, 1e-50);
}

---