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

procedure Filter is
   Verbose : boolean := true;
   N : size_t := 500 ;          -- Time Series length
   K : size_t := 51 ;           -- Kernel Size
   alpha : gsl.double_array := ( 0.5 , 3.0 , 10.0 );
   Status : Int ;

   procedure Test1 is
   -- 
   -- 24.5.1 Gaussian Example 1
   --
      rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.default) ;
      x, y1, y2, y3, k1, k2, k3 : access gsl.vector_double.gsl_vector ;
      gws : access gsl.filter.gsl_filter_gaussian_workspace :=
               gsl.filter.gaussian_alloc (K) ;
      sum : double := 0.0 ;
      ui : double ;
      kernelfile, datafile : File_Type ;
   begin
      Put_Line(enclosing_entity);
      x := gsl.vector_double.alloc(N);
      y1 := gsl.vector_double.alloc(N);
      y2 := gsl.vector_double.alloc(N);
      y3 := gsl.vector_double.alloc(N);
      k1 := gsl.vector_double.alloc(K);
      k2 := gsl.vector_double.alloc(K);
      k3 := gsl.vector_double.alloc(K);

      for i in 0..N-1
      loop
         ui := gsl.randist.gaussian(rng,1.0);
         sum := sum + ui ;
         gsl.vector_double.set(x, i , sum );
      end loop ;

      Status := gsl.filter.gaussian_kernel(alpha(0) , 0 , 0 , k1 );
      Status := gsl.filter.gaussian_kernel(alpha(1) , 0 , 0 , k2 );
      Status := gsl.filter.gaussian_kernel(alpha(2) , 0 , 0 , k3 );

      Status := gsl.filter.gaussian(gsl.filter.GSL_FILTER_END_PADVALUE,alpha(0),0,x,y1,gws) ;
      Status := gsl.filter.gaussian(gsl.filter.GSL_FILTER_END_PADVALUE,alpha(1),0,x,y2,gws) ;
      Status := gsl.filter.gaussian(gsl.filter.GSL_FILTER_END_PADVALUE,alpha(2),0,x,y3,gws) ;

      Ada.Long_Float_Text_IO.Default_aft := 4 ;
      Ada.Long_Float_Text_IO.Default_exp := 0 ;

      Create(kernelfile,Out_File,enclosing_entity & ".kernel.csv");
      Set_Output(kernelfile);
      for i in 0..K-1
      loop
         Put(Integer(i)); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(k1,i))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(k2,i))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(k3,i))); Put(" ; ");
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(kernelfile);

      Create(datafile,Out_File,enclosing_entity & ".data.csv");
      Set_Output(datafile);
      for i in 0..N-1
      loop
         Put(Integer(i)); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(x,i))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(y1,i))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(y2,i))); Put(" ; ");
         Put(Long_Float(gsl.vector_double.get(y3,i))); Put(" ; ");
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(datafile);
   end Test1 ;

   procedure Test2 is
   --
   -- 24.5.2 Gaussian Example 2
   --
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;

      x ,
      dxi ,
      y ,
      dy ,
      d2y : access gsl.vector_double.gsl_vector ;

      gws : access gsl.filter.gsl_filter_gaussian_workspace ;
      rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.default);
      xi, ei : double ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      N := 1000 ;
      x := gsl.vector_double.alloc(N) ;
      y := gsl.vector_double.alloc(N) ;
      dy := gsl.vector_double.alloc(N) ;
      d2y := gsl.vector_double.alloc(N) ;
      dxi := gsl.vector_double.alloc(N) ;
      gws := gsl.filter.gaussian_alloc (K) ;

      for i in 1..N
      loop
         if i > N/2
         then
            xi := 0.5;
         else 
            xi := 0.0 ;
         end if;
         ei := gsl.randist.gaussian(rng,0.1);
         gsl.vector_double.set(x,size_t(i-1),xi+ei);
      end loop ;

   for i in 1..N
      loop
         if i = 1
         then
            gsl.vector_double.set(dxi,i-1,
                                  gsl.vector_double.get(x,i) - gsl.vector_double.get(x,i-1));
         elsif i = n
         then
            gsl.vector_double.set(dxi,i-1,gsl.vector_double.get(x,i-1) - gsl.vector_double.get(x,i-2));
         else
            gsl.vector_double.set(dxi,i-1,0.5*(gsl.vector_double.get(x,i) - gsl.vector_double.get(x,i-2)));
         end if;

      end loop ;

      gsl.vector_double.WriteCSV(x,myname & ".data.csv");

      Status := gsl.filter.gaussian(gsl.filter.GSL_FILTER_END_PADVALUE,alpha(2),0,x,y,gws);
      Status := gsl.filter.gaussian(gsl.filter.GSL_FILTER_END_PADVALUE,alpha(2),1,x,dy,gws);
      Status := gsl.filter.gaussian(gsl.filter.GSL_FILTER_END_PADVALUE,alpha(2),2,x,d2y,gws);

      gsl.vector_double.WriteCSV( x , myname & ".filtered.csv" ,
                                 y , dy , d2y , dxi );

      gsl.vector_double.free(x);
      gsl.vector_double.free(y);
      gsl.vector_double.free(dy);
      gsl.vector_double.free(d2y);
      gsl.vector_double.free(dxi);
      gsl.rng.free(rng);
      gsl.filter.gaussian_free(gws);
   end Test2 ;
   procedure Test3 is
   -- 24.5.3 Square Wave Signal Example
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;

      gws : access gsl.filter.gsl_filter_gaussian_workspace ;
      rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.default);

      t,
      x,
      y_median,
      y_rmedian : access gsl.vector_double.gsl_vector ;

      medianws : access gsl.filter.gsl_filter_median_workspace ;
      rmedianws : access gsl.filter.gsl_filter_rmedian_workspace ;

      freq : double := 5.0 ;
      ti : double ;
      tmp : double ;
      xi, ei : double ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      N := 1000 ;
      K := 7 ;
      gws := gsl.filter.gaussian_alloc(K);
      t := gsl.vector_double.alloc(N);
      x := gsl.vector_double.alloc(N);
      y_median := gsl.vector_double.alloc(N) ;
      y_rmedian := gsl.vector_double.alloc(N) ;


      medianws := gsl.filter.median_alloc(K) ;
      rmedianws := gsl.filter.rmedian_alloc(K) ;
   
      for i in 1..N
      loop
         ti := double(i) / double(N-1) ;
         tmp := double(Ada.Numerics.Elementary_Functions.Sin( 2.0 * Ada.Numerics.PI * float(freq) * float(ti)));
         if tmp >= 0.0
         then
            xi := 1.0 ;
         else
            xi := -1.0 ;
         end if ;
         ei := gsl.randist.gaussian(rng,0.1);
         gsl.vector_double.set(t,i-1,ti);
         gsl.vector_double.set(x,i-1,xi+ei);
      end loop ;
      Status := gsl.filter.median(gsl.filter.GSL_FILTER_END_PADVALUE,x,y_median,medianws);
      Status := gsl.filter.rmedian(gsl.filter.GSL_FILTER_END_PADVALUE,x,y_rmedian,rmedianws);
      gsl.vector_double.WriteCSV(t,myname & ".csv" ,x,y_median,y_rmedian) ;

      gsl.vector_double.free(t);
      gsl.vector_double.free(x);
      gsl.vector_double.free(y_median) ;
      gsl.vector_double.free(y_rmedian) ;

      gsl.filter.median_free(medianws);
      gsl.filter.rmedian_free(rmedianws);

   end Test3 ;

   procedure Test4 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      t : double ;

      x,
      y,
      xmedian ,
      xsigma ,
      lupper,
      llower,
      watch : access gsl.vector_double.gsl_vector ;
      ioutlier : access gsl.vector_int.gsl_vector_int ;
      rng : access gsl.rng.gsl_rng := gsl.rng.alloc(gsl.rng.default) ;
      iws : access gsl.filter.gsl_filter_impulse_workspace ;

      xi, ei, u, outlier : double ;
      noutlier : aliased size_t ;
   begin
      -- 24.5.4 Impulse Detection Example
      if Verbose
      then
         Put_Line(myname);
      end if ;
      N := 1000 ;
      K := 25 ;
      t := 4.0 ;
      x := gsl.vector_double.alloc(N);
      y := gsl.vector_double.alloc(N);
      xmedian := gsl.vector_double.alloc(N);
      xsigma := gsl.vector_double.alloc(N);
      ioutlier := gsl.vector_int.alloc(N) ;
      lupper := gsl.vector_double.alloc(N) ;
      llower := gsl.vector_double.alloc(N) ;
      watch := gsl.vector_double.alloc(N) ;
      iws := gsl.filter.impulse_alloc(K) ;

      for i in 1..N
      loop
         u := gsl.rng.uniform(rng) ;
         ei := gsl.randist.gaussian(rng,2.0);
         xi := 10.0 * 
                  double(Ada.Numerics.Elementary_Functions.sin( 2.0 * 
                                                         Ada.Numerics.PI * 
                                                         Float(i) / Float (N))) ;
         if u < 0.01
         then
            outlier := 15.0 * double(gsl.math.GSL_SIGN(ei)) ;
         else
            outlier := 0.0;
         end if ;

         gsl.vector_double.set(x, i-1 , xi+ei+outlier);
      
      end loop ;
      Status := gsl.filter.impulse(gsl.filter.GSL_FILTER_END_TRUNCATE, 
                         gsl.filter.GSL_FILTER_SCALE_QN, t, x, y,
                         xmedian, xsigma, noutlier'access , ioutlier, iws );

      for i in 1..N
      loop
         gsl.vector_double.set(lupper,i-1,
                           gsl.vector_double.get(xmedian,i-1) + t * gsl.vector_double.get(xsigma,i-1));
         gsl.vector_double.set(llower,i-1,
                           gsl.vector_double.get(xmedian,i-1) - t * gsl.vector_double.get(xsigma,i-1));
         gsl.vector_double.set
           (watch,
            i - 1,
            gsl.vector_double.get(xmedian,i-1) * double(gsl.vector_int.get(ioutlier,i-1)) ); 
      end loop ;
      gsl.vector_double.WriteCSV
        (x,
         myname & ".csv",
         y , xmedian, xsigma, llower, lupper, watch );
      gsl.vector_int.WriteCSV
         (ioutlier, myname & ".outlier.csv" );

      gsl.vector_double.free(x);
      gsl.vector_double.free(y);
      gsl.vector_double.free(xmedian);
      gsl.vector_double.free(xsigma );
      gsl.vector_double.free(lupper);
      gsl.vector_double.free(llower);
      gsl.vector_int.free(ioutlier);
      gsl.rng.free(rng) ;
      gsl.filter.impulse_free(iws) ;
   end Test4 ;

   procedure Template is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      Set_Output(Standard_Output);
      Close(logfile);
   end Template ;
begin
   if Argument_Count > 0
   then
      N := size_t'Value( Argument(1) );
      if Argument_Count > 1
      then
         K := size_t'Value( Argument(2) );
      end if ;
   end if ;
   Template ;
   Test1 ;
   Test2 ;
   Test3 ;
   Test4 ;
end Filter;