with integer_io,Numbers_io;              use integer_io,Numbers_io;
with Communications_with_User;           use Communications_with_User;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Standard_Floating_Vectors_io;       use Standard_Floating_Vectors_io;
with Standard_Floating_Matrices;
with Standard_Floating_Matrices_io;      use Standard_Floating_Matrices_io;
with Standard_Complex_Matrices;
with Standard_Random_Matrices;           use Standard_Random_Matrices;
with Standard_Complex_VecMats_io;        use Standard_Complex_VecMats_io;
with Osculating_Planes;                  use Osculating_Planes;

package body Drivers_for_Input_Planes is

  function Finite ( dim : Bracket; fin_sum : natural ) return boolean is

  -- DESCRIPTION :
  --   Returns true if the co-dimensions in dim will lead to a finite
  --   number of solutions.                         

    sum : natural := 0;

  begin
    for i in dim'range loop
      sum := sum + dim(i);
    end loop;
    if sum = fin_sum
     then return true;
     else return false;
    end if;
  end Finite;

  function Read_Codimensions ( m,p,q : natural ) return Bracket is

    mpq : constant natural := m*p + q*(m+p);
    codim : Bracket(1..mpq);
    n : natural;

  begin
    loop
      put("Give number of intersection conditions : "); get(n);
      put("Give "); put(n,1); put(" co-dimensions (sum = ");
      put(mpq,1); put(") : ");
      for i in 1..n loop
        get(codim(i));
      end loop;
      for i in 1..n-1 loop
        put(codim(i),1); put(" + ");
      end loop;
      put(codim(n),1);
      if Finite(codim(1..n),mpq)
	   then put(" = "); put(mpq,1); put_line("  Finite #sols.");
            exit;
	   else put(" /= "); put(mpq,1);
            put_line("  Please try again.");
      end if;
    end loop;
    return codim(1..n);
  end Read_Codimensions;

  function Convert ( mat : Standard_Floating_Matrices.Matrix )
                   return Standard_Complex_Matrices.Matrix is

  -- DESCRIPTION :
  --   Converts a real matrix into a complex one.

    res : Standard_Complex_Matrices.Matrix(mat'range(1),mat'range(2));

  begin
    for i in mat'range(1) loop
      for j in mat'range(2) loop
        res(i,j) := Create(mat(i,j));
      end loop;
    end loop;
    return res;
  end Convert;

  function Random_Complex_Planes ( m,p : natural ) return VecMat is

    res : VecMat(1..m*p);
    n : constant natural := m+p;
    use Standard_Complex_Matrices;

  begin
    for i in res'range loop
      res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
    end loop;
    return res;
  end Random_Complex_Planes;

  function Random_Complex_Planes ( m,p : natural; k : Bracket ) return VecMat is

    res : VecMat(k'range);
    n : constant natural := m+p;
    use Standard_Complex_Matrices;

  begin
    for i in res'range loop
      res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m+1-k(i)));
    end loop;
    return res;
  end Random_Complex_Planes;

  function Random_Complex_Planes ( m,p,q : natural ) return VecMat is

    res : VecMat(1..(m*p+q*(m+p)));
    n : constant natural := m+p;
    use Standard_Complex_Matrices;

  begin
    for i in res'range loop
      res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
    end loop;
    return res;
  end Random_Complex_Planes;      

  function Random_Real_Planes ( m,p : natural ) return VecMat is

    res : VecMat(1..m*p);
    n : constant natural := m+p;
    fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
    cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);

  begin
    for i in res'range loop
      fltmat := Random_Orthogonal_Matrix(n,m);
      cmpmat := Convert(fltmat);
      res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
    end loop;
    return res;
  end Random_Real_Planes;

  function Random_Real_Planes ( m,p : natural; k : Bracket ) return VecMat is

    res : VecMat(k'range);
    n : constant natural := m+p;

  begin
    for i in res'range loop
      declare
        fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
        cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m+1-k(i));
      begin
        fltmat := Random_Orthogonal_Matrix(n,m+1-k(i));
        cmpmat := Convert(fltmat);
        res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
      end;
    end loop;
    return res;
  end Random_Real_Planes;

  function Random_Real_Planes ( m,p,q : natural ) return VecMat is

    res : VecMat(1..(m*p+q*(m+p)));
    n : constant natural := m+p;
    fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
    cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);

  begin
    for i in res'range loop
      fltmat := Random_Orthogonal_Matrix(n,m);
      cmpmat := Convert(fltmat);
      res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
    end loop;
    return res;
  end Random_Real_Planes;      

  function Equidistant_Interpolation_Points ( n : natural ) return Vector is

    res : Vector(1..n);
    s : double_float := Random;
    inc : constant double_float := 2.0/double_float(n);

  begin
    for i in res'range loop
      res(i) := s;
      s := s+inc;
      if s >= 1.0
       then s := s - 2.0;
      end if;
    end loop;
    return res;
  end Equidistant_Interpolation_Points;

  function Read_Interpolation_Points ( n : natural ) return Vector is

    res : Vector(1..n);

  begin
    new_line;
    put("Give "); put(n,1); put_line(" distinct real interpolation points : ");
    for i in res'range loop
      Read_Double_Float(res(i));
    end loop;       
    return res;
  end Read_Interpolation_Points;

  function Osculating_Input_Planes ( m,p : natural ) return VecMat is

    s : constant Vector := Equidistant_Interpolation_Points(m*p);

  begin
    return Osculating_Input_Planes(m,p,s);
  end Osculating_Input_Planes;

  function Osculating_Input_Planes
             ( m,p : natural; s : Vector ) return VecMat is

    res : VecMat(1..m*p);
    n : constant natural := m+p;
    realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);

  begin
    for i in res'range loop
      realmat := Orthogonal_Basis(n,m,s(i));
      res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
    end loop;
    return res;
  end Osculating_Input_Planes;

  function Osculating_Input_Planes
             ( m,p : natural; k : bracket ) return VecMat is

    s : constant Vector := Equidistant_Interpolation_Points(k'length);

  begin
    return Osculating_Input_Planes(m,p,k,s);
  end Osculating_Input_Planes;

  function Osculating_Input_Planes
             ( m,p : natural; k : bracket; s : Vector ) return VecMat is

    res : VecMat(k'range);
    n : constant natural := m+p;

  begin
    for i in res'range loop
      declare
        realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
      begin
        realmat := Orthogonal_Basis(n,m+1-k(i),s(i));
        res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
      end;
    end loop;
    return res;
  end Osculating_Input_Planes;

  function Osculating_Input_Planes ( m,p,q : natural ) return VecMat is

    dim : constant natural := m*p + q*(m+p);
    s : constant Vector := Equidistant_Interpolation_Points(dim);

  begin
    return Osculating_Input_Planes(m,p,q,s);
  end Osculating_Input_Planes;

  function Osculating_Input_Planes
             ( m,p,q : natural; s : Vector ) return VecMat is

    res : VecMat(1..m*p+q*(m+p));
    n : constant natural := m+p;
    realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);

  begin
    for i in res'range loop
      realmat := Orthogonal_Basis(n,m,s(i));
      res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
    end loop;
    return res;
  end Osculating_Input_Planes;

  function Read_Input_Planes ( m,p : natural ) return VecMat is

    res : VecMat(1..m*p);
    planesfile : file_type;
    n : constant natural := m+p;
    realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);

  begin
    new_line;
    put_line("Reading the name of the file with the input planes.");
    Read_Name_and_Open_File(planesfile);
    for i in res'range loop
      get(planesfile,realmat);
      realmat := Orthogonalize(realmat);
      res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
    end loop;
    Close(planesfile);
    return res;
  end Read_Input_Planes;

  function Read_Input_Planes ( m,p : natural; k : in Bracket ) return VecMat is

    res : VecMat(k'range);
    planesfile : file_type;
    n : constant natural := m+p;

  begin
    new_line;
    put_line("Reading the name of the file with the input planes.");
    Read_Name_and_Open_File(planesfile);
    for i in res'range loop
      declare
        realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
      begin
        get(planesfile,realmat);
        realmat := Orthogonalize(realmat);
        res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
      end;
    end loop;
    Close(planesfile);
    return res;
  end Read_Input_Planes;

  function Read_Input_Planes ( m,p,q : natural ) return VecMat is

    res : VecMat(1..m*p+q*(m+p));
    planesfile : file_type;
    n : constant natural := m+p;
    realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);

  begin
    new_line;
    put_line("Reading the name of the file with the input planes.");
    Read_Name_and_Open_File(planesfile);
    for i in res'range loop
      get(planesfile,realmat);
      realmat := Orthogonalize(realmat);
      res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
    end loop;
    Close(planesfile);
    return res;
  end Read_Input_Planes;

