Blame multifit/test_nonlinear.c

Packit 67cb25
/* multifit/test_nonlinear.c
Packit 67cb25
 * 
Packit 67cb25
 * Copyright (C) 2007, 2013, 2014 Brian Gough, Patrick Alken
Packit 67cb25
 * 
Packit 67cb25
 * This program is free software; you can redistribute it and/or modify
Packit 67cb25
 * it under the terms of the GNU General Public License as published by
Packit 67cb25
 * the Free Software Foundation; either version 3 of the License, or (at
Packit 67cb25
 * your option) any later version.
Packit 67cb25
 * 
Packit 67cb25
 * This program is distributed in the hope that it will be useful, but
Packit 67cb25
 * WITHOUT ANY WARRANTY; without even the implied warranty of
Packit 67cb25
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Packit 67cb25
 * General Public License for more details.
Packit 67cb25
 * 
Packit 67cb25
 * You should have received a copy of the GNU General Public License
Packit 67cb25
 * along with this program; if not, write to the Free Software
Packit 67cb25
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
Packit 67cb25
 */
Packit 67cb25
Packit 67cb25
typedef struct
Packit 67cb25
{
Packit 67cb25
  const char *name;
Packit 67cb25
  double *x0;       /* initial parameters (size p) */
Packit 67cb25
  double *sigma;
Packit 67cb25
  double *epsrel;   /* relative tolerance for solution checking */
Packit 67cb25
  size_t ntries;
Packit 67cb25
  void (*checksol) (const double x[], const double sumsq,
Packit 67cb25
                    const double epsrel, const char *sname,
Packit 67cb25
                    const char *pname);
Packit 67cb25
  gsl_multifit_function_fdf *fdf;
Packit 67cb25
} test_fdf_problem;
Packit 67cb25
Packit 67cb25
#include "test_bard.c"
Packit 67cb25
#include "test_beale.c"
Packit 67cb25
#include "test_biggs.c"
Packit 67cb25
#include "test_box.c"
Packit 67cb25
#include "test_boxbod.c"
Packit 67cb25
#include "test_brown1.c"
Packit 67cb25
#include "test_brown2.c"
Packit 67cb25
#include "test_brown3.c"
Packit 67cb25
#include "test_eckerle.c"
Packit 67cb25
#include "test_enso.c"
Packit 67cb25
#include "test_exp1.c"
Packit 67cb25
#include "test_gaussian.c"
Packit 67cb25
#include "test_hahn1.c"
Packit 67cb25
#include "test_helical.c"
Packit 67cb25
#include "test_jennrich.c"
Packit 67cb25
#include "test_kirby2.c"
Packit 67cb25
#include "test_kowalik.c"
Packit 67cb25
#include "test_lin1.c"
Packit 67cb25
#include "test_lin2.c"
Packit 67cb25
#include "test_lin3.c"
Packit 67cb25
#include "test_meyer.c"
Packit 67cb25
#include "test_meyerscal.c"
Packit 67cb25
#include "test_osborne.c"
Packit 67cb25
#include "test_penalty1.c"
Packit 67cb25
#include "test_penalty2.c"
Packit 67cb25
#include "test_powell1.c"
Packit 67cb25
#include "test_powell2.c"
Packit 67cb25
#include "test_powell3.c"
Packit 67cb25
#include "test_rat42.c"
Packit 67cb25
#include "test_rat43.c"
Packit 67cb25
#include "test_rosenbrock.c"
Packit 67cb25
#include "test_rosenbrocke.c"
Packit 67cb25
#include "test_roth.c"
Packit 67cb25
#include "test_thurber.c"
Packit 67cb25
#include "test_vardim.c"
Packit 67cb25
#include "test_watson.c"
Packit 67cb25
#include "test_wood.c"
Packit 67cb25
Packit 67cb25
#include "test_wnlin.c"
Packit 67cb25
Packit 67cb25
static void test_fdf(const gsl_multifit_fdfsolver_type * T,
Packit 67cb25
                     const double xtol, const double gtol,
Packit 67cb25
                     const double ftol, const double epsrel,
Packit 67cb25
                     const double x0_scale, test_fdf_problem *problem,
Packit 67cb25
                     const double *wts);
