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