Return to inner_normal_cones.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_Integer_Linear_Solvers;    use Standard_Integer_Linear_Solvers;
with Transformations;                    use Transformations;
with Integer_Vectors_Utilities;          use Integer_Vectors_Utilities;
with Lists_of_Vectors_Utilities;         use Lists_of_Vectors_Utilities;
with Transforming_Integer_Vector_Lists;  use Transforming_Integer_Vector_Lists;
with Vertices;                           use Vertices;

package body Inner_Normal_Cones is

-- NOTE ON THE IMPLEMENTATION :
--   There is a safety mode in which after each computation of the generators,
--   it is automatically tested whether the generators satisfy the 
--   inequalities of the normal cone.  If not, then the bug is reported and
--   a program_error is raised.

-- AUXILIAIRIES :

  function Lower ( a1,b1,a2,b2 : integer ) return boolean is

  -- DESCRIPTION :
  --   Returns true if a1/b1 < a2/b2, false otherwise. 

  -- REQUIRED : b1 /= 0 and b2 /= 0.

  begin
    if b1*b2 > 0
     then return (a1*b2 - a2*b1 < 0);
     else return (a1*b2 - a2*b1 > 0);
    end if;
  end Lower;

  function Higher ( a1,b1,a2,b2 : integer ) return boolean is

  -- DESCRIPTION :
  --   Returns true if a1/b1 > a2/b2, false otherwise. 

  -- REQUIRED : b1 /= 0 and b2 /= 0.

  begin
    if b1*b2 > 0
     then return (a1*b2 - a2*b1 > 0);
     else return (a1*b2 - a2*b1 < 0);
    end if;
  end Higher;

  function Compute_Facets ( l : List; x : Vector ) return Faces is

  -- DESCRIPTION :
  --   Returns all facets of conv(l) that all contain x.
  --   Checks whether x belongs to the list or not.

    res : Faces;
    n : constant natural := x'length;

  begin
    if Is_In(l,x)
     then res := Create(n-1,n,l,x);
     else declare
            wrk : List := l;
            lx : Link_to_Vector;
          begin
            lx := new Vector'(x);
            Construct(lx,wrk);
            res := Create(n-1,n,wrk,x);
          end;
    end if;
    return res;
  end Compute_Facets;

  procedure Shifted_Points ( l : in List; x : in Vector;
                             res,res_last : in out List ) is

  -- DESCRIPTION :
  --   Appends the list of shifted points w.r.t. x to the list res.

    tmp : List := l;

  begin
    tmp := l;
    while not Is_Null(tmp) loop
      Append_Diff(res,res_last,Head_Of(tmp).all-x);
      tmp := Tail_Of(tmp);
    end loop;
  end Shifted_Points;

  function In_Hull ( l : List; x : Vector ) return boolean is

  -- DESCRIPTION :
  --   Returns true if the vector x belongs the convex hull spanned by
  --   the points in l minus the point x itself.

    res : boolean;
    lx : List;

  begin
    Copy(l,lx);
    Remove(lx,x);
    res := Is_In_Hull(x,lx);
    Clear(lx);
    return res;
  end In_Hull;

