Return to floating_pruning_methods.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Stalift

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with Standard_Natural_Vectors;
with Standard_Floating_Matrices;         use Standard_Floating_Matrices;
with Standard_Floating_Linear_Solvers;   use Standard_Floating_Linear_Solvers;
with Floating_Linear_Inequalities;       use Floating_Linear_Inequalities;
with Lists_of_Floating_Vectors;          use Lists_of_Floating_Vectors;

package body Floating_Pruning_Methods is

-- GENERAL AUXILIARIES :

  procedure Normalize ( v : in out Standard_Floating_Vectors.Vector ) is

  -- DESCRIPTION : Divides every entry by v(v'last).

  begin
    for i in v'range loop
      v(i) := v(i)/v(v'last);
    end loop;
  end Normalize;

  function Convert ( fa : Face_Array ) return Array_of_Lists is

  -- DESCRIPTION :
  --   Converts the array of faces into an array of lists by
  --   converting the first element of each list of faces.

    res : Array_of_Lists(fa'range);

  begin
    for k in fa'range loop
      res(k) := Shallow_Create(fa(k).all);
    end loop;
    return res;
  end Convert;

-- AUXILIARIES FOR THE PRUNING ALGORITHMS :

  procedure Initialize ( n : in natural; mat : out Matrix; 
                         ipvt : out Standard_Natural_Vectors.Vector ) is

  -- DESCRIPTION :
  --   Initializes an n*(n+1) matrix with zeroes and the pivoting vector
  --   with the entries 1..n.

  begin
    for i in 1..n loop
      for j in 1..n+1 loop
        mat(i,j) := 0.0;
      end loop;
    end loop;
    for i in 1..n loop
      ipvt(i) := i;
    end loop;
  end Initialize;

  function Number_of_Inequalities
              ( mix : Vector; lifted : Array_of_Lists ) return natural is

  -- DESCRIPTION :
  --   Returns the maximal number of inequalities for pruning.

    res : natural := 0;

  begin
    for k in lifted'range loop
      res := res + Length_Of(lifted(k)) - mix(k) - 1;
    end loop;
    return res;
  end Number_of_Inequalities;

  procedure Ordered_Inequalities ( k : in natural; mat : in out Matrix ) is

  -- DESCRIPTION :
  --   Defines k inequalities mat(k,k) - mat(k+1,k) >= 0.

  begin
    for i in mat'first(1)..k loop
      for j in mat'range(2) loop
        mat(i,j) := 0.0;
      end loop;
      mat(i,i) := 1.0; mat(i,i+1) := -1.0;
    end loop;
  end Ordered_Inequalities;

  procedure Check_and_Update
                ( mic : in Face_Array; lifted : in Array_of_Lists;
                  m : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
                  tol : in double_float;
                  mixsub,mixsub_last : in out Mixed_Subdivision ) is

  -- DESCRIPTION :
  --   Computes the normal to the points in pts, by solving the 
  --   linear system defined by m and ipvt.
  --   If the computed normal is an inner normal w.r.t. the lifted points,
  --   then the mixed subdivision will be updated with a new cell.

    use Standard_Floating_Vectors;
    v : Standard_Floating_Vectors.Vector(m'range(2)) := Solve(m,tol,ipvt);
    pts : Array_of_Lists(mic'range);

  begin
    if abs(v(v'last)) > tol
     then Normalize(v);
          if v(v'last) < 0.0
           then Min(v);
          end if;
         -- pts := Convert(mic);
         -- Update(pts,v,mixsub,mixsub_last);
          declare
            cell : Mixed_Cell;
          begin
            cell.nor := new Standard_Floating_Vectors.Vector'(v);
            cell.pts := new Array_of_Lists'(Convert(mic));
            cell.sub := null;
            Append(mixsub,mixsub_last,cell);
          end;
    end if;
  end Check_and_Update;

  procedure Create_Equalities
                ( n : in natural; f : in Face; mat,ineq : in Matrix;
                  tol : in double_float;
                  ipvt : in Standard_Natural_Vectors.Vector;
                  row,rowineq : in natural;
                  newmat,newineq : in out Matrix;
                  newipvt : in out Standard_Natural_Vectors.Vector;
                  newrow,newrowineq : in out natural; fail : out boolean ) is

  -- DESCRIPTION :
  --   Creates new equalities and uses them to eliminate unknowns in the
  --   matrices of equalities and inequalities.  Failure is reported when
  --   a zero row is encountered.  On entry, all new* variables must be
  --   initialized with the corresponding *-ones.

    shi : Standard_Floating_Vectors.Vector(1..n+1) := f(f'first).all;
    fl : boolean := false;
    pivot : natural;

  begin
    for i in f'first+1..f'last loop
      newrow := newrow + 1;
      for j in f(i)'range loop
        newmat(newrow,j) := f(i)(j) - shi(j);
      end loop;
      Switch(newipvt,newrow,newmat);
      Upper_Triangulate(newrow,newmat,tol,newipvt,pivot);
      fl := (pivot = 0);
      exit when fl;
      Switch(newrow,pivot,ineq'first,rowineq,newineq);
    end loop;
    fail := fl;
  end Create_Equalities;

  procedure Complementary_Slackness
                  ( tableau : in out matrix;
                    tol : in double_float; feasible : out boolean ) is

    lastcol : constant integer := tableau'last(2)-1;
    rhs,sol : Standard_Floating_Vectors.Vector(tableau'range(1));
    columns : Vector(sol'range);

  begin
    for i in rhs'range loop
      rhs(i) := double_float(tableau(i,tableau'last(2)));
    end loop;
    Complementary_Slackness(tableau,lastcol,rhs,tol,sol,columns,feasible);
  end Complementary_Slackness;

  function Check_Feasibility ( k,m,n : natural; ineq : Matrix;
                               tol : double_float ) return boolean is
 
  -- DESCRIPTION :
  --   Returns true if -v(n+1) > 0 can be derived, false otherwise.

  -- ON ENTRY :
  --   k      current unknown that has been eliminated,
  --           for all i <= k : ineq(l,i) = 0, for l in ineq'first..m;
  --   m      number of inequalities;
  --   n      dimension;
  --   ineq   matrix of inequalities.

    tableau : Matrix(1..n-k+1,1..m+1);
    feasi : boolean;

  begin
    if m = 0
     then feasi := false;
     else for i in k+1..n+1 loop
            for j in 1..m loop
              tableau(i-k,j) := ineq(j,i);
            end loop;
            tableau(i-k,m+1) := 0.0;
          end loop;
          tableau(n-k+1,m+1) := -1.0;
          Complementary_Slackness(tableau,tol,feasi);
    end if;
    return feasi;
  end Check_Feasibility;

  procedure Update_Inequalities
               ( k,rowmat1,rowmat2,n : in natural;
                 mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
                 tol : in double_float;
                 rowineq : in out natural; ineq : in out Matrix;
                 lifted : in Array_of_Lists; mic : in out Face_Array ) is

  -- DESCRIPTION :
  --   The inequalities will be updated w.r.t. the equality
  --   constraints on the inner normal.

    tmp : List;
    pt,shi : Standard_Floating_Vectors.Link_to_Vector;

  begin
    for i in ineq'first..rowineq loop      -- update the old inequalities
      for j in rowmat1..rowmat2 loop
        Upper_Triangulate(j,mat,tol,i,ineq);
      end loop;
    end loop;
    shi := mic(k)(mic(k)'first);
    tmp := lifted(k);                      -- make new inequalities
    while not Is_Null(tmp) loop
      pt := Head_Of(tmp);
      if not Is_In(mic(k),pt.all)
       then rowineq := rowineq + 1;
            for j in pt'range loop
              ineq(rowineq,j) := pt(j) - shi(j);
            end loop;
            Switch(ipvt,rowineq,ineq);
            for i in 1..rowmat2 loop
              Upper_Triangulate(i,mat,tol,rowineq,ineq);
            end loop;
      end if;
      tmp := Tail_Of(tmp);
    end loop;
  end Update_Inequalities;

-- CONSTRUCTION WITH PRUNING :

  procedure Gen1_Create
               ( n : in natural; mix : in Vector; fa : in Array_of_Faces; 
                 lifted : in Array_of_Lists; tol : in double_float;
                 nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
                 mixsub : out Mixed_Subdivision ) is

    res,res_last : Mixed_Subdivision;
    accu : Face_Array(fa'range);
    ma : Matrix(1..n,1..n+1);
    ipvt : Standard_Natural_Vectors.Vector(1..n);
    ineqrows : natural;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in Matrix;
                   ipvt : in Standard_Natural_Vectors.Vector;
                   rowineq : in natural; ineq : in Matrix;
                   mic : in out Face_Array; continue : out boolean );

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

    -- ON ENTRY :
    --   k         index for current point set;
    --   row       number of rows already in the matrix;
    --   mat       matrix which determines the inner normal;
    --   ipvt      contains the pivoting information;
    --   rowineq   number of inequality constraints already in ineq;
    --   ineq      matrix for the inequality constraints on the
    --             inner normal v, of type <.,v> >= 0;
    --   mic       contains the current selected faces, up to k-1.

    -- ON RETURN :
    --   mic       updated selected faces.
    --   continue  indicates whether to continue the creation or not.

    procedure Process_Inequalities
                 ( k,rowmat1,rowmat2 : in natural;
                   mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
                   rowineq : in out natural; ineq : in out Matrix;
                   mic : in out Face_Array; cont : out boolean ) is

    -- DESCRIPTION :
    --   Updates inequalities and checks feasibility before proceeding.

      fl : boolean := false;
      tmp : List;
      pt,shi : Link_to_Vector;

    begin
      Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,tol,rowineq,ineq,
                          lifted,mic);
      if Check_Feasibility(rowmat2,rowineq,n,ineq,tol)
       then nbfail(k) := nbfail(k) + 1.0;
            cont := true;
       else nbsucc(k) := nbsucc(k) + 1.0;
            Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic,cont);
      end if;
    end Process_Inequalities;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in Matrix;
                   ipvt : in Standard_Natural_Vectors.Vector;
                   rowineq : in natural; ineq : in Matrix;
                   mic : in out Face_Array; continue : out boolean ) is

      old : Mixed_Subdivision := res_last;
      cont : boolean := true;
      tmpfa : Faces;

    begin
      if k > mic'last
       then Check_and_Update(mic,lifted,mat,ipvt,tol,res,res_last);
            if old /= res_last
             then Process(Head_Of(res_last),continue);
             else continue := true;
            end if;
       else tmpfa := fa(k);
            while not Is_Null(tmpfa) loop  -- enumerate faces of kth polytope
              mic(k) := Head_Of(tmpfa);
              declare                                 -- update the matrices
                fl : boolean;
                newipvt : Standard_Natural_Vectors.Vector(ipvt'range) := ipvt;
                newmat : Matrix(mat'range(1),mat'range(2)) := mat;
                newineq : Matrix(ineq'range(1),ineq'range(2)) := ineq;
                newrow : natural := row;
                newrowineq : natural := rowineq;
              begin
                Create_Equalities
                   (n,mic(k),mat,ineq,tol,ipvt,row,rowineq,newmat,newineq,
                    newipvt,newrow,newrowineq,fl);
                if fl
                 then nbfail(k) := nbfail(k) + 1.0;
                 else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
                                           newrowineq,newineq,mic,cont);
                end if;
              end;
              tmpfa := Tail_Of(tmpfa);
              exit when not cont;
            end loop;
            continue := cont;
      end if;
    end Compute_Mixed_Cells;

  begin
    Initialize(n,ma,ipvt);
    ineqrows := Number_of_Inequalities(mix,lifted);
    declare
      ineq : matrix(1..ineqrows,1..n+1);
      cont : boolean;
    begin
      ineq(1,1) := 0.0;
      Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu,cont);
    end;
    mixsub := res;
  end Gen1_Create;

  procedure Create
              ( n : in natural; mix : in Vector; fa : in Array_of_Faces; 
                lifted : in Array_of_Lists; tol : in double_float;
                nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
		mixsub : out Mixed_Subdivision ) is

    res,res_last : Mixed_Subdivision;
    accu : Face_Array(fa'range);
    ma : Matrix(1..n,1..n+1);
    ipvt : Standard_Natural_Vectors.Vector(1..n);
    ineqrows : natural;

    procedure Compute_Mixed_Cells
                 ( k,row : in natural; mat : in Matrix;
                   ipvt : in Standard_Natural_Vectors.Vector;
                   rowineq : in natural; ineq : in Matrix;
                   mic : in out Face_Array );

    -- DESCRIPTION :
    --   Backtrack recursive procedure to enumerate face-face combinations.

    -- ON ENTRY :
    --   k         index for current point set;
    --   row       number of rows already in the matrix;
    --   mat       matrix which determines the inner normal;
    --   ipvt      contains the pivoting information;
    --   rowineq   number of inequality constraints already in ineq;
    --   ineq      matrix for the inequality constraints on the
    --             inner normal v, of type <.,v> >= 0;
    --   mic       contains the current selected faces, up to k-1.

    -- ON RETURN :
    --   mic       updated selected faces.

    procedure Process_Inequalities
                 ( k,rowmat1,rowmat2 : in natural;
                   mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
                   rowineq : in out natural; ineq : in out matrix;
                   mic : in out Face_Array) is

    -- DESCRIPTION :
    --   Updates inequalities and checks feasibility before proceeding.

      tmp : List;
      pt,shi : Link_to_Vector;

    begin
      Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,tol,rowineq,ineq,
                          lifted,mic);
      if Check_Feasibility(rowmat2,rowineq,n,ineq,tol)
       then nbfail(k) := nbfail(k) + 1.0;
       else nbsucc(k) := nbsucc(k) + 1.0;
            Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic);
      end if;
    end Process_Inequalities;

    procedure Compute_Mixed_Cells 
                 ( k,row : in natural; mat : in Matrix;
                   ipvt : in Standard_Natural_Vectors.Vector;
                   rowineq : in natural; ineq : in Matrix;
                   mic : in out Face_Array ) is

      tmpfa : Faces;

    begin
      if k > mic'last
       then Check_and_Update(mic,lifted,mat,ipvt,tol,res,res_last);
       else tmpfa := fa(k);
            while not Is_Null(tmpfa) loop  -- enumerate faces of kth polytope
              mic(k) := Head_Of(tmpfa);
              declare                                      -- update matrices
                fl : boolean;
                newipvt : Standard_Natural_Vectors.Vector(ipvt'range) := ipvt;
                newmat : Matrix(mat'range(1),mat'range(2)) := mat;
                newineq : Matrix(ineq'range(1),ineq'range(2)) := ineq;
                newrow : natural := row;
                newrowineq : natural := rowineq;
              begin
                Create_Equalities
                    (n,mic(k),mat,ineq,tol,ipvt,row,rowineq,newmat,newineq,
                     newipvt,newrow,newrowineq,fl);
                if fl 
                 then nbfail(k) := nbfail(k) + 1.0;
                 else Process_Inequalities
                        (k,row+1,newrow,newmat,newipvt,newrowineq,newineq,mic);
                end if;
              end;
              tmpfa := Tail_Of(tmpfa);
            end loop;
      end if;
    end Compute_Mixed_Cells;

  begin
    Initialize(n,ma,ipvt);
    ineqrows := Number_of_Inequalities(mix,lifted);
    declare
      ineq : Matrix(1..ineqrows,1..n+1);
    begin
      ineq(1,1) := 0.0;
      Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu);
    end;
    mixsub := res;
  end Create;

end Floating_Pruning_Methods;