with Ada.Text_IO; use Ada.Text_IO; with Ada.Text_IO.C_Streams; with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO; with Ada.Integer_Text_IO; use Ada.Integer_Text_IO; with Interfaces.C; use Interfaces.C; with Interfaces.C.Strings; use Interfaces.C.Strings; with Ada.Strings.Unbounded; use Ada.Strings.Unbounded; with Ada.Command_Line; use Ada.Command_Line; with Ada.Numerics.Elementary_Functions ; use Ada.Numerics.Elementary_Functions ; with GNAT.Source_Info; use GNAT.Source_Info; with gsl; with gsl.vector_double; with gsl.matrix_double; with gsl.bspline; with gsl.rng; with gsl.randist; with gsl.version; procedure bspline is use gsl.double_text_Io ; use gsl.size_t_text_Io ; Verbose : boolean := true ; Status : Int ; procedure Test1 is -- 42.13.1 Example 1: Basis Splines and Uniform Knots myname : String := gnat.Source_Info.enclosing_entity ; logfile : File_Type ; procedure print_basis ( knotsfn : String ; splinefn : String ; order, nbreak : size_t ) is bspws : access gsl.bspline.gsl_bspline_workspace := gsl.bspline.alloc (order,nbreak); ncontrol : size_t := nbreak + order - 2; N : constant := 300 ; a, b, dx : double ; vec : access gsl.vector_double.gsl_vector ; u : double ; splinefile : File_Type ; begin a := 0.0 ; b := 1.0 ; dx := (b - a)/(double(N) - 1.0) ; vec := gsl.vector_double.alloc(ncontrol) ; Status := gsl.bspline.knots_uniform(a,b,bspws); gsl.vector_double.WriteCSV(bspws.knots,knotsfn); -- Status := gsl.vector_double.fprintf(bspws.knots,knotsfn & ".txt" , "%f"); Create(splinefile,Out_File,splinefn); Set_Output(splinefile); for i in 1..N loop u := double(i-1) * dx ; Status := gsl.bspline.eval(u,vec,bspws); Put(u) ; Put(" ; "); for j in 1..ncontrol loop Put( gsl.vector_double.get(vec,j-1) ); Put(" ; "); end loop ; New_Line ; end loop ; gsl.vector_double.free(vec) ; gsl.bspline.free(bspws); --Close(knotsfile); Set_Output(Standard_Output); Close(splinefile); end print_basis ; begin if Verbose then Put_Line(myname); end if ; print_basis(myname & ".linear.knots.csv" , myname & ".linear.spline.csv" , 2 , 11) ; print_basis(myname & ".quad.knots.csv" , myname & ".quad.spline.csv" , 3 , 11) ; print_basis(myname & ".cubic.knots.csv" , myname & ".cubic.spline.csv" , 4 , 11) ; print_basis(myname & ".quartic.knots.csv" , myname & ".quartic.spline.csv" , 5 , 11) ; end Test1 ; procedure Test2 is myname : String := gnat.Source_Info.enclosing_entity ; logfile : File_Type ; nbreak : constant := 6 ; spline_order : constant := 4 ; bspws : access gsl.bspline.gsl_bspline_workspace := gsl.bspline.alloc(spline_order,nbreak); ncontrol : size_t := nbreak + spline_order - 2; n : constant := 300; a : constant double := 0.0; b : constant double := 1.0; dx : constant double := (b - a) / (double(n) - 1.0); dB : access gsl.matrix_double.gsl_matrix := gsl.matrix_double.alloc(ncontrol, spline_order); lnum : Integer := 0; order_name : String := "0123" ; begin if Verbose then Put_Line(myname); end if ; Status := gsl.bspline.knots_uniform(a,b,bspws) ; gsl.vector_double.WriteCSV(bspws.knots,myname & ".knots.csv"); for i in 1..spline_order loop Create(logfile,Out_File,myname & "." & order_name(i) & ".csv"); Set_Output(logfile); for j in 1..n loop Status := gsl.bspline.deriv_eval( double(j-1) * dx , size_t(i-1) , dB , bspws ); lnum := lnum + 1 ; Put(lnum) ; Put(" ; "); Put(double(j-1) * dx ); Put (" ; "); for k in 1..ncontrol loop Put( gsl.matrix_double.get(dB,k-1,size_t(i-1))); Put(" ; "); end loop ; New_Line ; end loop ; Set_Output(Standard_Output); Close(logfile); end loop ; gsl.matrix_double.free(dB) ; gsl.bspline.free(bspws); end Test2 ; procedure Test3 is -- 2.13.3 Example 3: Least squares and breakpoints myname : String := gnat.Source_Info.enclosing_entity ; logfile : File_Type ; N : constant := 500 ; spline_order : constant := 4 ; a : constant double := 0.0 ; b : constant double := 15.0 ; sigma : constant double := 0.2 ; bspws40 : access gsl.bspline.gsl_bspline_workspace := gsl.bspline.alloc(spline_order,40); bspws10 : access gsl.bspline.gsl_bspline_workspace := gsl.bspline.alloc(spline_order,10); ncontrol40 : constant size_t := gsl.bspline.ncontrol(bspws40) ; ncontrol10 : constant size_t := gsl.bspline.ncontrol(bspws10) ; dof40 : constant size_t := N - ncontrol40; dof10 : constant size_t := N - ncontrol10; c40 : access gsl.vector_double.gsl_vector := gsl.vector_double.alloc(ncontrol40) ; c10 : access gsl.vector_double.gsl_vector := gsl.vector_double.alloc(ncontrol10) ; x : access gsl.vector_double.gsl_vector := gsl.vector_double.alloc(N) ; y : access gsl.vector_double.gsl_vector := gsl.vector_double.alloc(N) ; w : access gsl.vector_double.gsl_vector := gsl.vector_double.alloc(N) ; rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.env_setup); xi,yi,wi,dyi : double ; chisq40, chisq10 : aliased double ; begin if Verbose then Put_Line(myname); end if ; for i in 1..N loop xi := a + (b - a) / (double(N) - 1.0) * double(i-1) ; yi := double(Cos(Float(xi)) * Exp( - 0.1 * float(xi) )); dyi := gsl.randist.gaussian(rng,sigma); yi := yi + dyi ; gsl.vector_double.set(x,size_t(i-1),xi); gsl.vector_double.set(y,size_t(i-1),yi); gsl.vector_double.set(w,size_t(i-1),1.0/(sigma*sigma)); end loop ; gsl.vector_double.WriteCSV(x, myname & ".xyw.csv" , y , w ); Status := gsl.bspline.knots_uniform(a,b,bspws40); Status := gsl.bspline.knots_uniform(a,b,bspws10); Create(logfile,Out_File,myname & ".csv"); Set_Output(logfile); Set_Output(Standard_Output); Close(logfile); gsl.vector_double.free(c40); gsl.vector_double.free(c10); gsl.vector_double.free(x) ; gsl.vector_double.free(y) ; gsl.vector_double.free(w) ; gsl.bspline.free(bspws40); gsl.bspline.free(bspws10); end Test3 ; begin Ada.Long_Float_Text_IO.Default_exp := 0 ; Ada.Long_Float_Text_IO.Default_aft := 4 ; gsl.double_text_io.Default_aft := 4 ; gsl.double_text_Io.Default_exp := 0; gsl.version.Report ; Test1; Test2; end bspline ;