[PATCH] Add Greville abscissae functionality to B-splines

[Rhys Ulerich <rhys dot ulerich at gmail dot com>]

 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