-- test_lsfit5.adb

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Real_Arrays; -- types real, real_vector, real_matrix
use  Real_Arrays;
with Array5d;
use  Array5d;
with Integer_Arrays; -- types integer_vector
use  Integer_Arrays;
with Lsfit;
use  Lsfit;

procedure Test_Lsfit5 is
  Nx : Integer := 5;      -- number of data in each dimension
  Ny : Integer := 5;
  Nz : Integer := 5;
  Nt : Integer := 5;
  Nu : Integer := 5;
  Xg : Real_Vector(1..Nx);    -- x grid, does not need to be uniform
  Yg : Real_Vector(1..Ny);    -- y grid
  Zg : Real_Vector(1..Nz);    -- z grid
  Tg : Real_Vector(1..Nt);    -- t grid
  Ug : Real_Vector(1..Nu);    -- u grid
  Wg : Real_Matrix5d(1..Nx,1..Ny,1..Nz,1..Nt,1..Nu);
  Wa : Real_Matrix5d(1..Nx,1..Ny,1..Nz,1..Nt,1..Nu);
                     --computed solution at xyztu grid,
  Hx, Hy, Hz, Ht, Hu : real;
  xmin : Real := -1.0;
  xmax : Real :=  1.0;
  ymin : Real := -1.0;
  ymax : Real :=  1.0;
  zmin : Real := -1.0;
  zmax : Real :=  1.0;
  tmin : Real := -1.0;
  tmax : Real :=  1.0;
  umin : Real := -1.0;
  umax : Real :=  1.0;
  X, Y, Z, T, U, W, We : Real;
  Last : Boolean := False;
  Err, Maxerr : Real := 0.0;

  function ub(x:Real; y:Real; z:Real; t:Real; u:real) return real is
                                                      -- value
   begin
     return 1.0 + 2.0*x + 3.0*y + 4.0*z +5.0*t +6.0*u;
  end ub;

  function u1(x:Real; y:Real; z:Real; t:Real; u:real) return real is
                                                      -- value
   begin
     return 1.0 + 2.0*x + 3.0*y + 4.0*z +5.0*t +6.0*u;
  end u1;

  package Real_IO is new Float_IO(Real);
  use Real_IO;

begin
  Put_Line("test_lsfit5.adb");

  hx := (xmax-xmin)/real(nx-1);
  Put_Line("nx="&Integer'Image(nx)&", step="&Real'Image(hx));
  for i in 1..nx loop
    xg(i) := xmin+Real(i-1)*hx;
  end loop;
  hy := (ymax-ymin)/real(ny-1);
  Put_Line("ny="&Integer'Image(ny)&", step="&Real'Image(hy));
  for ii in 1..ny loop
    yg(ii) := ymin+Real(ii-1)*hy;
  end loop;
  hz := (zmax-zmin)/real(nz-1);
  Put_Line("nz="&Integer'Image(nz)&", step="& Real'Image(hz));
  for iii in 1..nz loop
    zg(iii) := zmin+Real(iii-1)*hz;
  end loop;
  ht := (tmax-tmin)/real(nt-1);
  Put_Line("nt="&Integer'Image(nt)&", step="& Real'Image(ht));
  for iiii in 1..nt loop
    tg(iiii) := tmin+Real(iiii-1)*ht;
  end loop;
  hu := (umax-umin)/real(nu-1);
  Put_Line("nu="&Integer'Image(nu)&", step="& Real'Image(hu));
  for iiiii in 1..nu loop
    ug(iiiii) := umin+Real(iiiii-1)*hu;
  end loop;

  Put_Line("power 1 data test 5D");
  Fit_Init(1, 5); -- first power, five independent variables
  Maxerr := 0.0;

  Last := False;
  for i in 1..nx loop
    X := Xg(i);
    for ii in 1..ny loop
      Y := Yg(ii);
      for iii in 1..nz loop
        Z := Zg(iii);
        for iiii in 1..nz loop
          T := Tg(iiii);
          for iiiii in 1..nu loop
            U := Ug(iiiii);
            W := U1(X,Y,Z,T,U);
            if i=nx and ii=ny and iii=nz and iiii=nt and iiiii=nu then
               Last := True;
            end if;
            Fit5_Data(X, Y, Z, T, U, W, Last);
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  -- built, now check
  for i in 1..nx loop
    X := Xg(i);
    for ii in 1..ny loop
      Y := Yg(ii);
      for iii in 1..nz loop
        Z := Zg(iii);
        for iiii in 1..nt loop
          T := Tg(iiii);
          for iiiii in 1..nu loop
            U := Ug(iiiii);
            W := U1(X,Y,Z,T,U);
            We := Eval5_Poly(X, Y, Z, T, U);
            Err := abs(W-We);
            if Err>Maxerr then Maxerr := Err; end if;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  Put_Line("power 1 data test 5D, maxerr="&Real'Image(Maxerr));
  New_Line;

  Put_Line("power 1 Bpoly test 5D");
  Maxerr := 0.0;
  for i in 1..nx loop
    for ii in 1..ny loop
      for iii in 1..nz loop
        for iiii in 1..nt loop
          for iiiii in 1..nu loop
            Wa(i,ii,iii,iiii,iiiii) :=
                               ub(xg(i),yg(ii),zg(iii),tg(iiii),ug(iiiii));
            if i=1 or i=nx or ii=1 or ii=ny or
               iii=1 or iii=nz or iiii=1 or iiii=nt or
               iiiii=1 or iiiii=nu then
              Wg(i,ii,iii,iiii,iiiii) := ub(xg(i),yg(ii),zg(iii),
                                            tg(iiii),ug(iiiii));
            else
              Wg(i,ii,iii,iiii,iiiii) := 0.0;
            end if;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  New_Line;

  Fit_Init(1, 5); -- first power, fifth dimension = 5 variables

  Fit5_Bpoly(Nx, Ny, Nz, Nt, Nu, Xg, Yg, Zg, Tg, Ug, Wg);

  -- test fit
  for i in 2..nx-1 loop
     X := Xg(i);
     for ii in 2..ny-1 loop
      Y := Yg(ii);
      for iii in 2..nz-1 loop
        Z := Zg(iii);
        for iiii in 2..nt-1 loop
          T := Tg(iiii);
          for iiiii in 2..nu-1 loop
            U := Ug(iiiii);
            W := Wa(i,ii,iii,iiii,iiiii);
            We := Eval5_Poly(X, Y, Z, T, U);
            Err := abs(W-We);
            if Err>Maxerr then Maxerr := Err; end if;
            Put("expected=");
            Put(Wa(i,ii,iii,iiii,iiiii),8,3,0);
            Put("  fit=");
            Put(We,8,3,0);
            Put("  error="&Real'Image(W-We));
            New_Line;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  Put_Line("power 1 Bpoly test 5D, maxerr="&Real'Image(Maxerr));
  New_Line;
  Put_Line("test_lsfit5.adb finished");
end Test_Lsfit5;