levenberg_marquardt.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 // Copyright (C) 2012 desire Nuentsa <desire.nuentsa_wakam@inria.fr
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 
12 // FIXME: These tests all check for hard-coded values. Ideally, parameters and start estimates should be randomized.
13 
14 
15 #include <stdio.h>
16 
17 #include "main.h"
18 #include <unsupported/Eigen/LevenbergMarquardt>
19 
20 // This disables some useless Warnings on MSVC.
21 // It is intended to be done for this test only.
22 #include <Eigen/src/Core/util/DisableStupidWarnings.h>
23 
24 using std::sqrt;
25 
26 // tolerance for chekcing number of iterations
27 #define LM_EVAL_COUNT_TOL 4/3
28 
29 struct lmder_functor : DenseFunctor<double>
30 {
31  lmder_functor(void): DenseFunctor<double>(3,15) {}
32  int operator()(const VectorXd &x, VectorXd &fvec) const
33  {
34  double tmp1, tmp2, tmp3;
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  for (int i = 0; i < values(); i++)
39  {
40  tmp1 = i+1;
41  tmp2 = 16 - i - 1;
42  tmp3 = (i>=8)? tmp2 : tmp1;
43  fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
44  }
45  return 0;
46  }
47 
48  int df(const VectorXd &x, MatrixXd &fjac) const
49  {
50  double tmp1, tmp2, tmp3, tmp4;
51  for (int i = 0; i < values(); i++)
52  {
53  tmp1 = i+1;
54  tmp2 = 16 - i - 1;
55  tmp3 = (i>=8)? tmp2 : tmp1;
56  tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
57  fjac(i,0) = -1;
58  fjac(i,1) = tmp1*tmp2/tmp4;
59  fjac(i,2) = tmp1*tmp3/tmp4;
60  }
61  return 0;
62  }
63 };
64 
65 void testLmder1()
66 {
67  int n=3, info;
68 
69  VectorXd x;
70 
71  /* the following starting values provide a rough fit. */
72  x.setConstant(n, 1.);
73 
74  // do the computation
75  lmder_functor functor;
76  LevenbergMarquardt<lmder_functor> lm(functor);
77  info = lm.lmder1(x);
78 
79  // check return value
80  VERIFY_IS_EQUAL(info, 1);
81  VERIFY_IS_EQUAL(lm.nfev(), 6);
82  VERIFY_IS_EQUAL(lm.njev(), 5);
83 
84  // check norm
85  VERIFY_IS_APPROX(lm.fvec().blueNorm(), 0.09063596);
86 
87  // check x
88  VectorXd x_ref(n);
89  x_ref << 0.08241058, 1.133037, 2.343695;
90  VERIFY_IS_APPROX(x, x_ref);
91 }
92 
93 void testLmder()
94 {
95  const int m=15, n=3;
96  int info;
97  double fnorm, covfac;
98  VectorXd x;
99 
100  /* the following starting values provide a rough fit. */
101  x.setConstant(n, 1.);
102 
103  // do the computation
104  lmder_functor functor;
105  LevenbergMarquardt<lmder_functor> lm(functor);
106  info = lm.minimize(x);
107 
108  // check return values
109  VERIFY_IS_EQUAL(info, 1);
110  VERIFY_IS_EQUAL(lm.nfev(), 6);
111  VERIFY_IS_EQUAL(lm.njev(), 5);
112 
113  // check norm
114  fnorm = lm.fvec().blueNorm();
115  VERIFY_IS_APPROX(fnorm, 0.09063596);
116 
117  // check x
118  VectorXd x_ref(n);
119  x_ref << 0.08241058, 1.133037, 2.343695;
120  VERIFY_IS_APPROX(x, x_ref);
121 
122  // check covariance
123  covfac = fnorm*fnorm/(m-n);
124  internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm
125 
126  MatrixXd cov_ref(n,n);
127  cov_ref <<
128  0.0001531202, 0.002869941, -0.002656662,
129  0.002869941, 0.09480935, -0.09098995,
130  -0.002656662, -0.09098995, 0.08778727;
131 
132 // std::cout << fjac*covfac << std::endl;
133 
134  MatrixXd cov;
135  cov = covfac*lm.matrixR().topLeftCorner<n,n>();
136  VERIFY_IS_APPROX( cov, cov_ref);
137  // TODO: why isn't this allowed ? :
138  // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
139 }
140 
141 struct lmdif_functor : DenseFunctor<double>
142 {
143  lmdif_functor(void) : DenseFunctor<double>(3,15) {}
144  int operator()(const VectorXd &x, VectorXd &fvec) const
145  {
146  int i;
147  double tmp1,tmp2,tmp3;
148  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,
149  3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0};
150 
151  assert(x.size()==3);
152  assert(fvec.size()==15);
153  for (i=0; i<15; i++)
154  {
155  tmp1 = i+1;
156  tmp2 = 15 - i;
157  tmp3 = tmp1;
158 
159  if (i >= 8) tmp3 = tmp2;
160  fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
161  }
162  return 0;
163  }
164 };
165 
166 void testLmdif1()
167 {
168  const int n=3;
169  int info;
170 
171  VectorXd x(n), fvec(15);
172 
173  /* the following starting values provide a rough fit. */
174  x.setConstant(n, 1.);
175 
176  // do the computation
177  lmdif_functor functor;
178  DenseIndex nfev;
179  info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev);
180 
181  // check return value
182  VERIFY_IS_EQUAL(info, 1);
183 // VERIFY_IS_EQUAL(nfev, 26);
184 
185  // check norm
186  functor(x, fvec);
187  VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596);
188 
189  // check x
190  VectorXd x_ref(n);
191  x_ref << 0.0824106, 1.1330366, 2.3436947;
192  VERIFY_IS_APPROX(x, x_ref);
193 
194 }
195 
196 void testLmdif()
197 {
198  const int m=15, n=3;
199  int info;
200  double fnorm, covfac;
201  VectorXd x(n);
202 
203  /* the following starting values provide a rough fit. */
204  x.setConstant(n, 1.);
205 
206  // do the computation
207  lmdif_functor functor;
208  NumericalDiff<lmdif_functor> numDiff(functor);
209  LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff);
210  info = lm.minimize(x);
211 
212  // check return values
213  VERIFY_IS_EQUAL(info, 1);
214 // VERIFY_IS_EQUAL(lm.nfev(), 26);
215 
216  // check norm
217  fnorm = lm.fvec().blueNorm();
218  VERIFY_IS_APPROX(fnorm, 0.09063596);
219 
220  // check x
221  VectorXd x_ref(n);
222  x_ref << 0.08241058, 1.133037, 2.343695;
223  VERIFY_IS_APPROX(x, x_ref);
224 
225  // check covariance
226  covfac = fnorm*fnorm/(m-n);
227  internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm
228 
229  MatrixXd cov_ref(n,n);
230  cov_ref <<
231  0.0001531202, 0.002869942, -0.002656662,
232  0.002869942, 0.09480937, -0.09098997,
233  -0.002656662, -0.09098997, 0.08778729;
234 
235 // std::cout << fjac*covfac << std::endl;
236 
237  MatrixXd cov;
238  cov = covfac*lm.matrixR().topLeftCorner<n,n>();
239  VERIFY_IS_APPROX( cov, cov_ref);
240  // TODO: why isn't this allowed ? :
241  // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
242 }
243 
244 struct chwirut2_functor : DenseFunctor<double>
245 {
246  chwirut2_functor(void) : DenseFunctor<double>(3,54) {}
247  static const double m_x[54];
248  static const double m_y[54];
249  int operator()(const VectorXd &b, VectorXd &fvec)
250  {
251  int i;
252 
253  assert(b.size()==3);
254  assert(fvec.size()==54);
255  for(i=0; i<54; i++) {
256  double x = m_x[i];
257  fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i];
258  }
259  return 0;
260  }
261  int df(const VectorXd &b, MatrixXd &fjac)
262  {
263  assert(b.size()==3);
264  assert(fjac.rows()==54);
265  assert(fjac.cols()==3);
266  for(int i=0; i<54; i++) {
267  double x = m_x[i];
268  double factor = 1./(b[1]+b[2]*x);
269  double e = exp(-b[0]*x);
270  fjac(i,0) = -x*e*factor;
271  fjac(i,1) = -e*factor*factor;
272  fjac(i,2) = -x*e*factor*factor;
273  }
274  return 0;
275  }
276 };
277 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};
278 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 };
279 
280 // http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml
281 void testNistChwirut2(void)
282 {
283  const int n=3;
284  LevenbergMarquardtSpace::Status info;
285 
286  VectorXd x(n);
287 
288  /*
289  * First try
290  */
291  x<< 0.1, 0.01, 0.02;
292  // do the computation
293  chwirut2_functor functor;
294  LevenbergMarquardt<chwirut2_functor> lm(functor);
295  info = lm.minimize(x);
296 
297  // check return value
298  VERIFY_IS_EQUAL(info, 1);
299 // VERIFY_IS_EQUAL(lm.nfev(), 10);
300  VERIFY_IS_EQUAL(lm.njev(), 8);
301  // check norm^2
302  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
303  // check x
304  VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
305  VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
306  VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
307 
308  /*
309  * Second try
310  */
311  x<< 0.15, 0.008, 0.010;
312  // do the computation
313  lm.resetParameters();
314  lm.setFtol(1.E6*NumTraits<double>::epsilon());
315  lm.setXtol(1.E6*NumTraits<double>::epsilon());
316  info = lm.minimize(x);
317 
318  // check return value
319  VERIFY_IS_EQUAL(info, 1);
320 // VERIFY_IS_EQUAL(lm.nfev(), 7);
321  VERIFY_IS_EQUAL(lm.njev(), 6);
322  // check norm^2
323  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
324  // check x
325  VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
326  VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
327  VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
328 }
329 
330 
331 struct misra1a_functor : DenseFunctor<double>
332 {
333  misra1a_functor(void) : DenseFunctor<double>(2,14) {}
334  static const double m_x[14];
335  static const double m_y[14];
336  int operator()(const VectorXd &b, VectorXd &fvec)
337  {
338  assert(b.size()==2);
339  assert(fvec.size()==14);
340  for(int i=0; i<14; i++) {
341  fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ;
342  }
343  return 0;
344  }
345  int df(const VectorXd &b, MatrixXd &fjac)
346  {
347  assert(b.size()==2);
348  assert(fjac.rows()==14);
349  assert(fjac.cols()==2);
350  for(int i=0; i<14; i++) {
351  fjac(i,0) = (1.-exp(-b[1]*m_x[i]));
352  fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i]));
353  }
354  return 0;
355  }
356 };
357 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};
358 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};
359 
360 // http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
361 void testNistMisra1a(void)
362 {
363  const int n=2;
364  int info;
365 
366  VectorXd x(n);
367 
368  /*
369  * First try
370  */
371  x<< 500., 0.0001;
372  // do the computation
373  misra1a_functor functor;
374  LevenbergMarquardt<misra1a_functor> lm(functor);
375  info = lm.minimize(x);
376 
377  // check return value
378  VERIFY_IS_EQUAL(info, 1);
379  VERIFY_IS_EQUAL(lm.nfev(), 19);
380  VERIFY_IS_EQUAL(lm.njev(), 15);
381  // check norm^2
382  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01);
383  // check x
384  VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
385  VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
386 
387  /*
388  * Second try
389  */
390  x<< 250., 0.0005;
391  // do the computation
392  info = lm.minimize(x);
393 
394  // check return value
395  VERIFY_IS_EQUAL(info, 1);
396  VERIFY_IS_EQUAL(lm.nfev(), 5);
397  VERIFY_IS_EQUAL(lm.njev(), 4);
398  // check norm^2
399  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01);
400  // check x
401  VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
402  VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
403 }
404 
405 struct hahn1_functor : DenseFunctor<double>
406 {
407  hahn1_functor(void) : DenseFunctor<double>(7,236) {}
408  static const double m_x[236];
409  int operator()(const VectorXd &b, VectorXd &fvec)
410  {
411  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 ,
412  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 ,
413 12.596E0 ,
414 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 };
415 
416  // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
417 
418  assert(b.size()==7);
419  assert(fvec.size()==236);
420  for(int i=0; i<236; i++) {
421  double x=m_x[i], xx=x*x, xxx=xx*x;
422  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];
423  }
424  return 0;
425  }
426 
427  int df(const VectorXd &b, MatrixXd &fjac)
428  {
429  assert(b.size()==7);
430  assert(fjac.rows()==236);
431  assert(fjac.cols()==7);
432  for(int i=0; i<236; i++) {
433  double x=m_x[i], xx=x*x, xxx=xx*x;
434  double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
435  fjac(i,0) = 1.*fact;
436  fjac(i,1) = x*fact;
437  fjac(i,2) = xx*fact;
438  fjac(i,3) = xxx*fact;
439  fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
440  fjac(i,4) = x*fact;
441  fjac(i,5) = xx*fact;
442  fjac(i,6) = xxx*fact;
443  }
444  return 0;
445  }
446 };
447 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 ,
448 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 ,
449 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};
450 
451 // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml
452 void testNistHahn1(void)
453 {
454  const int n=7;
455  int info;
456 
457  VectorXd x(n);
458 
459  /*
460  * First try
461  */
462  x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
463  // do the computation
464  hahn1_functor functor;
465  LevenbergMarquardt<hahn1_functor> lm(functor);
466  info = lm.minimize(x);
467 
468  // check return value
469  VERIFY_IS_EQUAL(info, 1);
470  VERIFY_IS_EQUAL(lm.nfev(), 11);
471  VERIFY_IS_EQUAL(lm.njev(), 10);
472  // check norm^2
473  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00);
474  // check x
475  VERIFY_IS_APPROX(x[0], 1.0776351733E+00);
476  VERIFY_IS_APPROX(x[1],-1.2269296921E-01);
477  VERIFY_IS_APPROX(x[2], 4.0863750610E-03);
478  VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06
479  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
480  VERIFY_IS_APPROX(x[5], 2.4053735503E-04);
481  VERIFY_IS_APPROX(x[6],-1.2314450199E-07);
482 
483  /*
484  * Second try
485  */
486  x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001;
487  // do the computation
488  info = lm.minimize(x);
489 
490  // check return value
491  VERIFY_IS_EQUAL(info, 1);
492 // VERIFY_IS_EQUAL(lm.nfev(), 11);
493  VERIFY_IS_EQUAL(lm.njev(), 10);
494  // check norm^2
495  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00);
496  // check x
497  VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00
498  VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01
499  VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03
500  VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06
501  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
502  VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04
503  VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07
504 
505 }
506 
507 struct misra1d_functor : DenseFunctor<double>
508 {
509  misra1d_functor(void) : DenseFunctor<double>(2,14) {}
510  static const double x[14];
511  static const double y[14];
512  int operator()(const VectorXd &b, VectorXd &fvec)
513  {
514  assert(b.size()==2);
515  assert(fvec.size()==14);
516  for(int i=0; i<14; i++) {
517  fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i];
518  }
519  return 0;
520  }
521  int df(const VectorXd &b, MatrixXd &fjac)
522  {
523  assert(b.size()==2);
524  assert(fjac.rows()==14);
525  assert(fjac.cols()==2);
526  for(int i=0; i<14; i++) {
527  double den = 1.+b[1]*x[i];
528  fjac(i,0) = b[1]*x[i] / den;
529  fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den;
530  }
531  return 0;
532  }
533 };
534 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};
535 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};
536 
537 // http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml
538 void testNistMisra1d(void)
539 {
540  const int n=2;
541  int info;
542 
543  VectorXd x(n);
544 
545  /*
546  * First try
547  */
548  x<< 500., 0.0001;
549  // do the computation
550  misra1d_functor functor;
551  LevenbergMarquardt<misra1d_functor> lm(functor);
552  info = lm.minimize(x);
553 
554  // check return value
555  VERIFY_IS_EQUAL(info, 1);
556  VERIFY_IS_EQUAL(lm.nfev(), 9);
557  VERIFY_IS_EQUAL(lm.njev(), 7);
558  // check norm^2
559  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02);
560  // check x
561  VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
562  VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
563 
564  /*
565  * Second try
566  */
567  x<< 450., 0.0003;
568  // do the computation
569  info = lm.minimize(x);
570 
571  // check return value
572  VERIFY_IS_EQUAL(info, 1);
573  VERIFY_IS_EQUAL(lm.nfev(), 4);
574  VERIFY_IS_EQUAL(lm.njev(), 3);
575  // check norm^2
576  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02);
577  // check x
578  VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
579  VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
580 }
581 
582 
583 struct lanczos1_functor : DenseFunctor<double>
584 {
585  lanczos1_functor(void) : DenseFunctor<double>(6,24) {}
586  static const double x[24];
587  static const double y[24];
588  int operator()(const VectorXd &b, VectorXd &fvec)
589  {
590  assert(b.size()==6);
591  assert(fvec.size()==24);
592  for(int i=0; i<24; i++)
593  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];
594  return 0;
595  }
596  int df(const VectorXd &b, MatrixXd &fjac)
597  {
598  assert(b.size()==6);
599  assert(fjac.rows()==24);
600  assert(fjac.cols()==6);
601  for(int i=0; i<24; i++) {
602  fjac(i,0) = exp(-b[1]*x[i]);
603  fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]);
604  fjac(i,2) = exp(-b[3]*x[i]);
605  fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]);
606  fjac(i,4) = exp(-b[5]*x[i]);
607  fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]);
608  }
609  return 0;
610  }
611 };
612 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 };
613 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 };
614 
615 // http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml
616 void testNistLanczos1(void)
617 {
618  const int n=6;
619  LevenbergMarquardtSpace::Status info;
620 
621  VectorXd x(n);
622 
623  /*
624  * First try
625  */
626  x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6;
627  // do the computation
628  lanczos1_functor functor;
629  LevenbergMarquardt<lanczos1_functor> lm(functor);
630  info = lm.minimize(x);
631 
632  // check return value
633  VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall);
634  VERIFY_IS_EQUAL(lm.nfev(), 79);
635  VERIFY_IS_EQUAL(lm.njev(), 72);
636  // check norm^2
637  VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25);
638  // check x
639  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
640  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
641  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
642  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
643  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
644  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
645 
646  /*
647  * Second try
648  */
649  x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3;
650  // do the computation
651  info = lm.minimize(x);
652 
653  // check return value
654  VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall);
655  VERIFY_IS_EQUAL(lm.nfev(), 9);
656  VERIFY_IS_EQUAL(lm.njev(), 8);
657  // check norm^2
658  VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25);
659  // check x
660  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
661  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
662  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
663  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
664  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
665  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
666 
667 }
668 
669 struct rat42_functor : DenseFunctor<double>
670 {
671  rat42_functor(void) : DenseFunctor<double>(3,9) {}
672  static const double x[9];
673  static const double y[9];
674  int operator()(const VectorXd &b, VectorXd &fvec)
675  {
676  assert(b.size()==3);
677  assert(fvec.size()==9);
678  for(int i=0; i<9; i++) {
679  fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i];
680  }
681  return 0;
682  }
683 
684  int df(const VectorXd &b, MatrixXd &fjac)
685  {
686  assert(b.size()==3);
687  assert(fjac.rows()==9);
688  assert(fjac.cols()==3);
689  for(int i=0; i<9; i++) {
690  double e = exp(b[1]-b[2]*x[i]);
691  fjac(i,0) = 1./(1.+e);
692  fjac(i,1) = -b[0]*e/(1.+e)/(1.+e);
693  fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e);
694  }
695  return 0;
696  }
697 };
698 const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 };
699 const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 };
700 
701 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml
702 void testNistRat42(void)
703 {
704  const int n=3;
705  LevenbergMarquardtSpace::Status info;
706 
707  VectorXd x(n);
708 
709  /*
710  * First try
711  */
712  x<< 100., 1., 0.1;
713  // do the computation
714  rat42_functor functor;
715  LevenbergMarquardt<rat42_functor> lm(functor);
716  info = lm.minimize(x);
717 
718  // check return value
719  VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
720  VERIFY_IS_EQUAL(lm.nfev(), 10);
721  VERIFY_IS_EQUAL(lm.njev(), 8);
722  // check norm^2
723  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00);
724  // check x
725  VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
726  VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
727  VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
728 
729  /*
730  * Second try
731  */
732  x<< 75., 2.5, 0.07;
733  // do the computation
734  info = lm.minimize(x);
735 
736  // check return value
737  VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
738  VERIFY_IS_EQUAL(lm.nfev(), 6);
739  VERIFY_IS_EQUAL(lm.njev(), 5);
740  // check norm^2
741  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00);
742  // check x
743  VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
744  VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
745  VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
746 }
747 
748 struct MGH10_functor : DenseFunctor<double>
749 {
750  MGH10_functor(void) : DenseFunctor<double>(3,16) {}
751  static const double x[16];
752  static const double y[16];
753  int operator()(const VectorXd &b, VectorXd &fvec)
754  {
755  assert(b.size()==3);
756  assert(fvec.size()==16);
757  for(int i=0; i<16; i++)
758  fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i];
759  return 0;
760  }
761  int df(const VectorXd &b, MatrixXd &fjac)
762  {
763  assert(b.size()==3);
764  assert(fjac.rows()==16);
765  assert(fjac.cols()==3);
766  for(int i=0; i<16; i++) {
767  double factor = 1./(x[i]+b[2]);
768  double e = exp(b[1]*factor);
769  fjac(i,0) = e;
770  fjac(i,1) = b[0]*factor*e;
771  fjac(i,2) = -b[1]*b[0]*factor*factor*e;
772  }
773  return 0;
774  }
775 };
776 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 };
777 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 };
778 
779 // http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
780 void testNistMGH10(void)
781 {
782  const int n=3;
783  LevenbergMarquardtSpace::Status info;
784 
785  VectorXd x(n);
786 
787  /*
788  * First try
789  */
790  x<< 2., 400000., 25000.;
791  // do the computation
792  MGH10_functor functor;
793  LevenbergMarquardt<MGH10_functor> lm(functor);
794  info = lm.minimize(x);
795  ++g_test_level;
796  VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
797  --g_test_level;
798  // was: VERIFY_IS_EQUAL(info, 1);
799 
800  // check norm^2
801  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
802  // check x
803  VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
804  VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
805  VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
806 
807  // check return value
808 
809  ++g_test_level;
810  VERIFY_IS_EQUAL(lm.nfev(), 284 );
811  VERIFY_IS_EQUAL(lm.njev(), 249 );
812  --g_test_level;
813  VERIFY(lm.nfev() < 284 * LM_EVAL_COUNT_TOL);
814  VERIFY(lm.njev() < 249 * LM_EVAL_COUNT_TOL);
815 
816  /*
817  * Second try
818  */
819  x<< 0.02, 4000., 250.;
820  // do the computation
821  info = lm.minimize(x);
822  ++g_test_level;
823  VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
824  // was: VERIFY_IS_EQUAL(info, 1);
825  --g_test_level;
826 
827  // check norm^2
828  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
829  // check x
830  VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
831  VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
832  VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
833 
834  // check return value
835  ++g_test_level;
836  VERIFY_IS_EQUAL(lm.nfev(), 126);
837  VERIFY_IS_EQUAL(lm.njev(), 116);
838  --g_test_level;
839  VERIFY(lm.nfev() < 126 * LM_EVAL_COUNT_TOL);
840  VERIFY(lm.njev() < 116 * LM_EVAL_COUNT_TOL);
841 }
842 
843 
844 struct BoxBOD_functor : DenseFunctor<double>
845 {
846  BoxBOD_functor(void) : DenseFunctor<double>(2,6) {}
847  static const double x[6];
848  int operator()(const VectorXd &b, VectorXd &fvec)
849  {
850  static const double y[6] = { 109., 149., 149., 191., 213., 224. };
851  assert(b.size()==2);
852  assert(fvec.size()==6);
853  for(int i=0; i<6; i++)
854  fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i];
855  return 0;
856  }
857  int df(const VectorXd &b, MatrixXd &fjac)
858  {
859  assert(b.size()==2);
860  assert(fjac.rows()==6);
861  assert(fjac.cols()==2);
862  for(int i=0; i<6; i++) {
863  double e = exp(-b[1]*x[i]);
864  fjac(i,0) = 1.-e;
865  fjac(i,1) = b[0]*x[i]*e;
866  }
867  return 0;
868  }
869 };
870 const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. };
871 
872 // http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml
873 void testNistBoxBOD(void)
874 {
875  const int n=2;
876  int info;
877 
878  VectorXd x(n);
879 
880  /*
881  * First try
882  */
883  x<< 1., 1.;
884  // do the computation
885  BoxBOD_functor functor;
886  LevenbergMarquardt<BoxBOD_functor> lm(functor);
887  lm.setFtol(1.E6*NumTraits<double>::epsilon());
888  lm.setXtol(1.E6*NumTraits<double>::epsilon());
889  lm.setFactor(10);
890  info = lm.minimize(x);
891 
892  // check norm^2
893  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03);
894  // check x
895  VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
896  VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
897 
898  // check return value
899  VERIFY_IS_EQUAL(info, 1);
900  VERIFY(lm.nfev() < 31); // 31
901  VERIFY(lm.njev() < 25); // 25
902 
903  /*
904  * Second try
905  */
906  x<< 100., 0.75;
907  // do the computation
908  lm.resetParameters();
909  lm.setFtol(NumTraits<double>::epsilon());
910  lm.setXtol( NumTraits<double>::epsilon());
911  info = lm.minimize(x);
912 
913  // check return value
914  VERIFY_IS_EQUAL(info, 1);
915  ++g_test_level;
916  VERIFY_IS_EQUAL(lm.nfev(), 16 );
917  VERIFY_IS_EQUAL(lm.njev(), 15 );
918  --g_test_level;
919  VERIFY(lm.nfev() < 16 * LM_EVAL_COUNT_TOL);
920  VERIFY(lm.njev() < 15 * LM_EVAL_COUNT_TOL);
921  // check norm^2
922  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03);
923  // check x
924  VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
925  VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
926 }
927 
928 struct MGH17_functor : DenseFunctor<double>
929 {
930  MGH17_functor(void) : DenseFunctor<double>(5,33) {}
931  static const double x[33];
932  static const double y[33];
933  int operator()(const VectorXd &b, VectorXd &fvec)
934  {
935  assert(b.size()==5);
936  assert(fvec.size()==33);
937  for(int i=0; i<33; i++)
938  fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i];
939  return 0;
940  }
941  int df(const VectorXd &b, MatrixXd &fjac)
942  {
943  assert(b.size()==5);
944  assert(fjac.rows()==33);
945  assert(fjac.cols()==5);
946  for(int i=0; i<33; i++) {
947  fjac(i,0) = 1.;
948  fjac(i,1) = exp(-b[3]*x[i]);
949  fjac(i,2) = exp(-b[4]*x[i]);
950  fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]);
951  fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]);
952  }
953  return 0;
954  }
955 };
956 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 };
957 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 };
958 
959 // http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml
960 void testNistMGH17(void)
961 {
962  const int n=5;
963  int info;
964 
965  VectorXd x(n);
966 
967  /*
968  * First try
969  */
970  x<< 50., 150., -100., 1., 2.;
971  // do the computation
972  MGH17_functor functor;
973  LevenbergMarquardt<MGH17_functor> lm(functor);
974  lm.setFtol(NumTraits<double>::epsilon());
975  lm.setXtol(NumTraits<double>::epsilon());
976  lm.setMaxfev(1000);
977  info = lm.minimize(x);
978 
979  // check norm^2
980  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05);
981  // check x
982  VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
983  VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
984  VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
985  VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
986  VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
987 
988  // check return value
989 // VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success)
990  VERIFY(lm.nfev() < 700 ); // 602
991  VERIFY(lm.njev() < 600 ); // 545
992 
993  /*
994  * Second try
995  */
996  x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02;
997  // do the computation
998  lm.resetParameters();
999  info = lm.minimize(x);
1000 
1001  // check return value
1002  VERIFY_IS_EQUAL(info, 1);
1003  VERIFY_IS_EQUAL(lm.nfev(), 18);
1004  VERIFY_IS_EQUAL(lm.njev(), 15);
1005  // check norm^2
1006  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05);
1007  // check x
1008  VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
1009  VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
1010  VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
1011  VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
1012  VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
1013 }
1014 
1015 struct MGH09_functor : DenseFunctor<double>
1016 {
1017  MGH09_functor(void) : DenseFunctor<double>(4,11) {}
1018  static const double _x[11];
1019  static const double y[11];
1020  int operator()(const VectorXd &b, VectorXd &fvec)
1021  {
1022  assert(b.size()==4);
1023  assert(fvec.size()==11);
1024  for(int i=0; i<11; i++) {
1025  double x = _x[i], xx=x*x;
1026  fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i];
1027  }
1028  return 0;
1029  }
1030  int df(const VectorXd &b, MatrixXd &fjac)
1031  {
1032  assert(b.size()==4);
1033  assert(fjac.rows()==11);
1034  assert(fjac.cols()==4);
1035  for(int i=0; i<11; i++) {
1036  double x = _x[i], xx=x*x;
1037  double factor = 1./(xx+x*b[2]+b[3]);
1038  fjac(i,0) = (xx+x*b[1]) * factor;
1039  fjac(i,1) = b[0]*x* factor;
1040  fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor;
1041  fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor;
1042  }
1043  return 0;
1044  }
1045 };
1046 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 };
1047 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 };
1048 
1049 // http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
1050 void testNistMGH09(void)
1051 {
1052  const int n=4;
1053  int info;
1054 
1055  VectorXd x(n);
1056 
1057  /*
1058  * First try
1059  */
1060  x<< 25., 39, 41.5, 39.;
1061  // do the computation
1062  MGH09_functor functor;
1063  LevenbergMarquardt<MGH09_functor> lm(functor);
1064  lm.setMaxfev(1000);
1065  info = lm.minimize(x);
1066 
1067  // check norm^2
1068  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04);
1069  // check x
1070  VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01
1071  VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01
1072  VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01
1073  VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01
1074  // check return value
1075  VERIFY_IS_EQUAL(info, 1);
1076  VERIFY(lm.nfev() < 510 ); // 490
1077  VERIFY(lm.njev() < 400 ); // 376
1078 
1079  /*
1080  * Second try
1081  */
1082  x<< 0.25, 0.39, 0.415, 0.39;
1083  // do the computation
1084  lm.resetParameters();
1085  info = lm.minimize(x);
1086 
1087  // check return value
1088  VERIFY_IS_EQUAL(info, 1);
1089  VERIFY_IS_EQUAL(lm.nfev(), 18);
1090  VERIFY_IS_EQUAL(lm.njev(), 16);
1091  // check norm^2
1092  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04);
1093  // check x
1094  VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01
1095  VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01
1096  VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01
1097  VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01
1098 }
1099 
1100 
1101 
1102 struct Bennett5_functor : DenseFunctor<double>
1103 {
1104  Bennett5_functor(void) : DenseFunctor<double>(3,154) {}
1105  static const double x[154];
1106  static const double y[154];
1107  int operator()(const VectorXd &b, VectorXd &fvec)
1108  {
1109  assert(b.size()==3);
1110  assert(fvec.size()==154);
1111  for(int i=0; i<154; i++)
1112  fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i];
1113  return 0;
1114  }
1115  int df(const VectorXd &b, MatrixXd &fjac)
1116  {
1117  assert(b.size()==3);
1118  assert(fjac.rows()==154);
1119  assert(fjac.cols()==3);
1120  for(int i=0; i<154; i++) {
1121  double e = pow(b[1]+x[i],-1./b[2]);
1122  fjac(i,0) = e;
1123  fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]);
1124  fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2];
1125  }
1126  return 0;
1127  }
1128 };
1129 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,
1130 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 };
1131 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
1132 ,-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 ,
1133 -31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 };
1134 
1135 // http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml
1136 void testNistBennett5(void)
1137 {
1138  const int n=3;
1139  int info;
1140 
1141  VectorXd x(n);
1142 
1143  /*
1144  * First try
1145  */
1146  x<< -2000., 50., 0.8;
1147  // do the computation
1148  Bennett5_functor functor;
1149  LevenbergMarquardt<Bennett5_functor> lm(functor);
1150  lm.setMaxfev(1000);
1151  info = lm.minimize(x);
1152 
1153  // check return value
1154  VERIFY_IS_EQUAL(info, 1);
1155  VERIFY_IS_EQUAL(lm.nfev(), 758);
1156  VERIFY_IS_EQUAL(lm.njev(), 744);
1157  // check norm^2
1158  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04);
1159  // check x
1160  VERIFY_IS_APPROX(x[0], -2.5235058043E+03);
1161  VERIFY_IS_APPROX(x[1], 4.6736564644E+01);
1162  VERIFY_IS_APPROX(x[2], 9.3218483193E-01);
1163  /*
1164  * Second try
1165  */
1166  x<< -1500., 45., 0.85;
1167  // do the computation
1168  lm.resetParameters();
1169  info = lm.minimize(x);
1170 
1171  // check return value
1172  VERIFY_IS_EQUAL(info, 1);
1173  VERIFY_IS_EQUAL(lm.nfev(), 203);
1174  VERIFY_IS_EQUAL(lm.njev(), 192);
1175  // check norm^2
1176  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04);
1177  // check x
1178  VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03
1179  VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01);
1180  VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01);
1181 }
1182 
1183 struct thurber_functor : DenseFunctor<double>
1184 {
1185  thurber_functor(void) : DenseFunctor<double>(7,37) {}
1186  static const double _x[37];
1187  static const double _y[37];
1188  int operator()(const VectorXd &b, VectorXd &fvec)
1189  {
1190  // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
1191  assert(b.size()==7);
1192  assert(fvec.size()==37);
1193  for(int i=0; i<37; i++) {
1194  double x=_x[i], xx=x*x, xxx=xx*x;
1195  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];
1196  }
1197  return 0;
1198  }
1199  int df(const VectorXd &b, MatrixXd &fjac)
1200  {
1201  assert(b.size()==7);
1202  assert(fjac.rows()==37);
1203  assert(fjac.cols()==7);
1204  for(int i=0; i<37; i++) {
1205  double x=_x[i], xx=x*x, xxx=xx*x;
1206  double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
1207  fjac(i,0) = 1.*fact;
1208  fjac(i,1) = x*fact;
1209  fjac(i,2) = xx*fact;
1210  fjac(i,3) = xxx*fact;
1211  fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
1212  fjac(i,4) = x*fact;
1213  fjac(i,5) = xx*fact;
1214  fjac(i,6) = xxx*fact;
1215  }
1216  return 0;
1217  }
1218 };
1219 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 };
1220 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};
1221 
1222 // http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml
1223 void testNistThurber(void)
1224 {
1225  const int n=7;
1226  int info;
1227 
1228  VectorXd x(n);
1229 
1230  /*
1231  * First try
1232  */
1233  x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ;
1234  // do the computation
1235  thurber_functor functor;
1236  LevenbergMarquardt<thurber_functor> lm(functor);
1237  lm.setFtol(1.E4*NumTraits<double>::epsilon());
1238  lm.setXtol(1.E4*NumTraits<double>::epsilon());
1239  info = lm.minimize(x);
1240 
1241  // check return value
1242  VERIFY_IS_EQUAL(info, 1);
1243  VERIFY_IS_EQUAL(lm.nfev(), 39);
1244  VERIFY_IS_EQUAL(lm.njev(), 36);
1245  // check norm^2
1246  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03);
1247  // check x
1248  VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
1249  VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
1250  VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
1251  VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
1252  VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
1253  VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
1254  VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
1255 
1256  /*
1257  * Second try
1258  */
1259  x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ;
1260  // do the computation
1261  lm.resetParameters();
1262  lm.setFtol(1.E4*NumTraits<double>::epsilon());
1263  lm.setXtol(1.E4*NumTraits<double>::epsilon());
1264  info = lm.minimize(x);
1265 
1266  // check return value
1267  VERIFY_IS_EQUAL(info, 1);
1268  VERIFY_IS_EQUAL(lm.nfev(), 29);
1269  VERIFY_IS_EQUAL(lm.njev(), 28);
1270  // check norm^2
1271  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03);
1272  // check x
1273  VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
1274  VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
1275  VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
1276  VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
1277  VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
1278  VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
1279  VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
1280 }
1281 
1282 struct rat43_functor : DenseFunctor<double>
1283 {
1284  rat43_functor(void) : DenseFunctor<double>(4,15) {}
1285  static const double x[15];
1286  static const double y[15];
1287  int operator()(const VectorXd &b, VectorXd &fvec)
1288  {
1289  assert(b.size()==4);
1290  assert(fvec.size()==15);
1291  for(int i=0; i<15; i++)
1292  fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i];
1293  return 0;
1294  }
1295  int df(const VectorXd &b, MatrixXd &fjac)
1296  {
1297  assert(b.size()==4);
1298  assert(fjac.rows()==15);
1299  assert(fjac.cols()==4);
1300  for(int i=0; i<15; i++) {
1301  double e = exp(b[1]-b[2]*x[i]);
1302  double power = -1./b[3];
1303  fjac(i,0) = pow(1.+e, power);
1304  fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.);
1305  fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.);
1306  fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power);
1307  }
1308  return 0;
1309  }
1310 };
1311 const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. };
1312 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 };
1313 
1314 // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml
1315 void testNistRat43(void)
1316 {
1317  const int n=4;
1318  int info;
1319 
1320  VectorXd x(n);
1321 
1322  /*
1323  * First try
1324  */
1325  x<< 100., 10., 1., 1.;
1326  // do the computation
1327  rat43_functor functor;
1328  LevenbergMarquardt<rat43_functor> lm(functor);
1329  lm.setFtol(1.E6*NumTraits<double>::epsilon());
1330  lm.setXtol(1.E6*NumTraits<double>::epsilon());
1331  info = lm.minimize(x);
1332 
1333  // check return value
1334  VERIFY_IS_EQUAL(info, 1);
1335  VERIFY_IS_EQUAL(lm.nfev(), 27);
1336  VERIFY_IS_EQUAL(lm.njev(), 20);
1337  // check norm^2
1338  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03);
1339  // check x
1340  VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
1341  VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
1342  VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
1343  VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
1344 
1345  /*
1346  * Second try
1347  */
1348  x<< 700., 5., 0.75, 1.3;
1349  // do the computation
1350  lm.resetParameters();
1351  lm.setFtol(1.E5*NumTraits<double>::epsilon());
1352  lm.setXtol(1.E5*NumTraits<double>::epsilon());
1353  info = lm.minimize(x);
1354 
1355  // check return value
1356  VERIFY_IS_EQUAL(info, 1);
1357  VERIFY_IS_EQUAL(lm.nfev(), 9);
1358  VERIFY_IS_EQUAL(lm.njev(), 8);
1359  // check norm^2
1360  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03);
1361  // check x
1362  VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
1363  VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
1364  VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
1365  VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
1366 }
1367 
1368 
1369 
1370 struct eckerle4_functor : DenseFunctor<double>
1371 {
1372  eckerle4_functor(void) : DenseFunctor<double>(3,35) {}
1373  static const double x[35];
1374  static const double y[35];
1375  int operator()(const VectorXd &b, VectorXd &fvec)
1376  {
1377  assert(b.size()==3);
1378  assert(fvec.size()==35);
1379  for(int i=0; i<35; i++)
1380  fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i];
1381  return 0;
1382  }
1383  int df(const VectorXd &b, MatrixXd &fjac)
1384  {
1385  assert(b.size()==3);
1386  assert(fjac.rows()==35);
1387  assert(fjac.cols()==3);
1388  for(int i=0; i<35; i++) {
1389  double b12 = b[1]*b[1];
1390  double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12);
1391  fjac(i,0) = e / b[1];
1392  fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12;
1393  fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12;
1394  }
1395  return 0;
1396  }
1397 };
1398 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};
1399 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 };
1400 
1401 // http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml
1402 void testNistEckerle4(void)
1403 {
1404  const int n=3;
1405  int info;
1406 
1407  VectorXd x(n);
1408 
1409  /*
1410  * First try
1411  */
1412  x<< 1., 10., 500.;
1413  // do the computation
1414  eckerle4_functor functor;
1415  LevenbergMarquardt<eckerle4_functor> lm(functor);
1416  info = lm.minimize(x);
1417 
1418  // check return value
1419  VERIFY_IS_EQUAL(info, 1);
1420  VERIFY_IS_EQUAL(lm.nfev(), 18);
1421  VERIFY_IS_EQUAL(lm.njev(), 15);
1422  // check norm^2
1423  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03);
1424  // check x
1425  VERIFY_IS_APPROX(x[0], 1.5543827178);
1426  VERIFY_IS_APPROX(x[1], 4.0888321754);
1427  VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
1428 
1429  /*
1430  * Second try
1431  */
1432  x<< 1.5, 5., 450.;
1433  // do the computation
1434  info = lm.minimize(x);
1435 
1436  // check return value
1437  VERIFY_IS_EQUAL(info, 1);
1438  VERIFY_IS_EQUAL(lm.nfev(), 7);
1439  VERIFY_IS_EQUAL(lm.njev(), 6);
1440  // check norm^2
1441  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03);
1442  // check x
1443  VERIFY_IS_APPROX(x[0], 1.5543827178);
1444  VERIFY_IS_APPROX(x[1], 4.0888321754);
1445  VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
1446 }
1447 
1448 void test_levenberg_marquardt()
1449 {
1450  // Tests using the examples provided by (c)minpack
1451  CALL_SUBTEST(testLmder1());
1452  CALL_SUBTEST(testLmder());
1453  CALL_SUBTEST(testLmdif1());
1454 // CALL_SUBTEST(testLmstr1());
1455 // CALL_SUBTEST(testLmstr());
1456  CALL_SUBTEST(testLmdif());
1457 
1458  // NIST tests, level of difficulty = "Lower"
1459  CALL_SUBTEST(testNistMisra1a());
1460  CALL_SUBTEST(testNistChwirut2());
1461 
1462  // NIST tests, level of difficulty = "Average"
1463  CALL_SUBTEST(testNistHahn1());
1464  CALL_SUBTEST(testNistMisra1d());
1465  CALL_SUBTEST(testNistMGH17());
1466  CALL_SUBTEST(testNistLanczos1());
1467 
1468 // // NIST tests, level of difficulty = "Higher"
1469  CALL_SUBTEST(testNistRat42());
1470  CALL_SUBTEST(testNistMGH10());
1471  CALL_SUBTEST(testNistBoxBOD());
1472 // CALL_SUBTEST(testNistMGH09());
1473  CALL_SUBTEST(testNistBennett5());
1474  CALL_SUBTEST(testNistThurber());
1475  CALL_SUBTEST(testNistRat43());
1476  CALL_SUBTEST(testNistEckerle4());
1477 }
lmder_functor::lmder_functor
lmder_functor(void)
Definition: levenberg_marquardt.cpp:31
BoxBOD_functor::x
static const double x[6]
Definition: levenberg_marquardt.cpp:847
rat42_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:684
rat43_functor
Definition: levenberg_marquardt.cpp:1282
hahn1_functor::hahn1_functor
hahn1_functor(void)
Definition: levenberg_marquardt.cpp:407
Bennett5_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:1107
MGH09_functor::MGH09_functor
MGH09_functor(void)
Definition: levenberg_marquardt.cpp:1017
misra1d_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:512
BoxBOD_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:857
lanczos1_functor
Definition: levenberg_marquardt.cpp:583
Eigen::half_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half pow(const half &a, const half &b)
Definition: Half.h:407
rat42_functor::rat42_functor
rat42_functor(void)
Definition: levenberg_marquardt.cpp:671
Bennett5_functor::y
static const double y[154]
Definition: levenberg_marquardt.cpp:1106
testLmder1
void testLmder1()
Definition: levenberg_marquardt.cpp:65
testNistRat42
void testNistRat42(void)
Definition: levenberg_marquardt.cpp:702
chwirut2_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:249
MGH09_functor::y
static const double y[11]
Definition: levenberg_marquardt.cpp:1019
misra1a_functor::m_y
static const double m_y[14]
Definition: levenberg_marquardt.cpp:335
thurber_functor::_y
static const double _y[37]
Definition: levenberg_marquardt.cpp:1187
thurber_functor
Definition: levenberg_marquardt.cpp:1183
exp
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:88
testLmdif1
void testLmdif1()
Definition: levenberg_marquardt.cpp:166
misra1a_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:336
MGH17_functor
Definition: levenberg_marquardt.cpp:928
MGH17_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:941
test_levenberg_marquardt
void test_levenberg_marquardt()
Definition: levenberg_marquardt.cpp:1448
eckerle4_functor
Definition: levenberg_marquardt.cpp:1370
BoxBOD_functor::BoxBOD_functor
BoxBOD_functor(void)
Definition: levenberg_marquardt.cpp:846
Bennett5_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:1115
chwirut2_functor
Definition: levenberg_marquardt.cpp:244
rat43_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:1295
misra1a_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:345
log
EIGEN_DEVICE_FUNC const LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:105
MGH10_functor::x
static const double x[16]
Definition: levenberg_marquardt.cpp:751
MGH10_functor
Definition: levenberg_marquardt.cpp:748
sqrt
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Definition: ArrayCwiseUnaryOps.h:152
eckerle4_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:1383
rat43_functor::x
static const double x[15]
Definition: levenberg_marquardt.cpp:1285
MGH10_functor::MGH10_functor
MGH10_functor(void)
Definition: levenberg_marquardt.cpp:750
eckerle4_functor::x
static const double x[35]
Definition: levenberg_marquardt.cpp:1373
testNistRat43
void testNistRat43(void)
Definition: levenberg_marquardt.cpp:1315
testNistHahn1
void testNistHahn1(void)
Definition: levenberg_marquardt.cpp:452
testNistThurber
void testNistThurber(void)
Definition: levenberg_marquardt.cpp:1223
eckerle4_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:1375
rat42_functor
Definition: levenberg_marquardt.cpp:669
MGH17_functor::y
static const double y[33]
Definition: levenberg_marquardt.cpp:932
lmdif_functor::operator()
int operator()(const VectorXd &x, VectorXd &fvec) const
Definition: levenberg_marquardt.cpp:144
misra1a_functor::misra1a_functor
misra1a_functor(void)
Definition: levenberg_marquardt.cpp:333
testLmdif
void testLmdif()
Definition: levenberg_marquardt.cpp:196
Bennett5_functor::x
static const double x[154]
Definition: levenberg_marquardt.cpp:1105
rat42_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:674
misra1d_functor::y
static const double y[14]
Definition: levenberg_marquardt.cpp:511
eckerle4_functor::y
static const double y[35]
Definition: levenberg_marquardt.cpp:1374
eckerle4_functor::eckerle4_functor
eckerle4_functor(void)
Definition: levenberg_marquardt.cpp:1372
testNistMGH09
void testNistMGH09(void)
Definition: levenberg_marquardt.cpp:1050
lmdif_functor::lmdif_functor
lmdif_functor(void)
Definition: levenberg_marquardt.cpp:143
MGH09_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:1020
MGH10_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:753
MGH10_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:761
rat42_functor::y
static const double y[9]
Definition: levenberg_marquardt.cpp:673
Bennett5_functor::Bennett5_functor
Bennett5_functor(void)
Definition: levenberg_marquardt.cpp:1104
misra1d_functor
Definition: levenberg_marquardt.cpp:507
MGH09_functor::_x
static const double _x[11]
Definition: levenberg_marquardt.cpp:1018
misra1d_functor::misra1d_functor
misra1d_functor(void)
Definition: levenberg_marquardt.cpp:509
testNistBennett5
void testNistBennett5(void)
Definition: levenberg_marquardt.cpp:1136
testNistChwirut2
void testNistChwirut2(void)
Definition: levenberg_marquardt.cpp:281
rat42_functor::x
static const double x[9]
Definition: levenberg_marquardt.cpp:672
chwirut2_functor::chwirut2_functor
chwirut2_functor(void)
Definition: levenberg_marquardt.cpp:246
misra1d_functor::x
static const double x[14]
Definition: levenberg_marquardt.cpp:510
hahn1_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:409
lanczos1_functor::y
static const double y[24]
Definition: levenberg_marquardt.cpp:587
rat43_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:1287
MGH17_functor::MGH17_functor
MGH17_functor(void)
Definition: levenberg_marquardt.cpp:930
misra1a_functor::m_x
static const double m_x[14]
Definition: levenberg_marquardt.cpp:334
hahn1_functor::m_x
static const double m_x[236]
Definition: levenberg_marquardt.cpp:408
lmdif_functor
Definition: levenberg_marquardt.cpp:141
testNistMisra1d
void testNistMisra1d(void)
Definition: levenberg_marquardt.cpp:538
hahn1_functor
Definition: levenberg_marquardt.cpp:405
testNistMisra1a
void testNistMisra1a(void)
Definition: levenberg_marquardt.cpp:361
testLmder
void testLmder()
Definition: levenberg_marquardt.cpp:93
rat43_functor::rat43_functor
rat43_functor(void)
Definition: levenberg_marquardt.cpp:1284
chwirut2_functor::m_x
static const double m_x[54]
Definition: levenberg_marquardt.cpp:247
lmder_functor::df
int df(const VectorXd &x, MatrixXd &fjac) const
Definition: levenberg_marquardt.cpp:48
Eigen::DenseIndex
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: Meta.h:25
MGH17_functor::x
static const double x[33]
Definition: levenberg_marquardt.cpp:931
thurber_functor::_x
static const double _x[37]
Definition: levenberg_marquardt.cpp:1186
lanczos1_functor::lanczos1_functor
lanczos1_functor(void)
Definition: levenberg_marquardt.cpp:585
MGH09_functor
Definition: levenberg_marquardt.cpp:1015
testNistBoxBOD
void testNistBoxBOD(void)
Definition: levenberg_marquardt.cpp:873
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: levenberg_marquardt.cpp:933
misra1a_functor
Definition: levenberg_marquardt.cpp:331
lmder_functor::operator()
int operator()(const VectorXd &x, VectorXd &fvec) const
Definition: levenberg_marquardt.cpp:32
lmder_functor
Definition: levenberg_marquardt.cpp:29
lanczos1_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:596
misra1d_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:521
chwirut2_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:261
Bennett5_functor
Definition: levenberg_marquardt.cpp:1102
BoxBOD_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:848
LM_EVAL_COUNT_TOL
#define LM_EVAL_COUNT_TOL
Definition: levenberg_marquardt.cpp:27
testNistEckerle4
void testNistEckerle4(void)
Definition: levenberg_marquardt.cpp:1402
thurber_functor::thurber_functor
thurber_functor(void)
Definition: levenberg_marquardt.cpp:1185
lanczos1_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:588
rat43_functor::y
static const double y[15]
Definition: levenberg_marquardt.cpp:1286
MGH10_functor::y
static const double y[16]
Definition: levenberg_marquardt.cpp:752
BoxBOD_functor
Definition: levenberg_marquardt.cpp:844
hahn1_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:427
Eigen::numext::b
EIGEN_DEVICE_FUNC const Scalar & b
Definition: SpecialFunctionsImpl.h:1556
MGH09_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:1030
thurber_functor::df
int df(const VectorXd &b, MatrixXd &fjac)
Definition: levenberg_marquardt.cpp:1199
testNistMGH17
void testNistMGH17(void)
Definition: levenberg_marquardt.cpp:960
thurber_functor::operator()
int operator()(const VectorXd &b, VectorXd &fvec)
Definition: levenberg_marquardt.cpp:1188
chwirut2_functor::m_y
static const double m_y[54]
Definition: levenberg_marquardt.cpp:248
testNistMGH10
void testNistMGH10(void)
Definition: levenberg_marquardt.cpp:780
testNistLanczos1
void testNistLanczos1(void)
Definition: levenberg_marquardt.cpp:616
lanczos1_functor::x
static const double x[24]
Definition: levenberg_marquardt.cpp:586

hebiros
Author(s): Xavier Artache , Matthew Tesch
autogenerated on Thu Sep 3 2020 04:08:20