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