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