* [PATCH] Add Greville abscissae functionality to B-splines @ 2009-06-14 17:02 Rhys Ulerich 2009-06-15 10:30 ` Brian Gough 0 siblings, 1 reply; 10+ messages in thread From: Rhys Ulerich @ 2009-06-14 17:02 UTC (permalink / raw) To: gsl-discuss This change adds computing Greville abscissae to the GSL B-spline routines. Updates to unit tests and documentation are included. The routines are written so that if the b-spline classes have lower continuity basis added later (i.e. by adding multiple knots per interior breakpoint), these should continue to do the right thing. --- bspline/Makefile.am | 2 + bspline/bspline.c | 24 +++++++++++++++ bspline/gsl_bspline.h | 2 + bspline/test.c | 80 +++++++++++++++++++++++++++++++++++++++++++++++++ doc/bspline.texi | 30 +++++++++++++++--- 5 files changed, 131 insertions(+), 7 deletions(-) diff --git a/bspline/Makefile.am b/bspline/Makefile.am index c96a3f5..b4179eb 100644 --- a/bspline/Makefile.am +++ b/bspline/Makefile.am @@ -12,6 +12,6 @@ check_PROGRAMS = test TESTS = $(check_PROGRAMS) -test_LDADD = libgslbspline.la ../linalg/libgsllinalg.la ../permutation/libgslpermutation.la ../blas/libgslblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../complex/libgslcomplex.la ../cblas/libgslcblas.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +test_LDADD = libgslbspline.la ../linalg/libgsllinalg.la ../permutation/libgslpermutation.la ../blas/libgslblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../complex/libgslcomplex.la ../cblas/libgslcblas.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la ../statistics/libgslstatistics.la test_SOURCES = test.c diff --git a/bspline/bspline.c b/bspline/bspline.c index f72ff27..d1e4470 100644 --- a/bspline/bspline.c +++ b/bspline/bspline.c @@ -21,6 +21,7 @@ #include <config.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_bspline.h> +#include <gsl/gsl_statistics.h> /* * This module contains routines related to calculating B-splines. @@ -207,6 +208,29 @@ 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))) + { + 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; + + return gsl_stats_mean(data, stride, w->k); +} + /* gsl_bspline_free() Free a gsl_bspline_workspace. diff --git a/bspline/gsl_bspline.h b/bspline/gsl_bspline.h index 04262dd..f30fefa 100644 --- a/bspline/gsl_bspline.h +++ b/bspline/gsl_bspline.h @@ -68,6 +68,8 @@ 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 gsl_bspline_knots(const gsl_vector *breakpts, gsl_bspline_workspace *w); diff --git a/bspline/test.c b/bspline/test.c index 2e233e6..0845b0b 100644 --- a/bspline/test.c +++ b/bspline/test.c @@ -267,5 +267,85 @@ main(int argc, char **argv) gsl_vector_free(breakpts); } + /* Check Greville abscissae functionality on a uniform k=3 */ + { + size_t i; /* looping */ + + /* Test parameters */ + const size_t k = 3; + const double bpoint_data[] = { 0.0, 2.0, 4.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 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); + 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]); + } + } + + /* Check Greville abscissae functionality on non-uniform k=3 */ + { + size_t i; /* looping */ + + /* Test parameters */ + const size_t k = 3; + 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, 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); + 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]); + } + } + + /* Check Greville abscissae functionality on non-uniform k=4 */ + { + size_t i; /* looping */ + + /* Test parameters */ + const size_t k = 4; + 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, 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 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); + 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]); + } + } + exit(gsl_test_summary()); } diff --git a/doc/bspline.texi b/doc/bspline.texi index d8d14a8..8df82be 100644 --- a/doc/bspline.texi +++ b/doc/bspline.texi @@ -16,6 +16,7 @@ bspline functions and related declarations. * Constructing the knots vector:: * Evaluation of B-spline basis functions:: * Evaluation of B-spline basis function derivatives:: +* Obtaining Greville abscissae for B-spline basis functions:: * Example programs for B-splines:: * References and Further Reading:: @end menu @@ -51,7 +52,7 @@ $$ @example B_(i,1)(x) = (1, t_i <= x < t_(i+1) (0, else -B_(i,k)(x) = [(x - t_i)/(t_(i+k-1) - t_i)] B_(i,k-1)(x) +B_(i,k)(x) = [(x - t_i)/(t_(i+k-1) - t_i)] B_(i,k-1)(x) + [(t_(i+k) - x)/(t_(i+k) - t_(i+1))] B_(i+1,k-1)(x) @end example @@ -155,10 +156,10 @@ that the @math{i}-th element is @c{$B_{(istart+i)}(x)$} The last element of @var{Bk} is @c{$B_{iend}(x)$} @math{B_(iend)(x)}. The vector @var{Bk} must be of length @math{k}. By returning only the nonzero basis functions, -this function +this function allows quantities involving linear combinations of the @math{B_i(x)} to be computed without unnecessary terms -(such linear combinations occur, for example, +(such linear combinations occur, for example, when evaluating an interpolated function). @end deftypefun @@ -175,8 +176,8 @@ This function returns the number of B-spline coefficients given by 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 the matrix @var{dB}. The @math{(i,j)}-th element of @var{dB} -is @math{d^jB_i(x)/dx^j}. The matrix @var{dB} must be -of size @math{n = nbreak + k - 2} by @math{nderiv + 1}. +is @math{d^jB_i(x)/dx^j}. The matrix @var{dB} must be +of size @math{n = nbreak + k - 2} by @math{nderiv + 1}. The value @math{n} may also be obtained by calling @code{gsl_bspline_ncoeffs}. Note that function evaluations are included as the zeroth order derivatives in @var{dB}. @@ -188,7 +189,7 @@ recurrence relation. @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}, gsl_bspline_deriv_workspace * @var{dw}) 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 the matrix @var{dB}. The +the position @var{x} and stores them in the matrix @var{dB}. The @math{(i,j)}-th element of @var{dB} is @c{$d^jB_{(istart+i)}(x)/dx^j$} @math{d^j/dx^j B_(istart+i)(x)}. The last row of @var{dB} contains @c{$d^jB_{iend}(x)/dx^j$} @@ -200,6 +201,23 @@ quantities involving linear combinations of the @math{B_i(x)} and their derivatives to be computed without unnecessary terms. @end deftypefun +@node Obtaining Greville abscissae for B-spline basis functions +@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. + +@deftypefun int gsl_bspline_greville_nabscissae (gsl_bspline_workspace * @var{w}) +Returns the number of Greville abscissae for the given spline basis. +@end deftypefun + +@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. +@end deftypefun + @node Example programs for B-splines @section Examples @cindex basis splines, examples ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-14 17:02 [PATCH] Add Greville abscissae functionality to B-splines Rhys Ulerich @ 2009-06-15 10:30 ` Brian Gough 2009-06-15 14:44 ` Rhys Ulerich 0 siblings, 1 reply; 10+ messages in thread From: Brian Gough @ 2009-06-15 10:30 UTC (permalink / raw) To: Rhys Ulerich; +Cc: gsl-discuss At Sun, 14 Jun 2009 12:02:10 -0500, Rhys Ulerich wrote: > This change adds computing Greville abscissae to the GSL B-spline > routines. Updates to unit tests and documentation are included. > The routines are written so that if the b-spline classes have lower > continuity basis added later (i.e. by adding multiple knots per > interior breakpoint), these should continue to do the right thing. Cool. I'll let Patrick take care of these. -- Brian Gough ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-15 10:30 ` Brian Gough @ 2009-06-15 14:44 ` Rhys Ulerich 2009-06-15 17:38 ` Patrick Alken 0 siblings, 1 reply; 10+ messages in thread From: Rhys Ulerich @ 2009-06-15 14:44 UTC (permalink / raw) To: Brian Gough; +Cc: gsl-discuss Thanks guys. FYI, I woke up this morning and realized I forgot a gsl_bspline_free call at the bottom of each of the unit test blocks in test.c. - Rhys On Mon, Jun 15, 2009 at 5:26 AM, Brian Gough<bjg@network-theory.co.uk> wrote: > At Sun, 14 Jun 2009 12:02:10 -0500, > Rhys Ulerich wrote: >> This change adds computing Greville abscissae to the GSL B-spline >> routines. Updates to unit tests and documentation are included. >> The routines are written so that if the b-spline classes have lower >> continuity basis added later (i.e. by adding multiple knots per >> interior breakpoint), these should continue to do the right thing. > > Cool. I'll let Patrick take care of these. > > -- > Brian Gough > > ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-15 14:44 ` Rhys Ulerich @ 2009-06-15 17:38 ` Patrick Alken 2009-06-15 18:14 ` Rhys Ulerich 0 siblings, 1 reply; 10+ messages in thread From: Patrick Alken @ 2009-06-15 17:38 UTC (permalink / raw) To: Rhys Ulerich; +Cc: gsl-discuss Nice work, I have added your patch to the repository and inserted the appropriate gsl_bspline_free's in test.c. Patrick On Mon, Jun 15, 2009 at 09:43:46AM -0500, Rhys Ulerich wrote: > Thanks guys. > > FYI, I woke up this morning and realized I forgot a gsl_bspline_free > call at the bottom of each of the unit test blocks in test.c. > > - Rhys > > On Mon, Jun 15, 2009 at 5:26 AM, Brian Gough<bjg@network-theory.co.uk> wrote: > > At Sun, 14 Jun 2009 12:02:10 -0500, > > Rhys Ulerich wrote: > >> Â This change adds computing Greville abscissae to the GSL B-spline > >> Â routines. Â Updates to unit tests and documentation are included. > >> Â The routines are written so that if the b-spline classes have lower > >> Â continuity basis added later (i.e. by adding multiple knots per > >> Â interior breakpoint), these should continue to do the right thing. > > > > Cool. I'll let Patrick take care of these. > > > > -- > > Brian Gough > > > > ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-15 17:38 ` Patrick Alken @ 2009-06-15 18:14 ` Rhys Ulerich 2009-06-15 18:31 ` Patrick Alken 0 siblings, 1 reply; 10+ messages in thread From: Rhys Ulerich @ 2009-06-15 18:14 UTC (permalink / raw) To: Patrick Alken; +Cc: gsl-discuss Hi Patrick, I noticed in the commit there were changes made to install.sh (http://git.savannah.gnu.org/gitweb/?p=gsl.git;a=commit;h=ccccc3cb7630ea43bdf365bcbaec71e31f8aea91). Were those deliberate? - Rhys On Mon, Jun 15, 2009 at 12:39 PM, Patrick Alken<patrick.alken@colorado.edu> wrote: > Nice work, I have added your patch to the repository and > inserted the appropriate gsl_bspline_free's in test.c. > > Patrick > > On Mon, Jun 15, 2009 at 09:43:46AM -0500, Rhys Ulerich wrote: >> Thanks guys. >> >> FYI, I woke up this morning and realized I forgot a gsl_bspline_free >> call at the bottom of each of the unit test blocks in test.c. >> >> - Rhys >> >> On Mon, Jun 15, 2009 at 5:26 AM, Brian Gough<bjg@network-theory.co.uk> wrote: >> > At Sun, 14 Jun 2009 12:02:10 -0500, >> > Rhys Ulerich wrote: >> >> This change adds computing Greville abscissae to the GSL B-spline >> >> routines. Updates to unit tests and documentation are included. >> >> The routines are written so that if the b-spline classes have lower >> >> continuity basis added later (i.e. by adding multiple knots per >> >> interior breakpoint), these should continue to do the right thing. >> > >> > Cool. I'll let Patrick take care of these. >> > >> > -- >> > Brian Gough >> > >> > > ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-15 18:14 ` Rhys Ulerich @ 2009-06-15 18:31 ` Patrick Alken 2009-06-17 21:54 ` Rhys Ulerich 0 siblings, 1 reply; 10+ messages in thread From: Patrick Alken @ 2009-06-15 18:31 UTC (permalink / raw) To: Rhys Ulerich; +Cc: gsl-discuss Oops, I must have ran some autoconf tools and changed that file. I think its fixed now. Patrick On Mon, Jun 15, 2009 at 01:13:32PM -0500, Rhys Ulerich wrote: > Hi Patrick, > > I noticed in the commit there were changes made to install.sh > (http://git.savannah.gnu.org/gitweb/?p=gsl.git;a=commit;h=ccccc3cb7630ea43bdf365bcbaec71e31f8aea91). > Were those deliberate? > > - Rhys > > On Mon, Jun 15, 2009 at 12:39 PM, Patrick > Alken<patrick.alken@colorado.edu> wrote: > > Nice work, I have added your patch to the repository and > > inserted the appropriate gsl_bspline_free's in test.c. > > > > Patrick > > > > On Mon, Jun 15, 2009 at 09:43:46AM -0500, Rhys Ulerich wrote: > >> Thanks guys. > >> > >> FYI, I woke up this morning and realized I forgot a gsl_bspline_free > >> call at the bottom of each of the unit test blocks in test.c. > >> > >> - Rhys > >> > >> On Mon, Jun 15, 2009 at 5:26 AM, Brian Gough<bjg@network-theory.co.uk> wrote: > >> > At Sun, 14 Jun 2009 12:02:10 -0500, > >> > Rhys Ulerich wrote: > >> >> Â This change adds computing Greville abscissae to the GSL B-spline > >> >> Â routines. Â Updates to unit tests and documentation are included. > >> >> Â The routines are written so that if the b-spline classes have lower > >> >> Â continuity basis added later (i.e. by adding multiple knots per > >> >> Â interior breakpoint), these should continue to do the right thing. > >> > > >> > Cool. I'll let Patrick take care of these. > >> > > >> > -- > >> > Brian Gough > >> > > >> > > > ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-15 18:31 ` Patrick Alken @ 2009-06-17 21:54 ` Rhys Ulerich 2009-06-17 22:11 ` Patrick Alken 2009-07-05 1:10 ` Patrick Alken 0 siblings, 2 replies; 10+ messages in thread From: Rhys Ulerich @ 2009-06-17 21:54 UTC (permalink / raw) To: Patrick Alken; +Cc: gsl-discuss [-- Attachment #1: Type: text/plain, Size: 1955 bytes --] Hi Patrick (and all), Attached is a patch to apply atop my last Greville abscissae contribution. I had badly misunderstood some of the details when I implemented it originally. This new patch corrects the implementation, units tests, and documentation. - Rhys On Mon, Jun 15, 2009 at 1:31 PM, Patrick Alken<patrick.alken@colorado.edu> wrote: > Oops, I must have ran some autoconf tools and changed that file. > I think its fixed now. > > Patrick > > On Mon, Jun 15, 2009 at 01:13:32PM -0500, Rhys Ulerich wrote: >> Hi Patrick, >> >> I noticed in the commit there were changes made to install.sh >> (http://git.savannah.gnu.org/gitweb/?p=gsl.git;a=commit;h=ccccc3cb7630ea43bdf365bcbaec71e31f8aea91). >> Were those deliberate? >> >> - Rhys >> >> On Mon, Jun 15, 2009 at 12:39 PM, Patrick >> Alken<patrick.alken@colorado.edu> wrote: >> > Nice work, I have added your patch to the repository and >> > inserted the appropriate gsl_bspline_free's in test.c. >> > >> > Patrick >> > >> > On Mon, Jun 15, 2009 at 09:43:46AM -0500, Rhys Ulerich wrote: >> >> Thanks guys. >> >> >> >> FYI, I woke up this morning and realized I forgot a gsl_bspline_free >> >> call at the bottom of each of the unit test blocks in test.c. >> >> >> >> - Rhys >> >> >> >> On Mon, Jun 15, 2009 at 5:26 AM, Brian Gough<bjg@network-theory.co.uk> wrote: >> >> > At Sun, 14 Jun 2009 12:02:10 -0500, >> >> > Rhys Ulerich wrote: >> >> >> This change adds computing Greville abscissae to the GSL B-spline >> >> >> routines. Updates to unit tests and documentation are included. >> >> >> The routines are written so that if the b-spline classes have lower >> >> >> continuity basis added later (i.e. by adding multiple knots per >> >> >> interior breakpoint), these should continue to do the right thing. >> >> > >> >> > Cool. I'll let Patrick take care of these. >> >> > >> >> > -- >> >> > Brian Gough >> >> > >> >> > >> > > [-- Attachment #2: 0001-Fixed-Greville-abscissae-implementation.patch --] [-- Type: text/x-diff, Size: 9601 bytes --] From 7ac6ed29cc302f00228732d4000bd324a2a95241 Mon Sep 17 00:00:00 2001 From: Rhys Ulerich <rhys.ulerich@gmail.com> 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-17 21:54 ` Rhys Ulerich @ 2009-06-17 22:11 ` Patrick Alken 2009-06-17 22:52 ` Rhys Ulerich 2009-07-05 1:10 ` Patrick Alken 1 sibling, 1 reply; 10+ messages in thread From: Patrick Alken @ 2009-06-17 22:11 UTC (permalink / raw) To: gsl-discuss I'll add that I've never used Greville abscissae before, but did do some reading about it when I looked at your patch. I also found a matlab code on it which I used to verify some of your test cases, which did look correct. Maybe it would be worthwhile to back out this patch from the git and do more testing on it first? I'll be out of town for 2 weeks and won't be able to look at your new patch until I get back. Patrick On Wed, Jun 17, 2009 at 04:54:21PM -0500, Rhys Ulerich wrote: > Hi Patrick (and all), > > Attached is a patch to apply atop my last Greville abscissae > contribution. I had badly misunderstood some of the details when I > implemented it originally. This new patch corrects the > implementation, units tests, and documentation. > > - Rhys > > On Mon, Jun 15, 2009 at 1:31 PM, Patrick > Alken<patrick.alken@colorado.edu> wrote: > > Oops, I must have ran some autoconf tools and changed that file. > > I think its fixed now. > > > > Patrick > > > > On Mon, Jun 15, 2009 at 01:13:32PM -0500, Rhys Ulerich wrote: > >> Hi Patrick, > >> > >> I noticed in the commit there were changes made to install.sh > >> (http://git.savannah.gnu.org/gitweb/?p=gsl.git;a=commit;h=ccccc3cb7630ea43bdf365bcbaec71e31f8aea91). > >> Â Were those deliberate? > >> > >> - Rhys > >> > >> On Mon, Jun 15, 2009 at 12:39 PM, Patrick > >> Alken<patrick.alken@colorado.edu> wrote: > >> > Nice work, I have added your patch to the repository and > >> > inserted the appropriate gsl_bspline_free's in test.c. > >> > > >> > Patrick > >> > > >> > On Mon, Jun 15, 2009 at 09:43:46AM -0500, Rhys Ulerich wrote: > >> >> Thanks guys. > >> >> > >> >> FYI, I woke up this morning and realized I forgot a gsl_bspline_free > >> >> call at the bottom of each of the unit test blocks in test.c. > >> >> > >> >> - Rhys > >> >> > >> >> On Mon, Jun 15, 2009 at 5:26 AM, Brian Gough<bjg@network-theory.co.uk> wrote: > >> >> > At Sun, 14 Jun 2009 12:02:10 -0500, > >> >> > Rhys Ulerich wrote: > >> >> >> Â This change adds computing Greville abscissae to the GSL B-spline > >> >> >> Â routines. Â Updates to unit tests and documentation are included. > >> >> >> Â The routines are written so that if the b-spline classes have lower > >> >> >> Â continuity basis added later (i.e. by adding multiple knots per > >> >> >> Â interior breakpoint), these should continue to do the right thing. > >> >> > > >> >> > Cool. I'll let Patrick take care of these. > >> >> > > >> >> > -- > >> >> > Brian Gough > >> >> > > >> >> > > >> > > > ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-17 22:11 ` Patrick Alken @ 2009-06-17 22:52 ` Rhys Ulerich 0 siblings, 0 replies; 10+ messages in thread From: Rhys Ulerich @ 2009-06-17 22:52 UTC (permalink / raw) To: Patrick Alken; +Cc: gsl-discuss Hi Patrick, I'm fine if you'd like to back out the previous patch until you can review the two applied together. The trap I fell into is that, confusingly, the indexing of the abscissa as it usually appears excludes the very first and very last knots in GSL's w->knots vector. Apparently this is a common convention I didn't know about until I ran into a friend in the hall today. Also, the abscissa location should be averaged over degree (w->km1) knots and not order (w->k) knots. > I also found a matlab code on it which I used to verify some of your > test cases, which did look correct. Odds are good that you can get the new unit tests will also look okay against that same Matlab code if you (a) throw out the first and last knot, and (b) use the degree and not the order. The new patch's documentation includes a reference to a 2005 paper by Johnson that gives both a definition and an example set of knots/abscissae for a C3 quartic (Figure 2) if you'd like another comparison point. Have a good trip, Rhys ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: [PATCH] Add Greville abscissae functionality to B-splines 2009-06-17 21:54 ` Rhys Ulerich 2009-06-17 22:11 ` Patrick Alken @ 2009-07-05 1:10 ` Patrick Alken 1 sibling, 0 replies; 10+ messages in thread From: Patrick Alken @ 2009-07-05 1:10 UTC (permalink / raw) To: Rhys Ulerich; +Cc: gsl-discuss I have committed this patch. Patrick On Wed, Jun 17, 2009 at 04:54:21PM -0500, Rhys Ulerich wrote: > Hi Patrick (and all), > > Attached is a patch to apply atop my last Greville abscissae > contribution. I had badly misunderstood some of the details when I > implemented it originally. This new patch corrects the > implementation, units tests, and documentation. > > - Rhys > > On Mon, Jun 15, 2009 at 1:31 PM, Patrick > Alken<patrick.alken@colorado.edu> wrote: > > Oops, I must have ran some autoconf tools and changed that file. > > I think its fixed now. > > > > Patrick > > > > On Mon, Jun 15, 2009 at 01:13:32PM -0500, Rhys Ulerich wrote: > >> Hi Patrick, > >> > >> I noticed in the commit there were changes made to install.sh > >> (http://git.savannah.gnu.org/gitweb/?p=gsl.git;a=commit;h=ccccc3cb7630ea43bdf365bcbaec71e31f8aea91). > >> Â Were those deliberate? > >> > >> - Rhys > >> > >> On Mon, Jun 15, 2009 at 12:39 PM, Patrick > >> Alken<patrick.alken@colorado.edu> wrote: > >> > Nice work, I have added your patch to the repository and > >> > inserted the appropriate gsl_bspline_free's in test.c. > >> > > >> > Patrick > >> > > >> > On Mon, Jun 15, 2009 at 09:43:46AM -0500, Rhys Ulerich wrote: > >> >> Thanks guys. > >> >> > >> >> FYI, I woke up this morning and realized I forgot a gsl_bspline_free > >> >> call at the bottom of each of the unit test blocks in test.c. > >> >> > >> >> - Rhys > >> >> > >> >> On Mon, Jun 15, 2009 at 5:26 AM, Brian Gough<bjg@network-theory.co.uk> wrote: > >> >> > At Sun, 14 Jun 2009 12:02:10 -0500, > >> >> > Rhys Ulerich wrote: > >> >> >> Â This change adds computing Greville abscissae to the GSL B-spline > >> >> >> Â routines. Â Updates to unit tests and documentation are included. > >> >> >> Â The routines are written so that if the b-spline classes have lower > >> >> >> Â continuity basis added later (i.e. by adding multiple knots per > >> >> >> Â interior breakpoint), these should continue to do the right thing. > >> >> > > >> >> > Cool. I'll let Patrick take care of these. > >> >> > > >> >> > -- > >> >> > Brian Gough > >> >> > > >> >> > > >> > > > ^ permalink raw reply [flat|nested] 10+ messages in thread
end of thread, other threads:[~2009-07-05 1:10 UTC | newest] Thread overview: 10+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2009-06-14 17:02 [PATCH] Add Greville abscissae functionality to B-splines Rhys Ulerich 2009-06-15 10:30 ` Brian Gough 2009-06-15 14:44 ` Rhys Ulerich 2009-06-15 17:38 ` Patrick Alken 2009-06-15 18:14 ` Rhys Ulerich 2009-06-15 18:31 ` Patrick Alken 2009-06-17 21:54 ` Rhys Ulerich 2009-06-17 22:11 ` Patrick Alken 2009-06-17 22:52 ` Rhys Ulerich 2009-07-05 1:10 ` Patrick Alken
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).