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.fit;
with gsl.multifit;
with gsl.rng ;
with gsl.randist;
with gsl.vector_double ;
with gsl.matrix_double ;

procedure Fit is
   use gsl.double_text_Io ;

   Verbose : boolean := true ;
   Status : Int ;

   procedure Test1 (weighted : boolean := true) is

      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;

-- 40.8.1 Simple Linear Regression Example

      x : gsl.double_array := ( 1970.0 , 1980.0 , 1990.0 , 2000.0 );
      y : gsl.double_array := (  12.0 ,   11.0 ,   14.0 ,   13.0 );
      w : gsl.double_array := (  0.1,  0.2,  0.3,  0.4 ) ;
      c0, c1 : aliased double ;
      cov00, cov01, cov11 : aliased double ;
      -- sumsq : aliased double ;
      chisq : aliased double ;
      N1 : constant := -30 ;
      N2 : constant := 130 ;

      xf, yf, yf_error : aliased double ;

   begin
      if Verbose
      then
         Put_Line(gnat.source_info.enclosing_entity);
      end if ;

      if weighted
      then 
         Put_Line("Weighted fit");
         Status := gsl.fit.wlinear(x , 1 , w , 1 , y , 1 , size_t(x'Length) , 
                c0'access , c1'access , 
                cov00'access , cov01'access , cov11'access , chisq'access ) ;
      else
         Status := gsl.fit.linear(x , 1 , y , 1 , size_t(x'Length) , 
                c0'access , c1'access , 
                cov00'access , cov01'access , cov11'access , chisq'access ) ;

      end if ;

      Put("c0 "); Put(Long_Float(c0)); Put(" c1 "); Put(Long_Float(c1)); New_Line ;
      Put(cov00); Put(" "); Put(cov01); Put(" "); Put(cov01); Put(" "); Put(cov11); New_Line ;
      Put("chisq "); Put(Long_Float(chisq)); New_Line ;

      if weighted
      then
         Create(logfile,Out_File,myname & ".w.csv");
      else
         Create(logfile,Out_File,myname & ".csv");
      end if ;
      Set_Output(logfile);
      for i in x'Range
      loop
         Put(x(i)); Put(" ; ");
         Put(y(i)); Put(" ; "); 
         Put(1.0/double(Sqrt(Float(w(i))))); Put(" ; "); 
         New_Line;
      end loop;
      Set_Output(Standard_Output);
      Close(logfile);

      Put_Line("Calculated");
      for i in x'Range
      loop
         Put(Long_FLoat(x(i))); Put(" ");
         Put(Long_Float(y(i))); Put(" ");
         Put(Long_Float(c0 + c1*x(i))); Put(" ");
         Put("error "); Put(Abs(Long_Float(y(i)-(c0+c1*x(i)))));
         New_Line;
      end loop;

      if weighted
      then 
         Create(logfile,Out_File,myname & ".w.est.csv");
      else
         Create(logfile,Out_File,myname & ".est.csv");
      end if ;

      Set_Output(logfile);

      for i in N1..N2
      loop
         xf := x(x'first) + (double(i)/100.0) * (x(x'Last) - x(x'first)) ;
         Put(xf); Put(" ; ");
         Status := gsl.fit.linear_est(double(xf) , c0 , c1 , 
                               cov00 , cov01 , cov11 , 
                               yf'access , yf_error'access);

            Put(yf); Put(" ; ");
            Put(yf_error); Put(" ; ");
            Put(yf-yf_error); Put(" ; ");
            Put(yf+yf_error) ; Put( " ; ");
            New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);
   end Test1 ;

   -- 40.8.2 Multi-parameter Linear Regression Example
   procedure Test2 is
      use gsl.double_elementary_functions ;
      N : constant := 19 ;
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      XData , cov : access gsl.matrix_double.gsl_matrix ;
      Y, W, c : access gsl.vector_double.gsl_vector ;

      ws : access gsl.multifit.gsl_multifit_linear_workspace ;
      chisq : aliased double ;
      xfrom, xto, xdelta, x : double ;

      yest, yerr : aliased double ;

      XCand : access gsl.vector_double.gsl_vector ;
      procedure CreateTestData is
         rng : access gsl.rng.gsl_rng ; 
         tdfile : File_Type ;

         y0, sigma, dy : double ;
         rownum : size_t ;
      begin
         rng := gsl.rng.alloc(gsl.rng.env_setup);
         xfrom := 0.0 ;
         xdelta := 0.1 ;
         xto := xfrom + double( N - 1 ) * xdelta ;
         Create( tdfile , Out_File , myname & ".td.csv");
         Set_Output(tdfile);
         x := xfrom ;
         rownum := 0 ;
         loop
            y0 := exp(x);
            sigma := 0.1 * y0 ;
            dy := gsl.randist.gaussian(rng,sigma);
            Put(x) ; Put(" ; ");
            Put(y0+dy); Put(" ; ");
            Put(sigma); Put(" ; ");
            New_Line;

            x := x + xdelta ;

            gsl.matrix_double.set(XData , rownum , 0 , 1.0 );
            gsl.matrix_double.set(XData , rownum , 1 , x );
            gsl.matrix_double.set(XData , rownum , 2 , x * x );
 
            gsl.vector_double.set (Y , rownum , y0+dy );
            gsl.vector_double.set (W , rownum , 1.0/(sigma*sigma)) ;
            rownum := rownum + 1 ;

            if x > xto
            then
               exit ;
            end if ; 

         end loop;
         Set_Output(Standard_Output);
         Close(tdfile);

         gsl.rng.free(rng);
      end CreateTestData;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;

      XData := gsl.matrix_double.alloc(N,3);
      cov := gsl.matrix_double.alloc(3,3) ;
      Y := gsl.vector_double.alloc(N);
      W := gsl.vector_double.alloc(N);
      c := gsl.vector_double.alloc(3);

      CreateTestData ;

      ws := gsl.multifit.linear_alloc(N,3);
      Status := gsl.multifit.wlinear (XData, W, Y, c, cov, chisq'Access, ws );

      Put("Chisq "); Put( chisq ) ; New_Line ;
      for i in 1..3
      loop
         for j in 1..3
         loop
            Put( gsl.matrix_double.get(cov,size_t(i-1),size_t(j-1)) ); Put (" ; ");
         end loop ;
         New_Line ;
      end loop;

      Create(logfile,Out_File,myname & ".est.csv");
      Set_Output(logfile);

      xfrom := 0.0 ;
      xdelta := 0.05 ;
      XCand := gsl.vector_double.alloc( size_t((xto-xfrom)/xdelta) );
      x := xfrom ;
      for n in 1..XCand.size
      loop
         gsl.vector_double.set(XCand,n-1,x);
         x := x + xdelta ;
      end loop ;
      declare
         xout, yout : double ;
         c0, c1, c2 : double ;
      begin
         c0 := gsl.vector_double.get(C,0);
         c1 := gsl.vector_double.get(C,1);
         c2 := gsl.vector_double.get(C,2);
        for n in 1..Integer(XCand.size)
         loop
            xout := gsl.vector_double.get(XCand,size_t(n-1)) ;
            Put( xout ); Put(" ; ");
            Put( c0 + c1 * xout + c2 * xout * xout ); Put(" ; ");
            New_Line;
         end loop ;
      end ;

      Set_Output(Standard_Output);
      Close(logfile);

      gsl.multifit.linear_free(ws);
      gsl.matrix_double.free(XData);
      gsl.matrix_double.free(cov);
      gsl.vector_double.free(Y);
      gsl.vector_double.free(W);
      gsl.vector_double.free(c);
   end Test2 ;


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;

   Test1(True);
   Test1(False);
   Test2;

end Fit;