mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-03-13 18:37:27 +08:00
beautify Functor for lmder : we now have f,df,debug methods
This commit is contained in:
Thomas Capricelli
parent
878f15b8a5
commit
acd757737a
@ -64,7 +64,7 @@ L20:
|
||||
/* evaluate the function at the starting point */
|
||||
/* and calculate its norm. */
|
||||
|
||||
iflag = Functor::f(x, fvec, fjac, 1);
|
||||
iflag = Functor::f(x, fvec);
|
||||
nfev = 1;
|
||||
if (iflag < 0) {
|
||||
goto L300;
|
||||
@ -82,7 +82,7 @@ L30:
|
||||
|
||||
/* calculate the jacobian matrix. */
|
||||
|
||||
iflag = Functor::f(x, fvec, fjac, 2);
|
||||
iflag = Functor::df(x, fjac);
|
||||
++njev;
|
||||
if (iflag < 0) {
|
||||
goto L300;
|
||||
@ -95,7 +95,7 @@ L30:
|
||||
}
|
||||
iflag = 0;
|
||||
if ((iter - 1) % nprint == 0) {
|
||||
iflag = Functor::f(x, fvec, fjac, 0);
|
||||
iflag = Functor::debug(x, fvec, fjac);
|
||||
}
|
||||
if (iflag < 0) {
|
||||
goto L300;
|
||||
@ -237,7 +237,7 @@ L200:
|
||||
|
||||
/* evaluate the function at x + p and calculate its norm. */
|
||||
|
||||
iflag = Functor::f(wa2, wa4, fjac, 1);
|
||||
iflag = Functor::f(wa2, wa4);
|
||||
++nfev;
|
||||
if (iflag < 0) {
|
||||
goto L300;
|
||||
@ -378,7 +378,7 @@ L300:
|
||||
}
|
||||
iflag = 0;
|
||||
if (nprint > 0) {
|
||||
iflag = Functor::f(x, fvec, fjac, 0);
|
||||
iflag = Functor::debug(x, fvec, fjac);
|
||||
}
|
||||
return info;
|
||||
|
||||
|
||||
@ -102,38 +102,35 @@ void testChkder()
|
||||
|
||||
|
||||
struct lmder_functor {
|
||||
static int f(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &x, VectorXd &fvec)
|
||||
{
|
||||
|
||||
/* subroutine fcn for lmder1 example. */
|
||||
|
||||
int i;
|
||||
double tmp1, tmp2, tmp3, tmp4;
|
||||
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 != 2)
|
||||
for (int i = 0; i < 15; i++)
|
||||
{
|
||||
for (i = 0; i < 15; i++)
|
||||
{
|
||||
tmp1 = i+1;
|
||||
tmp2 = 16 - i - 1;
|
||||
tmp3 = (i>=8)? tmp2 : tmp1;
|
||||
fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
|
||||
}
|
||||
tmp1 = i+1;
|
||||
tmp2 = 16 - i - 1;
|
||||
tmp3 = (i>=8)? tmp2 : tmp1;
|
||||
fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
|
||||
}
|
||||
else
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int df(const VectorXd &x, MatrixXd &fjac)
|
||||
{
|
||||
double tmp1, tmp2, tmp3, tmp4;
|
||||
for (int i = 0; i < 15; i++)
|
||||
{
|
||||
for ( i = 0; i < 15; i++)
|
||||
{
|
||||
tmp1 = i+1;
|
||||
tmp2 = 16 - i - 1;
|
||||
tmp3 = (i>=8)? tmp2 : tmp1;
|
||||
tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
|
||||
fjac(i,0) = -1;
|
||||
fjac(i,1) = tmp1*tmp2/tmp4;
|
||||
fjac(i,2) = tmp1*tmp3/tmp4;
|
||||
}
|
||||
tmp1 = i+1;
|
||||
tmp2 = 16 - i - 1;
|
||||
tmp3 = (i>=8)? tmp2 : tmp1;
|
||||
tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
|
||||
fjac(i,0) = -1;
|
||||
fjac(i,1) = tmp1*tmp2/tmp4;
|
||||
fjac(i,2) = tmp1*tmp3/tmp4;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
@ -598,35 +595,39 @@ void testLmdif()
|
||||
}
|
||||
|
||||
struct chwirut2_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double m_x[54];
|
||||
static const double m_y[54];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==3);
|
||||
assert(fvec.size()==54);
|
||||
for(i=0; i<54; i++) {
|
||||
double x = m_x[i];
|
||||
fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==3);
|
||||
assert(fjac.rows()==54);
|
||||
assert(fjac.cols()==3);
|
||||
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,0) = -x*e*factor;
|
||||
fjac(i,1) = -e*factor*factor;
|
||||
fjac(i,2) = -x*e*factor*factor;
|
||||
}
|
||||
for(int i=0; i<54; i++) {
|
||||
double x = m_x[i];
|
||||
double factor = 1./(b[1]+b[2]*x);
|
||||
double e = exp(-b[0]*x);
|
||||
fjac(i,0) = -x*e*factor;
|
||||
fjac(i,1) = -e*factor*factor;
|
||||
fjac(i,2) = -x*e*factor*factor;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double chwirut2_functor::m_x[54] = { 0.500E0, 1.000E0, 1.750E0, 3.750E0, 5.750E0, 0.875E0, 2.250E0, 3.250E0, 5.250E0, 0.750E0, 1.750E0, 2.750E0, 4.750E0, 0.625E0, 1.250E0, 2.250E0, 4.250E0, .500E0, 3.000E0, .750E0, 3.000E0, 1.500E0, 6.000E0, 3.000E0, 6.000E0, 1.500E0, 3.000E0, .500E0, 2.000E0, 4.000E0, .750E0, 2.000E0, 5.000E0, .750E0, 2.250E0, 3.750E0, 5.750E0, 3.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .500E0, 6.000E0, 3.000E0, .500E0, 2.750E0, .500E0, 1.750E0};
|
||||
const double chwirut2_functor::m_y[54] = { 92.9000E0 ,57.1000E0 ,31.0500E0 ,11.5875E0 ,8.0250E0 ,63.6000E0 ,21.4000E0 ,14.2500E0 ,8.4750E0 ,63.8000E0 ,26.8000E0 ,16.4625E0 ,7.1250E0 ,67.3000E0 ,41.0000E0 ,21.1500E0 ,8.1750E0 ,81.5000E0 ,13.1200E0 ,59.9000E0 ,14.6200E0 ,32.9000E0 ,5.4400E0 ,12.5600E0 ,5.4400E0 ,32.0000E0 ,13.9500E0 ,75.8000E0 ,20.0000E0 ,10.4200E0 ,59.5000E0 ,21.6700E0 ,8.5500E0 ,62.0000E0 ,20.2000E0 ,7.7600E0 ,3.7500E0 ,11.8100E0 ,54.7000E0 ,23.7000E0 ,11.5500E0 ,61.3000E0 ,17.7000E0 ,8.7400E0 ,59.2000E0 ,16.3000E0 ,8.6200E0 ,81.0000E0 ,4.8700E0 ,14.6200E0 ,81.7000E0 ,17.1700E0 ,81.3000E0 ,28.9000E0 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml
|
||||
void testNistChwirut2(void)
|
||||
@ -678,30 +679,32 @@ void testNistChwirut2(void)
|
||||
|
||||
|
||||
struct misra1a_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double m_x[14];
|
||||
static const double m_y[14];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==2);
|
||||
assert(fvec.size()==14);
|
||||
for(int i=0; i<14; i++) {
|
||||
fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==2);
|
||||
assert(fjac.rows()==14);
|
||||
assert(fjac.cols()==2);
|
||||
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,0) = (1.-exp(-b[1]*x[i]));
|
||||
fjac(i,1) = (b[0]*x[i]*exp(-b[1]*x[i]));
|
||||
}
|
||||
for(int i=0; i<14; i++) {
|
||||
fjac(i,0) = (1.-exp(-b[1]*m_x[i]));
|
||||
fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i]));
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double misra1a_functor::m_x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0};
|
||||
const double misra1a_functor::m_y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0};
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
|
||||
void testNistMisra1a(void)
|
||||
@ -749,41 +752,44 @@ void testNistMisra1a(void)
|
||||
}
|
||||
|
||||
struct hahn1_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double m_x[236];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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 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 };
|
||||
|
||||
// static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
|
||||
// static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
|
||||
|
||||
assert(b.size()==7);
|
||||
assert(fvec.size()==236);
|
||||
for(int i=0; i<236; i++) {
|
||||
double x=m_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) - m_y[i];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==7);
|
||||
assert(fjac.rows()==236);
|
||||
assert(fjac.cols()==7);
|
||||
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,0) = 1.*fact;
|
||||
fjac(i,1) = x*fact;
|
||||
fjac(i,2) = xx*fact;
|
||||
fjac(i,3) = xxx*fact;
|
||||
fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
|
||||
fjac(i,4) = x*fact;
|
||||
fjac(i,5) = xx*fact;
|
||||
fjac(i,6) = xxx*fact;
|
||||
}
|
||||
for(int i=0; i<236; i++) {
|
||||
double x=m_x[i], xx=x*x, xxx=xx*x;
|
||||
double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
|
||||
fjac(i,0) = 1.*fact;
|
||||
fjac(i,1) = x*fact;
|
||||
fjac(i,2) = xx*fact;
|
||||
fjac(i,3) = xxx*fact;
|
||||
fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
|
||||
fjac(i,4) = x*fact;
|
||||
fjac(i,5) = xx*fact;
|
||||
fjac(i,6) = xxx*fact;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , 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};
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml
|
||||
void testNistHahn1(void)
|
||||
@ -842,31 +848,33 @@ void testNistHahn1(void)
|
||||
}
|
||||
|
||||
struct misra1d_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[14];
|
||||
static const double y[14];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==2);
|
||||
assert(fvec.size()==14);
|
||||
for(int i=0; i<14; i++) {
|
||||
fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==2);
|
||||
assert(fjac.rows()==14);
|
||||
assert(fjac.cols()==2);
|
||||
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,0) = b[1]*x[i] / den;
|
||||
fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den;
|
||||
}
|
||||
for(int i=0; i<14; i++) {
|
||||
double den = 1.+b[1]*x[i];
|
||||
fjac(i,0) = b[1]*x[i] / den;
|
||||
fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double misra1d_functor::x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0};
|
||||
const double misra1d_functor::y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0};
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml
|
||||
void testNistMisra1d(void)
|
||||
@ -915,34 +923,35 @@ void testNistMisra1d(void)
|
||||
|
||||
|
||||
struct lanczos1_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[24];
|
||||
static const double y[24];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==6);
|
||||
assert(fvec.size()==24);
|
||||
for(int 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];
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==6);
|
||||
assert(fjac.rows()==24);
|
||||
assert(fjac.cols()==6);
|
||||
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,0) = exp(-b[1]*x[i]);
|
||||
fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]);
|
||||
fjac(i,2) = exp(-b[3]*x[i]);
|
||||
fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]);
|
||||
fjac(i,4) = exp(-b[5]*x[i]);
|
||||
fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]);
|
||||
}
|
||||
for(int i=0; i<24; i++) {
|
||||
fjac(i,0) = exp(-b[1]*x[i]);
|
||||
fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]);
|
||||
fjac(i,2) = exp(-b[3]*x[i]);
|
||||
fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]);
|
||||
fjac(i,4) = exp(-b[5]*x[i]);
|
||||
fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double lanczos1_functor::x[24] = { 0.000000000000E+00, 5.000000000000E-02, 1.000000000000E-01, 1.500000000000E-01, 2.000000000000E-01, 2.500000000000E-01, 3.000000000000E-01, 3.500000000000E-01, 4.000000000000E-01, 4.500000000000E-01, 5.000000000000E-01, 5.500000000000E-01, 6.000000000000E-01, 6.500000000000E-01, 7.000000000000E-01, 7.500000000000E-01, 8.000000000000E-01, 8.500000000000E-01, 9.000000000000E-01, 9.500000000000E-01, 1.000000000000E+00, 1.050000000000E+00, 1.100000000000E+00, 1.150000000000E+00 };
|
||||
const double lanczos1_functor::y[24] = { 2.513400000000E+00 ,2.044333373291E+00 ,1.668404436564E+00 ,1.366418021208E+00 ,1.123232487372E+00 ,9.268897180037E-01 ,7.679338563728E-01 ,6.388775523106E-01 ,5.337835317402E-01 ,4.479363617347E-01 ,3.775847884350E-01 ,3.197393199326E-01 ,2.720130773746E-01 ,2.324965529032E-01 ,1.996589546065E-01 ,1.722704126914E-01 ,1.493405660168E-01 ,1.300700206922E-01 ,1.138119324644E-01 ,1.000415587559E-01 ,8.833209084540E-02 ,7.833544019350E-02 ,6.976693743449E-02 ,6.239312536719E-02 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml
|
||||
void testNistLanczos1(void)
|
||||
@ -999,32 +1008,35 @@ void testNistLanczos1(void)
|
||||
}
|
||||
|
||||
struct rat42_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[9];
|
||||
static const double y[9];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==3);
|
||||
assert(fvec.size()==9);
|
||||
for(int i=0; i<9; i++) {
|
||||
fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==3);
|
||||
assert(fjac.rows()==9);
|
||||
assert(fjac.cols()==3);
|
||||
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,0) = 1./(1.+e);
|
||||
fjac(i,1) = -b[0]*e/(1.+e)/(1.+e);
|
||||
fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e);
|
||||
}
|
||||
for(int i=0; i<9; i++) {
|
||||
double e = exp(b[1]-b[2]*x[i]);
|
||||
fjac(i,0) = 1./(1.+e);
|
||||
fjac(i,1) = -b[0]*e/(1.+e)/(1.+e);
|
||||
fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 };
|
||||
const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml
|
||||
void testNistRat42(void)
|
||||
@ -1074,33 +1086,34 @@ void testNistRat42(void)
|
||||
}
|
||||
|
||||
struct MGH10_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[16];
|
||||
static const double y[16];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==3);
|
||||
assert(fvec.size()==16);
|
||||
for(int i=0; i<16; i++)
|
||||
fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i];
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==3);
|
||||
assert(fjac.rows()==16);
|
||||
assert(fjac.cols()==3);
|
||||
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,0) = e;
|
||||
fjac(i,1) = b[0]*factor*e;
|
||||
fjac(i,2) = -b[1]*b[0]*factor*factor*e;
|
||||
}
|
||||
for(int i=0; i<16; i++) {
|
||||
double factor = 1./(x[i]+b[2]);
|
||||
double e = exp(b[1]*factor);
|
||||
fjac(i,0) = e;
|
||||
fjac(i,1) = b[0]*factor*e;
|
||||
fjac(i,2) = -b[1]*b[0]*factor*factor*e;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double MGH10_functor::x[16] = { 5.000000E+01, 5.500000E+01, 6.000000E+01, 6.500000E+01, 7.000000E+01, 7.500000E+01, 8.000000E+01, 8.500000E+01, 9.000000E+01, 9.500000E+01, 1.000000E+02, 1.050000E+02, 1.100000E+02, 1.150000E+02, 1.200000E+02, 1.250000E+02 };
|
||||
const double MGH10_functor::y[16] = { 3.478000E+04, 2.861000E+04, 2.365000E+04, 1.963000E+04, 1.637000E+04, 1.372000E+04, 1.154000E+04, 9.744000E+03, 8.261000E+03, 7.030000E+03, 6.005000E+03, 5.147000E+03, 4.427000E+03, 3.820000E+03, 3.307000E+03, 2.872000E+03 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
|
||||
void testNistMGH10(void)
|
||||
@ -1151,31 +1164,31 @@ void testNistMGH10(void)
|
||||
|
||||
|
||||
struct BoxBOD_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[6];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==2);
|
||||
assert(fvec.size()==6);
|
||||
for(int i=0; i<6; i++)
|
||||
fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i];
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==2);
|
||||
assert(fjac.rows()==6);
|
||||
assert(fjac.cols()==2);
|
||||
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,0) = 1.-e;
|
||||
fjac(i,1) = b[0]*x[i]*e;
|
||||
}
|
||||
for(int i=0; i<6; i++) {
|
||||
double e = exp(-b[1]*x[i]);
|
||||
fjac(i,0) = 1.-e;
|
||||
fjac(i,1) = b[0]*x[i]*e;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml
|
||||
void testNistBoxBOD(void)
|
||||
@ -1225,33 +1238,34 @@ void testNistBoxBOD(void)
|
||||
}
|
||||
|
||||
struct MGH17_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[33];
|
||||
static const double y[33];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==5);
|
||||
assert(fvec.size()==33);
|
||||
for(int 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];
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==5);
|
||||
assert(fjac.rows()==33);
|
||||
assert(fjac.cols()==5);
|
||||
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,0) = 1.;
|
||||
fjac(i,1) = exp(-b[3]*x[i]);
|
||||
fjac(i,2) = exp(-b[4]*x[i]);
|
||||
fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]);
|
||||
fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]);
|
||||
}
|
||||
for(int i=0; i<33; i++) {
|
||||
fjac(i,0) = 1.;
|
||||
fjac(i,1) = exp(-b[3]*x[i]);
|
||||
fjac(i,2) = exp(-b[4]*x[i]);
|
||||
fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]);
|
||||
fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double MGH17_functor::x[33] = { 0.000000E+00, 1.000000E+01, 2.000000E+01, 3.000000E+01, 4.000000E+01, 5.000000E+01, 6.000000E+01, 7.000000E+01, 8.000000E+01, 9.000000E+01, 1.000000E+02, 1.100000E+02, 1.200000E+02, 1.300000E+02, 1.400000E+02, 1.500000E+02, 1.600000E+02, 1.700000E+02, 1.800000E+02, 1.900000E+02, 2.000000E+02, 2.100000E+02, 2.200000E+02, 2.300000E+02, 2.400000E+02, 2.500000E+02, 2.600000E+02, 2.700000E+02, 2.800000E+02, 2.900000E+02, 3.000000E+02, 3.100000E+02, 3.200000E+02 };
|
||||
const double MGH17_functor::y[33] = { 8.440000E-01, 9.080000E-01, 9.320000E-01, 9.360000E-01, 9.250000E-01, 9.080000E-01, 8.810000E-01, 8.500000E-01, 8.180000E-01, 7.840000E-01, 7.510000E-01, 7.180000E-01, 6.850000E-01, 6.580000E-01, 6.280000E-01, 6.030000E-01, 5.800000E-01, 5.580000E-01, 5.380000E-01, 5.220000E-01, 5.060000E-01, 4.900000E-01, 4.780000E-01, 4.670000E-01, 4.570000E-01, 4.480000E-01, 4.380000E-01, 4.310000E-01, 4.240000E-01, 4.200000E-01, 4.140000E-01, 4.110000E-01, 4.060000E-01 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml
|
||||
void testNistMGH17(void)
|
||||
@ -1309,35 +1323,37 @@ void testNistMGH17(void)
|
||||
}
|
||||
|
||||
struct MGH09_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double _x[11];
|
||||
static const double y[11];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==4);
|
||||
assert(fvec.size()==11);
|
||||
for(int 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];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==4);
|
||||
assert(fjac.rows()==11);
|
||||
assert(fjac.cols()==4);
|
||||
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,0) = (xx+x*b[1]) * factor;
|
||||
fjac(i,1) = b[0]*x* factor;
|
||||
fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor;
|
||||
fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor;
|
||||
}
|
||||
for(int i=0; i<11; i++) {
|
||||
double x = _x[i], xx=x*x;
|
||||
double factor = 1./(xx+x*b[2]+b[3]);
|
||||
fjac(i,0) = (xx+x*b[1]) * factor;
|
||||
fjac(i,1) = b[0]*x* factor;
|
||||
fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor;
|
||||
fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double MGH09_functor::_x[11] = { 4., 2., 1., 5.E-1 , 2.5E-01, 1.670000E-01, 1.250000E-01, 1.E-01, 8.330000E-02, 7.140000E-02, 6.250000E-02 };
|
||||
const double MGH09_functor::y[11] = { 1.957000E-01, 1.947000E-01, 1.735000E-01, 1.600000E-01, 8.440000E-02, 6.270000E-02, 4.560000E-02, 3.420000E-02, 3.230000E-02, 2.350000E-02, 2.460000E-02 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
|
||||
void testNistMGH09(void)
|
||||
@ -1393,33 +1409,33 @@ void testNistMGH09(void)
|
||||
|
||||
|
||||
struct Bennett5_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[154];
|
||||
static const double y[154];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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(b.size()==3);
|
||||
assert(fvec.size()==154);
|
||||
for(int i=0; i<154; i++)
|
||||
fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i];
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==3);
|
||||
assert(fjac.rows()==154);
|
||||
assert(fjac.cols()==3);
|
||||
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,0) = e;
|
||||
fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]);
|
||||
fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2];
|
||||
}
|
||||
for(int i=0; i<154; i++) {
|
||||
double e = pow(b[1]+x[i],-1./b[2]);
|
||||
fjac(i,0) = e;
|
||||
fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]);
|
||||
fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, 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 };
|
||||
const double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 ,-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 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml
|
||||
void testNistBennett5(void)
|
||||
@ -1469,40 +1485,42 @@ void testNistBennett5(void)
|
||||
}
|
||||
|
||||
struct thurber_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double _x[37];
|
||||
static const double _y[37];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
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++;
|
||||
// static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
|
||||
assert(b.size()==7);
|
||||
assert(fvec.size()==37);
|
||||
for(int 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];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==7);
|
||||
assert(fjac.rows()==37);
|
||||
assert(fjac.cols()==7);
|
||||
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,0) = 1.*fact;
|
||||
fjac(i,1) = x*fact;
|
||||
fjac(i,2) = xx*fact;
|
||||
fjac(i,3) = xxx*fact;
|
||||
fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
|
||||
fjac(i,4) = x*fact;
|
||||
fjac(i,5) = xx*fact;
|
||||
fjac(i,6) = xxx*fact;
|
||||
}
|
||||
for(int 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,0) = 1.*fact;
|
||||
fjac(i,1) = x*fact;
|
||||
fjac(i,2) = xx*fact;
|
||||
fjac(i,3) = xxx*fact;
|
||||
fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
|
||||
fjac(i,4) = x*fact;
|
||||
fjac(i,5) = xx*fact;
|
||||
fjac(i,6) = xxx*fact;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double thurber_functor::_x[37] = { -3.067E0, -2.981E0, -2.921E0, -2.912E0, -2.840E0, -2.797E0, -2.702E0, -2.699E0, -2.633E0, -2.481E0, -2.363E0, -2.322E0, -1.501E0, -1.460E0, -1.274E0, -1.212E0, -1.100E0, -1.046E0, -0.915E0, -0.714E0, -0.566E0, -0.545E0, -0.400E0, -0.309E0, -0.109E0, -0.103E0, 0.010E0, 0.119E0, 0.377E0, 0.790E0, 0.963E0, 1.006E0, 1.115E0, 1.572E0, 1.841E0, 2.047E0, 2.200E0 };
|
||||
const double thurber_functor::_y[37] = { 80.574E0, 84.248E0, 87.264E0, 87.195E0, 89.076E0, 89.608E0, 89.868E0, 90.101E0, 92.405E0, 95.854E0, 100.696E0, 101.060E0, 401.672E0, 390.724E0, 567.534E0, 635.316E0, 733.054E0, 759.087E0, 894.206E0, 990.785E0, 1090.109E0, 1080.914E0, 1122.643E0, 1178.351E0, 1260.531E0, 1273.514E0, 1288.339E0, 1327.543E0, 1353.863E0, 1414.509E0, 1425.208E0, 1421.384E0, 1442.962E0, 1464.350E0, 1468.705E0, 1447.894E0, 1457.628E0};
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml
|
||||
void testNistThurber(void)
|
||||
@ -1562,34 +1580,35 @@ void testNistThurber(void)
|
||||
}
|
||||
|
||||
struct rat43_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[15];
|
||||
static const double y[15];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
static const double x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. };
|
||||
static const double 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 };
|
||||
int i;
|
||||
|
||||
assert(b.size()==4);
|
||||
assert(fvec.size()==15);
|
||||
for(int i=0; i<15; i++)
|
||||
fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i];
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==4);
|
||||
assert(fjac.rows()==15);
|
||||
assert(fjac.cols()==4);
|
||||
if (iflag != 2) {// compute fvec at b
|
||||
for(i=0; i<15; i++) {
|
||||
fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i];
|
||||
}
|
||||
}
|
||||
else { // compute fjac at b
|
||||
for(i=0; i<15; i++) {
|
||||
double e = exp(b[1]-b[2]*x[i]);
|
||||
double power = -1./b[3];
|
||||
fjac(i,0) = pow(1.+e, power);
|
||||
fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.);
|
||||
fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.);
|
||||
fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power);
|
||||
}
|
||||
for(int i=0; i<15; i++) {
|
||||
double e = exp(b[1]-b[2]*x[i]);
|
||||
double power = -1./b[3];
|
||||
fjac(i,0) = pow(1.+e, power);
|
||||
fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.);
|
||||
fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.);
|
||||
fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. };
|
||||
const double rat43_functor::y[15] = { 16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml
|
||||
void testNistRat43(void)
|
||||
@ -1645,34 +1664,34 @@ void testNistRat43(void)
|
||||
|
||||
|
||||
struct eckerle4_functor {
|
||||
static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
|
||||
static const double x[35];
|
||||
static const double y[35];
|
||||
static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
|
||||
static int f(const VectorXd &b, VectorXd &fvec)
|
||||
{
|
||||
static const double 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};
|
||||
static const double 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 };
|
||||
|
||||
int i;
|
||||
|
||||
assert(b.size()==3);
|
||||
assert(fvec.size()==35);
|
||||
for(int i=0; i<35; i++)
|
||||
fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i];
|
||||
return 0;
|
||||
}
|
||||
static int df(const VectorXd &b, MatrixXd &fjac)
|
||||
{
|
||||
assert(b.size()==3);
|
||||
assert(fjac.rows()==35);
|
||||
assert(fjac.cols()==3);
|
||||
if (iflag != 2) {// compute fvec at b
|
||||
for(i=0; i<35; i++) {
|
||||
fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i];
|
||||
}
|
||||
}
|
||||
else { // compute fjac at b
|
||||
for(i=0; i<35; i++) {
|
||||
double b12 = b[1]*b[1];
|
||||
double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12);
|
||||
fjac(i,0) = e / b[1];
|
||||
fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12;
|
||||
fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12;
|
||||
}
|
||||
for(int i=0; i<35; i++) {
|
||||
double b12 = b[1]*b[1];
|
||||
double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12);
|
||||
fjac(i,0) = e / b[1];
|
||||
fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12;
|
||||
fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
const double eckerle4_functor::x[35] = { 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0};
|
||||
const double eckerle4_functor::y[35] = { 0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 };
|
||||
|
||||
// http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml
|
||||
void testNistEckerle4(void)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user