# brent.cpp

Go to the documentation of this file.
```00001 # include <cstdlib>
00002 # include <iostream>
00003 # include <cmath>
00004 # include <ctime>
00005
00006 using namespace std;
00007
00008 # include "brent.h"
00009
00010 namespace shogun
00011 {
00012
00013 //****************************************************************************80
00014
00015 double glomin ( double a, double b, double c, double m, double e, double t,
00016   func_base& f, double &x )
00017
00018 //****************************************************************************80
00019 //
00020 //  Purpose:
00021 //
00022 //    GLOMIN seeks a global minimum of a function F(X) in an interval [A,B].
00023 //
00024 //  Discussion:
00025 //
00026 //    This function assumes that F(X) is twice continuously differentiable
00027 //    over [A,B] and that F''(X) <= M for all X in [A,B].
00028 //
00029 //  Licensing:
00030 //
00032 //
00033 //  Modified:
00034 //
00035 //    06 May 2012
00036 //
00037 //  Author:
00038 //
00039 //    Original FORTRAN77 version by Richard Brent.
00040 //    C++ version by John Burkardt.
00041 //    Modifications by John Denker.
00042 //
00043 //  Reference:
00044 //
00045 //    Richard Brent,
00046 //    Algorithms for Minimization Without Derivatives,
00047 //    Dover, 2002,
00048 //    ISBN: 0-486-41998-3,
00049 //    LC: QA402.5.B74.
00050 //
00051 //  Parameters:
00052 //
00053 //    Input, double A, B, the endpoints of the interval.
00054 //    It must be the case that A < B.
00055 //
00056 //    Input, double C, an initial guess for the global
00057 //    minimizer.  If no good guess is known, C = A or B is acceptable.
00058 //
00059 //    Input, double M, the bound on the second derivative.
00060 //
00061 //    Input, double E, a positive tolerance, a bound for the
00062 //    absolute error in the evaluation of F(X) for any X in [A,B].
00063 //
00064 //    Input, double T, a positive error tolerance.
00065 //
00066 //    Input, func_base& F, a user-supplied c++ functor whose
00067 //    global minimum is being sought.  The input and output
00068 //    of F() are of type double.
00069 //
00070 //    Output, double &X, the estimated value of the abscissa
00071 //    for which F attains its global minimum value in [A,B].
00072 //
00073 //    Output, double GLOMIN, the value F(X).
00074 //
00075 {
00076   double a0;
00077   double a2;
00078   double a3;
00079   double d0;
00080   double d1;
00081   double d2;
00082   double h;
00083   int k;
00084   double m2;
00085   double macheps;
00086   double p;
00087   double q;
00088   double qs;
00089   double r;
00090   double s;
00091   double sc;
00092   double y;
00093   double y0;
00094   double y1;
00095   double y2;
00096   double y3;
00097   double yb;
00098   double z0;
00099   double z1;
00100   double z2;
00101
00102   a0 = b;
00103   x = a0;
00104   a2 = a;
00105   y0 = f ( b );
00106   yb = y0;
00107   y2 = f ( a );
00108   y = y2;
00109
00110   if ( y0 < y )
00111   {
00112     y = y0;
00113   }
00114   else
00115   {
00116     x = a;
00117   }
00118
00119   if ( m <= 0.0 || b <= a )
00120   {
00121     return y;
00122   }
00123
00124   macheps = r8_epsilon ( );
00125
00126   m2 = 0.5 * ( 1.0 + 16.0 * macheps ) * m;
00127
00128   if ( c <= a || b <= c )
00129   {
00130     sc = 0.5 * ( a + b );
00131   }
00132   else
00133   {
00134     sc = c;
00135   }
00136
00137   y1 = f ( sc );
00138   k = 3;
00139   d0 = a2 - sc;
00140   h = 9.0 / 11.0;
00141
00142   if ( y1 < y )
00143   {
00144     x = sc;
00145     y = y1;
00146   }
00147 //
00148 //  Loop.
00149 //
00150   for ( ; ; )
00151   {
00152     d1 = a2 - a0;
00153     d2 = sc - a0;
00154     z2 = b - a2;
00155     z0 = y2 - y1;
00156     z1 = y2 - y0;
00157     r = d1 * d1 * z0 - d0 * d0 * z1;
00158     p = r;
00159     qs = 2.0 * ( d0 * z1 - d1 * z0 );
00160     q = qs;
00161
00162     if ( k < 1000000 || y2 <= y )
00163     {
00164       for ( ; ; )
00165       {
00166         if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
00167           z2 * m2 * r * ( z2 * q - r ) )
00168         {
00169           a3 = a2 + r / q;
00170           y3 = f ( a3 );
00171
00172           if ( y3 < y )
00173           {
00174             x = a3;
00175             y = y3;
00176           }
00177         }
00178         k = ( ( 1611 * k ) % 1048576 );
00179         q = 1.0;
00180         r = ( b - a ) * 0.00001 * ( double ) ( k );
00181
00182         if ( z2 <= r )
00183         {
00184           break;
00185         }
00186       }
00187     }
00188     else
00189     {
00190       k = ( ( 1611 * k ) % 1048576 );
00191       q = 1.0;
00192       r = ( b - a ) * 0.00001 * ( double ) ( k );
00193
00194       while ( r < z2 )
00195       {
00196         if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) <
00197           z2 * m2 * r * ( z2 * q - r ) )
00198         {
00199           a3 = a2 + r / q;
00200           y3 = f ( a3 );
00201
00202           if ( y3 < y )
00203           {
00204             x = a3;
00205             y = y3;
00206           }
00207         }
00208         k = ( ( 1611 * k ) % 1048576 );
00209         q = 1.0;
00210         r = ( b - a ) * 0.00001 * ( double ) ( k );
00211       }
00212     }
00213
00214     r = m2 * d0 * d1 * d2;
00215     s = sqrt ( ( ( y2 - y ) + t ) / m2 );
00216     h = 0.5 * ( 1.0 + h );
00217     p = h * ( p + 2.0 * r * s );
00218     q = q + 0.5 * qs;
00219     r = - 0.5 * ( d0 + ( z0 + 2.01 * e ) / ( d0 * m2 ) );
00220
00221     if ( r < s || d0 < 0.0 )
00222     {
00223       r = a2 + s;
00224     }
00225     else
00226     {
00227       r = a2 + r;
00228     }
00229
00230     if ( 0.0 < p * q )
00231     {
00232       a3 = a2 + p / q;
00233     }
00234     else
00235     {
00236       a3 = r;
00237     }
00238
00239     for ( ; ; )
00240     {
00241       a3 = r8_max ( a3, r );
00242
00243       if ( b <= a3 )
00244       {
00245         a3 = b;
00246         y3 = yb;
00247       }
00248       else
00249       {
00250         y3 = f ( a3 );
00251       }
00252
00253       if ( y3 < y )
00254       {
00255         x = a3;
00256         y = y3;
00257       }
00258
00259       d0 = a3 - a2;
00260
00261       if ( a3 <= r )
00262       {
00263         break;
00264       }
00265
00266       p = 2.0 * ( y2 - y3 ) / ( m * d0 );
00267
00268       if ( ( 1.0 + 9.0 * macheps ) * d0 <= r8_abs ( p ) )
00269       {
00270         break;
00271       }
00272
00273       if ( 0.5 * m2 * ( d0 * d0 + p * p ) <= ( y2 - y ) + ( y3 - y ) + 2.0 * t )
00274       {
00275         break;
00276       }
00277       a3 = 0.5 * ( a2 + a3 );
00278       h = 0.9 * h;
00279     }
00280
00281     if ( b <= a3 )
00282     {
00283       break;
00284     }
00285
00286     a0 = sc;
00287     sc = a2;
00288     a2 = a3;
00289     y0 = y1;
00290     y1 = y2;
00291     y2 = y3;
00292   }
00293
00294   return y;
00295 }
00296 //****************************************************************************80
00297
00298 double local_min ( double a, double b, double t, func_base& f,
00299   double &x )
00300
00301 //****************************************************************************80
00302 //
00303 //  Purpose:
00304 //
00305 //    LOCAL_MIN seeks a local minimum of a function F(X) in an interval [A,B].
00306 //
00307 //  Discussion:
00308 //
00309 //    The method used is a combination of golden section search and
00310 //    successive parabolic interpolation.  Convergence is never much slower
00311 //    than that for a Fibonacci search.  If F has a continuous second
00312 //    derivative which is positive at the minimum (which is not at A or
00313 //    B), then convergence is superlinear, and usually of the order of
00315 //
00316 //    The values EPS and T define a tolerance TOL = EPS * abs ( X ) + T.
00317 //    F is never evaluated at two points closer than TOL.
00318 //
00319 //    If F is a unimodal function and the computed values of F are always
00320 //    unimodal when separated by at least SQEPS * abs ( X ) + (T/3), then
00321 //    LOCAL_MIN approximates the abscissa of the global minimum of F on the
00322 //    interval [A,B] with an error less than 3*SQEPS*abs(LOCAL_MIN)+T.
00323 //
00324 //    If F is not unimodal, then LOCAL_MIN may approximate a local, but
00325 //    perhaps non-global, minimum to the same accuracy.
00326 //
00327 //  Licensing:
00328 //
00330 //
00331 //  Modified:
00332 //
00333 //    17 July 2011
00334 //
00335 //  Author:
00336 //
00337 //    Original FORTRAN77 version by Richard Brent.
00338 //    C++ version by John Burkardt.
00339 //    Modifications by John Denker.
00340 //
00341 //  Reference:
00342 //
00343 //    Richard Brent,
00344 //    Algorithms for Minimization Without Derivatives,
00345 //    Dover, 2002,
00346 //    ISBN: 0-486-41998-3,
00347 //    LC: QA402.5.B74.
00348 //
00349 //  Parameters:
00350 //
00351 //    Input, double A, B, the endpoints of the interval.
00352 //
00353 //    Input, double T, a positive absolute error tolerance.
00354 //
00355 //    Input, func_base& F, a user-supplied c++ functor whose
00356 //    local minimum is being sought.  The input and output
00357 //    of F() are of type double.
00358 //
00359 //    Output, double &X, the estimated value of an abscissa
00360 //    for which F attains a local minimum value in [A,B].
00361 //
00362 //    Output, double LOCAL_MIN, the value F(X).
00363 //
00364 {
00365   double c;
00366   double d = 0.0;
00367   double e;
00368   double eps;
00369   double fu;
00370   double fv;
00371   double fw;
00372   double fx;
00373   double m;
00374   double p;
00375   double q;
00376   double r;
00377   double sa;
00378   double sb;
00379   double t2;
00380   double tol;
00381   double u;
00382   double v;
00383   double w;
00384 //
00385 //  C is the square of the inverse of the golden ratio.
00386 //
00387   c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
00388
00389   eps = sqrt ( r8_epsilon ( ) );
00390
00391   sa = a;
00392   sb = b;
00393   x = sa + c * ( b - a );
00394   w = x;
00395   v = w;
00396   e = 0.0;
00397   fx = f ( x );
00398   fw = fx;
00399   fv = fw;
00400
00401   for ( ; ; )
00402   {
00403     m = 0.5 * ( sa + sb ) ;
00404     tol = eps * r8_abs ( x ) + t;
00405     t2 = 2.0 * tol;
00406 //
00407 //  Check the stopping criterion.
00408 //
00409     if ( r8_abs ( x - m ) <= t2 - 0.5 * ( sb - sa ) )
00410     {
00411       break;
00412     }
00413 //
00414 //  Fit a parabola.
00415 //
00416     r = 0.0;
00417     q = r;
00418     p = q;
00419
00420     if ( tol < r8_abs ( e ) )
00421     {
00422       r = ( x - w ) * ( fx - fv );
00423       q = ( x - v ) * ( fx - fw );
00424       p = ( x - v ) * q - ( x - w ) * r;
00425       q = 2.0 * ( q - r );
00426       if ( 0.0 < q )
00427       {
00428         p = - p;
00429       }
00430       q = r8_abs ( q );
00431       r = e;
00432       e = d;
00433     }
00434
00435     if ( r8_abs ( p ) < r8_abs ( 0.5 * q * r ) &&
00436          q * ( sa - x ) < p &&
00437          p < q * ( sb - x ) )
00438     {
00439 //
00440 //  Take the parabolic interpolation step.
00441 //
00442       d = p / q;
00443       u = x + d;
00444 //
00445 //  F must not be evaluated too close to A or B.
00446 //
00447       if ( ( u - sa ) < t2 || ( sb - u ) < t2 )
00448       {
00449         if ( x < m )
00450         {
00451           d = tol;
00452         }
00453         else
00454         {
00455           d = - tol;
00456         }
00457       }
00458     }
00459 //
00460 //  A golden-section step.
00461 //
00462     else
00463     {
00464       if ( x < m )
00465       {
00466         e = sb - x;
00467       }
00468       else
00469       {
00470         e = sa - x;
00471       }
00472       d = c * e;
00473     }
00474 //
00475 //  F must not be evaluated too close to X.
00476 //
00477     if ( tol <= r8_abs ( d ) )
00478     {
00479       u = x + d;
00480     }
00481     else if ( 0.0 < d )
00482     {
00483       u = x + tol;
00484     }
00485     else
00486     {
00487       u = x - tol;
00488     }
00489
00490     fu = f ( u );
00491 //
00492 //  Update A, B, V, W, and X.
00493 //
00494     if ( fu <= fx )
00495     {
00496       if ( u < x )
00497       {
00498         sb = x;
00499       }
00500       else
00501       {
00502         sa = x;
00503       }
00504       v = w;
00505       fv = fw;
00506       w = x;
00507       fw = fx;
00508       x = u;
00509       fx = fu;
00510     }
00511     else
00512     {
00513       if ( u < x )
00514       {
00515         sa = u;
00516       }
00517       else
00518       {
00519         sb = u;
00520       }
00521
00522       if ( fu <= fw || w == x )
00523       {
00524         v = w;
00525         fv = fw;
00526         w = u;
00527         fw = fu;
00528       }
00529       else if ( fu <= fv || v == x || v== w )
00530       {
00531         v = u;
00532         fv = fu;
00533       }
00534     }
00535   }
00536   return fx;
00537 }
00538 //****************************************************************************80
00539
00540 double local_min_rc ( double &a, double &b, int &status, double value )
00541
00542 //****************************************************************************80
00543 //
00544 //  Purpose:
00545 //
00546 //    LOCAL_MIN_RC seeks a minimizer of a scalar function of a scalar variable.
00547 //
00548 //  Discussion:
00549 //
00550 //    This routine seeks an approximation to the point where a function
00551 //    F attains a minimum on the interval (A,B).
00552 //
00553 //    The method used is a combination of golden section search and
00554 //    successive parabolic interpolation.  Convergence is never much
00555 //    slower than that for a Fibonacci search.  If F has a continuous
00556 //    second derivative which is positive at the minimum (which is not
00557 //    at A or B), then convergence is superlinear, and usually of the
00558 //    order of about 1.324...
00559 //
00560 //    The routine is a revised version of the Brent local minimization
00561 //    algorithm, using reverse communication.
00562 //
00563 //    It is worth stating explicitly that this routine will NOT be
00564 //    able to detect a minimizer that occurs at either initial endpoint
00565 //    A or B.  If this is a concern to the user, then the user must
00566 //    either ensure that the initial interval is larger, or to check
00567 //    the function value at the returned minimizer against the values
00568 //    at either endpoint.
00569 //
00570 //  Licensing:
00571 //
00573 //
00574 //  Modified:
00575 //
00576 //    17 July 2011
00577 //
00578 //  Author:
00579 //
00580 //    John Burkardt
00581 //
00582 //  Reference:
00583 //
00584 //    Richard Brent,
00585 //    Algorithms for Minimization Without Derivatives,
00586 //    Dover, 2002,
00587 //    ISBN: 0-486-41998-3,
00588 //    LC: QA402.5.B74.
00589 //
00590 //    David Kahaner, Cleve Moler, Steven Nash,
00591 //    Numerical Methods and Software,
00592 //    Prentice Hall, 1989,
00593 //    ISBN: 0-13-627258-4,
00594 //    LC: TA345.K34.
00595 //
00596 //  Parameters
00597 //
00598 //    Input/output, double &A, &B.  On input, the left and right
00599 //    endpoints of the initial interval.  On output, the lower and upper
00600 //    bounds for an interval containing the minimizer.  It is required
00601 //    that A < B.
00602 //
00603 //    Input/output, int &STATUS, used to communicate between
00604 //    the user and the routine.  The user only sets STATUS to zero on the first
00605 //    call, to indicate that this is a startup call.  The routine returns STATUS
00606 //    positive to request that the function be evaluated at ARG, or returns
00607 //    STATUS as 0, to indicate that the iteration is complete and that
00608 //    ARG is the estimated minimizer.
00609 //
00610 //    Input, double VALUE, the function value at ARG, as requested
00611 //    by the routine on the previous call.
00612 //
00613 //    Output, double LOCAL_MIN_RC, the currently considered point.
00614 //    On return with STATUS positive, the user is requested to evaluate the
00615 //    function at this point, and return the value in VALUE.  On return with
00616 //    STATUS zero, this is the routine's estimate for the function minimizer.
00617 //
00618 //  Local parameters:
00619 //
00620 //    C is the squared inverse of the golden ratio.
00621 //
00622 //    EPS is the square root of the relative machine precision.
00623 //
00624 {
00625   static double arg;
00626   static double c;
00627   static double d;
00628   static double e;
00629   static double eps;
00630   static double fu;
00631   static double fv;
00632   static double fw;
00633   static double fx;
00634   static double midpoint;
00635   static double p;
00636   static double q;
00637   static double r;
00638   static double tol;
00639   static double tol1;
00640   static double tol2;
00641   static double u;
00642   static double v;
00643   static double w;
00644   static double x;
00645 //
00646 //  STATUS (INPUT) = 0, startup.
00647 //
00648   if ( status == 0 )
00649   {
00650     if ( b <= a )
00651     {
00652       cout << "\n";
00653       cout << "LOCAL_MIN_RC - Fatal error!\n";
00654       cout << "  A < B is required, but\n";
00655       cout << "  A = " << a << "\n";
00656       cout << "  B = " << b << "\n";
00657       status = -1;
00658       exit ( 1 );
00659     }
00660     c = 0.5 * ( 3.0 - sqrt ( 5.0 ) );
00661
00662     eps = sqrt ( r8_epsilon ( ) );
00663     tol = r8_epsilon ( );
00664
00665     v = a + c * ( b - a );
00666     w = v;
00667     x = v;
00668     e = 0.0;
00669
00670     status = 1;
00671     arg = x;
00672
00673     return arg;
00674   }
00675 //
00676 //  STATUS (INPUT) = 1, return with initial function value of FX.
00677 //
00678   else if ( status == 1 )
00679   {
00680     fx = value;
00681     fv = fx;
00682     fw = fx;
00683   }
00684 //
00685 //  STATUS (INPUT) = 2 or more, update the data.
00686 //
00687   else if ( 2 <= status )
00688   {
00689     fu = value;
00690
00691     if ( fu <= fx )
00692     {
00693       if ( x <= u )
00694       {
00695         a = x;
00696       }
00697       else
00698       {
00699         b = x;
00700       }
00701       v = w;
00702       fv = fw;
00703       w = x;
00704       fw = fx;
00705       x = u;
00706       fx = fu;
00707     }
00708     else
00709     {
00710       if ( u < x )
00711       {
00712         a = u;
00713       }
00714       else
00715       {
00716         b = u;
00717       }
00718
00719       if ( fu <= fw || w == x )
00720       {
00721         v = w;
00722         fv = fw;
00723         w = u;
00724         fw = fu;
00725       }
00726       else if ( fu <= fv || v == x || v == w )
00727       {
00728         v = u;
00729         fv = fu;
00730       }
00731     }
00732   }
00733 //
00734 //  Take the next step.
00735 //
00736   midpoint = 0.5 * ( a + b );
00737   tol1 = eps * r8_abs ( x ) + tol / 3.0;
00738   tol2 = 2.0 * tol1;
00739 //
00740 //  If the stopping criterion is satisfied, we can exit.
00741 //
00742   if ( r8_abs ( x - midpoint ) <= ( tol2 - 0.5 * ( b - a ) ) )
00743   {
00744     status = 0;
00745     return arg;
00746   }
00747 //
00748 //  Is golden-section necessary?
00749 //
00750   if ( r8_abs ( e ) <= tol1 )
00751   {
00752     if ( midpoint <= x )
00753     {
00754       e = a - x;
00755     }
00756     else
00757     {
00758       e = b - x;
00759     }
00760     d = c * e;
00761   }
00762 //
00763 //  Consider fitting a parabola.
00764 //
00765   else
00766   {
00767     r = ( x - w ) * ( fx - fv );
00768     q = ( x - v ) * ( fx - fw );
00769     p = ( x - v ) * q - ( x - w ) * r;
00770     q = 2.0 * ( q - r );
00771     if ( 0.0 < q )
00772     {
00773       p = - p;
00774     }
00775     q = r8_abs ( q );
00776     r = e;
00777     e = d;
00778 //
00779 //  Choose a golden-section step if the parabola is not advised.
00780 //
00781     if (
00782       ( r8_abs ( 0.5 * q * r ) <= r8_abs ( p ) ) ||
00783       ( p <= q * ( a - x ) ) ||
00784       ( q * ( b - x ) <= p ) )
00785     {
00786       if ( midpoint <= x )
00787       {
00788         e = a - x;
00789       }
00790       else
00791       {
00792         e = b - x;
00793       }
00794       d = c * e;
00795     }
00796 //
00797 //  Choose a parabolic interpolation step.
00798 //
00799     else
00800     {
00801       d = p / q;
00802       u = x + d;
00803
00804       if ( ( u - a ) < tol2 )
00805       {
00806         d = tol1 * r8_sign ( midpoint - x );
00807       }
00808
00809       if ( ( b - u ) < tol2 )
00810       {
00811         d = tol1 * r8_sign ( midpoint - x );
00812       }
00813     }
00814   }
00815 //
00816 //  F must not be evaluated too close to X.
00817 //
00818   if ( tol1 <= r8_abs ( d ) )
00819   {
00820     u = x + d;
00821   }
00822   if ( r8_abs ( d ) < tol1 )
00823   {
00824     u = x + tol1 * r8_sign ( d );
00825   }
00826 //
00827 //  Request value of F(U).
00828 //
00829   arg = u;
00830   status = status + 1;
00831
00832   return arg;
00833 }
00834 //****************************************************************************80
00835
00836 double r8_abs ( double x )
00837
00838 //****************************************************************************80
00839 //
00840 //  Purpose:
00841 //
00842 //    R8_ABS returns the absolute value of an R8.
00843 //
00844 //  Licensing:
00845 //
00847 //
00848 //  Modified:
00849 //
00850 //    07 May 2006
00851 //
00852 //  Author:
00853 //
00854 //    John Burkardt
00855 //
00856 //  Parameters:
00857 //
00858 //    Input, double X, the quantity whose absolute value is desired.
00859 //
00860 //    Output, double R8_ABS, the absolute value of X.
00861 //
00862 {
00863   double value;
00864
00865   if ( 0.0 <= x )
00866   {
00867     value = x;
00868   }
00869   else
00870   {
00871     value = - x;
00872   }
00873   return value;
00874 }
00875 //****************************************************************************80
00876
00877 double r8_epsilon ( )
00878
00879 //****************************************************************************80
00880 //
00881 //  Purpose:
00882 //
00883 //    R8_EPSILON returns the R8 round off unit.
00884 //
00885 //  Discussion:
00886 //
00887 //    R8_EPSILON is a number R which is a power of 2 with the property that,
00888 //    to the precision of the computer's arithmetic,
00889 //      1 < 1 + R
00890 //    but
00891 //      1 = ( 1 + R / 2 )
00892 //
00893 //  Licensing:
00894 //
00896 //
00897 //  Modified:
00898 //
00899 //    08 May 2006
00900 //
00901 //  Author:
00902 //
00903 //    John Burkardt
00904 //
00905 //  Parameters:
00906 //
00907 //    Output, double R8_EPSILON, the double precision round-off unit.
00908 //
00909 {
00910   double r;
00911
00912   r = 1.0;
00913
00914   while ( 1.0 < ( double ) ( 1.0 + r )  )
00915   {
00916     r = r / 2.0;
00917   }
00918
00919   return ( 2.0 * r );
00920 }
00921 //****************************************************************************80
00922
00923 double r8_max ( double x, double y )
00924
00925 //****************************************************************************80
00926 //
00927 //  Purpose:
00928 //
00929 //    R8_MAX returns the maximum of two R8's.
00930 //
00931 //  Licensing:
00932 //
00934 //
00935 //  Modified:
00936 //
00937 //    18 August 2004
00938 //
00939 //  Author:
00940 //
00941 //    John Burkardt
00942 //
00943 //  Parameters:
00944 //
00945 //    Input, double X, Y, the quantities to compare.
00946 //
00947 //    Output, double R8_MAX, the maximum of X and Y.
00948 //
00949 {
00950   double value;
00951
00952   if ( y < x )
00953   {
00954     value = x;
00955   }
00956   else
00957   {
00958     value = y;
00959   }
00960   return value;
00961 }
00962 //****************************************************************************80
00963
00964 double r8_sign ( double x )
00965
00966 //****************************************************************************80
00967 //
00968 //  Purpose:
00969 //
00970 //    R8_SIGN returns the sign of an R8.
00971 //
00972 //  Licensing:
00973 //
00975 //
00976 //  Modified:
00977 //
00978 //    18 October 2004
00979 //
00980 //  Author:
00981 //
00982 //    John Burkardt
00983 //
00984 //  Parameters:
00985 //
00986 //    Input, double X, the number whose sign is desired.
00987 //
00988 //    Output, double R8_SIGN, the sign of X.
00989 //
00990 {
00991   double value;
00992
00993   if ( x < 0.0 )
00994   {
00995     value = -1.0;
00996   }
00997   else
00998   {
00999     value = 1.0;
01000   }
01001   return value;
01002 }
01003 //****************************************************************************80
01004
01005 void timestamp ( )
01006
01007 //****************************************************************************80
01008 //
01009 //  Purpose:
01010 //
01011 //    TIMESTAMP prints the current YMDHMS date as a time stamp.
01012 //
01013 //  Example:
01014 //
01015 //    31 May 2001 09:45:54 AM
01016 //
01017 //  Licensing:
01018 //
01020 //
01021 //  Modified:
01022 //
01023 //    24 September 2003
01024 //
01025 //  Author:
01026 //
01027 //    John Burkardt
01028 //
01029 //  Parameters:
01030 //
01031 //    None
01032 //
01033 {
01034 # define TIME_SIZE 40
01035
01036   static char time_buffer[TIME_SIZE];
01037   const struct tm *tm;
01038   time_t now;
01039
01040   now = time ( NULL );
01041   tm = localtime ( &now );
01042
01043   strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
01044
01045   cout << time_buffer << "\n";
01046
01047   return;
01048 # undef TIME_SIZE
01049 }
01050 //****************************************************************************80
01051
01052 double zero ( double a, double b, double t, func_base& f )
01053
01054 //****************************************************************************80
01055 //
01056 //  Purpose:
01057 //
01058 //    ZERO seeks the root of a function F(X) in an interval [A,B].
01059 //
01060 //  Discussion:
01061 //
01062 //    The interval [A,B] must be a change of sign interval for F.
01063 //    That is, F(A) and F(B) must be of opposite signs.  Then
01064 //    assuming that F is continuous implies the existence of at least
01065 //    one value C between A and B for which F(C) = 0.
01066 //
01067 //    The location of the zero is determined to within an accuracy
01068 //    of 6 * MACHEPS * r8_abs ( C ) + 2 * T.
01069 //
01070 //  Licensing:
01071 //
01073 //
01074 //  Modified:
01075 //
01076 //    06 May 2012
01077 //
01078 //  Author:
01079 //
01080 //    Original FORTRAN77 version by Richard Brent.
01081 //    C++ version by John Burkardt.
01082 //    Modifications by John Denker.
01083 //
01084 //  Reference:
01085 //
01086 //    Richard Brent,
01087 //    Algorithms for Minimization Without Derivatives,
01088 //    Dover, 2002,
01089 //    ISBN: 0-486-41998-3,
01090 //    LC: QA402.5.B74.
01091 //
01092 //  Parameters:
01093 //
01094 //    Input, double A, B, the endpoints of the change of sign interval.
01095 //
01096 //    Input, double T, a positive error tolerance.
01097 //
01098 //    Input, func_base& F, the name of a user-supplied c++ functor
01099 //    whose zero is being sought.  The input and output
01100 //    of F() are of type double.
01101 //
01102 //    Output, double ZERO, the estimated value of a zero of
01103 //    the function F.
01104 //
01105 {
01106   double c;
01107   double d;
01108   double e;
01109   double fa;
01110   double fb;
01111   double fc;
01112   double m;
01113   double macheps;
01114   double p;
01115   double q;
01116   double r;
01117   double s;
01118   double sa;
01119   double sb;
01120   double tol;
01121 //
01122 //  Make local copies of A and B.
01123 //
01124   sa = a;
01125   sb = b;
01126   fa = f ( sa );
01127   fb = f ( sb );
01128
01129   c = sa;
01130   fc = fa;
01131   e = sb - sa;
01132   d = e;
01133
01134   macheps = r8_epsilon ( );
01135
01136   for ( ; ; )
01137   {
01138     if ( r8_abs ( fc ) < r8_abs ( fb ) )
01139     {
01140       sa = sb;
01141       sb = c;
01142       c = sa;
01143       fa = fb;
01144       fb = fc;
01145       fc = fa;
01146     }
01147
01148     tol = 2.0 * macheps * r8_abs ( sb ) + t;
01149     m = 0.5 * ( c - sb );
01150
01151     if ( r8_abs ( m ) <= tol || fb == 0.0 )
01152     {
01153       break;
01154     }
01155
01156     if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
01157     {
01158       e = m;
01159       d = e;
01160     }
01161     else
01162     {
01163       s = fb / fa;
01164
01165       if ( sa == c )
01166       {
01167         p = 2.0 * m * s;
01168         q = 1.0 - s;
01169       }
01170       else
01171       {
01172         q = fa / fc;
01173         r = fb / fc;
01174         p = s * ( 2.0 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
01175         q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
01176       }
01177
01178       if ( 0.0 < p )
01179       {
01180         q = - q;
01181       }
01182       else
01183       {
01184         p = - p;
01185       }
01186
01187       s = e;
01188       e = d;
01189
01190       if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
01191         p < r8_abs ( 0.5 * s * q ) )
01192       {
01193         d = p / q;
01194       }
01195       else
01196       {
01197         e = m;
01198         d = e;
01199       }
01200     }
01201     sa = sb;
01202     fa = fb;
01203
01204     if ( tol < r8_abs ( d ) )
01205     {
01206       sb = sb + d;
01207     }
01208     else if ( 0.0 < m )
01209     {
01210       sb = sb + tol;
01211     }
01212     else
01213     {
01214       sb = sb - tol;
01215     }
01216
01217     fb = f ( sb );
01218
01219     if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
01220     {
01221       c = sa;
01222       fc = fa;
01223       e = sb - sa;
01224       d = e;
01225     }
01226   }
01227   return sb;
01228 }
01229 //****************************************************************************80
01230
01231 void zero_rc ( double a, double b, double t, double &arg, int &status,
01232   double value )
01233
01234 //****************************************************************************80
01235 //
01236 //  Purpose:
01237 //
01238 //    ZERO_RC seeks the root of a function F(X) using reverse communication.
01239 //
01240 //  Discussion:
01241 //
01242 //    The interval [A,B] must be a change of sign interval for F.
01243 //    That is, F(A) and F(B) must be of opposite signs.  Then
01244 //    assuming that F is continuous implies the existence of at least
01245 //    one value C between A and B for which F(C) = 0.
01246 //
01247 //    The location of the zero is determined to within an accuracy
01248 //    of 6 * MACHEPS * r8_abs ( C ) + 2 * T.
01249 //
01250 //    The routine is a revised version of the Brent zero finder
01251 //    algorithm, using reverse communication.
01252 //
01253 //  Licensing:
01254 //
01256 //
01257 //  Modified:
01258 //
01259 //    17 July 2011
01260 //
01261 //  Author:
01262 //
01263 //    John Burkardt
01264 //
01265 //  Reference:
01266 //
01267 //    Richard Brent,
01268 //    Algorithms for Minimization Without Derivatives,
01269 //    Dover, 2002,
01270 //    ISBN: 0-486-41998-3,
01271 //    LC: QA402.5.B74.
01272 //
01273 //  Parameters:
01274 //
01275 //    Input, double A, B, the endpoints of the change of sign interval.
01276 //
01277 //    Input, double T, a positive error tolerance.
01278 //
01279 //    Output, double &ARG, the currently considered point.  The user
01280 //    does not need to initialize this value.  On return with STATUS positive,
01281 //    the user is requested to evaluate the function at ARG, and return
01282 //    the value in VALUE.  On return with STATUS zero, ARG is the routine's
01283 //    estimate for the function's zero.
01284 //
01285 //    Input/output, int &STATUS, used to communicate between
01286 //    the user and the routine.  The user only sets STATUS to zero on the first
01287 //    call, to indicate that this is a startup call.  The routine returns STATUS
01288 //    positive to request that the function be evaluated at ARG, or returns
01289 //    STATUS as 0, to indicate that the iteration is complete and that
01290 //    ARG is the estimated zero
01291 //
01292 //    Input, double VALUE, the function value at ARG, as requested
01293 //    by the routine on the previous call.
01294 //
01295 {
01296   static double c;
01297   static double d;
01298   static double e;
01299   static double fa;
01300   static double fb;
01301   static double fc;
01302   double m;
01303   static double macheps;
01304   double p;
01305   double q;
01306   double r;
01307   double s;
01308   static double sa;
01309   static double sb;
01310   double tol;
01311 //
01312 //  Input STATUS = 0.
01313 //  Initialize, request F(A).
01314 //
01315   if ( status == 0 )
01316   {
01317     macheps = r8_epsilon ( );
01318
01319     sa = a;
01320     sb = b;
01321     e = sb - sa;
01322     d = e;
01323
01324     status = 1;
01325     arg = a;
01326     return;
01327   }
01328 //
01329 //  Input STATUS = 1.
01330 //  Receive F(A), request F(B).
01331 //
01332   else if ( status == 1 )
01333   {
01334     fa = value;
01335     status = 2;
01336     arg = sb;
01337     return;
01338   }
01339 //
01340 //  Input STATUS = 2
01342 //
01343   else if ( status == 2 )
01344   {
01345     fb = value;
01346
01347     if ( 0.0 < fa * fb )
01348     {
01349       status = -1;
01350       return;
01351     }
01352     c = sa;
01353     fc = fa;
01354   }
01355   else
01356   {
01357     fb = value;
01358
01359     if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
01360     {
01361       c = sa;
01362       fc = fa;
01363       e = sb - sa;
01364       d = e;
01365     }
01366   }
01367 //
01368 //  Compute the next point at which a function value is requested.
01369 //
01370   if ( r8_abs ( fc ) < r8_abs ( fb ) )
01371   {
01372     sa = sb;
01373     sb = c;
01374     c = sa;
01375     fa = fb;
01376     fb = fc;
01377     fc = fa;
01378   }
01379
01380   tol = 2.0 * macheps * r8_abs ( sb ) + t;
01381   m = 0.5 * ( c - sb );
01382
01383   if ( r8_abs ( m ) <= tol || fb == 0.0 )
01384   {
01385     status = 0;
01386     arg = sb;
01387     return;
01388   }
01389
01390   if ( r8_abs ( e ) < tol || r8_abs ( fa ) <= r8_abs ( fb ) )
01391   {
01392     e = m;
01393     d = e;
01394   }
01395   else
01396   {
01397     s = fb / fa;
01398
01399     if ( sa == c )
01400     {
01401       p = 2.0 * m * s;
01402       q = 1.0 - s;
01403     }
01404     else
01405     {
01406       q = fa / fc;
01407       r = fb / fc;
01408       p = s * ( 2.0 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
01409       q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
01410     }
01411
01412     if ( 0.0 < p )
01413     {
01414       q = - q;
01415     }
01416     else
01417     {
01418       p = - p;
01419     }
01420     s = e;
01421     e = d;
01422
01423     if ( 2.0 * p < 3.0 * m * q - r8_abs ( tol * q ) &&
01424          p < r8_abs ( 0.5 * s * q ) )
01425     {
01426       d = p / q;
01427     }
01428     else
01429     {
01430       e = m;
01431       d = e;
01432     }
01433   }
01434
01435   sa = sb;
01436   fa = fb;
01437
01438   if ( tol < r8_abs ( d ) )
01439   {
01440     sb = sb + d;
01441   }
01442   else if ( 0.0 < m )
01443   {
01444     sb = sb + tol;
01445   }
01446   else
01447   {
01448     sb = sb - tol;
01449   }
01450
01451   arg = sb;
01452   status = status + 1;
01453
01454   return;
01455 }
01456
01457 // ======================================================================
01458 // === Simple wrapper functions
01459 // === for convenience and/or compatibility.
01460 //
01461 // === The three functions are the same as above,
01462 // === except that they take a plain function F
01463 // === instead of a c++ functor.  In all cases, the
01464 // === input and output of F() are of type double.
01465
01466 typedef double DoubleOfDouble (double);
01467
01468 class func_wrapper : public func_base {
01469   DoubleOfDouble* func;
01470 public:
01471   func_wrapper(DoubleOfDouble* f) {
01472     func = f;
01473   }
01474   virtual double operator() (double x){
01475     return func(x);
01476   }
01477 };
01478
01479 //****************************************************************************80
01480
01481 double glomin ( double a, double b, double c, double m, double e,
01482          double t, double f ( double x ), double &x ){
01483   func_wrapper foo(f);
01484   return glomin(a, b, c, m, e, t, foo, x);
01485 }
01486
01487 //****************************************************************************80
01488
01489 double local_min ( double a, double b, double t, double f ( double x ),
01490   double &x ){
01491   func_wrapper foo(f);
01492   return local_min(a, b, t, foo, x);
01493 }
01494
01495 //****************************************************************************80
01496
01497 double zero ( double a, double b, double t, double f ( double x ) ){
01498   func_wrapper foo(f);
01499   return zero(a, b, t, foo);
01500 }
01501
01502 // ======================================================================
01503 // Generally useful functor to evaluate a monic polynomial.
01504 // For details, see class definition in brent.hpp
01505
01506 double monicPoly::operator()(double x){
01507   double rslt(1);
01508   for (int ii = coeff.size()-1; ii >= 0; ii--){
01509     rslt *= x;
01510     rslt += coeff[ii];
01511   }
01512   return rslt;
01513 }
01514
01515 // Similarly, evaluate a general polynomial (not necessarily monic):
01516 double Poly::operator()(double x){
01517   double rslt(0);
01518   for (int ii = coeff.size()-1; ii >= 0; ii--){
01519     rslt *= x;
01520     rslt += coeff[ii];
01521   }
01522   return rslt;
01523 }
01524
01525 }
```

SHOGUN Machine Learning Toolbox - Documentation