with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
with Standard_Integer_Matrices; use Standard_Integer_Matrices;
with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
with Vertices; use Vertices;
with Frequency_Graph; use Frequency_Graph;
procedure Initial_Mixed_Cell
( n : in natural; mix : in Vector; pts : in Array_of_Lists;
mic : out Mixed_Cell; rest : in out Array_of_Lists ) is
res : Mixed_Cell;
acc,tmp,rest_last : List;
pt : Link_to_Vector;
expts : Array_of_Lists(pts'range);
perm : Vector(pts'range);
mind,minl,dims,lent,index : integer;
nullvec : constant Vector := (1..n => 0);
shiftvecs : VecVec(pts'range);
grap : Matrix(1..n,pts'range);
fail : boolean := false;
rowcnt,cnt : natural := 0;
mat : Matrix(1..n,1..n+1);
ipvt : Vector(1..n+1);
-- AUXILIAIRIES
procedure Add ( pt,shift : in Vector; l : in out List ) is
newpt : Link_to_Vector := new Vector'(pt);
begin
Add(newpt.all,shift);
Construct(newpt,l);
end Add;
procedure Add ( pt : in Vector; l : in out List ) is
newpt : Link_to_Vector := new Vector'(pt);
begin
Construct(newpt,l);
end Add;
procedure Fill ( l : in List; rowcnt : in out natural; m : in out Matrix ) is
-- DESCRIPTION :
-- Fills the matrix with the elements in l.
-- ON ENTRY :
-- l list of point vectors;
-- rowcnt m(m'first(1)..rowcnt) is already defined;
-- m matrix with m'range(2) = range of points in l.
-- rowcnt updated row counter;
-- m matrix with the nonzero points of l
tmp : List := l;
pt : Link_to_Vector;
begin
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
declare
nullvec : constant Vector(pt'range) := (pt'range => 0);
begin
if pt.all /= nullvec
then rowcnt := rowcnt + 1;
for k in pt'range loop
m(rowcnt,k) := pt(k);
end loop;
end if;
tmp := Tail_Of(tmp);
end;
end loop;
end Fill;
procedure Initialize ( l : in List; rowcnt : out natural;
m : in out Matrix; ipvt : in out Vector ) is
-- DESCRIPTION :
-- Given a list of linearly independent points, eventually with
-- 0 in l, a matrix will be constructed which contains these points
-- in upper triangular form.
cnt : natural := 0;
piv : natural;
begin
Fill(l,cnt,m);
for i in 1..cnt loop
m(i,m'last(2)) := i;
end loop;
for k in 1..cnt loop
Triangulate(k,m,ipvt,piv);
end loop;
rowcnt := cnt;
end Initialize;
procedure Linearly_Independent
( m : in out Matrix; rowcnt : in natural; ipvt : in Vector;
l : in List; len : out natural;
indep,indep_last : out List ) is
-- DESCRIPTION :
-- Computes the list of points which are linearly independent w.r.t.
-- the matrix m.
-- ON ENTRY :
-- m m(1..rowcnt,m'range(2)) is upper triangular,
-- can be used as work space;
-- rowcnt counter for the number of rows in m:
-- ipvt vector with the pivoting information;
-- l list of points to consider.
-- ON RETURN :
-- len length of the list indep;
-- indep the list of linearly independent points;
-- indep_last pointer to the last element of indep.
tmp,res,res_last : List;
pt : Link_to_Vector;
wipvt : Vector(ipvt'range) := ipvt;
piv,cnt : natural;
begin
if rowcnt < m'last(2)
then tmp := l; cnt := 0;
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
for i in pt'range loop
m(rowcnt+1,i) := pt(i);
end loop;
m(rowcnt+1,m'last(2)) := rowcnt+1;
Triangulate(rowcnt+1,m,wipvt,piv);
if m(rowcnt+1,rowcnt+1) /= 0
then cnt := cnt + 1;
Append(res,res_last,pt.all);
end if;
wipvt := ipvt;
tmp := Tail_Of(tmp);
end loop;
end if;
len := cnt;
indep := res; indep_last := res_last;
end Linearly_Independent;
procedure Construct_Basis ( l : in List; rowcnt : in out natural;
m : in out Matrix; ipvt : in out Vector ) is
-- DESCRIPTION :
-- Constructs a triangular basis for the vectors in l.
-- REQUIRED : dimensions of the vectors must match.
-- ON ENTRY :
-- l list of vector;
-- m triangular basis, stored in the rows of the matrix;
-- rowcnt index to the last significant row in m, could be 0,
-- when there is no initial basis to take into account;
-- ipvt initial pivoting vector, must be equal to the identity
-- permutation vector, when rowcnt = 0.
-- ON RETURN :
-- m triangular basis, stored in the rows of the matrix;
-- rowcnt index to the last significant row in m;
-- ipvt contains pivoting information.
tmp : List := l;
wipvt : Vector(ipvt'range);
pt : Link_to_Vector;
piv : natural;
begin
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
for i in pt'range loop
m(rowcnt+1,i) := pt(i);
end loop;
wipvt := ipvt;
m(rowcnt+1,m'last(2)) := rowcnt+1;
Triangulate(rowcnt+1,m,wipvt,piv);
if m(rowcnt+1,rowcnt+1) /= 0
then rowcnt := rowcnt + 1;
ipvt := wipvt;
end if;
exit when (rowcnt >= m'last(2));
tmp := Tail_Of(tmp);
end loop;
end Construct_Basis;
function Linearly_Independent ( l : List; x : Vector ) return boolean is
-- DESCRIPTION :
-- Returns true if the given vector x is linearly independent w.r.t.
-- the vectors in l.
basis : Matrix(1..Length_Of(l)+1,x'first..x'last+1);
ipvt : Vector(basis'range(2));
rowcnt : natural := 0;
piv : natural;
begin
for i in ipvt'range loop
ipvt(i) := i;
for j in basis'range(1) loop
basis(j,i) := 0;
end loop;
end loop;
Construct_Basis(l,rowcnt,basis,ipvt);
for i in x'range loop
basis(rowcnt+1,i) := x(i);
end loop;
basis(rowcnt+1,basis'last(2)) := rowcnt+1;
Triangulate(rowcnt+1,basis,ipvt,piv);
return (basis(rowcnt+1,rowcnt+1) /= 0);
end Linearly_Independent;
procedure Linearly_Independent
( l,xl : in List; indep,indep_last : in out List ) is
-- DESCRIPTION :
-- Returns those points in xl that are linearly independent from the
-- points in l.
begin
if not Is_Null(xl)
then
declare
x : Link_to_Vector := Head_Of(xl);
basis : Matrix(1..Length_Of(l)+1,x'first..x'last+1);
ipvt : Vector(basis'range(2));
rowcnt,len : natural := 0;
begin
for i in ipvt'range loop
ipvt(i) := i;
for j in basis'range(1) loop
basis(j,i) := 0;
end loop;
end loop;
Construct_Basis(l,rowcnt,basis,ipvt);
Linearly_Independent(basis,rowcnt,ipvt,xl,len,indep,indep_last);
end;
end if;
end Linearly_Independent;
function Construct_Rest ( l : Array_of_Lists; start : natural )
return List is
-- DESCRIPTION :
-- Collects all remaining points of l(i) into one single list,
-- for i in start..l'last.
tmp,res,res_last : List;
pt : Link_to_Vector;
begin
for i in start..l'last loop
tmp := l(i);
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
if not Is_In(res,pt)
then Append(res,res_last,pt.all);
end if;
tmp := Tail_Of(tmp);
end loop;
end loop;
return res;
end Construct_Rest;
procedure Sort_Co_Linearly_Independent
( l : in out List; al : in Array_of_Lists;
start : in natural ) is
-- DESCRIPTION :
-- Sorts the points in l, by putting those points that are linearly
-- independent w.r.t. the rest of the poins in al, in front of the list.
rest : List := Construct_Rest(al,start);
tmp,indep,indep_last : List;
pt : Link_to_Vector;
begin
-- put_line("The set of points in al just after Construct_Rest :"); put(al);
-- put_line("The rest of the points : "); put(rest);
Linearly_Independent(rest,l,indep,indep_last);
if Is_Null(indep)
then null; -- nothing to put in front of l
else tmp := l; -- append the other points of l to indep
-- put_line("Linearly co-independent points : "); put(indep);
-- put_line(" w.r.t. the rest : "); put(rest);
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
if not Is_In(indep,pt)
then Append(indep,indep_last,pt.all);
end if;
tmp := Tail_Of(tmp);
end loop;
Copy(indep,l); Deep_Clear(indep);
end if;
Deep_Clear(rest);
-- put_line("The list of points in al after Sort_Co : "); put(al);
end Sort_Co_Linearly_Independent;
procedure Incremental_Dimension
( m : in out Matrix; rowcnt : in natural; ipvt : in Vector;
l : in List; dim : out natural ) is
-- DESCRIPTION :
-- Computes the number of points which are linearly independent w.r.t.
-- the matrix m.
-- ON ENTRY :
-- m m(1..rowcnt,m'range(2)) is upper triangular,
-- can be used as work space;
-- rowcnt counter for the number of rows in m:
-- ipvt vector with the pivoting information;
-- l list of points to consider.
-- ON RETURN :
-- dim the number of linearly independent points.
tmp : List;
pt : Link_to_Vector;
cnt,piv : natural;
newipvt,wipvt : Vector(ipvt'range);
begin
wipvt := ipvt; newipvt := ipvt;
tmp := l; cnt := 0;
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
cnt := cnt + 1;
for i in pt'range loop
m(rowcnt+cnt,i) := pt(i);
end loop;
m(rowcnt+cnt,m'last(2)) := rowcnt+cnt;
Triangulate(rowcnt+cnt,m,newipvt,piv);
if m(rowcnt+cnt,rowcnt+cnt) /= 0
then wipvt := newipvt;
else newipvt := wipvt;
cnt := cnt - 1;
end if;
tmp := Tail_Of(tmp);
exit when ((rowcnt + cnt) >= m'last(1));
end loop;
dim := cnt;
end Incremental_Dimension;
procedure Incremental_Dimension
( m : in out Matrix; rowcnt : in natural; ipvt : in Vector;
l : in out List; dim,len : out natural ) is
-- DESCRIPTION :
-- Computes the increase in dimension by considering the points
-- in the list l, w.r.t. the matrix m.
-- ON ENTRY :
-- m m(1..rowcnt,m'range(2)) is upper triangular,
-- rowcnt counter for the number of rows in m:
-- ipvt vector with the pivoting information;
-- l list of points to consider.
-- ON RETURN :
-- l list of linearly independent points w.r.t. m;
-- dim increase in dimension;
-- len length of the list of linearly independent points
-- w.r.t. the rows in the matrix m.
-- ALGORITHM :
-- works in two stages:
-- 1. determination of all linearly independent points;
-- 2. determination of the dimension.
workm : Matrix(m'range(1),m'range(2));
indep,indep_last : List;
begin
-- Initialize workm :
-- put_line("The list to investigate : "); put(l);
for i in 1..rowcnt loop
for j in m'range(2) loop
workm(i,j) := m(i,j);
end loop;
end loop;
-- Determine linearly independent points :
Linearly_Independent(workm,rowcnt,ipvt,l,len,indep,indep_last);
-- Determine the incremental dimension :
Incremental_Dimension(workm,rowcnt,ipvt,indep,dim);
Copy(indep,l); Deep_Clear(indep);
end Incremental_Dimension;
procedure Next_Point ( acc,l : in List; pt : out Link_to_Vector;
fail : out boolean; rowcnt : in out natural;
m : in out Matrix; ipvt : in out Vector ) is
-- DESCRIPTION :
-- A new point out of l, not in the list acc will be chosen.
-- The point pt has to be linearly independent from the other,
-- already chosen points. Therefore, the upper triangular matrix m
-- will be used and properly updated.
res : Link_to_Vector;
tmp : List := l;
newipvt : Vector(ipvt'range) := ipvt;
done : boolean := false;
piv : natural;
begin
while not Is_Null(tmp) loop
res := Head_Of(tmp);
if not Is_In(acc,res)
then --put("Checking point "); put(res); new_line;
rowcnt := rowcnt + 1;
for i in res'range loop
m(rowcnt,i) := res(i);
end loop;
m(rowcnt,m'last(2)) := rowcnt;
Triangulate(rowcnt,m,newipvt,piv);
if m(rowcnt,rowcnt) = 0
then rowcnt := rowcnt - 1;
newipvt := ipvt;
else ipvt := newipvt;
pt := res;
done := true;
end if;
end if;
exit when done;
tmp := Tail_Of(tmp);
end loop;
fail := not done;
end Next_Point;
begin
-- INITIALIZE : compute for each list the basis points
res.nor := new Vector'(1..n+1 => 0); res.nor(n+1) := 1;
res.pts := new Array_of_Lists(mix'range);
for i in pts'range loop
expts(i) := Max_Extremal_Points(n,pts(i)); perm(i) := i;
-- put("Extremal points for component "); put(i,1); put_line(" :");
-- put(expts(i));
end loop;
-- INITIALIZE : order the lists according to occurencies
grap := Graph(n,expts);
-- put_line("The graph matrix of the extremal points : "); put(grap);
for i in expts'range loop
Sort(expts(i),grap);
-- put("Ordered extremal points for component "); put(i,1); put_line(" :");
-- put(expts(i));
end loop;
-- INITIALIZE : choose one anchor point for each component,
-- shift when necessary :
for i in expts'range loop
if Is_In(pts(i),nullvec)
then Add(nullvec,res.pts(i));
shiftvecs(i) := null;
else -- ADD LAST POINT = LEAST IMPORTANT
tmp := expts(i);
if not Is_Null(tmp)
then
while not Is_Null(Tail_Of(tmp)) loop
tmp := Tail_Of(tmp);
end loop;
pt := Head_Of(tmp);
Add(pt.all,res.pts(i));
shiftvecs(i) := new Vector'(pt.all);
-- put("The shift vector : "); put(shiftvecs(i)); new_line;
Shift(expts(i),shiftvecs(i));
-- put_line("The shifted list of points : "); put(expts(i));
end if;
end if;
end loop;
-- INITIALIZE : intermediate data structures mat and ipvt :
for i in ipvt'range loop
ipvt(i) := i;
end loop;
for i in mat'range(1) loop
for j in mat'range(2) loop
mat(i,j) := 0;
end loop;
end loop;
rowcnt := 0;
-- COMPUTE CELL : based on lists with extremal points
for i in expts'range loop -- MAIN LOOP
-- LOOK FOR SMALLEST SET :
index := i;
if i = 1
then mind := Length_Of(expts(i))-1;
for j in i+1..expts'last loop
dims := Length_Of(expts(j))-1;
if dims < mind
then index := j; mind := dims;
end if;
end loop;
else --put("Calling Incremental_Dimension for i = "); put(i,1); new_line;
Incremental_Dimension(mat,rowcnt,ipvt,expts(i),mind,minl);
--put(" dimension : "); put(mind,1);
--put(" length : " ); put(minl,1); new_line;
for j in i+1..expts'last loop
--put("Calling Incremental_Dimension for j = "); put(j,1); new_line;
Incremental_Dimension(mat,rowcnt,ipvt,expts(j),dims,lent);
--put(" dimension : "); put(dims,1);
--put(" length : " ); put(lent,1); new_line;
if (dims < mind) or else ((dims = mind) and then (lent < minl))
then index := j; mind := dims; minl := lent;
end if;
end loop;
end if;
-- put("The index : "); put(index,1); new_line;
-- put(" with minimal dimension : "); put(mind,1); new_line;
-- put("perm before permute : "); put(perm); new_line;
if index /= i
then tmp := expts(index); expts(index) := expts(i); expts(i) := tmp;
mind := perm(i); perm(i) := perm(index); perm(index) := mind;
end if;
-- SORT ACCORDING TO OCCURRENCIES AND TO CO-LINEARLY :
Sort(expts(i),grap,i-1,perm);
-- put_line("After first ordering : "); put(expts(i));
Sort_Co_Linearly_Independent(expts(i),expts,i+1);
-- put_line("The ordered list : "); put(expts(i));
-- put("perm after permute : "); put(perm); new_line;
-- CHOOSE THE POINTS :
if i = 1
then Add(nullvec,acc);
tmp := expts(i);
for j in 1..mix(perm(i)) loop
pt := Head_Of(tmp);
if pt.all = nullvec -- skip the null vector
then tmp := Tail_Of(tmp);
if not Is_Null(tmp)
then pt := Head_Of(tmp);
end if;
end if;
exit when Is_Null(tmp);
Add(pt.all,acc);
-- put_line("frequency graph before ignore :"); put(grap);
Ignore(grap,pt.all);
-- put_line("frequency graph after ignore :"); put(grap);
if shiftvecs(perm(i)) /= null
then Add(pt.all,shiftvecs(perm(i)).all,res.pts(perm(i)));
else Add(pt.all,res.pts(perm(i)));
end if;
tmp := Tail_Of(tmp);
exit when Is_Null(tmp);
end loop;
cnt := Length_Of(acc);
Initialize(acc,rowcnt,mat,ipvt);
-- put_line("The list acc : "); put(acc);
-- put_line("The matrix mat : "); put(mat,1,rowcnt);
-- put_line("The vector ipvt : "); put(ipvt); new_line;
fail := (cnt <= mix(perm(i)));
else for j in 1..mix(perm(i)) loop
Next_Point(acc,expts(i),pt,fail,rowcnt,mat,ipvt);
if not fail
then if shiftvecs(perm(i)) /= null
then Add(pt.all,shiftvecs(perm(i)).all,res.pts(perm(i)));
else Add(pt.all,res.pts(perm(i)));
end if;
Add(pt.all,acc); cnt := cnt + 1;
-- put_line("frequency graph before ignore :"); put(grap);
Ignore(grap,pt.all);
-- put_line("frequency graph after ignore :"); put(grap);
end if;
exit when fail;
end loop;
-- put_line("The list acc : "); put(acc);
-- put_line("The matrix mat : "); put(mat,1,rowcnt);
-- put_line("The vector ipvt : "); put(ipvt); new_line;
fail := (Length_Of(res.pts(perm(i))) <= mix(perm(i)));
end if;
exit when fail;
end loop; -- END OF MAIN LOOP
Deep_Clear(acc);
-- COMPUTE THE REST OF THE POINT LISTS :
if not fail
then for i in pts'range loop
tmp := pts(i);
rest_last := rest(i);
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
if not Is_In(res.pts(i),pt)
then Append(rest(i),rest_last,pt);
end if;
tmp := Tail_Of(tmp);
end loop;
end loop;
end if;
-- GIVE THE POINTS IN THE INITIAL SIMPLEX LIFTING VALUE 0 :
for i in res.pts'range loop
tmp := res.pts(i);
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
declare
lpt : Link_to_Vector;
begin
lpt := new Vector(1..n+1);
lpt(pt'range) := pt.all; lpt(n+1) := 0;
Set_Head(tmp,lpt);
end;
tmp := Tail_Of(tmp);
end loop;
end loop;
mic := res;
end Initial_Mixed_Cell;