Packit 67cb25
static void test_fdfridge(const gsl_multifit_fdfsolver_type * T,
Packit 67cb25
                          const double xtol, const double gtol,
Packit 67cb25
                          const double ftol, const double epsrel,
Packit 67cb25
                          const double x0_scale, test_fdf_problem *problem,
Packit 67cb25
                          const double *wts);
Packit 67cb25
static void test_fdf_checksol(const char *sname, const char *pname,
Packit 67cb25
                              const double epsrel,
Packit 67cb25
                              gsl_multifit_fdfsolver *s,
Packit 67cb25
                              test_fdf_problem *problem);
Packit 67cb25
static void test_scale_x0(gsl_vector *x0, const double scale);
Packit 67cb25
Packit 67cb25
/*
Packit 67cb25
 * These test problems are taken from
Packit 67cb25
 *
Packit 67cb25
 * H. B. Nielsen, UCTP test problems for unconstrained optimization,
Packit 67cb25
 * IMM Department of Mathematical Modeling, Tech. Report IMM-REP-2000-17,
Packit 67cb25
 * 2000.
Packit 67cb25
 */
Packit 67cb25
static test_fdf_problem *test_fdf_nielsen[] = {
Packit 67cb25
  &lin1_problem,       /* 1 */
Packit 67cb25
  &lin2_problem,       /* 2 */
Packit 67cb25
  &lin3_problem,       /* 3 */
Packit 67cb25
  &rosenbrock_problem, /* 4 */
Packit 67cb25
  &helical_problem,    /* 5 */
Packit 67cb25
  &powell1_problem,    /* 6 */
Packit 67cb25
  &roth_problem,       /* 7 */
Packit 67cb25
  &bard_problem,       /* 8 */
Packit 67cb25
  &kowalik_problem,    /* 9 */
Packit 67cb25
  &meyer_problem,      /* 10 */
Packit 67cb25
  &watson_problem,     /* 11 */
Packit 67cb25
  &box_problem,        /* 12 */
Packit 67cb25
  &jennrich_problem,   /* 13 */
Packit 67cb25
  &brown1_problem,     /* 14 */
Packit 67cb25
  &brown2_problem,     /* 16 */
Packit 67cb25
  &osborne_problem,    /* 17 */
Packit 67cb25
  &exp1_problem,       /* 18 */
Packit 67cb25
  &meyerscal_problem,  /* 20 */
Packit 67cb25
Packit 67cb25
  &powell2_problem,
Packit 67cb25
Packit 67cb25
  NULL
Packit 67cb25
};
Packit 67cb25
Packit 67cb25
/*
Packit 67cb25
 * These tests are from
Packit 67cb25
 *
Packit 67cb25
 * J. J. More, B. S. Garbow and K. E. Hillstrom, Testing
Packit 67cb25
 * Unconstrained Optimization Software, ACM Trans. Math. Soft.
Packit 67cb25
 * Vol 7, No 1, 1981.
Packit 67cb25
 *
Packit 67cb25
 * Many of these overlap with the Nielsen tests
Packit 67cb25
 */
