NonLinearOptimization.cpp (64597B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> 5 6 #include <stdio.h> 7 8 #include "main.h" 9 #include <unsupported/Eigen/NonLinearOptimization> 10 11 // This disables some useless Warnings on MSVC. 12 // It is intended to be done for this test only. 13 #include <Eigen/src/Core/util/DisableStupidWarnings.h> 14 15 // tolerance for chekcing number of iterations 16 #define LM_EVAL_COUNT_TOL 4/3 17 18 #define LM_CHECK_N_ITERS(SOLVER,NFEV,NJEV) { \ 19 ++g_test_level; \ 20 VERIFY_IS_EQUAL(SOLVER.nfev, NFEV); \ 21 VERIFY_IS_EQUAL(SOLVER.njev, NJEV); \ 22 --g_test_level; \ 23 VERIFY(SOLVER.nfev <= NFEV * LM_EVAL_COUNT_TOL); \ 24 VERIFY(SOLVER.njev <= NJEV * LM_EVAL_COUNT_TOL); \ 25 } 26 27 int fcn_chkder(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag) 28 { 29 /* subroutine fcn for chkder example. */ 30 31 int i; 32 assert(15 == fvec.size()); 33 assert(3 == x.size()); 34 double tmp1, tmp2, tmp3, tmp4; 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 39 if (iflag == 0) 40 return 0; 41 42 if (iflag != 2) 43 for (i=0; i<15; i++) { 44 tmp1 = i+1; 45 tmp2 = 16-i-1; 46 tmp3 = tmp1; 47 if (i >= 8) tmp3 = tmp2; 48 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 49 } 50 else { 51 for (i = 0; i < 15; i++) { 52 tmp1 = i+1; 53 tmp2 = 16-i-1; 54 55 /* error introduced into next statement for illustration. */ 56 /* corrected statement should read tmp3 = tmp1 . */ 57 58 tmp3 = tmp2; 59 if (i >= 8) tmp3 = tmp2; 60 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4=tmp4*tmp4; 61 fjac(i,0) = -1.; 62 fjac(i,1) = tmp1*tmp2/tmp4; 63 fjac(i,2) = tmp1*tmp3/tmp4; 64 } 65 } 66 return 0; 67 } 68 69 70 void testChkder() 71 { 72 const int m=15, n=3; 73 VectorXd x(n), fvec(m), xp, fvecp(m), err; 74 MatrixXd fjac(m,n); 75 VectorXi ipvt; 76 77 /* the following values should be suitable for */ 78 /* checking the jacobian matrix. */ 79 x << 9.2e-1, 1.3e-1, 5.4e-1; 80 81 internal::chkder(x, fvec, fjac, xp, fvecp, 1, err); 82 fcn_chkder(x, fvec, fjac, 1); 83 fcn_chkder(x, fvec, fjac, 2); 84 fcn_chkder(xp, fvecp, fjac, 1); 85 internal::chkder(x, fvec, fjac, xp, fvecp, 2, err); 86 87 fvecp -= fvec; 88 89 // check those 90 VectorXd fvec_ref(m), fvecp_ref(m), err_ref(m); 91 fvec_ref << 92 -1.181606, -1.429655, -1.606344, 93 -1.745269, -1.840654, -1.921586, 94 -1.984141, -2.022537, -2.468977, 95 -2.827562, -3.473582, -4.437612, 96 -6.047662, -9.267761, -18.91806; 97 fvecp_ref << 98 -7.724666e-09, -3.432406e-09, -2.034843e-10, 99 2.313685e-09, 4.331078e-09, 5.984096e-09, 100 7.363281e-09, 8.53147e-09, 1.488591e-08, 101 2.33585e-08, 3.522012e-08, 5.301255e-08, 102 8.26666e-08, 1.419747e-07, 3.19899e-07; 103 err_ref << 104 0.1141397, 0.09943516, 0.09674474, 105 0.09980447, 0.1073116, 0.1220445, 106 0.1526814, 1, 1, 107 1, 1, 1, 108 1, 1, 1; 109 110 VERIFY_IS_APPROX(fvec, fvec_ref); 111 VERIFY_IS_APPROX(fvecp, fvecp_ref); 112 VERIFY_IS_APPROX(err, err_ref); 113 } 114 115 // Generic functor 116 template<typename _Scalar, int NX=Dynamic, int NY=Dynamic> 117 struct Functor 118 { 119 typedef _Scalar Scalar; 120 enum { 121 InputsAtCompileTime = NX, 122 ValuesAtCompileTime = NY 123 }; 124 typedef Matrix<Scalar,InputsAtCompileTime,1> InputType; 125 typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType; 126 typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; 127 128 const int m_inputs, m_values; 129 130 Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} 131 Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {} 132 133 int inputs() const { return m_inputs; } 134 int values() const { return m_values; } 135 136 // you should define that in the subclass : 137 // void operator() (const InputType& x, ValueType* v, JacobianType* _j=0) const; 138 }; 139 140 struct lmder_functor : Functor<double> 141 { 142 lmder_functor(void): Functor<double>(3,15) {} 143 int operator()(const VectorXd &x, VectorXd &fvec) const 144 { 145 double tmp1, tmp2, tmp3; 146 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, 147 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 148 149 for (int i = 0; i < values(); i++) 150 { 151 tmp1 = i+1; 152 tmp2 = 16 - i - 1; 153 tmp3 = (i>=8)? tmp2 : tmp1; 154 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 155 } 156 return 0; 157 } 158 159 int df(const VectorXd &x, MatrixXd &fjac) const 160 { 161 double tmp1, tmp2, tmp3, tmp4; 162 for (int i = 0; i < values(); i++) 163 { 164 tmp1 = i+1; 165 tmp2 = 16 - i - 1; 166 tmp3 = (i>=8)? tmp2 : tmp1; 167 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 168 fjac(i,0) = -1; 169 fjac(i,1) = tmp1*tmp2/tmp4; 170 fjac(i,2) = tmp1*tmp3/tmp4; 171 } 172 return 0; 173 } 174 }; 175 176 void testLmder1() 177 { 178 int n=3, info; 179 180 VectorXd x; 181 182 /* the following starting values provide a rough fit. */ 183 x.setConstant(n, 1.); 184 185 // do the computation 186 lmder_functor functor; 187 LevenbergMarquardt<lmder_functor> lm(functor); 188 info = lm.lmder1(x); 189 190 // check return value 191 VERIFY_IS_EQUAL(info, 1); 192 LM_CHECK_N_ITERS(lm, 6, 5); 193 194 // check norm 195 VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 196 197 // check x 198 VectorXd x_ref(n); 199 x_ref << 0.08241058, 1.133037, 2.343695; 200 VERIFY_IS_APPROX(x, x_ref); 201 } 202 203 void testLmder() 204 { 205 const int m=15, n=3; 206 int info; 207 double fnorm, covfac; 208 VectorXd x; 209 210 /* the following starting values provide a rough fit. */ 211 x.setConstant(n, 1.); 212 213 // do the computation 214 lmder_functor functor; 215 LevenbergMarquardt<lmder_functor> lm(functor); 216 info = lm.minimize(x); 217 218 // check return values 219 VERIFY_IS_EQUAL(info, 1); 220 LM_CHECK_N_ITERS(lm, 6, 5); 221 222 // check norm 223 fnorm = lm.fvec.blueNorm(); 224 VERIFY_IS_APPROX(fnorm, 0.09063596); 225 226 // check x 227 VectorXd x_ref(n); 228 x_ref << 0.08241058, 1.133037, 2.343695; 229 VERIFY_IS_APPROX(x, x_ref); 230 231 // check covariance 232 covfac = fnorm*fnorm/(m-n); 233 internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 234 235 MatrixXd cov_ref(n,n); 236 cov_ref << 237 0.0001531202, 0.002869941, -0.002656662, 238 0.002869941, 0.09480935, -0.09098995, 239 -0.002656662, -0.09098995, 0.08778727; 240 241 // std::cout << fjac*covfac << std::endl; 242 243 MatrixXd cov; 244 cov = covfac*lm.fjac.topLeftCorner<n,n>(); 245 VERIFY_IS_APPROX( cov, cov_ref); 246 // TODO: why isn't this allowed ? : 247 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 248 } 249 250 struct hybrj_functor : Functor<double> 251 { 252 hybrj_functor(void) : Functor<double>(9,9) {} 253 254 int operator()(const VectorXd &x, VectorXd &fvec) 255 { 256 double temp, temp1, temp2; 257 const VectorXd::Index n = x.size(); 258 assert(fvec.size()==n); 259 for (VectorXd::Index k = 0; k < n; k++) 260 { 261 temp = (3. - 2.*x[k])*x[k]; 262 temp1 = 0.; 263 if (k) temp1 = x[k-1]; 264 temp2 = 0.; 265 if (k != n-1) temp2 = x[k+1]; 266 fvec[k] = temp - temp1 - 2.*temp2 + 1.; 267 } 268 return 0; 269 } 270 int df(const VectorXd &x, MatrixXd &fjac) 271 { 272 const VectorXd::Index n = x.size(); 273 assert(fjac.rows()==n); 274 assert(fjac.cols()==n); 275 for (VectorXd::Index k = 0; k < n; k++) 276 { 277 for (VectorXd::Index j = 0; j < n; j++) 278 fjac(k,j) = 0.; 279 fjac(k,k) = 3.- 4.*x[k]; 280 if (k) fjac(k,k-1) = -1.; 281 if (k != n-1) fjac(k,k+1) = -2.; 282 } 283 return 0; 284 } 285 }; 286 287 288 void testHybrj1() 289 { 290 const int n=9; 291 int info; 292 VectorXd x(n); 293 294 /* the following starting values provide a rough fit. */ 295 x.setConstant(n, -1.); 296 297 // do the computation 298 hybrj_functor functor; 299 HybridNonLinearSolver<hybrj_functor> solver(functor); 300 info = solver.hybrj1(x); 301 302 // check return value 303 VERIFY_IS_EQUAL(info, 1); 304 LM_CHECK_N_ITERS(solver, 11, 1); 305 306 // check norm 307 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 308 309 310 // check x 311 VectorXd x_ref(n); 312 x_ref << 313 -0.5706545, -0.6816283, -0.7017325, 314 -0.7042129, -0.701369, -0.6918656, 315 -0.665792, -0.5960342, -0.4164121; 316 VERIFY_IS_APPROX(x, x_ref); 317 } 318 319 void testHybrj() 320 { 321 const int n=9; 322 int info; 323 VectorXd x(n); 324 325 /* the following starting values provide a rough fit. */ 326 x.setConstant(n, -1.); 327 328 329 // do the computation 330 hybrj_functor functor; 331 HybridNonLinearSolver<hybrj_functor> solver(functor); 332 solver.diag.setConstant(n, 1.); 333 solver.useExternalScaling = true; 334 info = solver.solve(x); 335 336 // check return value 337 VERIFY_IS_EQUAL(info, 1); 338 LM_CHECK_N_ITERS(solver, 11, 1); 339 340 // check norm 341 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 342 343 344 // check x 345 VectorXd x_ref(n); 346 x_ref << 347 -0.5706545, -0.6816283, -0.7017325, 348 -0.7042129, -0.701369, -0.6918656, 349 -0.665792, -0.5960342, -0.4164121; 350 VERIFY_IS_APPROX(x, x_ref); 351 352 } 353 354 struct hybrd_functor : Functor<double> 355 { 356 hybrd_functor(void) : Functor<double>(9,9) {} 357 int operator()(const VectorXd &x, VectorXd &fvec) const 358 { 359 double temp, temp1, temp2; 360 const VectorXd::Index n = x.size(); 361 362 assert(fvec.size()==n); 363 for (VectorXd::Index k=0; k < n; k++) 364 { 365 temp = (3. - 2.*x[k])*x[k]; 366 temp1 = 0.; 367 if (k) temp1 = x[k-1]; 368 temp2 = 0.; 369 if (k != n-1) temp2 = x[k+1]; 370 fvec[k] = temp - temp1 - 2.*temp2 + 1.; 371 } 372 return 0; 373 } 374 }; 375 376 void testHybrd1() 377 { 378 int n=9, info; 379 VectorXd x(n); 380 381 /* the following starting values provide a rough solution. */ 382 x.setConstant(n, -1.); 383 384 // do the computation 385 hybrd_functor functor; 386 HybridNonLinearSolver<hybrd_functor> solver(functor); 387 info = solver.hybrd1(x); 388 389 // check return value 390 VERIFY_IS_EQUAL(info, 1); 391 VERIFY_IS_EQUAL(solver.nfev, 20); 392 393 // check norm 394 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 395 396 // check x 397 VectorXd x_ref(n); 398 x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121; 399 VERIFY_IS_APPROX(x, x_ref); 400 } 401 402 void testHybrd() 403 { 404 const int n=9; 405 int info; 406 VectorXd x; 407 408 /* the following starting values provide a rough fit. */ 409 x.setConstant(n, -1.); 410 411 // do the computation 412 hybrd_functor functor; 413 HybridNonLinearSolver<hybrd_functor> solver(functor); 414 solver.parameters.nb_of_subdiagonals = 1; 415 solver.parameters.nb_of_superdiagonals = 1; 416 solver.diag.setConstant(n, 1.); 417 solver.useExternalScaling = true; 418 info = solver.solveNumericalDiff(x); 419 420 // check return value 421 VERIFY_IS_EQUAL(info, 1); 422 VERIFY_IS_EQUAL(solver.nfev, 14); 423 424 // check norm 425 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 426 427 // check x 428 VectorXd x_ref(n); 429 x_ref << 430 -0.5706545, -0.6816283, -0.7017325, 431 -0.7042129, -0.701369, -0.6918656, 432 -0.665792, -0.5960342, -0.4164121; 433 VERIFY_IS_APPROX(x, x_ref); 434 } 435 436 struct lmstr_functor : Functor<double> 437 { 438 lmstr_functor(void) : Functor<double>(3,15) {} 439 int operator()(const VectorXd &x, VectorXd &fvec) 440 { 441 /* subroutine fcn for lmstr1 example. */ 442 double tmp1, tmp2, tmp3; 443 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, 444 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 445 446 assert(15==fvec.size()); 447 assert(3==x.size()); 448 449 for (int i=0; i<15; i++) 450 { 451 tmp1 = i+1; 452 tmp2 = 16 - i - 1; 453 tmp3 = (i>=8)? tmp2 : tmp1; 454 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 455 } 456 return 0; 457 } 458 int df(const VectorXd &x, VectorXd &jac_row, VectorXd::Index rownb) 459 { 460 assert(x.size()==3); 461 assert(jac_row.size()==x.size()); 462 double tmp1, tmp2, tmp3, tmp4; 463 464 VectorXd::Index i = rownb-2; 465 tmp1 = i+1; 466 tmp2 = 16 - i - 1; 467 tmp3 = (i>=8)? tmp2 : tmp1; 468 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 469 jac_row[0] = -1; 470 jac_row[1] = tmp1*tmp2/tmp4; 471 jac_row[2] = tmp1*tmp3/tmp4; 472 return 0; 473 } 474 }; 475 476 void testLmstr1() 477 { 478 const int n=3; 479 int info; 480 481 VectorXd x(n); 482 483 /* the following starting values provide a rough fit. */ 484 x.setConstant(n, 1.); 485 486 // do the computation 487 lmstr_functor functor; 488 LevenbergMarquardt<lmstr_functor> lm(functor); 489 info = lm.lmstr1(x); 490 491 // check return value 492 VERIFY_IS_EQUAL(info, 1); 493 LM_CHECK_N_ITERS(lm, 6, 5); 494 495 // check norm 496 VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 497 498 // check x 499 VectorXd x_ref(n); 500 x_ref << 0.08241058, 1.133037, 2.343695 ; 501 VERIFY_IS_APPROX(x, x_ref); 502 } 503 504 void testLmstr() 505 { 506 const int n=3; 507 int info; 508 double fnorm; 509 VectorXd x(n); 510 511 /* the following starting values provide a rough fit. */ 512 x.setConstant(n, 1.); 513 514 // do the computation 515 lmstr_functor functor; 516 LevenbergMarquardt<lmstr_functor> lm(functor); 517 info = lm.minimizeOptimumStorage(x); 518 519 // check return values 520 VERIFY_IS_EQUAL(info, 1); 521 LM_CHECK_N_ITERS(lm, 6, 5); 522 523 // check norm 524 fnorm = lm.fvec.blueNorm(); 525 VERIFY_IS_APPROX(fnorm, 0.09063596); 526 527 // check x 528 VectorXd x_ref(n); 529 x_ref << 0.08241058, 1.133037, 2.343695; 530 VERIFY_IS_APPROX(x, x_ref); 531 532 } 533 534 struct lmdif_functor : Functor<double> 535 { 536 lmdif_functor(void) : Functor<double>(3,15) {} 537 int operator()(const VectorXd &x, VectorXd &fvec) const 538 { 539 int i; 540 double tmp1,tmp2,tmp3; 541 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, 542 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; 543 544 assert(x.size()==3); 545 assert(fvec.size()==15); 546 for (i=0; i<15; i++) 547 { 548 tmp1 = i+1; 549 tmp2 = 15 - i; 550 tmp3 = tmp1; 551 552 if (i >= 8) tmp3 = tmp2; 553 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 554 } 555 return 0; 556 } 557 }; 558 559 void testLmdif1() 560 { 561 const int n=3; 562 int info; 563 564 VectorXd x(n), fvec(15); 565 566 /* the following starting values provide a rough fit. */ 567 x.setConstant(n, 1.); 568 569 // do the computation 570 lmdif_functor functor; 571 DenseIndex nfev = -1; // initialize to avoid maybe-uninitialized warning 572 info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev); 573 574 // check return value 575 VERIFY_IS_EQUAL(info, 1); 576 VERIFY_IS_EQUAL(nfev, 26); 577 578 // check norm 579 functor(x, fvec); 580 VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); 581 582 // check x 583 VectorXd x_ref(n); 584 x_ref << 0.0824106, 1.1330366, 2.3436947; 585 VERIFY_IS_APPROX(x, x_ref); 586 587 } 588 589 void testLmdif() 590 { 591 const int m=15, n=3; 592 int info; 593 double fnorm, covfac; 594 VectorXd x(n); 595 596 /* the following starting values provide a rough fit. */ 597 x.setConstant(n, 1.); 598 599 // do the computation 600 lmdif_functor functor; 601 NumericalDiff<lmdif_functor> numDiff(functor); 602 LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff); 603 info = lm.minimize(x); 604 605 // check return values 606 VERIFY_IS_EQUAL(info, 1); 607 VERIFY_IS_EQUAL(lm.nfev, 26); 608 609 // check norm 610 fnorm = lm.fvec.blueNorm(); 611 VERIFY_IS_APPROX(fnorm, 0.09063596); 612 613 // check x 614 VectorXd x_ref(n); 615 x_ref << 0.08241058, 1.133037, 2.343695; 616 VERIFY_IS_APPROX(x, x_ref); 617 618 // check covariance 619 covfac = fnorm*fnorm/(m-n); 620 internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 621 622 MatrixXd cov_ref(n,n); 623 cov_ref << 624 0.0001531202, 0.002869942, -0.002656662, 625 0.002869942, 0.09480937, -0.09098997, 626 -0.002656662, -0.09098997, 0.08778729; 627 628 // std::cout << fjac*covfac << std::endl; 629 630 MatrixXd cov; 631 cov = covfac*lm.fjac.topLeftCorner<n,n>(); 632 VERIFY_IS_APPROX( cov, cov_ref); 633 // TODO: why isn't this allowed ? : 634 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 635 } 636 637 struct chwirut2_functor : Functor<double> 638 { 639 chwirut2_functor(void) : Functor<double>(3,54) {} 640 static const double m_x[54]; 641 static const double m_y[54]; 642 int operator()(const VectorXd &b, VectorXd &fvec) 643 { 644 int i; 645 646 assert(b.size()==3); 647 assert(fvec.size()==54); 648 for(i=0; i<54; i++) { 649 double x = m_x[i]; 650 fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; 651 } 652 return 0; 653 } 654 int df(const VectorXd &b, MatrixXd &fjac) 655 { 656 assert(b.size()==3); 657 assert(fjac.rows()==54); 658 assert(fjac.cols()==3); 659 for(int i=0; i<54; i++) { 660 double x = m_x[i]; 661 double factor = 1./(b[1]+b[2]*x); 662 double e = exp(-b[0]*x); 663 fjac(i,0) = -x*e*factor; 664 fjac(i,1) = -e*factor*factor; 665 fjac(i,2) = -x*e*factor*factor; 666 } 667 return 0; 668 } 669 }; 670 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}; 671 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 }; 672 673 // http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml 674 void testNistChwirut2(void) 675 { 676 const int n=3; 677 int info; 678 679 VectorXd x(n); 680 681 /* 682 * First try 683 */ 684 x<< 0.1, 0.01, 0.02; 685 // do the computation 686 chwirut2_functor functor; 687 LevenbergMarquardt<chwirut2_functor> lm(functor); 688 info = lm.minimize(x); 689 690 // check return value 691 VERIFY_IS_EQUAL(info, 1); 692 LM_CHECK_N_ITERS(lm, 10, 8); 693 // check norm^2 694 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 695 // check x 696 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 697 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 698 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 699 700 /* 701 * Second try 702 */ 703 x<< 0.15, 0.008, 0.010; 704 // do the computation 705 lm.resetParameters(); 706 lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 707 lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 708 info = lm.minimize(x); 709 710 // check return value 711 VERIFY_IS_EQUAL(info, 1); 712 LM_CHECK_N_ITERS(lm, 7, 6); 713 // check norm^2 714 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 715 // check x 716 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 717 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 718 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 719 } 720 721 722 struct misra1a_functor : Functor<double> 723 { 724 misra1a_functor(void) : Functor<double>(2,14) {} 725 static const double m_x[14]; 726 static const double m_y[14]; 727 int operator()(const VectorXd &b, VectorXd &fvec) 728 { 729 assert(b.size()==2); 730 assert(fvec.size()==14); 731 for(int i=0; i<14; i++) { 732 fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; 733 } 734 return 0; 735 } 736 int df(const VectorXd &b, MatrixXd &fjac) 737 { 738 assert(b.size()==2); 739 assert(fjac.rows()==14); 740 assert(fjac.cols()==2); 741 for(int i=0; i<14; i++) { 742 fjac(i,0) = (1.-exp(-b[1]*m_x[i])); 743 fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); 744 } 745 return 0; 746 } 747 }; 748 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}; 749 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}; 750 751 // http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml 752 void testNistMisra1a(void) 753 { 754 const int n=2; 755 int info; 756 757 VectorXd x(n); 758 759 /* 760 * First try 761 */ 762 x<< 500., 0.0001; 763 // do the computation 764 misra1a_functor functor; 765 LevenbergMarquardt<misra1a_functor> lm(functor); 766 info = lm.minimize(x); 767 768 // check return value 769 VERIFY_IS_EQUAL(info, 1); 770 LM_CHECK_N_ITERS(lm, 19, 15); 771 // check norm^2 772 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 773 // check x 774 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 775 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 776 777 /* 778 * Second try 779 */ 780 x<< 250., 0.0005; 781 // do the computation 782 info = lm.minimize(x); 783 784 // check return value 785 VERIFY_IS_EQUAL(info, 1); 786 LM_CHECK_N_ITERS(lm, 5, 4); 787 // check norm^2 788 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 789 // check x 790 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 791 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 792 } 793 794 struct hahn1_functor : Functor<double> 795 { 796 hahn1_functor(void) : Functor<double>(7,236) {} 797 static const double m_x[236]; 798 int operator()(const VectorXd &b, VectorXd &fvec) 799 { 800 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 , 801 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 , 802 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 }; 803 804 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 805 806 assert(b.size()==7); 807 assert(fvec.size()==236); 808 for(int i=0; i<236; i++) { 809 double x=m_x[i], xx=x*x, xxx=xx*x; 810 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]; 811 } 812 return 0; 813 } 814 815 int df(const VectorXd &b, MatrixXd &fjac) 816 { 817 assert(b.size()==7); 818 assert(fjac.rows()==236); 819 assert(fjac.cols()==7); 820 for(int i=0; i<236; i++) { 821 double x=m_x[i], xx=x*x, xxx=xx*x; 822 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 823 fjac(i,0) = 1.*fact; 824 fjac(i,1) = x*fact; 825 fjac(i,2) = xx*fact; 826 fjac(i,3) = xxx*fact; 827 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 828 fjac(i,4) = x*fact; 829 fjac(i,5) = xx*fact; 830 fjac(i,6) = xxx*fact; 831 } 832 return 0; 833 } 834 }; 835 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 , 836 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 , 837 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}; 838 839 // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml 840 void testNistHahn1(void) 841 { 842 const int n=7; 843 int info; 844 845 VectorXd x(n); 846 847 /* 848 * First try 849 */ 850 x<< 10., -1., .05, -.00001, -.05, .001, -.000001; 851 // do the computation 852 hahn1_functor functor; 853 LevenbergMarquardt<hahn1_functor> lm(functor); 854 info = lm.minimize(x); 855 856 // check return value 857 VERIFY_IS_EQUAL(info, 1); 858 LM_CHECK_N_ITERS(lm, 11, 10); 859 // check norm^2 860 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 861 // check x 862 VERIFY_IS_APPROX(x[0], 1.0776351733E+00); 863 VERIFY_IS_APPROX(x[1],-1.2269296921E-01); 864 VERIFY_IS_APPROX(x[2], 4.0863750610E-03); 865 VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 866 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 867 VERIFY_IS_APPROX(x[5], 2.4053735503E-04); 868 VERIFY_IS_APPROX(x[6],-1.2314450199E-07); 869 870 /* 871 * Second try 872 */ 873 x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; 874 // do the computation 875 info = lm.minimize(x); 876 877 // check return value 878 VERIFY_IS_EQUAL(info, 1); 879 LM_CHECK_N_ITERS(lm, 11, 10); 880 // check norm^2 881 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 882 // check x 883 VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 884 VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 885 VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 886 VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 887 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 888 VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 889 VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 890 891 } 892 893 struct misra1d_functor : Functor<double> 894 { 895 misra1d_functor(void) : Functor<double>(2,14) {} 896 static const double x[14]; 897 static const double y[14]; 898 int operator()(const VectorXd &b, VectorXd &fvec) 899 { 900 assert(b.size()==2); 901 assert(fvec.size()==14); 902 for(int i=0; i<14; i++) { 903 fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; 904 } 905 return 0; 906 } 907 int df(const VectorXd &b, MatrixXd &fjac) 908 { 909 assert(b.size()==2); 910 assert(fjac.rows()==14); 911 assert(fjac.cols()==2); 912 for(int i=0; i<14; i++) { 913 double den = 1.+b[1]*x[i]; 914 fjac(i,0) = b[1]*x[i] / den; 915 fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; 916 } 917 return 0; 918 } 919 }; 920 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}; 921 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}; 922 923 // http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml 924 void testNistMisra1d(void) 925 { 926 const int n=2; 927 int info; 928 929 VectorXd x(n); 930 931 /* 932 * First try 933 */ 934 x<< 500., 0.0001; 935 // do the computation 936 misra1d_functor functor; 937 LevenbergMarquardt<misra1d_functor> lm(functor); 938 info = lm.minimize(x); 939 940 // check return value 941 VERIFY_IS_EQUAL(info, 3); 942 LM_CHECK_N_ITERS(lm, 9, 7); 943 // check norm^2 944 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 945 // check x 946 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 947 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 948 949 /* 950 * Second try 951 */ 952 x<< 450., 0.0003; 953 // do the computation 954 info = lm.minimize(x); 955 956 // check return value 957 VERIFY_IS_EQUAL(info, 1); 958 LM_CHECK_N_ITERS(lm, 4, 3); 959 // check norm^2 960 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 961 // check x 962 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 963 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 964 } 965 966 967 struct lanczos1_functor : Functor<double> 968 { 969 lanczos1_functor(void) : Functor<double>(6,24) {} 970 static const double x[24]; 971 static const double y[24]; 972 int operator()(const VectorXd &b, VectorXd &fvec) 973 { 974 assert(b.size()==6); 975 assert(fvec.size()==24); 976 for(int i=0; i<24; i++) 977 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]; 978 return 0; 979 } 980 int df(const VectorXd &b, MatrixXd &fjac) 981 { 982 assert(b.size()==6); 983 assert(fjac.rows()==24); 984 assert(fjac.cols()==6); 985 for(int i=0; i<24; i++) { 986 fjac(i,0) = exp(-b[1]*x[i]); 987 fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); 988 fjac(i,2) = exp(-b[3]*x[i]); 989 fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); 990 fjac(i,4) = exp(-b[5]*x[i]); 991 fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); 992 } 993 return 0; 994 } 995 }; 996 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 }; 997 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 }; 998 999 // http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml 1000 void testNistLanczos1(void) 1001 { 1002 const int n=6; 1003 int info; 1004 1005 VectorXd x(n); 1006 1007 /* 1008 * First try 1009 */ 1010 x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; 1011 // do the computation 1012 lanczos1_functor functor; 1013 LevenbergMarquardt<lanczos1_functor> lm(functor); 1014 info = lm.minimize(x); 1015 1016 // check return value 1017 VERIFY_IS_EQUAL(info, 2); 1018 LM_CHECK_N_ITERS(lm, 79, 72); 1019 // check norm^2 1020 std::cout.precision(30); 1021 std::cout << lm.fvec.squaredNorm() << "\n"; 1022 VERIFY(lm.fvec.squaredNorm() <= 1.4307867721E-25); 1023 // check x 1024 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1025 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1026 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1027 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1028 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1029 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1030 1031 /* 1032 * Second try 1033 */ 1034 x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; 1035 // do the computation 1036 info = lm.minimize(x); 1037 1038 // check return value 1039 VERIFY_IS_EQUAL(info, 2); 1040 LM_CHECK_N_ITERS(lm, 9, 8); 1041 // check norm^2 1042 VERIFY(lm.fvec.squaredNorm() <= 1.4307867721E-25); 1043 // check x 1044 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1045 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1046 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1047 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1048 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1049 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1050 1051 } 1052 1053 struct rat42_functor : Functor<double> 1054 { 1055 rat42_functor(void) : Functor<double>(3,9) {} 1056 static const double x[9]; 1057 static const double y[9]; 1058 int operator()(const VectorXd &b, VectorXd &fvec) 1059 { 1060 assert(b.size()==3); 1061 assert(fvec.size()==9); 1062 for(int i=0; i<9; i++) { 1063 fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; 1064 } 1065 return 0; 1066 } 1067 1068 int df(const VectorXd &b, MatrixXd &fjac) 1069 { 1070 assert(b.size()==3); 1071 assert(fjac.rows()==9); 1072 assert(fjac.cols()==3); 1073 for(int i=0; i<9; i++) { 1074 double e = exp(b[1]-b[2]*x[i]); 1075 fjac(i,0) = 1./(1.+e); 1076 fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); 1077 fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); 1078 } 1079 return 0; 1080 } 1081 }; 1082 const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; 1083 const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; 1084 1085 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml 1086 void testNistRat42(void) 1087 { 1088 const int n=3; 1089 int info; 1090 1091 VectorXd x(n); 1092 1093 /* 1094 * First try 1095 */ 1096 x<< 100., 1., 0.1; 1097 // do the computation 1098 rat42_functor functor; 1099 LevenbergMarquardt<rat42_functor> lm(functor); 1100 info = lm.minimize(x); 1101 1102 // check return value 1103 VERIFY_IS_EQUAL(info, 1); 1104 LM_CHECK_N_ITERS(lm, 10, 8); 1105 // check norm^2 1106 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00); 1107 // check x 1108 VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 1109 VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 1110 VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 1111 1112 /* 1113 * Second try 1114 */ 1115 x<< 75., 2.5, 0.07; 1116 // do the computation 1117 info = lm.minimize(x); 1118 1119 // check return value 1120 VERIFY_IS_EQUAL(info, 1); 1121 LM_CHECK_N_ITERS(lm, 6, 5); 1122 // check norm^2 1123 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00); 1124 // check x 1125 VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 1126 VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 1127 VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 1128 } 1129 1130 struct MGH10_functor : Functor<double> 1131 { 1132 MGH10_functor(void) : Functor<double>(3,16) {} 1133 static const double x[16]; 1134 static const double y[16]; 1135 int operator()(const VectorXd &b, VectorXd &fvec) 1136 { 1137 assert(b.size()==3); 1138 assert(fvec.size()==16); 1139 for(int i=0; i<16; i++) 1140 fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; 1141 return 0; 1142 } 1143 int df(const VectorXd &b, MatrixXd &fjac) 1144 { 1145 assert(b.size()==3); 1146 assert(fjac.rows()==16); 1147 assert(fjac.cols()==3); 1148 for(int i=0; i<16; i++) { 1149 double factor = 1./(x[i]+b[2]); 1150 double e = exp(b[1]*factor); 1151 fjac(i,0) = e; 1152 fjac(i,1) = b[0]*factor*e; 1153 fjac(i,2) = -b[1]*b[0]*factor*factor*e; 1154 } 1155 return 0; 1156 } 1157 }; 1158 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 }; 1159 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 }; 1160 1161 // http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml 1162 void testNistMGH10(void) 1163 { 1164 const int n=3; 1165 int info; 1166 1167 VectorXd x(n); 1168 1169 /* 1170 * First try 1171 */ 1172 x<< 2., 400000., 25000.; 1173 // do the computation 1174 MGH10_functor functor; 1175 LevenbergMarquardt<MGH10_functor> lm(functor); 1176 info = lm.minimize(x); 1177 1178 // check return value 1179 VERIFY_IS_EQUAL(info, 2); 1180 LM_CHECK_N_ITERS(lm, 284, 249); 1181 // check norm^2 1182 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1183 // check x 1184 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1185 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1186 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1187 1188 /* 1189 * Second try 1190 */ 1191 x<< 0.02, 4000., 250.; 1192 // do the computation 1193 info = lm.minimize(x); 1194 1195 // check return value 1196 VERIFY_IS_EQUAL(info, 3); 1197 LM_CHECK_N_ITERS(lm, 126, 116); 1198 // check norm^2 1199 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1200 // check x 1201 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1202 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1203 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1204 } 1205 1206 1207 struct BoxBOD_functor : Functor<double> 1208 { 1209 BoxBOD_functor(void) : Functor<double>(2,6) {} 1210 static const double x[6]; 1211 int operator()(const VectorXd &b, VectorXd &fvec) 1212 { 1213 static const double y[6] = { 109., 149., 149., 191., 213., 224. }; 1214 assert(b.size()==2); 1215 assert(fvec.size()==6); 1216 for(int i=0; i<6; i++) 1217 fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; 1218 return 0; 1219 } 1220 int df(const VectorXd &b, MatrixXd &fjac) 1221 { 1222 assert(b.size()==2); 1223 assert(fjac.rows()==6); 1224 assert(fjac.cols()==2); 1225 for(int i=0; i<6; i++) { 1226 double e = exp(-b[1]*x[i]); 1227 fjac(i,0) = 1.-e; 1228 fjac(i,1) = b[0]*x[i]*e; 1229 } 1230 return 0; 1231 } 1232 }; 1233 const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; 1234 1235 // http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml 1236 void testNistBoxBOD(void) 1237 { 1238 const int n=2; 1239 int info; 1240 1241 VectorXd x(n); 1242 1243 /* 1244 * First try 1245 */ 1246 x<< 1., 1.; 1247 // do the computation 1248 BoxBOD_functor functor; 1249 LevenbergMarquardt<BoxBOD_functor> lm(functor); 1250 lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1251 lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1252 lm.parameters.factor = 10.; 1253 info = lm.minimize(x); 1254 1255 // check return value 1256 VERIFY_IS_EQUAL(info, 1); 1257 LM_CHECK_N_ITERS(lm, 31, 25); 1258 // check norm^2 1259 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1260 // check x 1261 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1262 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1263 1264 /* 1265 * Second try 1266 */ 1267 x<< 100., 0.75; 1268 // do the computation 1269 lm.resetParameters(); 1270 lm.parameters.ftol = NumTraits<double>::epsilon(); 1271 lm.parameters.xtol = NumTraits<double>::epsilon(); 1272 info = lm.minimize(x); 1273 1274 // check return value 1275 VERIFY_IS_EQUAL(info, 1); 1276 LM_CHECK_N_ITERS(lm, 15, 14); 1277 // check norm^2 1278 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1279 // check x 1280 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1281 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1282 } 1283 1284 struct MGH17_functor : Functor<double> 1285 { 1286 MGH17_functor(void) : Functor<double>(5,33) {} 1287 static const double x[33]; 1288 static const double y[33]; 1289 int operator()(const VectorXd &b, VectorXd &fvec) 1290 { 1291 assert(b.size()==5); 1292 assert(fvec.size()==33); 1293 for(int i=0; i<33; i++) 1294 fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; 1295 return 0; 1296 } 1297 int df(const VectorXd &b, MatrixXd &fjac) 1298 { 1299 assert(b.size()==5); 1300 assert(fjac.rows()==33); 1301 assert(fjac.cols()==5); 1302 for(int i=0; i<33; i++) { 1303 fjac(i,0) = 1.; 1304 fjac(i,1) = exp(-b[3]*x[i]); 1305 fjac(i,2) = exp(-b[4]*x[i]); 1306 fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); 1307 fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); 1308 } 1309 return 0; 1310 } 1311 }; 1312 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 }; 1313 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 }; 1314 1315 // http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml 1316 void testNistMGH17(void) 1317 { 1318 const int n=5; 1319 int info; 1320 1321 VectorXd x(n); 1322 1323 /* 1324 * First try 1325 */ 1326 x<< 50., 150., -100., 1., 2.; 1327 // do the computation 1328 MGH17_functor functor; 1329 LevenbergMarquardt<MGH17_functor> lm(functor); 1330 lm.parameters.ftol = NumTraits<double>::epsilon(); 1331 lm.parameters.xtol = NumTraits<double>::epsilon(); 1332 lm.parameters.maxfev = 1000; 1333 info = lm.minimize(x); 1334 1335 // check norm^2 1336 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1337 // check x 1338 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1339 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1340 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1341 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1342 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1343 1344 // check return value 1345 VERIFY_IS_EQUAL(info, 2); 1346 LM_CHECK_N_ITERS(lm, 602, 545); 1347 1348 /* 1349 * Second try 1350 */ 1351 x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; 1352 // do the computation 1353 lm.resetParameters(); 1354 info = lm.minimize(x); 1355 1356 // check return value 1357 VERIFY_IS_EQUAL(info, 1); 1358 LM_CHECK_N_ITERS(lm, 18, 15); 1359 // check norm^2 1360 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1361 // check x 1362 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1363 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1364 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1365 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1366 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1367 } 1368 1369 struct MGH09_functor : Functor<double> 1370 { 1371 MGH09_functor(void) : Functor<double>(4,11) {} 1372 static const double _x[11]; 1373 static const double y[11]; 1374 int operator()(const VectorXd &b, VectorXd &fvec) 1375 { 1376 assert(b.size()==4); 1377 assert(fvec.size()==11); 1378 for(int i=0; i<11; i++) { 1379 double x = _x[i], xx=x*x; 1380 fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; 1381 } 1382 return 0; 1383 } 1384 int df(const VectorXd &b, MatrixXd &fjac) 1385 { 1386 assert(b.size()==4); 1387 assert(fjac.rows()==11); 1388 assert(fjac.cols()==4); 1389 for(int i=0; i<11; i++) { 1390 double x = _x[i], xx=x*x; 1391 double factor = 1./(xx+x*b[2]+b[3]); 1392 fjac(i,0) = (xx+x*b[1]) * factor; 1393 fjac(i,1) = b[0]*x* factor; 1394 fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; 1395 fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; 1396 } 1397 return 0; 1398 } 1399 }; 1400 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 }; 1401 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 }; 1402 1403 // http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml 1404 void testNistMGH09(void) 1405 { 1406 const int n=4; 1407 int info; 1408 1409 VectorXd x(n); 1410 1411 /* 1412 * First try 1413 */ 1414 x<< 25., 39, 41.5, 39.; 1415 // do the computation 1416 MGH09_functor functor; 1417 LevenbergMarquardt<MGH09_functor> lm(functor); 1418 lm.parameters.maxfev = 1000; 1419 info = lm.minimize(x); 1420 1421 // check return value 1422 VERIFY_IS_EQUAL(info, 1); 1423 LM_CHECK_N_ITERS(lm, 490, 376); 1424 // check norm^2 1425 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1426 // check x 1427 VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 1428 VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 1429 VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 1430 VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 1431 1432 /* 1433 * Second try 1434 */ 1435 x<< 0.25, 0.39, 0.415, 0.39; 1436 // do the computation 1437 lm.resetParameters(); 1438 info = lm.minimize(x); 1439 1440 // check return value 1441 VERIFY_IS_EQUAL(info, 1); 1442 LM_CHECK_N_ITERS(lm, 18, 16); 1443 // check norm^2 1444 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1445 // check x 1446 VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 1447 VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 1448 VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 1449 VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 1450 } 1451 1452 1453 1454 struct Bennett5_functor : Functor<double> 1455 { 1456 Bennett5_functor(void) : Functor<double>(3,154) {} 1457 static const double x[154]; 1458 static const double y[154]; 1459 int operator()(const VectorXd &b, VectorXd &fvec) 1460 { 1461 assert(b.size()==3); 1462 assert(fvec.size()==154); 1463 for(int i=0; i<154; i++) 1464 fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; 1465 return 0; 1466 } 1467 int df(const VectorXd &b, MatrixXd &fjac) 1468 { 1469 assert(b.size()==3); 1470 assert(fjac.rows()==154); 1471 assert(fjac.cols()==3); 1472 for(int i=0; i<154; i++) { 1473 double e = pow(b[1]+x[i],-1./b[2]); 1474 fjac(i,0) = e; 1475 fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); 1476 fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; 1477 } 1478 return 0; 1479 } 1480 }; 1481 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, 1482 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 }; 1483 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 1484 ,-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 , 1485 -31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; 1486 1487 // http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml 1488 void testNistBennett5(void) 1489 { 1490 const int n=3; 1491 int info; 1492 1493 VectorXd x(n); 1494 1495 /* 1496 * First try 1497 */ 1498 x<< -2000., 50., 0.8; 1499 // do the computation 1500 Bennett5_functor functor; 1501 LevenbergMarquardt<Bennett5_functor> lm(functor); 1502 lm.parameters.maxfev = 1000; 1503 info = lm.minimize(x); 1504 1505 // check return value 1506 VERIFY_IS_EQUAL(info, 1); 1507 LM_CHECK_N_ITERS(lm, 758, 744); 1508 // check norm^2 1509 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1510 // check x 1511 VERIFY_IS_APPROX(x[0], -2.5235058043E+03); 1512 VERIFY_IS_APPROX(x[1], 4.6736564644E+01); 1513 VERIFY_IS_APPROX(x[2], 9.3218483193E-01); 1514 /* 1515 * Second try 1516 */ 1517 x<< -1500., 45., 0.85; 1518 // do the computation 1519 lm.resetParameters(); 1520 info = lm.minimize(x); 1521 1522 // check return value 1523 VERIFY_IS_EQUAL(info, 1); 1524 LM_CHECK_N_ITERS(lm, 203, 192); 1525 // check norm^2 1526 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1527 // check x 1528 VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 1529 VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); 1530 VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); 1531 } 1532 1533 struct thurber_functor : Functor<double> 1534 { 1535 thurber_functor(void) : Functor<double>(7,37) {} 1536 static const double _x[37]; 1537 static const double _y[37]; 1538 int operator()(const VectorXd &b, VectorXd &fvec) 1539 { 1540 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 1541 assert(b.size()==7); 1542 assert(fvec.size()==37); 1543 for(int i=0; i<37; i++) { 1544 double x=_x[i], xx=x*x, xxx=xx*x; 1545 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]; 1546 } 1547 return 0; 1548 } 1549 int df(const VectorXd &b, MatrixXd &fjac) 1550 { 1551 assert(b.size()==7); 1552 assert(fjac.rows()==37); 1553 assert(fjac.cols()==7); 1554 for(int i=0; i<37; i++) { 1555 double x=_x[i], xx=x*x, xxx=xx*x; 1556 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 1557 fjac(i,0) = 1.*fact; 1558 fjac(i,1) = x*fact; 1559 fjac(i,2) = xx*fact; 1560 fjac(i,3) = xxx*fact; 1561 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 1562 fjac(i,4) = x*fact; 1563 fjac(i,5) = xx*fact; 1564 fjac(i,6) = xxx*fact; 1565 } 1566 return 0; 1567 } 1568 }; 1569 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 }; 1570 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}; 1571 1572 // http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml 1573 void testNistThurber(void) 1574 { 1575 const int n=7; 1576 int info; 1577 1578 VectorXd x(n); 1579 1580 /* 1581 * First try 1582 */ 1583 x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; 1584 // do the computation 1585 thurber_functor functor; 1586 LevenbergMarquardt<thurber_functor> lm(functor); 1587 lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1588 lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1589 info = lm.minimize(x); 1590 1591 // check return value 1592 VERIFY_IS_EQUAL(info, 1); 1593 LM_CHECK_N_ITERS(lm, 39,36); 1594 // check norm^2 1595 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1596 // check x 1597 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1598 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1599 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1600 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1601 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1602 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1603 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1604 1605 /* 1606 * Second try 1607 */ 1608 x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; 1609 // do the computation 1610 lm.resetParameters(); 1611 lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1612 lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1613 info = lm.minimize(x); 1614 1615 // check return value 1616 VERIFY_IS_EQUAL(info, 1); 1617 LM_CHECK_N_ITERS(lm, 29, 28); 1618 // check norm^2 1619 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1620 // check x 1621 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1622 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1623 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1624 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1625 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1626 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1627 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1628 } 1629 1630 struct rat43_functor : Functor<double> 1631 { 1632 rat43_functor(void) : Functor<double>(4,15) {} 1633 static const double x[15]; 1634 static const double y[15]; 1635 int operator()(const VectorXd &b, VectorXd &fvec) 1636 { 1637 assert(b.size()==4); 1638 assert(fvec.size()==15); 1639 for(int i=0; i<15; i++) 1640 fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; 1641 return 0; 1642 } 1643 int df(const VectorXd &b, MatrixXd &fjac) 1644 { 1645 assert(b.size()==4); 1646 assert(fjac.rows()==15); 1647 assert(fjac.cols()==4); 1648 for(int i=0; i<15; i++) { 1649 double e = exp(b[1]-b[2]*x[i]); 1650 double power = -1./b[3]; 1651 fjac(i,0) = pow(1.+e, power); 1652 fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); 1653 fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); 1654 fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); 1655 } 1656 return 0; 1657 } 1658 }; 1659 const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; 1660 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 }; 1661 1662 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml 1663 void testNistRat43(void) 1664 { 1665 const int n=4; 1666 int info; 1667 1668 VectorXd x(n); 1669 1670 /* 1671 * First try 1672 */ 1673 x<< 100., 10., 1., 1.; 1674 // do the computation 1675 rat43_functor functor; 1676 LevenbergMarquardt<rat43_functor> lm(functor); 1677 lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1678 lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1679 info = lm.minimize(x); 1680 1681 // check return value 1682 VERIFY_IS_EQUAL(info, 1); 1683 LM_CHECK_N_ITERS(lm, 27, 20); 1684 // check norm^2 1685 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1686 // check x 1687 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1688 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1689 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1690 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1691 1692 /* 1693 * Second try 1694 */ 1695 x<< 700., 5., 0.75, 1.3; 1696 // do the computation 1697 lm.resetParameters(); 1698 lm.parameters.ftol = 1.E5*NumTraits<double>::epsilon(); 1699 lm.parameters.xtol = 1.E5*NumTraits<double>::epsilon(); 1700 info = lm.minimize(x); 1701 1702 // check return value 1703 VERIFY_IS_EQUAL(info, 1); 1704 LM_CHECK_N_ITERS(lm, 9, 8); 1705 // check norm^2 1706 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1707 // check x 1708 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1709 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1710 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1711 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1712 } 1713 1714 1715 1716 struct eckerle4_functor : Functor<double> 1717 { 1718 eckerle4_functor(void) : Functor<double>(3,35) {} 1719 static const double x[35]; 1720 static const double y[35]; 1721 int operator()(const VectorXd &b, VectorXd &fvec) 1722 { 1723 assert(b.size()==3); 1724 assert(fvec.size()==35); 1725 for(int i=0; i<35; i++) 1726 fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; 1727 return 0; 1728 } 1729 int df(const VectorXd &b, MatrixXd &fjac) 1730 { 1731 assert(b.size()==3); 1732 assert(fjac.rows()==35); 1733 assert(fjac.cols()==3); 1734 for(int i=0; i<35; i++) { 1735 double b12 = b[1]*b[1]; 1736 double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); 1737 fjac(i,0) = e / b[1]; 1738 fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; 1739 fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; 1740 } 1741 return 0; 1742 } 1743 }; 1744 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}; 1745 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 }; 1746 1747 // http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml 1748 void testNistEckerle4(void) 1749 { 1750 const int n=3; 1751 int info; 1752 1753 VectorXd x(n); 1754 1755 /* 1756 * First try 1757 */ 1758 x<< 1., 10., 500.; 1759 // do the computation 1760 eckerle4_functor functor; 1761 LevenbergMarquardt<eckerle4_functor> lm(functor); 1762 info = lm.minimize(x); 1763 1764 // check return value 1765 VERIFY_IS_EQUAL(info, 1); 1766 LM_CHECK_N_ITERS(lm, 18, 15); 1767 // check norm^2 1768 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1769 // check x 1770 VERIFY_IS_APPROX(x[0], 1.5543827178); 1771 VERIFY_IS_APPROX(x[1], 4.0888321754); 1772 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1773 1774 /* 1775 * Second try 1776 */ 1777 x<< 1.5, 5., 450.; 1778 // do the computation 1779 info = lm.minimize(x); 1780 1781 // check return value 1782 VERIFY_IS_EQUAL(info, 1); 1783 LM_CHECK_N_ITERS(lm, 7, 6); 1784 // check norm^2 1785 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1786 // check x 1787 VERIFY_IS_APPROX(x[0], 1.5543827178); 1788 VERIFY_IS_APPROX(x[1], 4.0888321754); 1789 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1790 } 1791 1792 EIGEN_DECLARE_TEST(NonLinearOptimization) 1793 { 1794 // Tests using the examples provided by (c)minpack 1795 CALL_SUBTEST/*_1*/(testChkder()); 1796 CALL_SUBTEST/*_1*/(testLmder1()); 1797 CALL_SUBTEST/*_1*/(testLmder()); 1798 CALL_SUBTEST/*_2*/(testHybrj1()); 1799 CALL_SUBTEST/*_2*/(testHybrj()); 1800 CALL_SUBTEST/*_2*/(testHybrd1()); 1801 CALL_SUBTEST/*_2*/(testHybrd()); 1802 CALL_SUBTEST/*_3*/(testLmstr1()); 1803 CALL_SUBTEST/*_3*/(testLmstr()); 1804 CALL_SUBTEST/*_3*/(testLmdif1()); 1805 CALL_SUBTEST/*_3*/(testLmdif()); 1806 1807 // NIST tests, level of difficulty = "Lower" 1808 CALL_SUBTEST/*_4*/(testNistMisra1a()); 1809 CALL_SUBTEST/*_4*/(testNistChwirut2()); 1810 1811 // NIST tests, level of difficulty = "Average" 1812 CALL_SUBTEST/*_5*/(testNistHahn1()); 1813 CALL_SUBTEST/*_6*/(testNistMisra1d()); 1814 CALL_SUBTEST/*_7*/(testNistMGH17()); 1815 CALL_SUBTEST/*_8*/(testNistLanczos1()); 1816 1817 // // NIST tests, level of difficulty = "Higher" 1818 CALL_SUBTEST/*_9*/(testNistRat42()); 1819 // CALL_SUBTEST/*_10*/(testNistMGH10()); 1820 CALL_SUBTEST/*_11*/(testNistBoxBOD()); 1821 // CALL_SUBTEST/*_12*/(testNistMGH09()); 1822 CALL_SUBTEST/*_13*/(testNistBennett5()); 1823 CALL_SUBTEST/*_14*/(testNistThurber()); 1824 CALL_SUBTEST/*_15*/(testNistRat43()); 1825 CALL_SUBTEST/*_16*/(testNistEckerle4()); 1826 } 1827 1828 /* 1829 * Can be useful for debugging... 1830 printf("info, nfev : %d, %d\n", info, lm.nfev); 1831 printf("info, nfev, njev : %d, %d, %d\n", info, solver.nfev, solver.njev); 1832 printf("info, nfev : %d, %d\n", info, solver.nfev); 1833 printf("x[0] : %.32g\n", x[0]); 1834 printf("x[1] : %.32g\n", x[1]); 1835 printf("x[2] : %.32g\n", x[2]); 1836 printf("x[3] : %.32g\n", x[3]); 1837 printf("fvec.blueNorm() : %.32g\n", solver.fvec.blueNorm()); 1838 printf("fvec.blueNorm() : %.32g\n", lm.fvec.blueNorm()); 1839 1840 printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev); 1841 printf("fvec.squaredNorm() : %.13g\n", lm.fvec.squaredNorm()); 1842 std::cout << x << std::endl; 1843 std::cout.precision(9); 1844 std::cout << x[0] << std::endl; 1845 std::cout << x[1] << std::endl; 1846 std::cout << x[2] << std::endl; 1847 std::cout << x[3] << std::endl; 1848 */ 1849