Dear all, I found the following problem which I think is a bug on the construction of B-spline basis from non uniform breakpts
The constructor: gsl_bspline_knot(gsl_vector * breakpts, gsl_bspline_workspace *bw) does not *always* check correctly for increasing values in the supplied break points vector breakpts. Specifically (sometimes) it runs with no warning, and gives values for B-spline basis that are inconsistent (e.g. negative values for the basis functions B_i(x)). I attach: a) A plot of the calculated B-spline basis B_i(x) that were created for a knot vector with non increasing values. b) A program that demonstrates this. Specifically, I have two choices of breakpoints, one that the program compiles and runs and produces inconsistent results, and another that the program exits on run time, complaining for non increasing order of knot vector. GSL version 1.16 Operating system: Ubuntu 12.04 LTS Compiler version: g++ 4.7.3 File compiled with command: g++ -std=c++11 test_bspline_order.cpp -o test_bspline_order.xxx -lgsl -lgslcblas Thank you all for the great help you provide. All the best, Foivos
#include <iostream> #include <cstdlib> #include <cstdio> #include <vector> #include <gsl/gsl_math.h> #include <gsl/gsl_bspline.h> // Contains definitions for B_Splines based on GSL. using namespace std; int main(){ double xmin=0.0, xmax=15.0; //vector<double> breakpts {0.0,1.0,2.,3.,4.,5.,6.,7., xmax,8.,9.,10.,11.}; // With this compiles/runs vector<double> breakpts {0.0, 5. , 1.0,2.,3.,4.,5.,6.,7.,xmax,8.,9., 3,10.,11.}; // With this crashes due to non increasing order size_t k=5; gsl_vector * gsl_breakpts = gsl_vector_alloc(breakpts.size()); for(size_t i=0; i<breakpts.size(); i++) gsl_vector_set(gsl_breakpts,i,breakpts[i]); size_t nbreak = breakpts.size(); size_t ncoeffs = nbreak + k - 2; gsl_bspline_workspace *bw = gsl_bspline_alloc(k,nbreak); gsl_vector * B = gsl_vector_alloc(ncoeffs); // gsl_bspline_knots_uniform(xmin, xmax, bw); // Works perfect gsl_bspline_knots(gsl_breakpts,bw); // With this, does not always recognise non increasing order of breakpoints. for(double xi=xmin; xi<xmax; xi +=0.01) { cout<< xi<< " "; gsl_bspline_eval(xi, B, bw); for(size_t j=0; j<ncoeffs;j++) cout<< gsl_vector_get(B,j)<< " "; cout<<endl; } gsl_vector_free(gsl_breakpts); gsl_vector_free(B); gsl_bspline_free(bw); /* bspline_basis * mybasis = new bspline_basis(breakpts,k); for(double x=0.0; x<7.;x+=0.01) { cout<<x<< " "; for(size_t i=0; i< mybasis->get_ncoeffs(); i++) cout<<mybasis->get_Bix(i,x)<<" "; cout<<endl; } delete mybasis; -*/ }