NonLinearOptimization.cpp
Go to the documentation of this file.
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 
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 
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 
130  Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
131  Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
132 
133  int inputs() const { return m_inputs; }
134  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 {
142  lmder_functor(void): Functor<double>(3,15) {}
143  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 
159  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 
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 
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 {
252  hybrj_functor(void) : Functor<double>(9,9) {}
253 
254  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  }
270  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 
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 
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 {
356  hybrd_functor(void) : Functor<double>(9,9) {}
357  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 
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 
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 {
438  lmstr_functor(void) : Functor<double>(3,15) {}
439  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  }
458  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 
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 
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 {
536  lmdif_functor(void) : Functor<double>(3,15) {}
537  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 
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 
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 {
639  chwirut2_functor(void) : Functor<double>(3,54) {}
640  static const double m_x[54];
641  static const double m_y[54];
642  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  }
654  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
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 {
724  misra1a_functor(void) : Functor<double>(2,14) {}
725  static const double m_x[14];
726  static const double m_y[14];
727  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  }
736  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
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 {
796  hahn1_functor(void) : Functor<double>(7,236) {}
797  static const double m_x[236];
798  int operator()(const VectorXd &b, VectorXd &fvec)
799  {
800  static const double m_y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0 , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0 , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0 , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 ,
801  16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0 , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0 , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 12.596E0 ,
802 13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0 , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0 , 20.935E0 , 21.035E0 , 20.93E0 , 21.074E0 , 21.085E0 , 20.935E0 };
803 
804  // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
805 
806  assert(b.size()==7);
807  assert(fvec.size()==236);
808  for(int i=0; i<236; i++) {
809  double x=m_x[i], xx=x*x, xxx=xx*x;
810  fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - m_y[i];
811  }
812  return 0;
813  }
814 
815  int df(const VectorXd &b, MatrixXd &fjac)
816  {
817  assert(b.size()==7);
818  assert(fjac.rows()==236);
819  assert(fjac.cols()==7);
820  for(int i=0; i<236; i++) {
821  double x=m_x[i], xx=x*x, xxx=xx*x;
822  double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
823  fjac(i,0) = 1.*fact;
824  fjac(i,1) = x*fact;
825  fjac(i,2) = xx*fact;
826  fjac(i,3) = xxx*fact;
827  fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
828  fjac(i,4) = x*fact;
829  fjac(i,5) = xx*fact;
830  fjac(i,6) = xxx*fact;
831  }
832  return 0;
833  }
834 };
835 const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 ,
836 282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 ,
837 141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0 , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0 , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0};
838 
839 // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml
840 void testNistHahn1(void)
841 {
842  const int n=7;
843  int info;
844 
845  VectorXd x(n);
846 
847  /*
848  * First try
849  */
850  x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
851  // do the computation
852  hahn1_functor functor;
853  LevenbergMarquardt<hahn1_functor> lm(functor);
854  info = lm.minimize(x);
855 
856  // check return value
857  VERIFY_IS_EQUAL(info, 1);
858  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 {
895  misra1d_functor(void) : Functor<double>(2,14) {}
896  static const double x[14];
897  static const double y[14];
898  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  }
907  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
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 {
969  lanczos1_functor(void) : Functor<double>(6,24) {}
970  static const double x[24];
971  static const double y[24];
972  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  }
980  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
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 {
1055  rat42_functor(void) : Functor<double>(3,9) {}
1056  static const double x[9];
1057  static const double y[9];
1058  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 
1068  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
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 {
1132  MGH10_functor(void) : Functor<double>(3,16) {}
1133  static const double x[16];
1134  static const double y[16];
1135  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  }
1143  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
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 {
1209  BoxBOD_functor(void) : Functor<double>(2,6) {}
1210  static const double x[6];
1211  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  }
1220  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
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 {
1286  MGH17_functor(void) : Functor<double>(5,33) {}
1287  static const double x[33];
1288  static const double y[33];
1289  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  }
1297  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
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 {
1371  MGH09_functor(void) : Functor<double>(4,11) {}
1372  static const double _x[11];
1373  static const double y[11];
1374  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  }
1384  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
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 {
1456  Bennett5_functor(void) : Functor<double>(3,154) {}
1457  static const double x[154];
1458  static const double y[154];
1459  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  }
1467  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
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 {
1535  thurber_functor(void) : Functor<double>(7,37) {}
1536  static const double _x[37];
1537  static const double _y[37];
1538  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  }
1549  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
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 {
1632  rat43_functor(void) : Functor<double>(4,15) {}
1633  static const double x[15];
1634  static const double y[15];
1635  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  }
1643  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
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 {
1718  eckerle4_functor(void) : Functor<double>(3,35) {}
1719  static const double x[35];
1720  static const double y[35];
1721  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  }
1729  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
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 
1792 EIGEN_DECLARE_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 
hybrd_functor::operator()
int operator()(const VectorXd &x, VectorXd &fvec) const
Definition: NonLinearOptimization.cpp:357
Eigen::LevenbergMarquardt::fjac
JacobianType fjac
Definition: NonLinearOptimization/LevenbergMarquardt.h:108
DisableStupidWarnings.h
rat42_functor::y
static const double y[9]
Definition: levenberg_marquardt.cpp:673
testNistBoxBOD
void testNistBoxBOD(void)
Definition: NonLinearOptimization.cpp:1236
misra1a_functor::m_y
static const double m_y[14]
Definition: levenberg_marquardt.cpp:335
MGH17_functor::MGH17_functor
MGH17_functor(void)
Definition: NonLinearOptimization.cpp:1286
testHybrd
void testHybrd()
Definition: NonLinearOptimization.cpp:402
rat42_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1068
Functor::JacobianType
Matrix< Scalar, ValuesAtCompileTime, InputsAtCompileTime > JacobianType
Definition: NonLinearOptimization.cpp:126
eckerle4_functor::eckerle4_functor
eckerle4_functor(void)
Definition: NonLinearOptimization.cpp:1718
hybrj_functor::df
int df(const VectorXd &x, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:270
rat43_functor
Definition: levenberg_marquardt.cpp:1282
testNistMGH17
void testNistMGH17(void)
Definition: NonLinearOptimization.cpp:1316
misra1d_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:898
Eigen::LevenbergMarquardt::parameters
Parameters parameters
Definition: NonLinearOptimization/LevenbergMarquardt.h:106
rat43_functor::rat43_functor
rat43_functor(void)
Definition: NonLinearOptimization.cpp:1632
testNistRat43
void testNistRat43(void)
Definition: NonLinearOptimization.cpp:1663
thurber_functor::thurber_functor
thurber_functor(void)
Definition: NonLinearOptimization.cpp:1535
testNistMGH09
void testNistMGH09(void)
Definition: NonLinearOptimization.cpp:1404
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
lmder_functor::operator()
int operator()(const VectorXd &x, VectorXd &fvec) const
Definition: NonLinearOptimization.cpp:143
Eigen::PermutationMatrix::indices
const IndicesType & indices() const
Definition: PermutationMatrix.h:360
lmder_functor::lmder_functor
lmder_functor(void)
Definition: NonLinearOptimization.cpp:142
hahn1_functor::hahn1_functor
hahn1_functor(void)
Definition: NonLinearOptimization.cpp:796
rat42_functor::x
static const double x[9]
Definition: levenberg_marquardt.cpp:672
testLmdif1
void testLmdif1()
Definition: NonLinearOptimization.cpp:559
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
MGH17_functor::x
static const double x[33]
Definition: levenberg_marquardt.cpp:931
hybrj_functor::operator()
int operator()(const VectorXd &x, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:254
b
Scalar * b
Definition: benchVecAdd.cpp:17
eckerle4_functor::x
static const double x[35]
Definition: levenberg_marquardt.cpp:1373
testChkder
void testChkder()
Definition: NonLinearOptimization.cpp:70
thurber_functor::_y
static const double _y[37]
Definition: levenberg_marquardt.cpp:1187
chwirut2_functor::chwirut2_functor
chwirut2_functor(void)
Definition: NonLinearOptimization.cpp:639
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
BoxBOD_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1211
thurber_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1538
rat43_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1635
misra1a_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:736
MGH17_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1297
MGH10_functor
Definition: levenberg_marquardt.cpp:748
misra1a_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:727
MGH17_functor::y
static const double y[33]
Definition: levenberg_marquardt.cpp:932
lmstr_functor::operator()
int operator()(const VectorXd &x, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:439
hybrd_functor::hybrd_functor
hybrd_functor(void)
Definition: NonLinearOptimization.cpp:356
Bennett5_functor::x
static const double x[154]
Definition: levenberg_marquardt.cpp:1105
testNistMGH10
void testNistMGH10(void)
Definition: NonLinearOptimization.cpp:1162
MGH10_functor::MGH10_functor
MGH10_functor(void)
Definition: NonLinearOptimization.cpp:1132
rat42_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1058
Bennett5_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1467
lanczos1_functor::x
static const double x[24]
Definition: levenberg_marquardt.cpp:586
testNistChwirut2
void testNistChwirut2(void)
Definition: NonLinearOptimization.cpp:674
MGH10_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1135
MGH09_functor::MGH09_functor
MGH09_functor(void)
Definition: NonLinearOptimization.cpp:1371
Eigen::LevenbergMarquardt::lmstr1
LevenbergMarquardtSpace::Status lmstr1(FVectorType &x, const Scalar tol=sqrt_epsilon())
Definition: NonLinearOptimization/LevenbergMarquardt.h:363
chwirut2_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:642
rat42_functor
Definition: levenberg_marquardt.cpp:669
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
thurber_functor::_x
static const double _x[37]
Definition: levenberg_marquardt.cpp:1186
misra1a_functor::misra1a_functor
misra1a_functor(void)
Definition: NonLinearOptimization.cpp:724
Eigen::LevenbergMarquardt::fvec
FVectorType fvec
Definition: NonLinearOptimization/LevenbergMarquardt.h:107
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
Functor::InputsAtCompileTime
@ InputsAtCompileTime
Definition: NonLinearOptimization.cpp:121
Eigen::LevenbergMarquardt
Performs non linear optimization over a non-linear function, using a variant of the Levenberg Marquar...
Definition: LevenbergMarquardt/LevenbergMarquardt.h:110
Functor::Functor
Functor()
Definition: NonLinearOptimization.cpp:130
solver
BiCGSTAB< SparseMatrix< double > > solver
Definition: BiCGSTAB_simple.cpp:5
testNistMisra1a
void testNistMisra1a(void)
Definition: NonLinearOptimization.cpp:752
testLmdif
void testLmdif()
Definition: NonLinearOptimization.cpp:589
MGH09_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1374
chwirut2_functor::m_y
static const double m_y[54]
Definition: levenberg_marquardt.cpp:248
n
int n
Definition: BiCGSTAB_simple.cpp:1
epsilon
static double epsilon
Definition: testRot3.cpp:37
hahn1_functor::m_x
static const double m_x[236]
Definition: levenberg_marquardt.cpp:408
Functor::Scalar
_Scalar Scalar
Definition: NonLinearOptimization.cpp:119
Functor
Definition: NonLinearOptimization.cpp:117
Eigen::LevenbergMarquardt::resetParameters
void resetParameters()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:145
MGH10_functor::y
static const double y[16]
Definition: levenberg_marquardt.cpp:752
hybrj_functor::hybrj_functor
hybrj_functor(void)
Definition: NonLinearOptimization.cpp:252
lmstr_functor::df
int df(const VectorXd &x, VectorXd &jac_row, VectorXd::Index rownb)
Definition: NonLinearOptimization.cpp:458
hybrj_functor
Definition: NonLinearOptimization.cpp:250
testNistBennett5
void testNistBennett5(void)
Definition: NonLinearOptimization.cpp:1488
testLmstr
void testLmstr()
Definition: NonLinearOptimization.cpp:504
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
misra1d_functor::y
static const double y[14]
Definition: levenberg_marquardt.cpp:511
rat43_functor::y
static const double y[15]
Definition: levenberg_marquardt.cpp:1286
Bennett5_functor
Definition: levenberg_marquardt.cpp:1102
hybrd_functor
Definition: NonLinearOptimization.cpp:354
lmdif_functor::operator()
int operator()(const VectorXd &x, VectorXd &fvec) const
Definition: NonLinearOptimization.cpp:537
misra1d_functor::x
static const double x[14]
Definition: levenberg_marquardt.cpp:510
Bennett5_functor::y
static const double y[154]
Definition: levenberg_marquardt.cpp:1106
Eigen::NumericalDiff
Definition: NumericalDiff.h:36
testNistRat42
void testNistRat42(void)
Definition: NonLinearOptimization.cpp:1086
Eigen::LevenbergMarquardt::lmder1
LevenbergMarquardtSpace::Status lmder1(FVectorType &x, const Scalar tol=std::sqrt(NumTraits< Scalar >::epsilon()))
Definition: LevenbergMarquardt/LevenbergMarquardt.h:344
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
Functor::inputs
int inputs() const
Definition: NonLinearOptimization.cpp:133
testLmder
void testLmder()
Definition: NonLinearOptimization.cpp:203
BoxBOD_functor::x
static const double x[6]
Definition: levenberg_marquardt.cpp:847
misra1d_functor::misra1d_functor
misra1d_functor(void)
Definition: NonLinearOptimization.cpp:895
MGH10_functor::x
static const double x[16]
Definition: levenberg_marquardt.cpp:751
testHybrj
void testHybrj()
Definition: NonLinearOptimization.cpp:319
testNistThurber
void testNistThurber(void)
Definition: NonLinearOptimization.cpp:1573
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
misra1a_functor::m_x
static const double m_x[14]
Definition: levenberg_marquardt.cpp:334
ceres::pow
Jet< T, N > pow(const Jet< T, N > &f, double g)
Definition: jet.h:570
testNistHahn1
void testNistHahn1(void)
Definition: NonLinearOptimization.cpp:840
eckerle4_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1729
Functor::values
int values() const
Definition: NonLinearOptimization.cpp:134
MGH09_functor
Definition: levenberg_marquardt.cpp:1015
lmstr_functor
Definition: NonLinearOptimization.cpp:436
hahn1_functor
Definition: levenberg_marquardt.cpp:405
hahn1_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:815
Functor::InputType
Matrix< Scalar, InputsAtCompileTime, 1 > InputType
Definition: NonLinearOptimization.cpp:124
y
Scalar * y
Definition: level1_cplx_impl.h:124
rat43_functor::x
static const double x[15]
Definition: levenberg_marquardt.cpp:1285
lanczos1_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:972
BoxBOD_functor::BoxBOD_functor
BoxBOD_functor(void)
Definition: NonLinearOptimization.cpp:1209
Eigen::LevenbergMarquardt::permutation
PermutationMatrix< Dynamic, Dynamic > permutation
Definition: NonLinearOptimization/LevenbergMarquardt.h:109
E
DiscreteKey E(5, 2)
Eigen::DenseFunctor< double >::values
int values() const
Definition: LevenbergMarquardt/LevenbergMarquardt.h:59
Functor::Functor
Functor(int inputs, int values)
Definition: NonLinearOptimization.cpp:131
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
Bennett5_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1459
testHybrj1
void testHybrj1()
Definition: NonLinearOptimization.cpp:288
thurber_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1549
MGH10_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1143
lmder_functor
Definition: levenberg_marquardt.cpp:29
BoxBOD_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1220
lanczos1_functor
Definition: levenberg_marquardt.cpp:583
testHybrd1
void testHybrd1()
Definition: NonLinearOptimization.cpp:376
misra1d_functor
Definition: levenberg_marquardt.cpp:507
Eigen::internal::covar
void covar(Matrix< Scalar, Dynamic, Dynamic > &r, const VectorXi &ipvt, Scalar tol=std::sqrt(NumTraits< Scalar >::epsilon()))
Definition: LMcovar.h:20
MGH17_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1289
Functor::ValueType
Matrix< Scalar, ValuesAtCompileTime, 1 > ValueType
Definition: NonLinearOptimization.cpp:125
Eigen::LevenbergMarquardt::minimizeOptimumStorage
LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x)
Definition: NonLinearOptimization/LevenbergMarquardt.h:613
main.h
Eigen::DenseIndex
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: Meta.h:66
Functor::m_inputs
const int m_inputs
Definition: NonLinearOptimization.cpp:128
MGH09_functor::y
static const double y[11]
Definition: levenberg_marquardt.cpp:1019
testLmder1
void testLmder1()
Definition: NonLinearOptimization.cpp:176
eckerle4_functor::y
static const double y[35]
Definition: levenberg_marquardt.cpp:1374
lanczos1_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:980
BoxBOD_functor
Definition: levenberg_marquardt.cpp:844
lmdif_functor
Definition: levenberg_marquardt.cpp:141
MGH09_functor::_x
static const double _x[11]
Definition: levenberg_marquardt.cpp:1018
misra1d_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:907
Eigen::internal::chkder
void chkder(const Matrix< Scalar, Dynamic, 1 > &x, const Matrix< Scalar, Dynamic, 1 > &fvec, const Matrix< Scalar, Dynamic, Dynamic > &fjac, Matrix< Scalar, Dynamic, 1 > &xp, const Matrix< Scalar, Dynamic, 1 > &fvecp, int mode, Matrix< Scalar, Dynamic, 1 > &err)
Definition: chkder.h:9
rat43_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1643
eckerle4_functor
Definition: levenberg_marquardt.cpp:1370
testNistEckerle4
void testNistEckerle4(void)
Definition: NonLinearOptimization.cpp:1748
MGH09_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:1384
hahn1_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:798
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Eigen::HybridNonLinearSolver
Finds a zero of a system of n nonlinear functions in n variables by a modification of the Powell hybr...
Definition: HybridNonLinearSolver.h:43
lmstr_functor::lmstr_functor
lmstr_functor(void)
Definition: NonLinearOptimization.cpp:438
lanczos1_functor::y
static const double y[24]
Definition: levenberg_marquardt.cpp:587
chwirut2_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: NonLinearOptimization.cpp:654
chwirut2_functor::m_x
static const double m_x[54]
Definition: levenberg_marquardt.cpp:247
Bennett5_functor::Bennett5_functor
Bennett5_functor(void)
Definition: NonLinearOptimization.cpp:1456
fcn_chkder
int fcn_chkder(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag)
Definition: NonLinearOptimization.cpp:27
chwirut2_functor
Definition: levenberg_marquardt.cpp:244
Eigen::LevenbergMarquardt::minimize
LevenbergMarquardtSpace::Status minimize(FVectorType &x)
Definition: LevenbergMarquardt/LevenbergMarquardt.h:277
testNistLanczos1
void testNistLanczos1(void)
Definition: NonLinearOptimization.cpp:1000
rat42_functor::rat42_functor
rat42_functor(void)
Definition: NonLinearOptimization.cpp:1055
Functor::m_values
const int m_values
Definition: NonLinearOptimization.cpp:128
lmdif_functor::lmdif_functor
lmdif_functor(void)
Definition: NonLinearOptimization.cpp:536
Functor::ValuesAtCompileTime
@ ValuesAtCompileTime
Definition: NonLinearOptimization.cpp:122
misra1a_functor
Definition: levenberg_marquardt.cpp:331
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(NonLinearOptimization)
Definition: NonLinearOptimization.cpp:1792
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
testNistMisra1d
void testNistMisra1d(void)
Definition: NonLinearOptimization.cpp:924
MGH17_functor
Definition: levenberg_marquardt.cpp:928
testLmstr1
void testLmstr1()
Definition: NonLinearOptimization.cpp:476
eckerle4_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: NonLinearOptimization.cpp:1721
CALL_SUBTEST
#define CALL_SUBTEST(FUNC)
Definition: main.h:399
Eigen::LevenbergMarquardt::nfev
Index nfev
Definition: NonLinearOptimization/LevenbergMarquardt.h:110
thurber_functor
Definition: levenberg_marquardt.cpp:1183
VERIFY
#define VERIFY(a)
Definition: main.h:380
lmder_functor::df
int df(const VectorXd &x, MatrixXd &fjac) const
Definition: NonLinearOptimization.cpp:159
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
LM_CHECK_N_ITERS
#define LM_CHECK_N_ITERS(SOLVER, NFEV, NJEV)
Definition: NonLinearOptimization.cpp:18
lanczos1_functor::lanczos1_functor
lanczos1_functor(void)
Definition: NonLinearOptimization.cpp:969

gtsam
Author(s):
autogenerated on Tue Jun 25 2024 03:01:45