1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Thomas Capricelli <orzel (at) freehackers.org> 5 // Copyright (C) 2012 desire Nuentsa <desire.nuentsa_wakam (at) inria.fr 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 12 // FIXME: These tests all check for hard-coded values. Ideally, parameters and start estimates should be randomized. 13 14 15 #include <stdio.h> 16 17 #include "main.h" 18 #include <unsupported/Eigen/LevenbergMarquardt> 19 20 // This disables some useless Warnings on MSVC. 21 // It is intended to be done for this test only. 22 #include <Eigen/src/Core/util/DisableStupidWarnings.h> 23 24 using std::sqrt; 25 26 // tolerance for chekcing number of iterations 27 #define LM_EVAL_COUNT_TOL 4/3 28 29 struct lmder_functor : DenseFunctor<double> 30 { 31 lmder_functor(void): DenseFunctor<double>(3,15) {} 32 int operator()(const VectorXd &x, VectorXd &fvec) const 33 { 34 double tmp1, tmp2, tmp3; 35 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, 36 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 37 38 for (int i = 0; i < values(); i++) 39 { 40 tmp1 = i+1; 41 tmp2 = 16 - i - 1; 42 tmp3 = (i>=8)? tmp2 : tmp1; 43 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 44 } 45 return 0; 46 } 47 48 int df(const VectorXd &x, MatrixXd &fjac) const 49 { 50 double tmp1, tmp2, tmp3, tmp4; 51 for (int i = 0; i < values(); i++) 52 { 53 tmp1 = i+1; 54 tmp2 = 16 - i - 1; 55 tmp3 = (i>=8)? tmp2 : tmp1; 56 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 57 fjac(i,0) = -1; 58 fjac(i,1) = tmp1*tmp2/tmp4; 59 fjac(i,2) = tmp1*tmp3/tmp4; 60 } 61 return 0; 62 } 63 }; 64 65 void testLmder1() 66 { 67 int n=3, info; 68 69 VectorXd x; 70 71 /* the following starting values provide a rough fit. */ 72 x.setConstant(n, 1.); 73 74 // do the computation 75 lmder_functor functor; 76 LevenbergMarquardt<lmder_functor> lm(functor); 77 info = lm.lmder1(x); 78 79 // check return value 80 VERIFY_IS_EQUAL(info, 1); 81 VERIFY_IS_EQUAL(lm.nfev(), 6); 82 VERIFY_IS_EQUAL(lm.njev(), 5); 83 84 // check norm 85 VERIFY_IS_APPROX(lm.fvec().blueNorm(), 0.09063596); 86 87 // check x 88 VectorXd x_ref(n); 89 x_ref << 0.08241058, 1.133037, 2.343695; 90 VERIFY_IS_APPROX(x, x_ref); 91 } 92 93 void testLmder() 94 { 95 const int m=15, n=3; 96 int info; 97 double fnorm, covfac; 98 VectorXd x; 99 100 /* the following starting values provide a rough fit. */ 101 x.setConstant(n, 1.); 102 103 // do the computation 104 lmder_functor functor; 105 LevenbergMarquardt<lmder_functor> lm(functor); 106 info = lm.minimize(x); 107 108 // check return values 109 VERIFY_IS_EQUAL(info, 1); 110 VERIFY_IS_EQUAL(lm.nfev(), 6); 111 VERIFY_IS_EQUAL(lm.njev(), 5); 112 113 // check norm 114 fnorm = lm.fvec().blueNorm(); 115 VERIFY_IS_APPROX(fnorm, 0.09063596); 116 117 // check x 118 VectorXd x_ref(n); 119 x_ref << 0.08241058, 1.133037, 2.343695; 120 VERIFY_IS_APPROX(x, x_ref); 121 122 // check covariance 123 covfac = fnorm*fnorm/(m-n); 124 internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm 125 126 MatrixXd cov_ref(n,n); 127 cov_ref << 128 0.0001531202, 0.002869941, -0.002656662, 129 0.002869941, 0.09480935, -0.09098995, 130 -0.002656662, -0.09098995, 0.08778727; 131 132 // std::cout << fjac*covfac << std::endl; 133 134 MatrixXd cov; 135 cov = covfac*lm.matrixR().topLeftCorner<n,n>(); 136 VERIFY_IS_APPROX( cov, cov_ref); 137 // TODO: why isn't this allowed ? : 138 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 139 } 140 141 struct lmdif_functor : DenseFunctor<double> 142 { 143 lmdif_functor(void) : DenseFunctor<double>(3,15) {} 144 int operator()(const VectorXd &x, VectorXd &fvec) const 145 { 146 int i; 147 double tmp1,tmp2,tmp3; 148 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, 149 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; 150 151 assert(x.size()==3); 152 assert(fvec.size()==15); 153 for (i=0; i<15; i++) 154 { 155 tmp1 = i+1; 156 tmp2 = 15 - i; 157 tmp3 = tmp1; 158 159 if (i >= 8) tmp3 = tmp2; 160 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 161 } 162 return 0; 163 } 164 }; 165 166 void testLmdif1() 167 { 168 const int n=3; 169 int info; 170 171 VectorXd x(n), fvec(15); 172 173 /* the following starting values provide a rough fit. */ 174 x.setConstant(n, 1.); 175 176 // do the computation 177 lmdif_functor functor; 178 DenseIndex nfev; 179 info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev); 180 181 // check return value 182 VERIFY_IS_EQUAL(info, 1); 183 // VERIFY_IS_EQUAL(nfev, 26); 184 185 // check norm 186 functor(x, fvec); 187 VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); 188 189 // check x 190 VectorXd x_ref(n); 191 x_ref << 0.0824106, 1.1330366, 2.3436947; 192 VERIFY_IS_APPROX(x, x_ref); 193 194 } 195 196 void testLmdif() 197 { 198 const int m=15, n=3; 199 int info; 200 double fnorm, covfac; 201 VectorXd x(n); 202 203 /* the following starting values provide a rough fit. */ 204 x.setConstant(n, 1.); 205 206 // do the computation 207 lmdif_functor functor; 208 NumericalDiff<lmdif_functor> numDiff(functor); 209 LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff); 210 info = lm.minimize(x); 211 212 // check return values 213 VERIFY_IS_EQUAL(info, 1); 214 // VERIFY_IS_EQUAL(lm.nfev(), 26); 215 216 // check norm 217 fnorm = lm.fvec().blueNorm(); 218 VERIFY_IS_APPROX(fnorm, 0.09063596); 219 220 // check x 221 VectorXd x_ref(n); 222 x_ref << 0.08241058, 1.133037, 2.343695; 223 VERIFY_IS_APPROX(x, x_ref); 224 225 // check covariance 226 covfac = fnorm*fnorm/(m-n); 227 internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm 228 229 MatrixXd cov_ref(n,n); 230 cov_ref << 231 0.0001531202, 0.002869942, -0.002656662, 232 0.002869942, 0.09480937, -0.09098997, 233 -0.002656662, -0.09098997, 0.08778729; 234 235 // std::cout << fjac*covfac << std::endl; 236 237 MatrixXd cov; 238 cov = covfac*lm.matrixR().topLeftCorner<n,n>(); 239 VERIFY_IS_APPROX( cov, cov_ref); 240 // TODO: why isn't this allowed ? : 241 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 242 } 243 244 struct chwirut2_functor : DenseFunctor<double> 245 { 246 chwirut2_functor(void) : DenseFunctor<double>(3,54) {} 247 static const double m_x[54]; 248 static const double m_y[54]; 249 int operator()(const VectorXd &b, VectorXd &fvec) 250 { 251 int i; 252 253 assert(b.size()==3); 254 assert(fvec.size()==54); 255 for(i=0; i<54; i++) { 256 double x = m_x[i]; 257 fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; 258 } 259 return 0; 260 } 261 int df(const VectorXd &b, MatrixXd &fjac) 262 { 263 assert(b.size()==3); 264 assert(fjac.rows()==54); 265 assert(fjac.cols()==3); 266 for(int i=0; i<54; i++) { 267 double x = m_x[i]; 268 double factor = 1./(b[1]+b[2]*x); 269 double e = exp(-b[0]*x); 270 fjac(i,0) = -x*e*factor; 271 fjac(i,1) = -e*factor*factor; 272 fjac(i,2) = -x*e*factor*factor; 273 } 274 return 0; 275 } 276 }; 277 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}; 278 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 }; 279 280 // http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml 281 void testNistChwirut2(void) 282 { 283 const int n=3; 284 LevenbergMarquardtSpace::Status info; 285 286 VectorXd x(n); 287 288 /* 289 * First try 290 */ 291 x<< 0.1, 0.01, 0.02; 292 // do the computation 293 chwirut2_functor functor; 294 LevenbergMarquardt<chwirut2_functor> lm(functor); 295 info = lm.minimize(x); 296 297 // check return value 298 VERIFY_IS_EQUAL(info, 1); 299 // VERIFY_IS_EQUAL(lm.nfev(), 10); 300 VERIFY_IS_EQUAL(lm.njev(), 8); 301 // check norm^2 302 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); 303 // check x 304 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 305 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 306 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 307 308 /* 309 * Second try 310 */ 311 x<< 0.15, 0.008, 0.010; 312 // do the computation 313 lm.resetParameters(); 314 lm.setFtol(1.E6*NumTraits<double>::epsilon()); 315 lm.setXtol(1.E6*NumTraits<double>::epsilon()); 316 info = lm.minimize(x); 317 318 // check return value 319 VERIFY_IS_EQUAL(info, 1); 320 // VERIFY_IS_EQUAL(lm.nfev(), 7); 321 VERIFY_IS_EQUAL(lm.njev(), 6); 322 // check norm^2 323 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); 324 // check x 325 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 326 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 327 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 328 } 329 330 331 struct misra1a_functor : DenseFunctor<double> 332 { 333 misra1a_functor(void) : DenseFunctor<double>(2,14) {} 334 static const double m_x[14]; 335 static const double m_y[14]; 336 int operator()(const VectorXd &b, VectorXd &fvec) 337 { 338 assert(b.size()==2); 339 assert(fvec.size()==14); 340 for(int i=0; i<14; i++) { 341 fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; 342 } 343 return 0; 344 } 345 int df(const VectorXd &b, MatrixXd &fjac) 346 { 347 assert(b.size()==2); 348 assert(fjac.rows()==14); 349 assert(fjac.cols()==2); 350 for(int i=0; i<14; i++) { 351 fjac(i,0) = (1.-exp(-b[1]*m_x[i])); 352 fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); 353 } 354 return 0; 355 } 356 }; 357 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}; 358 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}; 359 360 // http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml 361 void testNistMisra1a(void) 362 { 363 const int n=2; 364 int info; 365 366 VectorXd x(n); 367 368 /* 369 * First try 370 */ 371 x<< 500., 0.0001; 372 // do the computation 373 misra1a_functor functor; 374 LevenbergMarquardt<misra1a_functor> lm(functor); 375 info = lm.minimize(x); 376 377 // check return value 378 VERIFY_IS_EQUAL(info, 1); 379 VERIFY_IS_EQUAL(lm.nfev(), 19); 380 VERIFY_IS_EQUAL(lm.njev(), 15); 381 // check norm^2 382 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); 383 // check x 384 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 385 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 386 387 /* 388 * Second try 389 */ 390 x<< 250., 0.0005; 391 // do the computation 392 info = lm.minimize(x); 393 394 // check return value 395 VERIFY_IS_EQUAL(info, 1); 396 VERIFY_IS_EQUAL(lm.nfev(), 5); 397 VERIFY_IS_EQUAL(lm.njev(), 4); 398 // check norm^2 399 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); 400 // check x 401 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 402 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 403 } 404 405 struct hahn1_functor : DenseFunctor<double> 406 { 407 hahn1_functor(void) : DenseFunctor<double>(7,236) {} 408 static const double m_x[236]; 409 int operator()(const VectorXd &b, VectorXd &fvec) 410 { 411 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 , 412 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 , 413 12.596E0 , 414 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 }; 415 416 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 417 418 assert(b.size()==7); 419 assert(fvec.size()==236); 420 for(int i=0; i<236; i++) { 421 double x=m_x[i], xx=x*x, xxx=xx*x; 422 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]; 423 } 424 return 0; 425 } 426 427 int df(const VectorXd &b, MatrixXd &fjac) 428 { 429 assert(b.size()==7); 430 assert(fjac.rows()==236); 431 assert(fjac.cols()==7); 432 for(int i=0; i<236; i++) { 433 double x=m_x[i], xx=x*x, xxx=xx*x; 434 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 435 fjac(i,0) = 1.*fact; 436 fjac(i,1) = x*fact; 437 fjac(i,2) = xx*fact; 438 fjac(i,3) = xxx*fact; 439 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 440 fjac(i,4) = x*fact; 441 fjac(i,5) = xx*fact; 442 fjac(i,6) = xxx*fact; 443 } 444 return 0; 445 } 446 }; 447 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 , 448 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 , 449 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}; 450 451 // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml 452 void testNistHahn1(void) 453 { 454 const int n=7; 455 int info; 456 457 VectorXd x(n); 458 459 /* 460 * First try 461 */ 462 x<< 10., -1., .05, -.00001, -.05, .001, -.000001; 463 // do the computation 464 hahn1_functor functor; 465 LevenbergMarquardt<hahn1_functor> lm(functor); 466 info = lm.minimize(x); 467 468 // check return value 469 VERIFY_IS_EQUAL(info, 1); 470 VERIFY_IS_EQUAL(lm.nfev(), 11); 471 VERIFY_IS_EQUAL(lm.njev(), 10); 472 // check norm^2 473 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); 474 // check x 475 VERIFY_IS_APPROX(x[0], 1.0776351733E+00); 476 VERIFY_IS_APPROX(x[1],-1.2269296921E-01); 477 VERIFY_IS_APPROX(x[2], 4.0863750610E-03); 478 VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 479 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 480 VERIFY_IS_APPROX(x[5], 2.4053735503E-04); 481 VERIFY_IS_APPROX(x[6],-1.2314450199E-07); 482 483 /* 484 * Second try 485 */ 486 x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; 487 // do the computation 488 info = lm.minimize(x); 489 490 // check return value 491 VERIFY_IS_EQUAL(info, 1); 492 // VERIFY_IS_EQUAL(lm.nfev(), 11); 493 VERIFY_IS_EQUAL(lm.njev(), 10); 494 // check norm^2 495 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); 496 // check x 497 VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 498 VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 499 VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 500 VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 501 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 502 VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 503 VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 504 505 } 506 507 struct misra1d_functor : DenseFunctor<double> 508 { 509 misra1d_functor(void) : DenseFunctor<double>(2,14) {} 510 static const double x[14]; 511 static const double y[14]; 512 int operator()(const VectorXd &b, VectorXd &fvec) 513 { 514 assert(b.size()==2); 515 assert(fvec.size()==14); 516 for(int i=0; i<14; i++) { 517 fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; 518 } 519 return 0; 520 } 521 int df(const VectorXd &b, MatrixXd &fjac) 522 { 523 assert(b.size()==2); 524 assert(fjac.rows()==14); 525 assert(fjac.cols()==2); 526 for(int i=0; i<14; i++) { 527 double den = 1.+b[1]*x[i]; 528 fjac(i,0) = b[1]*x[i] / den; 529 fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; 530 } 531 return 0; 532 } 533 }; 534 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}; 535 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}; 536 537 // http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml 538 void testNistMisra1d(void) 539 { 540 const int n=2; 541 int info; 542 543 VectorXd x(n); 544 545 /* 546 * First try 547 */ 548 x<< 500., 0.0001; 549 // do the computation 550 misra1d_functor functor; 551 LevenbergMarquardt<misra1d_functor> lm(functor); 552 info = lm.minimize(x); 553 554 // check return value 555 VERIFY_IS_EQUAL(info, 1); 556 VERIFY_IS_EQUAL(lm.nfev(), 9); 557 VERIFY_IS_EQUAL(lm.njev(), 7); 558 // check norm^2 559 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); 560 // check x 561 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 562 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 563 564 /* 565 * Second try 566 */ 567 x<< 450., 0.0003; 568 // do the computation 569 info = lm.minimize(x); 570 571 // check return value 572 VERIFY_IS_EQUAL(info, 1); 573 VERIFY_IS_EQUAL(lm.nfev(), 4); 574 VERIFY_IS_EQUAL(lm.njev(), 3); 575 // check norm^2 576 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); 577 // check x 578 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 579 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 580 } 581 582 583 struct lanczos1_functor : DenseFunctor<double> 584 { 585 lanczos1_functor(void) : DenseFunctor<double>(6,24) {} 586 static const double x[24]; 587 static const double y[24]; 588 int operator()(const VectorXd &b, VectorXd &fvec) 589 { 590 assert(b.size()==6); 591 assert(fvec.size()==24); 592 for(int i=0; i<24; i++) 593 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]; 594 return 0; 595 } 596 int df(const VectorXd &b, MatrixXd &fjac) 597 { 598 assert(b.size()==6); 599 assert(fjac.rows()==24); 600 assert(fjac.cols()==6); 601 for(int i=0; i<24; i++) { 602 fjac(i,0) = exp(-b[1]*x[i]); 603 fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); 604 fjac(i,2) = exp(-b[3]*x[i]); 605 fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); 606 fjac(i,4) = exp(-b[5]*x[i]); 607 fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); 608 } 609 return 0; 610 } 611 }; 612 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 }; 613 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 }; 614 615 // http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml 616 void testNistLanczos1(void) 617 { 618 const int n=6; 619 LevenbergMarquardtSpace::Status info; 620 621 VectorXd x(n); 622 623 /* 624 * First try 625 */ 626 x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; 627 // do the computation 628 lanczos1_functor functor; 629 LevenbergMarquardt<lanczos1_functor> lm(functor); 630 info = lm.minimize(x); 631 632 // check return value 633 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall); 634 VERIFY_IS_EQUAL(lm.nfev(), 79); 635 VERIFY_IS_EQUAL(lm.njev(), 72); 636 // check norm^2 637 VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25); 638 // check x 639 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 640 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 641 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 642 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 643 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 644 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 645 646 /* 647 * Second try 648 */ 649 x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; 650 // do the computation 651 info = lm.minimize(x); 652 653 // check return value 654 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall); 655 VERIFY_IS_EQUAL(lm.nfev(), 9); 656 VERIFY_IS_EQUAL(lm.njev(), 8); 657 // check norm^2 658 VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25); 659 // check x 660 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 661 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 662 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 663 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 664 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 665 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 666 667 } 668 669 struct rat42_functor : DenseFunctor<double> 670 { 671 rat42_functor(void) : DenseFunctor<double>(3,9) {} 672 static const double x[9]; 673 static const double y[9]; 674 int operator()(const VectorXd &b, VectorXd &fvec) 675 { 676 assert(b.size()==3); 677 assert(fvec.size()==9); 678 for(int i=0; i<9; i++) { 679 fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; 680 } 681 return 0; 682 } 683 684 int df(const VectorXd &b, MatrixXd &fjac) 685 { 686 assert(b.size()==3); 687 assert(fjac.rows()==9); 688 assert(fjac.cols()==3); 689 for(int i=0; i<9; i++) { 690 double e = exp(b[1]-b[2]*x[i]); 691 fjac(i,0) = 1./(1.+e); 692 fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); 693 fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); 694 } 695 return 0; 696 } 697 }; 698 const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; 699 const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; 700 701 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml 702 void testNistRat42(void) 703 { 704 const int n=3; 705 LevenbergMarquardtSpace::Status info; 706 707 VectorXd x(n); 708 709 /* 710 * First try 711 */ 712 x<< 100., 1., 0.1; 713 // do the computation 714 rat42_functor functor; 715 LevenbergMarquardt<rat42_functor> lm(functor); 716 info = lm.minimize(x); 717 718 // check return value 719 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 720 VERIFY_IS_EQUAL(lm.nfev(), 10); 721 VERIFY_IS_EQUAL(lm.njev(), 8); 722 // check norm^2 723 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); 724 // check x 725 VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 726 VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 727 VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 728 729 /* 730 * Second try 731 */ 732 x<< 75., 2.5, 0.07; 733 // do the computation 734 info = lm.minimize(x); 735 736 // check return value 737 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 738 VERIFY_IS_EQUAL(lm.nfev(), 6); 739 VERIFY_IS_EQUAL(lm.njev(), 5); 740 // check norm^2 741 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); 742 // check x 743 VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 744 VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 745 VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 746 } 747 748 struct MGH10_functor : DenseFunctor<double> 749 { 750 MGH10_functor(void) : DenseFunctor<double>(3,16) {} 751 static const double x[16]; 752 static const double y[16]; 753 int operator()(const VectorXd &b, VectorXd &fvec) 754 { 755 assert(b.size()==3); 756 assert(fvec.size()==16); 757 for(int i=0; i<16; i++) 758 fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; 759 return 0; 760 } 761 int df(const VectorXd &b, MatrixXd &fjac) 762 { 763 assert(b.size()==3); 764 assert(fjac.rows()==16); 765 assert(fjac.cols()==3); 766 for(int i=0; i<16; i++) { 767 double factor = 1./(x[i]+b[2]); 768 double e = exp(b[1]*factor); 769 fjac(i,0) = e; 770 fjac(i,1) = b[0]*factor*e; 771 fjac(i,2) = -b[1]*b[0]*factor*factor*e; 772 } 773 return 0; 774 } 775 }; 776 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 }; 777 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 }; 778 779 // http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml 780 void testNistMGH10(void) 781 { 782 const int n=3; 783 LevenbergMarquardtSpace::Status info; 784 785 VectorXd x(n); 786 787 /* 788 * First try 789 */ 790 x<< 2., 400000., 25000.; 791 // do the computation 792 MGH10_functor functor; 793 LevenbergMarquardt<MGH10_functor> lm(functor); 794 info = lm.minimize(x); 795 ++g_test_level; 796 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 797 --g_test_level; 798 // was: VERIFY_IS_EQUAL(info, 1); 799 800 // check norm^2 801 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); 802 // check x 803 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 804 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 805 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 806 807 // check return value 808 809 ++g_test_level; 810 VERIFY_IS_EQUAL(lm.nfev(), 284 ); 811 VERIFY_IS_EQUAL(lm.njev(), 249 ); 812 --g_test_level; 813 VERIFY(lm.nfev() < 284 * LM_EVAL_COUNT_TOL); 814 VERIFY(lm.njev() < 249 * LM_EVAL_COUNT_TOL); 815 816 /* 817 * Second try 818 */ 819 x<< 0.02, 4000., 250.; 820 // do the computation 821 info = lm.minimize(x); 822 ++g_test_level; 823 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 824 // was: VERIFY_IS_EQUAL(info, 1); 825 --g_test_level; 826 827 // check norm^2 828 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); 829 // check x 830 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 831 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 832 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 833 834 // check return value 835 ++g_test_level; 836 VERIFY_IS_EQUAL(lm.nfev(), 126); 837 VERIFY_IS_EQUAL(lm.njev(), 116); 838 --g_test_level; 839 VERIFY(lm.nfev() < 126 * LM_EVAL_COUNT_TOL); 840 VERIFY(lm.njev() < 116 * LM_EVAL_COUNT_TOL); 841 } 842 843 844 struct BoxBOD_functor : DenseFunctor<double> 845 { 846 BoxBOD_functor(void) : DenseFunctor<double>(2,6) {} 847 static const double x[6]; 848 int operator()(const VectorXd &b, VectorXd &fvec) 849 { 850 static const double y[6] = { 109., 149., 149., 191., 213., 224. }; 851 assert(b.size()==2); 852 assert(fvec.size()==6); 853 for(int i=0; i<6; i++) 854 fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; 855 return 0; 856 } 857 int df(const VectorXd &b, MatrixXd &fjac) 858 { 859 assert(b.size()==2); 860 assert(fjac.rows()==6); 861 assert(fjac.cols()==2); 862 for(int i=0; i<6; i++) { 863 double e = exp(-b[1]*x[i]); 864 fjac(i,0) = 1.-e; 865 fjac(i,1) = b[0]*x[i]*e; 866 } 867 return 0; 868 } 869 }; 870 const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; 871 872 // http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml 873 void testNistBoxBOD(void) 874 { 875 const int n=2; 876 int info; 877 878 VectorXd x(n); 879 880 /* 881 * First try 882 */ 883 x<< 1., 1.; 884 // do the computation 885 BoxBOD_functor functor; 886 LevenbergMarquardt<BoxBOD_functor> lm(functor); 887 lm.setFtol(1.E6*NumTraits<double>::epsilon()); 888 lm.setXtol(1.E6*NumTraits<double>::epsilon()); 889 lm.setFactor(10); 890 info = lm.minimize(x); 891 892 // check norm^2 893 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); 894 // check x 895 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 896 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 897 898 // check return value 899 VERIFY_IS_EQUAL(info, 1); 900 VERIFY(lm.nfev() < 31); // 31 901 VERIFY(lm.njev() < 25); // 25 902 903 /* 904 * Second try 905 */ 906 x<< 100., 0.75; 907 // do the computation 908 lm.resetParameters(); 909 lm.setFtol(NumTraits<double>::epsilon()); 910 lm.setXtol( NumTraits<double>::epsilon()); 911 info = lm.minimize(x); 912 913 // check return value 914 VERIFY_IS_EQUAL(info, 1); 915 ++g_test_level; 916 VERIFY_IS_EQUAL(lm.nfev(), 16 ); 917 VERIFY_IS_EQUAL(lm.njev(), 15 ); 918 --g_test_level; 919 VERIFY(lm.nfev() < 16 * LM_EVAL_COUNT_TOL); 920 VERIFY(lm.njev() < 15 * LM_EVAL_COUNT_TOL); 921 // check norm^2 922 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); 923 // check x 924 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 925 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 926 } 927 928 struct MGH17_functor : DenseFunctor<double> 929 { 930 MGH17_functor(void) : DenseFunctor<double>(5,33) {} 931 static const double x[33]; 932 static const double y[33]; 933 int operator()(const VectorXd &b, VectorXd &fvec) 934 { 935 assert(b.size()==5); 936 assert(fvec.size()==33); 937 for(int i=0; i<33; i++) 938 fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; 939 return 0; 940 } 941 int df(const VectorXd &b, MatrixXd &fjac) 942 { 943 assert(b.size()==5); 944 assert(fjac.rows()==33); 945 assert(fjac.cols()==5); 946 for(int i=0; i<33; i++) { 947 fjac(i,0) = 1.; 948 fjac(i,1) = exp(-b[3]*x[i]); 949 fjac(i,2) = exp(-b[4]*x[i]); 950 fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); 951 fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); 952 } 953 return 0; 954 } 955 }; 956 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 }; 957 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 }; 958 959 // http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml 960 void testNistMGH17(void) 961 { 962 const int n=5; 963 int info; 964 965 VectorXd x(n); 966 967 /* 968 * First try 969 */ 970 x<< 50., 150., -100., 1., 2.; 971 // do the computation 972 MGH17_functor functor; 973 LevenbergMarquardt<MGH17_functor> lm(functor); 974 lm.setFtol(NumTraits<double>::epsilon()); 975 lm.setXtol(NumTraits<double>::epsilon()); 976 lm.setMaxfev(1000); 977 info = lm.minimize(x); 978 979 // check norm^2 980 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); 981 // check x 982 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 983 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 984 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 985 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 986 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 987 988 // check return value 989 // VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success) 990 VERIFY(lm.nfev() < 700 ); // 602 991 VERIFY(lm.njev() < 600 ); // 545 992 993 /* 994 * Second try 995 */ 996 x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; 997 // do the computation 998 lm.resetParameters(); 999 info = lm.minimize(x); 1000 1001 // check return value 1002 VERIFY_IS_EQUAL(info, 1); 1003 VERIFY_IS_EQUAL(lm.nfev(), 18); 1004 VERIFY_IS_EQUAL(lm.njev(), 15); 1005 // check norm^2 1006 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); 1007 // check x 1008 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1009 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1010 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1011 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1012 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1013 } 1014 1015 struct MGH09_functor : DenseFunctor<double> 1016 { 1017 MGH09_functor(void) : DenseFunctor<double>(4,11) {} 1018 static const double _x[11]; 1019 static const double y[11]; 1020 int operator()(const VectorXd &b, VectorXd &fvec) 1021 { 1022 assert(b.size()==4); 1023 assert(fvec.size()==11); 1024 for(int i=0; i<11; i++) { 1025 double x = _x[i], xx=x*x; 1026 fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; 1027 } 1028 return 0; 1029 } 1030 int df(const VectorXd &b, MatrixXd &fjac) 1031 { 1032 assert(b.size()==4); 1033 assert(fjac.rows()==11); 1034 assert(fjac.cols()==4); 1035 for(int i=0; i<11; i++) { 1036 double x = _x[i], xx=x*x; 1037 double factor = 1./(xx+x*b[2]+b[3]); 1038 fjac(i,0) = (xx+x*b[1]) * factor; 1039 fjac(i,1) = b[0]*x* factor; 1040 fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; 1041 fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; 1042 } 1043 return 0; 1044 } 1045 }; 1046 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 }; 1047 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 }; 1048 1049 // http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml 1050 void testNistMGH09(void) 1051 { 1052 const int n=4; 1053 int info; 1054 1055 VectorXd x(n); 1056 1057 /* 1058 * First try 1059 */ 1060 x<< 25., 39, 41.5, 39.; 1061 // do the computation 1062 MGH09_functor functor; 1063 LevenbergMarquardt<MGH09_functor> lm(functor); 1064 lm.setMaxfev(1000); 1065 info = lm.minimize(x); 1066 1067 // check norm^2 1068 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); 1069 // check x 1070 VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 1071 VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 1072 VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 1073 VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 1074 // check return value 1075 VERIFY_IS_EQUAL(info, 1); 1076 VERIFY(lm.nfev() < 510 ); // 490 1077 VERIFY(lm.njev() < 400 ); // 376 1078 1079 /* 1080 * Second try 1081 */ 1082 x<< 0.25, 0.39, 0.415, 0.39; 1083 // do the computation 1084 lm.resetParameters(); 1085 info = lm.minimize(x); 1086 1087 // check return value 1088 VERIFY_IS_EQUAL(info, 1); 1089 VERIFY_IS_EQUAL(lm.nfev(), 18); 1090 VERIFY_IS_EQUAL(lm.njev(), 16); 1091 // check norm^2 1092 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); 1093 // check x 1094 VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 1095 VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 1096 VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 1097 VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 1098 } 1099 1100 1101 1102 struct Bennett5_functor : DenseFunctor<double> 1103 { 1104 Bennett5_functor(void) : DenseFunctor<double>(3,154) {} 1105 static const double x[154]; 1106 static const double y[154]; 1107 int operator()(const VectorXd &b, VectorXd &fvec) 1108 { 1109 assert(b.size()==3); 1110 assert(fvec.size()==154); 1111 for(int i=0; i<154; i++) 1112 fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; 1113 return 0; 1114 } 1115 int df(const VectorXd &b, MatrixXd &fjac) 1116 { 1117 assert(b.size()==3); 1118 assert(fjac.rows()==154); 1119 assert(fjac.cols()==3); 1120 for(int i=0; i<154; i++) { 1121 double e = pow(b[1]+x[i],-1./b[2]); 1122 fjac(i,0) = e; 1123 fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); 1124 fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; 1125 } 1126 return 0; 1127 } 1128 }; 1129 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, 1130 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 }; 1131 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 1132 ,-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 , 1133 -31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; 1134 1135 // http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml 1136 void testNistBennett5(void) 1137 { 1138 const int n=3; 1139 int info; 1140 1141 VectorXd x(n); 1142 1143 /* 1144 * First try 1145 */ 1146 x<< -2000., 50., 0.8; 1147 // do the computation 1148 Bennett5_functor functor; 1149 LevenbergMarquardt<Bennett5_functor> lm(functor); 1150 lm.setMaxfev(1000); 1151 info = lm.minimize(x); 1152 1153 // check return value 1154 VERIFY_IS_EQUAL(info, 1); 1155 VERIFY_IS_EQUAL(lm.nfev(), 758); 1156 VERIFY_IS_EQUAL(lm.njev(), 744); 1157 // check norm^2 1158 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); 1159 // check x 1160 VERIFY_IS_APPROX(x[0], -2.5235058043E+03); 1161 VERIFY_IS_APPROX(x[1], 4.6736564644E+01); 1162 VERIFY_IS_APPROX(x[2], 9.3218483193E-01); 1163 /* 1164 * Second try 1165 */ 1166 x<< -1500., 45., 0.85; 1167 // do the computation 1168 lm.resetParameters(); 1169 info = lm.minimize(x); 1170 1171 // check return value 1172 VERIFY_IS_EQUAL(info, 1); 1173 VERIFY_IS_EQUAL(lm.nfev(), 203); 1174 VERIFY_IS_EQUAL(lm.njev(), 192); 1175 // check norm^2 1176 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); 1177 // check x 1178 VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 1179 VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); 1180 VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); 1181 } 1182 1183 struct thurber_functor : DenseFunctor<double> 1184 { 1185 thurber_functor(void) : DenseFunctor<double>(7,37) {} 1186 static const double _x[37]; 1187 static const double _y[37]; 1188 int operator()(const VectorXd &b, VectorXd &fvec) 1189 { 1190 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 1191 assert(b.size()==7); 1192 assert(fvec.size()==37); 1193 for(int i=0; i<37; i++) { 1194 double x=_x[i], xx=x*x, xxx=xx*x; 1195 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]; 1196 } 1197 return 0; 1198 } 1199 int df(const VectorXd &b, MatrixXd &fjac) 1200 { 1201 assert(b.size()==7); 1202 assert(fjac.rows()==37); 1203 assert(fjac.cols()==7); 1204 for(int i=0; i<37; i++) { 1205 double x=_x[i], xx=x*x, xxx=xx*x; 1206 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 1207 fjac(i,0) = 1.*fact; 1208 fjac(i,1) = x*fact; 1209 fjac(i,2) = xx*fact; 1210 fjac(i,3) = xxx*fact; 1211 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 1212 fjac(i,4) = x*fact; 1213 fjac(i,5) = xx*fact; 1214 fjac(i,6) = xxx*fact; 1215 } 1216 return 0; 1217 } 1218 }; 1219 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 }; 1220 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}; 1221 1222 // http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml 1223 void testNistThurber(void) 1224 { 1225 const int n=7; 1226 int info; 1227 1228 VectorXd x(n); 1229 1230 /* 1231 * First try 1232 */ 1233 x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; 1234 // do the computation 1235 thurber_functor functor; 1236 LevenbergMarquardt<thurber_functor> lm(functor); 1237 lm.setFtol(1.E4*NumTraits<double>::epsilon()); 1238 lm.setXtol(1.E4*NumTraits<double>::epsilon()); 1239 info = lm.minimize(x); 1240 1241 // check return value 1242 VERIFY_IS_EQUAL(info, 1); 1243 VERIFY_IS_EQUAL(lm.nfev(), 39); 1244 VERIFY_IS_EQUAL(lm.njev(), 36); 1245 // check norm^2 1246 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); 1247 // check x 1248 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1249 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1250 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1251 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1252 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1253 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1254 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1255 1256 /* 1257 * Second try 1258 */ 1259 x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; 1260 // do the computation 1261 lm.resetParameters(); 1262 lm.setFtol(1.E4*NumTraits<double>::epsilon()); 1263 lm.setXtol(1.E4*NumTraits<double>::epsilon()); 1264 info = lm.minimize(x); 1265 1266 // check return value 1267 VERIFY_IS_EQUAL(info, 1); 1268 VERIFY_IS_EQUAL(lm.nfev(), 29); 1269 VERIFY_IS_EQUAL(lm.njev(), 28); 1270 // check norm^2 1271 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); 1272 // check x 1273 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1274 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1275 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1276 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1277 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1278 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1279 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1280 } 1281 1282 struct rat43_functor : DenseFunctor<double> 1283 { 1284 rat43_functor(void) : DenseFunctor<double>(4,15) {} 1285 static const double x[15]; 1286 static const double y[15]; 1287 int operator()(const VectorXd &b, VectorXd &fvec) 1288 { 1289 assert(b.size()==4); 1290 assert(fvec.size()==15); 1291 for(int i=0; i<15; i++) 1292 fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; 1293 return 0; 1294 } 1295 int df(const VectorXd &b, MatrixXd &fjac) 1296 { 1297 assert(b.size()==4); 1298 assert(fjac.rows()==15); 1299 assert(fjac.cols()==4); 1300 for(int i=0; i<15; i++) { 1301 double e = exp(b[1]-b[2]*x[i]); 1302 double power = -1./b[3]; 1303 fjac(i,0) = pow(1.+e, power); 1304 fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); 1305 fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); 1306 fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); 1307 } 1308 return 0; 1309 } 1310 }; 1311 const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; 1312 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 }; 1313 1314 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml 1315 void testNistRat43(void) 1316 { 1317 const int n=4; 1318 int info; 1319 1320 VectorXd x(n); 1321 1322 /* 1323 * First try 1324 */ 1325 x<< 100., 10., 1., 1.; 1326 // do the computation 1327 rat43_functor functor; 1328 LevenbergMarquardt<rat43_functor> lm(functor); 1329 lm.setFtol(1.E6*NumTraits<double>::epsilon()); 1330 lm.setXtol(1.E6*NumTraits<double>::epsilon()); 1331 info = lm.minimize(x); 1332 1333 // check return value 1334 VERIFY_IS_EQUAL(info, 1); 1335 VERIFY_IS_EQUAL(lm.nfev(), 27); 1336 VERIFY_IS_EQUAL(lm.njev(), 20); 1337 // check norm^2 1338 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); 1339 // check x 1340 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1341 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1342 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1343 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1344 1345 /* 1346 * Second try 1347 */ 1348 x<< 700., 5., 0.75, 1.3; 1349 // do the computation 1350 lm.resetParameters(); 1351 lm.setFtol(1.E5*NumTraits<double>::epsilon()); 1352 lm.setXtol(1.E5*NumTraits<double>::epsilon()); 1353 info = lm.minimize(x); 1354 1355 // check return value 1356 VERIFY_IS_EQUAL(info, 1); 1357 VERIFY_IS_EQUAL(lm.nfev(), 9); 1358 VERIFY_IS_EQUAL(lm.njev(), 8); 1359 // check norm^2 1360 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); 1361 // check x 1362 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1363 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1364 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1365 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1366 } 1367 1368 1369 1370 struct eckerle4_functor : DenseFunctor<double> 1371 { 1372 eckerle4_functor(void) : DenseFunctor<double>(3,35) {} 1373 static const double x[35]; 1374 static const double y[35]; 1375 int operator()(const VectorXd &b, VectorXd &fvec) 1376 { 1377 assert(b.size()==3); 1378 assert(fvec.size()==35); 1379 for(int i=0; i<35; i++) 1380 fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; 1381 return 0; 1382 } 1383 int df(const VectorXd &b, MatrixXd &fjac) 1384 { 1385 assert(b.size()==3); 1386 assert(fjac.rows()==35); 1387 assert(fjac.cols()==3); 1388 for(int i=0; i<35; i++) { 1389 double b12 = b[1]*b[1]; 1390 double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); 1391 fjac(i,0) = e / b[1]; 1392 fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; 1393 fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; 1394 } 1395 return 0; 1396 } 1397 }; 1398 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}; 1399 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 }; 1400 1401 // http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml 1402 void testNistEckerle4(void) 1403 { 1404 const int n=3; 1405 int info; 1406 1407 VectorXd x(n); 1408 1409 /* 1410 * First try 1411 */ 1412 x<< 1., 10., 500.; 1413 // do the computation 1414 eckerle4_functor functor; 1415 LevenbergMarquardt<eckerle4_functor> lm(functor); 1416 info = lm.minimize(x); 1417 1418 // check return value 1419 VERIFY_IS_EQUAL(info, 1); 1420 VERIFY_IS_EQUAL(lm.nfev(), 18); 1421 VERIFY_IS_EQUAL(lm.njev(), 15); 1422 // check norm^2 1423 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); 1424 // check x 1425 VERIFY_IS_APPROX(x[0], 1.5543827178); 1426 VERIFY_IS_APPROX(x[1], 4.0888321754); 1427 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1428 1429 /* 1430 * Second try 1431 */ 1432 x<< 1.5, 5., 450.; 1433 // do the computation 1434 info = lm.minimize(x); 1435 1436 // check return value 1437 VERIFY_IS_EQUAL(info, 1); 1438 VERIFY_IS_EQUAL(lm.nfev(), 7); 1439 VERIFY_IS_EQUAL(lm.njev(), 6); 1440 // check norm^2 1441 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); 1442 // check x 1443 VERIFY_IS_APPROX(x[0], 1.5543827178); 1444 VERIFY_IS_APPROX(x[1], 4.0888321754); 1445 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1446 } 1447 1448 void test_levenberg_marquardt() 1449 { 1450 // Tests using the examples provided by (c)minpack 1451 CALL_SUBTEST(testLmder1()); 1452 CALL_SUBTEST(testLmder()); 1453 CALL_SUBTEST(testLmdif1()); 1454 // CALL_SUBTEST(testLmstr1()); 1455 // CALL_SUBTEST(testLmstr()); 1456 CALL_SUBTEST(testLmdif()); 1457 1458 // NIST tests, level of difficulty = "Lower" 1459 CALL_SUBTEST(testNistMisra1a()); 1460 CALL_SUBTEST(testNistChwirut2()); 1461 1462 // NIST tests, level of difficulty = "Average" 1463 CALL_SUBTEST(testNistHahn1()); 1464 CALL_SUBTEST(testNistMisra1d()); 1465 CALL_SUBTEST(testNistMGH17()); 1466 CALL_SUBTEST(testNistLanczos1()); 1467 1468 // // NIST tests, level of difficulty = "Higher" 1469 CALL_SUBTEST(testNistRat42()); 1470 CALL_SUBTEST(testNistMGH10()); 1471 CALL_SUBTEST(testNistBoxBOD()); 1472 // CALL_SUBTEST(testNistMGH09()); 1473 CALL_SUBTEST(testNistBennett5()); 1474 CALL_SUBTEST(testNistThurber()); 1475 CALL_SUBTEST(testNistRat43()); 1476 CALL_SUBTEST(testNistEckerle4()); 1477 } 1478
