File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports / floating_face_enumerators.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:27 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD
Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with Standard_Floating_Vectors; use Standard_Floating_Vectors;
with Standard_Floating_Matrices; use Standard_Floating_Matrices;
with Givens_Rotations; use Givens_Rotations;
with Floating_Linear_Inequalities; use Floating_Linear_Inequalities;
package body Floating_Face_Enumerators is
-- AUXILIARIES :
function Is_In ( x : integer; v : Standard_Integer_Vectors.Vector )
return boolean is
-- DESCRIPTION :
-- Returns true if there exists an entry k in v, with v(k) = x.
begin
for k in v'range loop
if v(k) = x
then return true;
end if;
end loop;
return false;
end Is_In;
function Is_Zero ( v : Standard_Floating_Vectors.Vector; tol : double_float )
return boolean is
-- DESCRIPTION :
-- Returns true if abs(v(i)) < tol, for all i in v'range.
begin
for i in v'range loop
if abs(v(i)) > tol
then return false;
end if;
end loop;
return true;
end Is_Zero;
function In_Edge ( x,a,b : Standard_Floating_Vectors.Vector;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- Returns true if the vector x lies between a and b.
ba,xa : double_float;
piv : integer;
begin
for i in b'range loop -- look for first i: b(i) - a(i) /= 0
ba := b(i) - a(i);
if abs(ba) > tol
then piv := i; exit;
end if;
end loop;
if abs(ba) < tol
then return Equal(x,a); -- in this case b = a, so test x = a
else if ba < 0.0
then ba := -ba;
xa := x(piv) - b(piv);
else xa := x(piv) - a(piv);
end if;
if xa*ba >= 0.0 and then xa <= ba
then return true;
else return false;
end if;
end if;
end In_Edge;
function Convex_Decomposition
( k : natural; face : Standard_Integer_Vectors.Vector;
x : Standard_Floating_Vectors.Vector;
pts : Standard_Floating_VecVecs.VecVec;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- Decomposes the vector x in terms of the k+1 points of pts(face(*)).
-- Returns true if x can be written as a convex combination.
mat : Matrix(x'range,1..k);
rhs : Standard_Floating_Vectors.Vector(x'range) := x;
sol : Standard_Floating_Vectors.Vector(1..k) := (1..k => 0.0);
ipvt : Standard_Integer_Vectors.Vector(1..k);
sum : double_float := 0.0;
begin
for j in 1..k loop
for i in mat'range(1) loop
mat(i,j) := pts(face(j))(i) - pts(face(0))(i);
end loop;
ipvt(j) := j;
end loop;
Upper_Triangulate(mat,rhs,tol,ipvt);
for j in k+1..rhs'last loop
if abs(rhs(j)) > tol
then return false;
end if;
end loop;
Solve(mat,rhs,tol,sol);
for j in sol'range loop
if sol(j) < -tol
then return false;
else sum := sum + sol(j);
end if;
end loop;
return (abs(sum-1.0) < tol);
end Convex_Decomposition;
function In_Face ( k : in natural; face : Standard_Integer_Vectors.Vector;
x : Standard_Floating_Vectors.Vector;
pts : Standard_Floating_VecVecs.VecVec;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- Returns true if x lies in the given k-face, which contains entries
-- to its elements in pts. This means that we have to verify whether
-- the vector x can be written as a convex combination of the points
-- in the face.
begin
if k = 0
then return Equal(pts(face(face'first)).all,x);
elsif k = 1
then return In_Edge(x,pts(face(face'first)).all,
pts(face(face'first+1)).all,tol);
else return Convex_Decomposition(k,face,x,pts,tol);
end if;
end In_Face;
-- ELIMINATE TO MAKE UPPER TRIANGULAR :
procedure Eliminate ( pivot : in integer;
elim : in Standard_Floating_Vectors.Vector;
target : in out Standard_Floating_Vectors.Vector ) is
-- DESCRIPTION :
-- Makes target(pivot) = 0.0 by means of making a linear
-- combination of the vectors target and elim.
-- Note that target is viewed as an inequality, so that we make sure
-- to multiply it only with a positive number.
-- REQUIRED : abs(target(pivot)) > tol and abs(elim(pivot)) > tol.
a,b,tarfac,elifac : double_float;
begin
a := elim(pivot); b := target(pivot);
if a > 0.0
then tarfac := a; elifac := b;
else tarfac := -a; elifac := -b;
end if;
for j in target'range loop
target(j) := tarfac*target(j) - elifac*elim(j);
end loop;
end Eliminate;
procedure Eliminate ( l : in natural;
pivots : in Standard_Integer_Vectors.Vector;
elim : in Standard_Floating_VecVecs.VecVec;
tol : in double_float;
target : in out Standard_Floating_Vectors.Vector ) is
-- DESCRIPTION :
-- Makes target(pivots(i)) = 0 by means of making a linear
-- combination of the vectors target and elim(i), for i in 1..l.
-- REQUIRED : elim(i)(pivots(i))) > tol.
begin
for i in 1..l loop
if abs(target(pivots(i))) > tol
then Eliminate(pivots(i),elim(i).all,target);
end if;
end loop;
end Eliminate;
function Pivot_after_Elimination
( l,k : in natural; point : Standard_Floating_Vectors.Vector;
face,pivots : Standard_Integer_Vectors.Vector;
elim : Standard_Floating_VecVecs.VecVec; tol : double_float )
return natural is
-- DESCRIPTION :
-- Returns the first nonzero element of the given point after elimination
-- w.r.t. the entries in the face with lower index.
work : Standard_Floating_Vectors.Vector(point'range)
:= point - elim(1-k).all;
pivot : integer;
begin
for i in (face'first+1)..face'last loop
if (face(i) < l) and then (abs(work(pivots(i))) > tol)
then Eliminate(pivots(i),elim(i).all,work);
end if;
exit when (face(i) > l);
end loop;
pivot := 0;
for i in work'range loop
if abs(work(pivots(i))) > tol
then pivot := i;
end if;
exit when (pivot /= 0);
end loop;
return pivot;
end Pivot_after_Elimination;
procedure Update_Pivots
( point : in Standard_Floating_Vectors.Vector; l : in natural;
pivots : in out Standard_Integer_Vectors.Vector;
tol : in double_float; pivot : out natural ) is
-- DESCRIPTION :
-- Searches in the point(l..point'last) for the first nonzero entry
-- and updates the pivoting information.
temp,piv : integer;
begin
piv := 0;
for i in l..point'last loop
if abs(point(pivots(i))) > tol
then piv := i;
end if;
exit when (piv /= 0);
end loop;
if piv /= 0 and then (piv /= l)
then temp := pivots(l);
pivots(l) := pivots(piv);
pivots(piv) := temp;
end if;
pivot := piv;
end Update_Pivots;
procedure Update_Eliminator
( elim : in out Standard_Floating_VecVecs.VecVec;
l : in natural;
pivots : in out Standard_Integer_Vectors.Vector;
point : in Standard_Floating_Vectors.Vector;
tol : in double_float;
pivot : out natural ) is
-- DESCRIPTION :
-- Updates the vector of eliminators, by adding the lth elimination
-- equation. This procedure ensures the invariant condition on the
-- eliminator, which has to be upper triangular. If this cannot be
-- achieved, degeneracy is indicated by pivot = 0.
-- ON ENTRY :
-- elim vector of elimination equations: elim(i)(pivots(i)) /= 0
-- and for j in 1..(i-1) : elim(i)(pivots(j)) = 0,
-- for i in 1..(l-1), elim(0) contains the basis point;
-- l index of current elimination equation to be made;
-- pivots contains the pivoting information;
-- tol tolerance on the precision;
-- point new point to make the equation `point - elim(0) = 0'.
-- ON RETURN :
-- elim if not degen, then elim(l)(pivots(l)) /= 0 and
-- for i in 1..(l-1): elim(l)(pivots(i)) = 0;
-- pivots updated pivot information;
-- pivot the pivot that has been used for elim(l)(pivots(l)) /= 0;
-- piv = 0 when the new elimination equation elim(l)
-- became entirely zero after ensuring the invariant condition.
begin
elim(l) := new Standard_Floating_Vectors.Vector'(point - elim(0).all);
Eliminate(l-1,pivots,elim,tol,elim(l).all);
Update_Pivots(elim(l).all,l,pivots,tol,pivot);
end Update_Eliminator;
procedure Update_Eliminator_for_Sum
( elim : in out Standard_Floating_VecVecs.VecVec;
l : in natural;
pivots : in out Standard_Integer_Vectors.Vector;
point : in Standard_Floating_Vectors.Vector;
index : in natural;
tol : in double_float; pivot : out natural ) is
-- DESCRIPTION :
-- Updates the vector of eliminators, by adding the lth elimination
-- equation. This procedure ensures the invariant condition on the
-- eliminator, which has to be upper triangular. If this cannot be
-- achieved, degeneracy is indicated by pivot = 0.
-- ON ENTRY :
-- elim vector of elimination equations: elim(i)(pivots(i)) /= 0
-- and for j in 1..(i-1) : elim(i)(pivots(j)) = 0,
-- for i in 1..(l-1), elim(1-index) contains the basis point;
-- l index of current elimination equation to be made;
-- pivots contains the pivoting information;
-- point new point to make the equation `point - elim(1-index) = 0';
-- index indicates the current polytope;
-- tol tolerance for precision.
-- ON RETURN :
-- elim if not degen, then elim(l)(pivots(l)) /= 0 and
-- for i in 1..(l-1): elim(l)(pivots(i)) = 0;
-- pivots updated pivot information;
-- piv the pivot that has been used for elim(l)(pivots(l)) /= 0;
-- piv = 0 when the new elimination equation elim(l)
-- became entirely zero after ensuring the invariant condition.
begin
elim(l) := new Standard_Floating_Vectors.Vector'(point - elim(1-index).all);
Eliminate(l-1,pivots,elim,tol,elim(l).all);
Update_Pivots(elim(l).all,l,pivots,tol,pivot);
end Update_Eliminator_for_Sum;
-- CREATE TABLEAU OF INEQUALITIES :
procedure Create_Tableau_for_Vertices
( i,n : in integer;
pts : in Standard_Floating_VecVecs.VecVec;
tab : out Matrix;
rhs : out Standard_Floating_Vectors.Vector ) is
-- DESCRIPTION :
-- Creates a system of linear inequalities that is feasible when
-- the ith point is spanned by the other points in pts.
column : integer := tab'first(2);
begin
for j in pts'range loop
if j /= i
then for k in pts(j)'range loop
tab(k,column) := pts(j)(k);
end loop;
tab(tab'last(1),column) := 1.0;
column := column + 1;
end if;
end loop;
for k in pts(i)'range loop
rhs(k) := pts(i)(k);
end loop;
rhs(n+1) := 1.0;
end Create_Tableau_for_Vertices;
procedure Create_Tableau_for_Edges
( i,j,n,pivot : in integer;
elim : in Standard_Floating_Vectors.Vector;
pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float; tab : out matrix;
rhs : out Standard_Floating_Vectors.Vector;
lastcol : out integer; degenerate : out boolean ) is
-- DESCRIPTION :
-- Creates a system of linear inequalities that express the fact
-- that the points i and j span an edge of the convex hull.
-- Degenerate means that there exists a point on the edge spanned
-- by the points i and j that contains the tested edge.
column : integer := 1;
ineq : Standard_Floating_Vectors.Vector(1..n);
begin
degenerate := false;
for k in pts'range loop
if (k/=i) and then (k/=j)
then ineq := pts(k).all - pts(i).all; -- compute inequality
if abs(ineq(pivot)) > tol -- make ineq(pivot) = 0
then Eliminate(pivot,elim,ineq);
end if;
if Is_Zero(ineq,tol) -- check if degenerate
then if not In_Edge(pts(k).all,pts(i).all,pts(j).all,tol)
then degenerate := true; return;
end if;
else for l in 1..n loop -- fill in the column
if l < pivot
then tab(l,column) := ineq(l);
elsif l > pivot
then tab(l-1,column) := ineq(l);
end if;
end loop;
tab(tab'last(1),column) := 1.0;
column := column + 1;
end if;
end if;
end loop;
for k in tab'first(1)..(tab'last(1)-1) loop -- right hand side
rhs(k) := 0.0;
end loop;
rhs(tab'last(1)) := 1.0;
lastcol := column-1;
end Create_Tableau_for_Edges;
procedure Create_Tableau_for_Faces
( k,n : in natural;
face,pivots : in Standard_Integer_Vectors.Vector;
pts,elim : in Standard_Floating_VecVecs.VecVec;
tol : in double_float; tab : out Matrix;
rhs : out Standard_Floating_Vectors.Vector;
lastcol : out integer; degenerate : out boolean ) is
-- DESCRIPTION :
-- Creates a system of linear inequalities that express the fact
-- that the points pts(face(*)) span a face of the convex hull.
-- Degenerate means that there exists a point on the face spanned
-- by the points in the face that contains the tested face.
column : integer := 1;
ineq : Standard_Floating_Vectors.Vector(1..n);
begin
degenerate := false;
for l in pts'range loop
if not Is_In(l,face)
then ineq := pts(l).all - elim(0).all; -- new inequality
Eliminate(k,pivots,elim,tol,ineq); -- ineq(pivots(i)) = 0, i=1,2,..k
if Is_Zero(ineq,tol)
then
if not In_Face(k,face,pts(l).all,pts,tol)
and then l < face(face'last) -- lexicographic enumeration
and then (Pivot_after_Elimination
(l,1,pts(l).all,face,pivots,elim,tol) /= 0)
then degenerate := true; return;
end if;
else
for ll in (k+1)..n loop -- fill in the column
tab(ll-k,column) := ineq(pivots(ll));
end loop;
tab(tab'last(1),column) := 1.0;
column := column + 1;
end if;
end if;
end loop;
for l in tab'first(1)..(tab'last(1)-1) loop -- right hand side
rhs(l) := 0.0;
end loop;
rhs(tab'last(1)) := 1.0;
lastcol := column-1;
end Create_Tableau_for_Faces;
procedure Create_Tableau_for_Faces_of_Sum
( k,n,i,r : in natural;
ind,pivots : in Standard_Integer_Vectors.Vector;
pts,elim : Standard_Floating_VecVecs.VecVec;
tol : double_float;
face : in Standard_Integer_VecVecs.VecVec;
tab : out Matrix; rhs : out Standard_Floating_Vectors.Vector;
lastcol : out integer; degenerate : out boolean ) is
-- DESCRIPTION :
-- Creates the table of inequalities for determining whether the given
-- candidate face spans a face of the sum polytope.
-- ON ENTRY :
-- k dimension of the face on the sum;
-- n dimension of the points;
-- i current polytope;
-- r number of polytopes;
-- ind indicate the beginning vector in pts of each polytope;
-- etc...
column : integer := 1;
ineq : Standard_Floating_Vectors.Vector(1..n);
last : integer;
begin
degenerate := false;
for l1 in face'first..i loop
if l1 = r
then last := pts'last;
else last := ind(l1+1)-1;
end if;
for l2 in ind(l1)..last loop
if not Is_In(l2,face(l1).all)
then
ineq := pts(l2).all - elim(1-l1).all; -- new inequality
Eliminate(k,pivots,elim,tol,ineq); -- ineq(pivots(i)) = 0, i=1,2,..k
if Is_Zero(ineq,tol)
then
if not In_Face(face(l1)'length-1,face(l1).all,pts(l2).all,pts,tol)
and then face(l1)(face(l1)'first) <= l2
and then l2 < face(l1)(face(l1)'last)
-- lexicographic enumeration
and then (Pivot_after_Elimination
(l2,l1,pts(l2).all,face(l1).all,pivots,elim,tol) /= 0)
then degenerate := true; return;
end if;
else
for ll in (k+1)..n loop -- fill in the column
tab(ll-k,column) := ineq(pivots(ll));
end loop;
tab(tab'last(1),column) := 1.0;
column := column + 1;
end if;
end if;
end loop;
end loop;
for l in tab'first(1)..(tab'last(1)-1) loop -- right hand side
rhs(l) := 0.0;
end loop;
rhs(tab'last(1)) := 1.0;
lastcol := column-1;
end Create_Tableau_for_Faces_of_Sum;
procedure Create_Tableau_for_Lower_Edges
( i,j,n,pivot : in integer;
elim : in Standard_Floating_Vectors.Vector;
pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float; tab : out Matrix;
rhs : out Standard_Floating_Vectors.Vector;
lastcol : out integer; degenerate : out boolean ) is
-- DESCRIPTION :
-- Creates a system of linear inequalities that express the fact
-- that the points i and j spann an edge of the lower hull.
-- Degenerate means that there exists a point on the edge spanned
-- by the points i and j that contains the tested edge.
column : integer := 1;
ineq : Standard_Floating_Vectors.Vector(1..n);
begin
degenerate := false;
for k in pts'range loop
if (k/=i) and then (k/=j)
then ineq := pts(k).all - pts(i).all; -- compute inequality
if abs(ineq(pivot)) > tol -- make ineq(pivot) = 0
then Eliminate(pivot,elim,ineq);
end if;
if Is_Zero(ineq,tol) -- check if degenerate
then if not In_Edge(pts(k).all,pts(i).all,pts(j).all,tol)
then degenerate := true; return;
end if;
else for l in 1..n loop -- fill in the column
if l < pivot
then tab(l,column) := ineq(l);
elsif l > pivot
then tab(l-1,column) := ineq(l);
end if;
end loop;
column := column + 1;
end if;
end if;
end loop;
for k in tab'first(1)..(tab'last(1)-1) loop -- right hand side
rhs(k) := 0.0;
end loop;
rhs(tab'last(1)) := -1.0;
lastcol := column-1;
end Create_Tableau_for_Lower_Edges;
procedure Create_Tableau_for_Lower_Faces
( k,n : in natural;
face,pivots : in Standard_Integer_Vectors.Vector;
pts,elim : in Standard_Floating_VecVecs.VecVec;
tol : in double_float; tab : out Matrix;
rhs : out Standard_Floating_Vectors.Vector;
lastcol : out integer; degenerate : out boolean ) is
column : integer := 1;
ineq : Standard_Floating_Vectors.Vector(1..n);
begin
degenerate := false;
for l in pts'range loop
if not Is_In(l,face)
then ineq := pts(l).all - elim(0).all; -- new inequality
Eliminate(k,pivots,elim,tol,ineq); -- ineq(pivots(i)) = 0, i=1,2,..k
if Is_Zero(ineq,tol)
then
if not In_Face(k,face,pts(l).all,pts,tol)
and then l < face(face'last) -- lexicographic enumeration
and then (Pivot_after_Elimination
(l,1,pts(l).all,face,pivots,elim,tol) /= 0)
then degenerate := true; return;
end if;
else
for ll in (k+1)..n loop -- fill in the column
tab(ll-k,column) := ineq(pivots(ll));
end loop;
column := column + 1;
end if;
end if;
end loop;
for l in tab'first(1)..(tab'last(1)-1) loop -- right hand side
rhs(l) := 0.0;
end loop;
rhs(tab'last(1)) := -1.0;
lastcol := column-1;
end Create_Tableau_for_Lower_Faces;
-- DETERMINE WHETHER THE CANDIDATE IS VERTEX, SPANS AN EDGE OR A K-FACE :
function Is_Vertex ( i,m,n : integer;
pts : Standard_Floating_VecVecs.VecVec;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- This function determines whether the ith point is a vertex of
-- the m n-dimensional points in pts, by deciding on the feasibility
-- of a linear-inequality system.
tableau : Matrix(1..n+1,1..m);
rhs : Standard_Floating_Vectors.Vector(tableau'range(1));
sol : Standard_Floating_Vectors.Vector(tableau'range(1));
col : Standard_Integer_Vectors.Vector(tableau'range(1));
feasible : boolean;
begin
Create_Tableau_for_Vertices(i,n,pts,tableau,rhs);
Complementary_Slackness(tableau,rhs,tol,sol,col,feasible);
return not feasible;
end Is_Vertex;
function Is_Edge ( i,j,m,n : integer;
pts : Standard_Floating_VecVecs.VecVec;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- This function determines whether the ith and jth point span an
-- edge of the convex hull of m n-dimensional points in pts,
-- by deciding on the feasibility of a linear-inequality system.
elim : Standard_Floating_Vectors.Vector(1..n) := pts(i).all - pts(j).all;
pivot : integer;
begin
pivot := elim'first - 1;
for k in elim'range loop
if abs(elim(k)) > tol
then pivot := k;
end if;
exit when pivot >= elim'first;
end loop;
if pivot < elim'first
then return false;
else declare
tab : Matrix(1..n,1..m-1);
rhs : Standard_Floating_Vectors.Vector(tab'range(1));
sol : Standard_Floating_Vectors.Vector(tab'range(1));
col : Standard_Integer_Vectors.Vector(tab'range(1));
deg,feasible : boolean;
lst : integer;
begin
Create_Tableau_for_Edges(i,j,n,pivot,elim,pts,tol,tab,rhs,lst,deg);
if deg
then return false;
elsif lst = 0
then return true;
else Complementary_Slackness(tab,lst,rhs,tol,sol,col,feasible);
return not feasible;
end if;
end;
end if;
end Is_Edge;
function Is_Face ( k,n,m : natural;
elim,pts : Standard_Floating_VecVecs.VecVec;
pivots,face : Standard_Integer_Vectors.Vector;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- Applies Complementary Slackness to determine whether the given
-- candidate face is a face of the polytope.
-- ON ENTRY :
-- k dimension of the candidate face;
-- n dimension of the vector space;
-- m number of points which span the polytope;
-- elim elimination equations in upper triangular form;
-- pts the points which span the polytope;
-- pivots pivoting information for the elimination equations;
-- face entries of the points which span the candidate face;
-- tol tolerance on the precision.
nn : constant natural := n-k+1;
tab : Matrix(1..nn,1..(m-k));
rhs : Standard_Floating_Vectors.Vector(tab'range(1));
sol : Standard_Floating_Vectors.Vector(tab'range(1));
col : Standard_Integer_Vectors.Vector(tab'range(1));
deg,feasible : boolean;
lst : integer;
begin
-- put("The face being tested : "); put(face); new_line;
-- put("The pivots : "); put(pivots); new_line;
-- put_line("The elimination equations : ");
-- for i in 1..elim'last loop
-- put(elim(i)); put(" = 0 ");
-- end loop;
-- new_line;
if m - k <= 1
then return true;
else
Create_Tableau_for_Faces(k,n,face,pivots,pts,elim,tol,tab,rhs,lst,deg);
if deg
then return false;
elsif lst = 0
then return true;
else Complementary_Slackness(tab,lst,rhs,tol,sol,col,feasible);
return not feasible;
end if;
end if;
end Is_Face;
function Is_Face_of_Sum
( k,n,m,i,r : natural;
elim,pts : Standard_Floating_VecVecs.VecVec;
tol : double_float;
face : Standard_Integer_VecVecs.VecVec;
ind,pivots : Standard_Integer_Vectors.Vector )
return boolean is
-- DESCRIPTION :
-- Applies Complementary Slackness to determine whether the given
-- candidate face is a face of the polytope.
-- ON ENTRY :
-- k dimension of the candidate face;
-- n dimension of the vector space;
-- m number of total points which span the polytope;
-- i current polytope;
-- r number of different polytopes;
-- elim elimination equations in upper triangular form;
-- pts the points which span the polytope;
-- tol tolerance on precision;
-- face entries of the points which span the candidate face;
-- ind indicates the starting vector in pts of each polytope;
-- pivots pivoting information for the elimination equations.
nn : constant natural := n-k+1;
tab : Matrix(1..nn,1..(m-k));
rhs : Standard_Floating_Vectors.Vector(tab'range(1));
sol : Standard_Floating_Vectors.Vector(tab'range(1));
col : Standard_Integer_Vectors.Vector(tab'range(1));
deg,feasible : boolean;
lst : integer;
begin
if m - k <= 1
then return true;
else
Create_Tableau_for_Faces_of_Sum
(k,n,i,r,ind,pivots,pts,elim,tol,face,tab,rhs,lst,deg);
if deg
then return false;
elsif lst = 0
then return true;
else Complementary_Slackness(tab,lst,rhs,tol,sol,col,feasible);
return not feasible;
end if;
end if;
end Is_Face_of_Sum;
function Is_Lower_Edge ( i,j,m,n : integer;
pts : Standard_Floating_VecVecs.VecVec;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- This function determines whether the ith and jth point span an
-- edge of the lower convex hull of m n-dimensional points in pts,
-- by deciding on the feasibility of a linear-inequality system.
elim : Standard_Floating_Vectors.Vector(1..n) := pts(i).all - pts(j).all;
pivot : integer;
begin
pivot := elim'first - 1;
for k in elim'range loop
if abs(elim(k)) > tol
then pivot := k;
end if;
exit when pivot >= elim'first;
end loop;
if pivot < elim'first or else (pivot = elim'last)
then return false;
else declare
tab : Matrix(1..n-1,1..m-1);
rhs : Standard_Floating_Vectors.Vector(tab'range(1));
sol : Standard_Floating_Vectors.Vector(tab'range(1));
col : Standard_Integer_Vectors.Vector(tab'range(1));
deg,feasible : boolean;
lst : integer;
begin
Create_Tableau_for_Lower_Edges
(i,j,n,pivot,elim,pts,tol,tab,rhs,lst,deg);
if deg
then return false;
elsif lst = 0
then return true;
else Complementary_Slackness(tab,lst,rhs,tol,sol,col,feasible);
return not feasible;
end if;
end;
end if;
end Is_Lower_Edge;
function Is_Lower_Face
( k,n,m : in natural;
elim,pts : Standard_Floating_VecVecs.VecVec;
pivots,face : Standard_Integer_Vectors.Vector;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- Applies Complementary Slackness to determine whether the given
-- candidate face is a face of the lower hull of the polytope.
-- ON ENTRY :
-- k dimension of the candidate face;
-- n dimension of the vector space;
-- m number of points which span the polytope;
-- elim elimination equations in upper triangular form;
-- pts the points which span the polytope;
-- pivots pivoting information for the elimination equations;
-- face entries of the points which span the candidate face.
nn : constant natural := n-k;
tab : Matrix(1..nn,1..(m-k));
rhs : Standard_Floating_Vectors.Vector(tab'range(1));
sol : Standard_Floating_Vectors.Vector(tab'range(1));
col : Standard_Integer_Vectors.Vector(tab'range(1));
deg,feasible : boolean;
lst : integer;
begin
-- put("The pivots : "); put(pivots); new_line;
-- put_line("The elimination equations : ");
-- for i in 1..elim'last loop
-- put(elim(i)); put(" = 0 ");
-- end loop;
-- new_line;
if m - k <= 1
then return true;
else Create_Tableau_for_Lower_Faces
(k,n,face,pivots,pts,elim,tol,tab,rhs,lst,deg);
if deg
then return false;
elsif lst = 0
then return true;
else Complementary_Slackness(tab,lst,rhs,tol,sol,col,feasible);
return not feasible;
end if;
end if;
end Is_Lower_Face;
-- TARGET ROUTINES :
procedure Enumerate_Vertices ( pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float ) is
continue : boolean := true;
n : constant natural := pts(pts'first).all'length;
m : constant natural := pts'length - 1;
begin
for i in pts'range loop
if Is_Vertex(i,m,n,pts,tol)
then Process(i,continue);
end if;
exit when not continue;
end loop;
end Enumerate_Vertices;
procedure Enumerate_Edges ( pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float ) is
continue : boolean := true;
n : constant natural := pts(pts'first).all'length;
m : constant natural := pts'length;
procedure Candidate_Edges ( i,n : in integer ) is
begin
for j in (i+1)..pts'last loop -- enumerate all candidates
-- put("Verifying "); put(pts(i).all); put(" and");
-- put(pts(j).all); put(" :"); new_line;
if Is_Edge(i,j,m,n,pts,tol)
then Process(i,j,continue);
end if;
exit when not continue;
end loop;
end Candidate_Edges;
begin
for i in pts'first..(pts'last-1) loop
Candidate_Edges(i,n);
end loop;
end Enumerate_Edges;
procedure Enumerate_Lower_Edges
( pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float ) is
continue : boolean := true;
n : constant natural := pts(pts'first).all'length;
m : constant natural := pts'length;
procedure Candidate_Lower_Edges ( i,n : in integer ) is
begin
for j in (i+1)..pts'last loop -- enumerate all candidates
-- put("Verifying "); put(pts(i).all); put(" and");
-- put(pts(j).all); put(" :"); new_line;
if Is_Lower_Edge(i,j,m,n,pts,tol)
then Process(i,j,continue);
end if;
exit when not continue;
end loop;
end Candidate_Lower_Edges;
begin
for i in pts'first..(pts'last-1) loop
Candidate_Lower_Edges(i,n);
end loop;
end Enumerate_Lower_Edges;
procedure Enumerate_Faces ( k : in natural;
pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float ) is
m : constant natural := pts'length;
n : constant natural := pts(pts'first).all'length;
candidate : Standard_Integer_Vectors.Vector(0..k);
elim : Standard_Floating_VecVecs.VecVec(0..k);
continue : boolean := true;
procedure Candidate_Faces ( ipvt : in Standard_Integer_Vectors.Vector;
i,l : in integer ) is
-- DESCRIPTION :
-- Enumerates all candidate k-faces in lexicographic order.
piv : natural;
pivots : Standard_Integer_Vectors.Vector(1..n);
begin
if l > k
then if (k = m)
or else Is_Face(k,n,m,elim,pts,ipvt,candidate,tol)
then Process(candidate,continue);
end if;
else for j in (i+1)..pts'last loop
candidate(l) := j;
pivots := ipvt;
Update_Eliminator(elim,l,pivots,pts(j).all,tol,piv);
if piv /= 0
then Candidate_Faces(pivots,j,l+1);
end if;
Clear(elim(l));
exit when not continue;
end loop;
end if;
end Candidate_Faces;
begin
if k <= m
then declare
ipvt : Standard_Integer_Vectors.Vector(1..n);
begin
for i in ipvt'range loop
ipvt(i) := i;
end loop;
for i in pts'first..(pts'last-k) loop
candidate(0) := i;
elim(0) := pts(i); -- basis point
Candidate_Faces(ipvt,i,1);
end loop;
end;
end if;
end Enumerate_Faces;
procedure Enumerate_Lower_Faces
( k : in natural;
pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float ) is
m : constant natural := pts'length;
n : constant natural := pts(pts'first).all'length;
candidate : Standard_Integer_Vectors.Vector(0..k);
elim : Standard_Floating_VecVecs.VecVec(0..k);
continue : boolean := true;
procedure Candidate_Lower_Faces ( ipvt : in Standard_Integer_Vectors.Vector;
i,l : in integer ) is
-- DESCRIPTION :
-- Enumerates all candidate k-faces in lexicographic order.
piv : natural;
pivots : Standard_Integer_Vectors.Vector(1..n);
begin
if l > k
then --put_line("Testing the following candidate face :");
--for ii in candidate'first..candidate'last-1 loop
-- put(pts(candidate(ii))); put(" & ");
--end loop;
--put(pts(candidate(candidate'last))); new_line;
if (k = m) or else Is_Lower_Face(k,n,m,elim,pts,ipvt,candidate,tol)
then Process(candidate,continue);
end if;
else for j in (i+1)..pts'last loop
candidate(l) := j;
pivots := ipvt;
-- put("Picking "); put(pts(j));
-- put(" Pivots : "); put(pivots); new_line;
Update_Eliminator(elim,l,pivots,pts(j).all,tol,piv);
-- put(" update of eliminator piv = "); put(piv,1);
-- put(" Pivots : "); put(pivots); new_line;
if (piv /= 0) and (piv /= n)
then Candidate_Lower_Faces(pivots,j,l+1);
end if;
Clear(elim(l));
exit when not continue;
end loop;
end if;
end Candidate_Lower_Faces;
begin
if k <= m
then declare
ipvt : Standard_Integer_Vectors.Vector(1..n);
begin
for i in ipvt'range loop
ipvt(i) := i;
end loop;
for i in pts'first..(pts'last-k) loop
candidate(0) := i;
elim(0) := pts(i); -- basis point
Candidate_Lower_Faces(ipvt,i,1);
end loop;
end;
end if;
end Enumerate_Lower_Faces;
procedure Enumerate_Faces_of_Sum
( ind,typ : in Standard_Integer_Vectors.Vector;
k : in natural;
pts : in Standard_Floating_VecVecs.VecVec;
tol : in double_float ) is
m : constant natural := pts'length; -- number of points
n : constant natural := pts(pts'first)'length; -- dimension of vectors
r : constant natural := ind'length; -- number of polytopes
candidates : Standard_Integer_VecVecs.VecVec(1..r);
elim : Standard_Floating_VecVecs.VecVec(1-r..k);
pivots : Standard_Integer_Vectors.Vector(1..n);
continue : boolean := true;
lasti,sum : natural;
procedure Candidate_Faces_of_Sum
( ipvt : in Standard_Integer_Vectors.Vector;
i : in integer ) is
-- DESCRIPTION :
-- Enumerates all faces of the given type on the sum,
-- i indicates the current polytope.
procedure Candidate_Faces ( ipvt : in Standard_Integer_Vectors.Vector;
start,l : in integer ) is
-- DESCRIPTION :
-- Enumerates all candidate k-faces, with k = typ(i).
-- The parameter l indicates the current element of pts to be chosen.
-- The previously chosen point is given by start.
piv,last : natural;
begin
if i = r
then last := m;
else last := ind(i+1)-1;
end if;
if l > typ(i)
then
if (typ(i) = last-ind(i)+1)
or else Is_Face_of_Sum
(sum,n,last-i+1,i,r,elim,pts(pts'first..last),tol,
candidates,ind,ipvt)
then
if i = r
then Process(candidates,continue);
else Candidate_Faces_of_Sum(ipvt,i+1);
end if;
end if;
else
for j in (start+1)..(last-typ(i)+l) loop
candidates(i)(l) := j;
if l = 0
then
Candidate_Faces(ipvt,j,l+1);
else
pivots := ipvt;
Update_Eliminator_for_Sum
(elim,sum-typ(i)+l,pivots,pts(j).all,i,tol,piv);
if piv /= 0
then Candidate_Faces(pivots,j,l+1);
end if;
Clear(elim(sum-typ(i)+l));
end if;
exit when not continue;
end loop;
end if;
end Candidate_Faces;
begin
candidates(i) := new Standard_Integer_Vectors.Vector(0..typ(i));
if i = r
then lasti := pts'last;
else lasti := ind(i+1)-1;
end if;
sum := sum + typ(i);
for j in ind(i)..(lasti-typ(i)) loop
candidates(i)(0) := j;
elim(1-i) := pts(j);
Candidate_Faces(ipvt,j,1);
end loop;
sum := sum - typ(i);
-- for j in (sum+1)..(sum+typ(i)) loop
-- Clear(elim(j));
-- end loop;
Clear(candidates(i));
end Candidate_Faces_of_Sum;
begin
declare
ipvt : Standard_Integer_Vectors.Vector(1..n);
begin
for i in ipvt'range loop
ipvt(i) := i;
end loop;
sum := 0;
Candidate_Faces_of_Sum(ipvt,1);
end;
end Enumerate_Faces_of_Sum;
end Floating_Face_Enumerators;
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to floating_face_enumerators.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports