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