From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 11095 invoked by alias); 14 Jun 2009 17:02:27 -0000 Received: (qmail 9643 invoked by uid 22791); 14 Jun 2009 17:02:25 -0000 X-SWARE-Spam-Status: No, hits=0.4 required=5.0 tests=BAYES_20,J_CHICKENPOX_43,SPF_PASS,URI_NOVOWEL X-Spam-Check-By: sourceware.org Received: from qw-out-1920.google.com (HELO qw-out-1920.google.com) (74.125.92.147) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Sun, 14 Jun 2009 17:02:16 +0000 Received: by qw-out-1920.google.com with SMTP id 9so1748154qwj.24 for ; Sun, 14 Jun 2009 10:02:13 -0700 (PDT) Received: by 10.224.11.132 with SMTP id t4mr6401325qat.377.1244998933712; Sun, 14 Jun 2009 10:02:13 -0700 (PDT) Received: from localhost (cpe-70-112-187-78.austin.res.rr.com [70.112.187.78]) by mx.google.com with ESMTPS id 26sm6593613qwa.8.2009.06.14.10.02.11 (version=TLSv1/SSLv3 cipher=RC4-MD5); Sun, 14 Jun 2009 10:02:12 -0700 (PDT) Received: from [127.0.0.1] (localhost [127.0.0.1]) by localhost (Postfix) with ESMTP id 6E295C606D; Sun, 14 Jun 2009 12:02:10 -0500 (CDT) From: Rhys Ulerich Subject: [PATCH] Add Greville abscissae functionality to B-splines To: gsl-discuss@sourceware.org Date: Sun, 14 Jun 2009 17:02:00 -0000 Message-ID: <20090614165708.22998.15238.stgit@localhost> User-Agent: StGIT/0.14.2 MIME-Version: 1.0 Content-Type: text/plain; charset="utf-8" Content-Transfer-Encoding: 7bit Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2009-q2/txt/msg00016.txt.bz2 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 #include #include +#include /* * 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