C:/programs/UNCMIN/src/include/Fmin.h

Go to the documentation of this file.
00001 /*! \file Fmin.h
00002  
00003   \brief
00004   Univariate function minimization by Brent's method.
00005   Release 080309.
00006 
00007   A translation of Fmin.java by Steve Verrill to C++ 
00008   Version 0.000318 (March 18, 2000)
00009   by Brad Hanson (http://www.b-a-h.com/index.html)
00010   
00011   This software is available at 
00012   http://www.smallwaters.com/software/cpp/uncmin.html
00013   
00014   Fmin.java can be obtained from http://www1.fpl.fs.fed.us/optimization.html
00015   Information about Fmin.java is given below.
00016   
00017     Fmin.java copyright claim:
00018     This software is based on the public domain fmin routine.
00019     The FORTRAN version can be found at
00020     www.netlib.org
00021     This software was translated from the FORTRAN version
00022     to Java by a US government employee on official time.  
00023     Thus this software is also in the public domain.
00024     The translator's mail address is:
00025     Steve Verrill 
00026     USDA Forest Products Laboratory
00027     1 Gifford Pinchot Drive
00028     Madison, Wisconsin
00029     53705
00030     The translator's e-mail address is:
00031     steve@ws13.fpl.fs.fed.us
00032     
00033     (c) 2008, Werner Wothke, Fmin documentation.
00034     
00035   ***********************************************************************
00036   DISCLAIMER OF WARRANTIES:
00037   THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND. 
00038   THE TRANSLATOR DOES NOT WARRANT, GUARANTEE OR MAKE ANY REPRESENTATIONS 
00039   REGARDING THE SOFTWARE OR DOCUMENTATION IN TERMS OF THEIR CORRECTNESS, 
00040   RELIABILITY, CURRENTNESS, OR OTHERWISE. THE ENTIRE RISK AS TO 
00041   THE RESULTS AND PERFORMANCE OF THE SOFTWARE IS ASSUMED BY YOU. 
00042   IN NO CASE WILL ANY PARTY INVOLVED WITH THE CREATION OR DISTRIBUTION 
00043   OF THE SOFTWARE BE LIABLE FOR ANY DAMAGE THAT MAY RESULT FROM THE USE 
00044   OF THIS SOFTWARE.
00045   Sorry about that.
00046   ***********************************************************************
00047   History:
00048   Date        Translator        Changes
00049   3/09/08     Werner Wothke     Doxygen-compatible documentation.
00050   
00051   3/24/98     Steve Verrill     Translated
00052    */
00053   /**
00054    *
00055    *<p>
00056    *This class was translated by a statistician from the FORTRAN 
00057    *version of fmin.  It is NOT an official translation.  When 
00058    *public domain Java optimization routines become available 
00059    *from professional numerical analysts, then <b>THE CODE PRODUCED
00060    *BY THE NUMERICAL ANALYSTS SHOULD BE USED</b>.
00061    *
00062    *<p>
00063    *Meanwhile, if you have suggestions for improving this
00064    *code, please contact Steve Verrill at steve@ws13.fpl.fs.fed.us.
00065    *
00066    *@author Steve Verrill
00067    *@version .5 --- March 24, 1998
00068    * 
00069    */
00070 #include <limits>
00071 #include <cmath>
00072 /*!
00073   \brief 
00074   This function performs a uni-dimensional minimization by Brent's method.
00075   
00076   Brent's method combines a golden-section search and parabolic interpolation.  
00077   The following introductory comments are copied from the FORTRAN version.
00078   
00079   Convergence is never much slower than that for a Fibonacci search.  If f 
00080   has a continuous second derivative which is positive at the minimum (which 
00081   is not at ax or bx), then convergence is superlinear, and usually of the 
00082   order of about 1.324.
00083   
00084   The function f is never evaluated at two points closer together
00085   than eps*abs(fmin)+(tol/3), where eps is approximately the square
00086   root of the relative machine precision.  If f is a unimodal
00087   function and the computed values of f are always unimodal when
00088   separated by at least eps*abs(x)+(tol/3), then fmin approximates
00089   the abcissa of the global minimum of f on the interval (ax,bx) with
00090   an error less than 3*eps*abs(fmin)+tol.  If f is not unimodal,
00091   then fmin may approximate a local, but perhaps non-global, minimum to
00092   the same accuracy.
00093   
00094   This function subprogram is a slightly modified version of the
00095   Algol 60 procedure localmin given in Richard Brent, Algorithms For
00096   Minimization Without Derivatives, Prentice-Hall, Inc. (1973).
00097   
00098    This method is a translation from FORTRAN to Java of the Netlib 
00099    function fmin.  In the Netlib listing no author is given.
00100    
00101    Translated by Steve Verrill, March 24, 1998.
00102    
00103    Translated from Java to C++ by Brad Hanson (bradh@pobox.com)
00104         
00105   \section template_args Template Arguments
00106    
00107   \param T    Floating point type (e.g., double, float)
00108   \param FUNC Function object type giving function to minimize.
00109 
00110   \section func_args Function Arguments
00111    
00112   \param[in]  a   Left endpoint of initial interval
00113   \param[in]  b   Right endpoint of initial interval
00114   \param[in]  f   A class that overloads the function call operator
00115       [operator()()] to implement the function f(x) for any x
00116       in the interval (a,b)
00117   \param[in]  tol Desired length of the interval in which
00118       the minimum will be determined to lie
00119       (This should be greater than, roughly, 3.0e-8.)
00120  */
00121 template<class T, class FUNC>
00122 T Fmin(T a, T b, FUNC &f, T tol) 
00123 {
00124       T c,d,e,eps,xm,p,q,r,tol1,t2,
00125              u,v,w,fu,fv,fw,fx,x,tol3;
00126       c = .5*(3.0 - std::sqrt(5.0));
00127       d = 0.0;
00128       // 1.1102e-16 is machine precision
00129       eps = std::numeric_limits<T>::epsilon();
00130       tol1 = eps + 1.0;
00131       eps = std::sqrt(eps);
00132       v = a + c*(b-a);
00133       w = v;
00134       x = v;
00135       e = 0.0;
00136       fx = f(x);
00137       fv = fx;
00138       fw = fx;
00139       tol3 = tol/3.0;
00140       xm = .5*(a + b);
00141       tol1 = eps*std::fabs(x) + tol3;
00142       t2 = 2.0*tol1;
00143       // main loop
00144       while (std::fabs(x-xm) > (t2 - .5*(b-a))) {
00145          p = q = r = 0.0;
00146          if (std::fabs(e) > tol1) {
00147            // fit the parabola
00148             r = (x-w)*(fx-fv);
00149             q = (x-v)*(fx-fw);
00150             p = (x-v)*q - (x-w)*r;
00151             q = 2.0*(q-r);
00152             if (q > 0.0) {
00153                p = -p;
00154             } else {
00155                q = -q;
00156             }
00157             r = e;
00158             e = d;         
00159             // brace below corresponds to statement 50
00160          }
00161          if ((std::fabs(p) < std::fabs(.5*q*r)) &&
00162              (p > q*(a-x)) &&
00163              (p < q*(b-x))) {
00164            // a parabolic interpolation step
00165             d = p/q;
00166             u = x+d;
00167             // f must not be evaluated too close to a or b
00168             if (((u-a) < t2) || ((b-u) < t2)) {
00169                d = tol1;
00170                if (x >= xm) d = -d;
00171             }
00172             // brace below corresponds to statement 60
00173          } else {
00174            // a golden-section step
00175             if (x < xm) {
00176                e = b-x;
00177             } else {
00178                e = a-x;
00179             }
00180             d = c*e;
00181          }
00182          // f must not be evaluated too close to x
00183          if (std::fabs(d) >= tol1) {
00184             u = x+d;
00185          } else {
00186             if (d > 0.0) {
00187                u = x + tol1;
00188             } else {
00189                u = x - tol1;
00190             }
00191          }
00192          fu = f(u);
00193          // Update a, b, v, w, and x
00194          if (fx <= fu) {
00195             if (u < x) {
00196                a = u;
00197             } else {
00198                b = u;
00199             }
00200             // brace below corresponds to statement 140
00201          }
00202          if (fu <= fx) {
00203             if (u < x) {
00204                b = x;
00205             } else {
00206                a = x;
00207             }
00208             v = w;
00209             fv = fw;
00210             w = x;
00211             fw = fx;
00212             x = u;
00213             fx = fu;
00214             xm = .5*(a + b);
00215             tol1 = eps*std::fabs(x) + tol3;
00216             t2 = 2.0*tol1;
00217             // brace below corresponds to statement 170
00218          } else {
00219             if ((fu <= fw) || (w == x)) {
00220                v = w;
00221                fv = fw;
00222                w = u;
00223                fw = fu;
00224                xm = .5*(a + b);
00225                tol1 = eps*std::fabs(x) + tol3;
00226                t2 = 2.0*tol1;
00227             } else if ((fu > fv) && (v != x) && (v != w)) {
00228                xm = .5*(a + b);
00229                tol1 = eps*std::fabs(x) + tol3;
00230                t2 = 2.0*tol1;
00231             } else {
00232                v = u;
00233                fv = fu;
00234                xm = .5*(a + b);
00235                tol1 = eps*std::fabs(x) + tol3;
00236                t2 = 2.0*tol1;
00237             }
00238          }
00239          // brace below corresponds to statement 190
00240       }
00241       return x;
00242 }   
00243                                     

Generated on Sun Mar 9 14:31:58 2008 for UNCMIN by  doxygen 1.5.4