include/gsl/gsl_bspline.h

0001 /* bspline/gsl_bspline.h
0002  *
0003  * Copyright (C) 2018, 2019, 2020, 2021 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_BSPLINE_H__
0021 #define __GSL_BSPLINE_H__
0022
0023 #include <stdlib.h>
0024 #include <gsl/gsl_math.h>
0025 #include <gsl/gsl_inline.h>
0026 #include <gsl/gsl_vector.h>
0027 #include <gsl/gsl_matrix.h>
0028
0029 #undef __BEGIN_DECLS
0030 #undef __END_DECLS
0031 #ifdef __cplusplus
0032 # define __BEGIN_DECLS extern "C" {
0033 # define __END_DECLS }
0034 #else
0035 # define __BEGIN_DECLS /* empty */
0036 # define __END_DECLS /* empty */
0037 #endif
0038
0039 __BEGIN_DECLS
0040
0041 typedef struct
0042 {
0043     size_t spline_order;  /* spline order */
0044     size_t nbreak;        /* number of breakpoints */
0045     size_t ncontrol;      /* number of bspline control points (nbreak + spline_order - 2) */
0046
0047     gsl_vector *knots;    /* knots vector, length ncontrol + spline_order */
0048     gsl_vector *deltal;   /* left delta, length spline_order */
0049     gsl_vector *deltar;   /* right delta, length spline_order */
0050     gsl_vector *B;        /* temporary spline results, length spline_order */
0051     gsl_matrix *XTX;      /* stores diagonals of banded normal equations matrix, ncontrol-by-spline_order */
0052     gsl_matrix *R;        /* R factor for periodic least squares fitting, ncontrol-by-ncontrol */
0053     gsl_vector *work;     /* workspace, length 3*ncontrol */
0054
0055     /* bspline derivative parameters */
0056     gsl_matrix *A;        /* work matrix, spline_order-by-spline_order */
0057     gsl_matrix *dB;       /* temporary derivative results, spline_order-by-[2*(spline_order+1)] */
0058
0059     size_t icache;        /* cached index of current interval, in [0,n+k-2] */
0060 } gsl_bspline_workspace;
0061
0062 /* bspline.c */
0063
0064 gsl_bspline_workspace * gsl_bspline_alloc(const size_t k, const size_t nbreak);
0065
0066 gsl_bspline_workspace * gsl_bspline_alloc_ncontrol (const size_t spline_order, const size_t ncontrol);
0067
0068 void gsl_bspline_free(gsl_bspline_workspace *w);
0069
0070 size_t gsl_bspline_ncontrol(const gsl_bspline_workspace * w);
0071 size_t gsl_bspline_order(const gsl_bspline_workspace * w);
0072 size_t gsl_bspline_nbreak(const gsl_bspline_workspace * w);
0073 double gsl_bspline_breakpoint(const size_t i, const gsl_bspline_workspace * w);
0074 double gsl_bspline_greville_abscissa(const size_t i, const gsl_bspline_workspace *w);
0075
0076 int gsl_bspline_init_augment(const gsl_vector *breakpts, gsl_bspline_workspace *w);
0077
0078 int gsl_bspline_init_uniform(const double a, const double b, gsl_bspline_workspace *w);
0079
0080 int gsl_bspline_init_periodic (const double a, const double b, gsl_bspline_workspace * w);
0081
0082 int gsl_bspline_init (const gsl_vector * t, gsl_bspline_workspace * w);
0083
0084 int gsl_bspline_proj_rhs(const gsl_function * F, gsl_vector * y, gsl_bspline_workspace * w);
0085
0086 INLINE_DECL size_t gsl_bspline_find_interval (const double x, int *flag, gsl_bspline_workspace * w);
0087
0088 /* eval.c */
0089
0090 int gsl_bspline_calc(const double x, const gsl_vector * c,
0091                      double * result, gsl_bspline_workspace * w);
0092
0093 int gsl_bspline_calc_deriv(const double x, const gsl_vector * c, const size_t nderiv,
0094                            double * result, gsl_bspline_workspace * w);
0095
0096 int gsl_bspline_vector_calc(const double x, const gsl_matrix * c, gsl_vector * result,
0097                             gsl_bspline_workspace * w);
0098
0099 int gsl_bspline_vector_calc_deriv(const double x, const gsl_matrix * c, const size_t nderiv,
0100                                   gsl_vector * result, gsl_bspline_workspace * w);
0101
0102 int
0103 gsl_bspline_eval_basis(const double x, gsl_vector *B, gsl_bspline_workspace *w);
0104
0105 int
0106 gsl_bspline_basis(const double x, gsl_vector *Bk, size_t *istart,
0107                   gsl_bspline_workspace *w);
0108
0109 int
0110 gsl_bspline_eval_deriv_basis(const double x, const size_t nderiv,
0111                              gsl_matrix *dB, gsl_bspline_workspace *w);
0112
0113 int
0114 gsl_bspline_basis_deriv(const double x, const size_t nderiv,
0115                         gsl_matrix *dB, size_t *istart,
0116                         gsl_bspline_workspace *w);
0117
0118 int
0119 gsl_bspline_init_greville(const gsl_vector *abscissae,
0120                           gsl_bspline_workspace *w,
0121                           double *abserr);
0122
0123 /* gram.c */
0124
0125 int gsl_bspline_gram(const size_t nderiv, gsl_matrix * G, gsl_bspline_workspace * w);
0126
0127 int gsl_bspline_gram_interval(const double a, const double b, const size_t nderiv,
0128                               gsl_matrix * G, gsl_bspline_workspace * w);
0129
0130 int gsl_bspline_oprod(const size_t nderiv, const double x, gsl_matrix * A, gsl_bspline_workspace * w);
0131
0132 /* integ.c */
0133
0134 int gsl_bspline_calc_integ(const double a, const double b,
0135                            const gsl_vector * c, double * result,
0136                            gsl_bspline_workspace * w);
0137
0138 int gsl_bspline_basis_integ(const double a, const double b,
0139                             gsl_vector * bint, gsl_bspline_workspace * w);
0140
0141 /* interp.c */
0142
0143 int gsl_bspline_init_interp (const gsl_vector * x, gsl_bspline_workspace * w);
0144
0145 int gsl_bspline_init_hermite(const size_t nderiv, const gsl_vector * x, gsl_bspline_workspace * w);
0146
0147 int gsl_bspline_col_interp(const gsl_vector * tau, gsl_matrix * XB, gsl_bspline_workspace * w);
0148
0149 int gsl_bspline_interp_chermite(const gsl_vector * x, const gsl_vector * y,
0150                                 const gsl_vector * dy, gsl_vector * c,
0151                                 const gsl_bspline_workspace * w);
0152
0153 /* ls.c */
0154
0155 int gsl_bspline_lssolve(const gsl_vector * x, const gsl_vector * y, gsl_vector * c,
0156                         double * chisq, gsl_bspline_workspace * w);
0157
0158 int gsl_bspline_wlssolve(const gsl_vector * x, const gsl_vector * y, const gsl_vector * wts,
0159                          gsl_vector * c, double * chisq, gsl_bspline_workspace * w);
0160
0161 int gsl_bspline_lsnormal(const gsl_vector * x, const gsl_vector * y, const gsl_vector * wts,
0162                          gsl_vector * XTy, gsl_matrix * XTX, gsl_bspline_workspace * w);
0163
0164 int gsl_bspline_lsnormalm(const gsl_vector * x, const gsl_matrix * Y, const gsl_vector * wts,
0165                           gsl_matrix * XTY, gsl_matrix * XTX, gsl_bspline_workspace * w);
0166
0167 int gsl_bspline_plssolve(const gsl_vector * x, const gsl_vector * y,
0168                          gsl_vector * c, double * chisq, gsl_bspline_workspace * w);
0169
0170 int gsl_bspline_pwlssolve(const gsl_vector * x, const gsl_vector * y, const gsl_vector * wts,
0171                           gsl_vector * c, double * chisq, gsl_bspline_workspace * w);
0172
0173 int gsl_bspline_plsqr(const gsl_vector * x, const gsl_vector * y, const gsl_vector * wts,
0174                       gsl_matrix * R, gsl_vector * QTy, double * rnorm,
0175                       gsl_bspline_workspace * w);
0176
0177 int gsl_bspline_residuals(const gsl_vector * x, const gsl_vector * y, const gsl_vector * c,
0178                           gsl_vector * r, gsl_bspline_workspace * w);
0179
0180 int gsl_bspline_covariance(const gsl_matrix * XTX, gsl_matrix * cov, gsl_bspline_workspace * w);
0181
0182 int gsl_bspline_rcond(const gsl_matrix * XTX, double * rcond, gsl_bspline_workspace * w);
0183
0184 int gsl_bspline_err(const double x, const size_t nderiv,
0185                     const gsl_matrix * cov, double * err,
0186                     gsl_bspline_workspace * w);
0187
0188 /* future to be deprecated functions */
0189
0190 size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * w);
0191 int gsl_bspline_knots (const gsl_vector * breakpts, gsl_bspline_workspace * w);
0192 int gsl_bspline_knots_uniform (const double a, const double b, gsl_bspline_workspace * w);
0193 int gsl_bspline_eval (const double x, gsl_vector * B, gsl_bspline_workspace * w);
0194 int gsl_bspline_deriv_eval (const double x, const size_t nderiv,
0195                             gsl_matrix * dB, gsl_bspline_workspace * w);
0196 int gsl_bspline_knots_greville (const gsl_vector *abscissae,
0197                                 gsl_bspline_workspace *w,
0198                                 double *abserr);
0199
0200 /* end of future to be deprecated functions */
0201
0202 #ifdef HAVE_INLINE
0203
0204 /*
0205 gsl_bspline_find_interval()
0206   Find knot interval such that t_i <= x < t_{i + 1}
0207 where the t_i are knot values. The algorithm uses binary search
0208 since the knot array is sorted ascending.
0209
0210 Inputs: x    - x value
0211         flag - (output) error flag
0212         w    - bspline workspace
0213
0214 Return: i (index in w->knots corresponding to left limit of interval)
0215
0216 Notes: The error conditions are reported as follows:
0217
0218        Condition                        Return value        Flag
0219        ---------                        ------------        ----
0220        x < t_0                               0               -1
0221        t_i <= x < t_{i+1}                    i                0
0222        t_i < x = t_{i+1} = t_{n+k-1}         i                0
0223        t_i < t_{i+1} = t_{n+k-1} < x         i               +1
0224 */
0225
0226 INLINE_FUN
0227 size_t
0228 gsl_bspline_find_interval (const double x, int *flag, gsl_bspline_workspace * w)
0229 {
0230   const double *array = w->knots->data; /* stride is 1 */
0231   size_t ilo;
0232   size_t ihi = w->knots->size - 1;
0233
0234   /* check for quick return */
0235   if (array[w->icache] <= x && x < array[w->icache + 1])
0236     {
0237       *flag = 0;
0238       return w->icache;
0239     }
0240   else if (x < *array)
0241     {
0242       *flag = -1;
0243       return 0;
0244     }
0245   else if (x >= array[ihi])
0246     {
0247       *flag = (x > array[ihi]);
0248
0249       while (ihi > 0)
0250         {
0251           ihi--;
0252           if (array[ihi] < array[ihi + 1])
0253             break;
0254         }
0255
0256       return ihi;
0257     }
0258
0259   *flag = 0;
0260   ilo = 0;
0261
0262   while (ihi > ilo + 1)
0263     {
0264       size_t i = (ihi + ilo) >> 1;
0265       if(array[i] > x)
0266         ihi = i;
0267       else
0268         ilo = i;
0269     }
0270
0271   w->icache = ilo;
0272
0273   return ilo;
0274 }
0275
0276 #endif /* HAVE_INLINE */
0277
0278 __END_DECLS
0279
0280 #endif /* __GSL_BSPLINE_H__ */