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 | {
|
---|
25 | lmder_functor(void): DenseFunctor<double>(3,15) {}
|
---|
26 | 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 |
|
---|
42 | 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 |
|
---|
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 |
|
---|
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 | {
|
---|
137 | lmdif_functor(void) : DenseFunctor<double>(3,15) {}
|
---|
138 | 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 |
|
---|
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 |
|
---|
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 | {
|
---|
240 | chwirut2_functor(void) : DenseFunctor<double>(3,54) {}
|
---|
241 | static const double m_x[54];
|
---|
242 | static const double m_y[54];
|
---|
243 | 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 | }
|
---|
255 | 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
|
---|
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 | {
|
---|
327 | misra1a_functor(void) : DenseFunctor<double>(2,14) {}
|
---|
328 | static const double m_x[14];
|
---|
329 | static const double m_y[14];
|
---|
330 | 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 | }
|
---|
339 | 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
|
---|
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 | {
|
---|
401 | hahn1_functor(void) : DenseFunctor<double>(7,236) {}
|
---|
402 | static const double m_x[236];
|
---|
403 | 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 |
|
---|
421 | 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
|
---|
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 | {
|
---|
503 | misra1d_functor(void) : DenseFunctor<double>(2,14) {}
|
---|
504 | static const double x[14];
|
---|
505 | static const double y[14];
|
---|
506 | 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 | }
|
---|
515 | 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
|
---|
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 | {
|
---|
579 | lanczos1_functor(void) : DenseFunctor<double>(6,24) {}
|
---|
580 | static const double x[24];
|
---|
581 | static const double y[24];
|
---|
582 | 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 | }
|
---|
590 | 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
|
---|
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 | {
|
---|
665 | rat42_functor(void) : DenseFunctor<double>(3,9) {}
|
---|
666 | static const double x[9];
|
---|
667 | static const double y[9];
|
---|
668 | 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 |
|
---|
678 | 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
|
---|
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 | {
|
---|
744 | MGH10_functor(void) : DenseFunctor<double>(3,16) {}
|
---|
745 | static const double x[16];
|
---|
746 | static const double y[16];
|
---|
747 | 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 | }
|
---|
755 | 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
|
---|
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 | {
|
---|
823 | BoxBOD_functor(void) : DenseFunctor<double>(2,6) {}
|
---|
824 | static const double x[6];
|
---|
825 | 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 | }
|
---|
834 | 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
|
---|
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 | {
|
---|
902 | MGH17_functor(void) : DenseFunctor<double>(5,33) {}
|
---|
903 | static const double x[33];
|
---|
904 | static const double y[33];
|
---|
905 | 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 | }
|
---|
913 | 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
|
---|
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 | {
|
---|
988 | MGH09_functor(void) : DenseFunctor<double>(4,11) {}
|
---|
989 | static const double _x[11];
|
---|
990 | static const double y[11];
|
---|
991 | 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 | }
|
---|
1001 | 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
|
---|
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 | {
|
---|
1075 | Bennett5_functor(void) : DenseFunctor<double>(3,154) {}
|
---|
1076 | static const double x[154];
|
---|
1077 | static const double y[154];
|
---|
1078 | 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 | }
|
---|
1086 | 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
|
---|
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 | {
|
---|
1156 | thurber_functor(void) : DenseFunctor<double>(7,37) {}
|
---|
1157 | static const double _x[37];
|
---|
1158 | static const double _y[37];
|
---|
1159 | 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 | }
|
---|
1170 | 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
|
---|
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 | {
|
---|
1255 | rat43_functor(void) : DenseFunctor<double>(4,15) {}
|
---|
1256 | static const double x[15];
|
---|
1257 | static const double y[15];
|
---|
1258 | 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 | }
|
---|
1266 | 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
|
---|
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 | {
|
---|
1343 | eckerle4_functor(void) : DenseFunctor<double>(3,35) {}
|
---|
1344 | static const double x[35];
|
---|
1345 | static const double y[35];
|
---|
1346 | 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 | }
|
---|
1354 | 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
|
---|
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 |
|
---|
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 | }
|
---|