-- test_psimeq_dbg.adb
with text_io; use text_io;
with real_arrays;
use real_arrays;
with psimeq_dbg;

procedure Test_Psimeq_Dbg is
  rms_err, avg_err, Max_Err : real;
  Idata : Integer := 1; -- reset after each use of the data set
  N_Data : Integer := 20; -- number of data points
  X_Data : Real_Vector(1..N_Data);
  Y_Data : Real_Vector(1..N_Data);
  Yx : Real_Vector(1..2); -- for passing y,x
  Debug : boolean := false;

  function Data_Set  return Integer is
    -- returns 1 for data, 0 for end
    -- sets value of Yx(1) to Y
    -- for value of Yx(2) for X
  begin
    if Idata > N_Data then
      Idata := 1;
      return 0; -- ready for check
    end if;
    Yx(1) := Y_Data(Idata);
    Yx(2) := X_Data(Idata);
    Idata := Idata +1;
    return 1;
  end Data_Set;

  procedure fit_pn(C : out Real_Vector) is
    A : Real_Matrix(C'range,C'range);
    Y : Real_Vector(C'range);
    Xx, Yy : real;
    Pwr : Real_Vector(C'range);
  begin
    for I in C'range loop
      for J in C'range loop
        A(I,J) := 0.0;
      end loop;
      Y(I) := 0.0;
    end loop;
    while data_set=1 loop
      Yy := Yx(1);
      Xx := Yx(2);
      pwr(1) := 1.0;
      for I in 2..C'Last loop
        pwr(I) := pwr(I-1)*Xx;
      end loop;
      for I in C'range loop
        for J in C'range loop
          A(I,J) := A(I,J) + pwr(I)*pwr(J);
        end loop;
        Y(I) := Y(I) + Yy*pwr(I);
      end loop;
    end loop; -- while
    C := psimeq_dbg(A, Y);
    if Debug then
      for I in C'range loop
        Put_Line("C("&Integer'Image(I)&"):="&Real'Image(C(I)));
      end loop;
    end if;
  end Fit_Pn;

  procedure check_pn(C : Real_Vector) is
    x, y, ya, Diff : Real;
    Sumsq : Real := 0.0;
    Sum   : Real := 0.0;
    Maxe  : Real := 0.0;
    Xmin  : Real := 0.0;
    Xmax  : Real := 0.0;
    Ymin  : Real := 0.0;
    Ymax  : Real := 0.0;
    Xbad  : Real := 0.0;
    Ybad  : Real := 0.0;
    K, Imax : Integer :=0;
    N : Integer := C'Last;
  begin
    k := 0;
    while data_set=1 loop
      y := yx(1);
      x := yx(2);
      if k=0  then
        Xmin := x;
        Xmax := x;
        Ymin := y;
        Ymax := y;
        Imax := 1;
        Xbad := x;
        Ybad := y;
      end if;
      if x>Xmax then xmax:=x; end if;
      if x<Xmin then xmin:=x; end if;
      if y>Ymax then ymax:=y; end if;
      if y<Ymin then ymin:=y; end if;
      K := K +1;
      ya := C(n)*x;
      for I in reverse 2..N-1 loop
        ya := (C(i)+ya)*x;
      end loop;
      ya := ya + C(1);
      diff := abs(y-ya);
      if diff>Maxe then
        Maxe := diff;
        Imax := k;
        Xbad := x;
        Ybad := y;
      end if;
      sum := sum + diff;
      sumsq := sumsq + diff*diff;
    end loop; -- while
    if Debug then
       Put_Line("check_pn k:="&Integer'Image(K)&
                  ", xmin:="&Real'Image(Xmin)&
                  ", xmax:="&Real'Image(Xmax)&
                  ", ymin:="&Real'Image(Ymin)&
                  ", ymax:="&Real'Image(ymax));
    end if;
    max_err := maxe;
    avg_err := sum/Real(K);
    rms_err := sumsq/Real(K);
    if Debug then
       Put_Line("max="&Real'Image(Max_Err)&" at i="&Integer'Image(Imax)&
                ", xbad="&Real'Image(Xbad)&", ybad="&Real'Image(ybad));
    end if;
  end check_pn;

begin
  Put_Line("test_psimeq_dbg.adb");
  for I in 1..N_Data loop
    X_Data(i) := Real(I-1)/10.0;
  end loop;
  y_data(1) := 0.0;
  y_data(2) := 6.0;
  y_data(3) := 14.1;
  y_data(4) := 5.01;
  y_data(5) := 4.326;
  for I in 6..N_Data loop
    y_data(i) := y_data(i-1);
  end loop;
  y_data(n_data) := 0.0;
  if Debug then
    for I in 1..N_Data loop
      Put_Line("i:="&Integer'Image(I)&", x:="&Real'Image(x_data(i)));
      Put_Line(", y:="&Real'Image(y_data(i)));
    end loop;
  end if;
  -- sample polynomial least square fit, nth power
  for N in 4..5 loop -- 5 was n_data
    declare
      C : Real_Vector(1..N);
    begin
      Put_Line("fit data to "&Integer'Image(N-1)&" degree polynomial");
      fit_pn(C);
      check_pn(C);
      for I in 1..N loop
        Put_Line("C("&Integer'Image(I)&"):="&Real'Image(C(i)));
      end loop;
      Put_Line("max_err:="&Real'Image(Max_Err)&", avg_err:="&
                 Real'Image(Avg_Err));
      Put_Line("rms_err:="&Real'Image(Rms_Err));
      if Debug then
        New_Line;
      end if;
    end;
  end loop;
end Test_Psimeq_Dbg;