[PATCH] Add Greville abscissae functionality to B-splines
Rhys Ulerich
rhys.ulerich@gmail.com
Sun Jun 14 17:02:00 GMT 2009
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
More information about the Gsl-discuss
mailing list