Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:57:06 

0001 // -*- coding: utf-8
0002 // vim: set fileencoding=utf-8
0003 
0004 // This file is part of Eigen, a lightweight C++ template library
0005 // for linear algebra.
0006 //
0007 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
0008 //
0009 // This Source Code Form is subject to the terms of the Mozilla
0010 // Public License v. 2.0. If a copy of the MPL was not distributed
0011 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0012 
0013 #ifndef EIGEN_LEVENBERGMARQUARDT__H
0014 #define EIGEN_LEVENBERGMARQUARDT__H
0015 
0016 namespace Eigen { 
0017 
0018 namespace LevenbergMarquardtSpace {
0019     enum Status {
0020         NotStarted = -2,
0021         Running = -1,
0022         ImproperInputParameters = 0,
0023         RelativeReductionTooSmall = 1,
0024         RelativeErrorTooSmall = 2,
0025         RelativeErrorAndReductionTooSmall = 3,
0026         CosinusTooSmall = 4,
0027         TooManyFunctionEvaluation = 5,
0028         FtolTooSmall = 6,
0029         XtolTooSmall = 7,
0030         GtolTooSmall = 8,
0031         UserAsked = 9
0032     };
0033 }
0034 
0035 
0036 
0037 /**
0038   * \ingroup NonLinearOptimization_Module
0039   * \brief Performs non linear optimization over a non-linear function,
0040   * using a variant of the Levenberg Marquardt algorithm.
0041   *
0042   * Check wikipedia for more information.
0043   * http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
0044   */
0045 template<typename FunctorType, typename Scalar=double>
0046 class LevenbergMarquardt
0047 {
0048     static Scalar sqrt_epsilon()
0049     {
0050       using std::sqrt;
0051       return sqrt(NumTraits<Scalar>::epsilon());
0052     }
0053     
0054 public:
0055     LevenbergMarquardt(FunctorType &_functor)
0056         : functor(_functor) { nfev = njev = iter = 0;  fnorm = gnorm = 0.; useExternalScaling=false; }
0057 
0058     typedef DenseIndex Index;
0059     
0060     struct Parameters {
0061         Parameters()
0062             : factor(Scalar(100.))
0063             , maxfev(400)
0064             , ftol(sqrt_epsilon())
0065             , xtol(sqrt_epsilon())
0066             , gtol(Scalar(0.))
0067             , epsfcn(Scalar(0.)) {}
0068         Scalar factor;
0069         Index maxfev;   // maximum number of function evaluation
0070         Scalar ftol;
0071         Scalar xtol;
0072         Scalar gtol;
0073         Scalar epsfcn;
0074     };
0075 
0076     typedef Matrix< Scalar, Dynamic, 1 > FVectorType;
0077     typedef Matrix< Scalar, Dynamic, Dynamic > JacobianType;
0078 
0079     LevenbergMarquardtSpace::Status lmder1(
0080             FVectorType &x,
0081             const Scalar tol = sqrt_epsilon()
0082             );
0083 
0084     LevenbergMarquardtSpace::Status minimize(FVectorType &x);
0085     LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x);
0086     LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x);
0087 
0088     static LevenbergMarquardtSpace::Status lmdif1(
0089             FunctorType &functor,
0090             FVectorType &x,
0091             Index *nfev,
0092             const Scalar tol = sqrt_epsilon()
0093             );
0094 
0095     LevenbergMarquardtSpace::Status lmstr1(
0096             FVectorType  &x,
0097             const Scalar tol = sqrt_epsilon()
0098             );
0099 
0100     LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType  &x);
0101     LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType  &x);
0102     LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType  &x);
0103 
0104     void resetParameters(void) { parameters = Parameters(); }
0105 
0106     Parameters parameters;
0107     FVectorType  fvec, qtf, diag;
0108     JacobianType fjac;
0109     PermutationMatrix<Dynamic,Dynamic> permutation;
0110     Index nfev;
0111     Index njev;
0112     Index iter;
0113     Scalar fnorm, gnorm;
0114     bool useExternalScaling; 
0115 
0116     Scalar lm_param(void) { return par; }
0117 private:
0118     
0119     FunctorType &functor;
0120     Index n;
0121     Index m;
0122     FVectorType wa1, wa2, wa3, wa4;
0123 
0124     Scalar par, sum;
0125     Scalar temp, temp1, temp2;
0126     Scalar delta;
0127     Scalar ratio;
0128     Scalar pnorm, xnorm, fnorm1, actred, dirder, prered;
0129 
0130     LevenbergMarquardt& operator=(const LevenbergMarquardt&);
0131 };
0132 
0133 template<typename FunctorType, typename Scalar>
0134 LevenbergMarquardtSpace::Status
0135 LevenbergMarquardt<FunctorType,Scalar>::lmder1(
0136         FVectorType  &x,
0137         const Scalar tol
0138         )
0139 {
0140     n = x.size();
0141     m = functor.values();
0142 
0143     /* check the input parameters for errors. */
0144     if (n <= 0 || m < n || tol < 0.)
0145         return LevenbergMarquardtSpace::ImproperInputParameters;
0146 
0147     resetParameters();
0148     parameters.ftol = tol;
0149     parameters.xtol = tol;
0150     parameters.maxfev = 100*(n+1);
0151 
0152     return minimize(x);
0153 }
0154 
0155 
0156 template<typename FunctorType, typename Scalar>
0157 LevenbergMarquardtSpace::Status
0158 LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType  &x)
0159 {
0160     LevenbergMarquardtSpace::Status status = minimizeInit(x);
0161     if (status==LevenbergMarquardtSpace::ImproperInputParameters)
0162         return status;
0163     do {
0164         status = minimizeOneStep(x);
0165     } while (status==LevenbergMarquardtSpace::Running);
0166     return status;
0167 }
0168 
0169 template<typename FunctorType, typename Scalar>
0170 LevenbergMarquardtSpace::Status
0171 LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType  &x)
0172 {
0173     n = x.size();
0174     m = functor.values();
0175 
0176     wa1.resize(n); wa2.resize(n); wa3.resize(n);
0177     wa4.resize(m);
0178     fvec.resize(m);
0179     fjac.resize(m, n);
0180     if (!useExternalScaling)
0181         diag.resize(n);
0182     eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'");
0183     qtf.resize(n);
0184 
0185     /* Function Body */
0186     nfev = 0;
0187     njev = 0;
0188 
0189     /*     check the input parameters for errors. */
0190     if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
0191         return LevenbergMarquardtSpace::ImproperInputParameters;
0192 
0193     if (useExternalScaling)
0194         for (Index j = 0; j < n; ++j)
0195             if (diag[j] <= 0.)
0196                 return LevenbergMarquardtSpace::ImproperInputParameters;
0197 
0198     /*     evaluate the function at the starting point */
0199     /*     and calculate its norm. */
0200     nfev = 1;
0201     if ( functor(x, fvec) < 0)
0202         return LevenbergMarquardtSpace::UserAsked;
0203     fnorm = fvec.stableNorm();
0204 
0205     /*     initialize levenberg-marquardt parameter and iteration counter. */
0206     par = 0.;
0207     iter = 1;
0208 
0209     return LevenbergMarquardtSpace::NotStarted;
0210 }
0211 
0212 template<typename FunctorType, typename Scalar>
0213 LevenbergMarquardtSpace::Status
0214 LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType  &x)
0215 {
0216     using std::abs;
0217     using std::sqrt;
0218 
0219     eigen_assert(x.size()==n); // check the caller is not cheating us
0220 
0221     /* calculate the jacobian matrix. */
0222     Index df_ret = functor.df(x, fjac);
0223     if (df_ret<0)
0224         return LevenbergMarquardtSpace::UserAsked;
0225     if (df_ret>0)
0226         // numerical diff, we evaluated the function df_ret times
0227         nfev += df_ret;
0228     else njev++;
0229 
0230     /* compute the qr factorization of the jacobian. */
0231     wa2 = fjac.colwise().blueNorm();
0232     ColPivHouseholderQR<JacobianType> qrfac(fjac);
0233     fjac = qrfac.matrixQR();
0234     permutation = qrfac.colsPermutation();
0235 
0236     /* on the first iteration and if external scaling is not used, scale according */
0237     /* to the norms of the columns of the initial jacobian. */
0238     if (iter == 1) {
0239         if (!useExternalScaling)
0240             for (Index j = 0; j < n; ++j)
0241                 diag[j] = (wa2[j]==0.)? 1. : wa2[j];
0242 
0243         /* on the first iteration, calculate the norm of the scaled x */
0244         /* and initialize the step bound delta. */
0245         xnorm = diag.cwiseProduct(x).stableNorm();
0246         delta = parameters.factor * xnorm;
0247         if (delta == 0.)
0248             delta = parameters.factor;
0249     }
0250 
0251     /* form (q transpose)*fvec and store the first n components in */
0252     /* qtf. */
0253     wa4 = fvec;
0254     wa4.applyOnTheLeft(qrfac.householderQ().adjoint());
0255     qtf = wa4.head(n);
0256 
0257     /* compute the norm of the scaled gradient. */
0258     gnorm = 0.;
0259     if (fnorm != 0.)
0260         for (Index j = 0; j < n; ++j)
0261             if (wa2[permutation.indices()[j]] != 0.)
0262                 gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]]));
0263 
0264     /* test for convergence of the gradient norm. */
0265     if (gnorm <= parameters.gtol)
0266         return LevenbergMarquardtSpace::CosinusTooSmall;
0267 
0268     /* rescale if necessary. */
0269     if (!useExternalScaling)
0270         diag = diag.cwiseMax(wa2);
0271 
0272     do {
0273 
0274         /* determine the levenberg-marquardt parameter. */
0275         internal::lmpar2<Scalar>(qrfac, diag, qtf, delta, par, wa1);
0276 
0277         /* store the direction p and x + p. calculate the norm of p. */
0278         wa1 = -wa1;
0279         wa2 = x + wa1;
0280         pnorm = diag.cwiseProduct(wa1).stableNorm();
0281 
0282         /* on the first iteration, adjust the initial step bound. */
0283         if (iter == 1)
0284             delta = (std::min)(delta,pnorm);
0285 
0286         /* evaluate the function at x + p and calculate its norm. */
0287         if ( functor(wa2, wa4) < 0)
0288             return LevenbergMarquardtSpace::UserAsked;
0289         ++nfev;
0290         fnorm1 = wa4.stableNorm();
0291 
0292         /* compute the scaled actual reduction. */
0293         actred = -1.;
0294         if (Scalar(.1) * fnorm1 < fnorm)
0295             actred = 1. - numext::abs2(fnorm1 / fnorm);
0296 
0297         /* compute the scaled predicted reduction and */
0298         /* the scaled directional derivative. */
0299         wa3 = fjac.template triangularView<Upper>() * (qrfac.colsPermutation().inverse() *wa1);
0300         temp1 = numext::abs2(wa3.stableNorm() / fnorm);
0301         temp2 = numext::abs2(sqrt(par) * pnorm / fnorm);
0302         prered = temp1 + temp2 / Scalar(.5);
0303         dirder = -(temp1 + temp2);
0304 
0305         /* compute the ratio of the actual to the predicted */
0306         /* reduction. */
0307         ratio = 0.;
0308         if (prered != 0.)
0309             ratio = actred / prered;
0310 
0311         /* update the step bound. */
0312         if (ratio <= Scalar(.25)) {
0313             if (actred >= 0.)
0314                 temp = Scalar(.5);
0315             if (actred < 0.)
0316                 temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
0317             if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
0318                 temp = Scalar(.1);
0319             /* Computing MIN */
0320             delta = temp * (std::min)(delta, pnorm / Scalar(.1));
0321             par /= temp;
0322         } else if (!(par != 0. && ratio < Scalar(.75))) {
0323             delta = pnorm / Scalar(.5);
0324             par = Scalar(.5) * par;
0325         }
0326 
0327         /* test for successful iteration. */
0328         if (ratio >= Scalar(1e-4)) {
0329             /* successful iteration. update x, fvec, and their norms. */
0330             x = wa2;
0331             wa2 = diag.cwiseProduct(x);
0332             fvec = wa4;
0333             xnorm = wa2.stableNorm();
0334             fnorm = fnorm1;
0335             ++iter;
0336         }
0337 
0338         /* tests for convergence. */
0339         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
0340             return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
0341         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
0342             return LevenbergMarquardtSpace::RelativeReductionTooSmall;
0343         if (delta <= parameters.xtol * xnorm)
0344             return LevenbergMarquardtSpace::RelativeErrorTooSmall;
0345 
0346         /* tests for termination and stringent tolerances. */
0347         if (nfev >= parameters.maxfev)
0348             return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
0349         if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.)
0350             return LevenbergMarquardtSpace::FtolTooSmall;
0351         if (delta <= NumTraits<Scalar>::epsilon() * xnorm)
0352             return LevenbergMarquardtSpace::XtolTooSmall;
0353         if (gnorm <= NumTraits<Scalar>::epsilon())
0354             return LevenbergMarquardtSpace::GtolTooSmall;
0355 
0356     } while (ratio < Scalar(1e-4));
0357 
0358     return LevenbergMarquardtSpace::Running;
0359 }
0360 
0361 template<typename FunctorType, typename Scalar>
0362 LevenbergMarquardtSpace::Status
0363 LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
0364         FVectorType  &x,
0365         const Scalar tol
0366         )
0367 {
0368     n = x.size();
0369     m = functor.values();
0370 
0371     /* check the input parameters for errors. */
0372     if (n <= 0 || m < n || tol < 0.)
0373         return LevenbergMarquardtSpace::ImproperInputParameters;
0374 
0375     resetParameters();
0376     parameters.ftol = tol;
0377     parameters.xtol = tol;
0378     parameters.maxfev = 100*(n+1);
0379 
0380     return minimizeOptimumStorage(x);
0381 }
0382 
0383 template<typename FunctorType, typename Scalar>
0384 LevenbergMarquardtSpace::Status
0385 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType  &x)
0386 {
0387     n = x.size();
0388     m = functor.values();
0389 
0390     wa1.resize(n); wa2.resize(n); wa3.resize(n);
0391     wa4.resize(m);
0392     fvec.resize(m);
0393     // Only R is stored in fjac. Q is only used to compute 'qtf', which is
0394     // Q.transpose()*rhs. qtf will be updated using givens rotation,
0395     // instead of storing them in Q.
0396     // The purpose it to only use a nxn matrix, instead of mxn here, so
0397     // that we can handle cases where m>>n :
0398     fjac.resize(n, n);
0399     if (!useExternalScaling)
0400         diag.resize(n);
0401     eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'");
0402     qtf.resize(n);
0403 
0404     /* Function Body */
0405     nfev = 0;
0406     njev = 0;
0407 
0408     /*     check the input parameters for errors. */
0409     if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
0410         return LevenbergMarquardtSpace::ImproperInputParameters;
0411 
0412     if (useExternalScaling)
0413         for (Index j = 0; j < n; ++j)
0414             if (diag[j] <= 0.)
0415                 return LevenbergMarquardtSpace::ImproperInputParameters;
0416 
0417     /*     evaluate the function at the starting point */
0418     /*     and calculate its norm. */
0419     nfev = 1;
0420     if ( functor(x, fvec) < 0)
0421         return LevenbergMarquardtSpace::UserAsked;
0422     fnorm = fvec.stableNorm();
0423 
0424     /*     initialize levenberg-marquardt parameter and iteration counter. */
0425     par = 0.;
0426     iter = 1;
0427 
0428     return LevenbergMarquardtSpace::NotStarted;
0429 }
0430 
0431 
0432 template<typename FunctorType, typename Scalar>
0433 LevenbergMarquardtSpace::Status
0434 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType  &x)
0435 {
0436     using std::abs;
0437     using std::sqrt;
0438     
0439     eigen_assert(x.size()==n); // check the caller is not cheating us
0440 
0441     Index i, j;
0442     bool sing;
0443 
0444     /* compute the qr factorization of the jacobian matrix */
0445     /* calculated one row at a time, while simultaneously */
0446     /* forming (q transpose)*fvec and storing the first */
0447     /* n components in qtf. */
0448     qtf.fill(0.);
0449     fjac.fill(0.);
0450     Index rownb = 2;
0451     for (i = 0; i < m; ++i) {
0452         if (functor.df(x, wa3, rownb) < 0) return LevenbergMarquardtSpace::UserAsked;
0453         internal::rwupdt<Scalar>(fjac, wa3, qtf, fvec[i]);
0454         ++rownb;
0455     }
0456     ++njev;
0457 
0458     /* if the jacobian is rank deficient, call qrfac to */
0459     /* reorder its columns and update the components of qtf. */
0460     sing = false;
0461     for (j = 0; j < n; ++j) {
0462         if (fjac(j,j) == 0.)
0463             sing = true;
0464         wa2[j] = fjac.col(j).head(j).stableNorm();
0465     }
0466     permutation.setIdentity(n);
0467     if (sing) {
0468         wa2 = fjac.colwise().blueNorm();
0469         // TODO We have no unit test covering this code path, do not modify
0470         // until it is carefully tested
0471         ColPivHouseholderQR<JacobianType> qrfac(fjac);
0472         fjac = qrfac.matrixQR();
0473         wa1 = fjac.diagonal();
0474         fjac.diagonal() = qrfac.hCoeffs();
0475         permutation = qrfac.colsPermutation();
0476         // TODO : avoid this:
0477         for(Index ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors
0478 
0479         for (j = 0; j < n; ++j) {
0480             if (fjac(j,j) != 0.) {
0481                 sum = 0.;
0482                 for (i = j; i < n; ++i)
0483                     sum += fjac(i,j) * qtf[i];
0484                 temp = -sum / fjac(j,j);
0485                 for (i = j; i < n; ++i)
0486                     qtf[i] += fjac(i,j) * temp;
0487             }
0488             fjac(j,j) = wa1[j];
0489         }
0490     }
0491 
0492     /* on the first iteration and if external scaling is not used, scale according */
0493     /* to the norms of the columns of the initial jacobian. */
0494     if (iter == 1) {
0495         if (!useExternalScaling)
0496             for (j = 0; j < n; ++j)
0497                 diag[j] = (wa2[j]==0.)? 1. : wa2[j];
0498 
0499         /* on the first iteration, calculate the norm of the scaled x */
0500         /* and initialize the step bound delta. */
0501         xnorm = diag.cwiseProduct(x).stableNorm();
0502         delta = parameters.factor * xnorm;
0503         if (delta == 0.)
0504             delta = parameters.factor;
0505     }
0506 
0507     /* compute the norm of the scaled gradient. */
0508     gnorm = 0.;
0509     if (fnorm != 0.)
0510         for (j = 0; j < n; ++j)
0511             if (wa2[permutation.indices()[j]] != 0.)
0512                 gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]]));
0513 
0514     /* test for convergence of the gradient norm. */
0515     if (gnorm <= parameters.gtol)
0516         return LevenbergMarquardtSpace::CosinusTooSmall;
0517 
0518     /* rescale if necessary. */
0519     if (!useExternalScaling)
0520         diag = diag.cwiseMax(wa2);
0521 
0522     do {
0523 
0524         /* determine the levenberg-marquardt parameter. */
0525         internal::lmpar<Scalar>(fjac, permutation.indices(), diag, qtf, delta, par, wa1);
0526 
0527         /* store the direction p and x + p. calculate the norm of p. */
0528         wa1 = -wa1;
0529         wa2 = x + wa1;
0530         pnorm = diag.cwiseProduct(wa1).stableNorm();
0531 
0532         /* on the first iteration, adjust the initial step bound. */
0533         if (iter == 1)
0534             delta = (std::min)(delta,pnorm);
0535 
0536         /* evaluate the function at x + p and calculate its norm. */
0537         if ( functor(wa2, wa4) < 0)
0538             return LevenbergMarquardtSpace::UserAsked;
0539         ++nfev;
0540         fnorm1 = wa4.stableNorm();
0541 
0542         /* compute the scaled actual reduction. */
0543         actred = -1.;
0544         if (Scalar(.1) * fnorm1 < fnorm)
0545             actred = 1. - numext::abs2(fnorm1 / fnorm);
0546 
0547         /* compute the scaled predicted reduction and */
0548         /* the scaled directional derivative. */
0549         wa3 = fjac.topLeftCorner(n,n).template triangularView<Upper>() * (permutation.inverse() * wa1);
0550         temp1 = numext::abs2(wa3.stableNorm() / fnorm);
0551         temp2 = numext::abs2(sqrt(par) * pnorm / fnorm);
0552         prered = temp1 + temp2 / Scalar(.5);
0553         dirder = -(temp1 + temp2);
0554 
0555         /* compute the ratio of the actual to the predicted */
0556         /* reduction. */
0557         ratio = 0.;
0558         if (prered != 0.)
0559             ratio = actred / prered;
0560 
0561         /* update the step bound. */
0562         if (ratio <= Scalar(.25)) {
0563             if (actred >= 0.)
0564                 temp = Scalar(.5);
0565             if (actred < 0.)
0566                 temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
0567             if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
0568                 temp = Scalar(.1);
0569             /* Computing MIN */
0570             delta = temp * (std::min)(delta, pnorm / Scalar(.1));
0571             par /= temp;
0572         } else if (!(par != 0. && ratio < Scalar(.75))) {
0573             delta = pnorm / Scalar(.5);
0574             par = Scalar(.5) * par;
0575         }
0576 
0577         /* test for successful iteration. */
0578         if (ratio >= Scalar(1e-4)) {
0579             /* successful iteration. update x, fvec, and their norms. */
0580             x = wa2;
0581             wa2 = diag.cwiseProduct(x);
0582             fvec = wa4;
0583             xnorm = wa2.stableNorm();
0584             fnorm = fnorm1;
0585             ++iter;
0586         }
0587 
0588         /* tests for convergence. */
0589         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
0590             return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
0591         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
0592             return LevenbergMarquardtSpace::RelativeReductionTooSmall;
0593         if (delta <= parameters.xtol * xnorm)
0594             return LevenbergMarquardtSpace::RelativeErrorTooSmall;
0595 
0596         /* tests for termination and stringent tolerances. */
0597         if (nfev >= parameters.maxfev)
0598             return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
0599         if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.)
0600             return LevenbergMarquardtSpace::FtolTooSmall;
0601         if (delta <= NumTraits<Scalar>::epsilon() * xnorm)
0602             return LevenbergMarquardtSpace::XtolTooSmall;
0603         if (gnorm <= NumTraits<Scalar>::epsilon())
0604             return LevenbergMarquardtSpace::GtolTooSmall;
0605 
0606     } while (ratio < Scalar(1e-4));
0607 
0608     return LevenbergMarquardtSpace::Running;
0609 }
0610 
0611 template<typename FunctorType, typename Scalar>
0612 LevenbergMarquardtSpace::Status
0613 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType  &x)
0614 {
0615     LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x);
0616     if (status==LevenbergMarquardtSpace::ImproperInputParameters)
0617         return status;
0618     do {
0619         status = minimizeOptimumStorageOneStep(x);
0620     } while (status==LevenbergMarquardtSpace::Running);
0621     return status;
0622 }
0623 
0624 template<typename FunctorType, typename Scalar>
0625 LevenbergMarquardtSpace::Status
0626 LevenbergMarquardt<FunctorType,Scalar>::lmdif1(
0627         FunctorType &functor,
0628         FVectorType  &x,
0629         Index *nfev,
0630         const Scalar tol
0631         )
0632 {
0633     Index n = x.size();
0634     Index m = functor.values();
0635 
0636     /* check the input parameters for errors. */
0637     if (n <= 0 || m < n || tol < 0.)
0638         return LevenbergMarquardtSpace::ImproperInputParameters;
0639 
0640     NumericalDiff<FunctorType> numDiff(functor);
0641     // embedded LevenbergMarquardt
0642     LevenbergMarquardt<NumericalDiff<FunctorType>, Scalar > lm(numDiff);
0643     lm.parameters.ftol = tol;
0644     lm.parameters.xtol = tol;
0645     lm.parameters.maxfev = 200*(n+1);
0646 
0647     LevenbergMarquardtSpace::Status info = LevenbergMarquardtSpace::Status(lm.minimize(x));
0648     if (nfev)
0649         * nfev = lm.nfev;
0650     return info;
0651 }
0652 
0653 } // end namespace Eigen
0654 
0655 #endif // EIGEN_LEVENBERGMARQUARDT__H
0656 
0657 //vim: ai ts=4 sts=4 et sw=4
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!