17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 desire Nuentsa <desire.nuentsa_wakam@inria.fr 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// FIXME: These tests all check for hard-coded values. Ideally, parameters and start estimates should be randomized. 132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <stdio.h> 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include "main.h" 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <unsupported/Eigen/LevenbergMarquardt> 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This disables some useless Warnings on MSVC. 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// It is intended to be done for this test only. 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <Eigen/src/Core/util/DisableStupidWarnings.h> 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezusing std::sqrt; 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// tolerance for chekcing number of iterations 272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#define LM_EVAL_COUNT_TOL 4/3 282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct lmder_functor : DenseFunctor<double> 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lmder_functor(void): DenseFunctor<double>(3,15) {} 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &x, VectorXd &fvec) const 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double tmp1, tmp2, tmp3; 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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, 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (int i = 0; i < values(); i++) 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp1 = i+1; 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp2 = 16 - i - 1; 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp3 = (i>=8)? tmp2 : tmp1; 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &x, MatrixXd &fjac) const 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double tmp1, tmp2, tmp3, tmp4; 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (int i = 0; i < values(); i++) 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp1 = i+1; 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp2 = 16 - i - 1; 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp3 = (i>=8)? tmp2 : tmp1; 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = -1; 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = tmp1*tmp2/tmp4; 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = tmp1*tmp3/tmp4; 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testLmder1() 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int n=3, info; 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x; 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* the following starting values provide a rough fit. */ 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x.setConstant(n, 1.); 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lmder_functor functor; 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<lmder_functor> lm(functor); 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.lmder1(x); 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 6); 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 5); 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().blueNorm(), 0.09063596); 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x_ref(n); 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x_ref << 0.08241058, 1.133037, 2.343695; 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x, x_ref); 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testLmder() 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int m=15, n=3; 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double fnorm, covfac; 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x; 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* the following starting values provide a rough fit. */ 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x.setConstant(n, 1.); 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lmder_functor functor; 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<lmder_functor> lm(functor); 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return values 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 6); 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 5); 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fnorm = lm.fvec().blueNorm(); 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(fnorm, 0.09063596); 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x_ref(n); 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x_ref << 0.08241058, 1.133037, 2.343695; 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x, x_ref); 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check covariance 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez covfac = fnorm*fnorm/(m-n); 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixXd cov_ref(n,n); 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cov_ref << 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 0.0001531202, 0.002869941, -0.002656662, 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 0.002869941, 0.09480935, -0.09098995, 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez -0.002656662, -0.09098995, 0.08778727; 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// std::cout << fjac*covfac << std::endl; 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixXd cov; 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cov = covfac*lm.matrixR().topLeftCorner<n,n>(); 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX( cov, cov_ref); 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // TODO: why isn't this allowed ? : 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct lmdif_functor : DenseFunctor<double> 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lmdif_functor(void) : DenseFunctor<double>(3,15) {} 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &x, VectorXd &fvec) const 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int i; 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double tmp1,tmp2,tmp3; 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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, 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(x.size()==3); 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==15); 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (i=0; i<15; i++) 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp1 = i+1; 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp2 = 15 - i; 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp3 = tmp1; 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (i >= 8) tmp3 = tmp2; 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testLmdif1() 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=3; 1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n), fvec(15); 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* the following starting values provide a rough fit. */ 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x.setConstant(n, 1.); 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lmdif_functor functor; 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez DenseIndex nfev; 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev); 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// VERIFY_IS_EQUAL(nfev, 26); 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez functor(x, fvec); 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x_ref(n); 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x_ref << 0.0824106, 1.1330366, 2.3436947; 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x, x_ref); 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testLmdif() 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int m=15, n=3; 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double fnorm, covfac; 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* the following starting values provide a rough fit. */ 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x.setConstant(n, 1.); 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lmdif_functor functor; 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez NumericalDiff<lmdif_functor> numDiff(functor); 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff); 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return values 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// VERIFY_IS_EQUAL(lm.nfev(), 26); 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fnorm = lm.fvec().blueNorm(); 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(fnorm, 0.09063596); 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x_ref(n); 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x_ref << 0.08241058, 1.133037, 2.343695; 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x, x_ref); 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check covariance 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez covfac = fnorm*fnorm/(m-n); 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixXd cov_ref(n,n); 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cov_ref << 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 0.0001531202, 0.002869942, -0.002656662, 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 0.002869942, 0.09480937, -0.09098997, 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez -0.002656662, -0.09098997, 0.08778729; 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// std::cout << fjac*covfac << std::endl; 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixXd cov; 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez cov = covfac*lm.matrixR().topLeftCorner<n,n>(); 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX( cov, cov_ref); 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // TODO: why isn't this allowed ? : 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct chwirut2_functor : DenseFunctor<double> 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez chwirut2_functor(void) : DenseFunctor<double>(3,54) {} 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double m_x[54]; 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double m_y[54]; 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int i; 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==54); 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(i=0; i<54; i++) { 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x = m_x[i]; 2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==54); 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==3); 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<54; i++) { 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x = m_x[i]; 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double factor = 1./(b[1]+b[2]*x); 2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double e = exp(-b[0]*x); 2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = -x*e*factor; 2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = -e*factor*factor; 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = -x*e*factor*factor; 2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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}; 2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistChwirut2(void) 2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=3; 2842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LevenbergMarquardtSpace::Status info; 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 0.1, 0.01, 0.02; 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez chwirut2_functor functor; 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<chwirut2_functor> lm(functor); 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// VERIFY_IS_EQUAL(lm.nfev(), 10); 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 8); 3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); 3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 3067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 3077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 3097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 3117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 0.15, 0.008, 0.010; 3127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 3137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.resetParameters(); 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(1.E6*NumTraits<double>::epsilon()); 3157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol(1.E6*NumTraits<double>::epsilon()); 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 3197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// VERIFY_IS_EQUAL(lm.nfev(), 7); 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 6); 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 3267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 3287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct misra1a_functor : DenseFunctor<double> 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez misra1a_functor(void) : DenseFunctor<double>(2,14) {} 3347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double m_x[14]; 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double m_y[14]; 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==2); 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==14); 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<14; i++) { 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==2); 3487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==14); 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==2); 3507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<14; i++) { 3517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = (1.-exp(-b[1]*m_x[i])); 3527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); 3537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 3557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 3577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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}; 3587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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}; 3597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml 3617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistMisra1a(void) 3627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=2; 3647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 3657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 3677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 3697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 3707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 3717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 500., 0.0001; 3727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 3737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez misra1a_functor functor; 3747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<misra1a_functor> lm(functor); 3757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 3767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 3787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 3797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 19); 3807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 15); 3817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 3827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); 3837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 3847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 3857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 3907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 250., 0.0005; 3917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 3937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 3957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 3967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 5); 3977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 4); 3987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 3997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 4017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 4027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 4037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 4047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct hahn1_functor : DenseFunctor<double> 4067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 4077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez hahn1_functor(void) : DenseFunctor<double>(7,236) {} 4087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double m_x[236]; 4097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 4107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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 , 4127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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 , 4137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez12.596E0 , 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez13.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 }; 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==7); 4197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==236); 4207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<236; i++) { 4217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x=m_x[i], xx=x*x, xxx=xx*x; 4227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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]; 4237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 4257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 4287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==7); 4307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==236); 4317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==7); 4327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<236; i++) { 4337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x=m_x[i], xx=x*x, xxx=xx*x; 4347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 4357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = 1.*fact; 4367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = x*fact; 4377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = xx*fact; 4387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,3) = xxx*fact; 4397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 4407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,4) = x*fact; 4417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,5) = xx*fact; 4427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,6) = xxx*fact; 4437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 4457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 4477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 , 4487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez282.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 , 4497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez141.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}; 4507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml 4527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistHahn1(void) 4537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 4547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=7; 4557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 4567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 4587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 4607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 4617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 4627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 10., -1., .05, -.00001, -.05, .001, -.000001; 4637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 4647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez hahn1_functor functor; 4657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<hahn1_functor> lm(functor); 4667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 4677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 4697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 4707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 11); 4717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 10); 4727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 4737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); 4747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 4757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.0776351733E+00); 4767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1],-1.2269296921E-01); 4777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 4.0863750610E-03); 4787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 4797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 4807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[5], 2.4053735503E-04); 4817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[6],-1.2314450199E-07); 4827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 4847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 4857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 4867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; 4877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 4887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 4897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 4917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 4927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// VERIFY_IS_EQUAL(lm.nfev(), 11); 4937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 10); 4947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 4957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); 4967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 4977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 4987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 4997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 5007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 5017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 5027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 5037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 5047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 5067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct misra1d_functor : DenseFunctor<double> 5087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 5097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez misra1d_functor(void) : DenseFunctor<double>(2,14) {} 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[14]; 5117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[14]; 5127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 5137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==2); 5157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==14); 5167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<14; i++) { 5177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; 5187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 5207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 5227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==2); 5247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==14); 5257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==2); 5267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<14; i++) { 5277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double den = 1.+b[1]*x[i]; 5287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = b[1]*x[i] / den; 5297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; 5307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 5327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 5347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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}; 5357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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}; 5367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml 5387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistMisra1d(void) 5397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 5407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=2; 5417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 5427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 5447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 5467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 5477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 500., 0.0001; 5497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 5507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez misra1d_functor functor; 5517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<misra1d_functor> lm(functor); 5527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 5537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 5557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 5567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 9); 5577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 7); 5587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 5597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); 5607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 5617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 5627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 5637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 5657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 5667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 450., 0.0003; 5687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 5697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 5707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 5727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 5737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 4); 5747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 3); 5757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 5767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); 5777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 5787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 5797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 5807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 5817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct lanczos1_functor : DenseFunctor<double> 5847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 5857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lanczos1_functor(void) : DenseFunctor<double>(6,24) {} 5867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[24]; 5877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[24]; 5887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 5897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==6); 5917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==24); 5927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<24; i++) 5937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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]; 5947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 5957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 5977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==6); 5997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==24); 6007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==6); 6017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<24; i++) { 6027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = exp(-b[1]*x[i]); 6037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); 6047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = exp(-b[3]*x[i]); 6057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); 6067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,4) = exp(-b[5]*x[i]); 6077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); 6087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 6107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 6127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 6137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 6147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml 6167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistLanczos1(void) 6177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 6187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=6; 6192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LevenbergMarquardtSpace::Status info; 6207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 6227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 6247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 6257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 6267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; 6277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 6287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lanczos1_functor functor; 6297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<lanczos1_functor> lm(functor); 6307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 6317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 6332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall); 6347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 79); 6357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 72); 6367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 6372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25); 6387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 6397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 6407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 6417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 6427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 6437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 6447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 6457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 6477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 6487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 6497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; 6507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 6517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 6527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 6542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall); 6557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 9); 6567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 8); 6577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 6582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25); 6597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 6607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 6617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 6627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 6637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 6647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 6657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 6667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 6687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct rat42_functor : DenseFunctor<double> 6707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 6717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rat42_functor(void) : DenseFunctor<double>(3,9) {} 6727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[9]; 6737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[9]; 6747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 6757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 6777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==9); 6787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<9; i++) { 6797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; 6807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 6827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 6857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 6877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==9); 6887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==3); 6897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<9; i++) { 6907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double e = exp(b[1]-b[2]*x[i]); 6917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = 1./(1.+e); 6927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); 6937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); 6947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 6967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 6987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; 6997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; 7007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml 7027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistRat42(void) 7037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=3; 7052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LevenbergMarquardtSpace::Status info; 7067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 7087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 7107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 7117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 7127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 100., 1., 0.1; 7137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 7147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rat42_functor functor; 7157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<rat42_functor> lm(functor); 7167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 7177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 7192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 7207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 10); 7217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 8); 7227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 7237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); 7247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 7257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 7267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 7277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 7287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 7307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 7317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 7327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 75., 2.5, 0.07; 7337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 7347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 7357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 7372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 7387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 6); 7397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 5); 7407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 7417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); 7427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 7437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 7447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 7457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 7467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 7477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct MGH10_functor : DenseFunctor<double> 7497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MGH10_functor(void) : DenseFunctor<double>(3,16) {} 7517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[16]; 7527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[16]; 7537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 7547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 7567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==16); 7577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<16; i++) 7587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; 7597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 7607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 7627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 7647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==16); 7657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==3); 7667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<16; i++) { 7677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double factor = 1./(x[i]+b[2]); 7687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double e = exp(b[1]*factor); 7697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = e; 7707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = b[0]*factor*e; 7717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = -b[1]*b[0]*factor*factor*e; 7727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 7747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 7777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 7787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml 7807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistMGH10(void) 7817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=3; 7832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang LevenbergMarquardtSpace::Status info; 7847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 7867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 7887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 7897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 7907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 2., 400000., 25000.; 7917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 7927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MGH10_functor functor; 7937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<MGH10_functor> lm(functor); 7947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 7952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ++g_test_level; 7962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 7972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang --g_test_level; 7982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // was: VERIFY_IS_EQUAL(info, 1); 7997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 8017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); 8027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 8037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 8047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 8057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 8062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // check return value 8082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ++g_test_level; 8102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(lm.nfev(), 284 ); 8112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(lm.njev(), 249 ); 8122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang --g_test_level; 8132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.nfev() < 284 * LM_EVAL_COUNT_TOL); 8142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.njev() < 249 * LM_EVAL_COUNT_TOL); 8157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 8177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 8187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 8197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 0.02, 4000., 250.; 8207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 8217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 8222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ++g_test_level; 8232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall); 8242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // was: VERIFY_IS_EQUAL(info, 1); 8252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang --g_test_level; 8267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 8287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); 8297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 8307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 8317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 8327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 8332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // check return value 8352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ++g_test_level; 8362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(lm.nfev(), 126); 8372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(lm.njev(), 116); 8382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang --g_test_level; 8392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.nfev() < 126 * LM_EVAL_COUNT_TOL); 8402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.njev() < 116 * LM_EVAL_COUNT_TOL); 8417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 8427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct BoxBOD_functor : DenseFunctor<double> 8457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 8467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez BoxBOD_functor(void) : DenseFunctor<double>(2,6) {} 8477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[6]; 8487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 8497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 8507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[6] = { 109., 149., 149., 191., 213., 224. }; 8517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==2); 8527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==6); 8537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<6; i++) 8547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; 8557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 8567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 8577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 8587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 8597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==2); 8607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==6); 8617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==2); 8627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<6; i++) { 8637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double e = exp(-b[1]*x[i]); 8647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = 1.-e; 8657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = b[0]*x[i]*e; 8667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 8677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 8687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 8697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 8707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; 8717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml 8737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistBoxBOD(void) 8747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 8757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=2; 8767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 8777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 8797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 8817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 8827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 8837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 1., 1.; 8847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 8857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez BoxBOD_functor functor; 8867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<BoxBOD_functor> lm(functor); 8877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(1.E6*NumTraits<double>::epsilon()); 8887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol(1.E6*NumTraits<double>::epsilon()); 8897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFactor(10); 8907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 8917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 8937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); 8947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 8957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 8967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 8972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // check return value 8992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, 1); 9002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.nfev() < 31); // 31 9012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.njev() < 25); // 25 9027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 9047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 9057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 9067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 100., 0.75; 9077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 9087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.resetParameters(); 9097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(NumTraits<double>::epsilon()); 9107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol( NumTraits<double>::epsilon()); 9117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 9127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 9147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 9152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ++g_test_level; 9162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(lm.nfev(), 16 ); 9172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(lm.njev(), 15 ); 9182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang --g_test_level; 9192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.nfev() < 16 * LM_EVAL_COUNT_TOL); 9202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.njev() < 15 * LM_EVAL_COUNT_TOL); 9217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 9227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); 9237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 9247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 9257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 9267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 9277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct MGH17_functor : DenseFunctor<double> 9297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 9307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MGH17_functor(void) : DenseFunctor<double>(5,33) {} 9317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[33]; 9327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[33]; 9337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 9347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 9357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==5); 9367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==33); 9377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<33; i++) 9387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; 9397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 9407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 9417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 9427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 9437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==5); 9447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==33); 9457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==5); 9467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<33; i++) { 9477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = 1.; 9487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = exp(-b[3]*x[i]); 9497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = exp(-b[4]*x[i]); 9507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); 9517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); 9527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 9537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 9547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 9557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 9567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 9577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 9587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml 9607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistMGH17(void) 9617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 9627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=5; 9637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 9647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 9667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 9687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 9697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 9707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 50., 150., -100., 1., 2.; 9717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 9727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MGH17_functor functor; 9737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<MGH17_functor> lm(functor); 9747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(NumTraits<double>::epsilon()); 9757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol(NumTraits<double>::epsilon()); 9767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setMaxfev(1000); 9777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 9787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 9807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); 9817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 9827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 9837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 9847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 9857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 9867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 9872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // check return value 9892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success) 9902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.nfev() < 700 ); // 602 9912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.njev() < 600 ); // 545 9927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 9947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 9957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 9967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; 9977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 9987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.resetParameters(); 9997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 10007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 10027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 10037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 18); 10047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 15); 10057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 10067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); 10077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 10087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 10097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 10107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 10117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 10127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 10137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 10147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct MGH09_functor : DenseFunctor<double> 10167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 10177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MGH09_functor(void) : DenseFunctor<double>(4,11) {} 10187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double _x[11]; 10197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[11]; 10207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 10217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 10227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==4); 10237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==11); 10247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<11; i++) { 10257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x = _x[i], xx=x*x; 10267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; 10277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 10287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 10297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 10307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 10317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 10327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==4); 10337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==11); 10347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==4); 10357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<11; i++) { 10367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x = _x[i], xx=x*x; 10377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double factor = 1./(xx+x*b[2]+b[3]); 10387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = (xx+x*b[1]) * factor; 10397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = b[0]*x* factor; 10407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; 10417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; 10427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 10437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 10447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 10457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 10467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 10477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 10487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml 10507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistMGH09(void) 10517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 10527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=4; 10537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 10547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 10567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 10587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 10597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 10607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 25., 39, 41.5, 39.; 10617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 10627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MGH09_functor functor; 10637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<MGH09_functor> lm(functor); 10647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setMaxfev(1000); 10657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 10667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 10687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); 10697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 10707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 10717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 10727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 10737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 10742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // check return value 10752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_EQUAL(info, 1); 10762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.nfev() < 510 ); // 490 10772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY(lm.njev() < 400 ); // 376 10787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 10807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 10817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 10827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 0.25, 0.39, 0.415, 0.39; 10837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 10847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.resetParameters(); 10857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 10867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 10877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 10887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 10897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 18); 10907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 16); 10917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 10927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); 10937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 10947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 10957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 10967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 10977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 10987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 10997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct Bennett5_functor : DenseFunctor<double> 11037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 11047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Bennett5_functor(void) : DenseFunctor<double>(3,154) {} 11057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[154]; 11067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[154]; 11077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 11087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 11097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 11107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==154); 11117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<154; i++) 11127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; 11137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 11147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 11157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 11167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 11177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 11187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==154); 11197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==3); 11207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<154; i++) { 11217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double e = pow(b[1]+x[i],-1./b[2]); 11227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = e; 11237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); 11247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; 11257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 11267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 11277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 11287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 11297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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, 11307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez11.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 }; 11317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 11327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez,-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 , 11337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; 11347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml 11367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistBennett5(void) 11377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 11387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=3; 11397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 11407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 11427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 11447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 11457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 11467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< -2000., 50., 0.8; 11477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 11487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Bennett5_functor functor; 11497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<Bennett5_functor> lm(functor); 11507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setMaxfev(1000); 11517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 11527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 11547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 11557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 758); 11567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 744); 11577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 11587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); 11597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 11607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], -2.5235058043E+03); 11617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 4.6736564644E+01); 11627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 9.3218483193E-01); 11637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 11647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 11657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 11667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< -1500., 45., 0.85; 11677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 11687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.resetParameters(); 11697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 11707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 11727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 11737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 203); 11747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 192); 11757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 11767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); 11777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 11787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 11797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); 11807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); 11817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 11827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 11837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct thurber_functor : DenseFunctor<double> 11847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 11857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez thurber_functor(void) : DenseFunctor<double>(7,37) {} 11867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double _x[37]; 11877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double _y[37]; 11887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 11897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 11907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 11917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==7); 11927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==37); 11937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<37; i++) { 11947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x=_x[i], xx=x*x, xxx=xx*x; 11957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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]; 11967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 11977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 11987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 11997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 12007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 12017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==7); 12027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==37); 12037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==7); 12047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<37; i++) { 12057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double x=_x[i], xx=x*x, xxx=xx*x; 12067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 12077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = 1.*fact; 12087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = x*fact; 12097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = xx*fact; 12107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,3) = xxx*fact; 12117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 12127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,4) = x*fact; 12137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,5) = xx*fact; 12147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,6) = xxx*fact; 12157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 12167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 12177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 12187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 12197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 12207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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}; 12217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 12227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml 12237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistThurber(void) 12247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 12257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=7; 12267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 12277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 12287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 12297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 12307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 12317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 12327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 12337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; 12347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 12357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez thurber_functor functor; 12367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<thurber_functor> lm(functor); 12377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(1.E4*NumTraits<double>::epsilon()); 12387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol(1.E4*NumTraits<double>::epsilon()); 12397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 12407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 12417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 12427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 12437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 39); 12447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 36); 12457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 12467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); 12477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 12487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 12497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 12507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 12517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 12527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 12537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 12547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 12557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 12567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 12577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 12587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 12597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; 12607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 12617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.resetParameters(); 12627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(1.E4*NumTraits<double>::epsilon()); 12637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol(1.E4*NumTraits<double>::epsilon()); 12647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 12657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 12667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 12677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 12687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 29); 12697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 28); 12707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 12717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); 12727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 12737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 12747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 12757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 12767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 12777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 12787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 12797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 12807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 12817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 12827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct rat43_functor : DenseFunctor<double> 12837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 12847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rat43_functor(void) : DenseFunctor<double>(4,15) {} 12857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[15]; 12867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[15]; 12877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 12887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 12897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==4); 12907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==15); 12917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<15; i++) 12927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; 12937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 12947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 12957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 12967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 12977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==4); 12987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==15); 12997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==4); 13007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<15; i++) { 13017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double e = exp(b[1]-b[2]*x[i]); 13027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double power = -1./b[3]; 13037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = pow(1.+e, power); 13047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); 13057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); 13067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); 13077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 13087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 13097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 13107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 13117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; 13127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 13137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml 13157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistRat43(void) 13167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 13177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=4; 13187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 13197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 13217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 13237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 13247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 13257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 100., 10., 1., 1.; 13267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 13277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rat43_functor functor; 13287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<rat43_functor> lm(functor); 13297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(1.E6*NumTraits<double>::epsilon()); 13307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol(1.E6*NumTraits<double>::epsilon()); 13317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 13327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 13347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 13357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 27); 13367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 20); 13377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 13387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); 13397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 13407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 13417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 13427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 13437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 13447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 13467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 13477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 13487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 700., 5., 0.75, 1.3; 13497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 13507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.resetParameters(); 13517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setFtol(1.E5*NumTraits<double>::epsilon()); 13527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lm.setXtol(1.E5*NumTraits<double>::epsilon()); 13537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 13547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 13567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 13577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 9); 13587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 8); 13597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 13607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); 13617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 13627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 13637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 13647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 13657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 13667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 13677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 13707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct eckerle4_functor : DenseFunctor<double> 13717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 13727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eckerle4_functor(void) : DenseFunctor<double>(3,35) {} 13737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double x[35]; 13747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const double y[35]; 13757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int operator()(const VectorXd &b, VectorXd &fvec) 13767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 13777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 13787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fvec.size()==35); 13797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<35; i++) 13807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; 13817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 13827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 13837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int df(const VectorXd &b, MatrixXd &fjac) 13847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 13857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(b.size()==3); 13867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.rows()==35); 13877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez assert(fjac.cols()==3); 13887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<35; i++) { 13897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double b12 = b[1]*b[1]; 13907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); 13917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,0) = e / b[1]; 13927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; 13937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; 13947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 13957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return 0; 13967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 13977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 13987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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}; 13997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 }; 14007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml 14027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid testNistEckerle4(void) 14037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 14047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const int n=3; 14057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int info; 14067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorXd x(n); 14087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 14107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * First try 14117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 14127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 1., 10., 500.; 14137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 14147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eckerle4_functor functor; 14157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LevenbergMarquardt<eckerle4_functor> lm(functor); 14167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 14177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 14197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 14207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 18); 14217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 15); 14227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 14237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); 14247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 14257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.5543827178); 14267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 4.0888321754); 14277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 14287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* 14307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Second try 14317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 14327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x<< 1.5, 5., 450.; 14337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do the computation 14347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = lm.minimize(x); 14357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check return value 14377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(info, 1); 14387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.nfev(), 7); 14397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(lm.njev(), 6); 14407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check norm^2 14417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); 14427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check x 14437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[0], 1.5543827178); 14447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[1], 4.0888321754); 14457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 14467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 14477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid test_levenberg_marquardt() 14497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 14507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Tests using the examples provided by (c)minpack 14517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testLmder1()); 14527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testLmder()); 14537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testLmdif1()); 14547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// CALL_SUBTEST(testLmstr1()); 14557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// CALL_SUBTEST(testLmstr()); 14567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testLmdif()); 14577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // NIST tests, level of difficulty = "Lower" 14597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistMisra1a()); 14607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistChwirut2()); 14617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // NIST tests, level of difficulty = "Average" 14637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistHahn1()); 14647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistMisra1d()); 14657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistMGH17()); 14667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistLanczos1()); 14677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 14687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// // NIST tests, level of difficulty = "Higher" 14697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistRat42()); 14707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistMGH10()); 14717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistBoxBOD()); 14727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// CALL_SUBTEST(testNistMGH09()); 14737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistBennett5()); 14747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistThurber()); 14757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistRat43()); 14767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(testNistEckerle4()); 14777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1478