-- test_lsfit6.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 Array6d;
use  Array6d;
with Integer_Arrays; -- types integer_vector
use  Integer_Arrays;
with Lsfit;
use  Lsfit;

procedure Test_Lsfit6 is
  Nx : Integer := 4;      -- number of data in each dimension
  Ny : Integer := 4;
  Nz : Integer := 4;
  Nt : Integer := 4;
  Nu : Integer := 4;
  Nv : Integer := 4;
  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
  Vg : Real_Vector(1..Nv);    -- v grid
  Wg : Real_Matrix6d(1..Nx,1..Ny,1..Nz,1..Nt,1..Nu,1..Nv);
  Wa : Real_Matrix6d(1..Nx,1..Ny,1..Nz,1..Nt,1..Nu,1..Nv);
                     --computed solution at xyztuv grid,
  Hx, Hy, Hz, Ht, Hu, Hv : 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;
  vmin : Real := -1.0;
  vmax : Real :=  1.0;
  X, Y, Z, T, U, V, W, We : Real;
  Last : Boolean := False;
  Err, Maxerr : Real := 0.0;

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

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

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

begin
  Put_Line("test_lsfit6.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;
  hv := (vmax-vmin)/real(nv-1);
  Put_Line("nv="&Integer'Image(nv)&", step="& Real'Image(hv));
  for iiiiii in 1..nv loop
    vg(iiiiii) := vmin+Real(iiiiii-1)*hv;
  end loop;

  Put_Line("power data 1 test 6D");
  Fit_Init(1, 6); -- first power, six independent variables

  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);
            for iiiiii in 1..nv loop
              V := Vg(iiiiii);
              W := U1(X,Y,Z,T,U,V);
              if i=nx and ii=ny and iii=nz and iiii=nt and
                 iiiii=nu and iiiiii=nv then
                 Last := True;
              end if;
              Fit6_Data(X, Y, Z, T, U, V, W, Last);
            end loop;
          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);
            for iiiiii in 1..nv loop
              V := Vg(iiiii);
              W := U1(X,Y,Z,T,U,V);
              We := Eval6_Poly(X, Y, Z, T, U, V);
              Err := abs(W-We);
              if Err>Maxerr then Maxerr := Err; end if;
            end loop;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  Put_Line("power data 1 test 6D, maxerr="&Real'Image(Maxerr));
  New_Line;

  Put_Line("power Bpoly 1 test 6D");
  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
            for iiiiii in 1..nv loop
              Wa(i,ii,iii,iiii,iiiii,iiiiii) :=
                  ub(xg(i),yg(ii),zg(iii),tg(iiii),ug(iiiii),vg(iiiiii));
              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 or iiiiii=nv then
                Wg(i,ii,iii,iiii,iiiii,iiiiii) := ub(xg(i),yg(ii),zg(iii),
                                              tg(iiii),ug(iiiii),ug(iiiiii));
              else
                Wg(i,ii,iii,iiii,iiiii,iiiiii) := 0.0;
              end if;
            end loop;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  New_Line;

  Fit_Init(1, 6); -- first power, sixth dimension = 6 variables

  Fit6_Bpoly(Nx, Ny, Nz, Nt, Nu, Nv, Xg, Yg, Zg, Tg, Ug, Vg, 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);
            for iiiiii in 2..nv-1 loop
              V := Vg(iiiiii);
              W := Wa(i,ii,iii,iiii,iiiii,iiiiii);
              We := Eval6_Poly(X, Y, Z, T, U, V);
              Err := abs(W-We);
              if Err>Maxerr then Maxerr := Err; end if;
            end loop;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  Put_Line("power Bpoly 1 test 6D, maxerr="&Real'Image(Maxerr));
  New_Line;

  Put_Line("power data 6 test 6D, sin");
  Fit_Init(6, 6); -- sixth power, six independent variables

  Last := False;
  for i in 1..nx loop
    x := 0.04*real(i);
    for ii in 1..ny loop
      y := 0.04*real(ii);
      for iii in 1..nz loop
        z := 0.04*real(iii);
        for iiii in 1..nz loop
          t := 0.04*real(iiii);
          for iiiii in 1..nu loop
            u := 0.04*real(iiiii);
            for iiiiii in 1..nv loop
              v := 0.04*real(iiiiii);
              W := sin(X+Y+Z+T+U+V);
              if i=nx and ii=ny and iii=nz and iiii=nt and
                 iiiii=nu and iiiiii=nv then
                Last := True;
              end if;
              Fit6_Data(X, Y, Z, T, U, V, W, Last);
            end loop;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  -- built, now check
  for i in 1..nx loop
    x := 0.04*real(i);
    for ii in 1..ny loop
      y := 0.04*real(ii);
      for iii in 1..nz loop
        z := 0.04*real(iii);
        for iiii in 1..nt loop
          t := 0.04*real(iiii);
          for iiiii in 1..nu loop
            u := 0.04*real(iiiii);
            for iiiiii in 1..nv loop
              v := 0.04*real(iiiiii);
              W := sin(X+Y+Z+T+U+V);
              We := Eval6_Poly(X, Y, Z, T, U, V);
              Err := abs(W-We);
              if Err>Maxerr then Maxerr := Err; end if;
            end loop;
          end loop;
        end loop;
      end loop;
    end loop;
  end loop;
  Put_Line("power data 6 test on sin, maxerr="&Real'Image(Maxerr));


  Put_Line("test_lsfit6.adb finished");
end Test_Lsfit6;