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