NonLinearOptimization.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
00005 
00006 #include <stdio.h>
00007 
00008 #include "main.h"
00009 #include <unsupported/Eigen/NonLinearOptimization>
00010 
00011 // This disables some useless Warnings on MSVC.
00012 // It is intended to be done for this test only.
00013 #include <Eigen/src/Core/util/DisableStupidWarnings.h>
00014 
00015 int fcn_chkder(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag)
00016 {
00017     /*      subroutine fcn for chkder example. */
00018 
00019     int i;
00020     assert(15 ==  fvec.size());
00021     assert(3 ==  x.size());
00022     double tmp1, tmp2, tmp3, tmp4;
00023     static const double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
00024         3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
00025 
00026 
00027     if (iflag == 0)
00028         return 0;
00029 
00030     if (iflag != 2)
00031         for (i=0; i<15; i++) {
00032             tmp1 = i+1;
00033             tmp2 = 16-i-1;
00034             tmp3 = tmp1;
00035             if (i >= 8) tmp3 = tmp2;
00036             fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
00037         }
00038     else {
00039         for (i = 0; i < 15; i++) {
00040             tmp1 = i+1;
00041             tmp2 = 16-i-1;
00042 
00043             /* error introduced into next statement for illustration. */
00044             /* corrected statement should read    tmp3 = tmp1 . */
00045 
00046             tmp3 = tmp2;
00047             if (i >= 8) tmp3 = tmp2;
00048             tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4=tmp4*tmp4;
00049             fjac(i,0) = -1.;
00050             fjac(i,1) = tmp1*tmp2/tmp4;
00051             fjac(i,2) = tmp1*tmp3/tmp4;
00052         }
00053     }
00054     return 0;
00055 }
00056 
00057 
00058 void testChkder()
00059 {
00060   const int m=15, n=3;
00061   VectorXd x(n), fvec(m), xp, fvecp(m), err;
00062   MatrixXd fjac(m,n);
00063   VectorXi ipvt;
00064 
00065   /*      the following values should be suitable for */
00066   /*      checking the jacobian matrix. */
00067   x << 9.2e-1, 1.3e-1, 5.4e-1;
00068 
00069   internal::chkder(x, fvec, fjac, xp, fvecp, 1, err);
00070   fcn_chkder(x, fvec, fjac, 1);
00071   fcn_chkder(x, fvec, fjac, 2);
00072   fcn_chkder(xp, fvecp, fjac, 1);
00073   internal::chkder(x, fvec, fjac, xp, fvecp, 2, err);
00074 
00075   fvecp -= fvec;
00076 
00077   // check those
00078   VectorXd fvec_ref(m), fvecp_ref(m), err_ref(m);
00079   fvec_ref <<
00080       -1.181606, -1.429655, -1.606344,
00081       -1.745269, -1.840654, -1.921586,
00082       -1.984141, -2.022537, -2.468977,
00083       -2.827562, -3.473582, -4.437612,
00084       -6.047662, -9.267761, -18.91806;
00085   fvecp_ref <<
00086       -7.724666e-09, -3.432406e-09, -2.034843e-10,
00087       2.313685e-09,  4.331078e-09,  5.984096e-09,
00088       7.363281e-09,   8.53147e-09,  1.488591e-08,
00089       2.33585e-08,  3.522012e-08,  5.301255e-08,
00090       8.26666e-08,  1.419747e-07,   3.19899e-07;
00091   err_ref <<
00092       0.1141397,  0.09943516,  0.09674474,
00093       0.09980447,  0.1073116, 0.1220445,
00094       0.1526814, 1, 1,
00095       1, 1, 1,
00096       1, 1, 1;
00097 
00098   VERIFY_IS_APPROX(fvec, fvec_ref);
00099   VERIFY_IS_APPROX(fvecp, fvecp_ref);
00100   VERIFY_IS_APPROX(err, err_ref);
00101 }
00102 
00103 // Generic functor
00104 template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
00105 struct Functor
00106 {
00107   typedef _Scalar Scalar;
00108   enum {
00109     InputsAtCompileTime = NX,
00110     ValuesAtCompileTime = NY
00111   };
00112   typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
00113   typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
00114   typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
00115 
00116   const int m_inputs, m_values;
00117 
00118   Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
00119   Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
00120 
00121   int inputs() const { return m_inputs; }
00122   int values() const { return m_values; }
00123 
00124   // you should define that in the subclass :
00125 //  void operator() (const InputType& x, ValueType* v, JacobianType* _j=0) const;
00126 };
00127 
00128 struct lmder_functor : Functor<double>
00129 {
00130     lmder_functor(void): Functor<double>(3,15) {}
00131     int operator()(const VectorXd &x, VectorXd &fvec) const
00132     {
00133         double tmp1, tmp2, tmp3;
00134         static const double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
00135             3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
00136 
00137         for (int i = 0; i < values(); i++)
00138         {
00139             tmp1 = i+1;
00140             tmp2 = 16 - i - 1;
00141             tmp3 = (i>=8)? tmp2 : tmp1;
00142             fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
00143         }
00144         return 0;
00145     }
00146 
00147     int df(const VectorXd &x, MatrixXd &fjac) const
00148     {
00149         double tmp1, tmp2, tmp3, tmp4;
00150         for (int i = 0; i < values(); i++)
00151         {
00152             tmp1 = i+1;
00153             tmp2 = 16 - i - 1;
00154             tmp3 = (i>=8)? tmp2 : tmp1;
00155             tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
00156             fjac(i,0) = -1;
00157             fjac(i,1) = tmp1*tmp2/tmp4;
00158             fjac(i,2) = tmp1*tmp3/tmp4;
00159         }
00160         return 0;
00161     }
00162 };
00163 
00164 void testLmder1()
00165 {
00166   int n=3, info;
00167 
00168   VectorXd x;
00169 
00170   /* the following starting values provide a rough fit. */
00171   x.setConstant(n, 1.);
00172 
00173   // do the computation
00174   lmder_functor functor;
00175   LevenbergMarquardt<lmder_functor> lm(functor);
00176   info = lm.lmder1(x);
00177 
00178   // check return value
00179   VERIFY_IS_EQUAL(info, 1);
00180   VERIFY_IS_EQUAL(lm.nfev, 6);
00181   VERIFY_IS_EQUAL(lm.njev, 5);
00182 
00183   // check norm
00184   VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596);
00185 
00186   // check x
00187   VectorXd x_ref(n);
00188   x_ref << 0.08241058, 1.133037, 2.343695;
00189   VERIFY_IS_APPROX(x, x_ref);
00190 }
00191 
00192 void testLmder()
00193 {
00194   const int m=15, n=3;
00195   int info;
00196   double fnorm, covfac;
00197   VectorXd x;
00198 
00199   /* the following starting values provide a rough fit. */
00200   x.setConstant(n, 1.);
00201 
00202   // do the computation
00203   lmder_functor functor;
00204   LevenbergMarquardt<lmder_functor> lm(functor);
00205   info = lm.minimize(x);
00206 
00207   // check return values
00208   VERIFY_IS_EQUAL(info, 1);
00209   VERIFY_IS_EQUAL(lm.nfev, 6);
00210   VERIFY_IS_EQUAL(lm.njev, 5);
00211 
00212   // check norm
00213   fnorm = lm.fvec.blueNorm();
00214   VERIFY_IS_APPROX(fnorm, 0.09063596);
00215 
00216   // check x
00217   VectorXd x_ref(n);
00218   x_ref << 0.08241058, 1.133037, 2.343695;
00219   VERIFY_IS_APPROX(x, x_ref);
00220 
00221   // check covariance
00222   covfac = fnorm*fnorm/(m-n);
00223   internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm
00224 
00225   MatrixXd cov_ref(n,n);
00226   cov_ref <<
00227       0.0001531202,   0.002869941,  -0.002656662,
00228       0.002869941,    0.09480935,   -0.09098995,
00229       -0.002656662,   -0.09098995,    0.08778727;
00230 
00231 //  std::cout << fjac*covfac << std::endl;
00232 
00233   MatrixXd cov;
00234   cov =  covfac*lm.fjac.topLeftCorner<n,n>();
00235   VERIFY_IS_APPROX( cov, cov_ref);
00236   // TODO: why isn't this allowed ? :
00237   // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
00238 }
00239 
00240 struct hybrj_functor : Functor<double>
00241 {
00242     hybrj_functor(void) : Functor<double>(9,9) {}
00243 
00244     int operator()(const VectorXd &x, VectorXd &fvec)
00245     {
00246         double temp, temp1, temp2;
00247         const int n = x.size();
00248         assert(fvec.size()==n);
00249         for (int k = 0; k < n; k++)
00250         {
00251             temp = (3. - 2.*x[k])*x[k];
00252             temp1 = 0.;
00253             if (k) temp1 = x[k-1];
00254             temp2 = 0.;
00255             if (k != n-1) temp2 = x[k+1];
00256             fvec[k] = temp - temp1 - 2.*temp2 + 1.;
00257         }
00258         return 0;
00259     }
00260     int df(const VectorXd &x, MatrixXd &fjac)
00261     {
00262         const int n = x.size();
00263         assert(fjac.rows()==n);
00264         assert(fjac.cols()==n);
00265         for (int k = 0; k < n; k++)
00266         {
00267             for (int j = 0; j < n; j++)
00268                 fjac(k,j) = 0.;
00269             fjac(k,k) = 3.- 4.*x[k];
00270             if (k) fjac(k,k-1) = -1.;
00271             if (k != n-1) fjac(k,k+1) = -2.;
00272         }
00273         return 0;
00274     }
00275 };
00276 
00277 
00278 void testHybrj1()
00279 {
00280   const int n=9;
00281   int info;
00282   VectorXd x(n);
00283 
00284   /* the following starting values provide a rough fit. */
00285   x.setConstant(n, -1.);
00286 
00287   // do the computation
00288   hybrj_functor functor;
00289   HybridNonLinearSolver<hybrj_functor> solver(functor);
00290   info = solver.hybrj1(x);
00291 
00292   // check return value
00293   VERIFY_IS_EQUAL(info, 1);
00294   VERIFY_IS_EQUAL(solver.nfev, 11);
00295   VERIFY_IS_EQUAL(solver.njev, 1);
00296 
00297   // check norm
00298   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
00299 
00300 
00301 // check x
00302   VectorXd x_ref(n);
00303   x_ref <<
00304      -0.5706545,    -0.6816283,    -0.7017325,
00305      -0.7042129,     -0.701369,    -0.6918656,
00306      -0.665792,    -0.5960342,    -0.4164121;
00307   VERIFY_IS_APPROX(x, x_ref);
00308 }
00309 
00310 void testHybrj()
00311 {
00312   const int n=9;
00313   int info;
00314   VectorXd x(n);
00315 
00316   /* the following starting values provide a rough fit. */
00317   x.setConstant(n, -1.);
00318 
00319 
00320   // do the computation
00321   hybrj_functor functor;
00322   HybridNonLinearSolver<hybrj_functor> solver(functor);
00323   solver.diag.setConstant(n, 1.);
00324   solver.useExternalScaling = true;
00325   info = solver.solve(x);
00326 
00327   // check return value
00328   VERIFY_IS_EQUAL(info, 1);
00329   VERIFY_IS_EQUAL(solver.nfev, 11);
00330   VERIFY_IS_EQUAL(solver.njev, 1);
00331 
00332   // check norm
00333   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
00334 
00335 
00336 // check x
00337   VectorXd x_ref(n);
00338   x_ref <<
00339      -0.5706545,    -0.6816283,    -0.7017325,
00340      -0.7042129,     -0.701369,    -0.6918656,
00341      -0.665792,    -0.5960342,    -0.4164121;
00342   VERIFY_IS_APPROX(x, x_ref);
00343 
00344 }
00345 
00346 struct hybrd_functor : Functor<double>
00347 {
00348     hybrd_functor(void) : Functor<double>(9,9) {}
00349     int operator()(const VectorXd &x, VectorXd &fvec) const
00350     {
00351         double temp, temp1, temp2;
00352         const int n = x.size();
00353 
00354         assert(fvec.size()==n);
00355         for (int k=0; k < n; k++)
00356         {
00357             temp = (3. - 2.*x[k])*x[k];
00358             temp1 = 0.;
00359             if (k) temp1 = x[k-1];
00360             temp2 = 0.;
00361             if (k != n-1) temp2 = x[k+1];
00362             fvec[k] = temp - temp1 - 2.*temp2 + 1.;
00363         }
00364         return 0;
00365     }
00366 };
00367 
00368 void testHybrd1()
00369 {
00370   int n=9, info;
00371   VectorXd x(n);
00372 
00373   /* the following starting values provide a rough solution. */
00374   x.setConstant(n, -1.);
00375 
00376   // do the computation
00377   hybrd_functor functor;
00378   HybridNonLinearSolver<hybrd_functor> solver(functor);
00379   info = solver.hybrd1(x);
00380 
00381   // check return value
00382   VERIFY_IS_EQUAL(info, 1);
00383   VERIFY_IS_EQUAL(solver.nfev, 20);
00384 
00385   // check norm
00386   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
00387 
00388   // check x
00389   VectorXd x_ref(n);
00390   x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121;
00391   VERIFY_IS_APPROX(x, x_ref);
00392 }
00393 
00394 void testHybrd()
00395 {
00396   const int n=9;
00397   int info;
00398   VectorXd x;
00399 
00400   /* the following starting values provide a rough fit. */
00401   x.setConstant(n, -1.);
00402 
00403   // do the computation
00404   hybrd_functor functor;
00405   HybridNonLinearSolver<hybrd_functor> solver(functor);
00406   solver.parameters.nb_of_subdiagonals = 1;
00407   solver.parameters.nb_of_superdiagonals = 1;
00408   solver.diag.setConstant(n, 1.);
00409   solver.useExternalScaling = true;
00410   info = solver.solveNumericalDiff(x);
00411 
00412   // check return value
00413   VERIFY_IS_EQUAL(info, 1);
00414   VERIFY_IS_EQUAL(solver.nfev, 14);
00415 
00416   // check norm
00417   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
00418 
00419   // check x
00420   VectorXd x_ref(n);
00421   x_ref <<
00422       -0.5706545,    -0.6816283,    -0.7017325,
00423       -0.7042129,     -0.701369,    -0.6918656,
00424       -0.665792,    -0.5960342,    -0.4164121;
00425   VERIFY_IS_APPROX(x, x_ref);
00426 }
00427 
00428 struct lmstr_functor : Functor<double>
00429 {
00430     lmstr_functor(void) : Functor<double>(3,15) {}
00431     int operator()(const VectorXd &x, VectorXd &fvec)
00432     {
00433         /*  subroutine fcn for lmstr1 example. */
00434         double tmp1, tmp2, tmp3;
00435         static const double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
00436             3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
00437 
00438         assert(15==fvec.size());
00439         assert(3==x.size());
00440 
00441         for (int i=0; i<15; i++)
00442         {
00443             tmp1 = i+1;
00444             tmp2 = 16 - i - 1;
00445             tmp3 = (i>=8)? tmp2 : tmp1;
00446             fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
00447         }
00448         return 0;
00449     }
00450     int df(const VectorXd &x, VectorXd &jac_row, VectorXd::Index rownb)
00451     {
00452         assert(x.size()==3);
00453         assert(jac_row.size()==x.size());
00454         double tmp1, tmp2, tmp3, tmp4;
00455 
00456         int i = rownb-2;
00457         tmp1 = i+1;
00458         tmp2 = 16 - i - 1;
00459         tmp3 = (i>=8)? tmp2 : tmp1;
00460         tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
00461         jac_row[0] = -1;
00462         jac_row[1] = tmp1*tmp2/tmp4;
00463         jac_row[2] = tmp1*tmp3/tmp4;
00464         return 0;
00465     }
00466 };
00467 
00468 void testLmstr1()
00469 {
00470   const int n=3;
00471   int info;
00472 
00473   VectorXd x(n);
00474 
00475   /* the following starting values provide a rough fit. */
00476   x.setConstant(n, 1.);
00477 
00478   // do the computation
00479   lmstr_functor functor;
00480   LevenbergMarquardt<lmstr_functor> lm(functor);
00481   info = lm.lmstr1(x);
00482 
00483   // check return value
00484   VERIFY_IS_EQUAL(info, 1);
00485   VERIFY_IS_EQUAL(lm.nfev, 6);
00486   VERIFY_IS_EQUAL(lm.njev, 5);
00487 
00488   // check norm
00489   VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596);
00490 
00491   // check x
00492   VectorXd x_ref(n);
00493   x_ref << 0.08241058, 1.133037, 2.343695 ;
00494   VERIFY_IS_APPROX(x, x_ref);
00495 }
00496 
00497 void testLmstr()
00498 {
00499   const int n=3;
00500   int info;
00501   double fnorm;
00502   VectorXd x(n);
00503 
00504   /* the following starting values provide a rough fit. */
00505   x.setConstant(n, 1.);
00506 
00507   // do the computation
00508   lmstr_functor functor;
00509   LevenbergMarquardt<lmstr_functor> lm(functor);
00510   info = lm.minimizeOptimumStorage(x);
00511 
00512   // check return values
00513   VERIFY_IS_EQUAL(info, 1);
00514   VERIFY_IS_EQUAL(lm.nfev, 6);
00515   VERIFY_IS_EQUAL(lm.njev, 5);
00516 
00517   // check norm
00518   fnorm = lm.fvec.blueNorm();
00519   VERIFY_IS_APPROX(fnorm, 0.09063596);
00520 
00521   // check x
00522   VectorXd x_ref(n);
00523   x_ref << 0.08241058, 1.133037, 2.343695;
00524   VERIFY_IS_APPROX(x, x_ref);
00525 
00526 }
00527 
00528 struct lmdif_functor : Functor<double>
00529 {
00530     lmdif_functor(void) : Functor<double>(3,15) {}
00531     int operator()(const VectorXd &x, VectorXd &fvec) const
00532     {
00533         int i;
00534         double tmp1,tmp2,tmp3;
00535         static const double y[15]={1.4e-1,1.8e-1,2.2e-1,2.5e-1,2.9e-1,3.2e-1,3.5e-1,3.9e-1,
00536             3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0};
00537 
00538         assert(x.size()==3);
00539         assert(fvec.size()==15);
00540         for (i=0; i<15; i++)
00541         {
00542             tmp1 = i+1;
00543             tmp2 = 15 - i;
00544             tmp3 = tmp1;
00545 
00546             if (i >= 8) tmp3 = tmp2;
00547             fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
00548         }
00549         return 0;
00550     }
00551 };
00552 
00553 void testLmdif1()
00554 {
00555   const int n=3;
00556   int info;
00557 
00558   VectorXd x(n), fvec(15);
00559 
00560   /* the following starting values provide a rough fit. */
00561   x.setConstant(n, 1.);
00562 
00563   // do the computation
00564   lmdif_functor functor;
00565   DenseIndex nfev;
00566   info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev);
00567 
00568   // check return value
00569   VERIFY_IS_EQUAL(info, 1);
00570   VERIFY_IS_EQUAL(nfev, 26);
00571 
00572   // check norm
00573   functor(x, fvec);
00574   VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596);
00575 
00576   // check x
00577   VectorXd x_ref(n);
00578   x_ref << 0.0824106, 1.1330366, 2.3436947;
00579   VERIFY_IS_APPROX(x, x_ref);
00580 
00581 }
00582 
00583 void testLmdif()
00584 {
00585   const int m=15, n=3;
00586   int info;
00587   double fnorm, covfac;
00588   VectorXd x(n);
00589 
00590   /* the following starting values provide a rough fit. */
00591   x.setConstant(n, 1.);
00592 
00593   // do the computation
00594   lmdif_functor functor;
00595   NumericalDiff<lmdif_functor> numDiff(functor);
00596   LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff);
00597   info = lm.minimize(x);
00598 
00599   // check return values
00600   VERIFY_IS_EQUAL(info, 1);
00601   VERIFY_IS_EQUAL(lm.nfev, 26);
00602 
00603   // check norm
00604   fnorm = lm.fvec.blueNorm();
00605   VERIFY_IS_APPROX(fnorm, 0.09063596);
00606 
00607   // check x
00608   VectorXd x_ref(n);
00609   x_ref << 0.08241058, 1.133037, 2.343695;
00610   VERIFY_IS_APPROX(x, x_ref);
00611 
00612   // check covariance
00613   covfac = fnorm*fnorm/(m-n);
00614   internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm
00615 
00616   MatrixXd cov_ref(n,n);
00617   cov_ref <<
00618       0.0001531202,   0.002869942,  -0.002656662,
00619       0.002869942,    0.09480937,   -0.09098997,
00620       -0.002656662,   -0.09098997,    0.08778729;
00621 
00622 //  std::cout << fjac*covfac << std::endl;
00623 
00624   MatrixXd cov;
00625   cov =  covfac*lm.fjac.topLeftCorner<n,n>();
00626   VERIFY_IS_APPROX( cov, cov_ref);
00627   // TODO: why isn't this allowed ? :
00628   // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
00629 }
00630 
00631 struct chwirut2_functor : Functor<double>
00632 {
00633     chwirut2_functor(void) : Functor<double>(3,54) {}
00634     static const double m_x[54];
00635     static const double m_y[54];
00636     int operator()(const VectorXd &b, VectorXd &fvec)
00637     {
00638         int i;
00639 
00640         assert(b.size()==3);
00641         assert(fvec.size()==54);
00642         for(i=0; i<54; i++) {
00643             double x = m_x[i];
00644             fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i];
00645         }
00646         return 0;
00647     }
00648     int df(const VectorXd &b, MatrixXd &fjac)
00649     {
00650         assert(b.size()==3);
00651         assert(fjac.rows()==54);
00652         assert(fjac.cols()==3);
00653         for(int i=0; i<54; i++) {
00654             double x = m_x[i];
00655             double factor = 1./(b[1]+b[2]*x);
00656             double e = exp(-b[0]*x);
00657             fjac(i,0) = -x*e*factor;
00658             fjac(i,1) = -e*factor*factor;
00659             fjac(i,2) = -x*e*factor*factor;
00660         }
00661         return 0;
00662     }
00663 };
00664 const double chwirut2_functor::m_x[54] = { 0.500E0, 1.000E0, 1.750E0, 3.750E0, 5.750E0, 0.875E0, 2.250E0, 3.250E0, 5.250E0, 0.750E0, 1.750E0, 2.750E0, 4.750E0, 0.625E0, 1.250E0, 2.250E0, 4.250E0, .500E0, 3.000E0, .750E0, 3.000E0, 1.500E0, 6.000E0, 3.000E0, 6.000E0, 1.500E0, 3.000E0, .500E0, 2.000E0, 4.000E0, .750E0, 2.000E0, 5.000E0, .750E0, 2.250E0, 3.750E0, 5.750E0, 3.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .500E0, 6.000E0, 3.000E0, .500E0, 2.750E0, .500E0, 1.750E0};
00665 const double chwirut2_functor::m_y[54] = { 92.9000E0 ,57.1000E0 ,31.0500E0 ,11.5875E0 ,8.0250E0 ,63.6000E0 ,21.4000E0 ,14.2500E0 ,8.4750E0 ,63.8000E0 ,26.8000E0 ,16.4625E0 ,7.1250E0 ,67.3000E0 ,41.0000E0 ,21.1500E0 ,8.1750E0 ,81.5000E0 ,13.1200E0 ,59.9000E0 ,14.6200E0 ,32.9000E0 ,5.4400E0 ,12.5600E0 ,5.4400E0 ,32.0000E0 ,13.9500E0 ,75.8000E0 ,20.0000E0 ,10.4200E0 ,59.5000E0 ,21.6700E0 ,8.5500E0 ,62.0000E0 ,20.2000E0 ,7.7600E0 ,3.7500E0 ,11.8100E0 ,54.7000E0 ,23.7000E0 ,11.5500E0 ,61.3000E0 ,17.7000E0 ,8.7400E0 ,59.2000E0 ,16.3000E0 ,8.6200E0 ,81.0000E0 ,4.8700E0 ,14.6200E0 ,81.7000E0 ,17.1700E0 ,81.3000E0 ,28.9000E0  };
00666 
00667 // http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml
00668 void testNistChwirut2(void)
00669 {
00670   const int n=3;
00671   int info;
00672 
00673   VectorXd x(n);
00674 
00675   /*
00676    * First try
00677    */
00678   x<< 0.1, 0.01, 0.02;
00679   // do the computation
00680   chwirut2_functor functor;
00681   LevenbergMarquardt<chwirut2_functor> lm(functor);
00682   info = lm.minimize(x);
00683 
00684   // check return value
00685   VERIFY_IS_EQUAL(info, 1);
00686   VERIFY_IS_EQUAL(lm.nfev, 10);
00687   VERIFY_IS_EQUAL(lm.njev, 8);
00688   // check norm^2
00689   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02);
00690   // check x
00691   VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
00692   VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
00693   VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
00694 
00695   /*
00696    * Second try
00697    */
00698   x<< 0.15, 0.008, 0.010;
00699   // do the computation
00700   lm.resetParameters();
00701   lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
00702   lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
00703   info = lm.minimize(x);
00704 
00705   // check return value
00706   VERIFY_IS_EQUAL(info, 1);
00707   VERIFY_IS_EQUAL(lm.nfev, 7);
00708   VERIFY_IS_EQUAL(lm.njev, 6);
00709   // check norm^2
00710   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02);
00711   // check x
00712   VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
00713   VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
00714   VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
00715 }
00716 
00717 
00718 struct misra1a_functor : Functor<double>
00719 {
00720     misra1a_functor(void) : Functor<double>(2,14) {}
00721     static const double m_x[14];
00722     static const double m_y[14];
00723     int operator()(const VectorXd &b, VectorXd &fvec)
00724     {
00725         assert(b.size()==2);
00726         assert(fvec.size()==14);
00727         for(int i=0; i<14; i++) {
00728             fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ;
00729         }
00730         return 0;
00731     }
00732     int df(const VectorXd &b, MatrixXd &fjac)
00733     {
00734         assert(b.size()==2);
00735         assert(fjac.rows()==14);
00736         assert(fjac.cols()==2);
00737         for(int i=0; i<14; i++) {
00738             fjac(i,0) = (1.-exp(-b[1]*m_x[i]));
00739             fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i]));
00740         }
00741         return 0;
00742     }
00743 };
00744 const double misra1a_functor::m_x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0};
00745 const double misra1a_functor::m_y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0};
00746 
00747 // http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
00748 void testNistMisra1a(void)
00749 {
00750   const int n=2;
00751   int info;
00752 
00753   VectorXd x(n);
00754 
00755   /*
00756    * First try
00757    */
00758   x<< 500., 0.0001;
00759   // do the computation
00760   misra1a_functor functor;
00761   LevenbergMarquardt<misra1a_functor> lm(functor);
00762   info = lm.minimize(x);
00763 
00764   // check return value
00765   VERIFY_IS_EQUAL(info, 1);
00766   VERIFY_IS_EQUAL(lm.nfev, 19);
00767   VERIFY_IS_EQUAL(lm.njev, 15);
00768   // check norm^2
00769   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01);
00770   // check x
00771   VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
00772   VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
00773 
00774   /*
00775    * Second try
00776    */
00777   x<< 250., 0.0005;
00778   // do the computation
00779   info = lm.minimize(x);
00780 
00781   // check return value
00782   VERIFY_IS_EQUAL(info, 1);
00783   VERIFY_IS_EQUAL(lm.nfev, 5);
00784   VERIFY_IS_EQUAL(lm.njev, 4);
00785   // check norm^2
00786   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01);
00787   // check x
00788   VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
00789   VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
00790 }
00791 
00792 struct hahn1_functor : Functor<double>
00793 {
00794     hahn1_functor(void) : Functor<double>(7,236) {}
00795     static const double m_x[236];
00796     int operator()(const VectorXd &b, VectorXd &fvec)
00797     {
00798         static const double m_y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0  , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0  , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0  , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 , 16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0  , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0   , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 12.596E0 , 13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0  , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0  , 20.935E0 , 21.035E0 , 20.93E0  , 21.074E0 , 21.085E0 , 20.935E0 };
00799 
00800         //        int called=0; printf("call hahn1_functor with  iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
00801 
00802         assert(b.size()==7);
00803         assert(fvec.size()==236);
00804         for(int i=0; i<236; i++) {
00805             double x=m_x[i], xx=x*x, xxx=xx*x;
00806             fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - m_y[i];
00807         }
00808         return 0;
00809     }
00810 
00811     int df(const VectorXd &b, MatrixXd &fjac)
00812     {
00813         assert(b.size()==7);
00814         assert(fjac.rows()==236);
00815         assert(fjac.cols()==7);
00816         for(int i=0; i<236; i++) {
00817             double x=m_x[i], xx=x*x, xxx=xx*x;
00818             double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
00819             fjac(i,0) = 1.*fact;
00820             fjac(i,1) = x*fact;
00821             fjac(i,2) = xx*fact;
00822             fjac(i,3) = xxx*fact;
00823             fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
00824             fjac(i,4) = x*fact;
00825             fjac(i,5) = xx*fact;
00826             fjac(i,6) = xxx*fact;
00827         }
00828         return 0;
00829     }
00830 };
00831 const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0  , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , 282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 , 141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0  , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0  , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0};
00832 
00833 // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml
00834 void testNistHahn1(void)
00835 {
00836   const int  n=7;
00837   int info;
00838 
00839   VectorXd x(n);
00840 
00841   /*
00842    * First try
00843    */
00844   x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
00845   // do the computation
00846   hahn1_functor functor;
00847   LevenbergMarquardt<hahn1_functor> lm(functor);
00848   info = lm.minimize(x);
00849 
00850   // check return value
00851   VERIFY_IS_EQUAL(info, 1);
00852   VERIFY_IS_EQUAL(lm.nfev, 11);
00853   VERIFY_IS_EQUAL(lm.njev, 10);
00854   // check norm^2
00855   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00);
00856   // check x
00857   VERIFY_IS_APPROX(x[0], 1.0776351733E+00);
00858   VERIFY_IS_APPROX(x[1],-1.2269296921E-01);
00859   VERIFY_IS_APPROX(x[2], 4.0863750610E-03);
00860   VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06
00861   VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
00862   VERIFY_IS_APPROX(x[5], 2.4053735503E-04);
00863   VERIFY_IS_APPROX(x[6],-1.2314450199E-07);
00864 
00865   /*
00866    * Second try
00867    */
00868   x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001;
00869   // do the computation
00870   info = lm.minimize(x);
00871 
00872   // check return value
00873   VERIFY_IS_EQUAL(info, 1);
00874   VERIFY_IS_EQUAL(lm.nfev, 11);
00875   VERIFY_IS_EQUAL(lm.njev, 10);
00876   // check norm^2
00877   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00);
00878   // check x
00879   VERIFY_IS_APPROX(x[0], 1.077640); // should be :  1.0776351733E+00
00880   VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01
00881   VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03
00882   VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06
00883   VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
00884   VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04
00885   VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07
00886 
00887 }
00888 
00889 struct misra1d_functor : Functor<double>
00890 {
00891     misra1d_functor(void) : Functor<double>(2,14) {}
00892     static const double x[14];
00893     static const double y[14];
00894     int operator()(const VectorXd &b, VectorXd &fvec)
00895     {
00896         assert(b.size()==2);
00897         assert(fvec.size()==14);
00898         for(int i=0; i<14; i++) {
00899             fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i];
00900         }
00901         return 0;
00902     }
00903     int df(const VectorXd &b, MatrixXd &fjac)
00904     {
00905         assert(b.size()==2);
00906         assert(fjac.rows()==14);
00907         assert(fjac.cols()==2);
00908         for(int i=0; i<14; i++) {
00909             double den = 1.+b[1]*x[i];
00910             fjac(i,0) = b[1]*x[i] / den;
00911             fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den;
00912         }
00913         return 0;
00914     }
00915 };
00916 const double misra1d_functor::x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0};
00917 const double misra1d_functor::y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0};
00918 
00919 // http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml
00920 void testNistMisra1d(void)
00921 {
00922   const int n=2;
00923   int info;
00924 
00925   VectorXd x(n);
00926 
00927   /*
00928    * First try
00929    */
00930   x<< 500., 0.0001;
00931   // do the computation
00932   misra1d_functor functor;
00933   LevenbergMarquardt<misra1d_functor> lm(functor);
00934   info = lm.minimize(x);
00935 
00936   // check return value
00937   VERIFY_IS_EQUAL(info, 3);
00938   VERIFY_IS_EQUAL(lm.nfev, 9);
00939   VERIFY_IS_EQUAL(lm.njev, 7);
00940   // check norm^2
00941   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02);
00942   // check x
00943   VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
00944   VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
00945 
00946   /*
00947    * Second try
00948    */
00949   x<< 450., 0.0003;
00950   // do the computation
00951   info = lm.minimize(x);
00952 
00953   // check return value
00954   VERIFY_IS_EQUAL(info, 1);
00955   VERIFY_IS_EQUAL(lm.nfev, 4);
00956   VERIFY_IS_EQUAL(lm.njev, 3);
00957   // check norm^2
00958   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02);
00959   // check x
00960   VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
00961   VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
00962 }
00963 
00964 
00965 struct lanczos1_functor : Functor<double>
00966 {
00967     lanczos1_functor(void) : Functor<double>(6,24) {}
00968     static const double x[24];
00969     static const double y[24];
00970     int operator()(const VectorXd &b, VectorXd &fvec)
00971     {
00972         assert(b.size()==6);
00973         assert(fvec.size()==24);
00974         for(int i=0; i<24; i++)
00975             fvec[i] = b[0]*exp(-b[1]*x[i]) + b[2]*exp(-b[3]*x[i]) + b[4]*exp(-b[5]*x[i])  - y[i];
00976         return 0;
00977     }
00978     int df(const VectorXd &b, MatrixXd &fjac)
00979     {
00980         assert(b.size()==6);
00981         assert(fjac.rows()==24);
00982         assert(fjac.cols()==6);
00983         for(int i=0; i<24; i++) {
00984             fjac(i,0) = exp(-b[1]*x[i]);
00985             fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]);
00986             fjac(i,2) = exp(-b[3]*x[i]);
00987             fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]);
00988             fjac(i,4) = exp(-b[5]*x[i]);
00989             fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]);
00990         }
00991         return 0;
00992     }
00993 };
00994 const double lanczos1_functor::x[24] = { 0.000000000000E+00, 5.000000000000E-02, 1.000000000000E-01, 1.500000000000E-01, 2.000000000000E-01, 2.500000000000E-01, 3.000000000000E-01, 3.500000000000E-01, 4.000000000000E-01, 4.500000000000E-01, 5.000000000000E-01, 5.500000000000E-01, 6.000000000000E-01, 6.500000000000E-01, 7.000000000000E-01, 7.500000000000E-01, 8.000000000000E-01, 8.500000000000E-01, 9.000000000000E-01, 9.500000000000E-01, 1.000000000000E+00, 1.050000000000E+00, 1.100000000000E+00, 1.150000000000E+00 };
00995 const double lanczos1_functor::y[24] = { 2.513400000000E+00 ,2.044333373291E+00 ,1.668404436564E+00 ,1.366418021208E+00 ,1.123232487372E+00 ,9.268897180037E-01 ,7.679338563728E-01 ,6.388775523106E-01 ,5.337835317402E-01 ,4.479363617347E-01 ,3.775847884350E-01 ,3.197393199326E-01 ,2.720130773746E-01 ,2.324965529032E-01 ,1.996589546065E-01 ,1.722704126914E-01 ,1.493405660168E-01 ,1.300700206922E-01 ,1.138119324644E-01 ,1.000415587559E-01 ,8.833209084540E-02 ,7.833544019350E-02 ,6.976693743449E-02 ,6.239312536719E-02 };
00996 
00997 // http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml
00998 void testNistLanczos1(void)
00999 {
01000   const int n=6;
01001   int info;
01002 
01003   VectorXd x(n);
01004 
01005   /*
01006    * First try
01007    */
01008   x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6;
01009   // do the computation
01010   lanczos1_functor functor;
01011   LevenbergMarquardt<lanczos1_functor> lm(functor);
01012   info = lm.minimize(x);
01013 
01014   // check return value
01015   VERIFY_IS_EQUAL(info, 2);
01016   VERIFY_IS_EQUAL(lm.nfev, 79);
01017   VERIFY_IS_EQUAL(lm.njev, 72);
01018   // check norm^2
01019   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430899764097e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
01020   // check x
01021   VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
01022   VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
01023   VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
01024   VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
01025   VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
01026   VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
01027 
01028   /*
01029    * Second try
01030    */
01031   x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3;
01032   // do the computation
01033   info = lm.minimize(x);
01034 
01035   // check return value
01036   VERIFY_IS_EQUAL(info, 2);
01037   VERIFY_IS_EQUAL(lm.nfev, 9);
01038   VERIFY_IS_EQUAL(lm.njev, 8);
01039   // check norm^2
01040   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.428595533845e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
01041   // check x
01042   VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
01043   VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
01044   VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
01045   VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
01046   VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
01047   VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
01048 
01049 }
01050 
01051 struct rat42_functor : Functor<double>
01052 {
01053     rat42_functor(void) : Functor<double>(3,9) {}
01054     static const double x[9];
01055     static const double y[9];
01056     int operator()(const VectorXd &b, VectorXd &fvec)
01057     {
01058         assert(b.size()==3);
01059         assert(fvec.size()==9);
01060         for(int i=0; i<9; i++) {
01061             fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i];
01062         }
01063         return 0;
01064     }
01065 
01066     int df(const VectorXd &b, MatrixXd &fjac)
01067     {
01068         assert(b.size()==3);
01069         assert(fjac.rows()==9);
01070         assert(fjac.cols()==3);
01071         for(int i=0; i<9; i++) {
01072             double e = exp(b[1]-b[2]*x[i]);
01073             fjac(i,0) = 1./(1.+e);
01074             fjac(i,1) = -b[0]*e/(1.+e)/(1.+e);
01075             fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e);
01076         }
01077         return 0;
01078     }
01079 };
01080 const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 };
01081 const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 };
01082 
01083 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml
01084 void testNistRat42(void)
01085 {
01086   const int n=3;
01087   int info;
01088 
01089   VectorXd x(n);
01090 
01091   /*
01092    * First try
01093    */
01094   x<< 100., 1., 0.1;
01095   // do the computation
01096   rat42_functor functor;
01097   LevenbergMarquardt<rat42_functor> lm(functor);
01098   info = lm.minimize(x);
01099 
01100   // check return value
01101   VERIFY_IS_EQUAL(info, 1);
01102   VERIFY_IS_EQUAL(lm.nfev, 10);
01103   VERIFY_IS_EQUAL(lm.njev, 8);
01104   // check norm^2
01105   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00);
01106   // check x
01107   VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
01108   VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
01109   VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
01110 
01111   /*
01112    * Second try
01113    */
01114   x<< 75., 2.5, 0.07;
01115   // do the computation
01116   info = lm.minimize(x);
01117 
01118   // check return value
01119   VERIFY_IS_EQUAL(info, 1);
01120   VERIFY_IS_EQUAL(lm.nfev, 6);
01121   VERIFY_IS_EQUAL(lm.njev, 5);
01122   // check norm^2
01123   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00);
01124   // check x
01125   VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
01126   VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
01127   VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
01128 }
01129 
01130 struct MGH10_functor : Functor<double>
01131 {
01132     MGH10_functor(void) : Functor<double>(3,16) {}
01133     static const double x[16];
01134     static const double y[16];
01135     int operator()(const VectorXd &b, VectorXd &fvec)
01136     {
01137         assert(b.size()==3);
01138         assert(fvec.size()==16);
01139         for(int i=0; i<16; i++)
01140             fvec[i] =  b[0] * exp(b[1]/(x[i]+b[2])) - y[i];
01141         return 0;
01142     }
01143     int df(const VectorXd &b, MatrixXd &fjac)
01144     {
01145         assert(b.size()==3);
01146         assert(fjac.rows()==16);
01147         assert(fjac.cols()==3);
01148         for(int i=0; i<16; i++) {
01149             double factor = 1./(x[i]+b[2]);
01150             double e = exp(b[1]*factor);
01151             fjac(i,0) = e;
01152             fjac(i,1) = b[0]*factor*e;
01153             fjac(i,2) = -b[1]*b[0]*factor*factor*e;
01154         }
01155         return 0;
01156     }
01157 };
01158 const double MGH10_functor::x[16] = { 5.000000E+01, 5.500000E+01, 6.000000E+01, 6.500000E+01, 7.000000E+01, 7.500000E+01, 8.000000E+01, 8.500000E+01, 9.000000E+01, 9.500000E+01, 1.000000E+02, 1.050000E+02, 1.100000E+02, 1.150000E+02, 1.200000E+02, 1.250000E+02 };
01159 const double MGH10_functor::y[16] = { 3.478000E+04, 2.861000E+04, 2.365000E+04, 1.963000E+04, 1.637000E+04, 1.372000E+04, 1.154000E+04, 9.744000E+03, 8.261000E+03, 7.030000E+03, 6.005000E+03, 5.147000E+03, 4.427000E+03, 3.820000E+03, 3.307000E+03, 2.872000E+03 };
01160 
01161 // http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
01162 void testNistMGH10(void)
01163 {
01164   const int n=3;
01165   int info;
01166 
01167   VectorXd x(n);
01168 
01169   /*
01170    * First try
01171    */
01172   x<< 2., 400000., 25000.;
01173   // do the computation
01174   MGH10_functor functor;
01175   LevenbergMarquardt<MGH10_functor> lm(functor);
01176   info = lm.minimize(x);
01177 
01178   // check return value
01179   VERIFY_IS_EQUAL(info, 2); 
01180   VERIFY_IS_EQUAL(lm.nfev, 284 ); 
01181   VERIFY_IS_EQUAL(lm.njev, 249 ); 
01182   // check norm^2
01183   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01);
01184   // check x
01185   VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
01186   VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
01187   VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
01188 
01189   /*
01190    * Second try
01191    */
01192   x<< 0.02, 4000., 250.;
01193   // do the computation
01194   info = lm.minimize(x);
01195 
01196   // check return value
01197   VERIFY_IS_EQUAL(info, 3);
01198   VERIFY_IS_EQUAL(lm.nfev, 126);
01199   VERIFY_IS_EQUAL(lm.njev, 116);
01200   // check norm^2
01201   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01);
01202   // check x
01203   VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
01204   VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
01205   VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
01206 }
01207 
01208 
01209 struct BoxBOD_functor : Functor<double>
01210 {
01211     BoxBOD_functor(void) : Functor<double>(2,6) {}
01212     static const double x[6];
01213     int operator()(const VectorXd &b, VectorXd &fvec)
01214     {
01215         static const double y[6] = { 109., 149., 149., 191., 213., 224. };
01216         assert(b.size()==2);
01217         assert(fvec.size()==6);
01218         for(int i=0; i<6; i++)
01219             fvec[i] =  b[0]*(1.-exp(-b[1]*x[i])) - y[i];
01220         return 0;
01221     }
01222     int df(const VectorXd &b, MatrixXd &fjac)
01223     {
01224         assert(b.size()==2);
01225         assert(fjac.rows()==6);
01226         assert(fjac.cols()==2);
01227         for(int i=0; i<6; i++) {
01228             double e = exp(-b[1]*x[i]);
01229             fjac(i,0) = 1.-e;
01230             fjac(i,1) = b[0]*x[i]*e;
01231         }
01232         return 0;
01233     }
01234 };
01235 const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. };
01236 
01237 // http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml
01238 void testNistBoxBOD(void)
01239 {
01240   const int n=2;
01241   int info;
01242 
01243   VectorXd x(n);
01244 
01245   /*
01246    * First try
01247    */
01248   x<< 1., 1.;
01249   // do the computation
01250   BoxBOD_functor functor;
01251   LevenbergMarquardt<BoxBOD_functor> lm(functor);
01252   lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
01253   lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
01254   lm.parameters.factor = 10.;
01255   info = lm.minimize(x);
01256 
01257   // check return value
01258   VERIFY_IS_EQUAL(info, 1);
01259   VERIFY_IS_EQUAL(lm.nfev, 31);
01260   VERIFY_IS_EQUAL(lm.njev, 25);
01261   // check norm^2
01262   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
01263   // check x
01264   VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
01265   VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
01266 
01267   /*
01268    * Second try
01269    */
01270   x<< 100., 0.75;
01271   // do the computation
01272   lm.resetParameters();
01273   lm.parameters.ftol = NumTraits<double>::epsilon();
01274   lm.parameters.xtol = NumTraits<double>::epsilon();
01275   info = lm.minimize(x);
01276 
01277   // check return value
01278   VERIFY_IS_EQUAL(info, 1); 
01279   VERIFY_IS_EQUAL(lm.nfev, 15 ); 
01280   VERIFY_IS_EQUAL(lm.njev, 14 ); 
01281   // check norm^2
01282   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
01283   // check x
01284   VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
01285   VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
01286 }
01287 
01288 struct MGH17_functor : Functor<double>
01289 {
01290     MGH17_functor(void) : Functor<double>(5,33) {}
01291     static const double x[33];
01292     static const double y[33];
01293     int operator()(const VectorXd &b, VectorXd &fvec)
01294     {
01295         assert(b.size()==5);
01296         assert(fvec.size()==33);
01297         for(int i=0; i<33; i++)
01298             fvec[i] =  b[0] + b[1]*exp(-b[3]*x[i]) +  b[2]*exp(-b[4]*x[i]) - y[i];
01299         return 0;
01300     }
01301     int df(const VectorXd &b, MatrixXd &fjac)
01302     {
01303         assert(b.size()==5);
01304         assert(fjac.rows()==33);
01305         assert(fjac.cols()==5);
01306         for(int i=0; i<33; i++) {
01307             fjac(i,0) = 1.;
01308             fjac(i,1) = exp(-b[3]*x[i]);
01309             fjac(i,2) = exp(-b[4]*x[i]);
01310             fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]);
01311             fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]);
01312         }
01313         return 0;
01314     }
01315 };
01316 const double MGH17_functor::x[33] = { 0.000000E+00, 1.000000E+01, 2.000000E+01, 3.000000E+01, 4.000000E+01, 5.000000E+01, 6.000000E+01, 7.000000E+01, 8.000000E+01, 9.000000E+01, 1.000000E+02, 1.100000E+02, 1.200000E+02, 1.300000E+02, 1.400000E+02, 1.500000E+02, 1.600000E+02, 1.700000E+02, 1.800000E+02, 1.900000E+02, 2.000000E+02, 2.100000E+02, 2.200000E+02, 2.300000E+02, 2.400000E+02, 2.500000E+02, 2.600000E+02, 2.700000E+02, 2.800000E+02, 2.900000E+02, 3.000000E+02, 3.100000E+02, 3.200000E+02 };
01317 const double MGH17_functor::y[33] = { 8.440000E-01, 9.080000E-01, 9.320000E-01, 9.360000E-01, 9.250000E-01, 9.080000E-01, 8.810000E-01, 8.500000E-01, 8.180000E-01, 7.840000E-01, 7.510000E-01, 7.180000E-01, 6.850000E-01, 6.580000E-01, 6.280000E-01, 6.030000E-01, 5.800000E-01, 5.580000E-01, 5.380000E-01, 5.220000E-01, 5.060000E-01, 4.900000E-01, 4.780000E-01, 4.670000E-01, 4.570000E-01, 4.480000E-01, 4.380000E-01, 4.310000E-01, 4.240000E-01, 4.200000E-01, 4.140000E-01, 4.110000E-01, 4.060000E-01 };
01318 
01319 // http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml
01320 void testNistMGH17(void)
01321 {
01322   const int n=5;
01323   int info;
01324 
01325   VectorXd x(n);
01326 
01327   /*
01328    * First try
01329    */
01330   x<< 50., 150., -100., 1., 2.;
01331   // do the computation
01332   MGH17_functor functor;
01333   LevenbergMarquardt<MGH17_functor> lm(functor);
01334   lm.parameters.ftol = NumTraits<double>::epsilon();
01335   lm.parameters.xtol = NumTraits<double>::epsilon();
01336   lm.parameters.maxfev = 1000;
01337   info = lm.minimize(x);
01338 
01339   // check return value
01340   VERIFY_IS_EQUAL(info, 2); 
01341   VERIFY_IS_EQUAL(lm.nfev, 602 ); 
01342   VERIFY_IS_EQUAL(lm.njev, 545 ); 
01343   // check norm^2
01344   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
01345   // check x
01346   VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
01347   VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
01348   VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
01349   VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
01350   VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
01351 
01352   /*
01353    * Second try
01354    */
01355   x<< 0.5  ,1.5  ,-1   ,0.01 ,0.02;
01356   // do the computation
01357   lm.resetParameters();
01358   info = lm.minimize(x);
01359 
01360   // check return value
01361   VERIFY_IS_EQUAL(info, 1);
01362   VERIFY_IS_EQUAL(lm.nfev, 18);
01363   VERIFY_IS_EQUAL(lm.njev, 15);
01364   // check norm^2
01365   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
01366   // check x
01367   VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
01368   VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
01369   VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
01370   VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
01371   VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
01372 }
01373 
01374 struct MGH09_functor : Functor<double>
01375 {
01376     MGH09_functor(void) : Functor<double>(4,11) {}
01377     static const double _x[11];
01378     static const double y[11];
01379     int operator()(const VectorXd &b, VectorXd &fvec)
01380     {
01381         assert(b.size()==4);
01382         assert(fvec.size()==11);
01383         for(int i=0; i<11; i++) {
01384             double x = _x[i], xx=x*x;
01385             fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i];
01386         }
01387         return 0;
01388     }
01389     int df(const VectorXd &b, MatrixXd &fjac)
01390     {
01391         assert(b.size()==4);
01392         assert(fjac.rows()==11);
01393         assert(fjac.cols()==4);
01394         for(int i=0; i<11; i++) {
01395             double x = _x[i], xx=x*x;
01396             double factor = 1./(xx+x*b[2]+b[3]);
01397             fjac(i,0) = (xx+x*b[1]) * factor;
01398             fjac(i,1) = b[0]*x* factor;
01399             fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor;
01400             fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor;
01401         }
01402         return 0;
01403     }
01404 };
01405 const double MGH09_functor::_x[11] = { 4., 2., 1., 5.E-1 , 2.5E-01, 1.670000E-01, 1.250000E-01,  1.E-01, 8.330000E-02, 7.140000E-02, 6.250000E-02 };
01406 const double MGH09_functor::y[11] = { 1.957000E-01, 1.947000E-01, 1.735000E-01, 1.600000E-01, 8.440000E-02, 6.270000E-02, 4.560000E-02, 3.420000E-02, 3.230000E-02, 2.350000E-02, 2.460000E-02 };
01407 
01408 // http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
01409 void testNistMGH09(void)
01410 {
01411   const int n=4;
01412   int info;
01413 
01414   VectorXd x(n);
01415 
01416   /*
01417    * First try
01418    */
01419   x<< 25., 39, 41.5, 39.;
01420   // do the computation
01421   MGH09_functor functor;
01422   LevenbergMarquardt<MGH09_functor> lm(functor);
01423   lm.parameters.maxfev = 1000;
01424   info = lm.minimize(x);
01425 
01426   // check return value
01427   VERIFY_IS_EQUAL(info, 1); 
01428   VERIFY_IS_EQUAL(lm.nfev, 490 ); 
01429   VERIFY_IS_EQUAL(lm.njev, 376 ); 
01430   // check norm^2
01431   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04);
01432   // check x
01433   VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01
01434   VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01
01435   VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01
01436   VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01
01437 
01438   /*
01439    * Second try
01440    */
01441   x<< 0.25, 0.39, 0.415, 0.39;
01442   // do the computation
01443   lm.resetParameters();
01444   info = lm.minimize(x);
01445 
01446   // check return value
01447   VERIFY_IS_EQUAL(info, 1);
01448   VERIFY_IS_EQUAL(lm.nfev, 18);
01449   VERIFY_IS_EQUAL(lm.njev, 16);
01450   // check norm^2
01451   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04);
01452   // check x
01453   VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01
01454   VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01
01455   VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01
01456   VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01
01457 }
01458 
01459 
01460 
01461 struct Bennett5_functor : Functor<double>
01462 {
01463     Bennett5_functor(void) : Functor<double>(3,154) {}
01464     static const double x[154];
01465     static const double y[154];
01466     int operator()(const VectorXd &b, VectorXd &fvec)
01467     {
01468         assert(b.size()==3);
01469         assert(fvec.size()==154);
01470         for(int i=0; i<154; i++)
01471             fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i];
01472         return 0;
01473     }
01474     int df(const VectorXd &b, MatrixXd &fjac)
01475     {
01476         assert(b.size()==3);
01477         assert(fjac.rows()==154);
01478         assert(fjac.cols()==3);
01479         for(int i=0; i<154; i++) {
01480             double e = pow(b[1]+x[i],-1./b[2]);
01481             fjac(i,0) = e;
01482             fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]);
01483             fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2];
01484         }
01485         return 0;
01486     }
01487 };
01488 const double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, 11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 };
01489 const double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 ,-32.238602E0 ,-32.229900E0 ,-32.224098E0 ,-32.215401E0 ,-32.203800E0 ,-32.198002E0 ,-32.189400E0 ,-32.183601E0 ,-32.174900E0 ,-32.169102E0 ,-32.163300E0 ,-32.154598E0 ,-32.145901E0 ,-32.140099E0 ,-32.131401E0 ,-32.125599E0 ,-32.119801E0 ,-32.111198E0 ,-32.105400E0 ,-32.096699E0 ,-32.090900E0 ,-32.088001E0 ,-32.079300E0 ,-32.073502E0 ,-32.067699E0 ,-32.061901E0 ,-32.056099E0 ,-32.050301E0 ,-32.044498E0 ,-32.038799E0 ,-32.033001E0 ,-32.027199E0 ,-32.024300E0 ,-32.018501E0 ,-32.012699E0 ,-32.004002E0 ,-32.001099E0 ,-31.995300E0 ,-31.989500E0 ,-31.983700E0 ,-31.977900E0 ,-31.972099E0 ,-31.969299E0 ,-31.963501E0 ,-31.957701E0 ,-31.951900E0 ,-31.946100E0 ,-31.940300E0 ,-31.937401E0 ,-31.931601E0 ,-31.925800E0 ,-31.922899E0 ,-31.917101E0 ,-31.911301E0 ,-31.908400E0 ,-31.902599E0 ,-31.896900E0 ,-31.893999E0 ,-31.888201E0 ,-31.885300E0 ,-31.882401E0 ,-31.876600E0 ,-31.873699E0 ,-31.867901E0 ,-31.862101E0 ,-31.859200E0 ,-31.856300E0 ,-31.850500E0 ,-31.844700E0 ,-31.841801E0 ,-31.838900E0 ,-31.833099E0 ,-31.830200E0 ,-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 };
01490 
01491 // http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml
01492 void testNistBennett5(void)
01493 {
01494   const int  n=3;
01495   int info;
01496 
01497   VectorXd x(n);
01498 
01499   /*
01500    * First try
01501    */
01502   x<< -2000., 50., 0.8;
01503   // do the computation
01504   Bennett5_functor functor;
01505   LevenbergMarquardt<Bennett5_functor> lm(functor);
01506   lm.parameters.maxfev = 1000;
01507   info = lm.minimize(x);
01508 
01509   // check return value
01510   VERIFY_IS_EQUAL(info, 1);
01511   VERIFY_IS_EQUAL(lm.nfev, 758);
01512   VERIFY_IS_EQUAL(lm.njev, 744);
01513   // check norm^2
01514   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04);
01515   // check x
01516   VERIFY_IS_APPROX(x[0], -2.5235058043E+03);
01517   VERIFY_IS_APPROX(x[1], 4.6736564644E+01);
01518   VERIFY_IS_APPROX(x[2], 9.3218483193E-01);
01519   /*
01520    * Second try
01521    */
01522   x<< -1500., 45., 0.85;
01523   // do the computation
01524   lm.resetParameters();
01525   info = lm.minimize(x);
01526 
01527   // check return value
01528   VERIFY_IS_EQUAL(info, 1);
01529   VERIFY_IS_EQUAL(lm.nfev, 203);
01530   VERIFY_IS_EQUAL(lm.njev, 192);
01531   // check norm^2
01532   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04);
01533   // check x
01534   VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03
01535   VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01);
01536   VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01);
01537 }
01538 
01539 struct thurber_functor : Functor<double>
01540 {
01541     thurber_functor(void) : Functor<double>(7,37) {}
01542     static const double _x[37];
01543     static const double _y[37];
01544     int operator()(const VectorXd &b, VectorXd &fvec)
01545     {
01546         //        int called=0; printf("call hahn1_functor with  iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
01547         assert(b.size()==7);
01548         assert(fvec.size()==37);
01549         for(int i=0; i<37; i++) {
01550             double x=_x[i], xx=x*x, xxx=xx*x;
01551             fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - _y[i];
01552         }
01553         return 0;
01554     }
01555     int df(const VectorXd &b, MatrixXd &fjac)
01556     {
01557         assert(b.size()==7);
01558         assert(fjac.rows()==37);
01559         assert(fjac.cols()==7);
01560         for(int i=0; i<37; i++) {
01561             double x=_x[i], xx=x*x, xxx=xx*x;
01562             double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
01563             fjac(i,0) = 1.*fact;
01564             fjac(i,1) = x*fact;
01565             fjac(i,2) = xx*fact;
01566             fjac(i,3) = xxx*fact;
01567             fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
01568             fjac(i,4) = x*fact;
01569             fjac(i,5) = xx*fact;
01570             fjac(i,6) = xxx*fact;
01571         }
01572         return 0;
01573     }
01574 };
01575 const double thurber_functor::_x[37] = { -3.067E0, -2.981E0, -2.921E0, -2.912E0, -2.840E0, -2.797E0, -2.702E0, -2.699E0, -2.633E0, -2.481E0, -2.363E0, -2.322E0, -1.501E0, -1.460E0, -1.274E0, -1.212E0, -1.100E0, -1.046E0, -0.915E0, -0.714E0, -0.566E0, -0.545E0, -0.400E0, -0.309E0, -0.109E0, -0.103E0, 0.010E0, 0.119E0, 0.377E0, 0.790E0, 0.963E0, 1.006E0, 1.115E0, 1.572E0, 1.841E0, 2.047E0, 2.200E0 };
01576 const double thurber_functor::_y[37] = { 80.574E0, 84.248E0, 87.264E0, 87.195E0, 89.076E0, 89.608E0, 89.868E0, 90.101E0, 92.405E0, 95.854E0, 100.696E0, 101.060E0, 401.672E0, 390.724E0, 567.534E0, 635.316E0, 733.054E0, 759.087E0, 894.206E0, 990.785E0, 1090.109E0, 1080.914E0, 1122.643E0, 1178.351E0, 1260.531E0, 1273.514E0, 1288.339E0, 1327.543E0, 1353.863E0, 1414.509E0, 1425.208E0, 1421.384E0, 1442.962E0, 1464.350E0, 1468.705E0, 1447.894E0, 1457.628E0};
01577 
01578 // http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml
01579 void testNistThurber(void)
01580 {
01581   const int n=7;
01582   int info;
01583 
01584   VectorXd x(n);
01585 
01586   /*
01587    * First try
01588    */
01589   x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ;
01590   // do the computation
01591   thurber_functor functor;
01592   LevenbergMarquardt<thurber_functor> lm(functor);
01593   lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon();
01594   lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon();
01595   info = lm.minimize(x);
01596 
01597   // check return value
01598   VERIFY_IS_EQUAL(info, 1);
01599   VERIFY_IS_EQUAL(lm.nfev, 39);
01600   VERIFY_IS_EQUAL(lm.njev, 36);
01601   // check norm^2
01602   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03);
01603   // check x
01604   VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
01605   VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
01606   VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
01607   VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
01608   VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
01609   VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
01610   VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
01611 
01612   /*
01613    * Second try
01614    */
01615   x<< 1300 ,1500 ,500  ,75   ,1    ,0.4  ,0.05  ;
01616   // do the computation
01617   lm.resetParameters();
01618   lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon();
01619   lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon();
01620   info = lm.minimize(x);
01621 
01622   // check return value
01623   VERIFY_IS_EQUAL(info, 1);
01624   VERIFY_IS_EQUAL(lm.nfev, 29);
01625   VERIFY_IS_EQUAL(lm.njev, 28);
01626   // check norm^2
01627   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03);
01628   // check x
01629   VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
01630   VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
01631   VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
01632   VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
01633   VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
01634   VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
01635   VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
01636 }
01637 
01638 struct rat43_functor : Functor<double>
01639 {
01640     rat43_functor(void) : Functor<double>(4,15) {}
01641     static const double x[15];
01642     static const double y[15];
01643     int operator()(const VectorXd &b, VectorXd &fvec)
01644     {
01645         assert(b.size()==4);
01646         assert(fvec.size()==15);
01647         for(int i=0; i<15; i++)
01648             fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i];
01649         return 0;
01650     }
01651     int df(const VectorXd &b, MatrixXd &fjac)
01652     {
01653         assert(b.size()==4);
01654         assert(fjac.rows()==15);
01655         assert(fjac.cols()==4);
01656         for(int i=0; i<15; i++) {
01657             double e = exp(b[1]-b[2]*x[i]);
01658             double power = -1./b[3];
01659             fjac(i,0) = pow(1.+e, power);
01660             fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.);
01661             fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.);
01662             fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power);
01663         }
01664         return 0;
01665     }
01666 };
01667 const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. };
01668 const double rat43_functor::y[15] = { 16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41 };
01669 
01670 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml
01671 void testNistRat43(void)
01672 {
01673   const int n=4;
01674   int info;
01675 
01676   VectorXd x(n);
01677 
01678   /*
01679    * First try
01680    */
01681   x<< 100., 10., 1., 1.;
01682   // do the computation
01683   rat43_functor functor;
01684   LevenbergMarquardt<rat43_functor> lm(functor);
01685   lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
01686   lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
01687   info = lm.minimize(x);
01688 
01689   // check return value
01690   VERIFY_IS_EQUAL(info, 1);
01691   VERIFY_IS_EQUAL(lm.nfev, 27);
01692   VERIFY_IS_EQUAL(lm.njev, 20);
01693   // check norm^2
01694   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03);
01695   // check x
01696   VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
01697   VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
01698   VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
01699   VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
01700 
01701   /*
01702    * Second try
01703    */
01704   x<< 700., 5., 0.75, 1.3;
01705   // do the computation
01706   lm.resetParameters();
01707   lm.parameters.ftol = 1.E5*NumTraits<double>::epsilon();
01708   lm.parameters.xtol = 1.E5*NumTraits<double>::epsilon();
01709   info = lm.minimize(x);
01710 
01711   // check return value
01712   VERIFY_IS_EQUAL(info, 1);
01713   VERIFY_IS_EQUAL(lm.nfev, 9);
01714   VERIFY_IS_EQUAL(lm.njev, 8);
01715   // check norm^2
01716   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03);
01717   // check x
01718   VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
01719   VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
01720   VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
01721   VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
01722 }
01723 
01724 
01725 
01726 struct eckerle4_functor : Functor<double>
01727 {
01728     eckerle4_functor(void) : Functor<double>(3,35) {}
01729     static const double x[35];
01730     static const double y[35];
01731     int operator()(const VectorXd &b, VectorXd &fvec)
01732     {
01733         assert(b.size()==3);
01734         assert(fvec.size()==35);
01735         for(int i=0; i<35; i++)
01736             fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i];
01737         return 0;
01738     }
01739     int df(const VectorXd &b, MatrixXd &fjac)
01740     {
01741         assert(b.size()==3);
01742         assert(fjac.rows()==35);
01743         assert(fjac.cols()==3);
01744         for(int i=0; i<35; i++) {
01745             double b12 = b[1]*b[1];
01746             double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12);
01747             fjac(i,0) = e / b[1];
01748             fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12;
01749             fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12;
01750         }
01751         return 0;
01752     }
01753 };
01754 const double eckerle4_functor::x[35] = { 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0};
01755 const double eckerle4_functor::y[35] = { 0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 };
01756 
01757 // http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml
01758 void testNistEckerle4(void)
01759 {
01760   const int n=3;
01761   int info;
01762 
01763   VectorXd x(n);
01764 
01765   /*
01766    * First try
01767    */
01768   x<< 1., 10., 500.;
01769   // do the computation
01770   eckerle4_functor functor;
01771   LevenbergMarquardt<eckerle4_functor> lm(functor);
01772   info = lm.minimize(x);
01773 
01774   // check return value
01775   VERIFY_IS_EQUAL(info, 1);
01776   VERIFY_IS_EQUAL(lm.nfev, 18);
01777   VERIFY_IS_EQUAL(lm.njev, 15);
01778   // check norm^2
01779   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03);
01780   // check x
01781   VERIFY_IS_APPROX(x[0], 1.5543827178);
01782   VERIFY_IS_APPROX(x[1], 4.0888321754);
01783   VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
01784 
01785   /*
01786    * Second try
01787    */
01788   x<< 1.5, 5., 450.;
01789   // do the computation
01790   info = lm.minimize(x);
01791 
01792   // check return value
01793   VERIFY_IS_EQUAL(info, 1);
01794   VERIFY_IS_EQUAL(lm.nfev, 7);
01795   VERIFY_IS_EQUAL(lm.njev, 6);
01796   // check norm^2
01797   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03);
01798   // check x
01799   VERIFY_IS_APPROX(x[0], 1.5543827178);
01800   VERIFY_IS_APPROX(x[1], 4.0888321754);
01801   VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
01802 }
01803 
01804 void test_NonLinearOptimization()
01805 {
01806     // Tests using the examples provided by (c)minpack
01807     CALL_SUBTEST/*_1*/(testChkder());
01808     CALL_SUBTEST/*_1*/(testLmder1());
01809     CALL_SUBTEST/*_1*/(testLmder());
01810     CALL_SUBTEST/*_2*/(testHybrj1());
01811     CALL_SUBTEST/*_2*/(testHybrj());
01812     CALL_SUBTEST/*_2*/(testHybrd1());
01813     CALL_SUBTEST/*_2*/(testHybrd());
01814     CALL_SUBTEST/*_3*/(testLmstr1());
01815     CALL_SUBTEST/*_3*/(testLmstr());
01816     CALL_SUBTEST/*_3*/(testLmdif1());
01817     CALL_SUBTEST/*_3*/(testLmdif());
01818 
01819     // NIST tests, level of difficulty = "Lower"
01820     CALL_SUBTEST/*_4*/(testNistMisra1a());
01821     CALL_SUBTEST/*_4*/(testNistChwirut2());
01822 
01823     // NIST tests, level of difficulty = "Average"
01824     CALL_SUBTEST/*_5*/(testNistHahn1());
01825     CALL_SUBTEST/*_6*/(testNistMisra1d());
01826     CALL_SUBTEST/*_7*/(testNistMGH17());
01827     CALL_SUBTEST/*_8*/(testNistLanczos1());
01828 
01829     // NIST tests, level of difficulty = "Higher"
01830     CALL_SUBTEST/*_9*/(testNistRat42());
01831     CALL_SUBTEST/*_10*/(testNistMGH10());
01832     CALL_SUBTEST/*_11*/(testNistBoxBOD());
01833     CALL_SUBTEST/*_12*/(testNistMGH09());
01834     CALL_SUBTEST/*_13*/(testNistBennett5());
01835     CALL_SUBTEST/*_14*/(testNistThurber());
01836     CALL_SUBTEST/*_15*/(testNistRat43());
01837     CALL_SUBTEST/*_16*/(testNistEckerle4());
01838 }
01839 
01840 /*
01841  * Can be useful for debugging...
01842   printf("info, nfev : %d, %d\n", info, lm.nfev);
01843   printf("info, nfev, njev : %d, %d, %d\n", info, solver.nfev, solver.njev);
01844   printf("info, nfev : %d, %d\n", info, solver.nfev);
01845   printf("x[0] : %.32g\n", x[0]);
01846   printf("x[1] : %.32g\n", x[1]);
01847   printf("x[2] : %.32g\n", x[2]);
01848   printf("x[3] : %.32g\n", x[3]);
01849   printf("fvec.blueNorm() : %.32g\n", solver.fvec.blueNorm());
01850   printf("fvec.blueNorm() : %.32g\n", lm.fvec.blueNorm());
01851 
01852   printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev);
01853   printf("fvec.squaredNorm() : %.13g\n", lm.fvec.squaredNorm());
01854   std::cout << x << std::endl;
01855   std::cout.precision(9);
01856   std::cout << x[0] << std::endl;
01857   std::cout << x[1] << std::endl;
01858   std::cout << x[2] << std::endl;
01859   std::cout << x[3] << std::endl;
01860 */
01861 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:08