with Chebychev_Polynomials;              use Chebychev_Polynomials;
with Standard_Natural_Vectors;
with Standard_Floating_Vectors;          use Standard_Floating_Vectors;
with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Standard_Floating_Linear_Solvers;   use Standard_Floating_Linear_Solvers;
with Standard_Floating_QR_Decomposition; use Standard_Floating_QR_Decomposition;

package body Osculating_Planes is

  function Standard_Basis
                ( n,d : natural; s : double_float ) return Matrix is

    res : Matrix(1..n,1..d);
    acc : double_float;
    fac,lim : integer;

  begin
    for i in 1..d loop                     -- set the zeros and ones
      res(i,i) := 1.0;
      for j in (i+1)..d loop
        res(i,j) := 0.0;
      end loop;
    end loop;
    for j in 1..d loop                     -- set the powers of the s-values
      acc := s;
      for i in (j+1)..n loop
        res(i,j) := acc;
        acc := acc*s;
      end loop;
    end loop;
    for i in 3..n loop                  -- compute the factors from derivation
      fac := 1;
      if i-1 > d
       then lim := d;
       else lim := i-1;
      end if;
      for j in 2..lim loop
        fac := fac*(i+1-j);
        res(i,j) := double_float(fac)*res(i,j);
      end loop;
      if i <= d
       then res(i,i) := double_float(fac);
      end if;
    end loop;
    for j in 3..d loop                     -- scale the columns
      for i in (j+1)..n loop
        res(i,j) := res(i,j)/res(j,j);
      end loop;
      res(j,j) := 1.0;
    end loop;
    return res;
  end Standard_Basis;

  function Chebychev_Basis
                ( n,d : natural; s : double_float ) return Matrix is

    res : Matrix(1..n,1..d);
    lim : integer;

  begin
    for i in 1..d loop                     -- set the zeros and ones
      res(i,i) := 1.0;
      for j in (i+1)..d loop
        res(i,j) := 0.0;
      end loop;
    end loop;
    for i in 2..n loop                     -- evaluate and differentiate
      declare
        p : constant Vector := Chebychev_Polynomials.Create(i-1);
      begin
        res(i,1) := Eval(p,s);
        if i > d
         then lim := d;
         else lim := i;
        end if;
        for j in 2..lim loop
          declare
            dp : constant Vector := Chebychev_Polynomials.Diff(p,j-1);
          begin
            res(i,j) := Eval(dp,s);
          end;
        end loop;
      end;
    end loop;
    for j in 3..d loop                     -- scale the columns
      for i in (j+1)..n loop
        res(i,j) := res(i,j)/res(j,j);
      end loop;
      res(j,j) := 1.0;
    end loop;
    return res;
  end Chebychev_Basis;

  function Orthogonal_Basis
              ( n,d : natural; s : double_float ) return Matrix is

    res : Matrix(1..n,1..d) := Chebychev_Basis(n,d,s);
    wrk : Matrix(1..n,1..d);
    bas : Matrix(1..n,1..n);
    qra : Standard_Floating_Vectors.Vector(1..n) := (1..n => 0.0);
    pvt : Standard_Natural_Vectors.Vector(1..n) := (1..n => 0); 

  begin
    wrk := res;
    QRD(wrk,qra,pvt,false);
    for i in wrk'range(1) loop
      for j in wrk'range(2) loop
        bas(i,j) := wrk(i,j);
      end loop;
      for j in d+1..n loop
        bas(i,j) := 0.0;
	  end loop;
    end loop;
    Basis(bas,res); 
    for i in res'range(1) loop
      for j in res'range(2) loop
        res(i,j) := bas(i,j);
      end loop;
    end loop;
    return res;
  end Orthogonal_Basis;

  function Determinant
              ( mat : Matrix; rows : Standard_Natural_Vectors.Vector )
                return double_float is


  -- DESCRIPTION :
  --   Computes the determinant of the matrix obtained by selecting rows.

    res : double_float := 1.0;
    sqm : Matrix(rows'range,rows'range);
    piv : Standard_Natural_Vectors.Vector(rows'range);
    inf : natural;

  begin
    for i in rows'range loop
      piv(i) := i;
      for j in rows'range loop
        sqm(i,j) := mat(rows(i),j);
      end loop;
    end loop;
    lufac(sqm,rows'last,piv,inf);
    for i in rows'range loop
      res := res*sqm(i,i);
    end loop;
    for i in piv'range loop
      if piv(i) > i
       then res := -res;
      end if;
    end loop;
    return res;
  end Determinant;

  procedure Maximal_Minors ( n,d : in natural; mat : in Matrix;
                             min,max : out double_float ) is

  -- DESCRIPTION :
  --   Computes all maximal minors of a (nxd)-matrix mat, d < n.

    rows : Standard_Natural_Vectors.Vector(1..d);
    first : boolean := true;
    mindet,maxdet : double_float;

    procedure Select_Rows ( k,start : in natural ) is

      det : double_float;

    begin
      if k > d 
       then det := Determinant(mat,rows);
            det := abs(det);
            if first
             then mindet := det; maxdet := det; first := false;
             else if det > maxdet
                   then maxdet := det;
                   elsif det < mindet
                       then mindet := det;
                  end if;
            end if;
       else for j in start..n loop
              rows(k) := j;
              Select_Rows(k+1,j+1);
            end loop;
      end if;
    end Select_Rows;

  begin
    Select_Rows(1,1);
    min := mindet; max := maxdet;
  end Maximal_Minors;

  procedure Sampled_Chebychev_Basis
                  ( n,d,m : in natural; mat : out Matrix;
                    s,ratio : out double_float ) is                   

    rs,min,max,rat,bestrat,bests : double_float;
    first : boolean;
    themat : Matrix(1..n,1..d);

  begin
    first := true;
    for i in 1..m loop
      rs := Random;
      themat := Chebychev_Basis(n,d,rs); 
      Maximal_Minors(n,d,themat,min,max);
      rat := max/min; 
      if first
       then bestrat := rat; bests := rs; first := false;
       else if rat < bestrat
             then bestrat := rat; bests := rs;
                  mat := themat;
            end if;
      end if;
    end loop;
    ratio := bestrat;
    s := bests;
  end Sampled_Chebychev_Basis;

end Osculating_Planes;