File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift / vertices.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:29 2000 UTC (23 years, 10 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_Numbers; use Standard_Floating_Numbers;
with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
with Standard_Floating_Matrices; use Standard_Floating_Matrices;
with Dictionaries;
with Linear_Programming; use Linear_Programming;
with Integer_Face_Enumerators; use Integer_Face_Enumerators;
with Integer_Vectors_Utilities; use Integer_Vectors_Utilities;
with Lists_of_Vectors_Utilities; use Lists_of_Vectors_Utilities;
with Transformations; use Transformations;
with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
package body Vertices is
function Is_In_Hull ( point : Standard_Integer_Vectors.Vector; l : List )
return boolean is
fpt : Standard_Floating_Vectors.Vector(point'range);
begin
for i in point'range loop
fpt(i) := double_float(point(i));
end loop;
return Is_In_Hull(fpt,l);
end Is_In_Hull;
function Is_In_Hull ( point : Standard_Floating_Vectors.Vector; l : List )
return boolean is
-- ALGORITHM:
-- The following linear program will be solved:
--
-- min u_0 + u_1 + .. + u_n
--
-- l_1*p_1i + l_2*p_2i + .. + l_m*p_mi + u_i*q_i = q_i i=1,2,..,n
--
-- l_1 + l_2 + .. + l_m + u_0 = 1
--
-- to determine whether q belongs to the convex hull spanned by the
-- vectors p_1,p_2,..,p_m.
-- If all u_i are zero and all constraints are satisfied,
-- then q belongs to the convex hull.
-- CONSTANTS :
n : constant natural := point'last;
m : constant natural := 2*n+2; -- number of constraints
nb : constant natural := Length_Of(l)+n+1; -- number of unknowns
eps : constant double_float := 10.0**(-10);
-- VARIABLES :
dic : matrix(0..m,0..nb);
sol : Standard_Floating_Vectors.vector(1..nb);
inbas : Standard_Integer_Vectors.Vector(1..m) := (1..m => 0);
outbas : Standard_Integer_Vectors.Vector(1..nb) := (1..nb => 0);
nit : natural := 0;
feasi : boolean;
tmp : List;
pt : Link_to_Vector;
s : double_float;
begin
-- INITIALIZATION OF target :
for i in 0..(nb-n-1) loop
dic(0,i) := 0.0; -- sum of the lambda's
end loop;
for i in (nb-n)..nb loop
dic(0,i) := -1.0; -- sum of the mu's
end loop;
-- INITIALIZATION OF dic :
for i in 0..(nb-n) loop
dic(m-1,i) := 1.0; -- sum of the lambda's + mu_0
dic(m,i) := -1.0;
end loop;
for i in (nb-n+1)..nb loop
dic(m-1,i) := 0.0;
dic(m,i) := 0.0;
end loop;
for i in 1..n loop
for j in (nb-n)..nb loop
if i /= j-nb+n
then dic(i,j) := 0.0;
dic(i+n,j) := 0.0;
end if;
end loop;
end loop;
tmp := l;
for j in 1..(nb-n-1) loop
pt := Head_Of(tmp);
for i in pt'range loop
dic(i,j) := double_float(pt(i));
dic(i+n,j) := -double_float(pt(i));
end loop;
tmp := Tail_Of(tmp);
end loop;
for i in point'range loop
dic(i,0) := point(i);
dic(i+n,0) := -point(i);
dic(i,i+nb-n) := point(i);
dic(i+n,i+nb-n) := -point(i);
end loop;
-- SOLVE THE LINEAR PROGRAM :
Dictionaries.Init_Basis(inbas,outbas);
Dual_Simplex(dic,eps,inbas,outbas,nit,feasi);
if not feasi
then return false;
else
sol := Dictionaries.Primal_Solution(dic,inbas,outbas);
-- CHECK THE SOLUTION :
s := 0.0;
for i in 1..(nb-n-1) loop
s := s + sol(i);
end loop;
if abs(s - 1.0) > eps
then return false;
end if;
s := 0.0;
for i in (nb-n)..nb loop
s := s + sol(i);
end loop;
if abs(s) > eps
then return false;
end if;
return true;
end if;
end Is_In_Hull;
function Vertex_Points ( l : List ) return List is
result,result_last : List;
tmp,rest,origrest,rest_last : List;
pt : Link_to_Vector;
begin
if Is_Null(l) or else Is_Null(Tail_Of(l))
then return l;
else Copy(Tail_Of(l),rest);
origrest := rest;
rest_last := rest;
while not Is_Null(Tail_Of(rest_last)) loop
rest_last := Tail_Of(rest_last);
end loop;
tmp := l;
for i in 1..Length_Of(l) loop
pt := Head_Of(tmp);
if not Is_In_Hull(pt.all,rest)
then Append(result,result_last,pt.all);
Append(rest,rest_last,pt.all);
end if;
rest := Tail_Of(rest);
tmp := Tail_Of(tmp);
end loop;
Clear(origrest);
return result;
end if;
end Vertex_Points;
function Vertex_Points1 ( l : List ) return List is
len : constant natural := Length_Of(l);
pts : VecVec(1..len) := Shallow_Create(l);
result,result_last : List;
procedure Collect_Vertex ( i : in integer; cont : out boolean ) is
begin
Append(result,result_last,pts(i).all);
cont := true;
end Collect_Vertex;
procedure Enum_Vertices is new Enumerate_Vertices(Collect_Vertex);
begin
Enum_Vertices(pts);
return result;
end Vertex_Points1;
procedure Add ( pt : Link_to_Vector; l : in out List ) is
-- DESCRIPTION :
-- Constructs the point to the list, but without sharing.
newpt : Link_to_Vector := new Vector'(pt.all);
begin
Construct(newpt,l);
end Add;
function Extremal_Points ( l : List; v : Link_to_Vector ) return List is
min,max,sp : integer;
tmp,res : List;
minpt,maxpt,pt : Link_to_Vector;
begin
if Length_Of(l) <= 1
then Copy(l,res);
else pt := Head_Of(l);
min := pt*v; max := min;
minpt := pt; maxpt := pt;
tmp := Tail_Of(l);
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
sp := pt*v;
if sp > max
then max := sp; maxpt := pt;
elsif sp < min
then min := sp; minpt := pt;
end if;
tmp := Tail_Of(tmp);
end loop;
Add(minpt,res);
if min /= max
then Add(maxpt,res);
end if;
end if;
return res;
end Extremal_Points;
function Extremal_Points ( k,n : natural; exl,l : List ) return List is
res,tmp,nres : List;
iv : VecVec(1..k);
v : Link_to_Vector := new Vector(1..n);
sp : integer;
pt : Link_to_Vector;
done : boolean;
begin
tmp := exl;
for j in iv'range loop
iv(j) := Head_Of(tmp);
tmp := Tail_Of(tmp);
end loop;
v := Compute_Normal(iv);
sp := Head_Of(exl)*v;
nres := Extremal_Points(l,v);
tmp := nres;
res := exl;
done := false;
while not done and not Is_Null(tmp) loop
pt := Head_Of(tmp);
if pt*v /= sp
then Add(pt,res);
done := true;
else tmp := Tail_Of(tmp);
end if;
end loop;
Clear(v);
return res;
end Extremal_Points;
function Max_Extremal_Points ( k,n : natural; exl,l : List ) return List is
res,tmp,nres : List;
iv : VecVec(1..k);
v : Link_to_Vector := new Vector(1..n);
sp : integer;
pt : Link_to_Vector;
done : boolean;
begin
tmp := exl;
for j in iv'range loop
iv(j) := Head_Of(tmp);
tmp := Tail_Of(tmp);
end loop;
v := Compute_Normal(iv);
sp := Head_Of(exl)*v;
nres := Extremal_Points(l,v);
--put("The computed normal : "); put(v); new_line;
--put("The inner product : "); put(sp,1); new_line;
--put_line("The list of extremal points :"); put(nres);
tmp := nres;
Copy(exl,res);
done := false;
while not done and not Is_Null(tmp) loop
pt := Head_Of(tmp);
if pt*v /= sp
then Add(pt,res);
done := true;
else tmp := Tail_Of(tmp);
end if;
end loop;
Clear(nres);
if not done
then -- for all points x in l : <x,v> = sp
if n >= 2
then declare
i : integer := Pivot(v);
t,invt : Transfo;
texl,tl,tres : List;
begin
if i <= v'last
then t := Build_Transfo(v,i);
invt := Invert(t);
texl := Transform_and_Reduce(t,i,exl);
tl := Transform_and_Reduce(t,i,l);
tres := Max_Extremal_Points(k,n-1,texl,tl);
--put("The normal : "); put(v); new_line;
--put_line("The list of points : "); put(l);
--put_line("The list of extremal points :"); put(exl);
--put_line("The transformed list of points : "); put(tl);
--put_line("The transformed list of extremal points : ");
--put(texl);
--put_line("The computed transformed extremal points : ");
--put(tres);
--put("k : "); put(k,1); new_line;
Clear(t); Clear(texl); Clear(tl);
if Length_Of(tres) = k+1
then res := Insert_and_Transform(tres,i,sp,invt);
else Copy(exl,res);
end if;
--put_line("The computed extremal points : "); put(res);
Clear(tres); --responsible for segmentation violation...
Clear(invt);
else Copy(exl,res);
end if;
end;
else Copy(exl,res);
end if;
end if;
Clear(v);
--put_line("The new list of extremal points : "); put(res);
return res;
end Max_Extremal_Points;
function Extremal_Points ( n : natural; l : List ) return List is
res : List;
nor : Link_to_Vector := new Vector'(1..n => 1);
k : natural;
begin
res := Extremal_Points(l,nor);
Clear(nor);
if Length_Of(res) < 2
then return res;
else k := 2;
while k < n+1 loop
res := Extremal_Points(k,n,res,l);
exit when Length_Of(res) = k;
k := k+1;
end loop;
return res;
end if;
end Extremal_Points;
function Two_Extremal_Points ( n : natural; l : List ) return List is
-- DESCRIPTION :
-- Computes two extremal points of the list l.
res,tmp : List;
pt,minpt,maxpt : Link_to_Vector;
min,max : integer;
begin
if Length_Of(l) <= 2
then Copy(l,res);
else for i in 1..n loop
pt := Head_Of(l);
max := pt(i); min := max;
maxpt := pt; minpt := pt;
tmp := Tail_Of(l);
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
if pt(i) < min
then min := pt(i); minpt := pt;
elsif pt(i) > max
then max := pt(i); maxpt := pt;
end if;
tmp := Tail_Of(tmp);
end loop;
if Is_Null(res)
then Add(minpt,res);
if min /= max
then Add(maxpt,res);
end if;
else pt := Head_Of(res);
if pt.all /= minpt.all
then Add(minpt,res);
elsif pt.all /= maxpt.all
then Add(maxpt,res);
end if;
end if;
exit when (Length_Of(res) = 2);
end loop;
end if;
return res;
end Two_Extremal_Points;
function Max_Extremal_Points ( n : natural; l : List ) return List is
res,wl,wres,newres : List;
k : natural;
nullvec,shiftvec : Link_to_Vector;
-- shifted : boolean;
begin
if Length_Of(l) <= 2
then Copy(l,res);
else res := Two_Extremal_Points(n,l);
k := Length_Of(res);
if k = 2
then --nullvec := new Vector'(1..n => 0);
--if Is_In(nullvec,res)
-- then Copy(l,wl);
-- Copy(res,wres);
-- shifted := false;
-- else shiftvec := Head_Of(res);
-- wl := Shift(shiftvec,l);
-- wres := Shift(shiftvec,res);
-- shifted := true;
--end if;
--Clear(nullvec);
Copy(res,wres);
while k < n+1 loop
newres := Max_Extremal_Points(k,n,wres,l);
Copy(newres,wres);
exit when Length_Of(wres) = k;
k := k+1;
end loop;
res := wres;
-- Clear(res);
-- if shifted
-- then Min_Vector(shiftvec);
-- res := Shift(shiftvec,wres);
-- Clear(wres);
-- else Copy(wres,res);
-- end if;
-- Clear(wl);
end if;
end if;
return res;
end Max_Extremal_Points;
end Vertices;