htswanalysis/MACS/lib/gsl/gsl-1.11/multifit/test.c

   1 /* multifit/test.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* These tests are based on the NIST Statistical Reference Datasets
  21    See http://www.nist.gov/itl/div898/strd/index.html for more
  22    information. */
  23
  24 #include <config.h>
  25 #include <stdlib.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_test.h>
  28 #include <gsl/gsl_multifit.h>
  29 #include <gsl/gsl_multifit_nlin.h>
  30 #include <gsl/gsl_blas.h>
  31
  32 #include <gsl/gsl_ieee_utils.h>
  33
  34 #include "test_longley.c"
  35 #include "test_filip.c"
  36 #include "test_pontius.c"
  37 #include "test_brown.c"
  38 #include "test_enso.c"
  39 #include "test_kirby2.c"
  40 #include "test_hahn1.c"
  41 #include "test_nelson.c"
  42 #include "test_fn.c"
  43 #include "test_estimator.c"
  44
  45 void
  46 test_lmder (gsl_multifit_function_fdf * f, double x0[],
  47             double * X, double F[], double * cov);
  48
  49 void
  50 test_fdf (const char * name, gsl_multifit_function_fdf * f,
  51           double x0[], double x[], double sumsq,
  52           double sigma[]);
  53
  54 int
  55 main (void)
  56 {
  57   gsl_ieee_env_setup();
  58
  59   test_longley();
  60   test_filip();
  61   test_pontius();
  62   test_estimator();
  63
  64   {
  65     gsl_multifit_function_fdf f = make_fdf (&brown_f, &brown_df, &brown_fdf,
  66                                             brown_N, brown_P, 0);
  67
  68     test_lmder(&f, brown_x0, &brown_X[0][0], brown_F, &brown_cov[0][0]);
  69   }
  70
  71   {
  72     gsl_multifit_function_fdf f = make_fdf (&enso_f, &enso_df, &enso_fdf,
  73                                             enso_N, enso_P, 0);
  74
  75     test_fdf("nist-ENSO", &f, enso_x0, enso_x, enso_sumsq, enso_sigma);
  76   }
  77
  78   {
  79     gsl_multifit_function_fdf f = make_fdf (&kirby2_f, &kirby2_df, &kirby2_fdf,
  80                                             kirby2_N, kirby2_P, 0);
  81
  82     test_fdf("nist-kirby2", &f, kirby2_x0, kirby2_x, kirby2_sumsq, kirby2_sigma);
  83   }
  84
  85   {
  86     gsl_multifit_function_fdf f = make_fdf (&hahn1_f, &hahn1_df, &hahn1_fdf,
  87                                             hahn1_N, hahn1_P, 0);
  88
  89     test_fdf("nist-hahn1", &f, hahn1_x0, hahn1_x, hahn1_sumsq, hahn1_sigma);
  90   }
  91
  92 #ifdef JUNK
  93   {
  94     gsl_multifit_function_fdf f = make_fdf (&nelson_f, &nelson_df, &nelson_fdf,
  95                                             nelson_N, nelson_P, 0);
  96
  97     test_fdf("nist-nelson", &f, nelson_x0, nelson_x, nelson_sumsq, nelson_sigma);
  98   }
  99 #endif
 100
 101   /* now summarize the results */
 102
 103   exit (gsl_test_summary ());
 104 }
 105
 106
 107 void
 108 test_lmder (gsl_multifit_function_fdf * f, double x0[],
 109             double * X, double F[], double * cov)
 110 {
 111   const gsl_multifit_fdfsolver_type *T;
 112   gsl_multifit_fdfsolver *s;
 113
 114   const size_t n = f->n;
 115   const size_t p = f->p;
 116
 117   int status;
 118   size_t iter = 0, i;
 119
 120   gsl_vector_view x = gsl_vector_view_array (x0, p);
 121
 122   T = gsl_multifit_fdfsolver_lmsder;
 123   s = gsl_multifit_fdfsolver_alloc (T, n, p);
 124   gsl_multifit_fdfsolver_set (s, f, &x.vector);
 125
 126   do
 127     {
 128       status = gsl_multifit_fdfsolver_iterate (s);
 129
 130       for (i = 0 ; i < p; i++)
 131         {
 132           gsl_test_rel (gsl_vector_get (s->x, i), X[p*iter+i], 1e-5,
 133                         "lmsder, iter=%u, x%u", iter, i);
 134         }
 135
 136       gsl_test_rel (gsl_blas_dnrm2 (s->f), F[iter], 1e-5,
 137                     "lmsder, iter=%u, f", iter);
 138
 139       iter++;
 140     }
 141   while (iter < 20);
 142
 143   {
 144     size_t i, j;
 145     gsl_matrix * covar = gsl_matrix_alloc (4, 4);
 146     gsl_multifit_covar (s->J, 0.0, covar);
 147
 148     for (i = 0; i < 4; i++)
 149       {
 150         for (j = 0; j < 4; j++)
 151           {
 152             gsl_test_rel (gsl_matrix_get(covar,i,j), cov[i*p + j], 1e-7,
 153                           "gsl_multifit_covar cov(%d,%d)", i, j) ;
 154           }
 155       }
 156
 157     gsl_matrix_free (covar);
 158   }
 159
 160   gsl_multifit_fdfsolver_free (s);
 161
 162 }
 163
 164 void
 165 test_fdf (const char * name, gsl_multifit_function_fdf * f,
 166           double x0[], double x_final[],
 167           double f_sumsq, double sigma[])
 168 {
 169   const gsl_multifit_fdfsolver_type *T;
 170   gsl_multifit_fdfsolver *s;
 171
 172   const size_t n = f->n;
 173   const size_t p = f->p;
 174
 175   int status;
 176   size_t iter = 0;
 177
 178   gsl_vector_view x = gsl_vector_view_array (x0, p);
 179
 180   T = gsl_multifit_fdfsolver_lmsder;
 181   s = gsl_multifit_fdfsolver_alloc (T, n, p);
 182   gsl_multifit_fdfsolver_set (s, f, &x.vector);
 183
 184   do
 185     {
 186       status = gsl_multifit_fdfsolver_iterate (s);
 187
 188 #ifdef DEBUG
 189        printf("iter = %d  status = %d  |f| = %.18e x = \n",
 190          iter, status, gsl_blas_dnrm2 (s->f));
 191
 192          gsl_vector_fprintf(stdout, s->x, "%.8e");
 193 #endif
 194       status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-7);
 195
 196       iter++;
 197     }
 198   while (status == GSL_CONTINUE && iter < 1000);
 199
 200   {
 201     size_t i;
 202     gsl_matrix * covar = gsl_matrix_alloc (p, p);
 203     gsl_multifit_covar (s->J, 0.0, covar);
 204
 205     for (i = 0 ; i < p; i++)
 206       {
 207         gsl_test_rel (gsl_vector_get (s->x, i), x_final[i], 1e-5,
 208                       "%s, lmsder, x%u", name, i);
 209       }
 210
 211
 212     {
 213       double s2 = pow(gsl_blas_dnrm2 (s->f), 2.0);
 214
 215       gsl_test_rel (s2, f_sumsq, 1e-5, "%s, lmsder, |f|^2", name);
 216
 217       for (i = 0; i < p; i++)
 218         {
 219           double ei = sqrt(s2/(n-p))*sqrt(gsl_matrix_get(covar,i,i));
 220           gsl_test_rel (ei, sigma[i], 1e-4,
 221                         "%s, sigma(%d)", name, i) ;
 222         }
 223     }
 224
 225     gsl_matrix_free (covar);
 226   }
 227
 228   /* Check that there is no hidden state, restarting should
 229      produce identical results. */
 230
 231   {
 232     int status0, status1;
 233     size_t i;
 234     gsl_multifit_fdfsolver *t = gsl_multifit_fdfsolver_alloc (T, n, p);
 235     gsl_multifit_fdfsolver_set (t, f, &x.vector);
 236
 237     /* do a few extra iterations to stir things up */
 238
 239     gsl_multifit_fdfsolver_set (s, f, &x.vector);
 240
 241     for (i = 0; i < 3; i++)
 242       {
 243         gsl_multifit_fdfsolver_iterate (s);
 244       }
 245
 246     gsl_multifit_fdfsolver_set (s, f, &x.vector);
 247
 248     do
 249       {
 250         status0 = gsl_multifit_fdfsolver_iterate (s);
 251         status1 = gsl_multifit_fdfsolver_iterate (t);
 252
 253         gsl_test_int(status0, status1, "%s, lmsder status after set iter=%u", name, iter);
 254
 255         for (i = 0; i < p; i++) {
 256           double sxi = gsl_vector_get(s->x,i);
 257           double txi = gsl_vector_get(t->x,i);
 258 #ifdef DEBUG
 259           printf("%d %g %g\n", i, sxi, txi);
 260 #endif
 261           gsl_test_rel(sxi, txi, 1e-15, "%s, lmsder after set, %u/%u", name, iter, i);
 262         }
 263
 264 #ifdef DEBUG
 265         printf("iter = %d  status = %d  |f| = %.18e x = \n",
 266                iter, status, gsl_blas_dnrm2 (s->f));
 267
 268         gsl_vector_fprintf(stdout, s->x, "%.8e");
 269 #endif
 270         status0 = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-7);
 271         status1 = gsl_multifit_test_delta (t->dx, s->x, 0.0, 1e-7);
 272
 273         gsl_test_int(status0, status1, "%s, lmsder test delta status after set iter=%u", name, iter);
 274
 275         iter++;
 276       }
 277     while (status1 == GSL_CONTINUE && iter < 1000);
 278
 279     gsl_multifit_fdfsolver_free (t);
 280   }
 281
 282   gsl_multifit_fdfsolver_free (s);
 283 }