Re: RE : [Bug-gsl] GSL - bspline

[Patrick Alken]

Thank you for reporting this. I have fixed it in the latest CVS.

On Tue, Nov 27, 2007 at 10:30:24PM +0100, Christophe Cloquet wrote:
>    Dear Patrick,
> 
>    Thanks for your e-mail. I apologize for this extremely late answer. I
>    enclosed a sample code in this mail. I don't actually know if this is a
>    bug, or a bad implementation from my side.
> 
>    The code defines a cubic spline basis bw3, then a linear spline basis bw2
>    based on bw3. The problem arise when evaluating the splines in bw2 at the
>    last breakpoint. It is only a minor problem, and I could fix it for my
>    application. The final aim is to compute the second derivative of the bw3
>    basis, based on De Boor, Practical Guide to splines, 1978, p. 138.
> 
>    If something is unclear, please ask me. I'll try to answer as soon as
>    possible.
> 
> 
> 
>    Best regards,
> 
> 
> 
>    Christophe
> 
>    Here is the code :
>    
> ====================================================================================
> 
>    #include <stdlib.h>
>    #include <stdio.h>
>    #include <gsl/gsl_bspline.h>
> 
>    #define bugfixed // comment this line to see the effect of the bug on the
>    last
>                     // line of the table that is displayed on the screen
> 
>    int main(void) {
> 
>         int j, n;
>         double s;
> 
>        //bw3 is a cubic spline basis
>         size_t                 nbreaks  =10;
>         size_t                 nsplines =nbreaks+2;
>         size_t                 k        =4;
>         gsl_bspline_workspace *bw3;
>         gsl_vector            *breaks = gsl_vector_alloc(nbreaks);
>         for (j=0; j< nbreaks; j++) gsl_vector_set(breaks, j, (double)j);
> 
>         bw3            = gsl_bspline_alloc(k, nbreaks);  // the workspace for
>    cubical bspline (k=4)
>         gsl_vector *B3 =
>    gsl_vector_alloc(nsplines);
> 
>         gsl_bspline_knots(breaks, bw3);                  // knots
>    allocation
> 
>         //bw2 is a linear spline basis deriveed from bw3
>         double *t           = calloc(nsplines+4, sizeof(double));
>         gsl_vector *breaks2 = gsl_vector_alloc(nbreaks+4);
> 
>         for (j=0; j < nbreaks+4; j++)
>    gsl_vector_set(breaks2,j,gsl_vector_get(bw3->knots,j+2)); // sets the
>    breakpoints of the (k-2)order bspline, in order to have the **same knots**
> 
>         gsl_bspline_workspace *bw2 = gsl_bspline_alloc(k-2, nbreaks+4); //
>    quadric spline, with nbreaks+4 breakpoints,
>                                                                         //
>    hence nbreaks+4+2*2 knots
>         gsl_bspline_knots(breaks2, bw2);                                //
>    knots allocation
> 
>         size_t ncoeffs_2           = gsl_bspline_ncoeffs(bw2);
>         gsl_vector *B2             = gsl_vector_alloc(ncoeffs_2);
> 
>         printf("\n\n The following table shows the result of the evaluation
>    of the bsplines \n for several abscissa. \n Each row is the result for 1
>    breakpoint. \n Each column is the result for 1 B spline.\n\n");
> 
>         for(j=0; j<nbreaks; ++j) {
>             s = gsl_vector_get(breaks, j);
> 
>             gsl_bspline_eval(s, B2, bw2); // value of all the splines of
>    order 2 @ x=x[j].
>             gsl_bspline_eval(s, B3, bw3);
> 
>             //
>    ******************************************************************
>             // * for the evaluation of the linear spline at the last
>    breakpoint,*
>             // * I need to use the following
>    trick.                             *
>             // * Otherwise the evaluation for the last
>    breakpoints              *
>             // * gives 0 everywhere (last line of the table that is
>    displayed   *
>             // * on the
>    screen).                                                *
>             // * This problem only appears with B2 (the linear spline
>    basis)    *
>             // * and not with B3 (the cubic spline
>    basis)                       *
>             // * Is this a bug, or a bad implementation from my side
>    ?          *
>             //
>    ******************************************************************
> 
>             #ifdef bugfixed
>             if(j==nbreaks-1) {
>                          gsl_bspline_eval(s-1e-10, B2, bw2); // value of all
>    the splines of order 2 @ x=x[j].
>                          gsl_bspline_eval(s-1e-10, B3, bw3);
>                          }
>             #endif
> 
>             printf("bkpnt nr %d : ", j);
>             for (n=0; n<nsplines; n++) printf("%2.0f ",
>    gsl_vector_get(B2,n));
>             printf("\n");
>        }
>    }
>    ======================================================================
> 
>    I compile the code with the following lines :
> 
>    lcc -c -I"C:\progam files\lcc\include" -I"C:\Program
>    Files\GnuWin32\include" -g2  test_bspline.c
> 
>    lcclnk -L"C:\Program Files\GnuWin32\lib" -subsystem console -o
>    test_bspline.exe test_bspline.obj "C:\Program Files\GnuWin32\lib\libgsl.a"
>    "C:\Program Files\GnuWin32\lib\libgslcblas.a" "C:\Program
>    Files\GnuWin32\lib\bspline\bspline.obj"