with System ;
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.movstat ;
with gsl.rng ;
with gsl.randist ;
with gsl.vector_double ;
with math ;

with windowing ;

procedure Movstat is

   N : Integer := 500 ;        -- Length of time series
   K : Integer := 11 ;         -- Window Size

   NN : Integer := 1000 ;
   KK : Integer := 41 ;

   Status : Int ;
   function sign(x : double) return double is
   begin
      if x >= 0.0
      then
         return 1.0 ;
      else
         return -1.0 ;
      end if ;
   end sign ;

   procedure Test1 is

      ws : access gsl.movstat.gsl_movstat_workspace :=
                  gsl.movstat.alloc( size_t(K) );
      x : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(N) ) ;
      xmean : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc (size_t(N)) ;
      xmin : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(N) ) ;
      xmax : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc (size_t(N)) ;
      rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.default) ;
      ei : double ;
      xi : double ;
      csvfilename : String := gnat.source_info.enclosing_entity & ".csv" ;
      csvfile : File_Type ;
   begin
      Put_Line(gnat.source_info.enclosing_entity);
      Create(csvfile,Out_File,csvfilename);
      Set_Output(csvfile);
      for i in 1..N
      loop
         xi := double(cos(4.0 * math.M_PI * Float(i) / Float(N))) ;
         ei := gsl.randist.gaussian(rng,0.1);
         gsl.vector_double.set(x , size_t(i-1) , xi + ei );
         --Put(Long_Float(xi)); Put(" "); Put(Long_Float(ei)); Put(" "); New_Line;
      end loop ;
      Status := gsl.movstat.mean(gsl.movstat.GSL_MOVSTAT_END_PADVALUE,x , xmean, ws);
      Status := gsl.movstat.minmax(gsl.movstat.GSL_MOVSTAT_END_PADVALUE,x , xmin, xmax, ws);
      for i in 1..N
      loop
         Put(i) ; Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(x,size_t(i-1)))); Put("  ; ");
         Put(Long_Float(gsl.vector_double.get(xmean,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(xmin,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(xmax,size_t(i-1)))); Put(" ");
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(csvfile);
      gsl.vector_double.free(xmax);
      gsl.vector_double.free(xmin);
      gsl.vector_double.free(xmean);
      gsl.vector_double.free(x);
   end Test1 ;
   procedure Test2 is

      rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.default) ;
      ws : access gsl.movstat.gsl_movstat_workspace :=
                  gsl.movstat.alloc( size_t(KK) );

      x : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(NN) ) ;
      xmedian : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(NN) ) ;
      xmad : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(NN) ) ;
      xiqr : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(NN) ) ;
      xSn : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(NN) ) ;
      xQn : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(NN) ) ;
      xsd : access gsl.vector_double.gsl_vector :=
                  gsl.vector_double.alloc ( size_t(NN) ) ;
      csvfilename : String := gnat.source_info.enclosing_entity & ".csv" ;
      csvfile : File_Type ;
      sigma : gsl.double_array := (1.0, 5.0, 1.0, 3.0, 5.0) ;
      N_Sigma : array (1..5) of size_t := (200,450,600,850,1000) ;
      gi, u, outlier, xi : double ;
      varidx : size_t := 1 ;     -- variance changes
   begin
      Put_Line(gnat.source_info.enclosing_entity);
      Create(csvfile,Out_File,csvfilename);
      Set_Output(csvfile);
      for i in 1..NN
      loop
         gi := gsl.randist.gaussian(rng,sigma(Integer(varidx)-1)) ;
         u := gsl.rng.uniform(rng);
         if u < 0.01
         then
            outlier := 15.0 * sign(gi);
         else
            outlier := 0.0 ;
         end if ;
         xi := gi + outlier ;
         gsl.vector_double.set(x,size_t(i-1),xi);
         if size_t(i) = N_Sigma(Integer(varidx))
         then
            varidx := varidx + 1;
         end if ;
      end loop ;
      Status := gsl.movstat.mad(gsl.movstat.GSL_MOVSTAT_END_TRUNCATE, x, xmedian, xmad, ws);
      Status := gsl.movstat.qqr(gsl.movstat.GSL_MOVSTAT_END_TRUNCATE, x, 0.25, xiqr, ws);
      Status := gsl.movstat.Sn(gsl.movstat.GSL_MOVSTAT_END_TRUNCATE, x, xSn, ws);
      Status := gsl.movstat.Qn(gsl.movstat.GSL_MOVSTAT_END_TRUNCATE, x, xQn, ws);
      Status := gsl.movstat.sd(gsl.movstat.GSL_MOVSTAT_END_TRUNCATE, x, xsd, ws);
      Status := gsl.vector_double.scale(xiqr, 0.7413);

      varidx := 1;
      for i in 1..NN
      loop
         Put(i); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(x,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(sigma(Integer(varidx)-1))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(xmad,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(xiqr,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(xSn,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(xQn,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(xsd,size_t(i-1)))); Put(" ; ");
         New_Line;
         if size_t(i) = N_Sigma(Integer(varidx))
         then
            varidx := varidx + 1;
         end if ;
      end loop ;

      Set_Output(Standard_Output);
      Close(csvfile);
      gsl.vector_double.free(x);
      gsl.vector_double.free(xmedian);
      gsl.vector_double.free(xmad);
      gsl.vector_double.free(xiqr);
      gsl.vector_double.free(xSn);
      gsl.vector_double.free(xQn);
      gsl.vector_double.free(xsd);
      gsl.rng.free(rng);
      gsl.movstat.free(ws);
   end Test2 ;

   fn : aliased gsl.movstat.gsl_movstat_function ;

   procedure Test3 is
      proc : String := gnat.source_info.enclosing_entity ;
      N3 : constant := 1000 ;
      K3 : constant := 11 ;
      alpha : double := 0.1 ;
      x : access gsl.vector_double.gsl_vector := gsl.vector_double.alloc(N3) ;
      y : access gsl.vector_double.gsl_vector := gsl.vector_double.alloc(N3) ;
      rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.default);
      ws : access gsl.movstat.gsl_movstat_workspace :=
                  gsl.movstat.alloc( size_t(K3) );
      ui, outlier : double ;
      sum : double := 0.0 ;
      temp : double ;
      csvfile : File_Type ;
   begin
      Put_line(proc);
      for i in 1..N3
      loop
         ui := gsl.randist.gaussian(rng,1.0);
         temp := gsl.rng.uniform(rng);
         if temp < 0.01
         then
            outlier := sign(ui) * 10.0;
         else
            outlier := 0.0 ;
         end if;
         sum := sum + ui ;
         gsl.vector_double.set(x,size_t(i-1),sum+outlier);
      end loop;
      Create(csvfile,Out_File,proc & ".csv") ;
      Set_Output(csvfile);
      fn.c_function := windowing.windowingfn'access ;
      fn.params := alpha'Address ;
      Status := gsl.movstat.apply(gsl.movstat.GSL_MOVSTAT_END_PADVALUE , fn'access , x , y , ws );
      for i in 1..N3
      loop
         Put(i) ; Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(x,size_t(i-1)))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(y,size_t(i-1)))); Put(" ; ");
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(csvfile);
      gsl.vector_double.free(x);
      gsl.vector_double.free(y);
      gsl.movstat.free(ws);
      gsl.rng.free(rng);
   end Test3 ;
begin
   if Argument_Count > 0
   then
      N := Integer'Value( Argument(1) );
      if Argument_Count > 1
      then
         K := Integer'Value( Argument(2) );
      end if ;
   end if ;
   Ada.Long_Float_Text_IO.Default_aft := 4 ;
   Ada.Long_Float_Text_IO.Default_exp := 0 ;
   Test1 ;
   Test2 ;
   Test3 ;
end Movstat;