1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
5 // Copyright (C) 2012 desire Nuentsa <desire.nuentsa_wakam@inria.fr
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11
12 // FIXME: These tests all check for hard-coded values. Ideally, parameters and start estimates should be randomized.
13
14
15 #include <stdio.h>
16
17 #include "main.h"
18 #include <unsupported/Eigen/LevenbergMarquardt>
19
20 // This disables some useless Warnings on MSVC.
21 // It is intended to be done for this test only.
22 #include <Eigen/src/Core/util/DisableStupidWarnings.h>
23
24 using std::sqrt;
25
26 // tolerance for chekcing number of iterations
27 #define LM_EVAL_COUNT_TOL 4/3
28
29 struct lmder_functor : DenseFunctor<double>
30 {
lmder_functorlmder_functor31 lmder_functor(void): DenseFunctor<double>(3,15) {}
operator ()lmder_functor32 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
dflmder_functor48 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
testLmder1()65 void testLmder1()
66 {
67 int n=3, info;
68
69 VectorXd x;
70
71 /* the following starting values provide a rough fit. */
72 x.setConstant(n, 1.);
73
74 // do the computation
75 lmder_functor functor;
76 LevenbergMarquardt<lmder_functor> lm(functor);
77 info = lm.lmder1(x);
78
79 // check return value
80 VERIFY_IS_EQUAL(info, 1);
81 VERIFY_IS_EQUAL(lm.nfev(), 6);
82 VERIFY_IS_EQUAL(lm.njev(), 5);
83
84 // check norm
85 VERIFY_IS_APPROX(lm.fvec().blueNorm(), 0.09063596);
86
87 // check x
88 VectorXd x_ref(n);
89 x_ref << 0.08241058, 1.133037, 2.343695;
90 VERIFY_IS_APPROX(x, x_ref);
91 }
92
testLmder()93 void testLmder()
94 {
95 const int m=15, n=3;
96 int info;
97 double fnorm, covfac;
98 VectorXd x;
99
100 /* the following starting values provide a rough fit. */
101 x.setConstant(n, 1.);
102
103 // do the computation
104 lmder_functor functor;
105 LevenbergMarquardt<lmder_functor> lm(functor);
106 info = lm.minimize(x);
107
108 // check return values
109 VERIFY_IS_EQUAL(info, 1);
110 VERIFY_IS_EQUAL(lm.nfev(), 6);
111 VERIFY_IS_EQUAL(lm.njev(), 5);
112
113 // check norm
114 fnorm = lm.fvec().blueNorm();
115 VERIFY_IS_APPROX(fnorm, 0.09063596);
116
117 // check x
118 VectorXd x_ref(n);
119 x_ref << 0.08241058, 1.133037, 2.343695;
120 VERIFY_IS_APPROX(x, x_ref);
121
122 // check covariance
123 covfac = fnorm*fnorm/(m-n);
124 internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm
125
126 MatrixXd cov_ref(n,n);
127 cov_ref <<
128 0.0001531202, 0.002869941, -0.002656662,
129 0.002869941, 0.09480935, -0.09098995,
130 -0.002656662, -0.09098995, 0.08778727;
131
132 // std::cout << fjac*covfac << std::endl;
133
134 MatrixXd cov;
135 cov = covfac*lm.matrixR().topLeftCorner<n,n>();
136 VERIFY_IS_APPROX( cov, cov_ref);
137 // TODO: why isn't this allowed ? :
138 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
139 }
140
141 struct lmdif_functor : DenseFunctor<double>
142 {
lmdif_functorlmdif_functor143 lmdif_functor(void) : DenseFunctor<double>(3,15) {}
operator ()lmdif_functor144 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
testLmdif1()166 void testLmdif1()
167 {
168 const int n=3;
169 int info;
170
171 VectorXd x(n), fvec(15);
172
173 /* the following starting values provide a rough fit. */
174 x.setConstant(n, 1.);
175
176 // do the computation
177 lmdif_functor functor;
178 DenseIndex nfev;
179 info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev);
180
181 // check return value
182 VERIFY_IS_EQUAL(info, 1);
183 // VERIFY_IS_EQUAL(nfev, 26);
184
185 // check norm
186 functor(x, fvec);
187 VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596);
188
189 // check x
190 VectorXd x_ref(n);
191 x_ref << 0.0824106, 1.1330366, 2.3436947;
192 VERIFY_IS_APPROX(x, x_ref);
193
194 }
195
testLmdif()196 void testLmdif()
197 {
198 const int m=15, n=3;
199 int info;
200 double fnorm, covfac;
201 VectorXd x(n);
202
203 /* the following starting values provide a rough fit. */
204 x.setConstant(n, 1.);
205
206 // do the computation
207 lmdif_functor functor;
208 NumericalDiff<lmdif_functor> numDiff(functor);
209 LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff);
210 info = lm.minimize(x);
211
212 // check return values
213 VERIFY_IS_EQUAL(info, 1);
214 // VERIFY_IS_EQUAL(lm.nfev(), 26);
215
216 // check norm
217 fnorm = lm.fvec().blueNorm();
218 VERIFY_IS_APPROX(fnorm, 0.09063596);
219
220 // check x
221 VectorXd x_ref(n);
222 x_ref << 0.08241058, 1.133037, 2.343695;
223 VERIFY_IS_APPROX(x, x_ref);
224
225 // check covariance
226 covfac = fnorm*fnorm/(m-n);
227 internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm
228
229 MatrixXd cov_ref(n,n);
230 cov_ref <<
231 0.0001531202, 0.002869942, -0.002656662,
232 0.002869942, 0.09480937, -0.09098997,
233 -0.002656662, -0.09098997, 0.08778729;
234
235 // std::cout << fjac*covfac << std::endl;
236
237 MatrixXd cov;
238 cov = covfac*lm.matrixR().topLeftCorner<n,n>();
239 VERIFY_IS_APPROX( cov, cov_ref);
240 // TODO: why isn't this allowed ? :
241 // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
242 }
243
244 struct chwirut2_functor : DenseFunctor<double>
245 {
chwirut2_functorchwirut2_functor246 chwirut2_functor(void) : DenseFunctor<double>(3,54) {}
247 static const double m_x[54];
248 static const double m_y[54];
operator ()chwirut2_functor249 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 }
dfchwirut2_functor261 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
testNistChwirut2(void)281 void testNistChwirut2(void)
282 {
283 const int n=3;
284 LevenbergMarquardtSpace::Status info;
285
286 VectorXd x(n);
287
288 /*
289 * First try
290 */
291 x<< 0.1, 0.01, 0.02;
292 // do the computation
293 chwirut2_functor functor;
294 LevenbergMarquardt<chwirut2_functor> lm(functor);
295 info = lm.minimize(x);
296
297 // check return value
298 VERIFY_IS_EQUAL(info, 1);
299 // VERIFY_IS_EQUAL(lm.nfev(), 10);
300 VERIFY_IS_EQUAL(lm.njev(), 8);
301 // check norm^2
302 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
303 // check x
304 VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
305 VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
306 VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
307
308 /*
309 * Second try
310 */
311 x<< 0.15, 0.008, 0.010;
312 // do the computation
313 lm.resetParameters();
314 lm.setFtol(1.E6*NumTraits<double>::epsilon());
315 lm.setXtol(1.E6*NumTraits<double>::epsilon());
316 info = lm.minimize(x);
317
318 // check return value
319 VERIFY_IS_EQUAL(info, 1);
320 // VERIFY_IS_EQUAL(lm.nfev(), 7);
321 VERIFY_IS_EQUAL(lm.njev(), 6);
322 // check norm^2
323 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
324 // check x
325 VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
326 VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
327 VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
328 }
329
330
331 struct misra1a_functor : DenseFunctor<double>
332 {
misra1a_functormisra1a_functor333 misra1a_functor(void) : DenseFunctor<double>(2,14) {}
334 static const double m_x[14];
335 static const double m_y[14];
operator ()misra1a_functor336 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 }
dfmisra1a_functor345 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
testNistMisra1a(void)361 void testNistMisra1a(void)
362 {
363 const int n=2;
364 int info;
365
366 VectorXd x(n);
367
368 /*
369 * First try
370 */
371 x<< 500., 0.0001;
372 // do the computation
373 misra1a_functor functor;
374 LevenbergMarquardt<misra1a_functor> lm(functor);
375 info = lm.minimize(x);
376
377 // check return value
378 VERIFY_IS_EQUAL(info, 1);
379 VERIFY_IS_EQUAL(lm.nfev(), 19);
380 VERIFY_IS_EQUAL(lm.njev(), 15);
381 // check norm^2
382 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01);
383 // check x
384 VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
385 VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
386
387 /*
388 * Second try
389 */
390 x<< 250., 0.0005;
391 // do the computation
392 info = lm.minimize(x);
393
394 // check return value
395 VERIFY_IS_EQUAL(info, 1);
396 VERIFY_IS_EQUAL(lm.nfev(), 5);
397 VERIFY_IS_EQUAL(lm.njev(), 4);
398 // check norm^2
399 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01);
400 // check x
401 VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
402 VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
403 }
404
405 struct hahn1_functor : DenseFunctor<double>
406 {
hahn1_functorhahn1_functor407 hahn1_functor(void) : DenseFunctor<double>(7,236) {}
408 static const double m_x[236];
operator ()hahn1_functor409 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
dfhahn1_functor427 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
testNistHahn1(void)452 void testNistHahn1(void)
453 {
454 const int n=7;
455 int info;
456
457 VectorXd x(n);
458
459 /*
460 * First try
461 */
462 x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
463 // do the computation
464 hahn1_functor functor;
465 LevenbergMarquardt<hahn1_functor> lm(functor);
466 info = lm.minimize(x);
467
468 // check return value
469 VERIFY_IS_EQUAL(info, 1);
470 VERIFY_IS_EQUAL(lm.nfev(), 11);
471 VERIFY_IS_EQUAL(lm.njev(), 10);
472 // check norm^2
473 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00);
474 // check x
475 VERIFY_IS_APPROX(x[0], 1.0776351733E+00);
476 VERIFY_IS_APPROX(x[1],-1.2269296921E-01);
477 VERIFY_IS_APPROX(x[2], 4.0863750610E-03);
478 VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06
479 VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
480 VERIFY_IS_APPROX(x[5], 2.4053735503E-04);
481 VERIFY_IS_APPROX(x[6],-1.2314450199E-07);
482
483 /*
484 * Second try
485 */
486 x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001;
487 // do the computation
488 info = lm.minimize(x);
489
490 // check return value
491 VERIFY_IS_EQUAL(info, 1);
492 // VERIFY_IS_EQUAL(lm.nfev(), 11);
493 VERIFY_IS_EQUAL(lm.njev(), 10);
494 // check norm^2
495 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00);
496 // check x
497 VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00
498 VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01
499 VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03
500 VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06
501 VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
502 VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04
503 VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07
504
505 }
506
507 struct misra1d_functor : DenseFunctor<double>
508 {
misra1d_functormisra1d_functor509 misra1d_functor(void) : DenseFunctor<double>(2,14) {}
510 static const double x[14];
511 static const double y[14];
operator ()misra1d_functor512 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 }
dfmisra1d_functor521 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
testNistMisra1d(void)538 void testNistMisra1d(void)
539 {
540 const int n=2;
541 int info;
542
543 VectorXd x(n);
544
545 /*
546 * First try
547 */
548 x<< 500., 0.0001;
549 // do the computation
550 misra1d_functor functor;
551 LevenbergMarquardt<misra1d_functor> lm(functor);
552 info = lm.minimize(x);
553
554 // check return value
555 VERIFY_IS_EQUAL(info, 1);
556 VERIFY_IS_EQUAL(lm.nfev(), 9);
557 VERIFY_IS_EQUAL(lm.njev(), 7);
558 // check norm^2
559 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02);
560 // check x
561 VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
562 VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
563
564 /*
565 * Second try
566 */
567 x<< 450., 0.0003;
568 // do the computation
569 info = lm.minimize(x);
570
571 // check return value
572 VERIFY_IS_EQUAL(info, 1);
573 VERIFY_IS_EQUAL(lm.nfev(), 4);
574 VERIFY_IS_EQUAL(lm.njev(), 3);
575 // check norm^2
576 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02);
577 // check x
578 VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
579 VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
580 }
581
582
583 struct lanczos1_functor : DenseFunctor<double>
584 {
lanczos1_functorlanczos1_functor585 lanczos1_functor(void) : DenseFunctor<double>(6,24) {}
586 static const double x[24];
587 static const double y[24];
operator ()lanczos1_functor588 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 }
dflanczos1_functor596 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
testNistLanczos1(void)616 void testNistLanczos1(void)
617 {
618 const int n=6;
619 LevenbergMarquardtSpace::Status info;
620
621 VectorXd x(n);
622
623 /*
624 * First try
625 */
626 x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6;
627 // do the computation
628 lanczos1_functor functor;
629 LevenbergMarquardt<lanczos1_functor> lm(functor);
630 info = lm.minimize(x);
631
632 // check return value
633 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall);
634 VERIFY_IS_EQUAL(lm.nfev(), 79);
635 VERIFY_IS_EQUAL(lm.njev(), 72);
636 // check norm^2
637 VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25);
638 // check x
639 VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
640 VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
641 VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
642 VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
643 VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
644 VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
645
646 /*
647 * Second try
648 */
649 x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3;
650 // do the computation
651 info = lm.minimize(x);
652
653 // check return value
654 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeErrorTooSmall);
655 VERIFY_IS_EQUAL(lm.nfev(), 9);
656 VERIFY_IS_EQUAL(lm.njev(), 8);
657 // check norm^2
658 VERIFY(lm.fvec().squaredNorm() <= 1.4307867721E-25);
659 // check x
660 VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
661 VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
662 VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
663 VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
664 VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
665 VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
666
667 }
668
669 struct rat42_functor : DenseFunctor<double>
670 {
rat42_functorrat42_functor671 rat42_functor(void) : DenseFunctor<double>(3,9) {}
672 static const double x[9];
673 static const double y[9];
operator ()rat42_functor674 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
dfrat42_functor684 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
testNistRat42(void)702 void testNistRat42(void)
703 {
704 const int n=3;
705 LevenbergMarquardtSpace::Status info;
706
707 VectorXd x(n);
708
709 /*
710 * First try
711 */
712 x<< 100., 1., 0.1;
713 // do the computation
714 rat42_functor functor;
715 LevenbergMarquardt<rat42_functor> lm(functor);
716 info = lm.minimize(x);
717
718 // check return value
719 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
720 VERIFY_IS_EQUAL(lm.nfev(), 10);
721 VERIFY_IS_EQUAL(lm.njev(), 8);
722 // check norm^2
723 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00);
724 // check x
725 VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
726 VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
727 VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
728
729 /*
730 * Second try
731 */
732 x<< 75., 2.5, 0.07;
733 // do the computation
734 info = lm.minimize(x);
735
736 // check return value
737 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
738 VERIFY_IS_EQUAL(lm.nfev(), 6);
739 VERIFY_IS_EQUAL(lm.njev(), 5);
740 // check norm^2
741 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00);
742 // check x
743 VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
744 VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
745 VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
746 }
747
748 struct MGH10_functor : DenseFunctor<double>
749 {
MGH10_functorMGH10_functor750 MGH10_functor(void) : DenseFunctor<double>(3,16) {}
751 static const double x[16];
752 static const double y[16];
operator ()MGH10_functor753 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 }
dfMGH10_functor761 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
testNistMGH10(void)780 void testNistMGH10(void)
781 {
782 const int n=3;
783 LevenbergMarquardtSpace::Status info;
784
785 VectorXd x(n);
786
787 /*
788 * First try
789 */
790 x<< 2., 400000., 25000.;
791 // do the computation
792 MGH10_functor functor;
793 LevenbergMarquardt<MGH10_functor> lm(functor);
794 info = lm.minimize(x);
795 ++g_test_level;
796 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
797 --g_test_level;
798 // was: VERIFY_IS_EQUAL(info, 1);
799
800 // check norm^2
801 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
802 // check x
803 VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
804 VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
805 VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
806
807 // check return value
808
809 ++g_test_level;
810 VERIFY_IS_EQUAL(lm.nfev(), 284 );
811 VERIFY_IS_EQUAL(lm.njev(), 249 );
812 --g_test_level;
813 VERIFY(lm.nfev() < 284 * LM_EVAL_COUNT_TOL);
814 VERIFY(lm.njev() < 249 * LM_EVAL_COUNT_TOL);
815
816 /*
817 * Second try
818 */
819 x<< 0.02, 4000., 250.;
820 // do the computation
821 info = lm.minimize(x);
822 ++g_test_level;
823 VERIFY_IS_EQUAL(info, LevenbergMarquardtSpace::RelativeReductionTooSmall);
824 // was: VERIFY_IS_EQUAL(info, 1);
825 --g_test_level;
826
827 // check norm^2
828 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
829 // check x
830 VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
831 VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
832 VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
833
834 // check return value
835 ++g_test_level;
836 VERIFY_IS_EQUAL(lm.nfev(), 126);
837 VERIFY_IS_EQUAL(lm.njev(), 116);
838 --g_test_level;
839 VERIFY(lm.nfev() < 126 * LM_EVAL_COUNT_TOL);
840 VERIFY(lm.njev() < 116 * LM_EVAL_COUNT_TOL);
841 }
842
843
844 struct BoxBOD_functor : DenseFunctor<double>
845 {
BoxBOD_functorBoxBOD_functor846 BoxBOD_functor(void) : DenseFunctor<double>(2,6) {}
847 static const double x[6];
operator ()BoxBOD_functor848 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 }
dfBoxBOD_functor857 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
testNistBoxBOD(void)873 void testNistBoxBOD(void)
874 {
875 const int n=2;
876 int info;
877
878 VectorXd x(n);
879
880 /*
881 * First try
882 */
883 x<< 1., 1.;
884 // do the computation
885 BoxBOD_functor functor;
886 LevenbergMarquardt<BoxBOD_functor> lm(functor);
887 lm.setFtol(1.E6*NumTraits<double>::epsilon());
888 lm.setXtol(1.E6*NumTraits<double>::epsilon());
889 lm.setFactor(10);
890 info = lm.minimize(x);
891
892 // check norm^2
893 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03);
894 // check x
895 VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
896 VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
897
898 // check return value
899 VERIFY_IS_EQUAL(info, 1);
900 VERIFY(lm.nfev() < 31); // 31
901 VERIFY(lm.njev() < 25); // 25
902
903 /*
904 * Second try
905 */
906 x<< 100., 0.75;
907 // do the computation
908 lm.resetParameters();
909 lm.setFtol(NumTraits<double>::epsilon());
910 lm.setXtol( NumTraits<double>::epsilon());
911 info = lm.minimize(x);
912
913 // check return value
914 VERIFY_IS_EQUAL(info, 1);
915 ++g_test_level;
916 VERIFY_IS_EQUAL(lm.nfev(), 16 );
917 VERIFY_IS_EQUAL(lm.njev(), 15 );
918 --g_test_level;
919 VERIFY(lm.nfev() < 16 * LM_EVAL_COUNT_TOL);
920 VERIFY(lm.njev() < 15 * LM_EVAL_COUNT_TOL);
921 // check norm^2
922 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03);
923 // check x
924 VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
925 VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
926 }
927
928 struct MGH17_functor : DenseFunctor<double>
929 {
MGH17_functorMGH17_functor930 MGH17_functor(void) : DenseFunctor<double>(5,33) {}
931 static const double x[33];
932 static const double y[33];
operator ()MGH17_functor933 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 }
dfMGH17_functor941 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
testNistMGH17(void)960 void testNistMGH17(void)
961 {
962 const int n=5;
963 int info;
964
965 VectorXd x(n);
966
967 /*
968 * First try
969 */
970 x<< 50., 150., -100., 1., 2.;
971 // do the computation
972 MGH17_functor functor;
973 LevenbergMarquardt<MGH17_functor> lm(functor);
974 lm.setFtol(NumTraits<double>::epsilon());
975 lm.setXtol(NumTraits<double>::epsilon());
976 lm.setMaxfev(1000);
977 info = lm.minimize(x);
978
979 // check norm^2
980 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05);
981 // check x
982 VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
983 VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
984 VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
985 VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
986 VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
987
988 // check return value
989 // VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success)
990 VERIFY(lm.nfev() < 700 ); // 602
991 VERIFY(lm.njev() < 600 ); // 545
992
993 /*
994 * Second try
995 */
996 x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02;
997 // do the computation
998 lm.resetParameters();
999 info = lm.minimize(x);
1000
1001 // check return value
1002 VERIFY_IS_EQUAL(info, 1);
1003 VERIFY_IS_EQUAL(lm.nfev(), 18);
1004 VERIFY_IS_EQUAL(lm.njev(), 15);
1005 // check norm^2
1006 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05);
1007 // check x
1008 VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
1009 VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
1010 VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
1011 VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
1012 VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
1013 }
1014
1015 struct MGH09_functor : DenseFunctor<double>
1016 {
MGH09_functorMGH09_functor1017 MGH09_functor(void) : DenseFunctor<double>(4,11) {}
1018 static const double _x[11];
1019 static const double y[11];
operator ()MGH09_functor1020 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 }
dfMGH09_functor1030 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
testNistMGH09(void)1050 void testNistMGH09(void)
1051 {
1052 const int n=4;
1053 int info;
1054
1055 VectorXd x(n);
1056
1057 /*
1058 * First try
1059 */
1060 x<< 25., 39, 41.5, 39.;
1061 // do the computation
1062 MGH09_functor functor;
1063 LevenbergMarquardt<MGH09_functor> lm(functor);
1064 lm.setMaxfev(1000);
1065 info = lm.minimize(x);
1066
1067 // check norm^2
1068 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04);
1069 // check x
1070 VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01
1071 VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01
1072 VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01
1073 VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01
1074 // check return value
1075 VERIFY_IS_EQUAL(info, 1);
1076 VERIFY(lm.nfev() < 510 ); // 490
1077 VERIFY(lm.njev() < 400 ); // 376
1078
1079 /*
1080 * Second try
1081 */
1082 x<< 0.25, 0.39, 0.415, 0.39;
1083 // do the computation
1084 lm.resetParameters();
1085 info = lm.minimize(x);
1086
1087 // check return value
1088 VERIFY_IS_EQUAL(info, 1);
1089 VERIFY_IS_EQUAL(lm.nfev(), 18);
1090 VERIFY_IS_EQUAL(lm.njev(), 16);
1091 // check norm^2
1092 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04);
1093 // check x
1094 VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01
1095 VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01
1096 VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01
1097 VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01
1098 }
1099
1100
1101
1102 struct Bennett5_functor : DenseFunctor<double>
1103 {
Bennett5_functorBennett5_functor1104 Bennett5_functor(void) : DenseFunctor<double>(3,154) {}
1105 static const double x[154];
1106 static const double y[154];
operator ()Bennett5_functor1107 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 }
dfBennett5_functor1115 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
testNistBennett5(void)1136 void testNistBennett5(void)
1137 {
1138 const int n=3;
1139 int info;
1140
1141 VectorXd x(n);
1142
1143 /*
1144 * First try
1145 */
1146 x<< -2000., 50., 0.8;
1147 // do the computation
1148 Bennett5_functor functor;
1149 LevenbergMarquardt<Bennett5_functor> lm(functor);
1150 lm.setMaxfev(1000);
1151 info = lm.minimize(x);
1152
1153 // check return value
1154 VERIFY_IS_EQUAL(info, 1);
1155 VERIFY_IS_EQUAL(lm.nfev(), 758);
1156 VERIFY_IS_EQUAL(lm.njev(), 744);
1157 // check norm^2
1158 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04);
1159 // check x
1160 VERIFY_IS_APPROX(x[0], -2.5235058043E+03);
1161 VERIFY_IS_APPROX(x[1], 4.6736564644E+01);
1162 VERIFY_IS_APPROX(x[2], 9.3218483193E-01);
1163 /*
1164 * Second try
1165 */
1166 x<< -1500., 45., 0.85;
1167 // do the computation
1168 lm.resetParameters();
1169 info = lm.minimize(x);
1170
1171 // check return value
1172 VERIFY_IS_EQUAL(info, 1);
1173 VERIFY_IS_EQUAL(lm.nfev(), 203);
1174 VERIFY_IS_EQUAL(lm.njev(), 192);
1175 // check norm^2
1176 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04);
1177 // check x
1178 VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03
1179 VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01);
1180 VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01);
1181 }
1182
1183 struct thurber_functor : DenseFunctor<double>
1184 {
thurber_functorthurber_functor1185 thurber_functor(void) : DenseFunctor<double>(7,37) {}
1186 static const double _x[37];
1187 static const double _y[37];
operator ()thurber_functor1188 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 }
dfthurber_functor1199 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
testNistThurber(void)1223 void testNistThurber(void)
1224 {
1225 const int n=7;
1226 int info;
1227
1228 VectorXd x(n);
1229
1230 /*
1231 * First try
1232 */
1233 x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ;
1234 // do the computation
1235 thurber_functor functor;
1236 LevenbergMarquardt<thurber_functor> lm(functor);
1237 lm.setFtol(1.E4*NumTraits<double>::epsilon());
1238 lm.setXtol(1.E4*NumTraits<double>::epsilon());
1239 info = lm.minimize(x);
1240
1241 // check return value
1242 VERIFY_IS_EQUAL(info, 1);
1243 VERIFY_IS_EQUAL(lm.nfev(), 39);
1244 VERIFY_IS_EQUAL(lm.njev(), 36);
1245 // check norm^2
1246 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03);
1247 // check x
1248 VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
1249 VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
1250 VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
1251 VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
1252 VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
1253 VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
1254 VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
1255
1256 /*
1257 * Second try
1258 */
1259 x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ;
1260 // do the computation
1261 lm.resetParameters();
1262 lm.setFtol(1.E4*NumTraits<double>::epsilon());
1263 lm.setXtol(1.E4*NumTraits<double>::epsilon());
1264 info = lm.minimize(x);
1265
1266 // check return value
1267 VERIFY_IS_EQUAL(info, 1);
1268 VERIFY_IS_EQUAL(lm.nfev(), 29);
1269 VERIFY_IS_EQUAL(lm.njev(), 28);
1270 // check norm^2
1271 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03);
1272 // check x
1273 VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
1274 VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
1275 VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
1276 VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
1277 VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
1278 VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
1279 VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
1280 }
1281
1282 struct rat43_functor : DenseFunctor<double>
1283 {
rat43_functorrat43_functor1284 rat43_functor(void) : DenseFunctor<double>(4,15) {}
1285 static const double x[15];
1286 static const double y[15];
operator ()rat43_functor1287 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 }
dfrat43_functor1295 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
testNistRat43(void)1315 void testNistRat43(void)
1316 {
1317 const int n=4;
1318 int info;
1319
1320 VectorXd x(n);
1321
1322 /*
1323 * First try
1324 */
1325 x<< 100., 10., 1., 1.;
1326 // do the computation
1327 rat43_functor functor;
1328 LevenbergMarquardt<rat43_functor> lm(functor);
1329 lm.setFtol(1.E6*NumTraits<double>::epsilon());
1330 lm.setXtol(1.E6*NumTraits<double>::epsilon());
1331 info = lm.minimize(x);
1332
1333 // check return value
1334 VERIFY_IS_EQUAL(info, 1);
1335 VERIFY_IS_EQUAL(lm.nfev(), 27);
1336 VERIFY_IS_EQUAL(lm.njev(), 20);
1337 // check norm^2
1338 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03);
1339 // check x
1340 VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
1341 VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
1342 VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
1343 VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
1344
1345 /*
1346 * Second try
1347 */
1348 x<< 700., 5., 0.75, 1.3;
1349 // do the computation
1350 lm.resetParameters();
1351 lm.setFtol(1.E5*NumTraits<double>::epsilon());
1352 lm.setXtol(1.E5*NumTraits<double>::epsilon());
1353 info = lm.minimize(x);
1354
1355 // check return value
1356 VERIFY_IS_EQUAL(info, 1);
1357 VERIFY_IS_EQUAL(lm.nfev(), 9);
1358 VERIFY_IS_EQUAL(lm.njev(), 8);
1359 // check norm^2
1360 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03);
1361 // check x
1362 VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
1363 VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
1364 VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
1365 VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
1366 }
1367
1368
1369
1370 struct eckerle4_functor : DenseFunctor<double>
1371 {
eckerle4_functoreckerle4_functor1372 eckerle4_functor(void) : DenseFunctor<double>(3,35) {}
1373 static const double x[35];
1374 static const double y[35];
operator ()eckerle4_functor1375 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 }
dfeckerle4_functor1383 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
testNistEckerle4(void)1402 void testNistEckerle4(void)
1403 {
1404 const int n=3;
1405 int info;
1406
1407 VectorXd x(n);
1408
1409 /*
1410 * First try
1411 */
1412 x<< 1., 10., 500.;
1413 // do the computation
1414 eckerle4_functor functor;
1415 LevenbergMarquardt<eckerle4_functor> lm(functor);
1416 info = lm.minimize(x);
1417
1418 // check return value
1419 VERIFY_IS_EQUAL(info, 1);
1420 VERIFY_IS_EQUAL(lm.nfev(), 18);
1421 VERIFY_IS_EQUAL(lm.njev(), 15);
1422 // check norm^2
1423 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03);
1424 // check x
1425 VERIFY_IS_APPROX(x[0], 1.5543827178);
1426 VERIFY_IS_APPROX(x[1], 4.0888321754);
1427 VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
1428
1429 /*
1430 * Second try
1431 */
1432 x<< 1.5, 5., 450.;
1433 // do the computation
1434 info = lm.minimize(x);
1435
1436 // check return value
1437 VERIFY_IS_EQUAL(info, 1);
1438 VERIFY_IS_EQUAL(lm.nfev(), 7);
1439 VERIFY_IS_EQUAL(lm.njev(), 6);
1440 // check norm^2
1441 VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03);
1442 // check x
1443 VERIFY_IS_APPROX(x[0], 1.5543827178);
1444 VERIFY_IS_APPROX(x[1], 4.0888321754);
1445 VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
1446 }
1447
test_levenberg_marquardt()1448 void test_levenberg_marquardt()
1449 {
1450 // Tests using the examples provided by (c)minpack
1451 CALL_SUBTEST(testLmder1());
1452 CALL_SUBTEST(testLmder());
1453 CALL_SUBTEST(testLmdif1());
1454 // CALL_SUBTEST(testLmstr1());
1455 // CALL_SUBTEST(testLmstr());
1456 CALL_SUBTEST(testLmdif());
1457
1458 // NIST tests, level of difficulty = "Lower"
1459 CALL_SUBTEST(testNistMisra1a());
1460 CALL_SUBTEST(testNistChwirut2());
1461
1462 // NIST tests, level of difficulty = "Average"
1463 CALL_SUBTEST(testNistHahn1());
1464 CALL_SUBTEST(testNistMisra1d());
1465 CALL_SUBTEST(testNistMGH17());
1466 CALL_SUBTEST(testNistLanczos1());
1467
1468 // // NIST tests, level of difficulty = "Higher"
1469 CALL_SUBTEST(testNistRat42());
1470 CALL_SUBTEST(testNistMGH10());
1471 CALL_SUBTEST(testNistBoxBOD());
1472 // CALL_SUBTEST(testNistMGH09());
1473 CALL_SUBTEST(testNistBennett5());
1474 CALL_SUBTEST(testNistThurber());
1475 CALL_SUBTEST(testNistRat43());
1476 CALL_SUBTEST(testNistEckerle4());
1477 }
1478