1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Thomas Capricelli <orzel (at) freehackers.org> 5 // Copyright (C) 2012 desire Nuentsa <desire.nuentsa_wakam (at) inria.fr 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 12 #include <stdio.h> 13 14 #include "main.h" 15 #include <unsupported/Eigen/LevenbergMarquardt> 16 17 // This disables some useless Warnings on MSVC. 18 // It is intended to be done for this test only. 19 #include <Eigen/src/Core/util/DisableStupidWarnings.h> 20 21 using std::sqrt; 22 23 struct lmder_functor : DenseFunctor<double> 24 { 25 lmder_functor(void): DenseFunctor<double>(3,15) {} 26 int operator()(const VectorXd &x, VectorXd &fvec) const 27 { 28 double tmp1, tmp2, tmp3; 29 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, 30 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 31 32 for (int i = 0; i < values(); i++) 33 { 34 tmp1 = i+1; 35 tmp2 = 16 - i - 1; 36 tmp3 = (i>=8)? tmp2 : tmp1; 37 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 38 } 39 return 0; 40 } 41 42 int df(const VectorXd &x, MatrixXd &fjac) const 43 { 44 double tmp1, tmp2, tmp3, tmp4; 45 for (int i = 0; i < values(); i++) 46 { 47 tmp1 = i+1; 48 tmp2 = 16 - i - 1; 49 tmp3 = (i>=8)? tmp2 : tmp1; 50 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 51 fjac(i,0) = -1; 52 fjac(i,1) = tmp1*tmp2/tmp4; 53 fjac(i,2) = tmp1*tmp3/tmp4; 54 } 55 return 0; 56 } 57 }; 58 59 void testLmder1() 60 { 61 int n=3, info; 62 63 VectorXd x; 64 65 /* the following starting values provide a rough fit. */ 66 x.setConstant(n, 1.); 67 68 // do the computation 69 lmder_functor functor; 70 LevenbergMarquardt<lmder_functor> lm(functor); 71 info = lm.lmder1(x); 72 73 // check return value 74 VERIFY_IS_EQUAL(info, 1); 75 VERIFY_IS_EQUAL(lm.nfev(), 6); 76 VERIFY_IS_EQUAL(lm.njev(), 5); 77 78 // check norm 79 VERIFY_IS_APPROX(lm.fvec().blueNorm(), 0.09063596); 80 81 // check x 82 VectorXd x_ref(n); 83 x_ref << 0.08241058, 1.133037, 2.343695; 84 VERIFY_IS_APPROX(x, x_ref); 85 } 86 87 void testLmder() 88 { 89 const int m=15, n=3; 90 int info; 91 double fnorm, covfac; 92 VectorXd x; 93 94 /* the following starting values provide a rough fit. */ 95 x.setConstant(n, 1.); 96 97 // do the computation 98 lmder_functor functor; 99 LevenbergMarquardt<lmder_functor> lm(functor); 100 info = lm.minimize(x); 101 102 // check return values 103 VERIFY_IS_EQUAL(info, 1); 104 VERIFY_IS_EQUAL(lm.nfev(), 6); 105 VERIFY_IS_EQUAL(lm.njev(), 5); 106 107 // check norm 108 fnorm = lm.fvec().blueNorm(); 109 VERIFY_IS_APPROX(fnorm, 0.09063596); 110 111 // check x 112 VectorXd x_ref(n); 113 x_ref << 0.08241058, 1.133037, 2.343695; 114 VERIFY_IS_APPROX(x, x_ref); 115 116 // check covariance 117 covfac = fnorm*fnorm/(m-n); 118 internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm 119 120 MatrixXd cov_ref(n,n); 121 cov_ref << 122 0.0001531202, 0.002869941, -0.002656662, 123 0.002869941, 0.09480935, -0.09098995, 124 -0.002656662, -0.09098995, 0.08778727; 125 126 // std::cout << fjac*covfac << std::endl; 127 128 MatrixXd cov; 129 cov = covfac*lm.matrixR().topLeftCorner<n,n>(); 130 VERIFY_IS_APPROX( cov, cov_ref); 131 // TODO: why isn't this allowed ? : 132 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 133 } 134 135 struct lmdif_functor : DenseFunctor<double> 136 { 137 lmdif_functor(void) : DenseFunctor<double>(3,15) {} 138 int operator()(const VectorXd &x, VectorXd &fvec) const 139 { 140 int i; 141 double tmp1,tmp2,tmp3; 142 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, 143 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; 144 145 assert(x.size()==3); 146 assert(fvec.size()==15); 147 for (i=0; i<15; i++) 148 { 149 tmp1 = i+1; 150 tmp2 = 15 - i; 151 tmp3 = tmp1; 152 153 if (i >= 8) tmp3 = tmp2; 154 fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 155 } 156 return 0; 157 } 158 }; 159 160 void testLmdif1() 161 { 162 const int n=3; 163 int info; 164 165 VectorXd x(n), fvec(15); 166 167 /* the following starting values provide a rough fit. */ 168 x.setConstant(n, 1.); 169 170 // do the computation 171 lmdif_functor functor; 172 DenseIndex nfev; 173 info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev); 174 175 // check return value 176 VERIFY_IS_EQUAL(info, 1); 177 // VERIFY_IS_EQUAL(nfev, 26); 178 179 // check norm 180 functor(x, fvec); 181 VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); 182 183 // check x 184 VectorXd x_ref(n); 185 x_ref << 0.0824106, 1.1330366, 2.3436947; 186 VERIFY_IS_APPROX(x, x_ref); 187 188 } 189 190 void testLmdif() 191 { 192 const int m=15, n=3; 193 int info; 194 double fnorm, covfac; 195 VectorXd x(n); 196 197 /* the following starting values provide a rough fit. */ 198 x.setConstant(n, 1.); 199 200 // do the computation 201 lmdif_functor functor; 202 NumericalDiff<lmdif_functor> numDiff(functor); 203 LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff); 204 info = lm.minimize(x); 205 206 // check return values 207 VERIFY_IS_EQUAL(info, 1); 208 // VERIFY_IS_EQUAL(lm.nfev(), 26); 209 210 // check norm 211 fnorm = lm.fvec().blueNorm(); 212 VERIFY_IS_APPROX(fnorm, 0.09063596); 213 214 // check x 215 VectorXd x_ref(n); 216 x_ref << 0.08241058, 1.133037, 2.343695; 217 VERIFY_IS_APPROX(x, x_ref); 218 219 // check covariance 220 covfac = fnorm*fnorm/(m-n); 221 internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm 222 223 MatrixXd cov_ref(n,n); 224 cov_ref << 225 0.0001531202, 0.002869942, -0.002656662, 226 0.002869942, 0.09480937, -0.09098997, 227 -0.002656662, -0.09098997, 0.08778729; 228 229 // std::cout << fjac*covfac << std::endl; 230 231 MatrixXd cov; 232 cov = covfac*lm.matrixR().topLeftCorner<n,n>(); 233 VERIFY_IS_APPROX( cov, cov_ref); 234 // TODO: why isn't this allowed ? : 235 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 236 } 237 238 struct chwirut2_functor : DenseFunctor<double> 239 { 240 chwirut2_functor(void) : DenseFunctor<double>(3,54) {} 241 static const double m_x[54]; 242 static const double m_y[54]; 243 int operator()(const VectorXd &b, VectorXd &fvec) 244 { 245 int i; 246 247 assert(b.size()==3); 248 assert(fvec.size()==54); 249 for(i=0; i<54; i++) { 250 double x = m_x[i]; 251 fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; 252 } 253 return 0; 254 } 255 int df(const VectorXd &b, MatrixXd &fjac) 256 { 257 assert(b.size()==3); 258 assert(fjac.rows()==54); 259 assert(fjac.cols()==3); 260 for(int i=0; i<54; i++) { 261 double x = m_x[i]; 262 double factor = 1./(b[1]+b[2]*x); 263 double e = exp(-b[0]*x); 264 fjac(i,0) = -x*e*factor; 265 fjac(i,1) = -e*factor*factor; 266 fjac(i,2) = -x*e*factor*factor; 267 } 268 return 0; 269 } 270 }; 271 const double chwirut2_functor::m_x[54] = { 0.500E0, 1.000E0, 1.750E0, 3.750E0, 5.750E0, 0.875E0, 2.250E0, 3.250E0, 5.250E0, 0.750E0, 1.750E0, 2.750E0, 4.750E0, 0.625E0, 1.250E0, 2.250E0, 4.250E0, .500E0, 3.000E0, .750E0, 3.000E0, 1.500E0, 6.000E0, 3.000E0, 6.000E0, 1.500E0, 3.000E0, .500E0, 2.000E0, 4.000E0, .750E0, 2.000E0, 5.000E0, .750E0, 2.250E0, 3.750E0, 5.750E0, 3.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .500E0, 6.000E0, 3.000E0, .500E0, 2.750E0, .500E0, 1.750E0}; 272 const double chwirut2_functor::m_y[54] = { 92.9000E0 ,57.1000E0 ,31.0500E0 ,11.5875E0 ,8.0250E0 ,63.6000E0 ,21.4000E0 ,14.2500E0 ,8.4750E0 ,63.8000E0 ,26.8000E0 ,16.4625E0 ,7.1250E0 ,67.3000E0 ,41.0000E0 ,21.1500E0 ,8.1750E0 ,81.5000E0 ,13.1200E0 ,59.9000E0 ,14.6200E0 ,32.9000E0 ,5.4400E0 ,12.5600E0 ,5.4400E0 ,32.0000E0 ,13.9500E0 ,75.8000E0 ,20.0000E0 ,10.4200E0 ,59.5000E0 ,21.6700E0 ,8.5500E0 ,62.0000E0 ,20.2000E0 ,7.7600E0 ,3.7500E0 ,11.8100E0 ,54.7000E0 ,23.7000E0 ,11.5500E0 ,61.3000E0 ,17.7000E0 ,8.7400E0 ,59.2000E0 ,16.3000E0 ,8.6200E0 ,81.0000E0 ,4.8700E0 ,14.6200E0 ,81.7000E0 ,17.1700E0 ,81.3000E0 ,28.9000E0 }; 273 274 // http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml 275 void testNistChwirut2(void) 276 { 277 const int n=3; 278 int info; 279 280 VectorXd x(n); 281 282 /* 283 * First try 284 */ 285 x<< 0.1, 0.01, 0.02; 286 // do the computation 287 chwirut2_functor functor; 288 LevenbergMarquardt<chwirut2_functor> lm(functor); 289 info = lm.minimize(x); 290 291 // check return value 292 VERIFY_IS_EQUAL(info, 1); 293 // VERIFY_IS_EQUAL(lm.nfev(), 10); 294 VERIFY_IS_EQUAL(lm.njev(), 8); 295 // check norm^2 296 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); 297 // check x 298 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 299 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 300 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 301 302 /* 303 * Second try 304 */ 305 x<< 0.15, 0.008, 0.010; 306 // do the computation 307 lm.resetParameters(); 308 lm.setFtol(1.E6*NumTraits<double>::epsilon()); 309 lm.setXtol(1.E6*NumTraits<double>::epsilon()); 310 info = lm.minimize(x); 311 312 // check return value 313 VERIFY_IS_EQUAL(info, 1); 314 // VERIFY_IS_EQUAL(lm.nfev(), 7); 315 VERIFY_IS_EQUAL(lm.njev(), 6); 316 // check norm^2 317 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); 318 // check x 319 VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 320 VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 321 VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 322 } 323 324 325 struct misra1a_functor : DenseFunctor<double> 326 { 327 misra1a_functor(void) : DenseFunctor<double>(2,14) {} 328 static const double m_x[14]; 329 static const double m_y[14]; 330 int operator()(const VectorXd &b, VectorXd &fvec) 331 { 332 assert(b.size()==2); 333 assert(fvec.size()==14); 334 for(int i=0; i<14; i++) { 335 fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; 336 } 337 return 0; 338 } 339 int df(const VectorXd &b, MatrixXd &fjac) 340 { 341 assert(b.size()==2); 342 assert(fjac.rows()==14); 343 assert(fjac.cols()==2); 344 for(int i=0; i<14; i++) { 345 fjac(i,0) = (1.-exp(-b[1]*m_x[i])); 346 fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); 347 } 348 return 0; 349 } 350 }; 351 const double misra1a_functor::m_x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; 352 const double misra1a_functor::m_y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; 353 354 // http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml 355 void testNistMisra1a(void) 356 { 357 const int n=2; 358 int info; 359 360 VectorXd x(n); 361 362 /* 363 * First try 364 */ 365 x<< 500., 0.0001; 366 // do the computation 367 misra1a_functor functor; 368 LevenbergMarquardt<misra1a_functor> lm(functor); 369 info = lm.minimize(x); 370 371 // check return value 372 VERIFY_IS_EQUAL(info, 1); 373 VERIFY_IS_EQUAL(lm.nfev(), 19); 374 VERIFY_IS_EQUAL(lm.njev(), 15); 375 // check norm^2 376 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); 377 // check x 378 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 379 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 380 381 /* 382 * Second try 383 */ 384 x<< 250., 0.0005; 385 // do the computation 386 info = lm.minimize(x); 387 388 // check return value 389 VERIFY_IS_EQUAL(info, 1); 390 VERIFY_IS_EQUAL(lm.nfev(), 5); 391 VERIFY_IS_EQUAL(lm.njev(), 4); 392 // check norm^2 393 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); 394 // check x 395 VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 396 VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 397 } 398 399 struct hahn1_functor : DenseFunctor<double> 400 { 401 hahn1_functor(void) : DenseFunctor<double>(7,236) {} 402 static const double m_x[236]; 403 int operator()(const VectorXd &b, VectorXd &fvec) 404 { 405 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 , 406 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 , 407 12.596E0 , 408 13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0 , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0 , 20.935E0 , 21.035E0 , 20.93E0 , 21.074E0 , 21.085E0 , 20.935E0 }; 409 410 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 411 412 assert(b.size()==7); 413 assert(fvec.size()==236); 414 for(int i=0; i<236; i++) { 415 double x=m_x[i], xx=x*x, xxx=xx*x; 416 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]; 417 } 418 return 0; 419 } 420 421 int df(const VectorXd &b, MatrixXd &fjac) 422 { 423 assert(b.size()==7); 424 assert(fjac.rows()==236); 425 assert(fjac.cols()==7); 426 for(int i=0; i<236; i++) { 427 double x=m_x[i], xx=x*x, xxx=xx*x; 428 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 429 fjac(i,0) = 1.*fact; 430 fjac(i,1) = x*fact; 431 fjac(i,2) = xx*fact; 432 fjac(i,3) = xxx*fact; 433 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 434 fjac(i,4) = x*fact; 435 fjac(i,5) = xx*fact; 436 fjac(i,6) = xxx*fact; 437 } 438 return 0; 439 } 440 }; 441 const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , 442 282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 , 443 141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0 , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0 , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0}; 444 445 // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml 446 void testNistHahn1(void) 447 { 448 const int n=7; 449 int info; 450 451 VectorXd x(n); 452 453 /* 454 * First try 455 */ 456 x<< 10., -1., .05, -.00001, -.05, .001, -.000001; 457 // do the computation 458 hahn1_functor functor; 459 LevenbergMarquardt<hahn1_functor> lm(functor); 460 info = lm.minimize(x); 461 462 // check return value 463 VERIFY_IS_EQUAL(info, 1); 464 VERIFY_IS_EQUAL(lm.nfev(), 11); 465 VERIFY_IS_EQUAL(lm.njev(), 10); 466 // check norm^2 467 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); 468 // check x 469 VERIFY_IS_APPROX(x[0], 1.0776351733E+00); 470 VERIFY_IS_APPROX(x[1],-1.2269296921E-01); 471 VERIFY_IS_APPROX(x[2], 4.0863750610E-03); 472 VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 473 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 474 VERIFY_IS_APPROX(x[5], 2.4053735503E-04); 475 VERIFY_IS_APPROX(x[6],-1.2314450199E-07); 476 477 /* 478 * Second try 479 */ 480 x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; 481 // do the computation 482 info = lm.minimize(x); 483 484 // check return value 485 VERIFY_IS_EQUAL(info, 1); 486 // VERIFY_IS_EQUAL(lm.nfev(), 11); 487 VERIFY_IS_EQUAL(lm.njev(), 10); 488 // check norm^2 489 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); 490 // check x 491 VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 492 VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 493 VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 494 VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 495 VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 496 VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 497 VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 498 499 } 500 501 struct misra1d_functor : DenseFunctor<double> 502 { 503 misra1d_functor(void) : DenseFunctor<double>(2,14) {} 504 static const double x[14]; 505 static const double y[14]; 506 int operator()(const VectorXd &b, VectorXd &fvec) 507 { 508 assert(b.size()==2); 509 assert(fvec.size()==14); 510 for(int i=0; i<14; i++) { 511 fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; 512 } 513 return 0; 514 } 515 int df(const VectorXd &b, MatrixXd &fjac) 516 { 517 assert(b.size()==2); 518 assert(fjac.rows()==14); 519 assert(fjac.cols()==2); 520 for(int i=0; i<14; i++) { 521 double den = 1.+b[1]*x[i]; 522 fjac(i,0) = b[1]*x[i] / den; 523 fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; 524 } 525 return 0; 526 } 527 }; 528 const double misra1d_functor::x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; 529 const double misra1d_functor::y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; 530 531 // http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml 532 void testNistMisra1d(void) 533 { 534 const int n=2; 535 int info; 536 537 VectorXd x(n); 538 539 /* 540 * First try 541 */ 542 x<< 500., 0.0001; 543 // do the computation 544 misra1d_functor functor; 545 LevenbergMarquardt<misra1d_functor> lm(functor); 546 info = lm.minimize(x); 547 548 // check return value 549 VERIFY_IS_EQUAL(info, 1); 550 VERIFY_IS_EQUAL(lm.nfev(), 9); 551 VERIFY_IS_EQUAL(lm.njev(), 7); 552 // check norm^2 553 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); 554 // check x 555 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 556 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 557 558 /* 559 * Second try 560 */ 561 x<< 450., 0.0003; 562 // do the computation 563 info = lm.minimize(x); 564 565 // check return value 566 VERIFY_IS_EQUAL(info, 1); 567 VERIFY_IS_EQUAL(lm.nfev(), 4); 568 VERIFY_IS_EQUAL(lm.njev(), 3); 569 // check norm^2 570 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); 571 // check x 572 VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 573 VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 574 } 575 576 577 struct lanczos1_functor : DenseFunctor<double> 578 { 579 lanczos1_functor(void) : DenseFunctor<double>(6,24) {} 580 static const double x[24]; 581 static const double y[24]; 582 int operator()(const VectorXd &b, VectorXd &fvec) 583 { 584 assert(b.size()==6); 585 assert(fvec.size()==24); 586 for(int i=0; i<24; i++) 587 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]; 588 return 0; 589 } 590 int df(const VectorXd &b, MatrixXd &fjac) 591 { 592 assert(b.size()==6); 593 assert(fjac.rows()==24); 594 assert(fjac.cols()==6); 595 for(int i=0; i<24; i++) { 596 fjac(i,0) = exp(-b[1]*x[i]); 597 fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); 598 fjac(i,2) = exp(-b[3]*x[i]); 599 fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); 600 fjac(i,4) = exp(-b[5]*x[i]); 601 fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); 602 } 603 return 0; 604 } 605 }; 606 const double lanczos1_functor::x[24] = { 0.000000000000E+00, 5.000000000000E-02, 1.000000000000E-01, 1.500000000000E-01, 2.000000000000E-01, 2.500000000000E-01, 3.000000000000E-01, 3.500000000000E-01, 4.000000000000E-01, 4.500000000000E-01, 5.000000000000E-01, 5.500000000000E-01, 6.000000000000E-01, 6.500000000000E-01, 7.000000000000E-01, 7.500000000000E-01, 8.000000000000E-01, 8.500000000000E-01, 9.000000000000E-01, 9.500000000000E-01, 1.000000000000E+00, 1.050000000000E+00, 1.100000000000E+00, 1.150000000000E+00 }; 607 const double lanczos1_functor::y[24] = { 2.513400000000E+00 ,2.044333373291E+00 ,1.668404436564E+00 ,1.366418021208E+00 ,1.123232487372E+00 ,9.268897180037E-01 ,7.679338563728E-01 ,6.388775523106E-01 ,5.337835317402E-01 ,4.479363617347E-01 ,3.775847884350E-01 ,3.197393199326E-01 ,2.720130773746E-01 ,2.324965529032E-01 ,1.996589546065E-01 ,1.722704126914E-01 ,1.493405660168E-01 ,1.300700206922E-01 ,1.138119324644E-01 ,1.000415587559E-01 ,8.833209084540E-02 ,7.833544019350E-02 ,6.976693743449E-02 ,6.239312536719E-02 }; 608 609 // http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml 610 void testNistLanczos1(void) 611 { 612 const int n=6; 613 int info; 614 615 VectorXd x(n); 616 617 /* 618 * First try 619 */ 620 x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; 621 // do the computation 622 lanczos1_functor functor; 623 LevenbergMarquardt<lanczos1_functor> lm(functor); 624 info = lm.minimize(x); 625 626 // check return value 627 VERIFY_IS_EQUAL(info, 2); 628 VERIFY_IS_EQUAL(lm.nfev(), 79); 629 VERIFY_IS_EQUAL(lm.njev(), 72); 630 // check norm^2 631 // VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.430899764097e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 632 // check x 633 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 634 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 635 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 636 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 637 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 638 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 639 640 /* 641 * Second try 642 */ 643 x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; 644 // do the computation 645 info = lm.minimize(x); 646 647 // check return value 648 VERIFY_IS_EQUAL(info, 2); 649 VERIFY_IS_EQUAL(lm.nfev(), 9); 650 VERIFY_IS_EQUAL(lm.njev(), 8); 651 // check norm^2 652 // VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.428595533845e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 653 // check x 654 VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 655 VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 656 VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 657 VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 658 VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 659 VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 660 661 } 662 663 struct rat42_functor : DenseFunctor<double> 664 { 665 rat42_functor(void) : DenseFunctor<double>(3,9) {} 666 static const double x[9]; 667 static const double y[9]; 668 int operator()(const VectorXd &b, VectorXd &fvec) 669 { 670 assert(b.size()==3); 671 assert(fvec.size()==9); 672 for(int i=0; i<9; i++) { 673 fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; 674 } 675 return 0; 676 } 677 678 int df(const VectorXd &b, MatrixXd &fjac) 679 { 680 assert(b.size()==3); 681 assert(fjac.rows()==9); 682 assert(fjac.cols()==3); 683 for(int i=0; i<9; i++) { 684 double e = exp(b[1]-b[2]*x[i]); 685 fjac(i,0) = 1./(1.+e); 686 fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); 687 fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); 688 } 689 return 0; 690 } 691 }; 692 const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; 693 const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; 694 695 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml 696 void testNistRat42(void) 697 { 698 const int n=3; 699 int info; 700 701 VectorXd x(n); 702 703 /* 704 * First try 705 */ 706 x<< 100., 1., 0.1; 707 // do the computation 708 rat42_functor functor; 709 LevenbergMarquardt<rat42_functor> lm(functor); 710 info = lm.minimize(x); 711 712 // check return value 713 VERIFY_IS_EQUAL(info, 1); 714 VERIFY_IS_EQUAL(lm.nfev(), 10); 715 VERIFY_IS_EQUAL(lm.njev(), 8); 716 // check norm^2 717 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); 718 // check x 719 VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 720 VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 721 VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 722 723 /* 724 * Second try 725 */ 726 x<< 75., 2.5, 0.07; 727 // do the computation 728 info = lm.minimize(x); 729 730 // check return value 731 VERIFY_IS_EQUAL(info, 1); 732 VERIFY_IS_EQUAL(lm.nfev(), 6); 733 VERIFY_IS_EQUAL(lm.njev(), 5); 734 // check norm^2 735 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); 736 // check x 737 VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 738 VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 739 VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 740 } 741 742 struct MGH10_functor : DenseFunctor<double> 743 { 744 MGH10_functor(void) : DenseFunctor<double>(3,16) {} 745 static const double x[16]; 746 static const double y[16]; 747 int operator()(const VectorXd &b, VectorXd &fvec) 748 { 749 assert(b.size()==3); 750 assert(fvec.size()==16); 751 for(int i=0; i<16; i++) 752 fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; 753 return 0; 754 } 755 int df(const VectorXd &b, MatrixXd &fjac) 756 { 757 assert(b.size()==3); 758 assert(fjac.rows()==16); 759 assert(fjac.cols()==3); 760 for(int i=0; i<16; i++) { 761 double factor = 1./(x[i]+b[2]); 762 double e = exp(b[1]*factor); 763 fjac(i,0) = e; 764 fjac(i,1) = b[0]*factor*e; 765 fjac(i,2) = -b[1]*b[0]*factor*factor*e; 766 } 767 return 0; 768 } 769 }; 770 const double MGH10_functor::x[16] = { 5.000000E+01, 5.500000E+01, 6.000000E+01, 6.500000E+01, 7.000000E+01, 7.500000E+01, 8.000000E+01, 8.500000E+01, 9.000000E+01, 9.500000E+01, 1.000000E+02, 1.050000E+02, 1.100000E+02, 1.150000E+02, 1.200000E+02, 1.250000E+02 }; 771 const double MGH10_functor::y[16] = { 3.478000E+04, 2.861000E+04, 2.365000E+04, 1.963000E+04, 1.637000E+04, 1.372000E+04, 1.154000E+04, 9.744000E+03, 8.261000E+03, 7.030000E+03, 6.005000E+03, 5.147000E+03, 4.427000E+03, 3.820000E+03, 3.307000E+03, 2.872000E+03 }; 772 773 // http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml 774 void testNistMGH10(void) 775 { 776 const int n=3; 777 int info; 778 779 VectorXd x(n); 780 781 /* 782 * First try 783 */ 784 x<< 2., 400000., 25000.; 785 // do the computation 786 MGH10_functor functor; 787 LevenbergMarquardt<MGH10_functor> lm(functor); 788 info = lm.minimize(x); 789 790 // check return value 791 VERIFY_IS_EQUAL(info, 1); 792 VERIFY_IS_EQUAL(lm.nfev(), 284 ); 793 VERIFY_IS_EQUAL(lm.njev(), 249 ); 794 // check norm^2 795 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); 796 // check x 797 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 798 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 799 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 800 801 /* 802 * Second try 803 */ 804 x<< 0.02, 4000., 250.; 805 // do the computation 806 info = lm.minimize(x); 807 808 // check return value 809 VERIFY_IS_EQUAL(info, 1); 810 VERIFY_IS_EQUAL(lm.nfev(), 126); 811 VERIFY_IS_EQUAL(lm.njev(), 116); 812 // check norm^2 813 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); 814 // check x 815 VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 816 VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 817 VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 818 } 819 820 821 struct BoxBOD_functor : DenseFunctor<double> 822 { 823 BoxBOD_functor(void) : DenseFunctor<double>(2,6) {} 824 static const double x[6]; 825 int operator()(const VectorXd &b, VectorXd &fvec) 826 { 827 static const double y[6] = { 109., 149., 149., 191., 213., 224. }; 828 assert(b.size()==2); 829 assert(fvec.size()==6); 830 for(int i=0; i<6; i++) 831 fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; 832 return 0; 833 } 834 int df(const VectorXd &b, MatrixXd &fjac) 835 { 836 assert(b.size()==2); 837 assert(fjac.rows()==6); 838 assert(fjac.cols()==2); 839 for(int i=0; i<6; i++) { 840 double e = exp(-b[1]*x[i]); 841 fjac(i,0) = 1.-e; 842 fjac(i,1) = b[0]*x[i]*e; 843 } 844 return 0; 845 } 846 }; 847 const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; 848 849 // http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml 850 void testNistBoxBOD(void) 851 { 852 const int n=2; 853 int info; 854 855 VectorXd x(n); 856 857 /* 858 * First try 859 */ 860 x<< 1., 1.; 861 // do the computation 862 BoxBOD_functor functor; 863 LevenbergMarquardt<BoxBOD_functor> lm(functor); 864 lm.setFtol(1.E6*NumTraits<double>::epsilon()); 865 lm.setXtol(1.E6*NumTraits<double>::epsilon()); 866 lm.setFactor(10); 867 info = lm.minimize(x); 868 869 // check return value 870 VERIFY_IS_EQUAL(info, 1); 871 VERIFY_IS_EQUAL(lm.nfev(), 31); 872 VERIFY_IS_EQUAL(lm.njev(), 25); 873 // check norm^2 874 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); 875 // check x 876 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 877 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 878 879 /* 880 * Second try 881 */ 882 x<< 100., 0.75; 883 // do the computation 884 lm.resetParameters(); 885 lm.setFtol(NumTraits<double>::epsilon()); 886 lm.setXtol( NumTraits<double>::epsilon()); 887 info = lm.minimize(x); 888 889 // check return value 890 VERIFY_IS_EQUAL(info, 1); 891 VERIFY_IS_EQUAL(lm.nfev(), 15 ); 892 VERIFY_IS_EQUAL(lm.njev(), 14 ); 893 // check norm^2 894 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); 895 // check x 896 VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 897 VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 898 } 899 900 struct MGH17_functor : DenseFunctor<double> 901 { 902 MGH17_functor(void) : DenseFunctor<double>(5,33) {} 903 static const double x[33]; 904 static const double y[33]; 905 int operator()(const VectorXd &b, VectorXd &fvec) 906 { 907 assert(b.size()==5); 908 assert(fvec.size()==33); 909 for(int i=0; i<33; i++) 910 fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; 911 return 0; 912 } 913 int df(const VectorXd &b, MatrixXd &fjac) 914 { 915 assert(b.size()==5); 916 assert(fjac.rows()==33); 917 assert(fjac.cols()==5); 918 for(int i=0; i<33; i++) { 919 fjac(i,0) = 1.; 920 fjac(i,1) = exp(-b[3]*x[i]); 921 fjac(i,2) = exp(-b[4]*x[i]); 922 fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); 923 fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); 924 } 925 return 0; 926 } 927 }; 928 const double MGH17_functor::x[33] = { 0.000000E+00, 1.000000E+01, 2.000000E+01, 3.000000E+01, 4.000000E+01, 5.000000E+01, 6.000000E+01, 7.000000E+01, 8.000000E+01, 9.000000E+01, 1.000000E+02, 1.100000E+02, 1.200000E+02, 1.300000E+02, 1.400000E+02, 1.500000E+02, 1.600000E+02, 1.700000E+02, 1.800000E+02, 1.900000E+02, 2.000000E+02, 2.100000E+02, 2.200000E+02, 2.300000E+02, 2.400000E+02, 2.500000E+02, 2.600000E+02, 2.700000E+02, 2.800000E+02, 2.900000E+02, 3.000000E+02, 3.100000E+02, 3.200000E+02 }; 929 const double MGH17_functor::y[33] = { 8.440000E-01, 9.080000E-01, 9.320000E-01, 9.360000E-01, 9.250000E-01, 9.080000E-01, 8.810000E-01, 8.500000E-01, 8.180000E-01, 7.840000E-01, 7.510000E-01, 7.180000E-01, 6.850000E-01, 6.580000E-01, 6.280000E-01, 6.030000E-01, 5.800000E-01, 5.580000E-01, 5.380000E-01, 5.220000E-01, 5.060000E-01, 4.900000E-01, 4.780000E-01, 4.670000E-01, 4.570000E-01, 4.480000E-01, 4.380000E-01, 4.310000E-01, 4.240000E-01, 4.200000E-01, 4.140000E-01, 4.110000E-01, 4.060000E-01 }; 930 931 // http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml 932 void testNistMGH17(void) 933 { 934 const int n=5; 935 int info; 936 937 VectorXd x(n); 938 939 /* 940 * First try 941 */ 942 x<< 50., 150., -100., 1., 2.; 943 // do the computation 944 MGH17_functor functor; 945 LevenbergMarquardt<MGH17_functor> lm(functor); 946 lm.setFtol(NumTraits<double>::epsilon()); 947 lm.setXtol(NumTraits<double>::epsilon()); 948 lm.setMaxfev(1000); 949 info = lm.minimize(x); 950 951 // check return value 952 // VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success) 953 // VERIFY_IS_EQUAL(lm.nfev(), 602 ); 954 VERIFY_IS_EQUAL(lm.njev(), 545 ); 955 // check norm^2 956 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); 957 // check x 958 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 959 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 960 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 961 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 962 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 963 964 /* 965 * Second try 966 */ 967 x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; 968 // do the computation 969 lm.resetParameters(); 970 info = lm.minimize(x); 971 972 // check return value 973 VERIFY_IS_EQUAL(info, 1); 974 VERIFY_IS_EQUAL(lm.nfev(), 18); 975 VERIFY_IS_EQUAL(lm.njev(), 15); 976 // check norm^2 977 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); 978 // check x 979 VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 980 VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 981 VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 982 VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 983 VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 984 } 985 986 struct MGH09_functor : DenseFunctor<double> 987 { 988 MGH09_functor(void) : DenseFunctor<double>(4,11) {} 989 static const double _x[11]; 990 static const double y[11]; 991 int operator()(const VectorXd &b, VectorXd &fvec) 992 { 993 assert(b.size()==4); 994 assert(fvec.size()==11); 995 for(int i=0; i<11; i++) { 996 double x = _x[i], xx=x*x; 997 fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; 998 } 999 return 0; 1000 } 1001 int df(const VectorXd &b, MatrixXd &fjac) 1002 { 1003 assert(b.size()==4); 1004 assert(fjac.rows()==11); 1005 assert(fjac.cols()==4); 1006 for(int i=0; i<11; i++) { 1007 double x = _x[i], xx=x*x; 1008 double factor = 1./(xx+x*b[2]+b[3]); 1009 fjac(i,0) = (xx+x*b[1]) * factor; 1010 fjac(i,1) = b[0]*x* factor; 1011 fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; 1012 fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; 1013 } 1014 return 0; 1015 } 1016 }; 1017 const double MGH09_functor::_x[11] = { 4., 2., 1., 5.E-1 , 2.5E-01, 1.670000E-01, 1.250000E-01, 1.E-01, 8.330000E-02, 7.140000E-02, 6.250000E-02 }; 1018 const double MGH09_functor::y[11] = { 1.957000E-01, 1.947000E-01, 1.735000E-01, 1.600000E-01, 8.440000E-02, 6.270000E-02, 4.560000E-02, 3.420000E-02, 3.230000E-02, 2.350000E-02, 2.460000E-02 }; 1019 1020 // http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml 1021 void testNistMGH09(void) 1022 { 1023 const int n=4; 1024 int info; 1025 1026 VectorXd x(n); 1027 1028 /* 1029 * First try 1030 */ 1031 x<< 25., 39, 41.5, 39.; 1032 // do the computation 1033 MGH09_functor functor; 1034 LevenbergMarquardt<MGH09_functor> lm(functor); 1035 lm.setMaxfev(1000); 1036 info = lm.minimize(x); 1037 1038 // check return value 1039 VERIFY_IS_EQUAL(info, 1); 1040 VERIFY_IS_EQUAL(lm.nfev(), 490 ); 1041 VERIFY_IS_EQUAL(lm.njev(), 376 ); 1042 // check norm^2 1043 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); 1044 // check x 1045 VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 1046 VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 1047 VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 1048 VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 1049 1050 /* 1051 * Second try 1052 */ 1053 x<< 0.25, 0.39, 0.415, 0.39; 1054 // do the computation 1055 lm.resetParameters(); 1056 info = lm.minimize(x); 1057 1058 // check return value 1059 VERIFY_IS_EQUAL(info, 1); 1060 VERIFY_IS_EQUAL(lm.nfev(), 18); 1061 VERIFY_IS_EQUAL(lm.njev(), 16); 1062 // check norm^2 1063 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); 1064 // check x 1065 VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 1066 VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 1067 VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 1068 VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 1069 } 1070 1071 1072 1073 struct Bennett5_functor : DenseFunctor<double> 1074 { 1075 Bennett5_functor(void) : DenseFunctor<double>(3,154) {} 1076 static const double x[154]; 1077 static const double y[154]; 1078 int operator()(const VectorXd &b, VectorXd &fvec) 1079 { 1080 assert(b.size()==3); 1081 assert(fvec.size()==154); 1082 for(int i=0; i<154; i++) 1083 fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; 1084 return 0; 1085 } 1086 int df(const VectorXd &b, MatrixXd &fjac) 1087 { 1088 assert(b.size()==3); 1089 assert(fjac.rows()==154); 1090 assert(fjac.cols()==3); 1091 for(int i=0; i<154; i++) { 1092 double e = pow(b[1]+x[i],-1./b[2]); 1093 fjac(i,0) = e; 1094 fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); 1095 fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; 1096 } 1097 return 0; 1098 } 1099 }; 1100 const double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, 1101 11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 }; 1102 const double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 1103 ,-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 , 1104 -31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; 1105 1106 // http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml 1107 void testNistBennett5(void) 1108 { 1109 const int n=3; 1110 int info; 1111 1112 VectorXd x(n); 1113 1114 /* 1115 * First try 1116 */ 1117 x<< -2000., 50., 0.8; 1118 // do the computation 1119 Bennett5_functor functor; 1120 LevenbergMarquardt<Bennett5_functor> lm(functor); 1121 lm.setMaxfev(1000); 1122 info = lm.minimize(x); 1123 1124 // check return value 1125 VERIFY_IS_EQUAL(info, 1); 1126 VERIFY_IS_EQUAL(lm.nfev(), 758); 1127 VERIFY_IS_EQUAL(lm.njev(), 744); 1128 // check norm^2 1129 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); 1130 // check x 1131 VERIFY_IS_APPROX(x[0], -2.5235058043E+03); 1132 VERIFY_IS_APPROX(x[1], 4.6736564644E+01); 1133 VERIFY_IS_APPROX(x[2], 9.3218483193E-01); 1134 /* 1135 * Second try 1136 */ 1137 x<< -1500., 45., 0.85; 1138 // do the computation 1139 lm.resetParameters(); 1140 info = lm.minimize(x); 1141 1142 // check return value 1143 VERIFY_IS_EQUAL(info, 1); 1144 VERIFY_IS_EQUAL(lm.nfev(), 203); 1145 VERIFY_IS_EQUAL(lm.njev(), 192); 1146 // check norm^2 1147 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); 1148 // check x 1149 VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 1150 VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); 1151 VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); 1152 } 1153 1154 struct thurber_functor : DenseFunctor<double> 1155 { 1156 thurber_functor(void) : DenseFunctor<double>(7,37) {} 1157 static const double _x[37]; 1158 static const double _y[37]; 1159 int operator()(const VectorXd &b, VectorXd &fvec) 1160 { 1161 // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 1162 assert(b.size()==7); 1163 assert(fvec.size()==37); 1164 for(int i=0; i<37; i++) { 1165 double x=_x[i], xx=x*x, xxx=xx*x; 1166 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]; 1167 } 1168 return 0; 1169 } 1170 int df(const VectorXd &b, MatrixXd &fjac) 1171 { 1172 assert(b.size()==7); 1173 assert(fjac.rows()==37); 1174 assert(fjac.cols()==7); 1175 for(int i=0; i<37; i++) { 1176 double x=_x[i], xx=x*x, xxx=xx*x; 1177 double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 1178 fjac(i,0) = 1.*fact; 1179 fjac(i,1) = x*fact; 1180 fjac(i,2) = xx*fact; 1181 fjac(i,3) = xxx*fact; 1182 fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 1183 fjac(i,4) = x*fact; 1184 fjac(i,5) = xx*fact; 1185 fjac(i,6) = xxx*fact; 1186 } 1187 return 0; 1188 } 1189 }; 1190 const double thurber_functor::_x[37] = { -3.067E0, -2.981E0, -2.921E0, -2.912E0, -2.840E0, -2.797E0, -2.702E0, -2.699E0, -2.633E0, -2.481E0, -2.363E0, -2.322E0, -1.501E0, -1.460E0, -1.274E0, -1.212E0, -1.100E0, -1.046E0, -0.915E0, -0.714E0, -0.566E0, -0.545E0, -0.400E0, -0.309E0, -0.109E0, -0.103E0, 0.010E0, 0.119E0, 0.377E0, 0.790E0, 0.963E0, 1.006E0, 1.115E0, 1.572E0, 1.841E0, 2.047E0, 2.200E0 }; 1191 const double thurber_functor::_y[37] = { 80.574E0, 84.248E0, 87.264E0, 87.195E0, 89.076E0, 89.608E0, 89.868E0, 90.101E0, 92.405E0, 95.854E0, 100.696E0, 101.060E0, 401.672E0, 390.724E0, 567.534E0, 635.316E0, 733.054E0, 759.087E0, 894.206E0, 990.785E0, 1090.109E0, 1080.914E0, 1122.643E0, 1178.351E0, 1260.531E0, 1273.514E0, 1288.339E0, 1327.543E0, 1353.863E0, 1414.509E0, 1425.208E0, 1421.384E0, 1442.962E0, 1464.350E0, 1468.705E0, 1447.894E0, 1457.628E0}; 1192 1193 // http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml 1194 void testNistThurber(void) 1195 { 1196 const int n=7; 1197 int info; 1198 1199 VectorXd x(n); 1200 1201 /* 1202 * First try 1203 */ 1204 x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; 1205 // do the computation 1206 thurber_functor functor; 1207 LevenbergMarquardt<thurber_functor> lm(functor); 1208 lm.setFtol(1.E4*NumTraits<double>::epsilon()); 1209 lm.setXtol(1.E4*NumTraits<double>::epsilon()); 1210 info = lm.minimize(x); 1211 1212 // check return value 1213 VERIFY_IS_EQUAL(info, 1); 1214 VERIFY_IS_EQUAL(lm.nfev(), 39); 1215 VERIFY_IS_EQUAL(lm.njev(), 36); 1216 // check norm^2 1217 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); 1218 // check x 1219 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1220 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1221 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1222 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1223 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1224 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1225 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1226 1227 /* 1228 * Second try 1229 */ 1230 x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; 1231 // do the computation 1232 lm.resetParameters(); 1233 lm.setFtol(1.E4*NumTraits<double>::epsilon()); 1234 lm.setXtol(1.E4*NumTraits<double>::epsilon()); 1235 info = lm.minimize(x); 1236 1237 // check return value 1238 VERIFY_IS_EQUAL(info, 1); 1239 VERIFY_IS_EQUAL(lm.nfev(), 29); 1240 VERIFY_IS_EQUAL(lm.njev(), 28); 1241 // check norm^2 1242 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); 1243 // check x 1244 VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1245 VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1246 VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1247 VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1248 VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1249 VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1250 VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1251 } 1252 1253 struct rat43_functor : DenseFunctor<double> 1254 { 1255 rat43_functor(void) : DenseFunctor<double>(4,15) {} 1256 static const double x[15]; 1257 static const double y[15]; 1258 int operator()(const VectorXd &b, VectorXd &fvec) 1259 { 1260 assert(b.size()==4); 1261 assert(fvec.size()==15); 1262 for(int i=0; i<15; i++) 1263 fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; 1264 return 0; 1265 } 1266 int df(const VectorXd &b, MatrixXd &fjac) 1267 { 1268 assert(b.size()==4); 1269 assert(fjac.rows()==15); 1270 assert(fjac.cols()==4); 1271 for(int i=0; i<15; i++) { 1272 double e = exp(b[1]-b[2]*x[i]); 1273 double power = -1./b[3]; 1274 fjac(i,0) = pow(1.+e, power); 1275 fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); 1276 fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); 1277 fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); 1278 } 1279 return 0; 1280 } 1281 }; 1282 const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; 1283 const double rat43_functor::y[15] = { 16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41 }; 1284 1285 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml 1286 void testNistRat43(void) 1287 { 1288 const int n=4; 1289 int info; 1290 1291 VectorXd x(n); 1292 1293 /* 1294 * First try 1295 */ 1296 x<< 100., 10., 1., 1.; 1297 // do the computation 1298 rat43_functor functor; 1299 LevenbergMarquardt<rat43_functor> lm(functor); 1300 lm.setFtol(1.E6*NumTraits<double>::epsilon()); 1301 lm.setXtol(1.E6*NumTraits<double>::epsilon()); 1302 info = lm.minimize(x); 1303 1304 // check return value 1305 VERIFY_IS_EQUAL(info, 1); 1306 VERIFY_IS_EQUAL(lm.nfev(), 27); 1307 VERIFY_IS_EQUAL(lm.njev(), 20); 1308 // check norm^2 1309 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); 1310 // check x 1311 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1312 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1313 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1314 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1315 1316 /* 1317 * Second try 1318 */ 1319 x<< 700., 5., 0.75, 1.3; 1320 // do the computation 1321 lm.resetParameters(); 1322 lm.setFtol(1.E5*NumTraits<double>::epsilon()); 1323 lm.setXtol(1.E5*NumTraits<double>::epsilon()); 1324 info = lm.minimize(x); 1325 1326 // check return value 1327 VERIFY_IS_EQUAL(info, 1); 1328 VERIFY_IS_EQUAL(lm.nfev(), 9); 1329 VERIFY_IS_EQUAL(lm.njev(), 8); 1330 // check norm^2 1331 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); 1332 // check x 1333 VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1334 VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1335 VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1336 VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1337 } 1338 1339 1340 1341 struct eckerle4_functor : DenseFunctor<double> 1342 { 1343 eckerle4_functor(void) : DenseFunctor<double>(3,35) {} 1344 static const double x[35]; 1345 static const double y[35]; 1346 int operator()(const VectorXd &b, VectorXd &fvec) 1347 { 1348 assert(b.size()==3); 1349 assert(fvec.size()==35); 1350 for(int i=0; i<35; i++) 1351 fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; 1352 return 0; 1353 } 1354 int df(const VectorXd &b, MatrixXd &fjac) 1355 { 1356 assert(b.size()==3); 1357 assert(fjac.rows()==35); 1358 assert(fjac.cols()==3); 1359 for(int i=0; i<35; i++) { 1360 double b12 = b[1]*b[1]; 1361 double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); 1362 fjac(i,0) = e / b[1]; 1363 fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; 1364 fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; 1365 } 1366 return 0; 1367 } 1368 }; 1369 const double eckerle4_functor::x[35] = { 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0}; 1370 const double eckerle4_functor::y[35] = { 0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 }; 1371 1372 // http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml 1373 void testNistEckerle4(void) 1374 { 1375 const int n=3; 1376 int info; 1377 1378 VectorXd x(n); 1379 1380 /* 1381 * First try 1382 */ 1383 x<< 1., 10., 500.; 1384 // do the computation 1385 eckerle4_functor functor; 1386 LevenbergMarquardt<eckerle4_functor> lm(functor); 1387 info = lm.minimize(x); 1388 1389 // check return value 1390 VERIFY_IS_EQUAL(info, 1); 1391 VERIFY_IS_EQUAL(lm.nfev(), 18); 1392 VERIFY_IS_EQUAL(lm.njev(), 15); 1393 // check norm^2 1394 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); 1395 // check x 1396 VERIFY_IS_APPROX(x[0], 1.5543827178); 1397 VERIFY_IS_APPROX(x[1], 4.0888321754); 1398 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1399 1400 /* 1401 * Second try 1402 */ 1403 x<< 1.5, 5., 450.; 1404 // do the computation 1405 info = lm.minimize(x); 1406 1407 // check return value 1408 VERIFY_IS_EQUAL(info, 1); 1409 VERIFY_IS_EQUAL(lm.nfev(), 7); 1410 VERIFY_IS_EQUAL(lm.njev(), 6); 1411 // check norm^2 1412 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); 1413 // check x 1414 VERIFY_IS_APPROX(x[0], 1.5543827178); 1415 VERIFY_IS_APPROX(x[1], 4.0888321754); 1416 VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1417 } 1418 1419 void test_levenberg_marquardt() 1420 { 1421 // Tests using the examples provided by (c)minpack 1422 CALL_SUBTEST(testLmder1()); 1423 CALL_SUBTEST(testLmder()); 1424 CALL_SUBTEST(testLmdif1()); 1425 // CALL_SUBTEST(testLmstr1()); 1426 // CALL_SUBTEST(testLmstr()); 1427 CALL_SUBTEST(testLmdif()); 1428 1429 // NIST tests, level of difficulty = "Lower" 1430 CALL_SUBTEST(testNistMisra1a()); 1431 CALL_SUBTEST(testNistChwirut2()); 1432 1433 // NIST tests, level of difficulty = "Average" 1434 CALL_SUBTEST(testNistHahn1()); 1435 CALL_SUBTEST(testNistMisra1d()); 1436 CALL_SUBTEST(testNistMGH17()); 1437 CALL_SUBTEST(testNistLanczos1()); 1438 1439 // // NIST tests, level of difficulty = "Higher" 1440 CALL_SUBTEST(testNistRat42()); 1441 CALL_SUBTEST(testNistMGH10()); 1442 CALL_SUBTEST(testNistBoxBOD()); 1443 // CALL_SUBTEST(testNistMGH09()); 1444 CALL_SUBTEST(testNistBennett5()); 1445 CALL_SUBTEST(testNistThurber()); 1446 CALL_SUBTEST(testNistRat43()); 1447 CALL_SUBTEST(testNistEckerle4()); 1448 } 1449
