mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-21 07:19:46 +08:00
Code Issues Packages Projects Releases Wiki Activity
d8c671f475
eigen/unsupported/test/NonLinear.cpp
Thomas Capricelli d8c671f475 yet another (difficult) Nist test : Thurber
2009-08-14 17:46:28 +02:00

1785 lines
62 KiB
C++
Raw Blame History

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
#include <stdio.h>
#include "main.h"
#include <unsupported/Eigen/NonLinear>
int fcn_chkder(int /*m*/, int /*n*/, const double *x, double *fvec, double *fjac, int ldfjac, int iflag)
{
/* subroutine fcn for chkder example. */
int i;
double tmp1, tmp2, tmp3, tmp4;
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, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
if (iflag == 0)
{
/* insert print statements here when nprint is positive. */
return 0;
}
if (iflag != 2)
for (i=1; i<=15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3));
}
else
{
for (i = 1; i <= 15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
/* error introduced into next statement for illustration. */
/* corrected statement should read tmp3 = tmp1 . */
tmp3 = tmp2;
if (i > 8) tmp3 = tmp2;
tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4=tmp4*tmp4;
fjac[i-1+ ldfjac*(1-1)] = -1.;
fjac[i-1+ ldfjac*(2-1)] = tmp1*tmp2/tmp4;
fjac[i-1+ ldfjac*(3-1)] = tmp1*tmp3/tmp4;
}
}
return 0;
}
void testChkder()
{
int i, m, n, ldfjac;
double x[3], fvec[15], fjac[15*3], xp[3], fvecp[15],
err[15];
m = 15;
n = 3;
/* the following values should be suitable for */
/* checking the jacobian matrix. */
x[1-1] = 9.2e-1;
x[2-1] = 1.3e-1;
x[3-1] = 5.4e-1;
ldfjac = 15;
chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, 1, err);
fcn_chkder(m, n, x, fvec, fjac, ldfjac, 1);
fcn_chkder(m, n, x, fvec, fjac, ldfjac, 2);
fcn_chkder(m, n, xp, fvecp, fjac, ldfjac, 1);
chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, 2, err);
for (i=1; i<=m; i++)
{
fvecp[i-1] = fvecp[i-1] - fvec[i-1];
}
double fvec_ref[] = {
-1.181606, -1.429655, -1.606344,
-1.745269, -1.840654, -1.921586,
-1.984141, -2.022537, -2.468977,
-2.827562, -3.473582, -4.437612,
-6.047662, -9.267761, -18.91806
};
double fvecp_ref[] = {
-7.724666e-09, -3.432406e-09, -2.034843e-10,
2.313685e-09, 4.331078e-09, 5.984096e-09,
7.363281e-09, 8.53147e-09, 1.488591e-08,
2.33585e-08, 3.522012e-08, 5.301255e-08,
8.26666e-08, 1.419747e-07, 3.19899e-07
};
double err_ref[] = {
0.1141397, 0.09943516, 0.09674474,
0.09980447, 0.1073116, 0.1220445,
0.1526814, 1, 1,
1, 1, 1,
1, 1, 1
};
for (i=1; i<=m; i++) VERIFY_IS_APPROX(fvec[i-1], fvec_ref[i-1]);
for (i=1; i<=m; i++) VERIFY_IS_APPROX(fvecp[i-1], fvecp_ref[i-1]);
for (i=1; i<=m; i++) VERIFY_IS_APPROX(err[i-1], err_ref[i-1]);
}
struct lmder1_functor {
static int f(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjac, int ldfjac, int iflag)
{
/* subroutine fcn for lmder1 example. */
int i;
double tmp1, tmp2, tmp3, tmp4;
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, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
if (iflag != 2)
{
for (i = 1; i <= 15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3));
}
}
else
{
for ( i = 1; i <= 15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4 = tmp4*tmp4;
fjac[i-1 + ldfjac*(1-1)] = -1.;
fjac[i-1 + ldfjac*(2-1)] = tmp1*tmp2/tmp4;
fjac[i-1 + ldfjac*(3-1)] = tmp1*tmp3/tmp4;
}
}
return 0;
}
};
void testLmder1()
{
int m=15, n=3, info;
Eigen::VectorXd x(n), fvec(m);
VectorXi ipvt;
/* the following starting values provide a rough fit. */
x.setConstant(n, 1.);
// do the computation
info = ei_lmder1<lmder1_functor,double>(x, fvec, ipvt);
// check return value
VERIFY( 1 == info);
// check norm
VERIFY_IS_APPROX(fvec.norm(), 0.09063596);
// check x
VectorXd x_ref(n);
x_ref << 0.08241058, 1.133037, 2.343695;
VERIFY_IS_APPROX(x, x_ref);
}
struct lmder_functor {
static int f(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjac, int ldfjac, int iflag)
{
/* subroutine fcn for lmder example. */
int i;
double tmp1, tmp2, tmp3, tmp4;
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, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
if (iflag == 0)
{
/* insert print statements here when nprint is positive. */
return 0;
}
if (iflag != 2)
{
for (i=1; i <= 15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3));
}
}
else
{
for (i=1; i<=15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4 = tmp4*tmp4;
fjac[i-1 + ldfjac*(1-1)] = -1.;
fjac[i-1 + ldfjac*(2-1)] = tmp1*tmp2/tmp4;
fjac[i-1 + ldfjac*(3-1)] = tmp1*tmp3/tmp4;
};
}
return 0;
}
};
void testLmder()
{
const int m=15, n=3;
int info, nfev, njev;
double fnorm, covfac, covar_ftol;
Eigen::VectorXd x(n), fvec(m), diag(n);
Eigen::MatrixXd fjac;
VectorXi ipvt;
/* the following starting values provide a rough fit. */
x.setConstant(n, 1.);
// do the computation
info = ei_lmder<lmder_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return values
VERIFY( 1 == info);
VERIFY(nfev==6);
VERIFY(njev==5);
// check norm
fnorm = fvec.norm();
VERIFY_IS_APPROX(fnorm, 0.09063596);
// check x
VectorXd x_ref(n);
x_ref << 0.08241058, 1.133037, 2.343695;
VERIFY_IS_APPROX(x, x_ref);
// check covariance
covar_ftol = dpmpar(1);
covfac = fnorm*fnorm/(m-n);
Eigen::VectorXd wa(n);
covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
Eigen::MatrixXd cov_ref(n,n);
cov_ref <<
0.0001531202, 0.002869941, -0.002656662,
0.002869941, 0.09480935, -0.09098995,
-0.002656662, -0.09098995, 0.08778727;
// std::cout << fjac*covfac << std::endl;
Eigen::MatrixXd cov;
cov = covfac*fjac.corner<n,n>(TopLeft);
VERIFY_IS_APPROX( cov, cov_ref);
// TODO: why isn't this allowed ? :
// VERIFY_IS_APPROX( covfac*fjac.corner<n,n>(TopLeft) , cov_ref);
}
struct hybrj1_functor {
static int f(void * /*p*/, int n, const double *x, double *fvec, double *fjac, int ldfjac,
int iflag)
{
/* subroutine fcn for hybrj1 example. */
int j, k;
double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4;
if (iflag != 2)
{
for (k = 1; k <= n; k++)
{
temp = (three - two*x[k-1])*x[k-1];
temp1 = zero;
if (k != 1) temp1 = x[k-1-1];
temp2 = zero;
if (k != n) temp2 = x[k+1-1];
fvec[k-1] = temp - temp1 - two*temp2 + one;
}
}
else
{
for (k = 1; k <= n; k++)
{
for (j = 1; j <= n; j++)
{
fjac[k-1 + ldfjac*(j-1)] = zero;
}
fjac[k-1 + ldfjac*(k-1)] = three - four*x[k-1];
if (k != 1) fjac[k-1 + ldfjac*(k-1-1)] = -one;
if (k != n) fjac[k-1 + ldfjac*(k+1-1)] = -two;
}
}
return 0;
}
};
void testHybrj1()
{
const int n=9;
int info;
Eigen::VectorXd x(n), fvec, diag(n);
Eigen::MatrixXd fjac;
/* the following starting values provide a rough fit. */
x.setConstant(n, -1.);
// do the computation
info = ei_hybrj1<hybrj1_functor, double>(x,fvec, fjac);
// check return value
VERIFY( 1 == info);
// check norm
VERIFY_IS_APPROX(fvec.norm(), 1.192636e-08);
// check x
VectorXd x_ref(n);
x_ref <<
-0.5706545, -0.6816283, -0.7017325,
-0.7042129, -0.701369, -0.6918656,
-0.665792, -0.5960342, -0.4164121;
VERIFY_IS_APPROX(x, x_ref);
}
struct hybrj_functor {
static int f(void * /*p*/, int n, const double *x, double *fvec, double *fjac, int ldfjac,
int iflag)
{
/* subroutine fcn for hybrj example. */
int j, k;
double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4;
if (iflag == 0)
{
/* insert print statements here when nprint is positive. */
return 0;
}
if (iflag != 2)
{
for (k=1; k <= n; k++)
{
temp = (three - two*x[k-1])*x[k-1];
temp1 = zero;
if (k != 1) temp1 = x[k-1-1];
temp2 = zero;
if (k != n) temp2 = x[k+1-1];
fvec[k-1] = temp - temp1 - two*temp2 + one;
}
}
else
{
for (k = 1; k <= n; k++)
{
for (j=1; j <= n; j++)
{
fjac[k-1 + ldfjac*(j-1)] = zero;
}
fjac[k-1 + ldfjac*(k-1)] = three - four*x[k-1];
if (k != 1) fjac[k-1 + ldfjac*(k-1-1)] = -one;
if (k != n) fjac[k-1 + ldfjac*(k+1-1)] = -two;
}
}
return 0;
}
};
void testHybrj()
{
const int n=9;
int info, nfev, njev, mode;
Eigen::VectorXd x(n), fvec, diag(n), R, qtf;
Eigen::MatrixXd fjac;
/* the following starting values provide a rough fit. */
x.setConstant(n, -1.);
mode = 2;
diag.setConstant(n, 1.);
// do the computation
info = ei_hybrj<hybrj_functor, double>(x,fvec, nfev, njev, fjac, R, qtf, diag, mode);
// check return value
VERIFY( 1 == info);
VERIFY(nfev==11);
VERIFY(njev==1);
// check norm
VERIFY_IS_APPROX(fvec.norm(), 1.192636e-08);
// check x
VectorXd x_ref(n);
x_ref <<
-0.5706545, -0.6816283, -0.7017325,
-0.7042129, -0.701369, -0.6918656,
-0.665792, -0.5960342, -0.4164121;
VERIFY_IS_APPROX(x, x_ref);
}
struct hybrd1_functor {
static int f(void * /*p*/, int n, const double *x, double *fvec, int /*iflag*/)
{
/* subroutine fcn for hybrd1 example. */
int k;
double one=1, temp, temp1, temp2, three=3, two=2, zero=0;
for (k=1; k <= n; k++)
{
temp = (three - two*x[k-1])*x[k-1];
temp1 = zero;
if (k != 1) temp1 = x[k-1-1];
temp2 = zero;
if (k != n) temp2 = x[k+1-1];
fvec[k-1] = temp - temp1 - two*temp2 + one;
}
return 0;
}
};
void testHybrd1()
{
int n=9, info;
Eigen::VectorXd x(n), fvec(n);
/* the following starting values provide a rough solution. */
x.setConstant(n, -1.);
// do the computation
info = ei_hybrd1<hybrd1_functor,double>(x, fvec);
// check return value
VERIFY( 1 == info);
// check norm
VERIFY_IS_APPROX(fvec.norm(), 1.192636e-08);
// check x
VectorXd x_ref(n);
x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121;
VERIFY_IS_APPROX(x, x_ref);
}
struct hybrd_functor {
static int f(void * /*p*/, int n, const double *x, double *fvec, int iflag)
{
/* subroutine fcn for hybrd example. */
int k;
double one=1, temp, temp1, temp2, three=3, two=2, zero=0;
if (iflag == 0)
{
/* insert print statements here when nprint is positive. */
return 0;
}
for (k=1; k<=n; k++)
{
temp = (three - two*x[k-1])*x[k-1];
temp1 = zero;
if (k != 1) temp1 = x[k-1-1];
temp2 = zero;
if (k != n) temp2 = x[k+1-1];
fvec[k-1] = temp - temp1 - two*temp2 + one;
}
return 0;
}
};
void testHybrd()
{
const int n=9;
int info, nfev, ml, mu, mode;
Eigen::VectorXd x(n), fvec, diag(n), R, qtf;
Eigen::MatrixXd fjac;
/* the following starting values provide a rough fit. */
x.setConstant(n, -1.);
ml = 1;
mu = 1;
mode = 2;
diag.setConstant(n, 1.);
// do the computation
info = ei_hybrd<hybrd_functor, double>(x,fvec, nfev, fjac, R, qtf, diag, mode, ml, mu);
// check return value
VERIFY( 1 == info);
VERIFY(nfev==14);
// check norm
VERIFY_IS_APPROX(fvec.norm(), 1.192636e-08);
// check x
VectorXd x_ref(n);
x_ref <<
-0.5706545, -0.6816283, -0.7017325,
-0.7042129, -0.701369, -0.6918656,
-0.665792, -0.5960342, -0.4164121;
VERIFY_IS_APPROX(x, x_ref);
}
struct lmstr1_functor {
static int f(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjrow, int iflag)
{
/* subroutine fcn for lmstr1 example. */
int i;
double tmp1, tmp2, tmp3, tmp4;
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, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
if (iflag < 2)
{
for (i=1; i<=15; i++)
{
tmp1=i;
tmp2 = 16-i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3));
}
}
else
{
i = iflag - 1;
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4=tmp4*tmp4;
fjrow[1-1] = -1;
fjrow[2-1] = tmp1*tmp2/tmp4;
fjrow[3-1] = tmp1*tmp3/tmp4;
}
return 0;
}
};
void testLmstr1()
{
int m=15, n=3, info;
Eigen::VectorXd x(n), fvec(m);
VectorXi ipvt;
/* the following starting values provide a rough fit. */
x.setConstant(n, 1.);
// do the computation
info = ei_lmstr1<lmstr1_functor,double>(x, fvec, ipvt);
// check return value
VERIFY( 1 == info);
// check norm
VERIFY_IS_APPROX(fvec.norm(), 0.09063596);
// check x
VectorXd x_ref(n);
x_ref << 0.08241058, 1.133037, 2.343695 ;
VERIFY_IS_APPROX(x, x_ref);
}
struct lmstr_functor {
static int f(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, double *fjrow, int iflag)
{
/* subroutine fcn for lmstr example. */
int i;
double tmp1, tmp2, tmp3, tmp4;
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, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
if (iflag == 0)
{
/* insert print statements here when nprint is positive. */
return 0;
}
if (iflag < 2)
{
for (i = 1; i <= 15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3));
}
}
else
{
i = iflag - 1;
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
tmp4 = (x[2-1]*tmp2 + x[3-1]*tmp3); tmp4 = tmp4*tmp4;
fjrow[1-1] = -1.;
fjrow[2-1] = tmp1*tmp2/tmp4;
fjrow[3-1] = tmp1*tmp3/tmp4;
}
return 0;
}
};
void testLmstr()
{
const int m=15, n=3;
int info, nfev, njev;
double fnorm;
Eigen::VectorXd x(n), fvec(m), diag(n);
Eigen::MatrixXd fjac;
VectorXi ipvt;
/* the following starting values provide a rough fit. */
x.setConstant(n, 1.);
// do the computation
info = ei_lmstr<lmstr_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
Eigen::VectorXd wa(n);
// check return values
VERIFY( 1 == info);
VERIFY(nfev==6);
VERIFY(njev==5);
// check norm
fnorm = fvec.norm();
VERIFY_IS_APPROX(fnorm, 0.09063596);
// check x
VectorXd x_ref(n);
x_ref << 0.08241058, 1.133037, 2.343695;
VERIFY_IS_APPROX(x, x_ref);
}
struct lmdif1_functor {
static int f(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, int /*iflag*/)
{
/* function fcn for lmdif1 example */
int i;
double tmp1,tmp2,tmp3;
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,
3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0};
for (i=0; i<15; i++)
{
tmp1 = i+1;
tmp2 = 15 - i;
tmp3 = tmp1;
if (i >= 8) tmp3 = tmp2;
fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
}
return 0;
}
};
void testLmdif1()
{
int m=15, n=3, info;
Eigen::VectorXd x(n), fvec(m);
VectorXi iwa;
/* the following starting values provide a rough fit. */
x.setConstant(n, 1.);
// do the computation
info = ei_lmdif1<lmdif1_functor,double>(x, fvec, iwa);
// check return value
VERIFY( 1 == info);
// check norm
VERIFY_IS_APPROX(fvec.norm(), 0.09063596);
// check x
VectorXd x_ref(n);
x_ref << 0.0824106, 1.1330366, 2.3436947;
VERIFY_IS_APPROX(x, x_ref);
}
struct lmdif_functor {
static int f(void * /*p*/, int /*m*/, int /*n*/, const double *x, double *fvec, int iflag)
{
/* subroutine fcn for lmdif example. */
int i;
double tmp1, tmp2, tmp3;
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, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
if (iflag == 0)
{
/* insert print statements here when nprint is positive. */
return 0;
}
for (i = 1; i <= 15; i++)
{
tmp1 = i;
tmp2 = 16 - i;
tmp3 = tmp1;
if (i > 8) tmp3 = tmp2;
fvec[i-1] = y[i-1] - (x[1-1] + tmp1/(x[2-1]*tmp2 + x[3-1]*tmp3));
}
return 0;
}
};
void testLmdif()
{
const int m=15, n=3;
int info, nfev;
double fnorm, covfac, covar_ftol;
Eigen::VectorXd x(n), fvec(m), diag(n);
Eigen::MatrixXd fjac;
VectorXi ipvt;
/* the following starting values provide a rough fit. */
x.setConstant(n, 1.);
// do the computation
info = ei_lmdif<lmdif_functor, double>(x, fvec, nfev, fjac, ipvt, diag);
// check return values
VERIFY( 1 == info);
VERIFY(nfev==21);
// check norm
fnorm = fvec.norm();
VERIFY_IS_APPROX(fnorm, 0.09063596);
// check x
VectorXd x_ref(n);
x_ref << 0.08241058, 1.133037, 2.343695;
VERIFY_IS_APPROX(x, x_ref);
// check covariance
covar_ftol = dpmpar(1);
covfac = fnorm*fnorm/(m-n);
Eigen::VectorXd wa(n);
covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
Eigen::MatrixXd cov_ref(n,n);
cov_ref <<
0.0001531202, 0.002869942, -0.002656662,
0.002869942, 0.09480937, -0.09098997,
-0.002656662, -0.09098997, 0.08778729;
// std::cout << fjac*covfac << std::endl;
Eigen::MatrixXd cov;
cov = covfac*fjac.corner<n,n>(TopLeft);
VERIFY_IS_APPROX( cov, cov_ref);
// TODO: why isn't this allowed ? :
// VERIFY_IS_APPROX( covfac*fjac.corner<n,n>(TopLeft) , cov_ref);
}
struct chwirut2_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double _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};
static const double 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 };
int i;
assert(m==54);
assert(n==3);
assert(ldfjac==54);
if (iflag != 2) {// compute fvec at b
for(i=0; i<54; i++) {
double x = _x[i];
fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<54; i++) {
double x = _x[i];
double factor = 1./(b[1]+b[2]*x);
double e = exp(-b[0]*x);
fjac[i+ldfjac*0] = -x*e*factor;
fjac[i+ldfjac*1] = -e*factor*factor;
fjac[i+ldfjac*2] = -x*e*factor*factor;
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
void testNistChwirut2(void)
{
const int m=54, n=3;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 0.1, 0.01, 0.02;
// do the computation
info = ei_lmder<chwirut2_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 10 == nfev);
VERIFY( 8 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.1304802941E+02);
// check x
VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
/*
* Second try
*/
x<< 0.15, 0.008, 0.010;
// do the computation
info = ei_lmder<chwirut2_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag,
1, 100., 400, 1.E6*Eigen::machine_epsilon<double>(), 1.E6*Eigen::machine_epsilon<double>());
// check return value
VERIFY( 1 == info);
VERIFY( 7 == nfev);
VERIFY( 6 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.1304802941E+02);
// check x
VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
}
struct misra1a_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double 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};
static const double 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};
int i;
assert(m==14);
assert(n==2);
assert(ldfjac==14);
if (iflag != 2) {// compute fvec at b
for(i=0; i<14; i++) {
fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i] ;
}
}
else { // compute fjac at b
for(i=0; i<14; i++) {
fjac[i+ldfjac*0] = (1.-exp(-b[1]*x[i]));
fjac[i+ldfjac*1] = (b[0]*x[i]*exp(-b[1]*x[i]));
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
void testNistMisra1a(void)
{
const int m=14, n=2;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 500., 0.0001;
// do the computation
info = ei_lmder<misra1a_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 19 == nfev);
VERIFY( 15 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.2455138894E-01);
// check x
VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
/*
* Second try
*/
x<< 250., 0.0005;
// do the computation
info = ei_lmder<misra1a_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 5 == nfev);
VERIFY( 4 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.2455138894E-01);
// check x
VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
}
struct hahn1_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double _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};
static const double _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 };
int i;
// static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
assert(m==236);
assert(n==7);
assert(ldfjac==236);
if (iflag != 2) {// compute fvec at x
for(i=0; i<236; i++) {
double x=_x[i], xx=x*x, xxx=xx*x;
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];
}
}
else { // compute fjac at x
for(i=0; i<236; i++) {
double x=_x[i], xx=x*x, xxx=xx*x;
double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
fjac[i+ldfjac*0] = 1.*fact;
fjac[i+ldfjac*1] = x*fact;
fjac[i+ldfjac*2] = xx*fact;
fjac[i+ldfjac*3] = xxx*fact;
fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
fjac[i+ldfjac*4] = x*fact;
fjac[i+ldfjac*5] = xx*fact;
fjac[i+ldfjac*6] = xxx*fact;
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml
void testNistHahn1(void)
{
const int m=236, n=7;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
// do the computation
info = ei_lmder<hahn1_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 11== nfev);
VERIFY( 10== njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.5324382854E+00);
// check x
VERIFY_IS_APPROX(x[0], 1.0776351733E+00 );
VERIFY_IS_APPROX(x[1],-1.2269296921E-01 );
VERIFY_IS_APPROX(x[2], 4.0863750610E-03 );
VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06
VERIFY_IS_APPROX(x[4],-5.7609940901E-03 );
VERIFY_IS_APPROX(x[5], 2.4053735503E-04 );
VERIFY_IS_APPROX(x[6],-1.2314450199E-07 );
/*
* Second try
*/
x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001;
// do the computation
info = ei_lmder<hahn1_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 11 == nfev);
VERIFY( 10 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.5324382854E+00);
// check x
VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00
VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01
VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03
VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06
VERIFY_IS_APPROX(x[4],-5.7609940901E-03 );
VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04
VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07
}
struct misra1d_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double 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};
static const double 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};
int i;
assert(m==14);
assert(n==2);
assert(ldfjac==14);
if (iflag != 2) {// compute fvec at b
for(i=0; i<14; i++) {
fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<14; i++) {
double den = 1.+b[1]*x[i];
fjac[i+ldfjac*0] = b[1]*x[i] / den;
fjac[i+ldfjac*1] = b[0]*x[i]*(den-b[1]*x[i])/den/den;
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml
void testNistMisra1d(void)
{
const int m=14, n=2;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 500., 0.0001;
// do the computation
info = ei_lmder<misra1d_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 3 == info);
VERIFY( 9 == nfev);
VERIFY( 7 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.6419295283E-02);
// check x
VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
/*
* Second try
*/
x<< 450., 0.0003;
// do the computation
info = ei_lmder<misra1d_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 4 == nfev);
VERIFY( 3 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.6419295283E-02);
// check x
VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
}
struct lanczos1_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double 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 };
static const double 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 };
int i;
assert(m==24);
assert(n==6);
assert(ldfjac==24);
if (iflag != 2) {// compute fvec at b
for(i=0; i<24; i++) {
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];
}
}
else { // compute fjac at b
for(i=0; i<24; i++) {
fjac[i+ldfjac*0] = exp(-b[1]*x[i]);
fjac[i+ldfjac*1] = -b[0]*x[i]*exp(-b[1]*x[i]);
fjac[i+ldfjac*2] = exp(-b[3]*x[i]);
fjac[i+ldfjac*3] = -b[2]*x[i]*exp(-b[3]*x[i]);
fjac[i+ldfjac*4] = exp(-b[5]*x[i]);
fjac[i+ldfjac*5] = -b[4]*x[i]*exp(-b[5]*x[i]);
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml
void testNistLanczos1(void)
{
const int m=24, n=6;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6;
// do the computation
info = ei_lmder<lanczos1_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 2 == info);
VERIFY( 79 == nfev);
VERIFY( 72 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.4292388868910E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
// check x
VERIFY_IS_APPROX(x[0], 9.5100000027E-02 );
VERIFY_IS_APPROX(x[1], 1.0000000001E+00 );
VERIFY_IS_APPROX(x[2], 8.6070000013E-01 );
VERIFY_IS_APPROX(x[3], 3.0000000002E+00 );
VERIFY_IS_APPROX(x[4], 1.5575999998E+00 );
VERIFY_IS_APPROX(x[5], 5.0000000001E+00 );
/*
* Second try
*/
x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3;
// do the computation
info = ei_lmder<lanczos1_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 2 == info);
VERIFY( 9 == nfev);
VERIFY( 8 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.43049947737308E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
// check x
VERIFY_IS_APPROX(x[0], 9.5100000027E-02 );
VERIFY_IS_APPROX(x[1], 1.0000000001E+00 );
VERIFY_IS_APPROX(x[2], 8.6070000013E-01 );
VERIFY_IS_APPROX(x[3], 3.0000000002E+00 );
VERIFY_IS_APPROX(x[4], 1.5575999998E+00 );
VERIFY_IS_APPROX(x[5], 5.0000000001E+00 );
}
struct rat42_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 };
static const double y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 };
int i;
assert(m==9);
assert(n==3);
assert(ldfjac==9);
if (iflag != 2) {// compute fvec at b
for(i=0; i<9; i++) {
fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<9; i++) {
double e = exp(b[1]-b[2]*x[i]);
fjac[i+ldfjac*0] = 1./(1.+e);
fjac[i+ldfjac*1] = -b[0]*e/(1.+e)/(1.+e);
fjac[i+ldfjac*2] = +b[0]*e*x[i]/(1.+e)/(1.+e);
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml
void testNistRat42(void)
{
const int m=9, n=3;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 100., 1., 0.1;
// do the computation
info = ei_lmder<rat42_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 10 == nfev);
VERIFY( 8 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 8.0565229338E+00);
// check x
VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
/*
* Second try
*/
x<< 75., 2.5, 0.07;
// do the computation
info = ei_lmder<rat42_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 6 == nfev);
VERIFY( 5 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 8.0565229338E+00);
// check x
VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
}
struct MGH10_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double 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 };
static const double 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 };
int i;
assert(m==16);
assert(n==3);
assert(ldfjac==16);
if (iflag != 2) {// compute fvec at b
for(i=0; i<16; i++) {
fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<16; i++) {
double factor = 1./(x[i]+b[2]);
double e = exp(b[1]*factor);
fjac[i+ldfjac*0] = e;
fjac[i+ldfjac*1] = b[0]*factor*e;
fjac[i+ldfjac*2] = -b[1]*b[0]*factor*factor*e;
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
void testNistMGH10(void)
{
const int m=16, n=3;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 2., 400000., 25000.;
// do the computation
info = ei_lmder<MGH10_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 3 == info);
VERIFY( 285 == nfev);
VERIFY( 250 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 8.7945855171E+01);
// check x
VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
/*
* Second try
*/
x<< 0.02, 4000., 250.;
// do the computation
info = ei_lmder<MGH10_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 3 == info);
VERIFY( 126 == nfev);
VERIFY( 116 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 8.7945855171E+01);
// check x
VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
}
struct BoxBOD_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double x[6] = { 1., 2., 3., 5., 7., 10. };
static const double y[6] = { 109., 149., 149., 191., 213., 224. };
int i;
assert(m==6);
assert(n==2);
assert(ldfjac==6);
if (iflag != 2) {// compute fvec at b
for(i=0; i<6; i++) {
fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<6; i++) {
double e = exp(-b[1]*x[i]);
fjac[i+ldfjac*0] = 1.-e;
fjac[i+ldfjac*1] = b[0]*x[i]*e;
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml
void testNistBoxBOD(void)
{
const int m=6, n=2;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
#if 0
/*
* First try
*/
x<< 1., 1.;
// do the computation
info = ei_lmder<BoxBOD_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
1, 100., 14000, Eigen::machine_epsilon<double>(), Eigen::machine_epsilon<double>());
// check return value
printf("info=%d, f,j: %d, %d\n", info, nfev, njev);
printf("norm2 = %.50g\n", fvec.squaredNorm());
std::cout << x << std::endl;
VERIFY( 1 == info);
VERIFY( 10 == nfev);
VERIFY( 6 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.1680088766E+03);
// check x
VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
#endif
/*
* Second try
*/
x<< 100., 0.75;
// do the computation
info = ei_lmder<BoxBOD_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag,
1, 100., 14000, Eigen::machine_epsilon<double>(), Eigen::machine_epsilon<double>());
// check return value
VERIFY( 1 == info);
VERIFY( 16 == nfev);
VERIFY( 15 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.1680088766E+03);
// check x
VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
}
struct MGH17_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double 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 };
static const double 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 };
int i;
assert(m==33);
assert(n==5);
assert(ldfjac==33);
if (iflag != 2) {// compute fvec at b
for(i=0; i<33; i++) {
fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<33; i++) {
fjac[i+ldfjac*0] = 1.;
fjac[i+ldfjac*1] = exp(-b[3]*x[i]);
fjac[i+ldfjac*2] = exp(-b[4]*x[i]);
fjac[i+ldfjac*3] = -x[i]*b[1]*exp(-b[3]*x[i]);
fjac[i+ldfjac*4] = -x[i]*b[2]*exp(-b[4]*x[i]);
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml
void testNistMGH17(void)
{
const int m=33, n=5;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
#if 1
/*
* First try
*/
x<< 50., 150., -100., 1., 2.;
// do the computation
info = ei_lmder<MGH17_functor, double>(
x, fvec, nfev, njev, fjac, ipvt, diag,
1, 100., 5000, Eigen::machine_epsilon<double>(), Eigen::machine_epsilon<double>());
// check return value
VERIFY( 2 == info);
VERIFY( 604 == nfev);
VERIFY( 545 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.4648946975E-05);
// check x
VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
#endif
/*
* Second try
*/
x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02;
// do the computation
info = ei_lmder<MGH17_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 18 == nfev);
VERIFY( 15 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.4648946975E-05);
// check x
VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
}
struct MGH09_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double _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 };
static const double 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 };
int i;
assert(m==11);
assert(n==4);
assert(ldfjac==11);
if (iflag != 2) {// compute fvec at b
for(i=0; i<11; i++) {
double x = _x[i], xx=x*x;
fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<11; i++) {
double x = _x[i], xx=x*x;
double factor = 1./(xx+x*b[2]+b[3]);
fjac[i+ldfjac*0] = (xx+x*b[1]) * factor;
fjac[i+ldfjac*1] = b[0]*x* factor;
fjac[i+ldfjac*2] = - b[0]*(xx+x*b[1]) * x * factor * factor;
fjac[i+ldfjac*3] = - b[0]*(xx+x*b[1]) * factor * factor;
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
void testNistMGH09(void)
{
const int m=11, n=4;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 25., 39, 41.5, 39.;
// do the computation
info = ei_lmder<MGH09_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag,
1, 100., 5000);
// 1, 100., 5000, Eigen::machine_epsilon<double>(), Eigen::machine_epsilon<double>());
// check return value
VERIFY( 1 == info);
VERIFY( 487== nfev);
VERIFY( 378 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 3.0750560385E-04);
// check x
VERIFY_IS_APPROX(x[0], 0.19280590); // should be 1.9280693458E-01
VERIFY_IS_APPROX(x[1], 0.19130543); // should be 1.9128232873E-01
VERIFY_IS_APPROX(x[2], 0.12306085); // should be 1.2305650693E-01
VERIFY_IS_APPROX(x[3], 0.13607303); // should be 1.3606233068E-01
/*
* Second try
*/
x<< 0.25, 0.39, 0.415, 0.39;
// do the computation
info = ei_lmder<MGH09_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 18 == nfev);
VERIFY( 16 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 3.0750560385E-04);
// check x
VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01
VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01
VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01
VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01
}
struct Bennett5_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double 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 };
static const double 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 };
int i;
assert(m==154);
assert(n==3);
assert(ldfjac==154);
if (iflag != 2) {// compute fvec at b
for(i=0; i<154; i++) {
fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i];
}
}
else { // compute fjac at b
for(i=0; i<154; i++) {
double e = pow(b[1]+x[i],-1./b[2]);
fjac[i+ldfjac*0] = e;
fjac[i+ldfjac*1] = - b[0]*e/b[2]/(b[1]+x[i]);
fjac[i+ldfjac*2] = b[0]*e*log(b[1]+x[i])/b[2]/b[2];
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
void testNistBennett5(void)
{
const int m=154, n=3;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< -2000., 50., 0.8;
// do the computation
info = ei_lmder<Bennett5_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag,
1, 100., 5000);
// check return value
VERIFY( 1 == info);
VERIFY( 758 == nfev);
VERIFY( 744 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.2404744073E-04);
// check x
VERIFY_IS_APPROX(x[0], -2.5235058043E+03);
VERIFY_IS_APPROX(x[1], 4.6736564644E+01);
VERIFY_IS_APPROX(x[2], 9.3218483193E-01);
/*
* Second try
*/
x<< -1500., 45., 0.85;
// do the computation
info = ei_lmder<Bennett5_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 1 == info);
VERIFY( 203 == nfev);
VERIFY( 192 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.2404744073E-04);
// check x
VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03
VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01);
VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01);
}
struct thurber_functor {
static int f(void * /*p*/, int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
{
static const double _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 };
static const double _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};
int i;
// static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
assert(m==37);
assert(n==7);
assert(ldfjac==37);
if (iflag != 2) {// compute fvec at x
for(i=0; i<37; i++) {
double x=_x[i], xx=x*x, xxx=xx*x;
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];
}
}
else { // compute fjac at x
for(i=0; i<37; i++) {
double x=_x[i], xx=x*x, xxx=xx*x;
double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
fjac[i+ldfjac*0] = 1.*fact;
fjac[i+ldfjac*1] = x*fact;
fjac[i+ldfjac*2] = xx*fact;
fjac[i+ldfjac*3] = xxx*fact;
fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
fjac[i+ldfjac*4] = x*fact;
fjac[i+ldfjac*5] = xx*fact;
fjac[i+ldfjac*6] = xxx*fact;
}
}
return 0;
}
};
// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml
void testNistThurber(void)
{
const int m=37, n=7;
int info, nfev, njev;
Eigen::VectorXd x(n), fvec(m), diag;
Eigen::MatrixXd fjac;
VectorXi ipvt;
/*
* First try
*/
x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ;
// do the computation
info = ei_lmder<thurber_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag,
1, 100., 400, 1.E4*Eigen::machine_epsilon<double>(), 1.E4*Eigen::machine_epsilon<double>());
// check return value
VERIFY( 1 == info);
VERIFY( 39 == nfev);
VERIFY( 36== njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.6427082397E+03);
// check x
VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
/*
* Second try
*/
x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ;
// do the computation
info = ei_lmder<thurber_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag,
1, 100., 400, 1.E4*Eigen::machine_epsilon<double>(), 1.E4*Eigen::machine_epsilon<double>());
// check return value
VERIFY( 1 == info);
VERIFY( 29 == nfev);
VERIFY( 28 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 5.6427082397E+03);
// check x
VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
}
void test_NonLinear()
{
// NIST tests, level of difficulty = "Lower"
CALL_SUBTEST(testNistMisra1a());
CALL_SUBTEST(testNistChwirut2());
// NIST tests, level of difficulty = "Average"
CALL_SUBTEST(testNistHahn1());
CALL_SUBTEST(testNistMisra1d());
CALL_SUBTEST(testNistMGH17());
CALL_SUBTEST(testNistLanczos1());
// NIST tests, level of difficulty = "Higher"
CALL_SUBTEST(testNistRat42());
CALL_SUBTEST(testNistMGH10());
CALL_SUBTEST(testNistBoxBOD());
CALL_SUBTEST(testNistMGH09());
CALL_SUBTEST(testNistBennett5());
CALL_SUBTEST(testNistThurber());
CALL_SUBTEST(testChkder());
CALL_SUBTEST(testLmder1());
CALL_SUBTEST(testLmder());
CALL_SUBTEST(testHybrj1());
CALL_SUBTEST(testHybrj());
CALL_SUBTEST(testHybrd1());
CALL_SUBTEST(testHybrd());
CALL_SUBTEST(testLmstr1());
CALL_SUBTEST(testLmstr());
CALL_SUBTEST(testLmdif1());
CALL_SUBTEST(testLmdif());
}
Reference in New Issue View Git Blame Copy Permalink