/* 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);

    }

}