SHOGUN
v2.0.0
|
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 // 00031 // This code is distributed under the GNU LGPL license. 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 00314 // about 1.324.... 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 // 00329 // This code is distributed under the GNU LGPL license. 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 // 00572 // This code is distributed under the GNU LGPL license. 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 // 00846 // This code is distributed under the GNU LGPL license. 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 // 00895 // This code is distributed under the GNU LGPL license. 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 // 00933 // This code is distributed under the GNU LGPL license. 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 // 00974 // This code is distributed under the GNU LGPL license. 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 // 01019 // This code is distributed under the GNU LGPL license. 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 // 01072 // This code is distributed under the GNU LGPL license. 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 // 01255 // This code is distributed under the GNU LGPL license. 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 01341 // Receive F(B). 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 }