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 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 int fcn_chkder(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag) 16 { 17 /* subroutine fcn for chkder example. */ 18 19 int i; 20 assert(15 == fvec.size()); 21 assert(3 == x.size()); 22 double tmp1, tmp2, tmp3, tmp4; 23 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, 24 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 25 26 27 if (iflag == 0) 28 return 0; 29 30 if (iflag != 2) 31 for (i=0; i<15; i++) { 32 tmp1 = i+1; 33 tmp2 = 16-i-1; 34 tmp3 = tmp1; 35 if (i >= 8) tmp3 = tmp2; 36 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 37 } 38 else { 39 for (i = 0; i < 15; i++) { 40 tmp1 = i+1; 41 tmp2 = 16-i-1; 42 43 /* error introduced into next statement for illustration. */ 44 /* corrected statement should read tmp3 = tmp1 . */ 45 46 tmp3 = tmp2; 47 if (i >= 8) tmp3 = tmp2; 48 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4=tmp4*tmp4; 49 fjac(i,0) = -1.; 50 fjac(i,1) = tmp1*tmp2/tmp4; 51 fjac(i,2) = tmp1*tmp3/tmp4; 52 } 53 } 54 return 0; 55 } 56 57 58 void testChkder() 59 { 60 const int m=15, n=3; 61 VectorXd x(n), fvec(m), xp, fvecp(m), err; 62 MatrixXd fjac(m,n); 63 VectorXi ipvt; 64 65 /* the following values should be suitable for */ 66 /* checking the jacobian matrix. */ 67 x << 9.2e-1, 1.3e-1, 5.4e-1; 68 69 internal::chkder(x, fvec, fjac, xp, fvecp, 1, err); 70 fcn_chkder(x, fvec, fjac, 1); 71 fcn_chkder(x, fvec, fjac, 2); 72 fcn_chkder(xp, fvecp, fjac, 1); 73 internal::chkder(x, fvec, fjac, xp, fvecp, 2, err); 74 75 fvecp -= fvec; 76 77 // check those 78 VectorXd fvec_ref(m), fvecp_ref(m), err_ref(m); 79 fvec_ref << 80 -1.181606, -1.429655, -1.606344, 81 -1.745269, -1.840654, -1.921586, 82 -1.984141, -2.022537, -2.468977, 83 -2.827562, -3.473582, -4.437612, 84 -6.047662, -9.267761, -18.91806; 85 fvecp_ref << 86 -7.724666e-09, -3.432406e-09, -2.034843e-10, 87 2.313685e-09, 4.331078e-09, 5.984096e-09, 88 7.363281e-09, 8.53147e-09, 1.488591e-08, 89 2.33585e-08, 3.522012e-08, 5.301255e-08, 90 8.26666e-08, 1.419747e-07, 3.19899e-07; 91 err_ref << 92 0.1141397, 0.09943516, 0.09674474, 93 0.09980447, 0.1073116, 0.1220445, 94 0.1526814, 1, 1, 95 1, 1, 1, 96 1, 1, 1; 97 98 VERIFY_IS_APPROX(fvec, fvec_ref); 99 VERIFY_IS_APPROX(fvecp, fvecp_ref); 100 VERIFY_IS_APPROX(err, err_ref); 101 } 102 103 // Generic functor 104 template<typename _Scalar, int NX=Dynamic, int NY=Dynamic> 105 struct Functor 106 { 107 typedef _Scalar Scalar; 108 enum { 109 InputsAtCompileTime = NX, 110 ValuesAtCompileTime = NY 111 }; 112 typedef Matrix<Scalar,InputsAtCompileTime,1> InputType; 113 typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType; 114 typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; 115 116 const int m_inputs, m_values; 117 118 Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} 119 Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {} 120 121 int inputs() const { return m_inputs; } 122 int values() const { return m_values; } 123 124 // you should define that in the subclass : 125 // void operator() (const InputType& x, ValueType* v, JacobianType* _j=0) const; 126 }; 127 128 struct lmder_functor : Functor<double> 129 { 130 lmder_functor(void): Functor<double>(3,15) {} 131 int operator()(const VectorXd &x, VectorXd &fvec) const 132 { 133 double tmp1, tmp2, tmp3; 134 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, 135 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 136 137 for (int i = 0; i < values(); i++) 138 { 139 tmp1 = i+1; 140 tmp2 = 16 - i - 1; 141 tmp3 = (i>=8)? tmp2 : tmp1; 142 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 143 } 144 return 0; 145 } 146 147 int df(const VectorXd &x, MatrixXd &fjac) const 148 { 149 double tmp1, tmp2, tmp3, tmp4; 150 for (int i = 0; i < values(); i++) 151 { 152 tmp1 = i+1; 153 tmp2 = 16 - i - 1; 154 tmp3 = (i>=8)? tmp2 : tmp1; 155 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 156 fjac(i,0) = -1; 157 fjac(i,1) = tmp1*tmp2/tmp4; 158 fjac(i,2) = tmp1*tmp3/tmp4; 159 } 160 return 0; 161 } 162 }; 163 164 void testLmder1() 165 { 166 int n=3, info; 167 168 VectorXd x; 169 170 /* the following starting values provide a rough fit. */ 171 x.setConstant(n, 1.); 172 173 // do the computation 174 lmder_functor functor; 175 LevenbergMarquardt<lmder_functor> lm(functor); 176 info = lm.lmder1(x); 177 178 // check return value 179 VERIFY_IS_EQUAL(info, 1); 180 VERIFY_IS_EQUAL(lm.nfev, 6); 181 VERIFY_IS_EQUAL(lm.njev, 5); 182 183 // check norm 184 VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 185 186 // check x 187 VectorXd x_ref(n); 188 x_ref << 0.08241058, 1.133037, 2.343695; 189 VERIFY_IS_APPROX(x, x_ref); 190 } 191 192 void testLmder() 193 { 194 const int m=15, n=3; 195 int info; 196 double fnorm, covfac; 197 VectorXd x; 198 199 /* the following starting values provide a rough fit. */ 200 x.setConstant(n, 1.); 201 202 // do the computation 203 lmder_functor functor; 204 LevenbergMarquardt<lmder_functor> lm(functor); 205 info = lm.minimize(x); 206 207 // check return values 208 VERIFY_IS_EQUAL(info, 1); 209 VERIFY_IS_EQUAL(lm.nfev, 6); 210 VERIFY_IS_EQUAL(lm.njev, 5); 211 212 // check norm 213 fnorm = lm.fvec.blueNorm(); 214 VERIFY_IS_APPROX(fnorm, 0.09063596); 215 216 // check x 217 VectorXd x_ref(n); 218 x_ref << 0.08241058, 1.133037, 2.343695; 219 VERIFY_IS_APPROX(x, x_ref); 220 221 // check covariance 222 covfac = fnorm*fnorm/(m-n); 223 internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 224 225 MatrixXd cov_ref(n,n); 226 cov_ref << 227 0.0001531202, 0.002869941, -0.002656662, 228 0.002869941, 0.09480935, -0.09098995, 229 -0.002656662, -0.09098995, 0.08778727; 230 231 // std::cout << fjac*covfac << std::endl; 232 233 MatrixXd cov; 234 cov = covfac*lm.fjac.topLeftCorner<n,n>(); 235 VERIFY_IS_APPROX( cov, cov_ref); 236 // TODO: why isn't this allowed ? : 237 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 238 } 239 240 struct hybrj_functor : Functor<double> 241 { 242 hybrj_functor(void) : Functor<double>(9,9) {} 243 244 int operator()(const VectorXd &x, VectorXd &fvec) 245 { 246 double temp, temp1, temp2; 247 const int n = x.size(); 248 assert(fvec.size()==n); 249 for (int k = 0; k < n; k++) 250 { 251 temp = (3. - 2.*x[k])*x[k]; 252 temp1 = 0.; 253 if (k) temp1 = x[k-1]; 254 temp2 = 0.; 255 if (k != n-1) temp2 = x[k+1]; 256 fvec[k] = temp - temp1 - 2.*temp2 + 1.; 257 } 258 return 0; 259 } 260 int df(const VectorXd &x, MatrixXd &fjac) 261 { 262 const int n = x.size(); 263 assert(fjac.rows()==n); 264 assert(fjac.cols()==n); 265 for (int k = 0; k < n; k++) 266 { 267 for (int j = 0; j < n; j++) 268 fjac(k,j) = 0.; 269 fjac(k,k) = 3.- 4.*x[k]; 270 if (k) fjac(k,k-1) = -1.; 271 if (k != n-1) fjac(k,k+1) = -2.; 272 } 273 return 0; 274 } 275 }; 276 277 278 void testHybrj1() 279 { 280 const int n=9; 281 int info; 282 VectorXd x(n); 283 284 /* the following starting values provide a rough fit. */ 285 x.setConstant(n, -1.); 286 287 // do the computation 288 hybrj_functor functor; 289 HybridNonLinearSolver<hybrj_functor> solver(functor); 290 info = solver.hybrj1(x); 291 292 // check return value 293 VERIFY_IS_EQUAL(info, 1); 294 VERIFY_IS_EQUAL(solver.nfev, 11); 295 VERIFY_IS_EQUAL(solver.njev, 1); 296 297 // check norm 298 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 299 300 301 // check x 302 VectorXd x_ref(n); 303 x_ref << 304 -0.5706545, -0.6816283, -0.7017325, 305 -0.7042129, -0.701369, -0.6918656, 306 -0.665792, -0.5960342, -0.4164121; 307 VERIFY_IS_APPROX(x, x_ref); 308 } 309 310 void testHybrj() 311 { 312 const int n=9; 313 int info; 314 VectorXd x(n); 315 316 /* the following starting values provide a rough fit. */ 317 x.setConstant(n, -1.); 318 319 320 // do the computation 321 hybrj_functor functor; 322 HybridNonLinearSolver<hybrj_functor> solver(functor); 323 solver.diag.setConstant(n, 1.); 324 solver.useExternalScaling = true; 325 info = solver.solve(x); 326 327 // check return value 328 VERIFY_IS_EQUAL(info, 1); 329 VERIFY_IS_EQUAL(solver.nfev, 11); 330 VERIFY_IS_EQUAL(solver.njev, 1); 331 332 // check norm 333 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 334 335 336 // check x 337 VectorXd x_ref(n); 338 x_ref << 339 -0.5706545, -0.6816283, -0.7017325, 340 -0.7042129, -0.701369, -0.6918656, 341 -0.665792, -0.5960342, -0.4164121; 342 VERIFY_IS_APPROX(x, x_ref); 343 344 } 345 346 struct hybrd_functor : Functor<double> 347 { 348 hybrd_functor(void) : Functor<double>(9,9) {} 349 int operator()(const VectorXd &x, VectorXd &fvec) const 350 { 351 double temp, temp1, temp2; 352 const int n = x.size(); 353 354 assert(fvec.size()==n); 355 for (int k=0; k < n; k++) 356 { 357 temp = (3. - 2.*x[k])*x[k]; 358 temp1 = 0.; 359 if (k) temp1 = x[k-1]; 360 temp2 = 0.; 361 if (k != n-1) temp2 = x[k+1]; 362 fvec[k] = temp - temp1 - 2.*temp2 + 1.; 363 } 364 return 0; 365 } 366 }; 367 368 void testHybrd1() 369 { 370 int n=9, info; 371 VectorXd x(n); 372 373 /* the following starting values provide a rough solution. */ 374 x.setConstant(n, -1.); 375 376 // do the computation 377 hybrd_functor functor; 378 HybridNonLinearSolver<hybrd_functor> solver(functor); 379 info = solver.hybrd1(x); 380 381 // check return value 382 VERIFY_IS_EQUAL(info, 1); 383 VERIFY_IS_EQUAL(solver.nfev, 20); 384 385 // check norm 386 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 387 388 // check x 389 VectorXd x_ref(n); 390 x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121; 391 VERIFY_IS_APPROX(x, x_ref); 392 } 393 394 void testHybrd() 395 { 396 const int n=9; 397 int info; 398 VectorXd x; 399 400 /* the following starting values provide a rough fit. */ 401 x.setConstant(n, -1.); 402 403 // do the computation 404 hybrd_functor functor; 405 HybridNonLinearSolver<hybrd_functor> solver(functor); 406 solver.parameters.nb_of_subdiagonals = 1; 407 solver.parameters.nb_of_superdiagonals = 1; 408 solver.diag.setConstant(n, 1.); 409 solver.useExternalScaling = true; 410 info = solver.solveNumericalDiff(x); 411 412 // check return value 413 VERIFY_IS_EQUAL(info, 1); 414 VERIFY_IS_EQUAL(solver.nfev, 14); 415 416 // check norm 417 VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 418 419 // check x 420 VectorXd x_ref(n); 421 x_ref << 422 -0.5706545, -0.6816283, -0.7017325, 423 -0.7042129, -0.701369, -0.6918656, 424 -0.665792, -0.5960342, -0.4164121; 425 VERIFY_IS_APPROX(x, x_ref); 426 } 427 428 struct lmstr_functor : Functor<double> 429 { 430 lmstr_functor(void) : Functor<double>(3,15) {} 431 int operator()(const VectorXd &x, VectorXd &fvec) 432 { 433 /* subroutine fcn for lmstr1 example. */ 434 double tmp1, tmp2, tmp3; 435 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, 436 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 437 438 assert(15==fvec.size()); 439 assert(3==x.size()); 440 441 for (int i=0; i<15; i++) 442 { 443 tmp1 = i+1; 444 tmp2 = 16 - i - 1; 445 tmp3 = (i>=8)? tmp2 : tmp1; 446 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 447 } 448 return 0; 449 } 450 int df(const VectorXd &x, VectorXd &jac_row, VectorXd::Index rownb) 451 { 452 assert(x.size()==3); 453 assert(jac_row.size()==x.size()); 454 double tmp1, tmp2, tmp3, tmp4; 455 456 int i = rownb-2; 457 tmp1 = i+1; 458 tmp2 = 16 - i - 1; 459 tmp3 = (i>=8)? tmp2 : tmp1; 460 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 461 jac_row[0] = -1; 462 jac_row[1] = tmp1*tmp2/tmp4; 463 jac_row[2] = tmp1*tmp3/tmp4; 464 return 0; 465 } 466 }; 467 468 void testLmstr1() 469 { 470 const int n=3; 471 int info; 472 473 VectorXd x(n); 474 475 /* the following starting values provide a rough fit. */ 476 x.setConstant(n, 1.); 477 478 // do the computation 479 lmstr_functor functor; 480 LevenbergMarquardt<lmstr_functor> lm(functor); 481 info = lm.lmstr1(x); 482 483 // check return value 484 VERIFY_IS_EQUAL(info, 1); 485 VERIFY_IS_EQUAL(lm.nfev, 6); 486 VERIFY_IS_EQUAL(lm.njev, 5); 487 488 // check norm 489 VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 490 491 // check x 492 VectorXd x_ref(n); 493 x_ref << 0.08241058, 1.133037, 2.343695 ; 494 VERIFY_IS_APPROX(x, x_ref); 495 } 496 497 void testLmstr() 498 { 499 const int n=3; 500 int info; 501 double fnorm; 502 VectorXd x(n); 503 504 /* the following starting values provide a rough fit. */ 505 x.setConstant(n, 1.); 506 507 // do the computation 508 lmstr_functor functor; 509 LevenbergMarquardt<lmstr_functor> lm(functor); 510 info = lm.minimizeOptimumStorage(x); 511 512 // check return values 513 VERIFY_IS_EQUAL(info, 1); 514 VERIFY_IS_EQUAL(lm.nfev, 6); 515 VERIFY_IS_EQUAL(lm.njev, 5); 516 517 // check norm 518 fnorm = lm.fvec.blueNorm(); 519 VERIFY_IS_APPROX(fnorm, 0.09063596); 520 521 // check x 522 VectorXd x_ref(n); 523 x_ref << 0.08241058, 1.133037, 2.343695; 524 VERIFY_IS_APPROX(x, x_ref); 525 526 } 527 528 struct lmdif_functor : Functor<double> 529 { 530 lmdif_functor(void) : Functor<double>(3,15) {} 531 int operator()(const VectorXd &x, VectorXd &fvec) const 532 { 533 int i; 534 double tmp1,tmp2,tmp3; 535 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, 536 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; 537 538 assert(x.size()==3); 539 assert(fvec.size()==15); 540 for (i=0; i<15; i++) 541 { 542 tmp1 = i+1; 543 tmp2 = 15 - i; 544 tmp3 = tmp1; 545 546 if (i >= 8) tmp3 = tmp2; 547 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 548 } 549 return 0; 550 } 551 }; 552 553 void testLmdif1() 554 { 555 const int n=3; 556 int info; 557 558 VectorXd x(n), fvec(15); 559 560 /* the following starting values provide a rough fit. */ 561 x.setConstant(n, 1.); 562 563 // do the computation 564 lmdif_functor functor; 565 DenseIndex nfev; 566 info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev); 567 568 // check return value 569 VERIFY_IS_EQUAL(info, 1); 570 VERIFY_IS_EQUAL(nfev, 26); 571 572 // check norm 573 functor(x, fvec); 574 VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); 575 576 // check x 577 VectorXd x_ref(n); 578 x_ref << 0.0824106, 1.1330366, 2.3436947; 579 VERIFY_IS_APPROX(x, x_ref); 580 581 } 582 583 void testLmdif() 584 { 585 const int m=15, n=3; 586 int info; 587 double fnorm, covfac; 588 VectorXd x(n); 589 590 /* the following starting values provide a rough fit. */ 591 x.setConstant(n, 1.); 592 593 // do the computation 594 lmdif_functor functor; 595 NumericalDiff<lmdif_functor> numDiff(functor); 596 LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff); 597 info = lm.minimize(x); 598 599 // check return values 600 VERIFY_IS_EQUAL(info, 1); 601 VERIFY_IS_EQUAL(lm.nfev, 26); 602 603 // check norm 604 fnorm = lm.fvec.blueNorm(); 605 VERIFY_IS_APPROX(fnorm, 0.09063596); 606 607 // check x 608 VectorXd x_ref(n); 609 x_ref << 0.08241058, 1.133037, 2.343695; 610 VERIFY_IS_APPROX(x, x_ref); 611 612 // check covariance 613 covfac = fnorm*fnorm/(m-n); 614 internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 615 616 MatrixXd cov_ref(n,n); 617 cov_ref << 618 0.0001531202, 0.002869942, -0.002656662, 619 0.002869942, 0.09480937, -0.09098997, 620 -0.002656662, -0.09098997, 0.08778729; 621 622 // std::cout << fjac*covfac << std::endl; 623 624 MatrixXd cov; 625 cov = covfac*lm.fjac.topLeftCorner<n,n>(); 626 VERIFY_IS_APPROX( cov, cov_ref); 627 // TODO: why isn't this allowed ? : 628 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 629 } 630 631 struct chwirut2_functor : Functor<double> 632 { 633 chwirut2_functor(void) : Functor<double>(3,54) {} 634 static const double m_x[54]; 635 static const double m_y[54]; 636 int operator()(const VectorXd &b, VectorXd &fvec) 637 { 638 int i; 639 640 assert(b.size()==3); 641 assert(fvec.size()==54); 642 for(i=0; i<54; i++) { 643 double x = m_x[i]; 644 fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; 645 } 646 return 0; 647 } 648 int df(const VectorXd &b, MatrixXd &fjac) 649 { 650 assert(b.size()==3); 651 assert(fjac.rows()==54); 652 assert(fjac.cols()==3); 653 for(int i=0; i<54; i++) { 654 double x = m_x[i]; 655 double factor = 1./(b[1]+b[2]*x); 656 double e = exp(-b[0]*x); 657 fjac(i,0) = -x*e*factor; 658 fjac(i,1) = -e*factor*factor; 659 fjac(i,2) = -x*e*factor*factor; 660 } 661 return 0; 662 } 663 }; 664 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}; 665 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 }; 666 667 // http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml 668 void testNistChwirut2(void) 669 { 670 const int n=3; 671 int info; 672 673 VectorXd x(n); 674 675 /* 676 * First try 677 */ 678 x<< 0.1, 0.01, 0.02; 679 // do the computation 680 chwirut2_functor functor; 681 LevenbergMarquardt<chwirut2_functor> lm(functor); 682 info = lm.minimize(x); 683 684 // check return value 685 VERIFY_IS_EQUAL(info, 1); 686 VERIFY_IS_EQUAL(lm.nfev, 10); 687 VERIFY_IS_EQUAL(lm.njev, 8); 688 // check norm^2 689 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 690 // check x 691 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 692 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 693 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 694 695 /* 696 * Second try 697 */ 698 x<< 0.15, 0.008, 0.010; 699 // do the computation 700 lm.resetParameters(); 701 lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 702 lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 703 info = lm.minimize(x); 704 705 // check return value 706 VERIFY_IS_EQUAL(info, 1); 707 VERIFY_IS_EQUAL(lm.nfev, 7); 708 VERIFY_IS_EQUAL(lm.njev, 6); 709 // check norm^2 710 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 711 // check x 712 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 713 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 714 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 715 } 716 717 718 struct misra1a_functor : Functor<double> 719 { 720 misra1a_functor(void) : Functor<double>(2,14) {} 721 static const double m_x[14]; 722 static const double m_y[14]; 723 int operator()(const VectorXd &b, VectorXd &fvec) 724 { 725 assert(b.size()==2); 726 assert(fvec.size()==14); 727 for(int i=0; i<14; i++) { 728 fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; 729 } 730 return 0; 731 } 732 int df(const VectorXd &b, MatrixXd &fjac) 733 { 734 assert(b.size()==2); 735 assert(fjac.rows()==14); 736 assert(fjac.cols()==2); 737 for(int i=0; i<14; i++) { 738 fjac(i,0) = (1.-exp(-b[1]*m_x[i])); 739 fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); 740 } 741 return 0; 742 } 743 }; 744 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}; 745 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}; 746 747 // http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml 748 void testNistMisra1a(void) 749 { 750 const int n=2; 751 int info; 752 753 VectorXd x(n); 754 755 /* 756 * First try 757 */ 758 x<< 500., 0.0001; 759 // do the computation 760 misra1a_functor functor; 761 LevenbergMarquardt<misra1a_functor> lm(functor); 762 info = lm.minimize(x); 763 764 // check return value 765 VERIFY_IS_EQUAL(info, 1); 766 VERIFY_IS_EQUAL(lm.nfev, 19); 767 VERIFY_IS_EQUAL(lm.njev, 15); 768 // check norm^2 769 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 770 // check x 771 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 772 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 773 774 /* 775 * Second try 776 */ 777 x<< 250., 0.0005; 778 // do the computation 779 info = lm.minimize(x); 780 781 // check return value 782 VERIFY_IS_EQUAL(info, 1); 783 VERIFY_IS_EQUAL(lm.nfev, 5); 784 VERIFY_IS_EQUAL(lm.njev, 4); 785 // check norm^2 786 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 787 // check x 788 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 789 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 790 } 791 792 struct hahn1_functor : Functor<double> 793 { 794 hahn1_functor(void) : Functor<double>(7,236) {} 795 static const double m_x[236]; 796 int operator()(const VectorXd &b, VectorXd &fvec) 797 { 798 static const double m_y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0 , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0 , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0 , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 , 16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0 , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0 , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 12.596E0 , 13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0 , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0 , 20.935E0 , 21.035E0 , 20.93E0 , 21.074E0 , 21.085E0 , 20.935E0 }; 799 800 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 801 802 assert(b.size()==7); 803 assert(fvec.size()==236); 804 for(int i=0; i<236; i++) { 805 double x=m_x[i], xx=x*x, xxx=xx*x; 806 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]; 807 } 808 return 0; 809 } 810 811 int df(const VectorXd &b, MatrixXd &fjac) 812 { 813 assert(b.size()==7); 814 assert(fjac.rows()==236); 815 assert(fjac.cols()==7); 816 for(int i=0; i<236; i++) { 817 double x=m_x[i], xx=x*x, xxx=xx*x; 818 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 819 fjac(i,0) = 1.*fact; 820 fjac(i,1) = x*fact; 821 fjac(i,2) = xx*fact; 822 fjac(i,3) = xxx*fact; 823 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 824 fjac(i,4) = x*fact; 825 fjac(i,5) = xx*fact; 826 fjac(i,6) = xxx*fact; 827 } 828 return 0; 829 } 830 }; 831 const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , 282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 , 141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0 , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0 , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0}; 832 833 // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml 834 void testNistHahn1(void) 835 { 836 const int n=7; 837 int info; 838 839 VectorXd x(n); 840 841 /* 842 * First try 843 */ 844 x<< 10., -1., .05, -.00001, -.05, .001, -.000001; 845 // do the computation 846 hahn1_functor functor; 847 LevenbergMarquardt<hahn1_functor> lm(functor); 848 info = lm.minimize(x); 849 850 // check return value 851 VERIFY_IS_EQUAL(info, 1); 852 VERIFY_IS_EQUAL(lm.nfev, 11); 853 VERIFY_IS_EQUAL(lm.njev, 10); 854 // check norm^2 855 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 856 // check x 857 VERIFY_IS_APPROX(x[0], 1.0776351733E+00); 858 VERIFY_IS_APPROX(x[1],-1.2269296921E-01); 859 VERIFY_IS_APPROX(x[2], 4.0863750610E-03); 860 VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 861 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 862 VERIFY_IS_APPROX(x[5], 2.4053735503E-04); 863 VERIFY_IS_APPROX(x[6],-1.2314450199E-07); 864 865 /* 866 * Second try 867 */ 868 x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; 869 // do the computation 870 info = lm.minimize(x); 871 872 // check return value 873 VERIFY_IS_EQUAL(info, 1); 874 VERIFY_IS_EQUAL(lm.nfev, 11); 875 VERIFY_IS_EQUAL(lm.njev, 10); 876 // check norm^2 877 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 878 // check x 879 VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 880 VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 881 VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 882 VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 883 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 884 VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 885 VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 886 887 } 888 889 struct misra1d_functor : Functor<double> 890 { 891 misra1d_functor(void) : Functor<double>(2,14) {} 892 static const double x[14]; 893 static const double y[14]; 894 int operator()(const VectorXd &b, VectorXd &fvec) 895 { 896 assert(b.size()==2); 897 assert(fvec.size()==14); 898 for(int i=0; i<14; i++) { 899 fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; 900 } 901 return 0; 902 } 903 int df(const VectorXd &b, MatrixXd &fjac) 904 { 905 assert(b.size()==2); 906 assert(fjac.rows()==14); 907 assert(fjac.cols()==2); 908 for(int i=0; i<14; i++) { 909 double den = 1.+b[1]*x[i]; 910 fjac(i,0) = b[1]*x[i] / den; 911 fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; 912 } 913 return 0; 914 } 915 }; 916 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}; 917 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}; 918 919 // http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml 920 void testNistMisra1d(void) 921 { 922 const int n=2; 923 int info; 924 925 VectorXd x(n); 926 927 /* 928 * First try 929 */ 930 x<< 500., 0.0001; 931 // do the computation 932 misra1d_functor functor; 933 LevenbergMarquardt<misra1d_functor> lm(functor); 934 info = lm.minimize(x); 935 936 // check return value 937 VERIFY_IS_EQUAL(info, 3); 938 VERIFY_IS_EQUAL(lm.nfev, 9); 939 VERIFY_IS_EQUAL(lm.njev, 7); 940 // check norm^2 941 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 942 // check x 943 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 944 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 945 946 /* 947 * Second try 948 */ 949 x<< 450., 0.0003; 950 // do the computation 951 info = lm.minimize(x); 952 953 // check return value 954 VERIFY_IS_EQUAL(info, 1); 955 VERIFY_IS_EQUAL(lm.nfev, 4); 956 VERIFY_IS_EQUAL(lm.njev, 3); 957 // check norm^2 958 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 959 // check x 960 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 961 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 962 } 963 964 965 struct lanczos1_functor : Functor<double> 966 { 967 lanczos1_functor(void) : Functor<double>(6,24) {} 968 static const double x[24]; 969 static const double y[24]; 970 int operator()(const VectorXd &b, VectorXd &fvec) 971 { 972 assert(b.size()==6); 973 assert(fvec.size()==24); 974 for(int i=0; i<24; i++) 975 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]; 976 return 0; 977 } 978 int df(const VectorXd &b, MatrixXd &fjac) 979 { 980 assert(b.size()==6); 981 assert(fjac.rows()==24); 982 assert(fjac.cols()==6); 983 for(int i=0; i<24; i++) { 984 fjac(i,0) = exp(-b[1]*x[i]); 985 fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); 986 fjac(i,2) = exp(-b[3]*x[i]); 987 fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); 988 fjac(i,4) = exp(-b[5]*x[i]); 989 fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); 990 } 991 return 0; 992 } 993 }; 994 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 }; 995 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 }; 996 997 // http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml 998 void testNistLanczos1(void) 999 { 1000 const int n=6; 1001 int info; 1002 1003 VectorXd x(n); 1004 1005 /* 1006 * First try 1007 */ 1008 x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; 1009 // do the computation 1010 lanczos1_functor functor; 1011 LevenbergMarquardt<lanczos1_functor> lm(functor); 1012 info = lm.minimize(x); 1013 1014 // check return value 1015 VERIFY_IS_EQUAL(info, 2); 1016 VERIFY_IS_EQUAL(lm.nfev, 79); 1017 VERIFY_IS_EQUAL(lm.njev, 72); 1018 // check norm^2 1019 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430899764097e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 1020 // check x 1021 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1022 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1023 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1024 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1025 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1026 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1027 1028 /* 1029 * Second try 1030 */ 1031 x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; 1032 // do the computation 1033 info = lm.minimize(x); 1034 1035 // check return value 1036 VERIFY_IS_EQUAL(info, 2); 1037 VERIFY_IS_EQUAL(lm.nfev, 9); 1038 VERIFY_IS_EQUAL(lm.njev, 8); 1039 // check norm^2 1040 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.428595533845e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 1041 // check x 1042 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1043 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1044 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1045 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1046 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1047 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1048 1049 } 1050 1051 struct rat42_functor : Functor<double> 1052 { 1053 rat42_functor(void) : Functor<double>(3,9) {} 1054 static const double x[9]; 1055 static const double y[9]; 1056 int operator()(const VectorXd &b, VectorXd &fvec) 1057 { 1058 assert(b.size()==3); 1059 assert(fvec.size()==9); 1060 for(int i=0; i<9; i++) { 1061 fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; 1062 } 1063 return 0; 1064 } 1065 1066 int df(const VectorXd &b, MatrixXd &fjac) 1067 { 1068 assert(b.size()==3); 1069 assert(fjac.rows()==9); 1070 assert(fjac.cols()==3); 1071 for(int i=0; i<9; i++) { 1072 double e = exp(b[1]-b[2]*x[i]); 1073 fjac(i,0) = 1./(1.+e); 1074 fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); 1075 fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); 1076 } 1077 return 0; 1078 } 1079 }; 1080 const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; 1081 const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; 1082 1083 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml 1084 void testNistRat42(void) 1085 { 1086 const int n=3; 1087 int info; 1088 1089 VectorXd x(n); 1090 1091 /* 1092 * First try 1093 */ 1094 x<< 100., 1., 0.1; 1095 // do the computation 1096 rat42_functor functor; 1097 LevenbergMarquardt<rat42_functor> lm(functor); 1098 info = lm.minimize(x); 1099 1100 // check return value 1101 VERIFY_IS_EQUAL(info, 1); 1102 VERIFY_IS_EQUAL(lm.nfev, 10); 1103 VERIFY_IS_EQUAL(lm.njev, 8); 1104 // check norm^2 1105 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00); 1106 // check x 1107 VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 1108 VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 1109 VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 1110 1111 /* 1112 * Second try 1113 */ 1114 x<< 75., 2.5, 0.07; 1115 // do the computation 1116 info = lm.minimize(x); 1117 1118 // check return value 1119 VERIFY_IS_EQUAL(info, 1); 1120 VERIFY_IS_EQUAL(lm.nfev, 6); 1121 VERIFY_IS_EQUAL(lm.njev, 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 VERIFY_IS_EQUAL(lm.nfev, 284 ); 1181 VERIFY_IS_EQUAL(lm.njev, 249 ); 1182 // check norm^2 1183 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1184 // check x 1185 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1186 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1187 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1188 1189 /* 1190 * Second try 1191 */ 1192 x<< 0.02, 4000., 250.; 1193 // do the computation 1194 info = lm.minimize(x); 1195 1196 // check return value 1197 VERIFY_IS_EQUAL(info, 3); 1198 VERIFY_IS_EQUAL(lm.nfev, 126); 1199 VERIFY_IS_EQUAL(lm.njev, 116); 1200 // check norm^2 1201 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1202 // check x 1203 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1204 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1205 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1206 } 1207 1208 1209 struct BoxBOD_functor : Functor<double> 1210 { 1211 BoxBOD_functor(void) : Functor<double>(2,6) {} 1212 static const double x[6]; 1213 int operator()(const VectorXd &b, VectorXd &fvec) 1214 { 1215 static const double y[6] = { 109., 149., 149., 191., 213., 224. }; 1216 assert(b.size()==2); 1217 assert(fvec.size()==6); 1218 for(int i=0; i<6; i++) 1219 fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; 1220 return 0; 1221 } 1222 int df(const VectorXd &b, MatrixXd &fjac) 1223 { 1224 assert(b.size()==2); 1225 assert(fjac.rows()==6); 1226 assert(fjac.cols()==2); 1227 for(int i=0; i<6; i++) { 1228 double e = exp(-b[1]*x[i]); 1229 fjac(i,0) = 1.-e; 1230 fjac(i,1) = b[0]*x[i]*e; 1231 } 1232 return 0; 1233 } 1234 }; 1235 const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; 1236 1237 // http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml 1238 void testNistBoxBOD(void) 1239 { 1240 const int n=2; 1241 int info; 1242 1243 VectorXd x(n); 1244 1245 /* 1246 * First try 1247 */ 1248 x<< 1., 1.; 1249 // do the computation 1250 BoxBOD_functor functor; 1251 LevenbergMarquardt<BoxBOD_functor> lm(functor); 1252 lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1253 lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1254 lm.parameters.factor = 10.; 1255 info = lm.minimize(x); 1256 1257 // check return value 1258 VERIFY_IS_EQUAL(info, 1); 1259 VERIFY_IS_EQUAL(lm.nfev, 31); 1260 VERIFY_IS_EQUAL(lm.njev, 25); 1261 // check norm^2 1262 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1263 // check x 1264 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1265 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1266 1267 /* 1268 * Second try 1269 */ 1270 x<< 100., 0.75; 1271 // do the computation 1272 lm.resetParameters(); 1273 lm.parameters.ftol = NumTraits<double>::epsilon(); 1274 lm.parameters.xtol = NumTraits<double>::epsilon(); 1275 info = lm.minimize(x); 1276 1277 // check return value 1278 VERIFY_IS_EQUAL(info, 1); 1279 VERIFY_IS_EQUAL(lm.nfev, 15 ); 1280 VERIFY_IS_EQUAL(lm.njev, 14 ); 1281 // check norm^2 1282 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1283 // check x 1284 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1285 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1286 } 1287 1288 struct MGH17_functor : Functor<double> 1289 { 1290 MGH17_functor(void) : Functor<double>(5,33) {} 1291 static const double x[33]; 1292 static const double y[33]; 1293 int operator()(const VectorXd &b, VectorXd &fvec) 1294 { 1295 assert(b.size()==5); 1296 assert(fvec.size()==33); 1297 for(int i=0; i<33; i++) 1298 fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; 1299 return 0; 1300 } 1301 int df(const VectorXd &b, MatrixXd &fjac) 1302 { 1303 assert(b.size()==5); 1304 assert(fjac.rows()==33); 1305 assert(fjac.cols()==5); 1306 for(int i=0; i<33; i++) { 1307 fjac(i,0) = 1.; 1308 fjac(i,1) = exp(-b[3]*x[i]); 1309 fjac(i,2) = exp(-b[4]*x[i]); 1310 fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); 1311 fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); 1312 } 1313 return 0; 1314 } 1315 }; 1316 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 }; 1317 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 }; 1318 1319 // http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml 1320 void testNistMGH17(void) 1321 { 1322 const int n=5; 1323 int info; 1324 1325 VectorXd x(n); 1326 1327 /* 1328 * First try 1329 */ 1330 x<< 50., 150., -100., 1., 2.; 1331 // do the computation 1332 MGH17_functor functor; 1333 LevenbergMarquardt<MGH17_functor> lm(functor); 1334 lm.parameters.ftol = NumTraits<double>::epsilon(); 1335 lm.parameters.xtol = NumTraits<double>::epsilon(); 1336 lm.parameters.maxfev = 1000; 1337 info = lm.minimize(x); 1338 1339 // check return value 1340 VERIFY_IS_EQUAL(info, 2); 1341 VERIFY_IS_EQUAL(lm.nfev, 602 ); 1342 VERIFY_IS_EQUAL(lm.njev, 545 ); 1343 // check norm^2 1344 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1345 // check x 1346 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1347 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1348 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1349 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1350 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1351 1352 /* 1353 * Second try 1354 */ 1355 x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; 1356 // do the computation 1357 lm.resetParameters(); 1358 info = lm.minimize(x); 1359 1360 // check return value 1361 VERIFY_IS_EQUAL(info, 1); 1362 VERIFY_IS_EQUAL(lm.nfev, 18); 1363 VERIFY_IS_EQUAL(lm.njev, 15); 1364 // check norm^2 1365 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1366 // check x 1367 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1368 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1369 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1370 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1371 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1372 } 1373 1374 struct MGH09_functor : Functor<double> 1375 { 1376 MGH09_functor(void) : Functor<double>(4,11) {} 1377 static const double _x[11]; 1378 static const double y[11]; 1379 int operator()(const VectorXd &b, VectorXd &fvec) 1380 { 1381 assert(b.size()==4); 1382 assert(fvec.size()==11); 1383 for(int i=0; i<11; i++) { 1384 double x = _x[i], xx=x*x; 1385 fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; 1386 } 1387 return 0; 1388 } 1389 int df(const VectorXd &b, MatrixXd &fjac) 1390 { 1391 assert(b.size()==4); 1392 assert(fjac.rows()==11); 1393 assert(fjac.cols()==4); 1394 for(int i=0; i<11; i++) { 1395 double x = _x[i], xx=x*x; 1396 double factor = 1./(xx+x*b[2]+b[3]); 1397 fjac(i,0) = (xx+x*b[1]) * factor; 1398 fjac(i,1) = b[0]*x* factor; 1399 fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; 1400 fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; 1401 } 1402 return 0; 1403 } 1404 }; 1405 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 }; 1406 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 }; 1407 1408 // http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml 1409 void testNistMGH09(void) 1410 { 1411 const int n=4; 1412 int info; 1413 1414 VectorXd x(n); 1415 1416 /* 1417 * First try 1418 */ 1419 x<< 25., 39, 41.5, 39.; 1420 // do the computation 1421 MGH09_functor functor; 1422 LevenbergMarquardt<MGH09_functor> lm(functor); 1423 lm.parameters.maxfev = 1000; 1424 info = lm.minimize(x); 1425 1426 // check return value 1427 VERIFY_IS_EQUAL(info, 1); 1428 VERIFY_IS_EQUAL(lm.nfev, 490 ); 1429 VERIFY_IS_EQUAL(lm.njev, 376 ); 1430 // check norm^2 1431 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1432 // check x 1433 VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 1434 VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 1435 VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 1436 VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 1437 1438 /* 1439 * Second try 1440 */ 1441 x<< 0.25, 0.39, 0.415, 0.39; 1442 // do the computation 1443 lm.resetParameters(); 1444 info = lm.minimize(x); 1445 1446 // check return value 1447 VERIFY_IS_EQUAL(info, 1); 1448 VERIFY_IS_EQUAL(lm.nfev, 18); 1449 VERIFY_IS_EQUAL(lm.njev, 16); 1450 // check norm^2 1451 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1452 // check x 1453 VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 1454 VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 1455 VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 1456 VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 1457 } 1458 1459 1460 1461 struct Bennett5_functor : Functor<double> 1462 { 1463 Bennett5_functor(void) : Functor<double>(3,154) {} 1464 static const double x[154]; 1465 static const double y[154]; 1466 int operator()(const VectorXd &b, VectorXd &fvec) 1467 { 1468 assert(b.size()==3); 1469 assert(fvec.size()==154); 1470 for(int i=0; i<154; i++) 1471 fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; 1472 return 0; 1473 } 1474 int df(const VectorXd &b, MatrixXd &fjac) 1475 { 1476 assert(b.size()==3); 1477 assert(fjac.rows()==154); 1478 assert(fjac.cols()==3); 1479 for(int i=0; i<154; i++) { 1480 double e = pow(b[1]+x[i],-1./b[2]); 1481 fjac(i,0) = e; 1482 fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); 1483 fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; 1484 } 1485 return 0; 1486 } 1487 }; 1488 const double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, 11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 }; 1489 const double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 ,-32.238602E0 ,-32.229900E0 ,-32.224098E0 ,-32.215401E0 ,-32.203800E0 ,-32.198002E0 ,-32.189400E0 ,-32.183601E0 ,-32.174900E0 ,-32.169102E0 ,-32.163300E0 ,-32.154598E0 ,-32.145901E0 ,-32.140099E0 ,-32.131401E0 ,-32.125599E0 ,-32.119801E0 ,-32.111198E0 ,-32.105400E0 ,-32.096699E0 ,-32.090900E0 ,-32.088001E0 ,-32.079300E0 ,-32.073502E0 ,-32.067699E0 ,-32.061901E0 ,-32.056099E0 ,-32.050301E0 ,-32.044498E0 ,-32.038799E0 ,-32.033001E0 ,-32.027199E0 ,-32.024300E0 ,-32.018501E0 ,-32.012699E0 ,-32.004002E0 ,-32.001099E0 ,-31.995300E0 ,-31.989500E0 ,-31.983700E0 ,-31.977900E0 ,-31.972099E0 ,-31.969299E0 ,-31.963501E0 ,-31.957701E0 ,-31.951900E0 ,-31.946100E0 ,-31.940300E0 ,-31.937401E0 ,-31.931601E0 ,-31.925800E0 ,-31.922899E0 ,-31.917101E0 ,-31.911301E0 ,-31.908400E0 ,-31.902599E0 ,-31.896900E0 ,-31.893999E0 ,-31.888201E0 ,-31.885300E0 ,-31.882401E0 ,-31.876600E0 ,-31.873699E0 ,-31.867901E0 ,-31.862101E0 ,-31.859200E0 ,-31.856300E0 ,-31.850500E0 ,-31.844700E0 ,-31.841801E0 ,-31.838900E0 ,-31.833099E0 ,-31.830200E0 ,-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; 1490 1491 // http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml 1492 void testNistBennett5(void) 1493 { 1494 const int n=3; 1495 int info; 1496 1497 VectorXd x(n); 1498 1499 /* 1500 * First try 1501 */ 1502 x<< -2000., 50., 0.8; 1503 // do the computation 1504 Bennett5_functor functor; 1505 LevenbergMarquardt<Bennett5_functor> lm(functor); 1506 lm.parameters.maxfev = 1000; 1507 info = lm.minimize(x); 1508 1509 // check return value 1510 VERIFY_IS_EQUAL(info, 1); 1511 VERIFY_IS_EQUAL(lm.nfev, 758); 1512 VERIFY_IS_EQUAL(lm.njev, 744); 1513 // check norm^2 1514 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1515 // check x 1516 VERIFY_IS_APPROX(x[0], -2.5235058043E+03); 1517 VERIFY_IS_APPROX(x[1], 4.6736564644E+01); 1518 VERIFY_IS_APPROX(x[2], 9.3218483193E-01); 1519 /* 1520 * Second try 1521 */ 1522 x<< -1500., 45., 0.85; 1523 // do the computation 1524 lm.resetParameters(); 1525 info = lm.minimize(x); 1526 1527 // check return value 1528 VERIFY_IS_EQUAL(info, 1); 1529 VERIFY_IS_EQUAL(lm.nfev, 203); 1530 VERIFY_IS_EQUAL(lm.njev, 192); 1531 // check norm^2 1532 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1533 // check x 1534 VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 1535 VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); 1536 VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); 1537 } 1538 1539 struct thurber_functor : Functor<double> 1540 { 1541 thurber_functor(void) : Functor<double>(7,37) {} 1542 static const double _x[37]; 1543 static const double _y[37]; 1544 int operator()(const VectorXd &b, VectorXd &fvec) 1545 { 1546 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 1547 assert(b.size()==7); 1548 assert(fvec.size()==37); 1549 for(int i=0; i<37; i++) { 1550 double x=_x[i], xx=x*x, xxx=xx*x; 1551 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]; 1552 } 1553 return 0; 1554 } 1555 int df(const VectorXd &b, MatrixXd &fjac) 1556 { 1557 assert(b.size()==7); 1558 assert(fjac.rows()==37); 1559 assert(fjac.cols()==7); 1560 for(int i=0; i<37; i++) { 1561 double x=_x[i], xx=x*x, xxx=xx*x; 1562 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 1563 fjac(i,0) = 1.*fact; 1564 fjac(i,1) = x*fact; 1565 fjac(i,2) = xx*fact; 1566 fjac(i,3) = xxx*fact; 1567 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 1568 fjac(i,4) = x*fact; 1569 fjac(i,5) = xx*fact; 1570 fjac(i,6) = xxx*fact; 1571 } 1572 return 0; 1573 } 1574 }; 1575 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 }; 1576 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}; 1577 1578 // http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml 1579 void testNistThurber(void) 1580 { 1581 const int n=7; 1582 int info; 1583 1584 VectorXd x(n); 1585 1586 /* 1587 * First try 1588 */ 1589 x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; 1590 // do the computation 1591 thurber_functor functor; 1592 LevenbergMarquardt<thurber_functor> lm(functor); 1593 lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1594 lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1595 info = lm.minimize(x); 1596 1597 // check return value 1598 VERIFY_IS_EQUAL(info, 1); 1599 VERIFY_IS_EQUAL(lm.nfev, 39); 1600 VERIFY_IS_EQUAL(lm.njev, 36); 1601 // check norm^2 1602 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1603 // check x 1604 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1605 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1606 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1607 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1608 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1609 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1610 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1611 1612 /* 1613 * Second try 1614 */ 1615 x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; 1616 // do the computation 1617 lm.resetParameters(); 1618 lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1619 lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1620 info = lm.minimize(x); 1621 1622 // check return value 1623 VERIFY_IS_EQUAL(info, 1); 1624 VERIFY_IS_EQUAL(lm.nfev, 29); 1625 VERIFY_IS_EQUAL(lm.njev, 28); 1626 // check norm^2 1627 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1628 // check x 1629 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1630 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1631 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1632 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1633 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1634 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1635 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1636 } 1637 1638 struct rat43_functor : Functor<double> 1639 { 1640 rat43_functor(void) : Functor<double>(4,15) {} 1641 static const double x[15]; 1642 static const double y[15]; 1643 int operator()(const VectorXd &b, VectorXd &fvec) 1644 { 1645 assert(b.size()==4); 1646 assert(fvec.size()==15); 1647 for(int i=0; i<15; i++) 1648 fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; 1649 return 0; 1650 } 1651 int df(const VectorXd &b, MatrixXd &fjac) 1652 { 1653 assert(b.size()==4); 1654 assert(fjac.rows()==15); 1655 assert(fjac.cols()==4); 1656 for(int i=0; i<15; i++) { 1657 double e = exp(b[1]-b[2]*x[i]); 1658 double power = -1./b[3]; 1659 fjac(i,0) = pow(1.+e, power); 1660 fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); 1661 fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); 1662 fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); 1663 } 1664 return 0; 1665 } 1666 }; 1667 const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; 1668 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 }; 1669 1670 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml 1671 void testNistRat43(void) 1672 { 1673 const int n=4; 1674 int info; 1675 1676 VectorXd x(n); 1677 1678 /* 1679 * First try 1680 */ 1681 x<< 100., 10., 1., 1.; 1682 // do the computation 1683 rat43_functor functor; 1684 LevenbergMarquardt<rat43_functor> lm(functor); 1685 lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1686 lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1687 info = lm.minimize(x); 1688 1689 // check return value 1690 VERIFY_IS_EQUAL(info, 1); 1691 VERIFY_IS_EQUAL(lm.nfev, 27); 1692 VERIFY_IS_EQUAL(lm.njev, 20); 1693 // check norm^2 1694 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1695 // check x 1696 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1697 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1698 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1699 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1700 1701 /* 1702 * Second try 1703 */ 1704 x<< 700., 5., 0.75, 1.3; 1705 // do the computation 1706 lm.resetParameters(); 1707 lm.parameters.ftol = 1.E5*NumTraits<double>::epsilon(); 1708 lm.parameters.xtol = 1.E5*NumTraits<double>::epsilon(); 1709 info = lm.minimize(x); 1710 1711 // check return value 1712 VERIFY_IS_EQUAL(info, 1); 1713 VERIFY_IS_EQUAL(lm.nfev, 9); 1714 VERIFY_IS_EQUAL(lm.njev, 8); 1715 // check norm^2 1716 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1717 // check x 1718 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1719 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1720 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1721 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1722 } 1723 1724 1725 1726 struct eckerle4_functor : Functor<double> 1727 { 1728 eckerle4_functor(void) : Functor<double>(3,35) {} 1729 static const double x[35]; 1730 static const double y[35]; 1731 int operator()(const VectorXd &b, VectorXd &fvec) 1732 { 1733 assert(b.size()==3); 1734 assert(fvec.size()==35); 1735 for(int i=0; i<35; i++) 1736 fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; 1737 return 0; 1738 } 1739 int df(const VectorXd &b, MatrixXd &fjac) 1740 { 1741 assert(b.size()==3); 1742 assert(fjac.rows()==35); 1743 assert(fjac.cols()==3); 1744 for(int i=0; i<35; i++) { 1745 double b12 = b[1]*b[1]; 1746 double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); 1747 fjac(i,0) = e / b[1]; 1748 fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; 1749 fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; 1750 } 1751 return 0; 1752 } 1753 }; 1754 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}; 1755 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 }; 1756 1757 // http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml 1758 void testNistEckerle4(void) 1759 { 1760 const int n=3; 1761 int info; 1762 1763 VectorXd x(n); 1764 1765 /* 1766 * First try 1767 */ 1768 x<< 1., 10., 500.; 1769 // do the computation 1770 eckerle4_functor functor; 1771 LevenbergMarquardt<eckerle4_functor> lm(functor); 1772 info = lm.minimize(x); 1773 1774 // check return value 1775 VERIFY_IS_EQUAL(info, 1); 1776 VERIFY_IS_EQUAL(lm.nfev, 18); 1777 VERIFY_IS_EQUAL(lm.njev, 15); 1778 // check norm^2 1779 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1780 // check x 1781 VERIFY_IS_APPROX(x[0], 1.5543827178); 1782 VERIFY_IS_APPROX(x[1], 4.0888321754); 1783 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1784 1785 /* 1786 * Second try 1787 */ 1788 x<< 1.5, 5., 450.; 1789 // do the computation 1790 info = lm.minimize(x); 1791 1792 // check return value 1793 VERIFY_IS_EQUAL(info, 1); 1794 VERIFY_IS_EQUAL(lm.nfev, 7); 1795 VERIFY_IS_EQUAL(lm.njev, 6); 1796 // check norm^2 1797 VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1798 // check x 1799 VERIFY_IS_APPROX(x[0], 1.5543827178); 1800 VERIFY_IS_APPROX(x[1], 4.0888321754); 1801 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1802 } 1803 1804 void test_NonLinearOptimization() 1805 { 1806 // Tests using the examples provided by (c)minpack 1807 CALL_SUBTEST/*_1*/(testChkder()); 1808 CALL_SUBTEST/*_1*/(testLmder1()); 1809 CALL_SUBTEST/*_1*/(testLmder()); 1810 CALL_SUBTEST/*_2*/(testHybrj1()); 1811 CALL_SUBTEST/*_2*/(testHybrj()); 1812 CALL_SUBTEST/*_2*/(testHybrd1()); 1813 CALL_SUBTEST/*_2*/(testHybrd()); 1814 CALL_SUBTEST/*_3*/(testLmstr1()); 1815 CALL_SUBTEST/*_3*/(testLmstr()); 1816 CALL_SUBTEST/*_3*/(testLmdif1()); 1817 CALL_SUBTEST/*_3*/(testLmdif()); 1818 1819 // NIST tests, level of difficulty = "Lower" 1820 CALL_SUBTEST/*_4*/(testNistMisra1a()); 1821 CALL_SUBTEST/*_4*/(testNistChwirut2()); 1822 1823 // NIST tests, level of difficulty = "Average" 1824 CALL_SUBTEST/*_5*/(testNistHahn1()); 1825 CALL_SUBTEST/*_6*/(testNistMisra1d()); 1826 CALL_SUBTEST/*_7*/(testNistMGH17()); 1827 CALL_SUBTEST/*_8*/(testNistLanczos1()); 1828 1829 // NIST tests, level of difficulty = "Higher" 1830 CALL_SUBTEST/*_9*/(testNistRat42()); 1831 CALL_SUBTEST/*_10*/(testNistMGH10()); 1832 CALL_SUBTEST/*_11*/(testNistBoxBOD()); 1833 CALL_SUBTEST/*_12*/(testNistMGH09()); 1834 CALL_SUBTEST/*_13*/(testNistBennett5()); 1835 CALL_SUBTEST/*_14*/(testNistThurber()); 1836 CALL_SUBTEST/*_15*/(testNistRat43()); 1837 CALL_SUBTEST/*_16*/(testNistEckerle4()); 1838 } 1839 1840 /* 1841 * Can be useful for debugging... 1842 printf("info, nfev : %d, %d\n", info, lm.nfev); 1843 printf("info, nfev, njev : %d, %d, %d\n", info, solver.nfev, solver.njev); 1844 printf("info, nfev : %d, %d\n", info, solver.nfev); 1845 printf("x[0] : %.32g\n", x[0]); 1846 printf("x[1] : %.32g\n", x[1]); 1847 printf("x[2] : %.32g\n", x[2]); 1848 printf("x[3] : %.32g\n", x[3]); 1849 printf("fvec.blueNorm() : %.32g\n", solver.fvec.blueNorm()); 1850 printf("fvec.blueNorm() : %.32g\n", lm.fvec.blueNorm()); 1851 1852 printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev); 1853 printf("fvec.squaredNorm() : %.13g\n", lm.fvec.squaredNorm()); 1854 std::cout << x << std::endl; 1855 std::cout.precision(9); 1856 std::cout << x[0] << std::endl; 1857 std::cout << x[1] << std::endl; 1858 std::cout << x[2] << std::endl; 1859 std::cout << x[3] << std::endl; 1860 */ 1861 1862
