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.
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;
77  info = lm.lmder1(x);
78 
79  // check return value
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;
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 
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);
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 
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
282 {
283  const int n=3;
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;
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();
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 
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;
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;
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 
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;
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 
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
617 {
618  const int n=6;
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;
630  info = lm.minimize(x);
631 
632  // check return value
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
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;
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;
716  info = lm.minimize(x);
717 
718  // check return value
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
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;
784 
785  VectorXd x(n);
786 
787  /*
788  * First try
789  */
790  x<< 2., 400000., 25000.;
791  // do the computation
792  MGH10_functor functor;
794  info = lm.minimize(x);
795  ++g_test_level;
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;
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;
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();
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;
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;
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 
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
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;
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 
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
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;
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();
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;
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();
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 
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
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;
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 EIGEN_DECLARE_TEST(levenberg_marquardt)
1449 {
1450  // Tests using the examples provided by (c)minpack
1454 // CALL_SUBTEST(testLmstr1());
1455 // CALL_SUBTEST(testLmstr());
1457 
1458  // NIST tests, level of difficulty = "Lower"
1461 
1462  // NIST tests, level of difficulty = "Average"
1467 
1468 // // NIST tests, level of difficulty = "Higher"
1472 // CALL_SUBTEST(testNistMGH09());
1477 }
static const double x[6]
Matrix3f m
int df(const VectorXd &b, MatrixXd &fjac)
LevenbergMarquardtSpace::Status minimize(FVectorType &x)
int operator()(const VectorXd &b, VectorXd &fvec)
int operator()(const VectorXd &b, VectorXd &fvec)
int df(const VectorXd &b, MatrixXd &fjac)
static const double y[154]
void testLmder1()
void testNistRat42(void)
Scalar * y
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 testLmdif1()
int operator()(const VectorXd &b, VectorXd &fvec)
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)
static const double x[16]
int df(const VectorXd &b, MatrixXd &fjac)
static const double x[15]
static const double x[35]
void testNistRat43(void)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
void testNistHahn1(void)
void testNistThurber(void)
int operator()(const VectorXd &b, VectorXd &fvec)
static const double y[33]
EIGEN_DEVICE_FUNC const LogReturnType log() const
void testLmdif()
EIGEN_DECLARE_TEST(levenberg_marquardt)
static const double x[154]
int operator()(const VectorXd &b, VectorXd &fvec)
static const double y[14]
else if n * info
static const double y[35]
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
void testNistMGH09(void)
#define VERIFY_IS_APPROX(a, b)
int operator()(const VectorXd &b, VectorXd &fvec)
int operator()(const VectorXd &b, VectorXd &fvec)
int df(const VectorXd &b, MatrixXd &fjac)
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
static const double y[9]
const IndicesType & indices() const
static const double _x[11]
void testNistBennett5(void)
void testNistChwirut2(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.)
static const double y[24]
int operator()(const VectorXd &b, VectorXd &fvec)
static const double m_x[14]
static const double m_x[236]
void testNistMisra1d(void)
void testNistMisra1a(void)
void testLmder()
static const double m_x[54]
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: Meta.h:66
static const double x[33]
static const double _x[37]
void testNistBoxBOD(void)
void covar(Matrix< Scalar, Dynamic, Dynamic > &r, const VectorXi &ipvt, Scalar tol=std::sqrt(NumTraits< Scalar >::epsilon()))
Definition: LMcovar.h:20
DiscreteKey E(5, 2)
int operator()(const VectorXd &b, VectorXd &fvec)
#define CALL_SUBTEST(FUNC)
Definition: main.h:399
#define VERIFY(a)
Definition: main.h:380
int df(const VectorXd &b, MatrixXd &fjac)
int df(const VectorXd &b, MatrixXd &fjac)
int df(const VectorXd &b, MatrixXd &fjac)
int operator()(const VectorXd &b, VectorXd &fvec)
static int g_test_level
Definition: main.h:168
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
#define LM_EVAL_COUNT_TOL
void testNistEckerle4(void)
int operator()(const VectorXd &x, VectorXd &fvec) const
int operator()(const VectorXd &b, VectorXd &fvec)
static const double y[15]
int df(const VectorXd &x, MatrixXd &fjac) const
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]
int df(const VectorXd &b, MatrixXd &fjac)
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)
int df(const VectorXd &b, MatrixXd &fjac)
void testNistMGH17(void)
int operator()(const VectorXd &b, VectorXd &fvec)
static const double m_y[54]
void testNistMGH10(void)
void testNistLanczos1(void)
int operator()(const VectorXd &x, VectorXd &fvec) const
static const double x[24]


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:32