Return to ts_fvector.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with text_io,integer_io;                 use text_io,integer_io;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Standard_Integer_Vectors;
with Standard_Floating_Vectors;
with Standard_Integer_VecVecs;
with Standard_Integer_VecVecs_io;        use Standard_Integer_VecVecs_io;
with Standard_Floating_VecVecs;
with Standard_Floating_VecVecs_io;       use Standard_Floating_VecVecs_io;
with Standard_Random_VecVecs;            use Standard_Random_VecVecs;
with Face_Cardinalities;                 use Face_Cardinalities;

procedure ts_fvector is

-- DESCRIPTION :
--   Computes the f-vector of a polytope and checks the Euler-Poincare formula.

  procedure Write ( f : in Standard_Integer_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Writes the f-vector on screen.

    n : constant natural := f'last;

  begin
    put_line("The f-vector : ");
    put(" f(-1)     : "); put(f(-1)); new_line;
    put(" #vertices : "); put(f(0)); new_line;
    put(" #edges    : "); put(f(1)); new_line;
    for i in 2..(n-1) loop
      put(" #"); put(i,1); put("-faces  : "); put(f(i)); new_line;
    end loop;
    put(" f("); put(n,1); put(")      : "); put(f(n)); new_line;
  end Write;

  function Euler_Poincare ( f : Standard_Integer_Vectors.Vector )
                          return integer is

  -- DESCRIPTION :
  --   Computes the alternating sum: sum_{i in f'range} (-1)^(i) f(i).

    sum : integer := 0;
    pos : boolean := false;

  begin
    for i in f'range loop
      if pos
       then sum := sum + f(i);
       else sum := sum - f(i);
      end if;
      pos := not pos;
    end loop;
    return sum;
  end Euler_Poincare;

  procedure Integer_Interactive_Testing is

    n,m : natural;
    tol : constant double_float := 10.0**(-12);
    sum : integer;

  begin
    put("Give the dimension n : "); get(n);
    put("Give the number of points that span the polytope : "); get(m);
    declare
      pts : Standard_Integer_VecVecs.VecVec(1..m);
      f : Standard_Integer_Vectors.Vector(-1..n);
    begin
      put("Give "); put(m,1); put(" "); put(n,1);
      put_line("-dimensional integer vectors :");
      get(n,pts);
      put_line("Counting vertices, edges , .., facets, ...");
      f := fvector(pts); sum := Euler_Poincare(f);
      Write(f); 
      put("The result of the Euler-Poincare formula : "); put(sum,1);
      if sum /= 0
       then put_line("   BUG DISCOVERED !!!");
       else put_line("   OK.");
      end if;
    end;
  end Integer_Interactive_Testing;

  procedure Floating_Interactive_Testing is

    n,m : natural;
    tol : constant double_float := 10.0**(-12);
    sum : integer;

  begin
    put("Give the dimension n : "); get(n);
    put("Give the number of points that span the polytope : "); get(m);
    declare
      pts : Standard_Floating_VecVecs.VecVec(1..m);
      f : Standard_Integer_Vectors.Vector(-1..n);
    begin
      put("Give "); put(m,1); put(" "); put(n,1);
      put_line("-dimensional floating point vectors :");
      get(n,pts);
      put_line("Counting vertices, edges , .., facets, ...");
      f := fvector(pts); sum := Euler_Poincare(f);
      Write(f);
      put("The result of the Euler-Poincare formula : "); put(sum,1);
      if sum /= 0
       then put_line("   BUG DISCOVERED !!!");
       else put_line("   OK.");
      end if;
    end;
  end Floating_Interactive_Testing;

  function Floating_Random_Polytope 
             ( n,m : natural ) return Standard_Floating_VecVecs.VecVec is

  -- DESCRIPTION :
  --   Returns m randomly chosen n-dimensional floating-point vectors.

    res : Standard_Floating_VecVecs.VecVec(1..m) := Random_VecVec(n,m);

  begin
    return res;
  end Floating_Random_Polytope;

  function Integer_Random_Polytope
             ( n,m : natural; lower,upper : integer )
             return Standard_Integer_VecVecs.VecVec is

  -- DESCRIPTION :
  --   Returns m randomly chosen n-dimensional vectors,
  --   with integer entries between lower and upper.

    res : Standard_Integer_VecVecs.VecVec(1..m);
    done : boolean;

    function Is_In ( i : integer ) return boolean is

    -- DESCRIPTION :
    --   Returns true if the ith vector already occurs in res(1..i-1).

      use Standard_Integer_Vectors;

    begin
      for j in 1..i-1 loop
        if Equal(res(j).all,res(i).all)
         then return true;
        end if;
      end loop;
      return false;
    end Is_In;

  begin
    for i in 1..m loop
      res(i) := new Standard_Integer_Vectors.Vector(1..n);
      done := false;
      while not done loop
        for j in 1..n loop
          res(i)(j) := Random(lower,upper);
        end loop;
        done := not Is_In(i);
      end loop;
    end loop;
    return res;
  end Integer_Random_Polytope;

  procedure Floating_Automatic_Testing is

    n,m,times,cnt : natural;
    tol : constant double_float := 10.0**(-12);
    sum : integer;
    bug : boolean := false;

  begin
    put("Give the dimension n : "); get(n);
    put("Give the number of points that span the polytope : "); get(m);
    put("Give the number of testing cycles : "); get(times);
    declare
      pts : Standard_Floating_VecVecs.VecVec(1..m);
      f : Standard_Integer_Vectors.Vector(-1..n);
    begin
      for i in 1..times loop
        cnt := i;
        pts := Floating_Random_Polytope(n,m);
        f := fvector(pts); sum := Euler_Poincare(f);
        Write(f);
        put("The result of the Euler-Poincare formula : "); put(sum,1);
        if sum /= 0
         then put_line("   BUG DISCOVERED !!!"); bug := true;
              put_line("The generated random configuration is");
              put(pts);
         else put_line("   OK.");  bug := false;
        end if;
        Standard_Floating_VecVecs.Clear(pts);
        exit when bug;
      end loop;
    end;
    if not bug
     then put("No bugs found, with "); put(times,1);
          put_line(" generated cases tested.");
          put("Dimension : "); put(n,1);
          put(" and #points : "); put(m,1); put_line(".");
     else put("Bug found at case "); put(cnt,1); put_line(".");
    end if;
  end Floating_Automatic_Testing;

  procedure Integer_Automatic_Testing is

    n,m,times,cnt : natural;
    tol : constant double_float := 10.0**(-12);
    sum : integer;
    bug : boolean := false;
    lower,upper : integer;

  begin
    put("Give the dimension n : "); get(n);
    put("Give the number of points that span the polytope : "); get(m);
    put("Give lower bound on the entries : "); get(lower);
    put("Give upper bound on the entries : "); get(upper);
    put("Give the number of testing cycles : "); get(times);
    declare
      pts : Standard_Integer_VecVecs.VecVec(1..m);
      f : Standard_Integer_Vectors.Vector(-1..n);
    begin
      for i in 1..times loop
        cnt := i;
        pts := Integer_Random_Polytope(n,m,lower,upper);
        f := fvector(pts); sum := Euler_Poincare(f);
        Write(f);
        put("The result of the Euler-Poincare formula : "); put(sum,1);
        if sum /= 0
         then put_line("   BUG DISCOVERED !!!"); bug := true;
              put_line("The generated random configuration is");
              put(pts);
         else put_line("   OK.");  bug := false;
        end if;
        Standard_Integer_VecVecs.Clear(pts);
        exit when bug;
      end loop;
    end;
    if not bug
     then put("No bugs found, with "); put(times,1);
          put_line(" generated cases tested.");
          put("Dimension : "); put(n,1);
          put(" and #points : "); put(m,1); put_line(".");
     else put("Bug found at case "); put(cnt,1); put_line(".");
    end if;
  end Integer_Automatic_Testing;

  procedure Main is

    ans : character;

  begin
    new_line;
    put_line("Testing the face enumerators by computing f-vectors.");
    loop
      new_line;
      put_line("Choose one of the following :                  ");
      put_line("  0. Exit this program.                        ");
      put_line("  1. f-vector of given integer polytope.       ");
      put_line("  2. f-vector of given floating polytope.      ");
      put_line("  3. f-vector of random integer polytope.      ");
      put_line("  4. f-vector of random floating polytope.     ");
      put("Type 0,1,2,3, or 4 to select : "); get(ans);
      exit when (ans = '0');
      case ans is
        when '1' => Integer_Interactive_Testing;
        when '2' => Floating_Interactive_Testing;
        when '3' => Integer_Automatic_Testing;
        when '4' => Floating_Automatic_Testing;
        when others => null;
      end case;
    end loop;
  end Main;

begin 
  Main;
end ts_fvector;