Packit 67cb25
static test_fdf_problem *test_fdf_more[] = {
Packit 67cb25
  &rosenbrock_problem,   /* 1 */
Packit 67cb25
  &roth_problem,         /* 2 */
Packit 67cb25
  &powell3_problem,      /* 3 */
Packit 67cb25
  &brown3_problem,       /* 4 */
Packit 67cb25
  &beale_problem,        /* 5 */
Packit 67cb25
  &jennrich_problem,     /* 6 */
Packit 67cb25
  &helical_problem,      /* 7 */
Packit 67cb25
  &bard_problem,         /* 8 */
Packit 67cb25
  &gaussian_problem,     /* 9 */
Packit 67cb25
  &meyer_problem,        /* 10 */
Packit 67cb25
  &box_problem,          /* 12 */
Packit 67cb25
  &powell1_problem,      /* 13 */
Packit 67cb25
  &wood_problem,         /* 14 */
Packit 67cb25
  &kowalik_problem,      /* 15 */
Packit 67cb25
  &brown1_problem,       /* 16 */
Packit 67cb25
  &osborne_problem,      /* 17 */
Packit 67cb25
  &biggs_problem,        /* 18 */
Packit 67cb25
  &watson_problem,       /* 20 */
Packit 67cb25
  &rosenbrocke_problem,  /* 21 */
Packit 67cb25
  &penalty1_problem,     /* 23 */
Packit 67cb25
  &penalty2_problem,     /* 24 */
Packit 67cb25
  &vardim_problem,       /* 25 */
Packit 67cb25
  &brown2_problem,       /* 27 */
Packit 67cb25
  &lin1_problem,         /* 32 */
Packit 67cb25
  &lin2_problem,         /* 33 */
Packit 67cb25
  &lin3_problem,         /* 34 */
Packit 67cb25
Packit 67cb25
  NULL
Packit 67cb25
};
Packit 67cb25
Packit 67cb25
/* NIST test cases */
Packit 67cb25
static test_fdf_problem *test_fdf_nist[] = {
Packit 67cb25
  &kirby2_problem,
Packit 67cb25
  &hahn1_problem,
Packit 67cb25
  &enso_problem,
Packit 67cb25
  &thurber_problem,
Packit 67cb25
  &boxbod_problem,
Packit 67cb25
  &rat42_problem,
Packit 67cb25
  &eckerle_problem,
Packit 67cb25
  &rat43_problem,
Packit 67cb25
Packit 67cb25
  NULL
Packit 67cb25
};
Packit 67cb25
Packit 67cb25
static void
Packit 67cb25
test_nonlinear(void)
Packit 67cb25
{
Packit f2184d
  const double xtol = pow(GSL_DBL_EPSILON, 0.9) * 10.0;
Packit f2184d
  const double gtol = pow(GSL_DBL_EPSILON, 0.9) * 10.0;
Packit 67cb25
  const double ftol = 0.0;
Packit 67cb25
  size_t i, j;
Packit 67cb25
Packit 67cb25
  /* Nielsen tests */
Packit 67cb25
  for (i = 0; test_fdf_nielsen[i] != NULL; ++i)
Packit 67cb25
    {
Packit 67cb25
      test_fdf_problem *problem = test_fdf_nielsen[i];
Packit 67cb25
      double epsrel = *(problem->epsrel);
Packit 67cb25
      double scale = 1.0;
Packit 67cb25
Packit 67cb25
      for (j = 0; j < problem->ntries; ++j)
Packit 67cb25
        {
Packit 67cb25
          double eps_scale = epsrel * scale;
Packit 67cb25
Packit 67cb25
          test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                   eps_scale, scale, problem, NULL);
Packit 67cb25
          test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                        eps_scale, scale, problem, NULL);
Packit 67cb25
Packit 67cb25
          /* test finite difference Jacobian */
Packit 67cb25
          {
Packit 67cb25
            gsl_multifit_function_fdf fdf;
Packit 67cb25
            fdf.df = problem->fdf->df;
Packit 67cb25
            problem->fdf->df = NULL;
Packit 67cb25
            test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                     1.0e5 * eps_scale, 1.0, problem, NULL);
Packit 67cb25
            test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                          1.0e5 * eps_scale, 1.0, problem, NULL);
Packit 67cb25
            problem->fdf->df = fdf.df;
Packit 67cb25
          }
Packit 67cb25
Packit 67cb25
          scale *= 10.0;
Packit 67cb25
        }
Packit 67cb25
Packit 67cb25
      test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
Packit 67cb25
               10.0 * epsrel, 1.0, problem, NULL);
Packit 67cb25
    }
Packit 67cb25
Packit 67cb25
  /* More tests */
Packit 67cb25
  for (i = 0; test_fdf_more[i] != NULL; ++i)
