util/X86MAC64/include/gsl/gsl_multifit_nlinear.h

   1 /* multifit_nlinear/gsl_multifit_nlinear.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  * Copyright (C) 2015, 2016 Patrick Alken
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 #ifndef __GSL_MULTIFIT_NLINEAR_H__
  22 #define __GSL_MULTIFIT_NLINEAR_H__
  23
  24 #include <stdlib.h>
  25 #include <gsl/gsl_types.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_vector.h>
  28 #include <gsl/gsl_matrix.h>
  29 #include <gsl/gsl_permutation.h>
  30
  31 #undef __BEGIN_DECLS
  32 #undef __END_DECLS
  33 #ifdef __cplusplus
  34 # define __BEGIN_DECLS extern "C" {
  35 # define __END_DECLS }
  36 #else
  37 # define __BEGIN_DECLS /* empty */
  38 # define __END_DECLS /* empty */
  39 #endif
  40
  41 __BEGIN_DECLS
  42
  43 typedef enum
  44 {
  45   GSL_MULTIFIT_NLINEAR_FWDIFF,
  46   GSL_MULTIFIT_NLINEAR_CTRDIFF
  47 } gsl_multifit_nlinear_fdtype;
  48
  49 /* Definition of vector-valued functions and gradient with parameters
  50    based on gsl_vector */
  51
  52 typedef struct
  53 {
  54   int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  55   int (* df) (const gsl_vector * x, void * params, gsl_matrix * df);
  56   int (* fvv) (const gsl_vector * x, const gsl_vector * v, void * params,
  57                gsl_vector * fvv);
  58   size_t n;        /* number of functions */
  59   size_t p;        /* number of independent variables */
  60   void * params;   /* user parameters */
  61   size_t nevalf;   /* number of function evaluations */
  62   size_t nevaldf;  /* number of Jacobian evaluations */
  63   size_t nevalfvv; /* number of fvv evaluations */
  64 } gsl_multifit_nlinear_fdf;
  65
  66 /* trust region subproblem method */
  67 typedef struct
  68 {
  69   const char *name;
  70   void * (*alloc) (const void * params, const size_t n, const size_t p);
  71   int (*init) (const void * vtrust_state, void * vstate);
  72   int (*preloop) (const void * vtrust_state, void * vstate);
  73   int (*step) (const void * vtrust_state, const double delta,
  74                gsl_vector * dx, void * vstate);
  75   int (*preduction) (const void * vtrust_state, const gsl_vector * dx,
  76                      double * pred, void * vstate);
  77   void (*free) (void * vstate);
  78 } gsl_multifit_nlinear_trs;
  79
  80 /* scaling matrix specification */
  81 typedef struct
  82 {
  83   const char *name;
  84   int (*init) (const gsl_matrix * J, gsl_vector * diag);
  85   int (*update) (const gsl_matrix * J, gsl_vector * diag);
  86 } gsl_multifit_nlinear_scale;
  87
  88 /*
  89  * linear least squares solvers - there are three steps to
  90  * solving a least squares problem using a trust region
  91  * method:
  92  *
  93  * 1. init: called once per iteration when a new Jacobian matrix
  94  *          is computed; perform factorization of Jacobian (qr,svd)
  95  *          or form normal equations matrix (cholesky)
  96  * 2. presolve: called each time a new LM parameter value mu is available;
  97  *              used for cholesky method in order to factor
  98  *              the (J^T J + mu D^T D) matrix
  99  * 3. solve: solve the least square system for a given rhs
 100  */
 101 typedef struct
 102 {
 103   const char *name;
 104   void * (*alloc) (const size_t n, const size_t p);
 105   int (*init) (const void * vtrust_state, void * vstate);
 106   int (*presolve) (const double mu, const void * vtrust_state, void * vstate);
 107   int (*solve) (const gsl_vector * f, gsl_vector * x,
 108                 const void * vtrust_state, void * vstate);
 109   int (*rcond) (double * rcond, void * vstate);
 110   void (*free) (void * vstate);
 111 } gsl_multifit_nlinear_solver;
 112
 113 /* tunable parameters */
 114 typedef struct
 115 {
 116   const gsl_multifit_nlinear_trs *trs;        /* trust region subproblem method */
 117   const gsl_multifit_nlinear_scale *scale;    /* scaling method */
 118   const gsl_multifit_nlinear_solver *solver;  /* solver method */
 119   gsl_multifit_nlinear_fdtype fdtype;         /* finite difference method */
 120   double factor_up;                           /* factor for increasing trust radius */
 121   double factor_down;                         /* factor for decreasing trust radius */
 122   double avmax;                               /* max allowed |a|/|v| */
 123   double h_df;                                /* step size for finite difference Jacobian */
 124   double h_fvv;                               /* step size for finite difference fvv */
 125 } gsl_multifit_nlinear_parameters;
 126
 127 typedef struct
 128 {
 129   const char *name;
 130   void * (*alloc) (const gsl_multifit_nlinear_parameters * params,
 131                    const size_t n, const size_t p);
 132   int (*init) (void * state, const gsl_vector * wts,
 133                gsl_multifit_nlinear_fdf * fdf, const gsl_vector * x,
 134                gsl_vector * f, gsl_matrix * J, gsl_vector * g);
 135   int (*iterate) (void * state, const gsl_vector * wts,
 136                   gsl_multifit_nlinear_fdf * fdf, gsl_vector * x,
 137                   gsl_vector * f, gsl_matrix * J, gsl_vector * g,
 138                   gsl_vector * dx);
 139   int (*rcond) (double * rcond, void * state);
 140   double (*avratio) (void * state);
 141   void (*free) (void * state);
 142 } gsl_multifit_nlinear_type;
 143
 144 /* current state passed to low-level trust region algorithms */
 145 typedef struct
 146 {
 147   const gsl_vector * x;             /* parameter values x */
 148   const gsl_vector * f;             /* residual vector f(x) */
 149   const gsl_vector * g;             /* gradient J^T f */
 150   const gsl_matrix * J;             /* Jacobian J(x) */
 151   const gsl_vector * diag;          /* scaling matrix D */
 152   const gsl_vector * sqrt_wts;      /* sqrt(diag(W)) or NULL for unweighted */
 153   const double *mu;                 /* LM parameter */
 154   const gsl_multifit_nlinear_parameters * params;
 155   void *solver_state;               /* workspace for linear least squares solver */
 156   gsl_multifit_nlinear_fdf * fdf;
 157   double *avratio;                  /* |a| / |v| */
 158 } gsl_multifit_nlinear_trust_state;
 159
 160 typedef struct
 161 {
 162   const gsl_multifit_nlinear_type * type;
 163   gsl_multifit_nlinear_fdf * fdf ;
 164   gsl_vector * x;             /* parameter values x */
 165   gsl_vector * f;             /* residual vector f(x) */
 166   gsl_vector * dx;            /* step dx */
 167   gsl_vector * g;             /* gradient J^T f */
 168   gsl_matrix * J;             /* Jacobian J(x) */
 169   gsl_vector * sqrt_wts_work; /* sqrt(W) */
 170   gsl_vector * sqrt_wts;      /* ptr to sqrt_wts_work, or NULL if not using weights */
 171   size_t niter;               /* number of iterations performed */
 172   gsl_multifit_nlinear_parameters params;
 173   void *state;
 174 } gsl_multifit_nlinear_workspace;
 175
 176 gsl_multifit_nlinear_workspace *
 177 gsl_multifit_nlinear_alloc (const gsl_multifit_nlinear_type * T,
 178                             const gsl_multifit_nlinear_parameters * params,
 179                             size_t n, size_t p);
 180
 181 void gsl_multifit_nlinear_free (gsl_multifit_nlinear_workspace * w);
 182
 183 gsl_multifit_nlinear_parameters gsl_multifit_nlinear_default_parameters(void);
 184
 185 int
 186 gsl_multifit_nlinear_init (const gsl_vector * x,
 187                            gsl_multifit_nlinear_fdf * fdf,
 188                            gsl_multifit_nlinear_workspace * w);
 189
 190 int gsl_multifit_nlinear_winit (const gsl_vector * x,
 191                                 const gsl_vector * wts,
 192                                 gsl_multifit_nlinear_fdf * fdf,
 193                                 gsl_multifit_nlinear_workspace * w);
 194
 195 int
 196 gsl_multifit_nlinear_iterate (gsl_multifit_nlinear_workspace * w);
 197
 198 double
 199 gsl_multifit_nlinear_avratio (const gsl_multifit_nlinear_workspace * w);
 200
 201 int
 202 gsl_multifit_nlinear_driver (const size_t maxiter,
 203                              const double xtol,
 204                              const double gtol,
 205                              const double ftol,
 206                              void (*callback)(const size_t iter, void *params,
 207                                               const gsl_multifit_nlinear_workspace *w),
 208                              void *callback_params,
 209                              int *info,
 210                              gsl_multifit_nlinear_workspace * w);
 211
 212 gsl_matrix *
 213 gsl_multifit_nlinear_jac (const gsl_multifit_nlinear_workspace * w);
 214
 215 const char *
 216 gsl_multifit_nlinear_name (const gsl_multifit_nlinear_workspace * w);
 217
 218 gsl_vector *
 219 gsl_multifit_nlinear_position (const gsl_multifit_nlinear_workspace * w);
 220
 221 gsl_vector *
 222 gsl_multifit_nlinear_residual (const gsl_multifit_nlinear_workspace * w);
 223
 224 size_t
 225 gsl_multifit_nlinear_niter (const gsl_multifit_nlinear_workspace * w);
 226
 227 int
 228 gsl_multifit_nlinear_rcond (double *rcond, const gsl_multifit_nlinear_workspace * w);
 229
 230 const char *
 231 gsl_multifit_nlinear_trs_name (const gsl_multifit_nlinear_workspace * w);
 232
 233 int gsl_multifit_nlinear_eval_f(gsl_multifit_nlinear_fdf *fdf,
 234                                 const gsl_vector *x,
 235                                 const gsl_vector *swts,
 236                                 gsl_vector *y);
 237
 238 int gsl_multifit_nlinear_eval_df(const gsl_vector *x,
 239                                  const gsl_vector *f,
 240                                  const gsl_vector *swts,
 241                                  const double h,
 242                                  const gsl_multifit_nlinear_fdtype fdtype,
 243                                  gsl_multifit_nlinear_fdf *fdf,
 244                                  gsl_matrix *df, gsl_vector *work);
 245
 246 int
 247 gsl_multifit_nlinear_eval_fvv(const double h,
 248                               const gsl_vector *x,
 249                               const gsl_vector *v,
 250                               const gsl_vector *f,
 251                               const gsl_matrix *J,
 252                               const gsl_vector *swts,
 253                               gsl_multifit_nlinear_fdf *fdf,
 254                               gsl_vector *yvv, gsl_vector *work);
 255
 256 /* covar.c */
 257 int
 258 gsl_multifit_nlinear_covar (const gsl_matrix * J, const double epsrel,
 259                             gsl_matrix * covar);
 260
 261 /* convergence.c */
 262 int
 263 gsl_multifit_nlinear_test (const double xtol, const double gtol,
 264                            const double ftol, int *info,
 265                            const gsl_multifit_nlinear_workspace * w);
 266
 267 /* fdjac.c */
 268 int
 269 gsl_multifit_nlinear_df(const double h, const gsl_multifit_nlinear_fdtype fdtype,
 270                         const gsl_vector *x, const gsl_vector *wts,
 271                         gsl_multifit_nlinear_fdf *fdf,
 272                         const gsl_vector *f, gsl_matrix *J, gsl_vector *work);
 273
 274 /* fdfvv.c */
 275 int
 276 gsl_multifit_nlinear_fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
 277                            const gsl_vector *f, const gsl_matrix *J,
 278                            const gsl_vector *swts, gsl_multifit_nlinear_fdf *fdf,
 279                            gsl_vector *fvv, gsl_vector *work);
 280
 281 /* top-level algorithms */
 282 GSL_VAR const gsl_multifit_nlinear_type * gsl_multifit_nlinear_trust;
 283
 284 /* trust region subproblem methods */
 285 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lm;
 286 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lmaccel;
 287 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_dogleg;
 288 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_ddogleg;
 289 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_subspace2D;
 290
 291 /* scaling matrix strategies */
 292 GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_levenberg;
 293 GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_marquardt;
 294 GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_more;
 295
 296 /* linear solvers */
 297 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_cholesky;
 298 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_mcholesky;
 299 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_qr;
 300 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_svd;
 301
 302 __END_DECLS
 303
 304 #endif /* __GSL_MULTIFIT_NLINEAR_H__ */