This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
Re: [PATCH] Add Greville abscissae functionality to B-splines
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
>> >> >
>> >> >
>> >
>
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