[PATCH] Add Greville abscissae functionality to B-splines

[PATCH] Add Greville abscissae functionality to B-splines

[Rhys Ulerich]

 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

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Brian Gough]

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

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Rhys Ulerich]

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
>
>

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Patrick Alken]

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
> >
> >

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Rhys Ulerich]

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
>> >
>> >
>

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Patrick Alken]

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
> >> >
> >> >
> >

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Rhys Ulerich]

[-- 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

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Patrick Alken]

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
> >> >> >
> >> >> >
> >> >
> >

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Rhys Ulerich]

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

Re: [PATCH] Add Greville abscissae functionality to B-splines

[Patrick Alken]

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
> >> >> >
> >> >> >
> >> >
> >