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