include/gsl/gsl_multifit_nlinear.h

0001 /* multifit_nlinear/gsl_multifit_nlinear.h
0002  *
0003  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
0004  * Copyright (C) 2015, 2016 Patrick Alken
0005  *
0006  * This program is free software; you can redistribute it and/or modify
0007  * it under the terms of the GNU General Public License as published by
0008  * the Free Software Foundation; either version 3 of the License, or (at
0009  * your option) any later version.
0010  *
0011  * This program is distributed in the hope that it will be useful, but
0012  * WITHOUT ANY WARRANTY; without even the implied warranty of
0013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0014  * General Public License for more details.
0015  *
0016  * You should have received a copy of the GNU General Public License
0017  * along with this program; if not, write to the Free Software
0018  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0019  */
0020
0021 #ifndef __GSL_MULTIFIT_NLINEAR_H__
0022 #define __GSL_MULTIFIT_NLINEAR_H__
0023
0024 #include <stdlib.h>
0025 #include <gsl/gsl_types.h>
0026 #include <gsl/gsl_math.h>
0027 #include <gsl/gsl_vector.h>
0028 #include <gsl/gsl_matrix.h>
0029 #include <gsl/gsl_permutation.h>
0030
0031 #undef __BEGIN_DECLS
0032 #undef __END_DECLS
0033 #ifdef __cplusplus
0034 # define __BEGIN_DECLS extern "C" {
0035 # define __END_DECLS }
0036 #else
0037 # define __BEGIN_DECLS /* empty */
0038 # define __END_DECLS /* empty */
0039 #endif
0040
0041 __BEGIN_DECLS
0042
0043 typedef enum
0044 {
0045   GSL_MULTIFIT_NLINEAR_FWDIFF,
0046   GSL_MULTIFIT_NLINEAR_CTRDIFF
0047 } gsl_multifit_nlinear_fdtype;
0048
0049 /* Definition of vector-valued functions and gradient with parameters
0050    based on gsl_vector */
0051
0052 typedef struct
0053 {
0054   int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
0055   int (* df) (const gsl_vector * x, void * params, gsl_matrix * df);
0056   int (* fvv) (const gsl_vector * x, const gsl_vector * v, void * params,
0057                gsl_vector * fvv);
0058   size_t n;        /* number of functions */
0059   size_t p;        /* number of independent variables */
0060   void * params;   /* user parameters */
0061   size_t nevalf;   /* number of function evaluations */
0062   size_t nevaldf;  /* number of Jacobian evaluations */
0063   size_t nevalfvv; /* number of fvv evaluations */
0064 } gsl_multifit_nlinear_fdf;
0065
0066 /* trust region subproblem method */
0067 typedef struct
0068 {
0069   const char *name;
0070   void * (*alloc) (const void * params, const size_t n, const size_t p);
0071   int (*init) (const void * vtrust_state, void * vstate);
0072   int (*preloop) (const void * vtrust_state, void * vstate);
0073   int (*step) (const void * vtrust_state, const double delta,
0074                gsl_vector * dx, void * vstate);
0075   int (*preduction) (const void * vtrust_state, const gsl_vector * dx,
0076                      double * pred, void * vstate);
0077   void (*free) (void * vstate);
0078 } gsl_multifit_nlinear_trs;
0079
0080 /* scaling matrix specification */
0081 typedef struct
0082 {
0083   const char *name;
0084   int (*init) (const gsl_matrix * J, gsl_vector * diag);
0085   int (*update) (const gsl_matrix * J, gsl_vector * diag);
0086 } gsl_multifit_nlinear_scale;
0087
0088 /*
0089  * linear least squares solvers - there are three steps to
0090  * solving a least squares problem using a trust region
0091  * method:
0092  *
0093  * 1. init: called once per iteration when a new Jacobian matrix
0094  *          is computed; perform factorization of Jacobian (qr,svd)
0095  *          or form normal equations matrix (cholesky)
0096  * 2. presolve: called each time a new LM parameter value mu is available;
0097  *              used for cholesky method in order to factor
0098  *              the (J^T J + mu D^T D) matrix
0099  * 3. solve: solve the least square system for a given rhs
0100  */
0101 typedef struct
0102 {
0103   const char *name;
0104   void * (*alloc) (const size_t n, const size_t p);
0105   int (*init) (const void * vtrust_state, void * vstate);
0106   int (*presolve) (const double mu, const void * vtrust_state, void * vstate);
0107   int (*solve) (const gsl_vector * f, gsl_vector * x,
0108                 const void * vtrust_state, void * vstate);
0109   int (*rcond) (double * rcond, void * vstate);
0110   void (*free) (void * vstate);
0111 } gsl_multifit_nlinear_solver;
0112
0113 /* tunable parameters */
0114 typedef struct
0115 {
0116   const gsl_multifit_nlinear_trs *trs;        /* trust region subproblem method */
0117   const gsl_multifit_nlinear_scale *scale;    /* scaling method */
0118   const gsl_multifit_nlinear_solver *solver;  /* solver method */
0119   gsl_multifit_nlinear_fdtype fdtype;         /* finite difference method */
0120   double factor_up;                           /* factor for increasing trust radius */
0121   double factor_down;                         /* factor for decreasing trust radius */
0122   double avmax;                               /* max allowed |a|/|v| */
0123   double h_df;                                /* step size for finite difference Jacobian */
0124   double h_fvv;                               /* step size for finite difference fvv */
0125 } gsl_multifit_nlinear_parameters;
0126
0127 typedef struct
0128 {
0129   const char *name;
0130   void * (*alloc) (const gsl_multifit_nlinear_parameters * params,
0131                    const size_t n, const size_t p);
0132   int (*init) (void * state, const gsl_vector * wts,
0133                gsl_multifit_nlinear_fdf * fdf, const gsl_vector * x,
0134                gsl_vector * f, gsl_matrix * J, gsl_vector * g);
0135   int (*iterate) (void * state, const gsl_vector * wts,
0136                   gsl_multifit_nlinear_fdf * fdf, gsl_vector * x,
0137                   gsl_vector * f, gsl_matrix * J, gsl_vector * g,
0138                   gsl_vector * dx);
0139   int (*rcond) (double * rcond, void * state);
0140   double (*avratio) (void * state);
0141   void (*free) (void * state);
0142 } gsl_multifit_nlinear_type;
0143
0144 /* current state passed to low-level trust region algorithms */
0145 typedef struct
0146 {
0147   const gsl_vector * x;             /* parameter values x */
0148   const gsl_vector * f;             /* residual vector f(x) */
0149   const gsl_vector * g;             /* gradient J^T f */
0150   const gsl_matrix * J;             /* Jacobian J(x) */
0151   const gsl_vector * diag;          /* scaling matrix D */
0152   const gsl_vector * sqrt_wts;      /* sqrt(diag(W)) or NULL for unweighted */
0153   const double *mu;                 /* LM parameter */
0154   const gsl_multifit_nlinear_parameters * params;
0155   void *solver_state;               /* workspace for linear least squares solver */
0156   gsl_multifit_nlinear_fdf * fdf;
0157   double *avratio;                  /* |a| / |v| */
0158 } gsl_multifit_nlinear_trust_state;
0159
0160 typedef struct
0161 {
0162   const gsl_multifit_nlinear_type * type;
0163   gsl_multifit_nlinear_fdf * fdf ;
0164   gsl_vector * x;             /* parameter values x */
0165   gsl_vector * f;             /* residual vector f(x) */
0166   gsl_vector * dx;            /* step dx */
0167   gsl_vector * g;             /* gradient J^T f */
0168   gsl_matrix * J;             /* Jacobian J(x) */
0169   gsl_vector * sqrt_wts_work; /* sqrt(W) */
0170   gsl_vector * sqrt_wts;      /* ptr to sqrt_wts_work, or NULL if not using weights */
0171   size_t niter;               /* number of iterations performed */
0172   gsl_multifit_nlinear_parameters params;
0173   void *state;
0174 } gsl_multifit_nlinear_workspace;
0175
0176 gsl_multifit_nlinear_workspace *
0177 gsl_multifit_nlinear_alloc (const gsl_multifit_nlinear_type * T,
0178                             const gsl_multifit_nlinear_parameters * params,
0179                             size_t n, size_t p);
0180
0181 void gsl_multifit_nlinear_free (gsl_multifit_nlinear_workspace * w);
0182
0183 gsl_multifit_nlinear_parameters gsl_multifit_nlinear_default_parameters(void);
0184
0185 int
0186 gsl_multifit_nlinear_init (const gsl_vector * x,
0187                            gsl_multifit_nlinear_fdf * fdf,
0188                            gsl_multifit_nlinear_workspace * w);
0189
0190 int gsl_multifit_nlinear_winit (const gsl_vector * x,
0191                                 const gsl_vector * wts,
0192                                 gsl_multifit_nlinear_fdf * fdf,
0193                                 gsl_multifit_nlinear_workspace * w);
0194
0195 int
0196 gsl_multifit_nlinear_iterate (gsl_multifit_nlinear_workspace * w);
0197
0198 double
0199 gsl_multifit_nlinear_avratio (const gsl_multifit_nlinear_workspace * w);
0200
0201 int
0202 gsl_multifit_nlinear_driver (const size_t maxiter,
0203                              const double xtol,
0204                              const double gtol,
0205                              const double ftol,
0206                              void (*callback)(const size_t iter, void *params,
0207                                               const gsl_multifit_nlinear_workspace *w),
0208                              void *callback_params,
0209                              int *info,
0210                              gsl_multifit_nlinear_workspace * w);
0211
0212 gsl_matrix *
0213 gsl_multifit_nlinear_jac (const gsl_multifit_nlinear_workspace * w);
0214
0215 const char *
0216 gsl_multifit_nlinear_name (const gsl_multifit_nlinear_workspace * w);
0217
0218 gsl_vector *
0219 gsl_multifit_nlinear_position (const gsl_multifit_nlinear_workspace * w);
0220
0221 gsl_vector *
0222 gsl_multifit_nlinear_residual (const gsl_multifit_nlinear_workspace * w);
0223
0224 size_t
0225 gsl_multifit_nlinear_niter (const gsl_multifit_nlinear_workspace * w);
0226
0227 int
0228 gsl_multifit_nlinear_rcond (double *rcond, const gsl_multifit_nlinear_workspace * w);
0229
0230 const char *
0231 gsl_multifit_nlinear_trs_name (const gsl_multifit_nlinear_workspace * w);
0232
0233 int gsl_multifit_nlinear_eval_f(gsl_multifit_nlinear_fdf *fdf,
0234                                 const gsl_vector *x,
0235                                 const gsl_vector *swts,
0236                                 gsl_vector *y);
0237
0238 int gsl_multifit_nlinear_eval_df(const gsl_vector *x,
0239                                  const gsl_vector *f,
0240                                  const gsl_vector *swts,
0241                                  const double h,
0242                                  const gsl_multifit_nlinear_fdtype fdtype,
0243                                  gsl_multifit_nlinear_fdf *fdf,
0244                                  gsl_matrix *df, gsl_vector *work);
0245
0246 int
0247 gsl_multifit_nlinear_eval_fvv(const double h,
0248                               const gsl_vector *x,
0249                               const gsl_vector *v,
0250                               const gsl_vector *f,
0251                               const gsl_matrix *J,
0252                               const gsl_vector *swts,
0253                               gsl_multifit_nlinear_fdf *fdf,
0254                               gsl_vector *yvv, gsl_vector *work);
0255
0256 /* covar.c */
0257 int
0258 gsl_multifit_nlinear_covar (const gsl_matrix * J, const double epsrel,
0259                             gsl_matrix * covar);
0260
0261 /* convergence.c */
0262 int
0263 gsl_multifit_nlinear_test (const double xtol, const double gtol,
0264                            const double ftol, int *info,
0265                            const gsl_multifit_nlinear_workspace * w);
0266
0267 /* fdjac.c */
0268 int
0269 gsl_multifit_nlinear_df(const double h, const gsl_multifit_nlinear_fdtype fdtype,
0270                         const gsl_vector *x, const gsl_vector *wts,
0271                         gsl_multifit_nlinear_fdf *fdf,
0272                         const gsl_vector *f, gsl_matrix *J, gsl_vector *work);
0273
0274 /* fdfvv.c */
0275 int
0276 gsl_multifit_nlinear_fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
0277                            const gsl_vector *f, const gsl_matrix *J,
0278                            const gsl_vector *swts, gsl_multifit_nlinear_fdf *fdf,
0279                            gsl_vector *fvv, gsl_vector *work);
0280
0281 /* top-level algorithms */
0282 GSL_VAR const gsl_multifit_nlinear_type * gsl_multifit_nlinear_trust;
0283
0284 /* trust region subproblem methods */
0285 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lm;
0286 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_lmaccel;
0287 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_dogleg;
0288 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_ddogleg;
0289 GSL_VAR const gsl_multifit_nlinear_trs * gsl_multifit_nlinear_trs_subspace2D;
0290
0291 /* scaling matrix strategies */
0292 GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_levenberg;
0293 GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_marquardt;
0294 GSL_VAR const gsl_multifit_nlinear_scale * gsl_multifit_nlinear_scale_more;
0295
0296 /* linear solvers */
0297 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_cholesky;
0298 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_mcholesky;
0299 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_qr;
0300 GSL_VAR const gsl_multifit_nlinear_solver * gsl_multifit_nlinear_solver_svd;
0301
0302 __END_DECLS
0303
0304 #endif /* __GSL_MULTIFIT_NLINEAR_H__ */