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
6#include <stdio.h>
7
8#include "main.h"
9#include <unsupported/Eigen/NonLinearOptimization>
10
11// This disables some useless Warnings on MSVC.
12// It is intended to be done for this test only.
13#include <Eigen/src/Core/util/DisableStupidWarnings.h>
14
15int fcn_chkder(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag)
16{
17    /*      subroutine fcn for chkder example. */
18
19    int i;
20    assert(15 ==  fvec.size());
21    assert(3 ==  x.size());
22    double tmp1, tmp2, tmp3, tmp4;
23    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,
24        3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
25
26
27    if (iflag == 0)
28        return 0;
29
30    if (iflag != 2)
31        for (i=0; i<15; i++) {
32            tmp1 = i+1;
33            tmp2 = 16-i-1;
34            tmp3 = tmp1;
35            if (i >= 8) tmp3 = tmp2;
36            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
37        }
38    else {
39        for (i = 0; i < 15; i++) {
40            tmp1 = i+1;
41            tmp2 = 16-i-1;
42
43            /* error introduced into next statement for illustration. */
44            /* corrected statement should read    tmp3 = tmp1 . */
45
46            tmp3 = tmp2;
47            if (i >= 8) tmp3 = tmp2;
48            tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4=tmp4*tmp4;
49            fjac(i,0) = -1.;
50            fjac(i,1) = tmp1*tmp2/tmp4;
51            fjac(i,2) = tmp1*tmp3/tmp4;
52        }
53    }
54    return 0;
55}
56
57
58void testChkder()
59{
60  const int m=15, n=3;
61  VectorXd x(n), fvec(m), xp, fvecp(m), err;
62  MatrixXd fjac(m,n);
63  VectorXi ipvt;
64
65  /*      the following values should be suitable for */
66  /*      checking the jacobian matrix. */
67  x << 9.2e-1, 1.3e-1, 5.4e-1;
68
69  internal::chkder(x, fvec, fjac, xp, fvecp, 1, err);
70  fcn_chkder(x, fvec, fjac, 1);
71  fcn_chkder(x, fvec, fjac, 2);
72  fcn_chkder(xp, fvecp, fjac, 1);
73  internal::chkder(x, fvec, fjac, xp, fvecp, 2, err);
74
75  fvecp -= fvec;
76
77  // check those
78  VectorXd fvec_ref(m), fvecp_ref(m), err_ref(m);
79  fvec_ref <<
80      -1.181606, -1.429655, -1.606344,
81      -1.745269, -1.840654, -1.921586,
82      -1.984141, -2.022537, -2.468977,
83      -2.827562, -3.473582, -4.437612,
84      -6.047662, -9.267761, -18.91806;
85  fvecp_ref <<
86      -7.724666e-09, -3.432406e-09, -2.034843e-10,
87      2.313685e-09,  4.331078e-09,  5.984096e-09,
88      7.363281e-09,   8.53147e-09,  1.488591e-08,
89      2.33585e-08,  3.522012e-08,  5.301255e-08,
90      8.26666e-08,  1.419747e-07,   3.19899e-07;
91  err_ref <<
92      0.1141397,  0.09943516,  0.09674474,
93      0.09980447,  0.1073116, 0.1220445,
94      0.1526814, 1, 1,
95      1, 1, 1,
96      1, 1, 1;
97
98  VERIFY_IS_APPROX(fvec, fvec_ref);
99  VERIFY_IS_APPROX(fvecp, fvecp_ref);
100  VERIFY_IS_APPROX(err, err_ref);
101}
102
103// Generic functor
104template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
105struct Functor
106{
107  typedef _Scalar Scalar;
108  enum {
109    InputsAtCompileTime = NX,
110    ValuesAtCompileTime = NY
111  };
112  typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
113  typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
114  typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
115
116  const int m_inputs, m_values;
117
118  Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
119  Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
120
121  int inputs() const { return m_inputs; }
122  int values() const { return m_values; }
123
124  // you should define that in the subclass :
125//  void operator() (const InputType& x, ValueType* v, JacobianType* _j=0) const;
126};
127
128struct lmder_functor : Functor<double>
129{
130    lmder_functor(void): Functor<double>(3,15) {}
131    int operator()(const VectorXd &x, VectorXd &fvec) const
132    {
133        double tmp1, tmp2, tmp3;
134        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,
135            3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
136
137        for (int i = 0; i < values(); i++)
138        {
139            tmp1 = i+1;
140            tmp2 = 16 - i - 1;
141            tmp3 = (i>=8)? tmp2 : tmp1;
142            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
143        }
144        return 0;
145    }
146
147    int df(const VectorXd &x, MatrixXd &fjac) const
148    {
149        double tmp1, tmp2, tmp3, tmp4;
150        for (int i = 0; i < values(); i++)
151        {
152            tmp1 = i+1;
153            tmp2 = 16 - i - 1;
154            tmp3 = (i>=8)? tmp2 : tmp1;
155            tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
156            fjac(i,0) = -1;
157            fjac(i,1) = tmp1*tmp2/tmp4;
158            fjac(i,2) = tmp1*tmp3/tmp4;
159        }
160        return 0;
161    }
162};
163
164void testLmder1()
165{
166  int n=3, info;
167
168  VectorXd x;
169
170  /* the following starting values provide a rough fit. */
171  x.setConstant(n, 1.);
172
173  // do the computation
174  lmder_functor functor;
175  LevenbergMarquardt<lmder_functor> lm(functor);
176  info = lm.lmder1(x);
177
178  // check return value
179  VERIFY_IS_EQUAL(info, 1);
180  VERIFY_IS_EQUAL(lm.nfev, 6);
181  VERIFY_IS_EQUAL(lm.njev, 5);
182
183  // check norm
184  VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596);
185
186  // check x
187  VectorXd x_ref(n);
188  x_ref << 0.08241058, 1.133037, 2.343695;
189  VERIFY_IS_APPROX(x, x_ref);
190}
191
192void testLmder()
193{
194  const int m=15, n=3;
195  int info;
196  double fnorm, covfac;
197  VectorXd x;
198
199  /* the following starting values provide a rough fit. */
200  x.setConstant(n, 1.);
201
202  // do the computation
203  lmder_functor functor;
204  LevenbergMarquardt<lmder_functor> lm(functor);
205  info = lm.minimize(x);
206
207  // check return values
208  VERIFY_IS_EQUAL(info, 1);
209  VERIFY_IS_EQUAL(lm.nfev, 6);
210  VERIFY_IS_EQUAL(lm.njev, 5);
211
212  // check norm
213  fnorm = lm.fvec.blueNorm();
214  VERIFY_IS_APPROX(fnorm, 0.09063596);
215
216  // check x
217  VectorXd x_ref(n);
218  x_ref << 0.08241058, 1.133037, 2.343695;
219  VERIFY_IS_APPROX(x, x_ref);
220
221  // check covariance
222  covfac = fnorm*fnorm/(m-n);
223  internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm
224
225  MatrixXd cov_ref(n,n);
226  cov_ref <<
227      0.0001531202,   0.002869941,  -0.002656662,
228      0.002869941,    0.09480935,   -0.09098995,
229      -0.002656662,   -0.09098995,    0.08778727;
230
231//  std::cout << fjac*covfac << std::endl;
232
233  MatrixXd cov;
234  cov =  covfac*lm.fjac.topLeftCorner<n,n>();
235  VERIFY_IS_APPROX( cov, cov_ref);
236  // TODO: why isn't this allowed ? :
237  // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
238}
239
240struct hybrj_functor : Functor<double>
241{
242    hybrj_functor(void) : Functor<double>(9,9) {}
243
244    int operator()(const VectorXd &x, VectorXd &fvec)
245    {
246        double temp, temp1, temp2;
247        const int n = x.size();
248        assert(fvec.size()==n);
249        for (int k = 0; k < n; k++)
250        {
251            temp = (3. - 2.*x[k])*x[k];
252            temp1 = 0.;
253            if (k) temp1 = x[k-1];
254            temp2 = 0.;
255            if (k != n-1) temp2 = x[k+1];
256            fvec[k] = temp - temp1 - 2.*temp2 + 1.;
257        }
258        return 0;
259    }
260    int df(const VectorXd &x, MatrixXd &fjac)
261    {
262        const int n = x.size();
263        assert(fjac.rows()==n);
264        assert(fjac.cols()==n);
265        for (int k = 0; k < n; k++)
266        {
267            for (int j = 0; j < n; j++)
268                fjac(k,j) = 0.;
269            fjac(k,k) = 3.- 4.*x[k];
270            if (k) fjac(k,k-1) = -1.;
271            if (k != n-1) fjac(k,k+1) = -2.;
272        }
273        return 0;
274    }
275};
276
277
278void testHybrj1()
279{
280  const int n=9;
281  int info;
282  VectorXd x(n);
283
284  /* the following starting values provide a rough fit. */
285  x.setConstant(n, -1.);
286
287  // do the computation
288  hybrj_functor functor;
289  HybridNonLinearSolver<hybrj_functor> solver(functor);
290  info = solver.hybrj1(x);
291
292  // check return value
293  VERIFY_IS_EQUAL(info, 1);
294  VERIFY_IS_EQUAL(solver.nfev, 11);
295  VERIFY_IS_EQUAL(solver.njev, 1);
296
297  // check norm
298  VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
299
300
301// check x
302  VectorXd x_ref(n);
303  x_ref <<
304     -0.5706545,    -0.6816283,    -0.7017325,
305     -0.7042129,     -0.701369,    -0.6918656,
306     -0.665792,    -0.5960342,    -0.4164121;
307  VERIFY_IS_APPROX(x, x_ref);
308}
309
310void testHybrj()
311{
312  const int n=9;
313  int info;
314  VectorXd x(n);
315
316  /* the following starting values provide a rough fit. */
317  x.setConstant(n, -1.);
318
319
320  // do the computation
321  hybrj_functor functor;
322  HybridNonLinearSolver<hybrj_functor> solver(functor);
323  solver.diag.setConstant(n, 1.);
324  solver.useExternalScaling = true;
325  info = solver.solve(x);
326
327  // check return value
328  VERIFY_IS_EQUAL(info, 1);
329  VERIFY_IS_EQUAL(solver.nfev, 11);
330  VERIFY_IS_EQUAL(solver.njev, 1);
331
332  // check norm
333  VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
334
335
336// check x
337  VectorXd x_ref(n);
338  x_ref <<
339     -0.5706545,    -0.6816283,    -0.7017325,
340     -0.7042129,     -0.701369,    -0.6918656,
341     -0.665792,    -0.5960342,    -0.4164121;
342  VERIFY_IS_APPROX(x, x_ref);
343
344}
345
346struct hybrd_functor : Functor<double>
347{
348    hybrd_functor(void) : Functor<double>(9,9) {}
349    int operator()(const VectorXd &x, VectorXd &fvec) const
350    {
351        double temp, temp1, temp2;
352        const int n = x.size();
353
354        assert(fvec.size()==n);
355        for (int k=0; k < n; k++)
356        {
357            temp = (3. - 2.*x[k])*x[k];
358            temp1 = 0.;
359            if (k) temp1 = x[k-1];
360            temp2 = 0.;
361            if (k != n-1) temp2 = x[k+1];
362            fvec[k] = temp - temp1 - 2.*temp2 + 1.;
363        }
364        return 0;
365    }
366};
367
368void testHybrd1()
369{
370  int n=9, info;
371  VectorXd x(n);
372
373  /* the following starting values provide a rough solution. */
374  x.setConstant(n, -1.);
375
376  // do the computation
377  hybrd_functor functor;
378  HybridNonLinearSolver<hybrd_functor> solver(functor);
379  info = solver.hybrd1(x);
380
381  // check return value
382  VERIFY_IS_EQUAL(info, 1);
383  VERIFY_IS_EQUAL(solver.nfev, 20);
384
385  // check norm
386  VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
387
388  // check x
389  VectorXd x_ref(n);
390  x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121;
391  VERIFY_IS_APPROX(x, x_ref);
392}
393
394void testHybrd()
395{
396  const int n=9;
397  int info;
398  VectorXd x;
399
400  /* the following starting values provide a rough fit. */
401  x.setConstant(n, -1.);
402
403  // do the computation
404  hybrd_functor functor;
405  HybridNonLinearSolver<hybrd_functor> solver(functor);
406  solver.parameters.nb_of_subdiagonals = 1;
407  solver.parameters.nb_of_superdiagonals = 1;
408  solver.diag.setConstant(n, 1.);
409  solver.useExternalScaling = true;
410  info = solver.solveNumericalDiff(x);
411
412  // check return value
413  VERIFY_IS_EQUAL(info, 1);
414  VERIFY_IS_EQUAL(solver.nfev, 14);
415
416  // check norm
417  VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
418
419  // check x
420  VectorXd x_ref(n);
421  x_ref <<
422      -0.5706545,    -0.6816283,    -0.7017325,
423      -0.7042129,     -0.701369,    -0.6918656,
424      -0.665792,    -0.5960342,    -0.4164121;
425  VERIFY_IS_APPROX(x, x_ref);
426}
427
428struct lmstr_functor : Functor<double>
429{
430    lmstr_functor(void) : Functor<double>(3,15) {}
431    int operator()(const VectorXd &x, VectorXd &fvec)
432    {
433        /*  subroutine fcn for lmstr1 example. */
434        double tmp1, tmp2, tmp3;
435        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,
436            3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
437
438        assert(15==fvec.size());
439        assert(3==x.size());
440
441        for (int i=0; i<15; i++)
442        {
443            tmp1 = i+1;
444            tmp2 = 16 - i - 1;
445            tmp3 = (i>=8)? tmp2 : tmp1;
446            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
447        }
448        return 0;
449    }
450    int df(const VectorXd &x, VectorXd &jac_row, VectorXd::Index rownb)
451    {
452        assert(x.size()==3);
453        assert(jac_row.size()==x.size());
454        double tmp1, tmp2, tmp3, tmp4;
455
456        int i = rownb-2;
457        tmp1 = i+1;
458        tmp2 = 16 - i - 1;
459        tmp3 = (i>=8)? tmp2 : tmp1;
460        tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
461        jac_row[0] = -1;
462        jac_row[1] = tmp1*tmp2/tmp4;
463        jac_row[2] = tmp1*tmp3/tmp4;
464        return 0;
465    }
466};
467
468void testLmstr1()
469{
470  const int n=3;
471  int info;
472
473  VectorXd x(n);
474
475  /* the following starting values provide a rough fit. */
476  x.setConstant(n, 1.);
477
478  // do the computation
479  lmstr_functor functor;
480  LevenbergMarquardt<lmstr_functor> lm(functor);
481  info = lm.lmstr1(x);
482
483  // check return value
484  VERIFY_IS_EQUAL(info, 1);
485  VERIFY_IS_EQUAL(lm.nfev, 6);
486  VERIFY_IS_EQUAL(lm.njev, 5);
487
488  // check norm
489  VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596);
490
491  // check x
492  VectorXd x_ref(n);
493  x_ref << 0.08241058, 1.133037, 2.343695 ;
494  VERIFY_IS_APPROX(x, x_ref);
495}
496
497void testLmstr()
498{
499  const int n=3;
500  int info;
501  double fnorm;
502  VectorXd x(n);
503
504  /* the following starting values provide a rough fit. */
505  x.setConstant(n, 1.);
506
507  // do the computation
508  lmstr_functor functor;
509  LevenbergMarquardt<lmstr_functor> lm(functor);
510  info = lm.minimizeOptimumStorage(x);
511
512  // check return values
513  VERIFY_IS_EQUAL(info, 1);
514  VERIFY_IS_EQUAL(lm.nfev, 6);
515  VERIFY_IS_EQUAL(lm.njev, 5);
516
517  // check norm
518  fnorm = lm.fvec.blueNorm();
519  VERIFY_IS_APPROX(fnorm, 0.09063596);
520
521  // check x
522  VectorXd x_ref(n);
523  x_ref << 0.08241058, 1.133037, 2.343695;
524  VERIFY_IS_APPROX(x, x_ref);
525
526}
527
528struct lmdif_functor : Functor<double>
529{
530    lmdif_functor(void) : Functor<double>(3,15) {}
531    int operator()(const VectorXd &x, VectorXd &fvec) const
532    {
533        int i;
534        double tmp1,tmp2,tmp3;
535        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,
536            3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0};
537
538        assert(x.size()==3);
539        assert(fvec.size()==15);
540        for (i=0; i<15; i++)
541        {
542            tmp1 = i+1;
543            tmp2 = 15 - i;
544            tmp3 = tmp1;
545
546            if (i >= 8) tmp3 = tmp2;
547            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
548        }
549        return 0;
550    }
551};
552
553void testLmdif1()
554{
555  const int n=3;
556  int info;
557
558  VectorXd x(n), fvec(15);
559
560  /* the following starting values provide a rough fit. */
561  x.setConstant(n, 1.);
562
563  // do the computation
564  lmdif_functor functor;
565  DenseIndex nfev;
566  info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev);
567
568  // check return value
569  VERIFY_IS_EQUAL(info, 1);
570  VERIFY_IS_EQUAL(nfev, 26);
571
572  // check norm
573  functor(x, fvec);
574  VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596);
575
576  // check x
577  VectorXd x_ref(n);
578  x_ref << 0.0824106, 1.1330366, 2.3436947;
579  VERIFY_IS_APPROX(x, x_ref);
580
581}
582
583void testLmdif()
584{
585  const int m=15, n=3;
586  int info;
587  double fnorm, covfac;
588  VectorXd x(n);
589
590  /* the following starting values provide a rough fit. */
591  x.setConstant(n, 1.);
592
593  // do the computation
594  lmdif_functor functor;
595  NumericalDiff<lmdif_functor> numDiff(functor);
596  LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff);
597  info = lm.minimize(x);
598
599  // check return values
600  VERIFY_IS_EQUAL(info, 1);
601  VERIFY_IS_EQUAL(lm.nfev, 26);
602
603  // check norm
604  fnorm = lm.fvec.blueNorm();
605  VERIFY_IS_APPROX(fnorm, 0.09063596);
606
607  // check x
608  VectorXd x_ref(n);
609  x_ref << 0.08241058, 1.133037, 2.343695;
610  VERIFY_IS_APPROX(x, x_ref);
611
612  // check covariance
613  covfac = fnorm*fnorm/(m-n);
614  internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm
615
616  MatrixXd cov_ref(n,n);
617  cov_ref <<
618      0.0001531202,   0.002869942,  -0.002656662,
619      0.002869942,    0.09480937,   -0.09098997,
620      -0.002656662,   -0.09098997,    0.08778729;
621
622//  std::cout << fjac*covfac << std::endl;
623
624  MatrixXd cov;
625  cov =  covfac*lm.fjac.topLeftCorner<n,n>();
626  VERIFY_IS_APPROX( cov, cov_ref);
627  // TODO: why isn't this allowed ? :
628  // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
629}
630
631struct chwirut2_functor : Functor<double>
632{
633    chwirut2_functor(void) : Functor<double>(3,54) {}
634    static const double m_x[54];
635    static const double m_y[54];
636    int operator()(const VectorXd &b, VectorXd &fvec)
637    {
638        int i;
639
640        assert(b.size()==3);
641        assert(fvec.size()==54);
642        for(i=0; i<54; i++) {
643            double x = m_x[i];
644            fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i];
645        }
646        return 0;
647    }
648    int df(const VectorXd &b, MatrixXd &fjac)
649    {
650        assert(b.size()==3);
651        assert(fjac.rows()==54);
652        assert(fjac.cols()==3);
653        for(int i=0; i<54; i++) {
654            double x = m_x[i];
655            double factor = 1./(b[1]+b[2]*x);
656            double e = exp(-b[0]*x);
657            fjac(i,0) = -x*e*factor;
658            fjac(i,1) = -e*factor*factor;
659            fjac(i,2) = -x*e*factor*factor;
660        }
661        return 0;
662    }
663};
664const 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};
665const 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  };
666
667// http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml
668void testNistChwirut2(void)
669{
670  const int n=3;
671  int info;
672
673  VectorXd x(n);
674
675  /*
676   * First try
677   */
678  x<< 0.1, 0.01, 0.02;
679  // do the computation
680  chwirut2_functor functor;
681  LevenbergMarquardt<chwirut2_functor> lm(functor);
682  info = lm.minimize(x);
683
684  // check return value
685  VERIFY_IS_EQUAL(info, 1);
686  VERIFY_IS_EQUAL(lm.nfev, 10);
687  VERIFY_IS_EQUAL(lm.njev, 8);
688  // check norm^2
689  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02);
690  // check x
691  VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
692  VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
693  VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
694
695  /*
696   * Second try
697   */
698  x<< 0.15, 0.008, 0.010;
699  // do the computation
700  lm.resetParameters();
701  lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
702  lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
703  info = lm.minimize(x);
704
705  // check return value
706  VERIFY_IS_EQUAL(info, 1);
707  VERIFY_IS_EQUAL(lm.nfev, 7);
708  VERIFY_IS_EQUAL(lm.njev, 6);
709  // check norm^2
710  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02);
711  // check x
712  VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
713  VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
714  VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
715}
716
717
718struct misra1a_functor : Functor<double>
719{
720    misra1a_functor(void) : Functor<double>(2,14) {}
721    static const double m_x[14];
722    static const double m_y[14];
723    int operator()(const VectorXd &b, VectorXd &fvec)
724    {
725        assert(b.size()==2);
726        assert(fvec.size()==14);
727        for(int i=0; i<14; i++) {
728            fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ;
729        }
730        return 0;
731    }
732    int df(const VectorXd &b, MatrixXd &fjac)
733    {
734        assert(b.size()==2);
735        assert(fjac.rows()==14);
736        assert(fjac.cols()==2);
737        for(int i=0; i<14; i++) {
738            fjac(i,0) = (1.-exp(-b[1]*m_x[i]));
739            fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i]));
740        }
741        return 0;
742    }
743};
744const 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};
745const 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};
746
747// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
748void testNistMisra1a(void)
749{
750  const int n=2;
751  int info;
752
753  VectorXd x(n);
754
755  /*
756   * First try
757   */
758  x<< 500., 0.0001;
759  // do the computation
760  misra1a_functor functor;
761  LevenbergMarquardt<misra1a_functor> lm(functor);
762  info = lm.minimize(x);
763
764  // check return value
765  VERIFY_IS_EQUAL(info, 1);
766  VERIFY_IS_EQUAL(lm.nfev, 19);
767  VERIFY_IS_EQUAL(lm.njev, 15);
768  // check norm^2
769  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01);
770  // check x
771  VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
772  VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
773
774  /*
775   * Second try
776   */
777  x<< 250., 0.0005;
778  // do the computation
779  info = lm.minimize(x);
780
781  // check return value
782  VERIFY_IS_EQUAL(info, 1);
783  VERIFY_IS_EQUAL(lm.nfev, 5);
784  VERIFY_IS_EQUAL(lm.njev, 4);
785  // check norm^2
786  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01);
787  // check x
788  VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
789  VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
790}
791
792struct hahn1_functor : Functor<double>
793{
794    hahn1_functor(void) : Functor<double>(7,236) {}
795    static const double m_x[236];
796    int operator()(const VectorXd &b, VectorXd &fvec)
797    {
798        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 , 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 , 12.596E0 , 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 };
799
800        //        int called=0; printf("call hahn1_functor with  iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
801
802        assert(b.size()==7);
803        assert(fvec.size()==236);
804        for(int i=0; i<236; i++) {
805            double x=m_x[i], xx=x*x, xxx=xx*x;
806            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];
807        }
808        return 0;
809    }
810
811    int df(const VectorXd &b, MatrixXd &fjac)
812    {
813        assert(b.size()==7);
814        assert(fjac.rows()==236);
815        assert(fjac.cols()==7);
816        for(int i=0; i<236; i++) {
817            double x=m_x[i], xx=x*x, xxx=xx*x;
818            double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
819            fjac(i,0) = 1.*fact;
820            fjac(i,1) = x*fact;
821            fjac(i,2) = xx*fact;
822            fjac(i,3) = xxx*fact;
823            fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
824            fjac(i,4) = x*fact;
825            fjac(i,5) = xx*fact;
826            fjac(i,6) = xxx*fact;
827        }
828        return 0;
829    }
830};
831const 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 , 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 , 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};
832
833// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml
834void testNistHahn1(void)
835{
836  const int  n=7;
837  int info;
838
839  VectorXd x(n);
840
841  /*
842   * First try
843   */
844  x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
845  // do the computation
846  hahn1_functor functor;
847  LevenbergMarquardt<hahn1_functor> lm(functor);
848  info = lm.minimize(x);
849
850  // check return value
851  VERIFY_IS_EQUAL(info, 1);
852  VERIFY_IS_EQUAL(lm.nfev, 11);
853  VERIFY_IS_EQUAL(lm.njev, 10);
854  // check norm^2
855  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00);
856  // check x
857  VERIFY_IS_APPROX(x[0], 1.0776351733E+00);
858  VERIFY_IS_APPROX(x[1],-1.2269296921E-01);
859  VERIFY_IS_APPROX(x[2], 4.0863750610E-03);
860  VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06
861  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
862  VERIFY_IS_APPROX(x[5], 2.4053735503E-04);
863  VERIFY_IS_APPROX(x[6],-1.2314450199E-07);
864
865  /*
866   * Second try
867   */
868  x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001;
869  // do the computation
870  info = lm.minimize(x);
871
872  // check return value
873  VERIFY_IS_EQUAL(info, 1);
874  VERIFY_IS_EQUAL(lm.nfev, 11);
875  VERIFY_IS_EQUAL(lm.njev, 10);
876  // check norm^2
877  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00);
878  // check x
879  VERIFY_IS_APPROX(x[0], 1.077640); // should be :  1.0776351733E+00
880  VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01
881  VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03
882  VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06
883  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
884  VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04
885  VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07
886
887}
888
889struct misra1d_functor : Functor<double>
890{
891    misra1d_functor(void) : Functor<double>(2,14) {}
892    static const double x[14];
893    static const double y[14];
894    int operator()(const VectorXd &b, VectorXd &fvec)
895    {
896        assert(b.size()==2);
897        assert(fvec.size()==14);
898        for(int i=0; i<14; i++) {
899            fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i];
900        }
901        return 0;
902    }
903    int df(const VectorXd &b, MatrixXd &fjac)
904    {
905        assert(b.size()==2);
906        assert(fjac.rows()==14);
907        assert(fjac.cols()==2);
908        for(int i=0; i<14; i++) {
909            double den = 1.+b[1]*x[i];
910            fjac(i,0) = b[1]*x[i] / den;
911            fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den;
912        }
913        return 0;
914    }
915};
916const 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};
917const 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};
918
919// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml
920void testNistMisra1d(void)
921{
922  const int n=2;
923  int info;
924
925  VectorXd x(n);
926
927  /*
928   * First try
929   */
930  x<< 500., 0.0001;
931  // do the computation
932  misra1d_functor functor;
933  LevenbergMarquardt<misra1d_functor> lm(functor);
934  info = lm.minimize(x);
935
936  // check return value
937  VERIFY_IS_EQUAL(info, 3);
938  VERIFY_IS_EQUAL(lm.nfev, 9);
939  VERIFY_IS_EQUAL(lm.njev, 7);
940  // check norm^2
941  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02);
942  // check x
943  VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
944  VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
945
946  /*
947   * Second try
948   */
949  x<< 450., 0.0003;
950  // do the computation
951  info = lm.minimize(x);
952
953  // check return value
954  VERIFY_IS_EQUAL(info, 1);
955  VERIFY_IS_EQUAL(lm.nfev, 4);
956  VERIFY_IS_EQUAL(lm.njev, 3);
957  // check norm^2
958  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02);
959  // check x
960  VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
961  VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
962}
963
964
965struct lanczos1_functor : Functor<double>
966{
967    lanczos1_functor(void) : Functor<double>(6,24) {}
968    static const double x[24];
969    static const double y[24];
970    int operator()(const VectorXd &b, VectorXd &fvec)
971    {
972        assert(b.size()==6);
973        assert(fvec.size()==24);
974        for(int i=0; i<24; i++)
975            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];
976        return 0;
977    }
978    int df(const VectorXd &b, MatrixXd &fjac)
979    {
980        assert(b.size()==6);
981        assert(fjac.rows()==24);
982        assert(fjac.cols()==6);
983        for(int i=0; i<24; i++) {
984            fjac(i,0) = exp(-b[1]*x[i]);
985            fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]);
986            fjac(i,2) = exp(-b[3]*x[i]);
987            fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]);
988            fjac(i,4) = exp(-b[5]*x[i]);
989            fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]);
990        }
991        return 0;
992    }
993};
994const 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 };
995const 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 };
996
997// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml
998void testNistLanczos1(void)
999{
1000  const int n=6;
1001  int info;
1002
1003  VectorXd x(n);
1004
1005  /*
1006   * First try
1007   */
1008  x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6;
1009  // do the computation
1010  lanczos1_functor functor;
1011  LevenbergMarquardt<lanczos1_functor> lm(functor);
1012  info = lm.minimize(x);
1013
1014  // check return value
1015  VERIFY_IS_EQUAL(info, 2);
1016  VERIFY_IS_EQUAL(lm.nfev, 79);
1017  VERIFY_IS_EQUAL(lm.njev, 72);
1018  // check norm^2
1019  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430899764097e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
1020  // check x
1021  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
1022  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
1023  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
1024  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
1025  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
1026  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
1027
1028  /*
1029   * Second try
1030   */
1031  x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3;
1032  // do the computation
1033  info = lm.minimize(x);
1034
1035  // check return value
1036  VERIFY_IS_EQUAL(info, 2);
1037  VERIFY_IS_EQUAL(lm.nfev, 9);
1038  VERIFY_IS_EQUAL(lm.njev, 8);
1039  // check norm^2
1040  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.428595533845e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
1041  // check x
1042  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
1043  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
1044  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
1045  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
1046  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
1047  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
1048
1049}
1050
1051struct rat42_functor : Functor<double>
1052{
1053    rat42_functor(void) : Functor<double>(3,9) {}
1054    static const double x[9];
1055    static const double y[9];
1056    int operator()(const VectorXd &b, VectorXd &fvec)
1057    {
1058        assert(b.size()==3);
1059        assert(fvec.size()==9);
1060        for(int i=0; i<9; i++) {
1061            fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i];
1062        }
1063        return 0;
1064    }
1065
1066    int df(const VectorXd &b, MatrixXd &fjac)
1067    {
1068        assert(b.size()==3);
1069        assert(fjac.rows()==9);
1070        assert(fjac.cols()==3);
1071        for(int i=0; i<9; i++) {
1072            double e = exp(b[1]-b[2]*x[i]);
1073            fjac(i,0) = 1./(1.+e);
1074            fjac(i,1) = -b[0]*e/(1.+e)/(1.+e);
1075            fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e);
1076        }
1077        return 0;
1078    }
1079};
1080const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 };
1081const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 };
1082
1083// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml
1084void testNistRat42(void)
1085{
1086  const int n=3;
1087  int info;
1088
1089  VectorXd x(n);
1090
1091  /*
1092   * First try
1093   */
1094  x<< 100., 1., 0.1;
1095  // do the computation
1096  rat42_functor functor;
1097  LevenbergMarquardt<rat42_functor> lm(functor);
1098  info = lm.minimize(x);
1099
1100  // check return value
1101  VERIFY_IS_EQUAL(info, 1);
1102  VERIFY_IS_EQUAL(lm.nfev, 10);
1103  VERIFY_IS_EQUAL(lm.njev, 8);
1104  // check norm^2
1105  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00);
1106  // check x
1107  VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
1108  VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
1109  VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
1110
1111  /*
1112   * Second try
1113   */
1114  x<< 75., 2.5, 0.07;
1115  // do the computation
1116  info = lm.minimize(x);
1117
1118  // check return value
1119  VERIFY_IS_EQUAL(info, 1);
1120  VERIFY_IS_EQUAL(lm.nfev, 6);
1121  VERIFY_IS_EQUAL(lm.njev, 5);
1122  // check norm^2
1123  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00);
1124  // check x
1125  VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
1126  VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
1127  VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
1128}
1129
1130struct MGH10_functor : Functor<double>
1131{
1132    MGH10_functor(void) : Functor<double>(3,16) {}
1133    static const double x[16];
1134    static const double y[16];
1135    int operator()(const VectorXd &b, VectorXd &fvec)
1136    {
1137        assert(b.size()==3);
1138        assert(fvec.size()==16);
1139        for(int i=0; i<16; i++)
1140            fvec[i] =  b[0] * exp(b[1]/(x[i]+b[2])) - y[i];
1141        return 0;
1142    }
1143    int df(const VectorXd &b, MatrixXd &fjac)
1144    {
1145        assert(b.size()==3);
1146        assert(fjac.rows()==16);
1147        assert(fjac.cols()==3);
1148        for(int i=0; i<16; i++) {
1149            double factor = 1./(x[i]+b[2]);
1150            double e = exp(b[1]*factor);
1151            fjac(i,0) = e;
1152            fjac(i,1) = b[0]*factor*e;
1153            fjac(i,2) = -b[1]*b[0]*factor*factor*e;
1154        }
1155        return 0;
1156    }
1157};
1158const 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 };
1159const 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 };
1160
1161// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
1162void testNistMGH10(void)
1163{
1164  const int n=3;
1165  int info;
1166
1167  VectorXd x(n);
1168
1169  /*
1170   * First try
1171   */
1172  x<< 2., 400000., 25000.;
1173  // do the computation
1174  MGH10_functor functor;
1175  LevenbergMarquardt<MGH10_functor> lm(functor);
1176  info = lm.minimize(x);
1177
1178  // check return value
1179  VERIFY_IS_EQUAL(info, 2);
1180  VERIFY_IS_EQUAL(lm.nfev, 284 );
1181  VERIFY_IS_EQUAL(lm.njev, 249 );
1182  // check norm^2
1183  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01);
1184  // check x
1185  VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
1186  VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
1187  VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
1188
1189  /*
1190   * Second try
1191   */
1192  x<< 0.02, 4000., 250.;
1193  // do the computation
1194  info = lm.minimize(x);
1195
1196  // check return value
1197  VERIFY_IS_EQUAL(info, 3);
1198  VERIFY_IS_EQUAL(lm.nfev, 126);
1199  VERIFY_IS_EQUAL(lm.njev, 116);
1200  // check norm^2
1201  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01);
1202  // check x
1203  VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
1204  VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
1205  VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
1206}
1207
1208
1209struct BoxBOD_functor : Functor<double>
1210{
1211    BoxBOD_functor(void) : Functor<double>(2,6) {}
1212    static const double x[6];
1213    int operator()(const VectorXd &b, VectorXd &fvec)
1214    {
1215        static const double y[6] = { 109., 149., 149., 191., 213., 224. };
1216        assert(b.size()==2);
1217        assert(fvec.size()==6);
1218        for(int i=0; i<6; i++)
1219            fvec[i] =  b[0]*(1.-exp(-b[1]*x[i])) - y[i];
1220        return 0;
1221    }
1222    int df(const VectorXd &b, MatrixXd &fjac)
1223    {
1224        assert(b.size()==2);
1225        assert(fjac.rows()==6);
1226        assert(fjac.cols()==2);
1227        for(int i=0; i<6; i++) {
1228            double e = exp(-b[1]*x[i]);
1229            fjac(i,0) = 1.-e;
1230            fjac(i,1) = b[0]*x[i]*e;
1231        }
1232        return 0;
1233    }
1234};
1235const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. };
1236
1237// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml
1238void testNistBoxBOD(void)
1239{
1240  const int n=2;
1241  int info;
1242
1243  VectorXd x(n);
1244
1245  /*
1246   * First try
1247   */
1248  x<< 1., 1.;
1249  // do the computation
1250  BoxBOD_functor functor;
1251  LevenbergMarquardt<BoxBOD_functor> lm(functor);
1252  lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
1253  lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
1254  lm.parameters.factor = 10.;
1255  info = lm.minimize(x);
1256
1257  // check return value
1258  VERIFY_IS_EQUAL(info, 1);
1259  VERIFY_IS_EQUAL(lm.nfev, 31);
1260  VERIFY_IS_EQUAL(lm.njev, 25);
1261  // check norm^2
1262  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
1263  // check x
1264  VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
1265  VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
1266
1267  /*
1268   * Second try
1269   */
1270  x<< 100., 0.75;
1271  // do the computation
1272  lm.resetParameters();
1273  lm.parameters.ftol = NumTraits<double>::epsilon();
1274  lm.parameters.xtol = NumTraits<double>::epsilon();
1275  info = lm.minimize(x);
1276
1277  // check return value
1278  VERIFY_IS_EQUAL(info, 1);
1279  VERIFY_IS_EQUAL(lm.nfev, 15 );
1280  VERIFY_IS_EQUAL(lm.njev, 14 );
1281  // check norm^2
1282  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
1283  // check x
1284  VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
1285  VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
1286}
1287
1288struct MGH17_functor : Functor<double>
1289{
1290    MGH17_functor(void) : Functor<double>(5,33) {}
1291    static const double x[33];
1292    static const double y[33];
1293    int operator()(const VectorXd &b, VectorXd &fvec)
1294    {
1295        assert(b.size()==5);
1296        assert(fvec.size()==33);
1297        for(int i=0; i<33; i++)
1298            fvec[i] =  b[0] + b[1]*exp(-b[3]*x[i]) +  b[2]*exp(-b[4]*x[i]) - y[i];
1299        return 0;
1300    }
1301    int df(const VectorXd &b, MatrixXd &fjac)
1302    {
1303        assert(b.size()==5);
1304        assert(fjac.rows()==33);
1305        assert(fjac.cols()==5);
1306        for(int i=0; i<33; i++) {
1307            fjac(i,0) = 1.;
1308            fjac(i,1) = exp(-b[3]*x[i]);
1309            fjac(i,2) = exp(-b[4]*x[i]);
1310            fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]);
1311            fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]);
1312        }
1313        return 0;
1314    }
1315};
1316const 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 };
1317const 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 };
1318
1319// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml
1320void testNistMGH17(void)
1321{
1322  const int n=5;
1323  int info;
1324
1325  VectorXd x(n);
1326
1327  /*
1328   * First try
1329   */
1330  x<< 50., 150., -100., 1., 2.;
1331  // do the computation
1332  MGH17_functor functor;
1333  LevenbergMarquardt<MGH17_functor> lm(functor);
1334  lm.parameters.ftol = NumTraits<double>::epsilon();
1335  lm.parameters.xtol = NumTraits<double>::epsilon();
1336  lm.parameters.maxfev = 1000;
1337  info = lm.minimize(x);
1338
1339  // check return value
1340  VERIFY_IS_EQUAL(info, 2);
1341  VERIFY_IS_EQUAL(lm.nfev, 602 );
1342  VERIFY_IS_EQUAL(lm.njev, 545 );
1343  // check norm^2
1344  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
1345  // check x
1346  VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
1347  VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
1348  VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
1349  VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
1350  VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
1351
1352  /*
1353   * Second try
1354   */
1355  x<< 0.5  ,1.5  ,-1   ,0.01 ,0.02;
1356  // do the computation
1357  lm.resetParameters();
1358  info = lm.minimize(x);
1359
1360  // check return value
1361  VERIFY_IS_EQUAL(info, 1);
1362  VERIFY_IS_EQUAL(lm.nfev, 18);
1363  VERIFY_IS_EQUAL(lm.njev, 15);
1364  // check norm^2
1365  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
1366  // check x
1367  VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
1368  VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
1369  VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
1370  VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
1371  VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
1372}
1373
1374struct MGH09_functor : Functor<double>
1375{
1376    MGH09_functor(void) : Functor<double>(4,11) {}
1377    static const double _x[11];
1378    static const double y[11];
1379    int operator()(const VectorXd &b, VectorXd &fvec)
1380    {
1381        assert(b.size()==4);
1382        assert(fvec.size()==11);
1383        for(int i=0; i<11; i++) {
1384            double x = _x[i], xx=x*x;
1385            fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i];
1386        }
1387        return 0;
1388    }
1389    int df(const VectorXd &b, MatrixXd &fjac)
1390    {
1391        assert(b.size()==4);
1392        assert(fjac.rows()==11);
1393        assert(fjac.cols()==4);
1394        for(int i=0; i<11; i++) {
1395            double x = _x[i], xx=x*x;
1396            double factor = 1./(xx+x*b[2]+b[3]);
1397            fjac(i,0) = (xx+x*b[1]) * factor;
1398            fjac(i,1) = b[0]*x* factor;
1399            fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor;
1400            fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor;
1401        }
1402        return 0;
1403    }
1404};
1405const 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 };
1406const 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 };
1407
1408// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
1409void testNistMGH09(void)
1410{
1411  const int n=4;
1412  int info;
1413
1414  VectorXd x(n);
1415
1416  /*
1417   * First try
1418   */
1419  x<< 25., 39, 41.5, 39.;
1420  // do the computation
1421  MGH09_functor functor;
1422  LevenbergMarquardt<MGH09_functor> lm(functor);
1423  lm.parameters.maxfev = 1000;
1424  info = lm.minimize(x);
1425
1426  // check return value
1427  VERIFY_IS_EQUAL(info, 1);
1428  VERIFY_IS_EQUAL(lm.nfev, 490 );
1429  VERIFY_IS_EQUAL(lm.njev, 376 );
1430  // check norm^2
1431  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04);
1432  // check x
1433  VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01
1434  VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01
1435  VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01
1436  VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01
1437
1438  /*
1439   * Second try
1440   */
1441  x<< 0.25, 0.39, 0.415, 0.39;
1442  // do the computation
1443  lm.resetParameters();
1444  info = lm.minimize(x);
1445
1446  // check return value
1447  VERIFY_IS_EQUAL(info, 1);
1448  VERIFY_IS_EQUAL(lm.nfev, 18);
1449  VERIFY_IS_EQUAL(lm.njev, 16);
1450  // check norm^2
1451  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04);
1452  // check x
1453  VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01
1454  VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01
1455  VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01
1456  VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01
1457}
1458
1459
1460
1461struct Bennett5_functor : Functor<double>
1462{
1463    Bennett5_functor(void) : Functor<double>(3,154) {}
1464    static const double x[154];
1465    static const double y[154];
1466    int operator()(const VectorXd &b, VectorXd &fvec)
1467    {
1468        assert(b.size()==3);
1469        assert(fvec.size()==154);
1470        for(int i=0; i<154; i++)
1471            fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i];
1472        return 0;
1473    }
1474    int df(const VectorXd &b, MatrixXd &fjac)
1475    {
1476        assert(b.size()==3);
1477        assert(fjac.rows()==154);
1478        assert(fjac.cols()==3);
1479        for(int i=0; i<154; i++) {
1480            double e = pow(b[1]+x[i],-1./b[2]);
1481            fjac(i,0) = e;
1482            fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]);
1483            fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2];
1484        }
1485        return 0;
1486    }
1487};
1488const 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, 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 };
1489const 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 ,-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 ,-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 };
1490
1491// http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml
1492void testNistBennett5(void)
1493{
1494  const int  n=3;
1495  int info;
1496
1497  VectorXd x(n);
1498
1499  /*
1500   * First try
1501   */
1502  x<< -2000., 50., 0.8;
1503  // do the computation
1504  Bennett5_functor functor;
1505  LevenbergMarquardt<Bennett5_functor> lm(functor);
1506  lm.parameters.maxfev = 1000;
1507  info = lm.minimize(x);
1508
1509  // check return value
1510  VERIFY_IS_EQUAL(info, 1);
1511  VERIFY_IS_EQUAL(lm.nfev, 758);
1512  VERIFY_IS_EQUAL(lm.njev, 744);
1513  // check norm^2
1514  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04);
1515  // check x
1516  VERIFY_IS_APPROX(x[0], -2.5235058043E+03);
1517  VERIFY_IS_APPROX(x[1], 4.6736564644E+01);
1518  VERIFY_IS_APPROX(x[2], 9.3218483193E-01);
1519  /*
1520   * Second try
1521   */
1522  x<< -1500., 45., 0.85;
1523  // do the computation
1524  lm.resetParameters();
1525  info = lm.minimize(x);
1526
1527  // check return value
1528  VERIFY_IS_EQUAL(info, 1);
1529  VERIFY_IS_EQUAL(lm.nfev, 203);
1530  VERIFY_IS_EQUAL(lm.njev, 192);
1531  // check norm^2
1532  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04);
1533  // check x
1534  VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03
1535  VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01);
1536  VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01);
1537}
1538
1539struct thurber_functor : Functor<double>
1540{
1541    thurber_functor(void) : Functor<double>(7,37) {}
1542    static const double _x[37];
1543    static const double _y[37];
1544    int operator()(const VectorXd &b, VectorXd &fvec)
1545    {
1546        //        int called=0; printf("call hahn1_functor with  iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
1547        assert(b.size()==7);
1548        assert(fvec.size()==37);
1549        for(int i=0; i<37; i++) {
1550            double x=_x[i], xx=x*x, xxx=xx*x;
1551            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];
1552        }
1553        return 0;
1554    }
1555    int df(const VectorXd &b, MatrixXd &fjac)
1556    {
1557        assert(b.size()==7);
1558        assert(fjac.rows()==37);
1559        assert(fjac.cols()==7);
1560        for(int i=0; i<37; i++) {
1561            double x=_x[i], xx=x*x, xxx=xx*x;
1562            double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
1563            fjac(i,0) = 1.*fact;
1564            fjac(i,1) = x*fact;
1565            fjac(i,2) = xx*fact;
1566            fjac(i,3) = xxx*fact;
1567            fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
1568            fjac(i,4) = x*fact;
1569            fjac(i,5) = xx*fact;
1570            fjac(i,6) = xxx*fact;
1571        }
1572        return 0;
1573    }
1574};
1575const 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 };
1576const 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};
1577
1578// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml
1579void testNistThurber(void)
1580{
1581  const int n=7;
1582  int info;
1583
1584  VectorXd x(n);
1585
1586  /*
1587   * First try
1588   */
1589  x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ;
1590  // do the computation
1591  thurber_functor functor;
1592  LevenbergMarquardt<thurber_functor> lm(functor);
1593  lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon();
1594  lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon();
1595  info = lm.minimize(x);
1596
1597  // check return value
1598  VERIFY_IS_EQUAL(info, 1);
1599  VERIFY_IS_EQUAL(lm.nfev, 39);
1600  VERIFY_IS_EQUAL(lm.njev, 36);
1601  // check norm^2
1602  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03);
1603  // check x
1604  VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
1605  VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
1606  VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
1607  VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
1608  VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
1609  VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
1610  VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
1611
1612  /*
1613   * Second try
1614   */
1615  x<< 1300 ,1500 ,500  ,75   ,1    ,0.4  ,0.05  ;
1616  // do the computation
1617  lm.resetParameters();
1618  lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon();
1619  lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon();
1620  info = lm.minimize(x);
1621
1622  // check return value
1623  VERIFY_IS_EQUAL(info, 1);
1624  VERIFY_IS_EQUAL(lm.nfev, 29);
1625  VERIFY_IS_EQUAL(lm.njev, 28);
1626  // check norm^2
1627  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03);
1628  // check x
1629  VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
1630  VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
1631  VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
1632  VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
1633  VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
1634  VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
1635  VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
1636}
1637
1638struct rat43_functor : Functor<double>
1639{
1640    rat43_functor(void) : Functor<double>(4,15) {}
1641    static const double x[15];
1642    static const double y[15];
1643    int operator()(const VectorXd &b, VectorXd &fvec)
1644    {
1645        assert(b.size()==4);
1646        assert(fvec.size()==15);
1647        for(int i=0; i<15; i++)
1648            fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i];
1649        return 0;
1650    }
1651    int df(const VectorXd &b, MatrixXd &fjac)
1652    {
1653        assert(b.size()==4);
1654        assert(fjac.rows()==15);
1655        assert(fjac.cols()==4);
1656        for(int i=0; i<15; i++) {
1657            double e = exp(b[1]-b[2]*x[i]);
1658            double power = -1./b[3];
1659            fjac(i,0) = pow(1.+e, power);
1660            fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.);
1661            fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.);
1662            fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power);
1663        }
1664        return 0;
1665    }
1666};
1667const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. };
1668const 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 };
1669
1670// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml
1671void testNistRat43(void)
1672{
1673  const int n=4;
1674  int info;
1675
1676  VectorXd x(n);
1677
1678  /*
1679   * First try
1680   */
1681  x<< 100., 10., 1., 1.;
1682  // do the computation
1683  rat43_functor functor;
1684  LevenbergMarquardt<rat43_functor> lm(functor);
1685  lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
1686  lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
1687  info = lm.minimize(x);
1688
1689  // check return value
1690  VERIFY_IS_EQUAL(info, 1);
1691  VERIFY_IS_EQUAL(lm.nfev, 27);
1692  VERIFY_IS_EQUAL(lm.njev, 20);
1693  // check norm^2
1694  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03);
1695  // check x
1696  VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
1697  VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
1698  VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
1699  VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
1700
1701  /*
1702   * Second try
1703   */
1704  x<< 700., 5., 0.75, 1.3;
1705  // do the computation
1706  lm.resetParameters();
1707  lm.parameters.ftol = 1.E5*NumTraits<double>::epsilon();
1708  lm.parameters.xtol = 1.E5*NumTraits<double>::epsilon();
1709  info = lm.minimize(x);
1710
1711  // check return value
1712  VERIFY_IS_EQUAL(info, 1);
1713  VERIFY_IS_EQUAL(lm.nfev, 9);
1714  VERIFY_IS_EQUAL(lm.njev, 8);
1715  // check norm^2
1716  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03);
1717  // check x
1718  VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
1719  VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
1720  VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
1721  VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
1722}
1723
1724
1725
1726struct eckerle4_functor : Functor<double>
1727{
1728    eckerle4_functor(void) : Functor<double>(3,35) {}
1729    static const double x[35];
1730    static const double y[35];
1731    int operator()(const VectorXd &b, VectorXd &fvec)
1732    {
1733        assert(b.size()==3);
1734        assert(fvec.size()==35);
1735        for(int i=0; i<35; i++)
1736            fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i];
1737        return 0;
1738    }
1739    int df(const VectorXd &b, MatrixXd &fjac)
1740    {
1741        assert(b.size()==3);
1742        assert(fjac.rows()==35);
1743        assert(fjac.cols()==3);
1744        for(int i=0; i<35; i++) {
1745            double b12 = b[1]*b[1];
1746            double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12);
1747            fjac(i,0) = e / b[1];
1748            fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12;
1749            fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12;
1750        }
1751        return 0;
1752    }
1753};
1754const 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};
1755const 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 };
1756
1757// http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml
1758void testNistEckerle4(void)
1759{
1760  const int n=3;
1761  int info;
1762
1763  VectorXd x(n);
1764
1765  /*
1766   * First try
1767   */
1768  x<< 1., 10., 500.;
1769  // do the computation
1770  eckerle4_functor functor;
1771  LevenbergMarquardt<eckerle4_functor> lm(functor);
1772  info = lm.minimize(x);
1773
1774  // check return value
1775  VERIFY_IS_EQUAL(info, 1);
1776  VERIFY_IS_EQUAL(lm.nfev, 18);
1777  VERIFY_IS_EQUAL(lm.njev, 15);
1778  // check norm^2
1779  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03);
1780  // check x
1781  VERIFY_IS_APPROX(x[0], 1.5543827178);
1782  VERIFY_IS_APPROX(x[1], 4.0888321754);
1783  VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
1784
1785  /*
1786   * Second try
1787   */
1788  x<< 1.5, 5., 450.;
1789  // do the computation
1790  info = lm.minimize(x);
1791
1792  // check return value
1793  VERIFY_IS_EQUAL(info, 1);
1794  VERIFY_IS_EQUAL(lm.nfev, 7);
1795  VERIFY_IS_EQUAL(lm.njev, 6);
1796  // check norm^2
1797  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03);
1798  // check x
1799  VERIFY_IS_APPROX(x[0], 1.5543827178);
1800  VERIFY_IS_APPROX(x[1], 4.0888321754);
1801  VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
1802}
1803
1804void test_NonLinearOptimization()
1805{
1806    // Tests using the examples provided by (c)minpack
1807    CALL_SUBTEST/*_1*/(testChkder());
1808    CALL_SUBTEST/*_1*/(testLmder1());
1809    CALL_SUBTEST/*_1*/(testLmder());
1810    CALL_SUBTEST/*_2*/(testHybrj1());
1811    CALL_SUBTEST/*_2*/(testHybrj());
1812    CALL_SUBTEST/*_2*/(testHybrd1());
1813    CALL_SUBTEST/*_2*/(testHybrd());
1814    CALL_SUBTEST/*_3*/(testLmstr1());
1815    CALL_SUBTEST/*_3*/(testLmstr());
1816    CALL_SUBTEST/*_3*/(testLmdif1());
1817    CALL_SUBTEST/*_3*/(testLmdif());
1818
1819    // NIST tests, level of difficulty = "Lower"
1820    CALL_SUBTEST/*_4*/(testNistMisra1a());
1821    CALL_SUBTEST/*_4*/(testNistChwirut2());
1822
1823    // NIST tests, level of difficulty = "Average"
1824    CALL_SUBTEST/*_5*/(testNistHahn1());
1825    CALL_SUBTEST/*_6*/(testNistMisra1d());
1826    CALL_SUBTEST/*_7*/(testNistMGH17());
1827    CALL_SUBTEST/*_8*/(testNistLanczos1());
1828
1829    // NIST tests, level of difficulty = "Higher"
1830    CALL_SUBTEST/*_9*/(testNistRat42());
1831    CALL_SUBTEST/*_10*/(testNistMGH10());
1832    CALL_SUBTEST/*_11*/(testNistBoxBOD());
1833    CALL_SUBTEST/*_12*/(testNistMGH09());
1834    CALL_SUBTEST/*_13*/(testNistBennett5());
1835    CALL_SUBTEST/*_14*/(testNistThurber());
1836    CALL_SUBTEST/*_15*/(testNistRat43());
1837    CALL_SUBTEST/*_16*/(testNistEckerle4());
1838}
1839
1840/*
1841 * Can be useful for debugging...
1842  printf("info, nfev : %d, %d\n", info, lm.nfev);
1843  printf("info, nfev, njev : %d, %d, %d\n", info, solver.nfev, solver.njev);
1844  printf("info, nfev : %d, %d\n", info, solver.nfev);
1845  printf("x[0] : %.32g\n", x[0]);
1846  printf("x[1] : %.32g\n", x[1]);
1847  printf("x[2] : %.32g\n", x[2]);
1848  printf("x[3] : %.32g\n", x[3]);
1849  printf("fvec.blueNorm() : %.32g\n", solver.fvec.blueNorm());
1850  printf("fvec.blueNorm() : %.32g\n", lm.fvec.blueNorm());
1851
1852  printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev);
1853  printf("fvec.squaredNorm() : %.13g\n", lm.fvec.squaredNorm());
1854  std::cout << x << std::endl;
1855  std::cout.precision(9);
1856  std::cout << x[0] << std::endl;
1857  std::cout << x[1] << std::endl;
1858  std::cout << x[2] << std::endl;
1859  std::cout << x[3] << std::endl;
1860*/
1861
1862