/* multifit_nlinear/test_fdf.c
*
* Copyright (C) 2007, 2013, 2014, 2015, 2016 Brian Gough, Patrick Alken
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
typedef struct
{
const char *name;
double *x0; /* initial parameters (size p) */
double *weights; /* data weights */
double *sigma;
double *epsrel; /* relative tolerance for solution checking */
void (*checksol) (const double x[], const double sumsq,
const double epsrel, const char *sname,
const char *pname);
gsl_multifit_nlinear_fdf *fdf;
} test_fdf_problem;
#include "test_bard.c"
#include "test_beale.c"
#include "test_biggs.c"
#include "test_box.c"
#include "test_boxbod.c"
#include "test_brown1.c"
#include "test_brown2.c"
#include "test_brown3.c"
#include "test_eckerle.c"
#include "test_enso.c"
#include "test_exp1.c"
#include "test_gaussian.c"
#include "test_hahn1.c"
#include "test_helical.c"
#include "test_jennrich.c"
#include "test_kirby2.c"
#include "test_kowalik.c"
#include "test_lin1.c"
#include "test_lin2.c"
#include "test_lin3.c"
#include "test_meyer.c"
#include "test_meyerscal.c"
#include "test_osborne.c"
#include "test_penalty1.c"
#include "test_penalty2.c"
#include "test_powell1.c"
#include "test_powell2.c"
#include "test_powell3.c"
#include "test_rat42.c"
#include "test_rat43.c"
#include "test_rosenbrock.c"
#include "test_rosenbrocke.c"
#include "test_roth.c"
#include "test_thurber.c"
#include "test_vardim.c"
#include "test_watson.c"
#include "test_wood.c"
#include "test_wnlin.c"
static void test_fdf(const gsl_multifit_nlinear_type * T,
const gsl_multifit_nlinear_parameters * params,
const double xtol, const double gtol,
const double ftol,
const double epsrel,
test_fdf_problem *problem);
static void test_fdf_checksol(const char *sname, const char *pname,
const double epsrel,
gsl_multifit_nlinear_workspace *s,
test_fdf_problem *problem);
/*
* FIXME: some test problems are disabled since they fail on certain
* solvers. Known failures are:
*
* Method test-problem
* ====== ============
* dogleg thurbera
* dogleg rat43a
* all boxboda
*/
static test_fdf_problem *test_problems[] = {
/*
* These test problems are taken from
*
* H. B. Nielsen, UCTP test problems for unconstrained optimization,
* IMM Department of Mathematical Modeling, Tech. Report
* IMM-REP-2000-17, 2000.
*/
&lin1_problem, /* 1 */
&lin2_problem, /* 2 */
&lin3_problem, /* 3 */
&rosenbrock_problem, /* 4 */
&helical_problem, /* 5 */
&powell1_problem, /* 6 */
&roth_problem, /* 7 */
&bard_problem, /* 8 */
&kowalik_problem, /* 9 */
&meyer_problem, /* 10 */
&watson_problem, /* 11 */
&box_problem, /* 12 */
&jennrich_problem, /* 13 */
&brown1_problem, /* 14 */
&brown2_problem, /* 16 */
&osborne_problem, /* 17 */
&exp1_problem, /* 18 */
&meyerscal_problem, /* 20 */
&powell2_problem,
/*
* These tests are from
*
* J. J. More, B. S. Garbow and K. E. Hillstrom, Testing
* Unconstrained Optimization Software, ACM Trans. Math. Soft.
* Vol 7, No 1, 1981.
*
* Many of these overlap with the Nielsen tests
*/
&rosenbrock_problem, /* 1 */
&roth_problem, /* 2 */
&powell3_problem, /* 3 */
&brown3_problem, /* 4 */
&beale_problem, /* 5 */
&jennrich_problem, /* 6 */
&helical_problem, /* 7 */
&bard_problem, /* 8 */
&gaussian_problem, /* 9 */
&meyer_problem, /* 10 */
&box_problem, /* 12 */
&powell1_problem, /* 13 */
&wood_problem, /* 14 */
&kowalik_problem, /* 15 */
&brown1_problem, /* 16 */
&osborne_problem, /* 17 */
&biggs_problem, /* 18 */
&watson_problem, /* 20 */
&rosenbrocke_problem, /* 21 */
&penalty1_problem, /* 23 */
&penalty2_problem, /* 24 */
&vardim_problem, /* 25 */
&brown2_problem, /* 27 */
&lin1_problem, /* 32 */
&lin2_problem, /* 33 */
&lin3_problem, /* 34 */
/* NIST test cases */
&kirby2a_problem,
&kirby2b_problem,
&hahn1a_problem,
&hahn1b_problem,
&ensoa_problem,
&ensob_problem,
/*&thurbera_problem,*/
&thurberb_problem,
/*&boxboda_problem,*/
&boxbodb_problem,
&rat42a_problem,
&rat42b_problem,
&eckerlea_problem,
&eckerleb_problem,
/*&rat43a_problem,*/
&rat43b_problem,
/* weighted test cases */
&wnlin_problem1,
&wnlin_problem2,
NULL
};
static void
test_fdf_main(const gsl_multifit_nlinear_parameters * params)
{
const double xtol = pow(GSL_DBL_EPSILON, 0.9);
const double gtol = pow(GSL_DBL_EPSILON, 0.9);
const double ftol = 0.0;
size_t i;
for (i = 0; test_problems[i] != NULL; ++i)
{
test_fdf_problem *problem = test_problems[i];
double epsrel = *(problem->epsrel);
gsl_multifit_nlinear_fdf fdf;
test_fdf(gsl_multifit_nlinear_trust, params, xtol, gtol, ftol,
epsrel, problem);
/* test finite difference Jacobian
* XXX: watson problem doesn't work with forward differences */
if (problem != &watson_problem)
{
fdf.df = problem->fdf->df;
problem->fdf->df = NULL;
test_fdf(gsl_multifit_nlinear_trust, params, xtol, gtol, ftol,
1.0e3 * epsrel, problem);
problem->fdf->df = fdf.df;
}
#if 0 /*XXX: box3d test fails on MacOS here */
if (params->trs == gsl_multifit_nlinear_trs_lmaccel && problem->fdf->fvv != NULL)
{
/* test finite difference second directional derivative */
fdf.fvv = problem->fdf->fvv;
problem->fdf->fvv = NULL;
test_fdf(gsl_multifit_nlinear_trust, params, xtol, gtol, ftol,
epsrel / params->h_fvv, problem);
problem->fdf->fvv = fdf.fvv;
}
#endif
}
}
/*
test_fdf()
Test a weighted nonlinear least squares problem
Inputs: T - solver to use
params - solver parameters
xtol - tolerance in x
gtol - tolerance in gradient
ftol - tolerance in residual vector
epsrel - relative error tolerance in solution
problem - contains the nonlinear problem and solution point
*/
static void
test_fdf(const gsl_multifit_nlinear_type * T,
const gsl_multifit_nlinear_parameters * params,
const double xtol, const double gtol, const double ftol,
const double epsrel,
test_fdf_problem *problem)
{
gsl_multifit_nlinear_fdf *fdf = problem->fdf;
const size_t n = fdf->n;
const size_t p = fdf->p;
const size_t max_iter = 2500;
gsl_vector *x0 = gsl_vector_alloc(p);
gsl_vector_view x0v = gsl_vector_view_array(problem->x0, p);
gsl_multifit_nlinear_workspace *w =
gsl_multifit_nlinear_alloc (T, params, n, p);
const char *pname = problem->name;
char buf[2048];
char sname[2048];
int status, info;
sprintf(buf, "%s/%s/scale=%s/solver=%s",
gsl_multifit_nlinear_name(w),
params->trs->name,
params->scale->name,
params->solver->name);
if (problem->fdf->df == NULL)
{
if (params->fdtype == GSL_MULTIFIT_NLINEAR_FWDIFF)
strcat(buf, "/fdjac,forward");
else
strcat(buf, "/fdjac,center");
}
if (problem->fdf->fvv == NULL)
{
strcat(buf, "/fdfvv");
}
strcpy(sname, buf);
gsl_vector_memcpy(x0, &x0v.vector);
if (problem->weights != NULL)
{
gsl_vector_const_view wv = gsl_vector_const_view_array(problem->weights, n);
gsl_multifit_nlinear_winit(x0, &wv.vector, fdf, w);
}
else
gsl_multifit_nlinear_init(x0, fdf, w);
status = gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
NULL, NULL, &info, w);
gsl_test(status, "%s/%s did not converge, status=%s",
sname, pname, gsl_strerror(status));
/* check solution */
test_fdf_checksol(sname, pname, epsrel, w, problem);
if (problem->weights == NULL)
{
/* test again with weighting matrix W = I */
gsl_vector *wv = gsl_vector_alloc(n);
sprintf(sname, "%s/weighted", buf);
gsl_vector_memcpy(x0, &x0v.vector);
gsl_vector_set_all(wv, 1.0);
gsl_multifit_nlinear_winit(x0, wv, fdf, w);
status = gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
NULL, NULL, &info, w);
gsl_test(status, "%s/%s did not converge, status=%s",
sname, pname, gsl_strerror(status));
test_fdf_checksol(sname, pname, epsrel, w, problem);
gsl_vector_free(wv);
}
gsl_multifit_nlinear_free(w);
gsl_vector_free(x0);
}
static void
test_fdf_checksol(const char *sname, const char *pname,
const double epsrel,
gsl_multifit_nlinear_workspace *w,
test_fdf_problem *problem)
{
gsl_multifit_nlinear_fdf *fdf = problem->fdf;
const double *sigma = problem->sigma;
gsl_vector *f = gsl_multifit_nlinear_residual(w);
gsl_vector *x = gsl_multifit_nlinear_position(w);
double sumsq;
/* check solution vector x and sumsq = ||f||^2 */
gsl_blas_ddot(f, f, &sumsq);
(problem->checksol)(x->data, sumsq, epsrel, sname, pname);
/* check variances */
if (sigma)
{
const size_t n = fdf->n;
const size_t p = fdf->p;
size_t i;
gsl_matrix * covar = gsl_matrix_alloc (p, p);
gsl_matrix *J = gsl_multifit_nlinear_jac (w);
gsl_multifit_nlinear_covar (J, 0.0, covar);
for (i = 0; i < p; i++)
{
double ei = sqrt(sumsq/(n-p))*sqrt(gsl_matrix_get(covar,i,i));
gsl_test_rel (ei, sigma[i], epsrel,
"%s/%s, sigma(%d)", sname, pname, i) ;
}
gsl_matrix_free (covar);
}
}