From 7ac6ed29cc302f00228732d4000bd324a2a95241 Mon Sep 17 00:00:00 2001 From: Rhys Ulerich Date: Wed, 17 Jun 2009 16:27:27 -0500 Subject: [PATCH] Fixed Greville abscissae implementation Repairs the previous, broken implementation of Greville abscissae. I had misunderstood the definition because it was implied that the first and last knots are excluded prior to the averaging process. That changes both the number and values of the locations. Updates to unit tests and documentation included. --- bspline/bspline.c | 21 +++++++++-------- bspline/gsl_bspline.h | 1 - bspline/test.c | 56 ++++++++++++++++++++++++++++++++++++------------ doc/bspline.texi | 24 ++++++++++++++------ 4 files changed, 70 insertions(+), 32 deletions(-) diff --git a/bspline/bspline.c b/bspline/bspline.c index d1e4470..b538e01 100644 --- a/bspline/bspline.c +++ b/bspline/bspline.c @@ -208,27 +208,28 @@ gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w) return gsl_vector_get (w->knots, j); } -/* Return the number of Greville abscissae for this basis */ -size_t -gsl_bspline_greville_nabscissae(gsl_bspline_workspace *w) -{ - return w->knots->size - w->km1; -} - /* Return the location of the i-th Greville abscissa */ double gsl_bspline_greville_abscissa(size_t i, gsl_bspline_workspace *w) { #if GSL_RANGE_CHECK - if (GSL_RANGE_COND(i >= gsl_bspline_greville_nabscissae(w))) + if (GSL_RANGE_COND(i >= gsl_bspline_ncoeffs(w))) { GSL_ERROR_VAL ("Greville abscissa index out of range", GSL_EINVAL, 0); } #endif const size_t stride = w->knots->stride; - const double * data = w->knots->data + i*stride; + size_t km1 = w->km1; + double * data = w->knots->data + (i+1)*stride; + + if (km1 == 0) + { + /* Return interval midpoints in degenerate k = 1 case*/ + km1 = 2; + data -= stride; + } - return gsl_stats_mean(data, stride, w->k); + return gsl_stats_mean(data, stride, km1); } /* diff --git a/bspline/gsl_bspline.h b/bspline/gsl_bspline.h index f30fefa..7b89a10 100644 --- a/bspline/gsl_bspline.h +++ b/bspline/gsl_bspline.h @@ -68,7 +68,6 @@ 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_greville_nabscissae(gsl_bspline_workspace *w); double gsl_bspline_greville_abscissa(size_t i, gsl_bspline_workspace *w); int diff --git a/bspline/test.c b/bspline/test.c index 9256a77..885fb87 100644 --- a/bspline/test.c +++ b/bspline/test.c @@ -267,25 +267,53 @@ main(int argc, char **argv) gsl_vector_free(breakpts); } - /* Check Greville abscissae functionality on a uniform k=3 */ + /* Check Greville abscissae functionality on a non-uniform k=1 */ { size_t i; /* looping */ /* Test parameters */ - const size_t k = 3; - const double bpoint_data[] = { 0.0, 2.0, 4.0 }; + const size_t k = 1; + const double bpoint_data[] = { 0.0, 0.2, 0.5, 0.75, 1.0 }; + const size_t nbreak = sizeof(bpoint_data)/sizeof(bpoint_data[0]); + + /* Expected results */ + const double abscissae_data[] = { 0.1, 0.35, 0.625, 0.875 }; + const size_t nabscissae = sizeof(abscissae_data)/sizeof(abscissae_data[0]); + + gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak); + gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak); + gsl_bspline_knots((const gsl_vector *) &bpoints, w); + + gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w), + "b-spline k=%d number of abscissae", k); + for (i = 0; i < nabscissae; ++i) + { + gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON, + "b-spline k=%d Greville abscissa #%d at x = %f", k, i, abscissae_data[i]); + } + + gsl_bspline_free(w); + } + + /* Check Greville abscissae functionality on a non-uniform k=2 */ + { + size_t i; /* looping */ + + /* Test parameters */ + const size_t k = 2; + const double bpoint_data[] = { 0.0, 0.2, 0.5, 0.75, 1.0 }; const size_t nbreak = sizeof(bpoint_data)/sizeof(bpoint_data[0]); /* Expected results */ - const double abscissae_data[] = { 0.0, 2.0/3.0, 2.0, 10.0/3.0, 4.0 }; + const double abscissae_data[] = { 0.0, 0.2, 0.5, 0.75, 1.0 }; const size_t nabscissae = sizeof(abscissae_data)/sizeof(abscissae_data[0]); gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak); gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak); gsl_bspline_knots((const gsl_vector *) &bpoints, w); - gsl_test_int(gsl_bspline_greville_nabscissae(w), nabscissae, - "b-spline k=%d number of Greville abscissae is %d", k, nabscissae); + gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w), + "b-spline k=%d number of abscissae", k); for (i = 0; i < nabscissae; ++i) { gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON, @@ -305,16 +333,16 @@ main(int argc, char **argv) const size_t nbreak = sizeof(bpoint_data)/sizeof(bpoint_data[0]); /* Expected results */ - const double abscissae_data[] = { 0.0, 1.0/15.0, 7.0/30.0, - 29.0/60.0, 3.0/4.0, 11.0/12.0, 1.0 }; + const double abscissae_data[] = { 0.0, 1.0/10.0, 7.0/20.0, + 5.0/ 8.0, 7.0/ 8.0, 1.0 }; const size_t nabscissae = sizeof(abscissae_data)/sizeof(abscissae_data[0]); gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak); gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak); gsl_bspline_knots((const gsl_vector *) &bpoints, w); - gsl_test_int(gsl_bspline_greville_nabscissae(w), nabscissae, - "b-spline k=%d number of Greville abscissae is %d", k, nabscissae); + gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w), + "b-spline k=%d number of abscissae", k); for (i = 0; i < nabscissae; ++i) { gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON, @@ -334,16 +362,16 @@ main(int argc, char **argv) const size_t nbreak = sizeof(bpoint_data)/sizeof(bpoint_data[0]); /* Expected results */ - const double abscissae_data[] = { 0.0, 1.0/20.0, 7.0/40.0, 29.0/80.0, - 49.0/80.0, 13.0/16.0, 15.0/16.0, 1.0 }; + const double abscissae_data[] = { 0.0, 1.0/15.0, 7.0/30.0, 29.0/60.0, + 3.0/ 4.0, 11.0/12.0, 1.0 }; const size_t nabscissae = sizeof(abscissae_data)/sizeof(abscissae_data[0]); gsl_vector_const_view bpoints = gsl_vector_const_view_array(bpoint_data, nbreak); gsl_bspline_workspace *w = gsl_bspline_alloc(k, nbreak); gsl_bspline_knots((const gsl_vector *) &bpoints, w); - gsl_test_int(gsl_bspline_greville_nabscissae(w), nabscissae, - "b-spline k=%d number of Greville abscissae is %d", k, nabscissae); + gsl_test_int(nabscissae, gsl_bspline_ncoeffs(w), + "b-spline k=%d number of abscissae", k); for (i = 0; i < nabscissae; ++i) { gsl_test_abs(gsl_bspline_greville_abscissa(i, w), abscissae_data[i], 2*k*GSL_DBL_EPSILON, diff --git a/doc/bspline.texi b/doc/bspline.texi index 8df82be..eb64c47 100644 --- a/doc/bspline.texi +++ b/doc/bspline.texi @@ -205,17 +205,19 @@ their derivatives to be computed without unnecessary terms. @section Greville abscissae @cindex basis splines, Greville abscissae -The Greville abscissae are defined to be the mean location of @math{k} -consecutive knots in the knot vector for each basis spline of order @math{k}. -They are often used in B-spline collocation applications. +The Greville abscissae are defined to be the mean location of @math{k-1} +consecutive knots in the knot vector for each basis spline function of order +@math{k}. Note that the first and last knots in the knot vector are excluded +when applying this definition; consequently there are +@code{gsl_bspline_ncoeffs} Greville abscissa. They are often used in B-spline +collocation applications and may also be called Marsden-Schoenberg points. -@deftypefun int gsl_bspline_greville_nabscissae (gsl_bspline_workspace * @var{w}) -Returns the number of Greville abscissae for the given spline basis. -@end deftypefun +The above definition is undefined for @math{k=1}. The implementation chooses +to return interval midpoints in the degenerate @math{k=1} case. @deftypefun double gsl_bspline_greville_abscissa (size_t @var{i}, gsl_bspline_workspace *@var{w}); Returns the location of the @math{i}-th Greville abscissa for the given spline -basis. +basis. Here, @math{i = 0}, ..., @code{gsl_bspline_ncoeffs(w)}. @end deftypefun @node Example programs for B-splines @@ -255,6 +257,14 @@ C. de Boor, @cite{A Practical Guide to Splines} (1978), Springer-Verlag, ISBN 0-387-90356-9. @end itemize +Further information of Greville abscissae and B-spline collocation +can be found in the following paper, + +@itemize @asis +Richard W. Johnson, Higher order B-spline collocation at the Greville +abscissae. @cite{Applied Numerical Mathematics}. vol.@: 52, 2005, 63--75. +@end itemize + @noindent A large collection of B-spline routines is available in the @sc{pppack} library available at @uref{http://www.netlib.org/pppack}, -- 1.5.4.3