Packit 67cb25
    {
Packit 67cb25
      test_fdf_problem *problem = test_fdf_more[i];
Packit 67cb25
      double epsrel = *(problem->epsrel);
Packit 67cb25
      double scale = 1.0;
Packit 67cb25
Packit 67cb25
      for (j = 0; j < problem->ntries; ++j)
Packit 67cb25
        {
Packit 67cb25
          double eps_scale = epsrel * scale;
Packit 67cb25
Packit 67cb25
          test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                   eps_scale, scale, problem, NULL);
Packit 67cb25
          test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                        eps_scale, scale, problem, NULL);
Packit 67cb25
Packit 67cb25
          /* test finite difference Jacobian */
Packit 67cb25
          {
Packit 67cb25
            gsl_multifit_function_fdf fdf;
Packit 67cb25
            fdf.df = problem->fdf->df;
Packit 67cb25
            problem->fdf->df = NULL;
Packit 67cb25
            test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                     1.0e5 * eps_scale, 1.0, problem, NULL);
Packit 67cb25
            test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                          1.0e5 * eps_scale, 1.0, problem, NULL);
Packit 67cb25
            problem->fdf->df = fdf.df;
Packit 67cb25
          }
Packit 67cb25
Packit 67cb25
          scale *= 10.0;
Packit 67cb25
        }
Packit 67cb25
Packit 67cb25
      test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
Packit 67cb25
               10.0 * epsrel, 1.0, problem, NULL);
Packit 67cb25
    }
Packit 67cb25
Packit 67cb25
  /* NIST tests */
Packit 67cb25
  for (i = 0; test_fdf_nist[i] != NULL; ++i)
Packit 67cb25
    {
Packit 67cb25
      test_fdf_problem *problem = test_fdf_nist[i];
Packit 67cb25
      double epsrel = *(problem->epsrel);
Packit 67cb25
      double scale = 1.0;
Packit 67cb25
Packit 67cb25
      for (j = 0; j < problem->ntries; ++j)
Packit 67cb25
        {
Packit 67cb25
          double eps_scale = epsrel * scale;
Packit 67cb25
Packit 67cb25
          test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                   eps_scale, scale, problem, NULL);
Packit 67cb25
          test_fdf(gsl_multifit_fdfsolver_lmder, xtol, gtol, ftol,
Packit 67cb25
                   eps_scale, scale, problem, NULL);
Packit 67cb25
Packit 67cb25
          /* test finite difference Jacobian */
Packit 67cb25
          {
Packit 67cb25
            gsl_multifit_function_fdf fdf;
Packit 67cb25
            fdf.df = problem->fdf->df;
Packit 67cb25
            problem->fdf->df = NULL;
Packit 67cb25
            test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                     eps_scale, 1.0, problem, NULL);
Packit 67cb25
            test_fdf(gsl_multifit_fdfsolver_lmder, xtol, gtol, ftol,
Packit 67cb25
                     eps_scale, scale, problem, NULL);
Packit 67cb25
            problem->fdf->df = fdf.df;
Packit 67cb25
          }
Packit 67cb25
Packit 67cb25
          scale *= 10.0;
Packit 67cb25
        }
Packit 67cb25
Packit 67cb25
      test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
Packit 67cb25
               epsrel, 1.0, problem, NULL);
Packit 67cb25
    }
Packit 67cb25
Packit 67cb25
  /* test weighted nonlinear least squares */
Packit 67cb25
Packit 67cb25
  /* internal weighting in _f and _df functions */
Packit 67cb25
  test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
           wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
Packit 67cb25
  test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
Packit 67cb25
           wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
Packit 67cb25
  test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
Packit 67cb25
  test_fdfridge(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
Packit 67cb25
                wnlin_epsrel, 1.0, &wnlin_problem1, NULL);
Packit 67cb25
Packit 67cb25
  /* weighting through fdfsolver_wset */
