File : fouriers-io.adb

with Messages, Reps.Ops, Reps.IO;
use Messages;

pragma Elaborate_All (Messages, Reps.IO, Reps.Ops);

package body Fouriers.IO is

  package Txt_Int_IO is new Ada.Text_IO.Integer_IO(Integer);
  use Txt_Int_IO;

  function IMax(I1,I2: Integer) return Integer renames Integer'Max;

  procedure EffDeg(H: in Fourier; Dm,Dn: out Integer) is
  begin
    Dm := -1;
    Dn := -1;
    for M in Freq1 loop
      for N in Nmin(M) .. Ncut(M) loop
        if not IsZero(H.C(M,N)) then
          Dm := Imax(Abs(M),Dm);
          Dn := Imax(N,Dn);
        end if;
      end loop;
    end loop;
  end EffDeg;

  procedure Show(S: in String; R: in Weights; Decimal: in Boolean := True) is
    use Reps.IO;
  begin
    Put(Current_Output,S & " ( ");
    Txt_Put(Current_Output,R.Q,Decimal);
    Put(" , ");
    Txt_Put(Current_Output,R.P,Decimal);
    Put_Line(" )");
  end Show;

  procedure Show(L: in Degrees) is
  begin
    Put("  \ M");
    for M in Freq1 loop
      Put(Current_Output,M,3,10);
    end loop;
    New_Line;
    Put("  N \_");
    for M in Freq1 loop
      Put("___");
    end loop;
    New_Line;
    for N in reverse Freq2 loop
      Put(Current_Output,N,3,10);
      Put(" |");
      for M in Freq1 loop
        if Kcut(M,N)<0 then
          Put("   ");
        elsif L(M,N)<0 then
          Put("  +");
        else
          Put(Current_Output,L(M,N),3,10);
        end if;
      end loop;
      New_Line;
    end loop;
    New_Line;
  end Show;

  procedure ShowIndices(J: in Parity; M: in Freq1; N: in Freq2; K: in Power) is
  begin
    if J=0 then
      Put(Current_Output,"cos(");
    else
      Put(Current_Output,"sin(");
    end if;
    Put(Current_Output,M,3,10);
    Put(Current_Output,"x+");
    Put(Current_Output,N,2,10);
    Put(Current_Output,"y) z**");
    Put(Current_Output,K,2,10);
    Put(Current_Output,"  ");
  end ShowIndices;

  procedure ShowIndices(J: in Parity; M: in Freq1; N: in Freq2) is
  begin
    if J=0 then
      Put(Current_Output,"cos(");
    else
      Put(Current_Output,"sin(");
    end if;
    Put(Current_Output,M,3,10);
    Put(Current_Output,"x+");
    Put(Current_Output,N,2,10);
    Put(Current_Output,"y) E(z)   ");
  end ShowIndices;

  procedure ShowIndices(I: in FreqE; K: in Power) is
  begin
    Put(Current_Output,"E(dist >= ");
    Put(Current_Output,I,2,10);
    if K>0 then
      Put(Current_Output,") z**");
      Put(Current_Output,K,2,10);
      Put(Current_Output,"  ");
    else
      Put(Current_Output,")       ");
    end if;
  end ShowIndices;

  procedure ShowIndices(I: in FreqE) is
  begin
    Put(Current_Output,"E(dist >= ");
    Put(Current_Output,I,2,10);
    Put(Current_Output,") error  ");
  end ShowIndices;

  procedure Show(H: in Fourier; Cut: in Rep := Zero) is
    C: Scalar;
  begin
    Put(Current_Output,"J=");
    Put(Current_Output,H.J,1,10);
    Show(" R=",H.R);
    if not Poly(H,True) then
      Show("Poly(H) is False");
    end if;

    for M in Freq1 loop
      for N in Nmin(M) .. Ncut(M) loop
        for K in 0 .. Deg(M,N,H) loop
          C := Component(M,N,K,H);
          if SupAbs(C) > Cut then
            if K>Kcut(M,N) then
              Put(Current_Output,Kcut(M,N),2,10);
              Put(" ??? ");
            end if;
            ShowIndices(H.J,M,N,K);
            Show(C);
          end if;
        end loop;
        C := ErrComp(H.C(M,N));
        if SupAbs(C) > Cut then
          ShowIndices(H.J,M,N);
          Show(C);
        end if;
      end loop;
    end loop;
    New_Line;
    if Poly(H) then return; end if;

    for I in FreqE loop
      for K in 0 .. Pdeg loop
        C := Component(K,H.E(I));
        if SupAbs(C) > Cut then
          ShowIndices(I,K);
          Show(C);
        end if;
      end loop;
      C := ErrComp(H.E(I));
      if SupAbs(C) > Cut then
        ShowIndices(I);
        Show(C);
      end if;
     end loop;
    New_Line;
  end Show;

  procedure Read(Name: in String; J: in Parity; R: in Weights; H: out Fourier; Decimal: in Boolean := True) is
    use Reps.IO, Components_IO;
    F: File_Type;
    GetPoly: Boolean;
    Warn,Truncated: Boolean := False;
    S,Sz,SzE,Dh,KcutMN: Integer;
    Rt: Weights;
    P: Taylor;
  begin
    Put_Line("Reading " & Name);
    SetZero(J,R,H);

    Open(F,In_File,Name);
    Get(F,S);
    UnJoin(S,SzE,Dh,Sz);
    GetPoly := S<10000;
    if not GetPoly then
      Truncated := Dh<Dho;
      H.P := False;
      Txt_Get(F,Rt.Q,Decimal);
      Txt_Get(F,Rt.P,Decimal);
      H.R := Min(Rt,R);
    end if;

    for M in -Sz .. Sz loop
      for N in Nmin(M) .. Sz loop
        Get(R.P,F,P,H.P,GetPoly,Decimal);
        KcutMN := Kcut(M,N);
        if KcutMN >= 0 then
          Copy(P,H.C(M,N));
          Warn := Warn or Deg(P)>KcutMN;
        elsif Deg(P) >= 0 then
          Truncated := True;
        end if;
      end loop;
    end loop;

    if not GetPoly then
      for I in 0 .. SzE loop
        Get(R.P,F,P,H.P,GetPoly,Decimal);
        if I <= SizeE then
          Copy(P,H.E(I));
        elsif Deg(P) >= 0 then
          Truncated := True;
        end if;
      end loop;
    end if;
    Close(F);

    H.P := H.P or Trunc;
    if H.P then
      H.R := R;
    elsif H.R /= R then
      Show("Reduced domain to ",H.R);
    end if;
    if Truncated then
      Put_Line("Truncated input");
      if not Trunc then
        Message("Fouriers.IO.Read",Index_Error);
      end if;
    elsif Warn then
      Put_Line("Some degrees exceed Dmax");
    end if;
  end Read;

  procedure Write(Name: in String; H: in Fourier; Decimal: in Boolean := True) is
    use Reps.IO, Components_IO;
    F: File_Type;
    PutPoly: constant Boolean := Poly(H,True);
    Dm,Dn,Sz,SzE: Integer;
    P0: Taylor;
  begin
    Put_Line("Writing H to " & Name);
    Show("H.R=",H.R);

    EffDeg(H,Dm,Dn);
    Sz := Imax(Dm,Dn);
    SzE := EffSizeE(H);
    SetZero(P0);
    Create(F,Out_File,Name);

    if PutPoly then
      Put(F,Sz,3);
    else
      Put(F,Join(SzE,Dho,Imax(Sz,0)),7);
      New_Line(F);
      Txt_Put(F,H.R.Q,Decimal);
      New_Line(F);
      Txt_Put(F,H.R.P,Decimal);
    end if;
    New_Line(F);

    for M in -Sz .. Sz loop
      for N in Nmin(M) .. Sz loop
        if Kcut(M,N)>=0 then
          Put(F,H.C(M,N),PutPoly,Decimal);
        else
          Put(F,P0,PutPoly,Decimal);
        end if;
      end loop;
    end loop;
    if Sz<0 then
      Put(F,P0,PutPoly,Decimal);
    end if;
    if not PutPoly then
      for I in 0 .. SzE loop
        Put(F,H.E(I),PutPoly,Decimal);
      end loop;
    end if;

    Close(F);
  end Write;

  procedure EffDeg(H: in Sobolev; Dm,Dn: out Integer) is
  begin
    Dm := -1;
    Dn := -1;
    for M in Freq1 loop
      for N in Nmin(M) .. Ncut(M) loop
        if H.C(M,N) /= ScalZero then
          Dm := Imax(Abs(M),Dm);
          Dn := Imax(N,Dn);
        end if;
      end loop;
    end loop;
  end EffDeg;

  procedure Show(H: in Sobolev; Cut: in Rep := Zero) is
    C: Scalar;
  begin
    Put(Current_Output,"J=");
    Put(Current_Output,H.J,1,10);
    Show(" R=",H.R);

    for M in Freq1 loop
      for N in Nmin(M) .. Ncut(M) loop
        C := Component(M,N,H);
        if SupAbs(C) > Cut then
          ShowIndices(H.J,M,N,0);
          Show(C);
        end if;
      end loop;
    end loop;
    New_Line;
    if Trunc then return; end if;

    for I in FreqE loop
      C := H.E(I);
      if SupAbs(C) > Cut then
        ShowIndices(I);
        Show(C);
      end if;
     end loop;
    New_Line;
  end Show;

  procedure ShowNorms(H: in Sobolev; Cut: in Rep := Zero) is
    C,E: Scalar := ScalZero;
  begin
    Put(Current_Output,"J=");
    Put(Current_Output,H.J,1,10);
    Show(" R=",H.R);

    for M in Freq1 loop
      for N in Nmin(M) .. Ncut(M) loop
        C := Exp(H.R*Log(ScalOne+Abs(Scal(M)+Rep(-N)*InvGold)));
        C := Abs(Component(M,N,H)*C);
        E := E+C;
        if Sup(C) > Cut then
          ShowIndices(H.J,M,N,0);
          Show(" Norm ",Sup(C));
        end if;
      end loop;
    end loop;
    New_Line;

    if not Trunc then
      for I in FreqE loop
        C := Abs(H.E(I));
        E := E+C;
        if Sup(C) > Cut then
          ShowIndices(I);
          Show(" Norm ",Sup(C));
        end if;
      end loop;
      New_Line;
    end if;
    New_Line;
    Show(" Norm total ",Sup(E));
    New_Line;
  end ShowNorms;

  procedure Read(Name: in String; J: in Parity; R: in Radius; H: out Sobolev; Decimal: in Boolean := True) is
    use Reps.IO, Components_IO;
    F: File_Type;
    GetPoly: Boolean;
    Truncated: Boolean := False;
    S,Sz,SzE,Dh: Integer;
    Rt: Radius;
    C: Scalar;
  begin
    Put_Line("Reading " & Name);
    SetZero(J,R,H);

    Open(F,In_File,Name);
    Get(F,S);
    UnJoin(S,SzE,Dh,Sz);
    GetPoly := S<10000;
    if not GetPoly then
      Txt_Get(F,Rt,Decimal);
      H.R := Radius'Min(Rt,R);
    end if;

    for M in -Sz .. Sz loop
      for N in Nmin(M) .. Sz loop
        Get(F,C,Decimal);
        if Kcut(M,N) >= 0 then
          H.C(M,N) := C;
        elsif C /= ScalZero then
          Truncated := True;
        end if;
      end loop;
    end loop;

    if not GetPoly then
      for I in 0 .. SzE loop
        Get(F,C,Decimal);
        if I <= SizeE then
          H.E(I) := C;
        elsif C /= ScalZero then
          Truncated := True;
        end if;
      end loop;
    end if;
    Close(F);

    if GetPoly or Trunc then
      H.R := R;
    elsif Rt<R then
      H.R := Rt;
      Show("Reduced domain to ",Rt);
    end if;
    if Truncated then
      Put_Line("Truncated input");
      if not Trunc then
        Message("Fouriers.IO.Read",Index_Error);
      end if;
    end if;
  end Read;

  procedure Write(Name: in String; H: in Sobolev; Decimal: in Boolean := True) is
    use Reps.IO, Components_IO;
    PutPoly: Boolean;
    F: File_Type;
    Dm,Dn,Sz,SzE: Integer;
  begin
    Put_Line("Writing H to " & Name);
    Show("H.R=",H.R);

    EffDeg(H,Dm,Dn);
    Sz := Imax(Dm,Dn);
    SzE := EffSizeE(H);
    PutPoly := Trunc or else SzE<0;
    Create(F,Out_File,Name);

    if PutPoly then
      Put(F,Sz,3);
    else
      Put(F,Join(SzE,0,Sz),7);
      New_Line(F);
      Txt_Put(F,H.R,Decimal);
    end if;
    New_Line(F);

    for M in -Sz .. Sz loop
      for N in Nmin(M) .. Sz loop
        if Kcut(M,N)>=0 then
          Put(F,H.C(M,N),Decimal);
          New_Line(F);
        else
          Put(F,ScalZero,Decimal);
          New_Line(F);
        end if;
      end loop;
    end loop;
    if not PutPoly then
      for I in 0 .. SzE loop
        Put(F,H.E(I),Decimal);
        New_Line(F);
      end loop;
    end if;

    Close(F);
  end Write;

end Fouriers.IO;