//=========================================================================== // GoTools Core - SINTEF Geometry Tools Core library, version 2.0.1 // // Copyright (C) 2000-2007, 2010 SINTEF ICT, Applied Mathematics, Norway. // // 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 version 2 of the License. // // 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., // 59 Temple Place - Suite 330, // Boston, MA 02111-1307, USA. // // Contact information: E-mail: tor.dokken@sintef.no // SINTEF ICT, Department of Applied Mathematics, // P.O. Box 124 Blindern, // 0314 Oslo, Norway. // // Other licenses are also available for this software, notably licenses // for: // - Building commercial software. // - Building software whose source code you wish to keep private. //=========================================================================== 00015 #ifndef _GENERAL_FUNCTIONMINIMIZER_IMPLEMENTATION_H 00016 #define _GENERAL_FUNCTIONMINIMIZER_IMPLEMENTATION_H 00017 00018 #include <limits> 00019 #include <assert.h> 00020 #include "GoTools/utils/errormacros.h" 00021 #include <algorithm> 00022 #include <iostream> // @@debug 00023 00024 namespace Go { 00025 00026 //=========================================================================== 00027 template<class Functor> 00028 const double FunctionMinimizer<Functor>::root_machine_precision_ = 00029 sqrt(numerical_tolerance(1)); 00030 00031 template<class Functor> 00032 const double FunctionMinimizer<Functor>::rel_tol_ = 1.0e-6; 00033 00034 template<class Functor> 00035 const double FunctionMinimizer<Functor>::perturbation_ = 1.0e-10; 00036 00037 template<class Functor> 00038 const double FunctionMinimizer<Functor>::golden_ratio_ = 0.61803399; 00039 00040 template<class Functor> 00041 const double FunctionMinimizer<Functor>::default_partition_ = 00042 1.0 / 20.0; // default partition ratio of an interval 00043 00044 template<class Functor> 00045 const double FunctionMinimizer<Functor>::tiny_ = 1.0e-20; // a really tiny number 00046 00047 template<class Functor> 00048 const int FunctionMinimizer<Functor>::max_iter_ = 100; 00049 //=========================================================================== 00050 00051 //=========================================================================== 00052 template<class Functor> 00053 void minimise_conjugated_gradient(FunctionMinimizer<Functor>& dfmin) 00054 //=========================================================================== 00055 { 00056 const double TOL = std::numeric_limits<double>::epsilon(); //1.0e-8; 00057 const double EPS = 1.0e-10; 00058 // minimising the 'dfmin' function using conjugated gradients. 00059 const int N = dfmin.numPars(); 00060 Point gradient(N), old_gradient(N), dir(N); 00061 dfmin.grad(old_gradient); 00062 dir = -old_gradient; 00063 double old_val = dfmin.fval(); 00064 while(true) { 00065 00066 // make sure direction is not uphill (is this already guaranteed??) 00067 // and truncating it if we are at the border of the domain 00068 if (dir * old_gradient > 0) { 00069 dir *= -1; 00070 } 00071 for (int i = 0; i < N; ++i) { 00072 if ((dfmin.atMin(i) && dir[i] < 0) || (dfmin.atMax(i) && dir[i] > 0)) { 00073 dir[i] = 0; 00074 // conjugate gradient value is broken since we modified dir. 00075 // Restarting cycle. 00076 } 00077 } 00078 00079 // minimize along this direction 00080 bool hit_domain_edge = false; 00081 double new_val = dfmin.minimize(dir, hit_domain_edge); 00082 if (2.0 * fabs(new_val - old_val) <= TOL * (fabs(new_val) + fabs(old_val)+ EPS)) { 00083 // we have reached a minimum 00084 break; 00085 } else { 00086 old_val = new_val; 00087 } 00088 00089 // choose new direction 00090 dfmin.grad(gradient); 00091 if (!hit_domain_edge) { 00092 // we are still in a conjugated gradient cycle. Choose new direction 00093 // using conjugated gradients (Polak-Ribiere variant) 00094 00095 Point diff = gradient - old_gradient; 00096 double factor = gradient * diff / old_gradient.length2(); 00097 Point dir_saved = dir; 00098 dir *= factor; 00099 dir -= gradient; 00100 00101 if (dir * old_gradient > 0) { 00102 dir *= -1; 00103 } 00104 00105 bool on_boundary = false; 00106 for (int i = 0; i != N; ++i) { 00107 if ((dfmin.atMin(i) && dir[i] < 0) || (dfmin.atMax(i) && dir[i] > 0)) { 00108 on_boundary = true; 00109 gradient[i] = old_gradient[i] = 0; 00110 dir_saved[i] = 0; // we will have to recalculated 'dir' based on this vector 00111 } 00112 } 00113 00114 if (on_boundary) { 00115 // the direction we choose will take us out of the domain. Reduce problem 00116 // to conjugated gradient in a lower dimension. 00117 dir = dir_saved; 00118 diff = gradient - old_gradient; 00119 factor = gradient * diff / old_gradient.length2(); 00120 dir *= factor; 00121 dir -= gradient; 00122 } 00123 old_gradient = gradient; 00124 } else { 00125 // We ran into an edge of the domain. Re-initialising conjugate 00126 // gradient method using steepest descent. 00127 dir = - gradient; 00128 for (int i = 0; i != N; ++i) { 00129 if ((dfmin.atMin(i) && dir[i] < 0) || (dfmin.atMax(i) && dir[i] > 0)) { 00130 dir[i] = 0; 00131 gradient[i] = 0; 00132 } 00133 } 00134 old_gradient = gradient; 00135 } 00136 } 00137 } 00138 00139 00140 //=========================================================================== 00141 template<class Functor> 00142 inline FunctionMinimizer<Functor>:: 00143 FunctionMinimizer(int num_param, const Functor& fun, const double* const seed, double tol) 00144 : fun_(fun), par_(seed, seed + num_param), 00145 at_min_(num_param), at_max_(num_param), cached_value_updated_(false), 00146 cached_grad_updated_(false), cached_grad_(num_param), 00147 minimization_tol_(tol > root_machine_precision_ ? tol : root_machine_precision_), 00148 param_tol_(num_param) 00149 //=========================================================================== 00150 { 00151 for (int i = 0; i < num_param; ++i) { 00152 param_tol_[i] = rel_tol_ * (fun_.maxPar(i) - fun_.minPar(i)); 00153 } 00154 00155 checkBorder(); // determine at_min_ and at_max_ 00156 resetCache(); 00157 } 00158 00159 //=========================================================================== 00160 template<class Functor> 00161 inline void FunctionMinimizer<Functor>::resetCache() const 00162 //=========================================================================== 00163 { 00164 cached_value_updated_ = false; 00165 cached_grad_updated_ = false; 00166 } 00167 00168 //=========================================================================== 00169 template<class Functor> 00170 inline void FunctionMinimizer<Functor>::checkBorder() 00171 //=========================================================================== 00172 { 00173 int n = numPars(); 00174 for (int i = 0; i < n; ++i) { 00175 at_min_[i] = par_[i] < (fun_.minPar(i) + param_tol_[i]); 00176 at_max_[i] = par_[i] > (fun_.maxPar(i) - param_tol_[i]); 00177 } 00178 } 00179 00180 //=========================================================================== 00181 template<class Functor> 00182 inline double FunctionMinimizer<Functor>::minimize(const Point& dir, 00183 bool& hit_domain_edge) 00184 //=========================================================================== 00185 { 00186 static bool rerun = false; // 'true' if we have recursively entered this function 00187 // flipping dir if it conflicts with gradient of function (will this ever happen??) 00188 Point oriented_dir(numPars()); 00189 if (!orientDirection(dir, oriented_dir)) { 00190 // gradient of function was exactly perpendicular to 'dir'. We are already at 00191 // a minimum along this direction 00192 return fval(); 00193 } 00194 // deciding maximum possible steplength in the given direction 00195 double max_step = determineMaxSteplength(oriented_dir); 00196 if (max_step < tiny_) { 00197 // already on boundary in this direction. 00198 return fval(); 00199 } 00200 // bracketing minimum in the given direction 00201 double bracket[3]; // startpoint, midpoint and endpoint of bracketing interval 00202 double fval_brak[3]; // function value at the three bracketing points 00203 00204 if (bracketInterval(oriented_dir, max_step, bracket, fval_brak)) { 00205 // bracketing of interval successful. Applying Brent 00206 hit_domain_edge = false; 00207 return linminBrent(oriented_dir, bracket, fval_brak); 00208 } 00209 // We were not able to bracket the interval. This may either be because the search 00210 // interval grew all the way to 'max_step' without being able to bracket a minimum, 00211 // or that the search interval shrank and collapsed in a singularity, indicating 00212 // that although we know this to be a direction of descent, numerically the function 00213 // is not able to produce a smaller value in the neighbourhood of the departing point. 00214 if (fval_brak[1] >= fval_brak[0]) { 00215 ASSERT(bracket[2] < max_step); 00216 // the interval has collapsed onto the starting point. We consider this to be 00217 // a minimum (its directional derivative, although not zero, is so small that 00218 // it does not modify the function within machine precision). 00219 return fval(); 00220 } 00221 ASSERT(bracket[2] == max_step); 00222 // If we got here, we have maximally expanded the interval without being able to 00223 // verify that we have bracketed a minimum. The conclusion is that the minimum 00224 // is located very close to max_step. Our strategy now is to find outwhat happens 00225 // at the end of the interval. If the gradient is positive, we know that we have 00226 // just passed the minimum, and we can search backwards from the endpoint. If not, 00227 // we conclude that this is a minimum. We check if this point gives a smaller 00228 // function value than where we stand now, and if that is the case, we move there 00229 // before returning. 00230 Point end_grad(numPars()); 00231 Point temp= par_ + max_step * oriented_dir; 00232 grad(temp, end_grad); 00233 double max_step_val = fval(temp); 00234 00235 if (scalarProductSign(end_grad, oriented_dir) == 1) { 00236 // derivative is POSITIVE at 'max_step'. There MUST be a minimum that we 00237 // missed. Let us search backwards for it. 00238 if (!rerun) { 00239 moveUV(oriented_dir, max_step); // moving to end of interval 00240 rerun = true; 00241 double tmp = minimize(dir, hit_domain_edge); // search backwds 00242 rerun = false; 00243 return tmp; 00244 } 00245 MESSAGE("Warning: unable to pinpoint exact minimum in FunctionMinimizer::minimize()"); 00246 } 00247 // the derivative is NEGATIVE at 'max_step', and we conclude that this is a minimum 00248 // as good as any. If it gives a smaller function value than where we stand now, 00249 // we will move there. Otherwise, we stay put and consider ourselves already on 00250 // a minimum. 00251 if (max_step_val < fval()) { 00252 hit_domain_edge = true; 00253 moveUV(oriented_dir, max_step); 00254 return max_step_val; 00255 } 00256 // unable to find any point better than where we are currently standing 00257 hit_domain_edge = false; 00258 return fval(); 00259 } 00260 00261 //=========================================================================== 00262 template<class Functor> inline 00263 int FunctionMinimizer<Functor>::scalarProductSign(const Point& p1, const Point& p2) 00264 //=========================================================================== 00265 { 00266 double prod = p1 * p2; 00267 double l1 = *std::max_element(p1.begin(), p1.end()); 00268 double l2 = *std::max_element(p2.begin(), p2.end()); 00269 // find the numerical tolerance related to the point coordinate of the largest 00270 // magnitude (its error dwarfs the errors of the others). 00271 double tol = 10 * numerical_tolerance(l1 > l2 ? l1 : l2); 00272 00273 if (prod < -tol) return -1; 00274 if (prod > tol ) return 1; 00275 00276 // indiscernable from 0 00277 return 0; 00278 } 00279 00280 //=========================================================================== 00281 template<class Functor> inline bool 00282 FunctionMinimizer<Functor>::orientDirection(const Point& dir, Point& result) 00283 //=========================================================================== 00284 { 00285 Point start_grad(numPars()); 00286 grad(start_grad); 00287 00288 result = dir / dir.lengthInf(); 00289 00290 int sign = scalarProductSign(start_grad, result); 00291 00292 if (sign == 1) { 00293 result *= -1; 00294 } 00295 return (sign != 0); 00296 } 00297 00298 //=========================================================================== 00299 template<class Functor> 00300 inline double FunctionMinimizer<Functor>::linminBrent(const Point& dir, 00301 const double* bracket, 00302 const double* fval_brak) 00303 //=========================================================================== 00304 { 00305 const int MAXITER = 100; 00306 double a = bracket[0]; double fa = fval_brak[0]; 00307 double b = bracket[1]; double fb = fval_brak[1]; 00308 double c = bracket[2]; double fc = fval_brak[2]; 00309 00310 if (!(a <= b && b <= c)) { 00311 THROW("Error: Assertion 'a <= b && b <= c' failed."); 00312 } 00313 if (!(fa >= fb && fb <= fc)) { 00314 THROW("Error: Assertion 'fa >= fb && fb <= fc' failed."); 00315 } 00316 00317 double b2, fb2, b3, fb3; 00318 if (fa < fc) { 00319 b2 = a; fb2 = fa; 00320 b3 = c; fb3 = fc; 00321 } else { 00322 b2 = c; fb2 = fc; 00323 b3 = a; fb3 = fa; 00324 } 00325 00326 double u = b; 00327 double last_step = std::numeric_limits<double>::max(); 00328 double before_last_step = last_step; 00329 double fu; 00330 int iter; 00331 // a step of 'delta' moves u and v by respectively delta * dir[0] and delta * dir[1]. 00332 // therefore, if we want a precision of (minimization_tol_ * central bracket value) 00333 // in each parameter, we must use the tolerance below: 00334 const double central_value = fabs(b) > perturbation_ ? fabs(b) : perturbation_; 00335 double tol_1 = minimization_tol_ * central_value; 00336 // Since a, b, c and u are "fictive" function arguments (the real argument in R^n being 00337 // 'par_ + u * dir'), we must check if the minimum tolerance of this last expression 00338 // exceeds the tolerance we have chosen. In that case, increase the tolerance to this 00339 // new value. 00340 double tol_2 = numerical_tolerance(par_, dir); 00341 const double tol = tol_1 > tol_2 ? tol_1 : tol_2; 00342 00343 // three below variables to speed up convergence when suspecting to be really close 00344 // to the function minimum 00345 bool last_step_was_tol = false; 00346 bool last_step_was_positive = false; 00347 double btemp = b2; 00348 for (iter = 0; iter < MAXITER; ++iter) { 00349 // ------------------------------------- 00350 // STEP 1: coming up with a point to try 00351 // ------------------------------------- 00352 00353 // trying parabolic interpolation 00354 bool use_parabolic = parabolicFit(b, fb, b2, fb2, b3, fb3, 0.5 * before_last_step, u); 00355 before_last_step = last_step; 00356 00357 if (use_parabolic) { 00358 // so far, it seems like we can use a parabolic fit. But is it inside the interval? 00359 use_parabolic = (u > a + tol && u < c - tol); // useful only if inside the interval 00360 } 00361 // here, we should know whether the parabolic fit is useful or not 00362 if (!use_parabolic) { 00363 // we must come up with another point to try. Resorting to golden search 00364 const double rdiff = c - b; 00365 const double ldiff = b - a; 00366 if (rdiff < ldiff) { 00367 u = a + golden_ratio_ * ldiff; 00368 } else { 00369 u = c - golden_ratio_ * rdiff; 00370 } 00371 } 00372 // assuring a minimum steplength 00373 if (fabs(u - b) < tol) { 00374 if (last_step_was_tol && btemp == b) { // exact equality is justified here! 00375 u = last_step_was_positive ? b - tol : b + tol; 00376 } else { 00377 u = (u > b) ? b + tol : b - tol; 00378 } 00379 btemp = b; 00380 last_step_was_positive = (u > b); 00381 last_step_was_tol = true; 00382 00383 } else { 00384 last_step_was_tol = false; 00385 } 00386 last_step = fabs(b - u); 00387 00388 // -------------------------------------------------------- 00389 // STEP 2: trying the point and deciding what to do with it 00390 // -------------------------------------------------------- 00391 00392 // evaluating function in u (and evaluating gradient if premature end is allowed) 00393 00394 fu = fval(par_ + u * dir); 00395 adjustBrackets(a, fa, b, fb, c, fc, b2, fb2, b3, fb3, u, fu); 00396 00397 // end criteria 00398 if (fabs(b - (0.5 * (c + a))) < 2 * tol - 0.5 * (c - a)) { 00399 break; 00400 } 00401 } 00402 if (iter == MAXITER) { 00403 MESSAGE("Failed to converge properly in linminBrent."); 00404 //throw runtime_error("Failed to converge in linminBrent."); 00405 } 00406 00407 // setting current point at the found minimum 00408 moveUV(dir, b); 00409 00410 return fb; 00411 } 00412 00413 //=========================================================================== 00414 template<class Functor> 00415 inline void FunctionMinimizer<Functor>::adjustBrackets(double& a, double& fa, 00416 double& b, double& fb, 00417 double& c, double& fc, 00418 double& b2, double& fb2, 00419 double& b3, double& fb3, 00420 double u, double fu) 00421 //=========================================================================== 00422 { 00423 if (fu < fb) { 00424 // narrowing down brackets 00425 if (u < b) { 00426 c = b; fc = fb; 00427 } else { 00428 a = b; fa = fb; 00429 } 00430 // updating minima 00431 b3 = b2; fb3 = fb2; 00432 b2 = b; fb2 = fb; 00433 b = u; fb = fu; 00434 } else { 00435 // we did not find a better minima, but we can still narrow down brackets 00436 if (u < b) { 00437 a = u; fa = fu; 00438 } else { 00439 c = u; fc = fu; 00440 } 00441 00442 // checking if previous minima should be updated 00443 if (fu < fb2) { 00444 b3 = b2; fb3 = fb2; 00445 b2 = u; fb2 = fu; 00446 } else if (fu < fb3) { 00447 b3 = u; fb3 = fu; 00448 } 00449 } 00450 } 00451 00452 //=========================================================================== 00453 template<class Functor> inline double 00454 FunctionMinimizer<Functor>::determineMaxSteplength(const Point& dir) const 00455 //=========================================================================== 00456 { 00457 const int n = numPars(); 00458 std::vector<double> max_uv(n); 00459 00460 for (int i = 0; i < n; ++i) { 00461 if (dir[i] > param_tol_[i]) { 00462 double delta = fun_.maxPar(i) - par_[i]; 00463 max_uv[i] = delta / dir[i]; // a positive value 00464 } else if (dir[i] < -param_tol_[i]) { 00465 double delta = fun_.minPar(i) - par_[i]; 00466 max_uv[i] = delta / dir[i]; // a positive value 00467 } else { 00468 max_uv[i] = std::numeric_limits<double>::max(); // largest machine-repr. number 00469 } 00470 } 00471 return *(min_element(max_uv.begin(), max_uv.end())); 00472 } 00473 00474 //=========================================================================== 00475 template<class Functor> inline double 00476 FunctionMinimizer<Functor>::parabolicEstimate(double p2, double fp2, const Point& dir) 00477 //=========================================================================== 00478 { 00479 Point gradient(numPars()); 00480 grad(gradient); 00481 double directional_deriv = gradient * dir; 00482 return -directional_deriv * p2 * p2 / (2 * (fp2 - fval() - p2 * directional_deriv)); 00483 } 00484 00485 //=========================================================================== 00486 template<class Functor> inline bool 00487 FunctionMinimizer<Functor>::bracketInterval(const Point& dir, // points to 2-array 00488 const double max_step, 00489 //const double oper_tol, // operating tolerance 00490 double* bracket, // points to 3-array 00491 double* fval_brak) // points to 3-array 00492 //=========================================================================== 00493 { 00494 // We suppose that 'dir' is a direciton of descent - anything else is an error! 00495 // we want to find three collinear points (a, b, c) on the line with direction 'dir' 00496 // and passing through the point (curU(), curV()) such that f(a) > f(b) < f(c). 00497 00498 if (!(max_step > 0)) { 00499 THROW("Error: Assertion 'max_step > 0' failed."); 00500 } 00501 //ASSERT(oper_tol > 0); 00502 //ASSERT(max_step > oper_tol); 00503 00504 double& a = bracket[0]; // the three points we want to determine 00505 double& b = bracket[1]; 00506 double& c = bracket[2]; 00507 double& fa = fval_brak[0]; 00508 double& fb = fval_brak[1]; 00509 double& fc = fval_brak[2]; // function value in these points 00510 00511 a = 0; 00512 fa = fval(); // current point (curU(), curV()) 00513 00514 b = default_partition_ * max_step; 00515 fb = fval(par_ + b * dir); 00516 00517 // determining tolerance based on numerical precision of machine 00518 Point g; grad(g); 00519 double num_tolerance = numerical_tolerance(par_, dir, fa, dir * g); 00520 00521 if (fb >= fa) { 00522 // we know there is a minimum between a and b, since 'dir' is assumed to be a 00523 // direction of descent. We now only have to bisect the interval sufficiently 00524 // many times towards 'a' in order to bracket the minimum 00525 //std::cout << "Fallen into slow bracketing." << std::endl; 00526 do { 00527 c = b; fc = fb; 00528 b = parabolicEstimate(c, fc, dir); 00529 fb = fval(par_ + b * dir); 00530 if (!(b > a && b < c)) { 00531 THROW("Error: Assertion 'b > a && b < c' failed."); 00532 } 00533 } while (fb >= fa && fabs(b-a) > num_tolerance); 00534 00535 if (fb >= fa) { 00536 // we did not succeed in bracketing a minimum between a and c, due to 00537 // numerical issues (although 'dir' is direction of descent, we were 00538 // unable to find a smaller function value within a distance from the 00539 // startpoint bigger than the numerical tolerance). 00540 return false; 00541 } 00542 if (!(fb < fa && fb <= fc)) { 00543 THROW("Error: Assertion 'fb < fa && fb <= fc' failed."); 00544 } 00545 return true; 00546 } 00547 00548 // if we got here, (fb < fa) -- we have (probably) not yet arrived at the minimum value 00549 // the below code is based on ideas from the book "Numerical Recipes in C" for bracketing 00550 // a minimum. 00551 c = b + (1 + golden_ratio_) * (b - a); 00552 c = (c < max_step) ? c : max_step; 00553 fc = fval(par_ + c * dir); 00554 00555 //std::cout << "Entering normal bracketing." << std::endl; 00556 // we here know that a < b < c <= max_step, and that fb < fa 00557 while (fb > fc) { 00558 // coming up with suggestion for new point, based on parabolic interp. of a, b and c. 00559 double fu, u; 00560 parabolicFit(a, fa, b, fb, c, fc, max_step - b, u); 00561 00562 double interval_length = c - a; 00563 if (fabs(u-b) > interval_length) { 00564 u = (u > b) ? b + interval_length : b - interval_length; 00565 } 00566 u = (u < max_step) ? u : max_step; // not really necessary, but safeguard against roundoff 00567 00568 // deciding what to do with new, estimated point 00569 if (b < u && u < c) { 00570 fu = fval(par_ + u * dir); 00571 if (fu < fc) { // minimum between b and c 00572 a = b; fa = fb; 00573 b = u; fb = fu; 00574 return true; 00575 } else if (fu > fb) { // minimum between a and u 00576 c = u; fc = fu; 00577 return true; 00578 } 00579 // no help in parabolic fit 00580 u = c + (1 + golden_ratio_) * (c - b); 00581 u = (u < max_step) ? u : max_step; 00582 } else if (c < u) { 00583 fu = fval(par_ + u * dir); 00584 if (fu < fc) { 00585 shift3(b, c, u, u + (1 + golden_ratio_) * (u - c)); 00586 u = (u < max_step) ? u : max_step; 00587 shift3(fb, fc, fu, fval(par_ + u * dir)); 00588 } 00589 } else { // reject parabolic u 00590 u = c + (1 + golden_ratio_) * (c - b); 00591 u = (u < max_step) ? u : max_step; 00592 fu = fval(par_ + u * dir); 00593 } 00594 shift3(a, b, c, u); 00595 shift3(fa, fb, fc, fu); 00596 if (c == max_step && fc <= fb) { // yes, equality check of doubles is justified here 00597 return false; 00598 } 00599 // if (fabs(c - max_step) < minimization_tol_ && fc <= fb) { 00600 // return false; 00601 // } 00602 } 00603 // if we got here, then fa > fb < fc 00604 if (!(fb < fa && fb <= fc)) { 00605 THROW("Error: Assertion 'fb < fa && fb <= fc' failed."); 00606 } 00607 return true; 00608 } 00609 00610 //=========================================================================== 00611 template<class Functor> 00612 inline bool FunctionMinimizer<Functor>::parabolicFit(double a, double fa, 00613 double b, double fb, 00614 double c, double fc, 00615 double max_from_b, double&u) 00616 //=========================================================================== 00617 { 00618 // returns 'true' if it could determine a minimum on the parabola containing 00619 // the points (a, fa), (b, fb) and (c, fc). If no such minimum could be found 00620 // (ie. the function is linear), then it returns false. 00621 if (!(max_from_b > 0)) { 00622 THROW("Error: Assertion 'max_from_b > 0' failed."); 00623 } 00624 00625 double tmp1 = (b - a) * (fb - fc); 00626 double tmp2 = (b - c) * (fb - fa); 00627 double denom = 2 * (tmp1 - tmp2); // denominator 00628 u = (b - a) * tmp1 - (b - c) * tmp2; // nominator 00629 00630 if (fabs(u) < fabs(denom) * max_from_b) { 00631 u /= denom; 00632 u = b - u; 00633 return true; 00634 } else { 00635 // we would have to step farther than allowed. Set to max and return false. 00636 u = (u * denom > 0) ? b - max_from_b : b + max_from_b; 00637 return false; 00638 } 00639 // never get here, but keep compiler happy 00640 return true; 00641 } 00642 00643 //=========================================================================== 00644 template<class Functor> 00645 inline double FunctionMinimizer<Functor>::numerical_tolerance(double c) 00646 //=========================================================================== 00647 { 00648 if (c == double(0)) { 00649 return tiny_; 00650 } 00651 static const double eps = std::numeric_limits<double>::epsilon(); 00652 return fabs(c) * eps; 00653 } 00654 00655 //=========================================================================== 00656 template<class Functor> 00657 inline double FunctionMinimizer<Functor>::numerical_tolerance(const Point& x, 00658 const Point& dir) 00659 //=========================================================================== 00660 { 00661 // we want to determine the smallest positive double value 'a' such that 00662 // (x + a * dir) != x. If we define tol(x_i) to be the numerical precision 00663 // around the i-component of x, and dir_i to be the i-component of y, 00664 // then we search for the value: 00665 // minimize_over_i (tol(x_i) / |dir_i|) 00666 00667 double res = std::numeric_limits<double>::max(); 00668 for (int i = 0; i < x.size(); ++i) { 00669 if (dir[i] != double(0.0)) { 00670 double temp = numerical_tolerance(x[i]) / fabs(dir[i]); 00671 res = res < temp ? res : temp; 00672 //double temp = fabs(x[i]) * std::numeric_limits<double>::epsilon(); //tol(x_i) 00673 //temp /= fabs(dir[i]); // |dir_i| 00674 //res = res < temp ? res : temp; 00675 } 00676 } 00677 return res; 00678 } 00679 00680 //=========================================================================== 00681 template<class Functor> 00682 inline double FunctionMinimizer<Functor>::numerical_tolerance(const Point& x, 00683 const Point& dir, 00684 const double fx, 00685 const double dfx) 00686 //=========================================================================== 00687 { 00688 // we want to determine the smallest positive double value 'a' such that 00689 // -> (x + a * dir) != x 00690 // and such that 00691 // -> f(x + a * dir) != f(x). 00692 // To satisfy the last criteria, we use the function value in x, given by 00693 // 'fx', as well as the directional derivative, which is given by 'dfx'. 00694 00695 00696 double a1 = numerical_tolerance(x, dir); 00697 00698 double a2 = (fabs(fx) + perturbation_) * std::numeric_limits<double>::epsilon(); 00699 a2 /= fabs(dfx); 00700 00701 return a1 > a2 ? a1 : a2; 00702 } 00703 00704 //=========================================================================== 00705 template<class Functor> 00706 inline double FunctionMinimizer<Functor>::fval() const 00707 //=========================================================================== 00708 { 00709 if (!cached_value_updated_) { 00710 cached_value_ = fun_(par_.begin()); 00711 cached_value_updated_ = true; 00712 } 00713 return cached_value_; 00714 } 00715 00716 //=========================================================================== 00717 template<class Functor> 00718 inline double FunctionMinimizer<Functor>::fval(const Point& par) const 00719 //=========================================================================== 00720 { 00721 return fun_(par.begin()); 00722 } 00723 00724 //=========================================================================== 00725 template<class Functor> 00726 inline void FunctionMinimizer<Functor>::grad(Point& result) const 00727 //=========================================================================== 00728 { 00729 if (!cached_grad_updated_) { 00730 fun_.grad(par_.begin(), cached_grad_.begin()); 00731 cached_grad_updated_ = true; 00732 } 00733 result = cached_grad_; 00734 } 00735 00736 //=========================================================================== 00737 template<class Functor> 00738 inline void FunctionMinimizer<Functor>::grad(const Point& param, Point& result) const 00739 //=========================================================================== 00740 { 00741 fun_.grad(param.begin(), result.begin()); 00742 } 00743 00744 //=========================================================================== 00745 template<class Functor> 00746 inline void FunctionMinimizer<Functor>::moveUV(const Point& dir, double multiplier) 00747 //=========================================================================== 00748 { 00749 par_ += dir * multiplier; 00750 checkBorder(); 00751 resetCache(); 00752 } 00753 00754 00755 }; // end namespace Go 00756 00757 #endif // _GENERAL_FUNCTIONMINIMIZER_IMPLEMENTATION_H 00758