Packit 67cb25
  test_fdf(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
           wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
Packit 67cb25
  test_fdf(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
Packit 67cb25
           wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
Packit 67cb25
  test_fdfridge(gsl_multifit_fdfsolver_lmsder, xtol, gtol, ftol,
Packit 67cb25
                wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
Packit 67cb25
  test_fdfridge(gsl_multifit_fdfsolver_lmniel, xtol, gtol, ftol,
Packit 67cb25
                wnlin_epsrel, 1.0, &wnlin_problem2, wnlin_W);
Packit 67cb25
}
Packit 67cb25
Packit 67cb25
/*
Packit 67cb25
test_fdf()
Packit 67cb25
  Test a weighted nonlinear least squares problem
Packit 67cb25
Packit 67cb25
Inputs: T        - solver to use
Packit 67cb25
        xtol     - tolerance in x
Packit 67cb25
        gtol     - tolerance in gradient
Packit 67cb25
        ftol     - tolerance in residual vector
Packit 67cb25
        epsrel   - relative error tolerance in solution
Packit 67cb25
        x0_scale - to test robustness against starting points,
Packit 67cb25
                   the standard starting point in 'problem' is
Packit 67cb25
                   multiplied by this scale factor:
Packit 67cb25
                   x0 <- x0 * x0_scale
Packit 67cb25
                   If x0 = 0, then all components of x0 are set to
Packit 67cb25
                   x0_scale
Packit 67cb25
        problem  - contains the nonlinear problem and solution point
Packit 67cb25
        wts      - weight vector (NULL for unweighted)
Packit 67cb25
*/
Packit 67cb25
Packit 67cb25
static void
Packit 67cb25
test_fdf(const gsl_multifit_fdfsolver_type * T, const double xtol,
Packit 67cb25
         const double gtol, const double ftol,
Packit 67cb25
         const double epsrel, const double x0_scale,
Packit 67cb25
         test_fdf_problem *problem,
Packit 67cb25
         const double *wts)
Packit 67cb25
{
Packit 67cb25
  gsl_multifit_function_fdf *fdf = problem->fdf;
Packit 67cb25
  const size_t n = fdf->n;
Packit 67cb25
  const size_t p = fdf->p;
Packit 67cb25
  const size_t max_iter = 1500;
Packit 67cb25
  gsl_vector *x0 = gsl_vector_alloc(p);
Packit 67cb25
  gsl_vector_view x0v = gsl_vector_view_array(problem->x0, p);
Packit 67cb25
  gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc (T, n, p);
Packit 67cb25
  const char *pname = problem->name;
Packit 67cb25
  char sname[2048];
Packit 67cb25
  int status, info;
Packit 67cb25
Packit 67cb25
  sprintf(sname, "%s/scale=%g%s",
Packit 67cb25
    gsl_multifit_fdfsolver_name(s), x0_scale,
Packit 67cb25
    problem->fdf->df ? "" : "/fdiff");
Packit 67cb25
Packit 67cb25
  /* scale starting point x0 */
Packit 67cb25
  gsl_vector_memcpy(x0, &x0v.vector);
Packit 67cb25
  test_scale_x0(x0, x0_scale);
Packit 67cb25
Packit 67cb25
  if (wts)
Packit 67cb25
    {
Packit 67cb25
      gsl_vector_const_view wv = gsl_vector_const_view_array(wts, n);
Packit 67cb25
      gsl_multifit_fdfsolver_wset(s, fdf, x0, &wv.vector);
Packit 67cb25
    }
Packit 67cb25
  else
Packit 67cb25
    gsl_multifit_fdfsolver_set(s, fdf, x0);
Packit 67cb25
Packit 67cb25
  status = gsl_multifit_fdfsolver_driver(s, max_iter, xtol, gtol,
Packit 67cb25
                                         ftol, &info;;
Packit 67cb25
  gsl_test(status, "%s/%s did not converge, status=%s",
Packit 67cb25
           sname, pname, gsl_strerror(status));
Packit 67cb25
Packit 67cb25
  /* check solution */
Packit 67cb25
  test_fdf_checksol(sname, pname, epsrel, s, problem);
Packit 67cb25
Packit 67cb25
  if (wts == NULL)
Packit 67cb25
    {
Packit 67cb25
      /* test again with weighting matrix W = I */
Packit 67cb25
      gsl_vector *wv = gsl_vector_alloc(n);
Packit 67cb25
Packit 67cb25
      sprintf(sname, "%s/scale=%g%s/weights",
Packit 67cb25
        gsl_multifit_fdfsolver_name(s), x0_scale,
Packit 67cb25
        problem->fdf->df ? "" : "/fdiff");
Packit 67cb25
Packit 67cb25
      gsl_vector_memcpy(x0, &x0v.vector);
Packit 67cb25
      test_scale_x0(x0, x0_scale);
Packit 67cb25
Packit 67cb25
      gsl_vector_set_all(wv, 1.0);
Packit 67cb25
      gsl_multifit_fdfsolver_wset(s, fdf, x0, wv);
Packit 67cb25
  
Packit 67cb25
      status = gsl_multifit_fdfsolver_driver(s, max_iter, xtol, gtol,
Packit 67cb25
                                             ftol, &info;;
Packit 67cb25
      gsl_test(status, "%s/%s did not converge, status=%s",
Packit 67cb25
               sname, pname, gsl_strerror(status));
Packit 67cb25
Packit 67cb25
      test_fdf_checksol(sname, pname, epsrel, s, problem);
Packit 67cb25
Packit 67cb25
      gsl_vector_free(wv);
Packit 67cb25
    }
Packit 67cb25
Packit 67cb25
  gsl_multifit_fdfsolver_free(s);
Packit 67cb25
  gsl_vector_free(x0);
Packit 67cb25
}
Packit 67cb25
Packit 67cb25
/*
Packit 67cb25
test_fdfridge()
Packit 67cb25
  Test a nonlinear least squares problem
Packit 67cb25
Packit 67cb25
Inputs: T        - solver to use
Packit 67cb25
        xtol     - tolerance in x
Packit 67cb25
        gtol     - tolerance in gradient
Packit 67cb25
        ftol     - tolerance in residual vector
Packit 67cb25
        epsrel   - relative error tolerance in solution
Packit 67cb25
        x0_scale - to test robustness against starting points,
Packit 67cb25
                   the standard starting point in 'problem' is
Packit 67cb25
                   multiplied by this scale factor:
Packit 67cb25
                   x0 <- x0 * x0_scale
Packit 67cb25
                   If x0 = 0, then all components of x0 are set to
Packit 67cb25
                   x0_scale
Packit 67cb25
        problem  - contains the nonlinear problem and solution point
Packit 67cb25
        wts      - weight vector
Packit 67cb25
*/
Packit 67cb25
Packit 67cb25
static void
Packit 67cb25
test_fdfridge(const gsl_multifit_fdfsolver_type * T, const double xtol,
Packit 67cb25
              const double gtol, const double ftol,
Packit 67cb25
              const double epsrel, const double x0_scale,
Packit 67cb25
              test_fdf_problem *problem, const double *wts)
Packit 67cb25
{
Packit 67cb25
  gsl_multifit_function_fdf *fdf = problem->fdf;
Packit 67cb25
  const size_t n = fdf->n;
Packit 67cb25
  const size_t p = fdf->p;
Packit 67cb25
  const size_t max_iter = 1500;
Packit 67cb25
  gsl_vector *x0 = gsl_vector_alloc(p);
Packit 67cb25
  gsl_vector_view x0v = gsl_vector_view_array(problem->x0, p);
Packit 67cb25
  gsl_multifit_fdfridge *w = gsl_multifit_fdfridge_alloc (T, n, p);
Packit 67cb25
  const char *pname = problem->name;
Packit 67cb25
  char sname[2048];
Packit 67cb25
  int status, info;
Packit 67cb25
  double lambda = 0.0;
Packit 67cb25
Packit 67cb25
  sprintf(sname, "ridge/%s", gsl_multifit_fdfridge_name(w));
Packit 67cb25
Packit 67cb25
  /* scale starting point x0 */
Packit 67cb25
  gsl_vector_memcpy(x0, &x0v.vector);
Packit 67cb25
  test_scale_x0(x0, x0_scale);
Packit 67cb25
Packit 67cb25
  /* test undamped case with lambda = 0 */
Packit 67cb25
  if (wts)
Packit 67cb25
    {
Packit 67cb25
      gsl_vector_const_view wv = gsl_vector_const_view_array(wts, n);
Packit 67cb25
      gsl_multifit_fdfridge_wset(w, fdf, x0, lambda, &wv.vector);
Packit 67cb25
    }
Packit 67cb25
  else
Packit 67cb25
    gsl_multifit_fdfridge_set(w, fdf, x0, lambda);
Packit 67cb25
Packit 67cb25
  status = gsl_multifit_fdfridge_driver(w, max_iter, xtol, gtol,
Packit 67cb25
                                        ftol, &info;;
Packit 67cb25
  gsl_test(status, "%s/%s did not converge, status=%s",
Packit 67cb25
           sname, pname, gsl_strerror(status));
Packit 67cb25
Packit 67cb25
  /* check solution */
Packit 67cb25
  test_fdf_checksol(sname, pname, epsrel, w->s, problem);
Packit 67cb25
Packit 67cb25
  /* test for self consisent solution with L = \lambda I */
Packit 67cb25
  {
Packit 67cb25
    const double eps = 1.0e-10;
Packit 67cb25
    gsl_matrix *L = gsl_matrix_calloc(p, p);
Packit 67cb25
    gsl_vector_view diag = gsl_matrix_diagonal(L);
Packit 67cb25
    gsl_multifit_fdfridge *w2 = gsl_multifit_fdfridge_alloc (T, n, p);
Packit 67cb25
    gsl_vector *y0 = gsl_vector_alloc(p);
Packit 67cb25
    size_t i;
Packit 67cb25
Packit 67cb25
    /* pick some value for lambda and set L = \lambda I */
Packit 67cb25
    lambda = 5.0;
Packit 67cb25
    gsl_vector_set_all(&diag.vector, lambda);
Packit 67cb25
Packit 67cb25
    /* scale initial vector */
Packit 67cb25
    gsl_vector_memcpy(x0, &x0v.vector);
Packit 67cb25
    test_scale_x0(x0, x0_scale);
Packit 67cb25
    gsl_vector_memcpy(y0, x0);
Packit 67cb25
Packit 67cb25
    if (wts)
Packit 67cb25
      {
Packit 67cb25
        gsl_vector_const_view wv = gsl_vector_const_view_array(wts, n);
Packit 67cb25
        gsl_multifit_fdfridge_wset(w, fdf, x0, lambda, &wv.vector);
Packit 67cb25
        gsl_multifit_fdfridge_wset3(w2, fdf, y0, L, &wv.vector);
Packit 67cb25
      }
Packit 67cb25
    else
Packit 67cb25
      {
Packit 67cb25
        gsl_multifit_fdfridge_set(w, fdf, x0, lambda);
Packit 67cb25
        gsl_multifit_fdfridge_set3(w2, fdf, y0, L);
Packit 67cb25
      }
Packit 67cb25
Packit 67cb25
    /* solve with scalar lambda routine */
Packit 67cb25
    status = gsl_multifit_fdfridge_driver(w, max_iter, xtol, gtol,
Packit 67cb25
                                          ftol, &info;;
Packit 67cb25
    gsl_test(status, "%s/lambda/%s did not converge, status=%s",
Packit 67cb25
             sname, pname, gsl_strerror(status));
Packit 67cb25
Packit 67cb25
    /* solve with general matrix routine */
Packit 67cb25
    status = gsl_multifit_fdfridge_driver(w2, max_iter, xtol, gtol,
Packit 67cb25
                                          ftol, &info;;
Packit 67cb25
    gsl_test(status, "%s/L/%s did not converge, status=%s",
Packit 67cb25
             sname, pname, gsl_strerror(status));
Packit 67cb25
Packit 67cb25
    /* test x = y */
Packit 67cb25
    for (i = 0; i < p; ++i)
Packit 67cb25
      {
Packit 67cb25
        double xi = gsl_vector_get(w->s->x, i);
Packit 67cb25
        double yi = gsl_vector_get(w2->s->x, i);
Packit 67cb25
Packit 67cb25
        if (fabs(xi) < eps)
Packit 67cb25
          {
Packit 67cb25
            gsl_test_abs(yi, xi, eps, "%s/%s ridge lambda=%g i=%zu",
Packit 67cb25
                         sname, pname, lambda, i);
Packit 67cb25
          }
Packit 67cb25
        else
Packit 67cb25
          {
Packit 67cb25
            gsl_test_rel(yi, xi, eps, "%s/%s ridge lambda=%g i=%zu",
Packit 67cb25
                         sname, pname, lambda, i);
Packit 67cb25
          }
Packit 67cb25
      }
Packit 67cb25
Packit 67cb25
    gsl_matrix_free(L);
Packit 67cb25
    gsl_vector_free(y0);
Packit 67cb25
    gsl_multifit_fdfridge_free(w2);
Packit 67cb25
  }
Packit 67cb25
Packit 67cb25
  gsl_multifit_fdfridge_free(w);
Packit 67cb25
  gsl_vector_free(x0);
Packit 67cb25
}
Packit 67cb25
Packit 67cb25
static void
Packit 67cb25
test_fdf_checksol(const char *sname, const char *pname,
Packit 67cb25
                  const double epsrel, gsl_multifit_fdfsolver *s,
Packit 67cb25
                  test_fdf_problem *problem)
Packit 67cb25
{
Packit 67cb25
  gsl_multifit_function_fdf *fdf = problem->fdf;
Packit 67cb25
  const double *sigma = problem->sigma;
Packit 67cb25
  gsl_vector *f = gsl_multifit_fdfsolver_residual(s);
Packit 67cb25
  gsl_vector *x = gsl_multifit_fdfsolver_position(s);
Packit 67cb25
  double sumsq;
Packit 67cb25
Packit 67cb25
  /* check solution vector x and sumsq = ||f||^2 */
Packit 67cb25
  gsl_blas_ddot(f, f, &sumsq);
Packit 67cb25
  (problem->checksol)(x->data, sumsq, epsrel, sname, pname);
Packit 67cb25
Packit 67cb25
#if 1
Packit 67cb25
  /* check variances */
Packit 67cb25
  if (sigma)
Packit 67cb25
    {
Packit 67cb25
      const size_t n = fdf->n;
Packit 67cb25
      const size_t p = fdf->p;
Packit 67cb25
      size_t i;
Packit 67cb25
      gsl_matrix * J = gsl_matrix_alloc(n, p);
Packit 67cb25
      gsl_matrix * covar = gsl_matrix_alloc (p, p);
Packit 67cb25
Packit 67cb25
      gsl_multifit_fdfsolver_jac (s, J);
Packit 67cb25
      gsl_multifit_covar(J, 0.0, covar);
Packit 67cb25
Packit 67cb25
      for (i = 0; i < p; i++) 
Packit 67cb25
        {
Packit 67cb25
          double ei = sqrt(sumsq/(n-p))*sqrt(gsl_matrix_get(covar,i,i));
Packit 67cb25
          gsl_test_rel (ei, sigma[i], epsrel, 
Packit 67cb25
                        "%s/%s, sigma(%d)", sname, pname, i) ;
Packit 67cb25
        }
Packit 67cb25
Packit 67cb25
      gsl_matrix_free (J);
Packit 67cb25
      gsl_matrix_free (covar);
Packit 67cb25
    }
Packit 67cb25
#endif
Packit 67cb25
}
Packit 67cb25
Packit 67cb25
static void
Packit 67cb25
test_scale_x0(gsl_vector *x0, const double scale)
Packit 67cb25
{
Packit 67cb25
  double nx = gsl_blas_dnrm2(x0);
Packit 67cb25
Packit 67cb25
  if (nx == 0.0)
Packit 67cb25
    gsl_vector_set_all(x0, scale);
Packit 67cb25
  else
Packit 67cb25
    gsl_vector_scale(x0, scale);
Packit 67cb25
} /* test_scale_x0() */