Move bspline derivative data into separate workspace. From: Rhys Ulerich Required for binary compatibility with existing code. --- bspline/TODO | 3 bspline/bspline.c | 558 +++++++++++++++++++++++++++---------------------- bspline/gsl_bspline.h | 66 ++++-- bspline/test.c | 20 +- doc/bspline.texi | 21 ++ 5 files changed, 389 insertions(+), 279 deletions(-) diff --git a/bspline/TODO b/bspline/TODO index 41e9bf0..37a52c3 100644 --- a/bspline/TODO +++ b/bspline/TODO @@ -1,5 +1,8 @@ Add functions: +add bspline derivatve code to examples/bspline.c or produce a +derivative-related example. + gsl_bspline_smooth to fit smoothing splines to data more efficiently than the standard least squares inversion (see pppack l2appr and smooth.spline() from GNU R) diff --git a/bspline/bspline.c b/bspline/bspline.c index 45bd628..f286707 100644 --- a/bspline/bspline.c +++ b/bspline/bspline.c @@ -36,12 +36,12 @@ * */ +static inline size_t +bspline_find_interval(const double x, int *flag, gsl_bspline_workspace *w); -static inline size_t bspline_find_interval(const double x, int *flag, - gsl_bspline_workspace *w); - -static inline int bspline_process_interval_for_eval( - const double x, size_t *i, int flag, gsl_bspline_workspace *w); +static inline int +bspline_process_interval_for_eval(const double x, size_t *i, int flag, + gsl_bspline_workspace *w); static void bspline_pppack_bsplvb(const gsl_vector *t, @@ -68,7 +68,7 @@ bspline_pppack_bsplvd(const gsl_vector *t, /* gsl_bspline_alloc() Allocate space for a bspline workspace. The size of the -workspace is O(5k + nbreak + 2k^2) +workspace is O(5k + nbreak) Inputs: k - spline order (cubic = 4) nbreak - number of breakpoints @@ -91,8 +91,7 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak) { gsl_bspline_workspace *w; - w = (gsl_bspline_workspace *) - calloc(1, sizeof(gsl_bspline_workspace)); + w = (gsl_bspline_workspace *) calloc(1, sizeof(gsl_bspline_workspace)); if (w == 0) { GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM); @@ -108,7 +107,8 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak) if (w->knots == 0) { gsl_bspline_free(w); - GSL_ERROR_NULL("failed to allocate space for knots vector", GSL_ENOMEM); + GSL_ERROR_NULL("failed to allocate space for knots vector", + GSL_ENOMEM); } w->deltal = gsl_vector_alloc(k); @@ -116,58 +116,120 @@ gsl_bspline_alloc(const size_t k, const size_t nbreak) if (!w->deltal || !w->deltar) { gsl_bspline_free(w); - GSL_ERROR_NULL("failed to allocate space for delta vectors", GSL_ENOMEM); + GSL_ERROR_NULL("failed to allocate space for delta vectors", + GSL_ENOMEM); } w->A = gsl_matrix_alloc(k, k); if (w->A == 0) { gsl_bspline_free(w); - GSL_ERROR_NULL("failed to allocate space for derivative work matrix", GSL_ENOMEM); + GSL_ERROR_NULL("failed to allocate space for derivative work matrix", + GSL_ENOMEM); } w->B = gsl_vector_alloc(k); if (w->B == 0) { gsl_bspline_free(w); - GSL_ERROR_NULL("failed to allocate space for temporary spline vector", GSL_ENOMEM); + GSL_ERROR_NULL( + "failed to allocate space for temporary spline vector", + GSL_ENOMEM); } w->dB = gsl_matrix_alloc(k, k + 1); if (w->dB == 0) { gsl_bspline_free(w); - GSL_ERROR_NULL("failed to allocate space for temporary derivative matrix", GSL_ENOMEM); + GSL_ERROR_NULL( + "failed to allocate space for temporary derivative matrix", + GSL_ENOMEM); } - return (w); } } /* gsl_bspline_alloc() */ +/* +gsl_bspline_deriv_alloc() + Allocate space for a bspline derivative workspace. The size of the +workspace is O(2k^2) + +Inputs: bspline - gsl_bspline_workspace of the associated bspline + +Return: pointer to workspace +*/ + +gsl_bspline_deriv_workspace * +gsl_bspline_deriv_alloc(gsl_bspline_workspace * bspline) +{ + if (bspline == 0) + { + GSL_ERROR_NULL("bspline must not be null", GSL_EINVAL); + } + else if (bspline->k == 0) + { + /* Not strictly necessary, but a sanity check on the bspline */ + GSL_ERROR_NULL("k must be at least 1", GSL_EINVAL); + } + else + { + gsl_bspline_deriv_workspace *w; + + w = (gsl_bspline_deriv_workspace *) calloc(1, + sizeof(gsl_bspline_deriv_workspace)); + if (w == 0) + { + GSL_ERROR_NULL("failed to allocate space for workspace", GSL_ENOMEM); + } + + w->A = gsl_matrix_alloc(bspline->k, bspline->k); + if (w->A == 0) + { + gsl_bspline_deriv_free(w); + GSL_ERROR_NULL("failed to allocate space for derivative work matrix", + GSL_ENOMEM); + } + + w->dB = gsl_matrix_alloc(bspline->k, bspline->k + 1); + if (w->dB == 0) + { + gsl_bspline_deriv_free(w); + GSL_ERROR_NULL( + "failed to allocate space for temporary derivative matrix", + GSL_ENOMEM); + } + + w->bspline = bspline; + + return (w); + } +} /* gsl_bspline_deriv_alloc() */ + /* Return number of coefficients */ size_t -gsl_bspline_ncoeffs (gsl_bspline_workspace * w) +gsl_bspline_ncoeffs(gsl_bspline_workspace * w) { return w->n; } /* Return order */ size_t -gsl_bspline_order (gsl_bspline_workspace * w) +gsl_bspline_order(gsl_bspline_workspace * w) { return w->k; } /* Return number of breakpoints */ size_t -gsl_bspline_nbreak (gsl_bspline_workspace * w) +gsl_bspline_nbreak(gsl_bspline_workspace * w) { return w->nbreak; } +/* Return the location of the i-th breakpoint*/ double -gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w) +gsl_bspline_breakpoint(size_t i, gsl_bspline_workspace * w) { size_t j = i + w->k - 1; return gsl_vector_get(w->knots, j); @@ -175,7 +237,8 @@ gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w) /* gsl_bspline_free() - Free a bspline workspace + Free a gsl_bspline_workspace. +Any associated gsl_bspline_deriv_workspace should be freed beforehand. Inputs: w - workspace to free @@ -194,28 +257,46 @@ gsl_bspline_free(gsl_bspline_workspace *w) if (w->deltar) gsl_vector_free(w->deltar); - if (w->A) - gsl_matrix_free(w->A); - if (w->B) gsl_vector_free(w->B); + free(w); +} /* gsl_bspline_free() */ + +/* +gsl_bspline_deriv_free() + Free a gsl_bspline_deriv_workspace. +The associated gsl_bspline_workspace should be freed afterwards. + +Inputs: w - workspace to free + +Return: none +*/ + +void +gsl_bspline_deriv_free(gsl_bspline_deriv_workspace *w) +{ + w->bspline = 0; + + if (w->A) + gsl_matrix_free(w->A); + if (w->dB) gsl_matrix_free(w->dB); free(w); -} /* gsl_bspline_free() */ +} /* gsl_bspline_deriv_free() */ /* gsl_bspline_knots() Compute the knots from the given breakpoints: -knots(1:k) = breakpts(1) -knots(k+1:k+l-1) = breakpts(i), i = 2 .. l -knots(n+1:n+k) = breakpts(l + 1) + knots(1:k) = breakpts(1) + knots(k+1:k+l-1) = breakpts(i), i = 2 .. l + knots(n+1:n+k) = breakpts(l + 1) where l is the number of polynomial pieces (l = nbreak - 1) and -n = k + l - 1 + n = k + l - 1 (using matlab syntax for the arrays) The repeated knots at the beginning and end of the interval @@ -244,8 +325,7 @@ gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w) for (i = 1; i < w->l; ++i) { - gsl_vector_set(w->knots, w->k - 1 + i, - gsl_vector_get(breakpts, i)); + gsl_vector_set(w->knots, w->k - 1 + i, gsl_vector_get(breakpts, i)); } for (i = w->n; i < w->n + w->k; ++i) @@ -323,8 +403,7 @@ Notes: The w->knots vector must be initialized prior to calling */ int -gsl_bspline_eval(const double x, gsl_vector *B, - gsl_bspline_workspace *w) +gsl_bspline_eval(const double x, gsl_vector *B, gsl_bspline_workspace *w) { if (B->size != w->n) { @@ -358,7 +437,6 @@ gsl_bspline_eval(const double x, gsl_vector *B, } } /* gsl_bspline_eval() */ - /* gsl_bspline_eval_nonzero() Evaluate all non-zero B-spline functions at point x. @@ -368,7 +446,7 @@ Always B_i(x) = 0 for i < istart and for i > iend. Inputs: x - point at which to evaluate splines Bk - (output) where to store B-spline values (length k) istart - (output) B-spline function index of - first non-zero basis for given x + first non-zero basis for given x iend - (output) B-spline function index of last non-zero basis for given x. This is also the knot index corresponding to x. @@ -381,18 +459,15 @@ Notes: 1) the w->knots vector must be initialized before calling 2) On output, B contains - [B_{istart,k}, B_{istart+1,k}, - ..., B_{iend-1,k}, B_{iend,k}] + [B_{istart,k}, B_{istart+1,k}, + ..., B_{iend-1,k}, B_{iend,k}] evaluated at the given x. */ int -gsl_bspline_eval_nonzero(const double x, - gsl_vector *Bk, - size_t *istart, - size_t *iend, - gsl_bspline_workspace *w) +gsl_bspline_eval_nonzero(const double x, gsl_vector *Bk, size_t *istart, + size_t *iend, gsl_bspline_workspace *w) { if (Bk->size != w->k) { @@ -400,9 +475,9 @@ gsl_bspline_eval_nonzero(const double x, } else { - size_t i; /* spline index */ - size_t j; /* looping */ - int flag = 0; /* interval search flag */ + size_t i; /* spline index */ + size_t j; /* looping */ + int flag = 0; /* interval search flag */ int error = 0; /* error flag */ i = bspline_find_interval(x, &flag, w); @@ -415,27 +490,26 @@ gsl_bspline_eval_nonzero(const double x, *istart = i - w->k + 1; *iend = i; - bspline_pppack_bsplvb(w->knots, w->k, 1, x, *iend, - &j, w->deltal, w->deltar, Bk); + bspline_pppack_bsplvb(w->knots, w->k, 1, x, *iend, &j, w->deltal, + w->deltar, Bk); return GSL_SUCCESS; } } /* gsl_bspline_eval_nonzero() */ /* -gsl_bspline_eval_deriv() +gsl_bspline_deriv_eval() Evaluate d^j/dx^j B_i(x) for all i, 0 <= j <= nderiv. -This is a wrapper function for gsl_bspline_eval_deriv_nonzero() +This is a wrapper function for gsl_bspline_deriv_eval_nonzero() which formats the output in a nice way. Inputs: x - point for evaluation - nderiv - number of derivatives to compute, - inclusive. + nderiv - number of derivatives to compute, inclusive. dB - (output) where to store d^j/dx^j B_i(x) values. the size of this matrix is (n = nbreak + k - 2 = l + k - 1 = w->n) by (nderiv + 1) - w - bspline workspace + w - bspline derivative workspace Return: success or error @@ -446,30 +520,29 @@ Notes: 1) The w->knots vector must be initialized prior to calling */ int -gsl_bspline_eval_deriv(const double x, - const size_t nderiv, - gsl_matrix *dB, - gsl_bspline_workspace *w) +gsl_bspline_deriv_eval(const double x, const size_t nderiv, gsl_matrix *dB, + gsl_bspline_deriv_workspace *w) { - if (dB->size1 != w->n) + if (dB->size1 != w->bspline->n) { GSL_ERROR("dB matrix first dimension not of length n", GSL_EBADLEN); } - else if (dB->size2 < nderiv+1) + else if (dB->size2 < nderiv + 1) { - GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", GSL_EBADLEN); + GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", + GSL_EBADLEN); } else { - size_t i; /* looping */ - size_t j; /* looping */ - size_t istart; /* first non-zero spline for x */ - size_t iend; /* last non-zero spline for x, knot for x*/ - int error; /* error handling */ + size_t i; /* looping */ + size_t j; /* looping */ + size_t istart; /* first non-zero spline for x */ + size_t iend; /* last non-zero spline for x, knot for x*/ + int error; /* error handling */ /* find all non-zero d^j/dx^j B_i(x) values */ - error = gsl_bspline_eval_deriv_nonzero( - x, nderiv, w->dB, &istart, &iend, w); + error = gsl_bspline_deriv_eval_nonzero(x, nderiv, w->dB, &istart, &iend, + w); if (error) { return error; @@ -484,16 +557,16 @@ gsl_bspline_eval_deriv(const double x, for (i = istart; i <= iend; ++i) gsl_matrix_set(dB, i, j, gsl_matrix_get(w->dB, i - istart, j)); - for (i = iend + 1; i < w->n; ++i) + for (i = iend + 1; i < w->bspline->n; ++i) gsl_matrix_set(dB, i, j, 0.0); } return GSL_SUCCESS; } -} /* gsl_bspline_eval_deriv() */ +} /* gsl_bspline_deriv_eval() */ /* -gsl_bspline_eval_deriv_nonzero() +gsl_bspline_deriv_eval_nonzero() At point x evaluate all requested, non-zero B-spline function derivatives and store them in dB. These are the d^j/dx^j B_i(x) with i in [istart, iend] and j in [0, nderiv]. @@ -508,7 +581,7 @@ Inputs: x - point at which to evaluate splines iend - (output) B-spline function index of last non-zero basis for given x. This is also the knot index corresponding to x. - w - bspline workspace + w - bspline derivative workspace Return: success or error @@ -517,14 +590,14 @@ Notes: 1) the w->knots vector must be initialized before calling 2) On output, dB contains - [[B_{istart, k}, ..., d^nderiv/dx^nderiv B_{istart ,k}], - [B_{istart+1,k}, ..., d^nderiv/dx^nderiv B_{istart+1,k}], - ... - [B_{iend-1, k}, ..., d^nderiv/dx^nderiv B_{iend-1, k}], - [B_{iend, k}, ..., d^nderiv/dx^nderiv B_{iend, k}]] + [[B_{istart, k}, ..., d^nderiv/dx^nderiv B_{istart ,k}], + [B_{istart+1,k}, ..., d^nderiv/dx^nderiv B_{istart+1,k}], + ... + [B_{iend-1, k}, ..., d^nderiv/dx^nderiv B_{iend-1, k}], + [B_{iend, k}, ..., d^nderiv/dx^nderiv B_{iend, k}]] - evaluated at x. B_{istart, k} is stored in dB(0,0). - Each additional column contains an additional derivative. + evaluated at x. B_{istart, k} is stored in dB(0,0). + Each additional column contains an additional derivative. 3) Note that the zero-th column of the result contains the 0th derivative, which is simply a function evaluation. @@ -533,48 +606,46 @@ Notes: 1) the w->knots vector must be initialized before calling */ int -gsl_bspline_eval_deriv_nonzero(const double x, - const size_t nderiv, - gsl_matrix *dB, - size_t *istart, - size_t *iend, - gsl_bspline_workspace *w) +gsl_bspline_deriv_eval_nonzero(const double x, const size_t nderiv, + gsl_matrix *dB, size_t *istart, size_t *iend, + gsl_bspline_deriv_workspace *w) { - if (dB->size1 != w->k) + if (dB->size1 != w->bspline->k) { GSL_ERROR("dB matrix first dimension not of length k", GSL_EBADLEN); } - else if (dB->size2 < nderiv+1) + else if (dB->size2 < nderiv + 1) { - GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", GSL_EBADLEN); + GSL_ERROR("dB matrix second dimension must be at least length nderiv+1", + GSL_EBADLEN); } else { size_t i; /* spline index */ size_t j; /* looping */ - int flag = 0; /* interval search flag */ + int flag = 0; /* interval search flag */ int error = 0; /* error flag */ size_t min_nderivk; - i = bspline_find_interval(x, &flag, w); - error = bspline_process_interval_for_eval(x, &i, flag, w); + i = bspline_find_interval(x, &flag, w->bspline); + error = bspline_process_interval_for_eval(x, &i, flag, w->bspline); if (error) { return error; } - *istart = i - w->k + 1; + *istart = i - w->bspline->k + 1; *iend = i; - bspline_pppack_bsplvd(w->knots, w->k, x, *iend, - w->deltal, w->deltar, w->A, dB, nderiv); + bspline_pppack_bsplvd(w->bspline->knots, w->bspline->k, x, *iend, + w->bspline->deltal, w->bspline->deltar, w->A, dB, nderiv); /* An order k b-spline has at most k-1 nonzero derivatives - so we need to zero all requested higher order derivatives */ - min_nderivk = GSL_MIN_INT(nderiv, w->k - 1); + so we need to zero all requested higher order derivatives */ + min_nderivk = GSL_MIN_INT(nderiv, w->bspline->k - 1); for (j = min_nderivk + 1; j <= nderiv; ++j) { - for (i = 0; i < w->k; ++i) + for (i = 0; i < w->bspline->k; ++i) { gsl_matrix_set(dB, i, j, 0.0); } @@ -582,8 +653,7 @@ gsl_bspline_eval_deriv_nonzero(const double x, return GSL_SUCCESS; } -} /* gsl_bspline_eval_deriv_nonzero() */ - +} /* gsl_bspline_deriv_eval_nonzero() */ /**************************************** * INTERNAL ROUTINES * @@ -591,10 +661,7 @@ gsl_bspline_eval_deriv_nonzero(const double x, /* bspline_find_interval() - Find knot interval such that - - t_i <= x < t_{i + 1} - + Find knot interval such that t_i <= x < t_{i + 1} where the t_i are knot values. Inputs: x - x value @@ -605,17 +672,16 @@ Return: i (index in w->knots corresponding to left limit of interval) Notes: The error conditions are reported as follows: -Condition Return value Flag ---------- ------------ ---- -x < t_0 0 -1 -t_i <= x < t_{i+1} i 0 -t_i < x = t_{i+1} = t_{n+k-1} i 0 -t_{n+k-1} < x l+k-1 +1 + Condition Return value Flag + --------- ------------ ---- + x < t_0 0 -1 + t_i <= x < t_{i+1} i 0 + t_i < x = t_{i+1} = t_{n+k-1} i 0 + t_{n+k-1} < x l+k-1 +1 */ static inline size_t -bspline_find_interval(const double x, int *flag, - gsl_bspline_workspace *w) +bspline_find_interval(const double x, int *flag, gsl_bspline_workspace *w) { size_t i; @@ -639,8 +705,8 @@ bspline_find_interval(const double x, int *flag, if (ti <= x && x < tip1) break; - if (ti < x && x == tip1 && - tip1 == gsl_vector_get(w->knots, w->k + w->l - 1)) + if (ti < x && x == tip1 && tip1 == gsl_vector_get(w->knots, w->k + w->l + - 1)) break; } @@ -652,17 +718,16 @@ bspline_find_interval(const double x, int *flag, return i; } /* bspline_find_interval() */ - /* bspline_process_interval_for_eval() - Consumes an x location, left knot from bspline_find_interval, flag + Consumes an x location, left knot from bspline_find_interval, flag from bspline_find_interval, and a workspace. Checks that x lies within the splines' knots, enforces some endpoint continuity requirements, and avoids divide by zero errors in the underlying bspline_pppack_* functions. - */ -static inline int bspline_process_interval_for_eval( - const double x, size_t *i, const int flag, gsl_bspline_workspace - *w) +*/ +static inline int +bspline_process_interval_for_eval(const double x, size_t *i, const int flag, + gsl_bspline_workspace *w) { if (flag == -1) { @@ -689,7 +754,6 @@ static inline int bspline_process_interval_for_eval( return GSL_SUCCESS; } /* bspline_process_interval_for_eval */ - /******************************************************************** * PPPACK ROUTINES * @@ -705,66 +769,71 @@ bspline_pppack_bsplvb() jout = max( jhigh , (j+1)*(index-1) ) with knot sequence t. Parameters: - t - knot sequence, of length left + jout , assumed to be - nondecreasing. assumption t(left).lt.t(left + 1). - division by zero will result if t(left) = t(left+1) - jhigh - - index - integers which determine the order jout = max(jhigh, - (j+1)*(index-1)) of the b-splines whose values at x - are to be returned. index is used to avoid - recalculations when several columns of the triangular - array of b-spline values are needed (e.g., in bsplpp - or in bsplvd ). - precisely, if index = 1 , - the calculation starts from scratch and the entire - triangular array of b-spline values of orders - 1,2,...,jhigh is generated order by order , i.e., - column by column . - if index = 2 , - only the b-spline values of order j+1, j+2, ..., jout - are generated, the assumption being that biatx, j, - deltal, deltar are, on entry, as they were on exit - at the previous call. - in particular, if jhigh = 0, then jout = j+1, i.e., - just the next column of b-spline values is generated. - x - the point at which the b-splines are to be evaluated. - left - an integer chosen (usually) so that - t(left) .le. x .le. t(left+1). - j - (output) a working scalar for indexing - deltal - (output) a working area which must be of length at - least jout - deltar - (output) a working area which must be of length at - least jout - biatx - (output) array of length jout, with biatx(i) - containing the value at x of the polynomial of order - jout which agrees with the b-spline b(left-jout+i,jout,t) - on the interval (t(left), t(left+1)) . + t - knot sequence, of length left + jout , assumed to be + nondecreasing. assumption t(left).lt.t(left + 1). + division by zero will result if t(left) = t(left+1) + jhigh - + index - integers which determine the order jout = max(jhigh, + (j+1)*(index-1)) of the b-splines whose values at x + are to be returned. index is used to avoid + recalculations when several columns of the triangular + array of b-spline values are needed (e.g., in bsplpp + or in bsplvd ). precisely, + + if index = 1 , + the calculation starts from scratch and the entire + triangular array of b-spline values of orders + 1,2,...,jhigh is generated order by order , i.e., + column by column . + + if index = 2 , + only the b-spline values of order j+1, j+2, ..., jout + are generated, the assumption being that biatx, j, + deltal, deltar are, on entry, as they were on exit + at the previous call. + + in particular, if jhigh = 0, then jout = j+1, i.e., + just the next column of b-spline values is generated. + x - the point at which the b-splines are to be evaluated. + left - an integer chosen (usually) so that + t(left) .le. x .le. t(left+1). + j - (output) a working scalar for indexing + deltal - (output) a working area which must be of length at least jout + deltar - (output) a working area which must be of length at least jout + biatx - (output) array of length jout, with biatx(i) + containing the value at x of the polynomial of order + jout which agrees with the b-spline b(left-jout+i,jout,t) + on the interval (t(left), t(left+1)) . Method: - the recurrence relation - - x - t(i) t(i+j+1) - x - b(i,j+1)(x) = -----------b(i,j)(x) + ---------------b(i+1,j)(x) - t(i+j)-t(i) t(i+j+1)-t(i+1) - - is used (repeatedly) to generate the (j+1)-vector b(left-j,j+1)(x), - ...,b(left,j+1)(x) from the j-vector b(left-j+1,j)(x),..., - b(left,j)(x), storing the new values in biatx over the old. the - facts that - b(i,1) = 1 if t(i) .le. x .lt. t(i+1) - and that - b(i,j)(x) = 0 unless t(i) .le. x .lt. t(i+j) - are used. the particular organization of the calculations follows - algorithm (8) in chapter x of [1]. + the recurrence relation + + x - t(i) t(i+j+1) - x + b(i,j+1)(x) = -----------b(i,j)(x) + ---------------b(i+1,j)(x) + t(i+j)-t(i) t(i+j+1)-t(i+1) + + is used (repeatedly) to generate the (j+1)-vector b(left-j,j+1)(x), + ...,b(left,j+1)(x) from the j-vector b(left-j+1,j)(x),..., + b(left,j)(x), storing the new values in biatx over the old. the + facts that + + b(i,1) = 1 if t(i) .le. x .lt. t(i+1) + + and that + + b(i,j)(x) = 0 unless t(i) .le. x .lt. t(i+j) + + are used. the particular organization of the calculations follows + algorithm (8) in chapter x of [1]. Notes: - (1) This is a direct translation of PPPACK's bsplvb routine with - j, deltal, deltar rewritten as input parameters and - utilizing zero-based indexing. + (1) This is a direct translation of PPPACK's bsplvb routine with + j, deltal, deltar rewritten as input parameters and + utilizing zero-based indexing. - (2) This routine contains no error checking. Please use routines - like gsl_bspline_eval(). + (2) This routine contains no error checking. Please use routines + like gsl_bspline_eval(). */ static void @@ -778,7 +847,7 @@ bspline_pppack_bsplvb(const gsl_vector *t, gsl_vector *deltar, gsl_vector *biatx) { - size_t i; /* looping */ + size_t i; /* looping */ double saved; double term; @@ -797,12 +866,10 @@ bspline_pppack_bsplvb(const gsl_vector *t, for (i = 0; i <= *j; ++i) { - term = gsl_vector_get(biatx, i) / - (gsl_vector_get(deltar, i) + - gsl_vector_get(deltal, *j - i)); + term = gsl_vector_get(biatx, i) / (gsl_vector_get(deltar, i) + + gsl_vector_get(deltal, *j - i)); - gsl_vector_set(biatx, i, - saved + gsl_vector_get(deltar, i) * term); + gsl_vector_set(biatx, i, saved + gsl_vector_get(deltar, i) * term); saved = gsl_vector_get(deltal, *j - i) * term; } @@ -815,48 +882,47 @@ bspline_pppack_bsplvb(const gsl_vector *t, /* bspline_pppack_bsplvd() -calculates value and derivs of all b-splines which do not vanish at x + calculates value and derivs of all b-splines which do not vanish at x Parameters: - t - the knot array, of length left+k (at least) - k - the order of the b-splines to be evaluated - x - the point at which these values are sought - left - an integer indicating the left endpoint of the interval - of interest. the k b-splines whose support contains the - interval (t(left), t(left+1)) are to be considered. - it is assumed that t(left) .lt. t(left+1) - division by zero will result otherwise (in bsplvb). - also, the output is as advertised only if - t(left) .le. x .le. t(left+1) . - deltal - a working area which must be of length at least k - deltar - a working area which must be of length at least k - a - an array of order (k,k), to contain b-coeffs of the - derivatives of a certain order of the k b-splines - of interest. - dbiatx - an array of order (k,nderiv). its entry (i,m) contains - value of (m)th derivative of (left-k+i)-th b-spline - of order k for knot sequence t, i=1,...,k, - m=0,...,nderiv. - nderiv - an integer indicating that values of b-splines and - their derivatives up to AND INCLUDING the nderiv-th - are asked for. (nderiv is replaced internally by the - integer mhigh in (1,k) closest to it.) + t - the knot array, of length left+k (at least) + k - the order of the b-splines to be evaluated + x - the point at which these values are sought + left - an integer indicating the left endpoint of the interval + of interest. the k b-splines whose support contains the + interval (t(left), t(left+1)) are to be considered. + it is assumed that t(left) .lt. t(left+1) + division by zero will result otherwise (in bsplvb). + also, the output is as advertised only if + t(left) .le. x .le. t(left+1) . + deltal - a working area which must be of length at least k + deltar - a working area which must be of length at least k + a - an array of order (k,k), to contain b-coeffs of the + derivatives of a certain order of the k b-splines + of interest. + dbiatx - an array of order (k,nderiv). its entry (i,m) contains + value of (m)th derivative of (left-k+i)-th b-spline + of order k for knot sequence t, i=1,...,k, m=0,...,nderiv. + nderiv - an integer indicating that values of b-splines and + their derivatives up to AND INCLUDING the nderiv-th + are asked for. (nderiv is replaced internally by the + integer mhigh in (1,k) closest to it.) Method: - values at x of all the relevant b-splines of order k,k-1,..., k+1-nderiv - are generated via bsplvb and stored temporarily in dbiatx. then, the - b-coeffs of the required derivatives of the b-splines of interest are - generated by differencing, each from the preceeding one of lower order, - and combined with the values of b-splines of corresponding order in - dbiatx to produce the desired values . + values at x of all the relevant b-splines of order k,k-1,..., k+1-nderiv + are generated via bsplvb and stored temporarily in dbiatx. then, the + b-coeffs of the required derivatives of the b-splines of interest are + generated by differencing, each from the preceeding one of lower order, + and combined with the values of b-splines of corresponding order in + dbiatx to produce the desired values . Notes: - (1) This is a direct translation of PPPACK's bsplvd routine with - deltal, deltar rewritten as input parameters (to later feed them - to bspline_pppack_bsplvb) and utilizing zero-based indexing. + (1) This is a direct translation of PPPACK's bsplvd routine with + deltal, deltar rewritten as input parameters (to later feed them + to bspline_pppack_bsplvb) and utilizing zero-based indexing. - (2) This routine contains no error checking. + (2) This routine contains no error checking. */ static void @@ -877,36 +943,36 @@ bspline_pppack_bsplvd(const gsl_vector *t, gsl_vector_view dbcol = gsl_matrix_column(dbiatx, 0); mhigh = GSL_MIN_INT(nderiv, k - 1); - bspline_pppack_bsplvb(t, k-mhigh, 1, x, left, - &bsplvb_j, deltal, deltar, &dbcol.vector); + bspline_pppack_bsplvb(t, k - mhigh, 1, x, left, &bsplvb_j, deltal, deltar, + &dbcol.vector); if (mhigh > 0) { /* the first column of dbiatx always contains the b-spline - values for the current order. these are stored in column - k-current order before bsplvb is called to put values - for the next higher order on top of it. */ + values for the current order. these are stored in column + k-current order before bsplvb is called to put values + for the next higher order on top of it. */ ideriv = mhigh; for (m = 1; m <= mhigh; ++m) { - for (j = ideriv, jp1mid = 0; j < (int)k; ++j, ++jp1mid) + for (j = ideriv, jp1mid = 0; j < (int) k; ++j, ++jp1mid) { - gsl_matrix_set(dbiatx, j, ideriv, - gsl_matrix_get(dbiatx, jp1mid, 0)); + gsl_matrix_set(dbiatx, j, ideriv, gsl_matrix_get(dbiatx, jp1mid, + 0)); } --ideriv; - bspline_pppack_bsplvb(t, k-ideriv, 2, x, left, - &bsplvb_j, deltal, deltar, &dbcol.vector); + bspline_pppack_bsplvb(t, k - ideriv, 2, x, left, &bsplvb_j, deltal, + deltar, &dbcol.vector); } /* at this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) - for i=j,...,k-1 and j=0,...,mhigh. in particular, the - first column of dbiatx is already in final form. to obtain - corresponding derivatives of b-splines in subsequent columns, - generate their b-repr. by differencing, then evaluate at x. */ + for i=j,...,k-1 and j=0,...,mhigh. in particular, the + first column of dbiatx is already in final form. to obtain + corresponding derivatives of b-splines in subsequent columns, + generate their b-repr. by differencing, then evaluate at x. */ jlow = 0; - for (i = 0; i < (int)k; ++i) + for (i = 0; i < (int) k; ++i) { - for (j = jlow; j < (int)k; ++j) + for (j = jlow; j < (int) k; ++j) { gsl_matrix_set(a, j, i, 0.0); } @@ -915,7 +981,7 @@ bspline_pppack_bsplvd(const gsl_vector *t, } /* at this point, a(.,j) contains the b-coeffs for the j-th of the - k b-splines of interest here. */ + k b-splines of interest here. */ for (m = 1; m <= mhigh; ++m) { kmm = k - m; @@ -924,40 +990,38 @@ bspline_pppack_bsplvd(const gsl_vector *t, i = k - 1; /* for j=1,...,k, construct b-coeffs of (m)th derivative - of b-splines from those for preceding derivative by - differencing and store again in a(.,j) . the fact that - a(i,j) = 0 for i .lt. j is used. */ + of b-splines from those for preceding derivative by + differencing and store again in a(.,j) . the fact that + a(i,j) = 0 for i .lt. j is used. */ for (ldummy = 0; ldummy < kmm; ++ldummy) { - factor = fkmm/( - gsl_vector_get(t, il + kmm) - gsl_vector_get(t, il) - ); + factor = fkmm / (gsl_vector_get(t, il + kmm) - gsl_vector_get(t, + il)); /* the assumption that t(left).lt.t(left+1) makes - denominator in factor nonzero. */ - for (j = 0; j <=i; ++j) + denominator in factor nonzero. */ + for (j = 0; j <= i; ++j) { - gsl_matrix_set(a, i, j, factor*( - gsl_matrix_get(a, i, j) - gsl_matrix_get(a, i-1, j))); + gsl_matrix_set(a, i, j, factor * (gsl_matrix_get(a, i, j) + - gsl_matrix_get(a, i - 1, j))); } --il; --i; } /* for i=1,...,k, combine b-coeffs a(.,i) with b-spline values - stored in dbiatx(.,m) to get value of (m)th derivative - of i-th b-spline (of interest here) at x, and store in - dbiatx(i,m). storage of this value over the value of a - b-spline of order m there is safe since the remaining - b-spline derivatives of the same order do not use this - value due to the fact that a(j,i) = 0 for j .lt. i . */ - for (i = 0; i < (int)k; ++i) + stored in dbiatx(.,m) to get value of (m)th derivative + of i-th b-spline (of interest here) at x, and store in + dbiatx(i,m). storage of this value over the value of a + b-spline of order m there is safe since the remaining + b-spline derivatives of the same order do not use this + value due to the fact that a(j,i) = 0 for j .lt. i . */ + for (i = 0; i < (int) k; ++i) { sum = 0; jlow = GSL_MAX_INT(i, m); - for (j = jlow; j < (int)k; ++j) + for (j = jlow; j < (int) k; ++j) { - sum += gsl_matrix_get(a, j, i) - * gsl_matrix_get(dbiatx, j, m); + sum += gsl_matrix_get(a, j, i) * gsl_matrix_get(dbiatx, j, m); } gsl_matrix_set(dbiatx, i, m, sum); } diff --git a/bspline/gsl_bspline.h b/bspline/gsl_bspline.h index 8a2ffc5..f7f8175 100644 --- a/bspline/gsl_bspline.h +++ b/bspline/gsl_bspline.h @@ -39,41 +39,53 @@ __BEGIN_DECLS typedef struct - { - size_t k; /* spline order */ - size_t km1; /* k - 1 (polynomial order) */ - size_t l; /* number of polynomial pieces on interval */ +{ + size_t k; /* spline order */ + size_t km1; /* k - 1 (polynomial order) */ + size_t l; /* number of polynomial pieces on interval */ size_t nbreak; /* number of breakpoints (l + 1) */ - size_t n; /* number of bspline basis functions (l + k - 1) */ + size_t n; /* number of bspline basis functions (l + k - 1) */ - gsl_vector *knots; /* knots vector */ + gsl_vector *knots; /* knots vector */ gsl_vector *deltal; /* left delta */ gsl_vector *deltar; /* right delta */ - gsl_matrix *A; /* work matrix */ - gsl_vector *B; /* temporary spline results */ - gsl_matrix *dB; /* temporary derivative results */ - } gsl_bspline_workspace; + gsl_matrix *A; /* work matrix */ + gsl_vector *B; /* temporary spline results */ + gsl_matrix *dB; /* temporary derivative results */ +} gsl_bspline_workspace; + +typedef struct +{ + gsl_bspline_workspace *bspline; /* associated bspline */ + gsl_matrix *A; /* work matrix */ + gsl_matrix *dB; /* temporary derivative results */ +} gsl_bspline_deriv_workspace; gsl_bspline_workspace * gsl_bspline_alloc(const size_t k, const size_t nbreak); -void gsl_bspline_free(gsl_bspline_workspace *w); +void +gsl_bspline_free(gsl_bspline_workspace *w); -size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * w); -size_t gsl_bspline_order (gsl_bspline_workspace * w); -size_t gsl_bspline_nbreak (gsl_bspline_workspace * w); -double gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w); +size_t +gsl_bspline_ncoeffs(gsl_bspline_workspace * w); +size_t +gsl_bspline_order(gsl_bspline_workspace * w); +size_t +gsl_bspline_nbreak(gsl_bspline_workspace * w); +double +gsl_bspline_breakpoint(size_t i, gsl_bspline_workspace * w); int gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w); -int gsl_bspline_knots_uniform(const double a, const double b, - gsl_bspline_workspace *w); +int +gsl_bspline_knots_uniform(const double a, + const double b, + gsl_bspline_workspace *w); int -gsl_bspline_eval(const double x, - gsl_vector *B, - gsl_bspline_workspace *w); +gsl_bspline_eval(const double x, gsl_vector *B, gsl_bspline_workspace *w); int gsl_bspline_eval_nonzero(const double x, @@ -82,19 +94,25 @@ gsl_bspline_eval_nonzero(const double x, size_t *iend, gsl_bspline_workspace *w); +gsl_bspline_deriv_workspace * +gsl_bspline_deriv_alloc(gsl_bspline_workspace * bspline); + +void +gsl_bspline_deriv_free(gsl_bspline_deriv_workspace *w); + int -gsl_bspline_eval_deriv(const double x, +gsl_bspline_deriv_eval(const double x, const size_t nderiv, gsl_matrix *dB, - gsl_bspline_workspace *w); + gsl_bspline_deriv_workspace *w); int -gsl_bspline_eval_deriv_nonzero(const double x, +gsl_bspline_deriv_eval_nonzero(const double x, const size_t nderiv, gsl_matrix *dB, size_t *istart, size_t *iend, - gsl_bspline_workspace *w); + gsl_bspline_deriv_workspace *w); __END_DECLS diff --git a/bspline/test.c b/bspline/test.c index 393fbf3..5b511d3 100644 --- a/bspline/test.c +++ b/bspline/test.c @@ -26,7 +26,7 @@ #include void -test_bspline(gsl_bspline_workspace * bw) +test_bspline(gsl_bspline_workspace * bw, gsl_bspline_deriv_workspace * dbw) { gsl_vector *B; gsl_matrix *dB; @@ -68,7 +68,7 @@ test_bspline(gsl_bspline_workspace * bw) { double xi = a + (b - a) * (i / (n - 1.0)); gsl_bspline_eval(xi, B, bw); - gsl_bspline_eval_deriv(xi, 0, dB, bw); + gsl_bspline_deriv_eval(xi, 0, dB, dbw); for (j = 0; j < ncoeffs; j++) { @@ -100,8 +100,10 @@ main(int argc, char **argv) { double a = -1.23 * order, b = 45.6 * order; gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints); + gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw); gsl_bspline_knots_uniform(a, b, bw); - test_bspline(bw); + test_bspline(bw, dbw); + gsl_bspline_deriv_free(dbw); gsl_bspline_free(bw); } } @@ -113,6 +115,7 @@ main(int argc, char **argv) { double a = -1.23 * order, b = 45.6 * order; gsl_bspline_workspace *bw = gsl_bspline_alloc(order, breakpoints); + gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw); gsl_vector *k = gsl_vector_alloc(breakpoints); for (i = 0; i < breakpoints; i++) { @@ -122,8 +125,9 @@ main(int argc, char **argv) gsl_vector_set(k, i, x); }; gsl_bspline_knots(k, bw); - test_bspline(bw); + test_bspline(bw, dbw); gsl_vector_free(k); + gsl_bspline_deriv_free(dbw); gsl_bspline_free(bw); } } @@ -143,6 +147,7 @@ main(int argc, char **argv) }; gsl_bspline_workspace *bw = gsl_bspline_alloc(2, 3); + gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw); gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw), gsl_bspline_order(bw) + 1); @@ -158,7 +163,7 @@ main(int argc, char **argv) /* Initialize dB with poison to ensure we overwrite it */ gsl_matrix_set_all(dB, GSL_NAN); - gsl_bspline_eval_deriv(xloc[i], gsl_bspline_order(bw), dB, bw); + gsl_bspline_deriv_eval(xloc[i], gsl_bspline_order(bw), dB, dbw); for (j = 0; j < gsl_bspline_ncoeffs(bw) ; ++j) { /* check basis function 1st deriv */ @@ -178,6 +183,7 @@ main(int argc, char **argv) } gsl_matrix_free(dB); + gsl_bspline_deriv_free(dbw); gsl_bspline_free(bw); gsl_vector_free(breakpts); } @@ -213,6 +219,7 @@ main(int argc, char **argv) }; gsl_bspline_workspace *bw = gsl_bspline_alloc(3, 5); + gsl_bspline_deriv_workspace *dbw = gsl_bspline_deriv_alloc(bw); gsl_matrix *dB = gsl_matrix_alloc(gsl_bspline_ncoeffs(bw), gsl_bspline_order(bw) + 1); @@ -229,7 +236,7 @@ main(int argc, char **argv) { /* Initialize dB with poison to ensure we overwrite it */ gsl_matrix_set_all(dB, GSL_NAN); - gsl_bspline_eval_deriv(xloc[i], gsl_bspline_order(bw), dB, bw); + gsl_bspline_deriv_eval(xloc[i], gsl_bspline_order(bw), dB, dbw); /* check basis function evaluation */ for (j = 0; j < gsl_bspline_ncoeffs(bw); ++j) @@ -255,6 +262,7 @@ main(int argc, char **argv) } gsl_matrix_free(dB); + gsl_bspline_deriv_free(dbw); gsl_bspline_free(bw); gsl_vector_free(breakpts); } diff --git a/doc/bspline.texi b/doc/bspline.texi index 36b3dce..b2cdb1b 100644 --- a/doc/bspline.texi +++ b/doc/bspline.texi @@ -78,6 +78,8 @@ be readily obtained from a least-squares fit. @section Initializing the B-splines solver @cindex basis splines, initializing +Using B-splines requires a gsl_bspline_workspace: + @deftypefun {gsl_bspline_workspace *} gsl_bspline_alloc (const size_t @var{k}, const size_t @var{nbreak}) This function allocates a workspace for computing B-splines of order @var{k}. The number of breakpoints is given by @var{nbreak}. This @@ -88,6 +90,21 @@ are specified by @math{k = 4}. The size of the workspace is @deftypefun void gsl_bspline_free (gsl_bspline_workspace * @var{w}) This function frees the memory associated with the workspace @var{w}. +Any associated gsl_bspline_deriv_workspace should be freed beforehand. +@end deftypefun + +Evaluation of B-spline basis function derivatives additionally requires +a gsl_bspline_deriv_workspace: + +@deftypefun {gsl_bspline_deriv_workspace *} gsl_bspline_deriv_alloc (gsl_bspline_workspace * bspline) +This function allocates a workspace for computing the B-spline basis +function derivatives for @var{bspline}. The size of the workspace is +@math{O(2k^2)}. +@end deftypefun + +@deftypefun void gsl_bspline_deriv_free (gsl_bspline_deriv_workspace * @var{w}) +This function frees the memory associated with the workspace @var{w}. +The associated gsl_bspline_workspace should be freed afterwards. @end deftypefun @node Constructing the knots vector @@ -140,7 +157,7 @@ This function returns the number of B-spline coefficients given by @section Evaluation of B-spline derivatives @cindex basis splines, derivatives -@deftypefun int gsl_bspline_eval_deriv (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, gsl_bspline_workspace * @var{w}) +@deftypefun int gsl_bspline_deriv_eval (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, gsl_bspline_workspace * @var{w}) This function evaluates all B-spline basis function derivatives of orders @math{0} through @math{nderiv} (inclusive) at the position @var{x} and stores them in @var{dB}. The @math{(i,j)}th element of @var{dB} @@ -153,7 +170,7 @@ at once than to compute them individually, due to the nature of the defining recurrence relation. @end deftypefun -@deftypefun int gsl_bspline_eval_deriv_nonzero (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, size_t * @var{istart}, size_t * @var{iend}, gsl_bspline_workspace * @var{w}) +@deftypefun int gsl_bspline_deriv_eval_nonzero (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, size_t * @var{istart}, size_t * @var{iend}, gsl_bspline_workspace * @var{w}) This function evaluates all potentially nonzero B-spline basis function derivatives of orders @math{0} through @math{nderiv} (inclusive) at the position @var{x} and stores them in @var{dB}. The @math{(i,j)}th