#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;
-*/


}