include/gsl/gsl_multilarge_nlinear.h

0001 /* multilarge_nlinear/gsl_multilarge_nlinear.h
0002  *
0003  * Copyright (C) 2015, 2016 Patrick Alken
0004  *
0005  * This program is free software; you can redistribute it and/or modify
0006  * it under the terms of the GNU General Public License as published by
0007  * the Free Software Foundation; either version 3 of the License, or (at
0008  * your option) any later version.
0009  *
0010  * This program is distributed in the hope that it will be useful, but
0011  * WITHOUT ANY WARRANTY; without even the implied warranty of
0012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0013  * General Public License for more details.
0014  *
0015  * You should have received a copy of the GNU General Public License
0016  * along with this program; if not, write to the Free Software
0017  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0018  */
0019
0020 #ifndef __GSL_MULTILARGE_NLINEAR_H__
0021 #define __GSL_MULTILARGE_NLINEAR_H__
0022
0023 #include <stdlib.h>
0024 #include <gsl/gsl_types.h>
0025 #include <gsl/gsl_math.h>
0026 #include <gsl/gsl_vector.h>
0027 #include <gsl/gsl_matrix.h>
0028 #include <gsl/gsl_permutation.h>
0029 #include <gsl/gsl_blas.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_MULTILARGE_NLINEAR_FWDIFF,
0046   GSL_MULTILARGE_NLINEAR_CTRDIFF
0047 } gsl_multilarge_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) (CBLAS_TRANSPOSE_t TransJ, const gsl_vector * x,
0056               const gsl_vector * u, void * params, gsl_vector * v,
0057               gsl_matrix * JTJ);
0058   int (* fvv) (const gsl_vector * x, const gsl_vector * v, void * params,
0059                gsl_vector * fvv);
0060   size_t n;        /* number of functions */
0061   size_t p;        /* number of independent variables */
0062   void * params;   /* user parameters */
0063   size_t nevalf;   /* number of function evaluations */
0064   size_t nevaldfu; /* number of Jacobian matrix-vector evaluations */
0065   size_t nevaldf2; /* number of Jacobian J^T J evaluations */
0066   size_t nevalfvv; /* number of fvv evaluations */
0067 } gsl_multilarge_nlinear_fdf;
0068
0069 /* trust region subproblem method */
0070 typedef struct
0071 {
0072   const char *name;
0073   void * (*alloc) (const void * params, const size_t n, const size_t p);
0074   int (*init) (const void * vtrust_state, void * vstate);
0075   int (*preloop) (const void * vtrust_state, void * vstate);
0076   int (*step) (const void * vtrust_state, const double delta,
0077                gsl_vector * dx, void * vstate);
0078   int (*preduction) (const void * vtrust_state, const gsl_vector * dx,
0079                      double * pred, void * vstate);
0080   void (*free) (void * vstate);
0081 } gsl_multilarge_nlinear_trs;
0082
0083 /* scaling matrix specification */
0084 typedef struct
0085 {
0086   const char *name;
0087   int (*init) (const gsl_matrix * JTJ, gsl_vector * diag);
0088   int (*update) (const gsl_matrix * JTJ, gsl_vector * diag);
0089 } gsl_multilarge_nlinear_scale;
0090
0091 /*
0092  * linear least squares solvers - there are three steps to
0093  * solving a least squares problem using a direct method:
0094  *
0095  * 1. init: called once per iteration when a new Jacobian matrix
0096  *          is required; form normal equations matrix J^T J
0097  * 2. presolve: called each time a new LM parameter value mu is available;
0098  *              used for cholesky method in order to factor
0099  *              the (J^T J + mu D^T D) matrix
0100  * 3. solve: solve the least square system for a given rhs
0101  */
0102 typedef struct
0103 {
0104   const char *name;
0105   void * (*alloc) (const size_t n, const size_t p);
0106   int (*init) (const void * vtrust_state, void * vstate);
0107   int (*presolve) (const double mu, const void * vtrust_state, void * vstate);
0108   int (*solve) (const gsl_vector * g, gsl_vector * x,
0109                 const void * vtrust_state, void * vstate);
0110   int (*rcond) (double * rcond, const gsl_matrix * JTJ, void * vstate);
0111   int (*covar) (const gsl_matrix * JTJ, gsl_matrix * covar, void * vstate);
0112   void (*free) (void * vstate);
0113 } gsl_multilarge_nlinear_solver;
0114
0115 /* tunable parameters */
0116 typedef struct
0117 {
0118   const gsl_multilarge_nlinear_trs *trs;       /* trust region subproblem method */
0119   const gsl_multilarge_nlinear_scale *scale;   /* scaling method */
0120   const gsl_multilarge_nlinear_solver *solver; /* solver method */
0121   gsl_multilarge_nlinear_fdtype fdtype;        /* finite difference method */
0122   double factor_up;                            /* factor for increasing trust radius */
0123   double factor_down;                          /* factor for decreasing trust radius */
0124   double avmax;                                /* max allowed |a|/|v| */
0125   double h_df;                                 /* step size for finite difference Jacobian */
0126   double h_fvv;                                /* step size for finite difference fvv */
0127   size_t max_iter;                             /* maximum iterations for trs method */
0128   double tol;                                  /* tolerance for solving trs */
0129 } gsl_multilarge_nlinear_parameters;
0130
0131 typedef struct
0132 {
0133   const char *name;
0134   void * (*alloc) (const gsl_multilarge_nlinear_parameters * params,
0135                    const size_t n, const size_t p);
0136   int (*init) (void * state, const gsl_vector * wts,
0137                gsl_multilarge_nlinear_fdf * fdf, const gsl_vector * x,
0138                gsl_vector * f, gsl_vector * g, gsl_matrix * JTJ);
0139   int (*iterate) (void * state, const gsl_vector * wts,
0140                   gsl_multilarge_nlinear_fdf * fdf, gsl_vector * x,
0141                   gsl_vector * f, gsl_vector * g, gsl_matrix * JTJ,
0142                   gsl_vector * dx);
0143   int (*rcond) (double * rcond, const gsl_matrix * JTJ, void * state);
0144   int (*covar) (const gsl_matrix * JTJ, gsl_matrix * covar, void * state);
0145   double (*avratio) (void * state);
0146   void (*free) (void * state);
0147 } gsl_multilarge_nlinear_type;
0148
0149 /* current state passed to low-level trust region algorithms */
0150 typedef struct
0151 {
0152   const gsl_vector * x;             /* parameter values x */
0153   const gsl_vector * f;             /* residual vector f(x) */
0154   const gsl_vector * g;             /* gradient J^T f */
0155   const gsl_matrix * JTJ;           /* matrix J^T J */
0156   const gsl_vector * diag;          /* scaling matrix D */
0157   const gsl_vector * sqrt_wts;      /* sqrt(diag(W)) or NULL for unweighted */
0158   const double *mu;                 /* LM parameter */
0159   const gsl_multilarge_nlinear_parameters * params;
0160   void *solver_state;               /* workspace for direct least squares solver */
0161   gsl_multilarge_nlinear_fdf * fdf;
0162   double *avratio;                  /* |a| / |v| */
0163 } gsl_multilarge_nlinear_trust_state;
0164
0165 typedef struct
0166 {
0167   const gsl_multilarge_nlinear_type * type;
0168   gsl_multilarge_nlinear_fdf * fdf ;
0169   gsl_vector * x;             /* parameter values x */
0170   gsl_vector * f;             /* residual vector f(x) */
0171   gsl_vector * dx;            /* step dx */
0172   gsl_vector * g;             /* gradient J^T f */
0173   gsl_matrix * JTJ;           /* matrix J^T J */
0174   gsl_vector * sqrt_wts_work; /* sqrt(W) */
0175   gsl_vector * sqrt_wts;      /* ptr to sqrt_wts_work, or NULL if not using weights */
0176   size_t n;                   /* number of residuals */
0177   size_t p;                   /* number of parameters */
0178   size_t niter;               /* number of iterations performed */
0179   gsl_multilarge_nlinear_parameters params;
0180   void *state;
0181 } gsl_multilarge_nlinear_workspace;
0182
0183 gsl_multilarge_nlinear_workspace *
0184 gsl_multilarge_nlinear_alloc (const gsl_multilarge_nlinear_type * T,
0185                               const gsl_multilarge_nlinear_parameters * params,
0186                               size_t n, size_t p);
0187
0188 void gsl_multilarge_nlinear_free (gsl_multilarge_nlinear_workspace * w);
0189
0190 gsl_multilarge_nlinear_parameters gsl_multilarge_nlinear_default_parameters(void);
0191
0192 int
0193 gsl_multilarge_nlinear_init (const gsl_vector * x,
0194                              gsl_multilarge_nlinear_fdf * fdf,
0195                              gsl_multilarge_nlinear_workspace * w);
0196
0197 int gsl_multilarge_nlinear_winit (const gsl_vector * x,
0198                                   const gsl_vector * wts,
0199                                   gsl_multilarge_nlinear_fdf * fdf,
0200                                   gsl_multilarge_nlinear_workspace * w);
0201
0202 int
0203 gsl_multilarge_nlinear_iterate (gsl_multilarge_nlinear_workspace * w);
0204
0205 double
0206 gsl_multilarge_nlinear_avratio (const gsl_multilarge_nlinear_workspace * w);
0207
0208 int
0209 gsl_multilarge_nlinear_rcond (double * rcond, const gsl_multilarge_nlinear_workspace * w);
0210
0211 int
0212 gsl_multilarge_nlinear_covar (gsl_matrix * covar, gsl_multilarge_nlinear_workspace * w);
0213
0214 int
0215 gsl_multilarge_nlinear_driver (const size_t maxiter,
0216                                const double xtol,
0217                                const double gtol,
0218                                const double ftol,
0219                                void (*callback)(const size_t iter, void *params,
0220                                                 const gsl_multilarge_nlinear_workspace *w),
0221                                void *callback_params,
0222                                int *info,
0223                                gsl_multilarge_nlinear_workspace * w);
0224
0225 const char *
0226 gsl_multilarge_nlinear_name (const gsl_multilarge_nlinear_workspace * w);
0227
0228 gsl_vector *
0229 gsl_multilarge_nlinear_position (const gsl_multilarge_nlinear_workspace * w);
0230
0231 gsl_vector *
0232 gsl_multilarge_nlinear_residual (const gsl_multilarge_nlinear_workspace * w);
0233
0234 gsl_vector *
0235 gsl_multilarge_nlinear_step (const gsl_multilarge_nlinear_workspace * w);
0236
0237 size_t
0238 gsl_multilarge_nlinear_niter (const gsl_multilarge_nlinear_workspace * w);
0239
0240 const char *
0241 gsl_multilarge_nlinear_trs_name (const gsl_multilarge_nlinear_workspace * w);
0242
0243 int gsl_multilarge_nlinear_eval_f(gsl_multilarge_nlinear_fdf *fdf,
0244                                   const gsl_vector *x,
0245                                   const gsl_vector *swts,
0246                                   gsl_vector *y);
0247
0248 int
0249 gsl_multilarge_nlinear_eval_df(const CBLAS_TRANSPOSE_t TransJ,
0250                                const gsl_vector *x,
0251                                const gsl_vector *f,
0252                                const gsl_vector *u,
0253                                const gsl_vector *swts,
0254                                const double h,
0255                                const gsl_multilarge_nlinear_fdtype fdtype,
0256                                gsl_multilarge_nlinear_fdf *fdf,
0257                                gsl_vector *v,
0258                                gsl_matrix *JTJ,
0259                                gsl_vector *work);
0260
0261 int
0262 gsl_multilarge_nlinear_eval_fvv(const double h,
0263                                 const gsl_vector *x,
0264                                 const gsl_vector *v,
0265                                 const gsl_vector *f,
0266                                 const gsl_vector *swts,
0267                                 gsl_multilarge_nlinear_fdf *fdf,
0268                                 gsl_vector *yvv,
0269                                 gsl_vector *work);
0270
0271 /* convergence.c */
0272 int
0273 gsl_multilarge_nlinear_test (const double xtol, const double gtol,
0274                              const double ftol, int *info,
0275                              const gsl_multilarge_nlinear_workspace * w);
0276
0277 /* fdjac.c */
0278 int
0279 gsl_multilarge_nlinear_df(const double h, const gsl_multilarge_nlinear_fdtype fdtype,
0280                           const gsl_vector *x, const gsl_vector *wts,
0281                           gsl_multilarge_nlinear_fdf *fdf,
0282                           const gsl_vector *f, gsl_matrix *J, gsl_vector *work);
0283
0284 /* fdfvv.c */
0285 int
0286 gsl_multilarge_nlinear_fdfvv(const double h, const gsl_vector *x, const gsl_vector *v,
0287                              const gsl_vector *f, const gsl_matrix *J,
0288                              const gsl_vector *swts, gsl_multilarge_nlinear_fdf *fdf,
0289                              gsl_vector *fvv, gsl_vector *work);
0290
0291 /* top-level algorithms */
0292 GSL_VAR const gsl_multilarge_nlinear_type * gsl_multilarge_nlinear_trust;
0293
0294 /* trust region subproblem methods */
0295 GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_lm;
0296 GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_lmaccel;
0297 GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_dogleg;
0298 GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_ddogleg;
0299 GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_subspace2D;
0300 GSL_VAR const gsl_multilarge_nlinear_trs * gsl_multilarge_nlinear_trs_cgst;
0301
0302 /* scaling matrix strategies */
0303 GSL_VAR const gsl_multilarge_nlinear_scale * gsl_multilarge_nlinear_scale_levenberg;
0304 GSL_VAR const gsl_multilarge_nlinear_scale * gsl_multilarge_nlinear_scale_marquardt;
0305 GSL_VAR const gsl_multilarge_nlinear_scale * gsl_multilarge_nlinear_scale_more;
0306
0307 /* linear solvers */
0308 GSL_VAR const gsl_multilarge_nlinear_solver * gsl_multilarge_nlinear_solver_cholesky;
0309 GSL_VAR const gsl_multilarge_nlinear_solver * gsl_multilarge_nlinear_solver_mcholesky;
0310 GSL_VAR const gsl_multilarge_nlinear_solver * gsl_multilarge_nlinear_solver_none;
0311
0312 __END_DECLS
0313
0314 #endif /* __GSL_MULTILARGE_NLINEAR_H__ */