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


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:04