-- test_lsfit4.adb  using stuff from  pde_nl44.adb  for test

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 Array4d;
use  Array4d;
with Integer_Arrays; -- types integer_vector
use  Integer_Arrays;
with Lsfit;
use  Lsfit;

procedure Test_Lsfit4 is
  Nx : Integer := 5;      -- number of data in each dimension
  Ny : Integer := 5;
  Nz : Integer := 5;
  Nt : 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);    -- z grid
  Ug : Real_Matrix4d(1..Nx,1..Ny,1..Nz,1..Nt); --computed solution at xyz grid,
  Ua : Real_matrix4d(1..nx,1..ny,1..Nz,1..nt); --known includes boundary
  Hx, Hy, Hz, Ht : 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;
  X, Y, Z, T, U, Ue : Real;
  Last : Boolean := False;
  Err, Maxerr : Real := 0.0;

  function ub(x:Real; y:Real; z:Real; t:Real) return real is -- boundary
  begin
     return x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*t*t*t*t+
            5.0*x*x*x*y+6.0*y*y*y*z+7.0*z*z*z*t+8.0*t*t*t*x+
            9.0*x + 10.0*y + 11.0*z + 12.0*t + 13.0;
  end ub;

  function U1(X:Real; Y:Real; Z:Real; T:Real) return Real is -- data
  begin
    return 1.0 + 2.0*X + 3.0*Y + 4.0*Z + 5.0*T;
  end U1;

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

begin
  Put_Line("test_lsfit4.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;

  Fit_Init(1, 4); -- first power, four 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);
          U := U1(X,Y,Z,T);
          if i=nx and ii=ny and iii=nz and iiii=nt then
             Last := True;
          end if;
          Fit4_Data(X, Y, Z, T, U, Last);
        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);
          U := U1(X,Y,Z,T);
          Ue := Eval4_Poly(X, Y, Z, T);
          Err := abs(U-Ue);
          if Err>Maxerr then Maxerr := Err; end if;
        end loop;
      end loop;
    end loop;
  end loop;
  Put_Line("power 1 test, maxerr="&Real'Image(Maxerr));

  for i in 1..nx loop
    for ii in 1..ny loop
      for iii in 1..nz loop
        for iiii in 1..nt loop
          Ua(i,ii,iii,iiii) := ub(xg(i),yg(ii),zg(iii),tg(iiii));
          Put("xg=");
          Put(xg(i),3,4,0);
          Put(", yg=");
          Put(yg(ii),3,4,0);
          Put(", zg=");
          Put(zg(iii),3,4,0);
          Put(", tg=");
          Put(tg(iiii),3,4,0);
          Put(", Ua("&Integer'Image(i)&","&Integer'Image(ii)&
              ","&Integer'Image(iii)&","&Integer'Image(iiii));
          Put(Ua(i,ii,iii,iiii),3,5,0);
          New_Line;
          if i=1 or i=nx or ii=1 or ii=ny or
             iii=1 or iii=nz or iiii=1 or iiii=nt then
            ug(i,ii,iii,iiii) := ub(xg(i),yg(ii),zg(iii),tg(iiii));
          else
            ug(i,ii,iii,iiii) := 0.0;
          end if;
        end loop;
      end loop;
    end loop;
  end loop;
  New_Line;

  Fit_Init(4,4); -- max power, dimension

  Fit4_Bpoly(Nx, Ny, Nz, Nt, Xg, Yg, Zg, Tg, Ug);

  -- 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);
          U := Eval4_Poly(X, Y, Z, T);
          Put("expected=");
          Put(Ua(i,ii,iii,iiii),8,3,0);
          Put("  fit=");
          Put(U,8,3,0);
          Put("  error="&Real'Image(Ua(i,ii,iii,iiii)-U));
          New_Line;
        end loop;
      end loop;
    end loop;
  end loop;
  New_Line;
  Put_Line("test_lsfit4.adb finished");
end Test_Lsfit4;