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