-- MAIN INTERACTIVE DRIVERS :

  function Select_Input_Choice return character is

  -- DESCRIPTION :
  --   Displays the menu to obtain a choice for the kind of input.

    res : character;

  begin
    new_line;
    put_line("MENU for generating real input planes : ");
    put_line("  0. Random real planes at equidistant interpolation points.");
    put_line("  1. Generate input planes osculating a rational normal curve.");
    put_line("  2. Interactively give s-values for the "
                               & "osculating input planes.");
    put_line("  3. Give the name of the file with input planes.");
    put("Type 0, 1, 2, or 3 to select : "); Ask_Alternative(res,"0123");
    return res;
  end Select_Input_Choice;

  procedure Driver_for_Input_Planes
              ( file : in file_type; m,p : in natural; planes : out VecMat ) is

    input_choice : character := Select_Input_Choice;
    dim : constant natural := m*p;
    svals : Vector(1..dim);

  begin
    case input_choice is
      when '0' => planes := Random_Real_Planes(m,p);
      when '1' => planes := Osculating_Input_Planes(m,p);
      when '2' => svals := Read_Interpolation_Points(dim);
                  planes := Osculating_Input_Planes(m,p,svals);
      when '3' => planes := Read_Input_Planes(m,p);
      when others => null;
    end case;
  end Driver_for_Input_Planes;

  procedure Driver_for_Input_Planes
              ( file : in file_type; m,p : in natural; k : in Bracket;
                planes : out VecMat ) is

    input_choice : character := Select_Input_Choice;
    svals : Vector(k'range);

  begin
    case input_choice is
      when '0' => planes := Random_Real_Planes(m,p,k);
      when '1' => planes := Osculating_Input_Planes(m,p,k);
      when '2' => svals := Read_Interpolation_Points(k'length);
                  planes := Osculating_Input_Planes(m,p,k,svals);
      when '3' => planes := Read_Input_Planes(m,p,k);
      when others => null;
    end case;
  end Driver_for_Input_Planes;

  procedure Driver_for_Input_Planes
              ( file : in file_type; m,p,q : in natural;
                s : out Vector; planes : out VecMat ) is

    input_choice : character := Select_Input_Choice;
    dim : constant natural := m*p + q*(m+p);
    svals : Vector(1..dim);

  begin
    case input_choice is
      when '0' => svals := Equidistant_Interpolation_Points(dim);
                  s := svals;
                  planes := Random_Real_Planes(m,p,q);
      when '1' => svals := Equidistant_Interpolation_Points(dim);
                  s := svals;
                  planes := Osculating_Input_Planes(m,p,q,svals);
      when '2' => svals := Read_Interpolation_Points(dim);
                  s := svals;
                  planes := Osculating_Input_Planes(m,p,q,svals);
      when '3' => svals := Read_Interpolation_Points(dim);
                  s := svals;
                  planes := Read_Input_Planes(m,p,q);
      when others => null;
    end case;
    put_line(file,"The interpolation points : ");
    put_line(file,svals);
    put_line(file,"The target planes : ");
    put(file,planes,3);
  end Driver_for_Input_Planes;

end Drivers_for_Input_Planes;