File: [local] / OpenXM_contrib / PHC / Ada / Continuation / valipoco.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:22 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 integer_io; use integer_io;
with Timing_Package; use Timing_Package;
with File_Scanning; use File_Scanning;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
with Standard_Integer_Vectors;
with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
with Standard_Integer_VecVecs;
with Standard_Floating_Vectors;
with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
with Standard_Floating_VecVecs;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
procedure valipoco ( pocofile,resultfile : in file_type ) is
timer : Timing_Widget;
procedure put ( file : in file_type;
v : in Standard_Floating_Vectors.Vector;
fore,aft,exp : in natural ) is
begin
for i in v'range loop
put(file,v(i),fore,aft,exp);
end loop;
end put;
-- SCANNING THE INFORMATION FROM FILE :
procedure Scan_System
( file : in file_type; p : out Link_to_Poly_Sys ) is
-- DESCRIPTION :
-- Scans the file for the target system.
found : boolean := false;
lp : Link_to_Poly_Sys;
begin
Scan_and_Skip(file,"TARGET SYSTEM",found);
if found
then get(file,lp); p := lp;
end if;
end Scan_System;
procedure Scan_Dimensions
( file : in file_type; n,npaths : out natural ) is
-- DESCRIPTION :
-- Scans the input file for the dimension and the number of paths.
found : boolean := false;
tmp : integer := 0;
begin
Scan_and_Skip(file,"START SOLUTIONS",found);
if not found
then n := 0; npaths := 0;
else get(file,tmp); npaths := tmp;
get(file,tmp); n := tmp;
end if;
end Scan_Dimensions;
procedure Scan_Data
( infile,outfile : in file_type; n,npaths : in natural;
nv : out Standard_Floating_VecVecs.VecVec;
em : out Standard_Integer_Vectors.Vector;
ev,rv : out Standard_Floating_Vectors.Vector ) is
-- DESCRIPTION :
-- Scans for the computed directions nv, the errors ev,
-- the residuals rv and the estimated multiplicities em.
-- REQUIRED : nv'range = em'range = ev'range = rv'range = 1..npaths.
-- ON ENTRY :
-- infile file with input data from poco;
-- outfile file to write results to;
-- n dimension of the vectors along the paths;
-- npaths number of paths.
-- ON RETURN :
-- nv the computed path directions;
-- em estimated multiplicities;
-- rv residuals.
v : Standard_Floating_VecVecs.VecVec(1..npaths);
m : Standard_Integer_Vectors.Vector(1..npaths);
e,r : Standard_Floating_Vectors.Vector(1..npaths);
found : boolean := false;
begin
for i in v'range loop -- scan for normals
Scan_and_Skip(infile,"Computed direction",found);
exit when not found;
v(i) := new Standard_Floating_Vectors.Vector(1..n);
get(infile,v(i).all);
Scan(infile,"error :",found);
exit when not found;
get(infile,e(i));
end loop;
for i in m'range loop -- scan for normals and residuals
Scan(infile,"m :",found);
exit when not found;
get(infile,m(i));
Scan(infile,"res :",found);
exit when not found;
get(infile,r(i));
end loop;
put_line(outfile,
"COMPUTED DIRECTIONS, ESTIMATED MULTIPLICITY AND RESIDUALS : ");
for i in v'range loop
put(outfile,i,1); put(outfile," :"); put(outfile,v(i).all,3,3,3);
put(outfile," err : "); put(outfile,e(i),3,3,3);
put(outfile," m : "); put(outfile,m(i),1);
put(outfile," res : "); put(outfile,r(i),2,3,3);
new_line(outfile);
end loop;
nv := v; em := m; ev := e; rv := r;
end Scan_Data;
-- PROCESSING THE DATA :
function Lattice_Error ( x : Standard_Floating_Vectors.Vector )
return double_float is
-- DESCRIPTION :
-- Returns the sum of the differences of the components in x with
-- the closest corresponding integer.
res : double_float := 0.0;
begin
for i in x'range loop
res := res + abs(x(i) - double_float(integer(x(i))));
end loop;
return res;
end Lattice_Error;
function Lattice_Errors ( v : Standard_Floating_VecVecs.VecVec )
return Standard_Floating_Vectors.Vector is
res : Standard_Floating_Vectors.Vector(v'range);
begin
for i in v'range loop
res(i) := Lattice_Error(v(i).all);
end loop;
return res;
end Lattice_Errors;
function Maximum ( v : Standard_Floating_Vectors.Vector )
return double_float is
-- DESCRIPTION :
-- Returns the largest component in absolute value.
max : double_float := abs(v(v'first));
begin
for i in v'first+1..v'last loop
if abs(v(i)) > max
then max := abs(v(i));
end if;
end loop;
return max;
end Maximum;
procedure Scale ( v : in out Standard_Floating_VecVecs.VecVec;
e : in out Standard_Floating_Vectors.Vector ) is
-- DESCRIPTION :
-- Divides every vector by its largest element in absolute value,
-- with as well the corresponding errors.
use Standard_Floating_Vectors;
tol : double_float := 10.0**(-8);
max : double_float;
begin
for i in v'range loop
max := Maximum(v(i).all);
if max > tol
then v(i).all := (1.0/max)*v(i).all;
e(i) := e(i)/max;
end if;
end loop;
end Scale;
-- COMPUTE FREQUENCY TABLE :
function Equal ( v1,v2 : Standard_Floating_Vectors.Vector;
tol : double_float ) return boolean is
-- DESCRIPTION :
-- Returns true if abs(v1(i)-v2(i)) <= tol, for i in v1'range=v2'range.
begin
for i in v1'range loop
if abs(v1(i)-v2(i)) > tol
then return false;
end if;
end loop;
return true;
end Equal;
procedure Update_Frequency
( vpath : in out Standard_Integer_VecVecs.VecVec;
i,pathno : in natural ) is
-- DESCRIPTION :
-- A path, with pathno, has as path direction i.
freq : Standard_Integer_Vectors.Vector(vpath(i)'first..vpath(i)'last+1);
begin
freq(vpath(i)'range) := vpath(i).all;
freq(freq'last) := pathno;
Standard_Integer_Vectors.Clear(vpath(i));
vpath(i) := new Standard_Integer_Vectors.Vector'(freq);
end Update_Frequency;
procedure Update_Frequency_Table
( freqv : in out Standard_Floating_VecVecs.VecVec;
vpath : in out Standard_Integer_VecVecs.VecVec;
cnt : in out natural; tol : in double_float;
v : in Standard_Floating_Vectors.Vector;
pathno : in natural ) is
-- DESCRIPTION :
-- Scans the current frequency table for the vector v.
done : boolean := false;
begin
for i in 1..cnt loop
if Equal(v,freqv(i).all,tol)
then Update_Frequency(vpath,i,pathno); done := true;
end if;
exit when done;
end loop;
if not done
then cnt := cnt + 1;
freqv(cnt) := new Standard_Floating_Vectors.Vector'(v);
vpath(cnt) := new Standard_Integer_Vectors.Vector'(1..1 => pathno);
end if;
end Update_Frequency_Table;
procedure Frequency_Table
( v : in Standard_Floating_VecVecs.VecVec;
e,r : in Standard_Floating_Vectors.Vector;
tol : in double_float;
freqv : out Standard_Floating_VecVecs.Link_to_VecVec;
vpath : out Standard_Integer_VecVecs.Link_to_VecVec ) is
-- DESCRIPTION :
-- Computes the frequency table of path directions.
-- Only those directions for which the error is smaller than tol
-- and the corresponding residual is higher than tol are considered.
-- ON ENTRY :
-- v scaled computed path directions;
-- e errors on the computed path directions;
-- r residuals of the solutions at the end of the paths;
-- tol tolerance for precision.
-- ON RETURN :
-- freqv different path directions;
-- vpath for each direction, vector of paths with same direction.
cnt : natural := 0;
freqdirs : Standard_Floating_VecVecs.VecVec(v'range);
freqpath : Standard_Integer_VecVecs.VecVec(v'range);
begin
for i in v'range loop
if r(i) > tol and then e(i) < tol
then Update_Frequency_Table(freqdirs,freqpath,cnt,tol,v(i).all,i);
end if;
end loop;
freqv := new Standard_Floating_VecVecs.VecVec'(freqdirs(1..cnt));
vpath := new Standard_Integer_VecVecs.VecVec'(freqpath(1..cnt));
end Frequency_Table;
function Average_Error ( pathdir : Standard_Integer_Vectors.Vector;
e : Standard_Floating_Vectors.Vector )
return double_float is
-- DESCRIPTION :
-- Returns the average error over all paths in pathdir.
cnt : constant natural := pathdir'length;
res : double_float := 0.0;
begin
for i in pathdir'range loop
res := res + e(pathdir(i));
end loop;
res := res/double_float(cnt);
return res;
end Average_Error;
procedure Write_Frequency_Table
( file : in file_type;
freqv : in Standard_Floating_VecVecs.VecVec;
vpath : in Standard_Integer_VecVecs.VecVec;
m : in Standard_Integer_Vectors.Vector;
e : in Standard_Floating_Vectors.Vector ) is
-- DESCRIPTION :
-- Writes the frequency table on file, together with estimated m and error.
begin
put_line(file,"FREQUENCY TABLE OF PATH DIRECTIONS :");
for i in freqv'range loop
put(file,"v("); put(file,i,1); put(file,") : ");
put(file,freqv(i).all,3,3,3);
put(file," m : "); put(file,m(vpath(i)(1)),1);
put(file," err : "); put(file,e(vpath(i)(1)),3,3,3);
put(file," avgerr : "); put(file,Average_Error(vpath(i).all,e),3,3,3);
put(file," laterr : "); put(file,Lattice_Error(freqv(i).all),3,3,3);
new_line(file);
put(file," with "); put(file,vpath(i)'last,1);
put(file," paths : "); put(file,vpath(i)); new_line(file);
end loop;
end Write_Frequency_Table;
-- COMPUTING FACES OF POLYNOMIAL SYSTEMS :
function Product ( d : Degrees; v : Standard_Floating_Vectors.Vector )
return double_float is
-- DESCRIPTION : Returns <d,v>.
res : double_float := 0.0;
begin
for i in d'range loop
res := res + double_float(d(i))*v(i);
end loop;
return res;
end Product;
function Minimal_Support ( p : Poly; v : Standard_Floating_Vectors.Vector )
return double_float is
min : double_float := 0.0;
first : boolean := true;
procedure Scan_Term ( t : in Term; cont : out boolean ) is
ip : double_float := Product(t.dg,v);
begin
if first
then min := ip; first := false;
elsif ip < min
then min := ip;
end if;
cont := true;
end Scan_Term;
procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
begin
Scan_Terms(p);
return min;
end Minimal_Support;
function Face ( p : Poly; tol,ip : double_float;
v : Standard_Floating_Vectors.Vector ) return Poly is
-- DESCRIPTION :
-- Returns those terms t for which abs(<t.dg,v> - ip) <= tol.
res : Poly := Null_Poly;
procedure Visit_Term ( t : in Term; cont : out boolean ) is
begin
if abs(Product(t.dg,v) - ip) <= tol
then Add(res,t);
end if;
cont := true;
end Visit_Term;
procedure Visit_Terms is new Visiting_Iterator(Visit_Term);
begin
Visit_Terms(p);
return res;
end Face;
function Face ( p : Poly; tol : double_float;
v : Standard_Floating_Vectors.Vector )
return Poly is
-- DESCRIPTION :
-- Returns the face of the polynomial with outer normal v.
ip : double_float := Minimal_Support(p,v);
begin
return Face(p,tol,ip,v);
end Face;
function Face ( p : Poly_Sys; tol : double_float;
v : Standard_Floating_Vectors.Vector )
return Poly_Sys is
-- DESCRIPTION :
-- Returns the face of the system with outer normal v.
res : Poly_Sys(p'range);
begin
for i in p'range loop
res(i) := Face(p(i),tol,v);
end loop;
return res;
end Face;
procedure Face_Systems
( file : in file_type; p : in Poly_Sys; tol : in double_float;
v : in Standard_Floating_VecVecs.VecVec ) is
-- DESCRIPTION :
-- Computes the faces and the number of paths that diverged to it.
begin
put_line(file,"FACE SYSTEMS :");
for i in v'range loop
put(file,"path direction ");
put(file,i,1); put_line(file," :");
put(file,Face(p,tol,v(i).all));
end loop;
end Face_Systems;
-- MAIN PROCEDURE :
procedure Main ( infile,outfile : in file_type ) is
tol : constant double_float := 10.0**(-2);
lp : Link_to_Poly_Sys;
n,npaths : natural;
begin
Scan_System(infile,lp);
put_line(outfile,"TARGET SYSTEM : "); put(outfile,lp.all);
Scan_Dimensions(infile,n,npaths);
put(outfile," n : "); put(outfile,n,1);
put(outfile," and #paths : "); put(outfile,npaths,1); new_line(outfile);
declare
v : Standard_Floating_VecVecs.VecVec(1..npaths);
m : Standard_Integer_Vectors.Vector(1..npaths);
e,r,le : Standard_Floating_Vectors.Vector(1..npaths);
freqv : Standard_Floating_VecVecs.Link_to_VecVec;
vpath : Standard_Integer_VecVecs.Link_to_VecVec;
begin
Scan_Data(infile,outfile,n,npaths,v,m,e,r);
le := Lattice_Errors(v);
-- Scale(v,e); -- avoid because, then fractions may arise
Frequency_Table(v,e,r,tol,freqv,vpath);
Write_Frequency_Table(outfile,freqv.all,vpath.all,m,e);
Face_Systems(outfile,lp.all,tol,freqv.all);
end;
end Main;
begin
tstart(timer);
Main(pocofile,resultfile);
tstop(timer);
new_line(resultfile);
print_times(resultfile,timer,"Validation stage of polyhedral end game");
new_line(resultfile);
end valipoco;