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