![]() |
Eigen-unsupported
3.3.3
|
00001 // -*- coding: utf-8 00002 // vim: set fileencoding=utf-8 00003 00004 // This file is part of Eigen, a lightweight C++ template library 00005 // for linear algebra. 00006 // 00007 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> 00008 // 00009 // This Source Code Form is subject to the terms of the Mozilla 00010 // Public License v. 2.0. If a copy of the MPL was not distributed 00011 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00012 00013 #ifndef EIGEN_LEVENBERGMARQUARDT__H 00014 #define EIGEN_LEVENBERGMARQUARDT__H 00015 00016 namespace Eigen { 00017 00018 namespace LevenbergMarquardtSpace { 00019 enum Status { 00020 NotStarted = -2, 00021 Running = -1, 00022 ImproperInputParameters = 0, 00023 RelativeReductionTooSmall = 1, 00024 RelativeErrorTooSmall = 2, 00025 RelativeErrorAndReductionTooSmall = 3, 00026 CosinusTooSmall = 4, 00027 TooManyFunctionEvaluation = 5, 00028 FtolTooSmall = 6, 00029 XtolTooSmall = 7, 00030 GtolTooSmall = 8, 00031 UserAsked = 9 00032 }; 00033 } 00034 00035 00036 00045 template<typename FunctorType, typename Scalar=double> 00046 class LevenbergMarquardt 00047 { 00048 static Scalar sqrt_epsilon() 00049 { 00050 using std::sqrt; 00051 return sqrt(NumTraits<Scalar>::epsilon()); 00052 } 00053 00054 public: 00055 LevenbergMarquardt(FunctorType &_functor) 00056 : functor(_functor) { nfev = njev = iter = 0; fnorm = gnorm = 0.; useExternalScaling=false; } 00057 00058 typedef DenseIndex Index; 00059 00060 struct Parameters { 00061 Parameters() 00062 : factor(Scalar(100.)) 00063 , maxfev(400) 00064 , ftol(sqrt_epsilon()) 00065 , xtol(sqrt_epsilon()) 00066 , gtol(Scalar(0.)) 00067 , epsfcn(Scalar(0.)) {} 00068 Scalar factor; 00069 Index maxfev; // maximum number of function evaluation 00070 Scalar ftol; 00071 Scalar xtol; 00072 Scalar gtol; 00073 Scalar epsfcn; 00074 }; 00075 00076 typedef Matrix< Scalar, Dynamic, 1 > FVectorType; 00077 typedef Matrix< Scalar, Dynamic, Dynamic > JacobianType; 00078 00079 LevenbergMarquardtSpace::Status lmder1( 00080 FVectorType &x, 00081 const Scalar tol = sqrt_epsilon() 00082 ); 00083 00084 LevenbergMarquardtSpace::Status minimize(FVectorType &x); 00085 LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x); 00086 LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x); 00087 00088 static LevenbergMarquardtSpace::Status lmdif1( 00089 FunctorType &functor, 00090 FVectorType &x, 00091 Index *nfev, 00092 const Scalar tol = sqrt_epsilon() 00093 ); 00094 00095 LevenbergMarquardtSpace::Status lmstr1( 00096 FVectorType &x, 00097 const Scalar tol = sqrt_epsilon() 00098 ); 00099 00100 LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x); 00101 LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType &x); 00102 LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType &x); 00103 00104 void resetParameters(void) { parameters = Parameters(); } 00105 00106 Parameters parameters; 00107 FVectorType fvec, qtf, diag; 00108 JacobianType fjac; 00109 PermutationMatrix<Dynamic,Dynamic> permutation; 00110 Index nfev; 00111 Index njev; 00112 Index iter; 00113 Scalar fnorm, gnorm; 00114 bool useExternalScaling; 00115 00116 Scalar lm_param(void) { return par; } 00117 private: 00118 00119 FunctorType &functor; 00120 Index n; 00121 Index m; 00122 FVectorType wa1, wa2, wa3, wa4; 00123 00124 Scalar par, sum; 00125 Scalar temp, temp1, temp2; 00126 Scalar delta; 00127 Scalar ratio; 00128 Scalar pnorm, xnorm, fnorm1, actred, dirder, prered; 00129 00130 LevenbergMarquardt& operator=(const LevenbergMarquardt&); 00131 }; 00132 00133 template<typename FunctorType, typename Scalar> 00134 LevenbergMarquardtSpace::Status 00135 LevenbergMarquardt<FunctorType,Scalar>::lmder1( 00136 FVectorType &x, 00137 const Scalar tol 00138 ) 00139 { 00140 n = x.size(); 00141 m = functor.values(); 00142 00143 /* check the input parameters for errors. */ 00144 if (n <= 0 || m < n || tol < 0.) 00145 return LevenbergMarquardtSpace::ImproperInputParameters; 00146 00147 resetParameters(); 00148 parameters.ftol = tol; 00149 parameters.xtol = tol; 00150 parameters.maxfev = 100*(n+1); 00151 00152 return minimize(x); 00153 } 00154 00155 00156 template<typename FunctorType, typename Scalar> 00157 LevenbergMarquardtSpace::Status 00158 LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType &x) 00159 { 00160 LevenbergMarquardtSpace::Status status = minimizeInit(x); 00161 if (status==LevenbergMarquardtSpace::ImproperInputParameters) 00162 return status; 00163 do { 00164 status = minimizeOneStep(x); 00165 } while (status==LevenbergMarquardtSpace::Running); 00166 return status; 00167 } 00168 00169 template<typename FunctorType, typename Scalar> 00170 LevenbergMarquardtSpace::Status 00171 LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType &x) 00172 { 00173 n = x.size(); 00174 m = functor.values(); 00175 00176 wa1.resize(n); wa2.resize(n); wa3.resize(n); 00177 wa4.resize(m); 00178 fvec.resize(m); 00179 fjac.resize(m, n); 00180 if (!useExternalScaling) 00181 diag.resize(n); 00182 eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'"); 00183 qtf.resize(n); 00184 00185 /* Function Body */ 00186 nfev = 0; 00187 njev = 0; 00188 00189 /* check the input parameters for errors. */ 00190 if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.) 00191 return LevenbergMarquardtSpace::ImproperInputParameters; 00192 00193 if (useExternalScaling) 00194 for (Index j = 0; j < n; ++j) 00195 if (diag[j] <= 0.) 00196 return LevenbergMarquardtSpace::ImproperInputParameters; 00197 00198 /* evaluate the function at the starting point */ 00199 /* and calculate its norm. */ 00200 nfev = 1; 00201 if ( functor(x, fvec) < 0) 00202 return LevenbergMarquardtSpace::UserAsked; 00203 fnorm = fvec.stableNorm(); 00204 00205 /* initialize levenberg-marquardt parameter and iteration counter. */ 00206 par = 0.; 00207 iter = 1; 00208 00209 return LevenbergMarquardtSpace::NotStarted; 00210 } 00211 00212 template<typename FunctorType, typename Scalar> 00213 LevenbergMarquardtSpace::Status 00214 LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType &x) 00215 { 00216 using std::abs; 00217 using std::sqrt; 00218 00219 eigen_assert(x.size()==n); // check the caller is not cheating us 00220 00221 /* calculate the jacobian matrix. */ 00222 Index df_ret = functor.df(x, fjac); 00223 if (df_ret<0) 00224 return LevenbergMarquardtSpace::UserAsked; 00225 if (df_ret>0) 00226 // numerical diff, we evaluated the function df_ret times 00227 nfev += df_ret; 00228 else njev++; 00229 00230 /* compute the qr factorization of the jacobian. */ 00231 wa2 = fjac.colwise().blueNorm(); 00232 ColPivHouseholderQR<JacobianType> qrfac(fjac); 00233 fjac = qrfac.matrixQR(); 00234 permutation = qrfac.colsPermutation(); 00235 00236 /* on the first iteration and if external scaling is not used, scale according */ 00237 /* to the norms of the columns of the initial jacobian. */ 00238 if (iter == 1) { 00239 if (!useExternalScaling) 00240 for (Index j = 0; j < n; ++j) 00241 diag[j] = (wa2[j]==0.)? 1. : wa2[j]; 00242 00243 /* on the first iteration, calculate the norm of the scaled x */ 00244 /* and initialize the step bound delta. */ 00245 xnorm = diag.cwiseProduct(x).stableNorm(); 00246 delta = parameters.factor * xnorm; 00247 if (delta == 0.) 00248 delta = parameters.factor; 00249 } 00250 00251 /* form (q transpose)*fvec and store the first n components in */ 00252 /* qtf. */ 00253 wa4 = fvec; 00254 wa4.applyOnTheLeft(qrfac.householderQ().adjoint()); 00255 qtf = wa4.head(n); 00256 00257 /* compute the norm of the scaled gradient. */ 00258 gnorm = 0.; 00259 if (fnorm != 0.) 00260 for (Index j = 0; j < n; ++j) 00261 if (wa2[permutation.indices()[j]] != 0.) 00262 gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]])); 00263 00264 /* test for convergence of the gradient norm. */ 00265 if (gnorm <= parameters.gtol) 00266 return LevenbergMarquardtSpace::CosinusTooSmall; 00267 00268 /* rescale if necessary. */ 00269 if (!useExternalScaling) 00270 diag = diag.cwiseMax(wa2); 00271 00272 do { 00273 00274 /* determine the levenberg-marquardt parameter. */ 00275 internal::lmpar2<Scalar>(qrfac, diag, qtf, delta, par, wa1); 00276 00277 /* store the direction p and x + p. calculate the norm of p. */ 00278 wa1 = -wa1; 00279 wa2 = x + wa1; 00280 pnorm = diag.cwiseProduct(wa1).stableNorm(); 00281 00282 /* on the first iteration, adjust the initial step bound. */ 00283 if (iter == 1) 00284 delta = (std::min)(delta,pnorm); 00285 00286 /* evaluate the function at x + p and calculate its norm. */ 00287 if ( functor(wa2, wa4) < 0) 00288 return LevenbergMarquardtSpace::UserAsked; 00289 ++nfev; 00290 fnorm1 = wa4.stableNorm(); 00291 00292 /* compute the scaled actual reduction. */ 00293 actred = -1.; 00294 if (Scalar(.1) * fnorm1 < fnorm) 00295 actred = 1. - numext::abs2(fnorm1 / fnorm); 00296 00297 /* compute the scaled predicted reduction and */ 00298 /* the scaled directional derivative. */ 00299 wa3 = fjac.template triangularView<Upper>() * (qrfac.colsPermutation().inverse() *wa1); 00300 temp1 = numext::abs2(wa3.stableNorm() / fnorm); 00301 temp2 = numext::abs2(sqrt(par) * pnorm / fnorm); 00302 prered = temp1 + temp2 / Scalar(.5); 00303 dirder = -(temp1 + temp2); 00304 00305 /* compute the ratio of the actual to the predicted */ 00306 /* reduction. */ 00307 ratio = 0.; 00308 if (prered != 0.) 00309 ratio = actred / prered; 00310 00311 /* update the step bound. */ 00312 if (ratio <= Scalar(.25)) { 00313 if (actred >= 0.) 00314 temp = Scalar(.5); 00315 if (actred < 0.) 00316 temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred); 00317 if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) 00318 temp = Scalar(.1); 00319 /* Computing MIN */ 00320 delta = temp * (std::min)(delta, pnorm / Scalar(.1)); 00321 par /= temp; 00322 } else if (!(par != 0. && ratio < Scalar(.75))) { 00323 delta = pnorm / Scalar(.5); 00324 par = Scalar(.5) * par; 00325 } 00326 00327 /* test for successful iteration. */ 00328 if (ratio >= Scalar(1e-4)) { 00329 /* successful iteration. update x, fvec, and their norms. */ 00330 x = wa2; 00331 wa2 = diag.cwiseProduct(x); 00332 fvec = wa4; 00333 xnorm = wa2.stableNorm(); 00334 fnorm = fnorm1; 00335 ++iter; 00336 } 00337 00338 /* tests for convergence. */ 00339 if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm) 00340 return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall; 00341 if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.) 00342 return LevenbergMarquardtSpace::RelativeReductionTooSmall; 00343 if (delta <= parameters.xtol * xnorm) 00344 return LevenbergMarquardtSpace::RelativeErrorTooSmall; 00345 00346 /* tests for termination and stringent tolerances. */ 00347 if (nfev >= parameters.maxfev) 00348 return LevenbergMarquardtSpace::TooManyFunctionEvaluation; 00349 if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.) 00350 return LevenbergMarquardtSpace::FtolTooSmall; 00351 if (delta <= NumTraits<Scalar>::epsilon() * xnorm) 00352 return LevenbergMarquardtSpace::XtolTooSmall; 00353 if (gnorm <= NumTraits<Scalar>::epsilon()) 00354 return LevenbergMarquardtSpace::GtolTooSmall; 00355 00356 } while (ratio < Scalar(1e-4)); 00357 00358 return LevenbergMarquardtSpace::Running; 00359 } 00360 00361 template<typename FunctorType, typename Scalar> 00362 LevenbergMarquardtSpace::Status 00363 LevenbergMarquardt<FunctorType,Scalar>::lmstr1( 00364 FVectorType &x, 00365 const Scalar tol 00366 ) 00367 { 00368 n = x.size(); 00369 m = functor.values(); 00370 00371 /* check the input parameters for errors. */ 00372 if (n <= 0 || m < n || tol < 0.) 00373 return LevenbergMarquardtSpace::ImproperInputParameters; 00374 00375 resetParameters(); 00376 parameters.ftol = tol; 00377 parameters.xtol = tol; 00378 parameters.maxfev = 100*(n+1); 00379 00380 return minimizeOptimumStorage(x); 00381 } 00382 00383 template<typename FunctorType, typename Scalar> 00384 LevenbergMarquardtSpace::Status 00385 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType &x) 00386 { 00387 n = x.size(); 00388 m = functor.values(); 00389 00390 wa1.resize(n); wa2.resize(n); wa3.resize(n); 00391 wa4.resize(m); 00392 fvec.resize(m); 00393 // Only R is stored in fjac. Q is only used to compute 'qtf', which is 00394 // Q.transpose()*rhs. qtf will be updated using givens rotation, 00395 // instead of storing them in Q. 00396 // The purpose it to only use a nxn matrix, instead of mxn here, so 00397 // that we can handle cases where m>>n : 00398 fjac.resize(n, n); 00399 if (!useExternalScaling) 00400 diag.resize(n); 00401 eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'"); 00402 qtf.resize(n); 00403 00404 /* Function Body */ 00405 nfev = 0; 00406 njev = 0; 00407 00408 /* check the input parameters for errors. */ 00409 if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.) 00410 return LevenbergMarquardtSpace::ImproperInputParameters; 00411 00412 if (useExternalScaling) 00413 for (Index j = 0; j < n; ++j) 00414 if (diag[j] <= 0.) 00415 return LevenbergMarquardtSpace::ImproperInputParameters; 00416 00417 /* evaluate the function at the starting point */ 00418 /* and calculate its norm. */ 00419 nfev = 1; 00420 if ( functor(x, fvec) < 0) 00421 return LevenbergMarquardtSpace::UserAsked; 00422 fnorm = fvec.stableNorm(); 00423 00424 /* initialize levenberg-marquardt parameter and iteration counter. */ 00425 par = 0.; 00426 iter = 1; 00427 00428 return LevenbergMarquardtSpace::NotStarted; 00429 } 00430 00431 00432 template<typename FunctorType, typename Scalar> 00433 LevenbergMarquardtSpace::Status 00434 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType &x) 00435 { 00436 using std::abs; 00437 using std::sqrt; 00438 00439 eigen_assert(x.size()==n); // check the caller is not cheating us 00440 00441 Index i, j; 00442 bool sing; 00443 00444 /* compute the qr factorization of the jacobian matrix */ 00445 /* calculated one row at a time, while simultaneously */ 00446 /* forming (q transpose)*fvec and storing the first */ 00447 /* n components in qtf. */ 00448 qtf.fill(0.); 00449 fjac.fill(0.); 00450 Index rownb = 2; 00451 for (i = 0; i < m; ++i) { 00452 if (functor.df(x, wa3, rownb) < 0) return LevenbergMarquardtSpace::UserAsked; 00453 internal::rwupdt<Scalar>(fjac, wa3, qtf, fvec[i]); 00454 ++rownb; 00455 } 00456 ++njev; 00457 00458 /* if the jacobian is rank deficient, call qrfac to */ 00459 /* reorder its columns and update the components of qtf. */ 00460 sing = false; 00461 for (j = 0; j < n; ++j) { 00462 if (fjac(j,j) == 0.) 00463 sing = true; 00464 wa2[j] = fjac.col(j).head(j).stableNorm(); 00465 } 00466 permutation.setIdentity(n); 00467 if (sing) { 00468 wa2 = fjac.colwise().blueNorm(); 00469 // TODO We have no unit test covering this code path, do not modify 00470 // until it is carefully tested 00471 ColPivHouseholderQR<JacobianType> qrfac(fjac); 00472 fjac = qrfac.matrixQR(); 00473 wa1 = fjac.diagonal(); 00474 fjac.diagonal() = qrfac.hCoeffs(); 00475 permutation = qrfac.colsPermutation(); 00476 // TODO : avoid this: 00477 for(Index ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors 00478 00479 for (j = 0; j < n; ++j) { 00480 if (fjac(j,j) != 0.) { 00481 sum = 0.; 00482 for (i = j; i < n; ++i) 00483 sum += fjac(i,j) * qtf[i]; 00484 temp = -sum / fjac(j,j); 00485 for (i = j; i < n; ++i) 00486 qtf[i] += fjac(i,j) * temp; 00487 } 00488 fjac(j,j) = wa1[j]; 00489 } 00490 } 00491 00492 /* on the first iteration and if external scaling is not used, scale according */ 00493 /* to the norms of the columns of the initial jacobian. */ 00494 if (iter == 1) { 00495 if (!useExternalScaling) 00496 for (j = 0; j < n; ++j) 00497 diag[j] = (wa2[j]==0.)? 1. : wa2[j]; 00498 00499 /* on the first iteration, calculate the norm of the scaled x */ 00500 /* and initialize the step bound delta. */ 00501 xnorm = diag.cwiseProduct(x).stableNorm(); 00502 delta = parameters.factor * xnorm; 00503 if (delta == 0.) 00504 delta = parameters.factor; 00505 } 00506 00507 /* compute the norm of the scaled gradient. */ 00508 gnorm = 0.; 00509 if (fnorm != 0.) 00510 for (j = 0; j < n; ++j) 00511 if (wa2[permutation.indices()[j]] != 0.) 00512 gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]])); 00513 00514 /* test for convergence of the gradient norm. */ 00515 if (gnorm <= parameters.gtol) 00516 return LevenbergMarquardtSpace::CosinusTooSmall; 00517 00518 /* rescale if necessary. */ 00519 if (!useExternalScaling) 00520 diag = diag.cwiseMax(wa2); 00521 00522 do { 00523 00524 /* determine the levenberg-marquardt parameter. */ 00525 internal::lmpar<Scalar>(fjac, permutation.indices(), diag, qtf, delta, par, wa1); 00526 00527 /* store the direction p and x + p. calculate the norm of p. */ 00528 wa1 = -wa1; 00529 wa2 = x + wa1; 00530 pnorm = diag.cwiseProduct(wa1).stableNorm(); 00531 00532 /* on the first iteration, adjust the initial step bound. */ 00533 if (iter == 1) 00534 delta = (std::min)(delta,pnorm); 00535 00536 /* evaluate the function at x + p and calculate its norm. */ 00537 if ( functor(wa2, wa4) < 0) 00538 return LevenbergMarquardtSpace::UserAsked; 00539 ++nfev; 00540 fnorm1 = wa4.stableNorm(); 00541 00542 /* compute the scaled actual reduction. */ 00543 actred = -1.; 00544 if (Scalar(.1) * fnorm1 < fnorm) 00545 actred = 1. - numext::abs2(fnorm1 / fnorm); 00546 00547 /* compute the scaled predicted reduction and */ 00548 /* the scaled directional derivative. */ 00549 wa3 = fjac.topLeftCorner(n,n).template triangularView<Upper>() * (permutation.inverse() * wa1); 00550 temp1 = numext::abs2(wa3.stableNorm() / fnorm); 00551 temp2 = numext::abs2(sqrt(par) * pnorm / fnorm); 00552 prered = temp1 + temp2 / Scalar(.5); 00553 dirder = -(temp1 + temp2); 00554 00555 /* compute the ratio of the actual to the predicted */ 00556 /* reduction. */ 00557 ratio = 0.; 00558 if (prered != 0.) 00559 ratio = actred / prered; 00560 00561 /* update the step bound. */ 00562 if (ratio <= Scalar(.25)) { 00563 if (actred >= 0.) 00564 temp = Scalar(.5); 00565 if (actred < 0.) 00566 temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred); 00567 if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) 00568 temp = Scalar(.1); 00569 /* Computing MIN */ 00570 delta = temp * (std::min)(delta, pnorm / Scalar(.1)); 00571 par /= temp; 00572 } else if (!(par != 0. && ratio < Scalar(.75))) { 00573 delta = pnorm / Scalar(.5); 00574 par = Scalar(.5) * par; 00575 } 00576 00577 /* test for successful iteration. */ 00578 if (ratio >= Scalar(1e-4)) { 00579 /* successful iteration. update x, fvec, and their norms. */ 00580 x = wa2; 00581 wa2 = diag.cwiseProduct(x); 00582 fvec = wa4; 00583 xnorm = wa2.stableNorm(); 00584 fnorm = fnorm1; 00585 ++iter; 00586 } 00587 00588 /* tests for convergence. */ 00589 if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm) 00590 return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall; 00591 if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.) 00592 return LevenbergMarquardtSpace::RelativeReductionTooSmall; 00593 if (delta <= parameters.xtol * xnorm) 00594 return LevenbergMarquardtSpace::RelativeErrorTooSmall; 00595 00596 /* tests for termination and stringent tolerances. */ 00597 if (nfev >= parameters.maxfev) 00598 return LevenbergMarquardtSpace::TooManyFunctionEvaluation; 00599 if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.) 00600 return LevenbergMarquardtSpace::FtolTooSmall; 00601 if (delta <= NumTraits<Scalar>::epsilon() * xnorm) 00602 return LevenbergMarquardtSpace::XtolTooSmall; 00603 if (gnorm <= NumTraits<Scalar>::epsilon()) 00604 return LevenbergMarquardtSpace::GtolTooSmall; 00605 00606 } while (ratio < Scalar(1e-4)); 00607 00608 return LevenbergMarquardtSpace::Running; 00609 } 00610 00611 template<typename FunctorType, typename Scalar> 00612 LevenbergMarquardtSpace::Status 00613 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType &x) 00614 { 00615 LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x); 00616 if (status==LevenbergMarquardtSpace::ImproperInputParameters) 00617 return status; 00618 do { 00619 status = minimizeOptimumStorageOneStep(x); 00620 } while (status==LevenbergMarquardtSpace::Running); 00621 return status; 00622 } 00623 00624 template<typename FunctorType, typename Scalar> 00625 LevenbergMarquardtSpace::Status 00626 LevenbergMarquardt<FunctorType,Scalar>::lmdif1( 00627 FunctorType &functor, 00628 FVectorType &x, 00629 Index *nfev, 00630 const Scalar tol 00631 ) 00632 { 00633 Index n = x.size(); 00634 Index m = functor.values(); 00635 00636 /* check the input parameters for errors. */ 00637 if (n <= 0 || m < n || tol < 0.) 00638 return LevenbergMarquardtSpace::ImproperInputParameters; 00639 00640 NumericalDiff<FunctorType> numDiff(functor); 00641 // embedded LevenbergMarquardt 00642 LevenbergMarquardt<NumericalDiff<FunctorType>, Scalar > lm(numDiff); 00643 lm.parameters.ftol = tol; 00644 lm.parameters.xtol = tol; 00645 lm.parameters.maxfev = 200*(n+1); 00646 00647 LevenbergMarquardtSpace::Status info = LevenbergMarquardtSpace::Status(lm.minimize(x)); 00648 if (nfev) 00649 * nfev = lm.nfev; 00650 return info; 00651 } 00652 00653 } // end namespace Eigen 00654 00655 #endif // EIGEN_LEVENBERGMARQUARDT__H 00656 00657 //vim: ai ts=4 sts=4 et sw=4