-- AUXILIAIRIES TO CONSTRUCT THE NORMALS :

  procedure Normal ( m : in out Matrix; v : out Vector ) is

  -- DESCRIPTION :
  --   Computes the normal to the vectors stored in the rows of the matrix m.

  -- REQUIRED : the matrix m is square!

    res : Vector(m'range(2));

  begin
    Upper_Triangulate(m); Scale(m);
    res := (res'range => 0);
    Solve0(m,res);
    v := res;
  end Normal;

  procedure Orientate_Normal ( l : in List; f : in Face; 
                               normal : in out Vector ) is

  -- DESCRIPTION :
  --   Orientates the normal of the face such that it becomes an inner normal
  --   w.r.t. the points in the list l.

    tmp : List := l;
    ip : constant integer := f(f'first).all*normal;
    pt : Vector(Head_Of(l)'range);
    done : boolean := false;
    ptip : integer;

  begin
    while not Is_Null(tmp) loop
      pt := Head_Of(tmp).all;
      if not Is_In(f,pt)
       then ptip := pt*normal;
            if ptip /= ip
             then if ptip < ip
                   then normal := -normal;
                  end if;
                  done := true;
            end if;
      end if;
      exit when done;
      tmp := Tail_Of(tmp);
    end loop;
  end Orientate_Normal;

  function Inner_Normal ( l : List; facet : Face ) return Vector is

  -- DESCRIPTION :
  --   Returns the inner normal to the facet.

    res,fst : Vector(facet(facet'first)'range);
    m : Matrix(facet'first+1..facet'last,facet(facet'first)'range);

  begin
    fst := facet(facet'first).all;
    for i in m'first(1)..m'last(1) loop
      for j in m'range(2) loop
        m(i,j) := facet(i)(j) - fst(j);
      end loop;
    end loop;
    Normal(m,res);
    Orientate_Normal(l,facet,res);
    return res;
  end Inner_Normal;

  procedure Inner_Normals ( l : in List; facets : in Faces;
                            iv,iv_last : in out List ) is

  -- DESCRIPTION :
  --   Returns the list of all inner normals to the facets of conv(l).
  --   If there is only one facet, then both `inner' and `outer' normal
  --   will be returned, because there is nothing to orientate with.
  --   This means that also the negation of the normal will be in the list iv.

    tmp : Faces := facets;

  begin
    while not Is_Null(tmp) loop
      declare
        f : Face := Head_Of(tmp);
        v : constant Vector := Inner_Normal(l,Head_Of(tmp));
      begin
        Append(iv,iv_last,v);
        if Length_Of(facets) = 1
         then Append(iv,iv_last,-v);
        end if;
      end;
      tmp := Tail_Of(tmp);
    end loop;
  end Inner_Normals;

  function Number_of_Rows ( l : List ) return natural is

  -- DESCRIPTION :
  --   Returns Maximum(dimension,Length_Of(l)).

  -- REQUIRED : not Is_Null(l).

    n : constant natural := Head_Of(l)'length;
    m : constant natural := Length_Of(l);

  begin
    if n > m
     then return n;
     else return m;
    end if;
  end Number_of_Rows;

  function Normal ( l : List ) return Vector is

  -- DESCRIPTION :
  --   Return a normal to the vectors in the list.

  -- REQUIRED : Length_Of(l) <= n = Head_Of(l)'length
  --   or the points in l all lie on the same facet.

    fst : Link_to_Vector := Head_Of(l);
    m : Matrix(1..Number_of_Rows(l),fst'range);
    cnt : natural := 0;
    res : Vector(fst'range);
    tmp : List := Tail_Of(l);

  begin
    while not Is_Null(tmp) loop
      res := Head_Of(tmp).all;
      cnt := cnt+1;
      for i in res'range loop
        m(cnt,i) := res(i) - fst(i);
      end loop;
      tmp := Tail_Of(tmp);
    end loop;
    for i in cnt+1..m'last(1) loop
      for j in m'range(2) loop
        m(i,j) := 0;
      end loop;
    end loop;
    Normal(m,res);
    return res;
  end Normal;

  procedure Normals ( l : in List; x : in Vector; iv,iv_last : in out List ) is

  -- DESCRIPTION :
  --   Computes a list of all normals to the vectors in the list.

  -- REQUIRED : Length_Of(l) <= n = x'length.
  --   or all points in l lie on the same facet.

    nor : constant Vector := Normal(l);
    piv : natural := Pivot(nor);
    pivnor : Vector(nor'range);   -- normal with pivnor(piv) > 0

  begin
    if piv <= nor'last
     then if x'length > 2
           then if nor(piv) < 0
                 then pivnor := -nor;
                 else pivnor := nor;
                end if;
                declare
                  t : Transfo := Build_Transfo(pivnor,piv);
                  tx : constant Vector := Reduce(t*x,piv);
                  tt : Transfo := Transpose(t);
                  tl : List := Transform_and_Reduce(t,piv,l);
                begin
                  if Length_Of(tl) <= nor'length-1
                   then Normals(tl,tx,iv,iv_last);
                   else declare
                          facets : Faces := Compute_Facets(tl,tx);
                        begin
                          if Length_Of(facets) > 1
                           then Inner_Normals(tl,facets,iv,iv_last);
                           else Normals(tl,tx,iv,iv_last);
                          end if;
                        end;
                  end if;
                  Insert_and_Transform(iv,piv,0,tt);
                  iv_last := Pointer_to_Last(iv);
                  Clear(t); Deep_Clear(tl); Clear(tt);
                end;
          end if;
          Append(iv,iv_last,nor); Append(iv,iv_last,-nor);
    end if;
  end Normals;

-- CONSTRUCTORS FOR PRIMAL REPRESENTATION :

  function Generators ( l : List; facets : Faces; x : Vector ) return List is

    res,res_last : List;

  begin
    if Length_Of(facets) > 1
     then Inner_Normals(l,facets,res,res_last);     -- compute inner normals
     elsif In_Hull(l,x)
         then return res;                           -- empty normal cone
         else Normals(l,x,res,res_last);            -- compute normals
              if Length_Of(l) = 2
               then Shifted_Points(l,x,res,res_last);
              end if;
    end if;
   -- SAFETY MODE :
   -- if not Contained_in_Cone(l,x,res)
   --  then put_line("Bug raised for computing generators for list"); put(l);
   --       put(" and the point : "); put(x); new_line;
   --       raise PROGRAM_ERROR;
   -- end if;
    return res;
  end Generators;

  function Generators ( l : List; x : Vector ) return List is

    res,res_last : List;
    n : constant natural := x'length;
    facets : Faces;

  begin
    if In_Hull(l,x)
     then return res;                                -- empty normal cone
     else if Length_Of(l) > n
           then facets := Compute_Facets(l,x);
          end if;
          if Length_Of(facets) > 1
           then res := Generators(l,facets,x);
           else Normals(l,x,res,res_last);
                if Length_Of(l) = 2
                 then Shifted_Points(l,x,res,res_last);
                end if;
          end if;
         -- SAFETY MODE :
         -- if not Contained_in_Cone(l,x,res)
         --  then put_line("Bug raised for computing generators for list");
         --       put(l); put(" and the point : "); put(x); new_line;
         --       raise PROGRAM_ERROR;
         -- end if;
          return res;
    end if;
  end Generators;

-- CONSTRUCTOR FOR DUAL REPRESENTATION :

  function Inner_Normal_Cone ( l : List; x : Vector ) return Matrix is

    res : Matrix(x'range,1..Length_Of(l)-1);
       -- CAUTION : Is_In(l,x) must hold !!
    y : Vector(x'range);
    cnt : natural := 0;
    tmp : List := l;

  begin
    while not Is_Null(tmp) loop
      y := Head_Of(tmp).all;
      if not Equal(x,y)                  -- avoid trivial inequality
       then cnt := cnt+1;
            for i in x'range loop
              res(i,cnt) := y(i) - x(i);
            end loop;
      end if;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Inner_Normal_Cone;

  function Included_Normal_Cone ( gx : List; x : Vector ) return Matrix is

    res : Matrix(x'first-1..x'last,1..Length_Of(gx));
    tmp : List := gx;
    v : Vector(x'range);

  begin
    for j in res'range(2) loop
      v := Head_Of(tmp).all;
      res(res'first(1),j) := x*v;
      for i in res'first(1)+1..res'last(1) loop
        res(i,j) := v(i);
      end loop;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Included_Normal_Cone;

-- PRIMITIVE SELECTORS :

  function Evaluate ( m : Matrix; i : integer; v : Vector ) return integer is

    sum : integer := 0;

  begin
    for j in v'range loop
      sum := sum + m(j,i)*v(j);
    end loop;
    return sum;
  end Evaluate;

  function Satisfies ( m : Matrix; i : integer; v : Vector ) return boolean is
  begin
    if m'first(1) = v'first-1
     then return (Evaluate(m,i,v) >= m(0,i));
     else return (Evaluate(m,i,v) >= 0);
    end if;
  end Satisfies;

  function Strictly_Satisfies ( m : Matrix; i : integer; v : Vector )
                              return boolean is
  begin
    if m'first(1) = v'first-1
     then return (Evaluate(m,i,v) > m(0,i));
     else return (Evaluate(m,i,v) > 0);
    end if;
  end Strictly_Satisfies;

  function Satisfies ( m : Matrix; v : Vector ) return boolean is
  begin
    for i in m'range(2) loop
      if not Satisfies(m,i,v)
       then return false;
      end if;
    end loop;
    return true;
  end Satisfies;

  function Strictly_Satisfies ( m : Matrix; v : Vector ) return boolean is
  begin
    for i in m'range(2) loop
      if not Strictly_Satisfies(m,i,v)
       then return false;
      end if;
    end loop;
    return true;
  end Strictly_Satisfies;

  function Satisfies ( m : Matrix; v : List ) return boolean is

    tmp : List := v;

  begin
    while not Is_Null(tmp) loop
      if not Satisfies(m,Head_Of(tmp).all)
       then return false;
       else tmp := Tail_Of(tmp);
      end if;
    end loop;
    return true;
  end Satisfies;

  function Strictly_Satisfies ( m : Matrix; v : List ) return boolean is

    tmp : List := v;

  begin
    while not Is_Null(tmp) loop
      if not Strictly_Satisfies(m,Head_Of(tmp).all)
       then return false;
       else tmp := Tail_Of(tmp);
      end if;
    end loop;
    return true;
  end Strictly_Satisfies;

-- SECONDARY SELECTORS :

  function Satisfy ( m : Matrix; l : List ) return List is

    res,res_last,tmp : List;

  begin
    tmp := l;
    while not Is_Null(tmp) loop
      if Satisfies(m,Head_Of(tmp).all)
       then Append(res,res_last,Head_Of(tmp).all);
      end if;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Satisfy;

  function Strictly_Satisfy ( m : Matrix; l : List ) return List is

    res,res_last,tmp : List;

  begin
    tmp := l;
    while not Is_Null(tmp) loop
      if Strictly_Satisfies(m,Head_Of(tmp).all)
       then Append(res,res_last,Head_Of(tmp).all);
      end if;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Strictly_Satisfy;

  function Contained_in_Cone ( l : List; x : Vector; v : List )
                             return boolean is

    ineq : Matrix(x'range,1..Length_Of(l)-1) := Inner_Normal_Cone(l,x);

  begin
    return Satisfies(ineq,v);
  end Contained_in_Cone;

  function Strictly_Contained_in_Cone ( l : List; x : Vector; v : List )
                                      return boolean is

    ineq : Matrix(x'range,1..Length_Of(l)-1) := Inner_Normal_Cone(l,x);

  begin
    return Strictly_Satisfies(ineq,v);
  end Strictly_Contained_in_Cone;

  function Contained_in_Cone ( l,v : List ) return boolean is

    tmp : List := l;

  begin
    while not Is_Null(tmp) loop
      if Contained_in_Cone(l,Head_Of(tmp).all,v)
       then -- put("true for vector : "); put(Head_Of(tmp).all); new_line;
            return true;
       else tmp := Tail_Of(tmp);
      end if;
    end loop;
    return false;
  end Contained_in_Cone;

  function Strictly_Contained_in_Cone ( l,v : List ) return boolean is

    tmp : List := l;

  begin
    while not Is_Null(tmp) loop
      if Strictly_Contained_in_Cone(l,Head_Of(tmp).all,v)
       then return true;
       else tmp := Tail_Of(tmp);
      end if;
    end loop;
    return false;
  end Strictly_Contained_in_Cone;

-- CONCERNING THE UNION OF CONES :

  function In_Union ( v1,v2 : Vector; k1,k2 : Matrix ) return boolean is

  -- ALGORITHM : consider x(t) = (1-t)*v1 + t*v2, for t in [0,1]
  --   and compute t1: smallest t such that Satisfies(k1,x(t))
  --   and t2: largest t such that Satisfies(k2,x(t)).
  --   Then t1 >= t2 makes the test returning true, false otherwise.

    diff : constant Vector := v1-v2;
    t1a,t1b,t2a,t2b,a,b : integer;          -- t1 = t1a/t1b and t2 = t2a/t2b

  begin
    t1a := 1; t1b := 1;                            -- initially t1 = 1 = 1/1
    for i in k1'range(2) loop
      b := Evaluate(k1,i,diff);
      if b /= 0
       then a := Evaluate(k1,i,v1);
            if Lower(a,b,t1a,t1b)
             then t1a := a; t1b := b;
            end if;
      end if;
    end loop;
    t2a := 0; t2b := 1;                            -- initially t2 = 0 = 0/1
    for i in k2'range(2) loop
      b := Evaluate(k2,i,diff);
      if b /= 0
       then a := Evaluate(k2,i,v1);
            if Higher(a,b,t2a,t2b)
             then t2a := a; t2b := b;
            end if;
      end if;
    end loop;
   -- put("In_Union for "); put(v1); put(" and "); put(v2);
   -- put(" : t1 = "); put(t1a,1); put("/"); put(t1b,1);
   -- put(" : t2 = "); put(t2a,1); put("/"); put(t2b,1); new_line;
   -- put_line("k1 : "); put(k1); put_line("k2 : "); put(k2);
    return not Lower(t1a,t1b,t2a,t2b);
  end In_Union;

  function In_Union ( v1,v2 : List; k1,k2 : Matrix ) return boolean is

    tmpv1,tmpv2 : List;

  begin
    tmpv1 := v1;
    while not Is_Null(tmpv1) loop
      tmpv2 := v2;
      while not Is_Null(tmpv2) loop
        if not In_Union(Head_Of(tmpv1).all,Head_Of(tmpv2).all,k1,k2)
         then return false;
         else tmpv2 := Tail_Of(tmpv2);
        end if;
      end loop;
      tmpv1 := Tail_Of(tmpv1);
    end loop;
    return true;
  end In_Union;

end Inner_Normal_Cones;