Annotation of OpenXM_contrib/PHC/Ada/Continuation/valipoco.adb, Revision 1.1
1.1 ! maekawa 1: with integer_io; use integer_io;
! 2: with Timing_Package; use Timing_Package;
! 3: with File_Scanning; use File_Scanning;
! 4: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 5: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
! 6: with Standard_Integer_Vectors;
! 7: with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
! 8: with Standard_Integer_VecVecs;
! 9: with Standard_Floating_Vectors;
! 10: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
! 11: with Standard_Floating_VecVecs;
! 12: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 13: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
! 14: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
! 15:
! 16: procedure valipoco ( pocofile,resultfile : in file_type ) is
! 17:
! 18: timer : Timing_Widget;
! 19:
! 20: procedure put ( file : in file_type;
! 21: v : in Standard_Floating_Vectors.Vector;
! 22: fore,aft,exp : in natural ) is
! 23: begin
! 24: for i in v'range loop
! 25: put(file,v(i),fore,aft,exp);
! 26: end loop;
! 27: end put;
! 28:
! 29: -- SCANNING THE INFORMATION FROM FILE :
! 30:
! 31: procedure Scan_System
! 32: ( file : in file_type; p : out Link_to_Poly_Sys ) is
! 33:
! 34: -- DESCRIPTION :
! 35: -- Scans the file for the target system.
! 36:
! 37: found : boolean := false;
! 38: lp : Link_to_Poly_Sys;
! 39:
! 40: begin
! 41: Scan_and_Skip(file,"TARGET SYSTEM",found);
! 42: if found
! 43: then get(file,lp); p := lp;
! 44: end if;
! 45: end Scan_System;
! 46:
! 47: procedure Scan_Dimensions
! 48: ( file : in file_type; n,npaths : out natural ) is
! 49:
! 50: -- DESCRIPTION :
! 51: -- Scans the input file for the dimension and the number of paths.
! 52:
! 53: found : boolean := false;
! 54: tmp : integer := 0;
! 55:
! 56: begin
! 57: Scan_and_Skip(file,"START SOLUTIONS",found);
! 58: if not found
! 59: then n := 0; npaths := 0;
! 60: else get(file,tmp); npaths := tmp;
! 61: get(file,tmp); n := tmp;
! 62: end if;
! 63: end Scan_Dimensions;
! 64:
! 65: procedure Scan_Data
! 66: ( infile,outfile : in file_type; n,npaths : in natural;
! 67: nv : out Standard_Floating_VecVecs.VecVec;
! 68: em : out Standard_Integer_Vectors.Vector;
! 69: ev,rv : out Standard_Floating_Vectors.Vector ) is
! 70:
! 71: -- DESCRIPTION :
! 72: -- Scans for the computed directions nv, the errors ev,
! 73: -- the residuals rv and the estimated multiplicities em.
! 74:
! 75: -- REQUIRED : nv'range = em'range = ev'range = rv'range = 1..npaths.
! 76:
! 77: -- ON ENTRY :
! 78: -- infile file with input data from poco;
! 79: -- outfile file to write results to;
! 80: -- n dimension of the vectors along the paths;
! 81: -- npaths number of paths.
! 82:
! 83: -- ON RETURN :
! 84: -- nv the computed path directions;
! 85: -- em estimated multiplicities;
! 86: -- rv residuals.
! 87:
! 88: v : Standard_Floating_VecVecs.VecVec(1..npaths);
! 89: m : Standard_Integer_Vectors.Vector(1..npaths);
! 90: e,r : Standard_Floating_Vectors.Vector(1..npaths);
! 91: found : boolean := false;
! 92:
! 93: begin
! 94: for i in v'range loop -- scan for normals
! 95: Scan_and_Skip(infile,"Computed direction",found);
! 96: exit when not found;
! 97: v(i) := new Standard_Floating_Vectors.Vector(1..n);
! 98: get(infile,v(i).all);
! 99: Scan(infile,"error :",found);
! 100: exit when not found;
! 101: get(infile,e(i));
! 102: end loop;
! 103: for i in m'range loop -- scan for normals and residuals
! 104: Scan(infile,"m :",found);
! 105: exit when not found;
! 106: get(infile,m(i));
! 107: Scan(infile,"res :",found);
! 108: exit when not found;
! 109: get(infile,r(i));
! 110: end loop;
! 111: put_line(outfile,
! 112: "COMPUTED DIRECTIONS, ESTIMATED MULTIPLICITY AND RESIDUALS : ");
! 113: for i in v'range loop
! 114: put(outfile,i,1); put(outfile," :"); put(outfile,v(i).all,3,3,3);
! 115: put(outfile," err : "); put(outfile,e(i),3,3,3);
! 116: put(outfile," m : "); put(outfile,m(i),1);
! 117: put(outfile," res : "); put(outfile,r(i),2,3,3);
! 118: new_line(outfile);
! 119: end loop;
! 120: nv := v; em := m; ev := e; rv := r;
! 121: end Scan_Data;
! 122:
! 123: -- PROCESSING THE DATA :
! 124:
! 125: function Lattice_Error ( x : Standard_Floating_Vectors.Vector )
! 126: return double_float is
! 127:
! 128: -- DESCRIPTION :
! 129: -- Returns the sum of the differences of the components in x with
! 130: -- the closest corresponding integer.
! 131:
! 132: res : double_float := 0.0;
! 133:
! 134: begin
! 135: for i in x'range loop
! 136: res := res + abs(x(i) - double_float(integer(x(i))));
! 137: end loop;
! 138: return res;
! 139: end Lattice_Error;
! 140:
! 141: function Lattice_Errors ( v : Standard_Floating_VecVecs.VecVec )
! 142: return Standard_Floating_Vectors.Vector is
! 143:
! 144: res : Standard_Floating_Vectors.Vector(v'range);
! 145:
! 146: begin
! 147: for i in v'range loop
! 148: res(i) := Lattice_Error(v(i).all);
! 149: end loop;
! 150: return res;
! 151: end Lattice_Errors;
! 152:
! 153: function Maximum ( v : Standard_Floating_Vectors.Vector )
! 154: return double_float is
! 155:
! 156: -- DESCRIPTION :
! 157: -- Returns the largest component in absolute value.
! 158:
! 159: max : double_float := abs(v(v'first));
! 160:
! 161: begin
! 162: for i in v'first+1..v'last loop
! 163: if abs(v(i)) > max
! 164: then max := abs(v(i));
! 165: end if;
! 166: end loop;
! 167: return max;
! 168: end Maximum;
! 169:
! 170: procedure Scale ( v : in out Standard_Floating_VecVecs.VecVec;
! 171: e : in out Standard_Floating_Vectors.Vector ) is
! 172:
! 173: -- DESCRIPTION :
! 174: -- Divides every vector by its largest element in absolute value,
! 175: -- with as well the corresponding errors.
! 176:
! 177: use Standard_Floating_Vectors;
! 178: tol : double_float := 10.0**(-8);
! 179: max : double_float;
! 180:
! 181: begin
! 182: for i in v'range loop
! 183: max := Maximum(v(i).all);
! 184: if max > tol
! 185: then v(i).all := (1.0/max)*v(i).all;
! 186: e(i) := e(i)/max;
! 187: end if;
! 188: end loop;
! 189: end Scale;
! 190:
! 191: -- COMPUTE FREQUENCY TABLE :
! 192:
! 193: function Equal ( v1,v2 : Standard_Floating_Vectors.Vector;
! 194: tol : double_float ) return boolean is
! 195:
! 196: -- DESCRIPTION :
! 197: -- Returns true if abs(v1(i)-v2(i)) <= tol, for i in v1'range=v2'range.
! 198:
! 199: begin
! 200: for i in v1'range loop
! 201: if abs(v1(i)-v2(i)) > tol
! 202: then return false;
! 203: end if;
! 204: end loop;
! 205: return true;
! 206: end Equal;
! 207:
! 208: procedure Update_Frequency
! 209: ( vpath : in out Standard_Integer_VecVecs.VecVec;
! 210: i,pathno : in natural ) is
! 211:
! 212: -- DESCRIPTION :
! 213: -- A path, with pathno, has as path direction i.
! 214:
! 215: freq : Standard_Integer_Vectors.Vector(vpath(i)'first..vpath(i)'last+1);
! 216:
! 217: begin
! 218: freq(vpath(i)'range) := vpath(i).all;
! 219: freq(freq'last) := pathno;
! 220: Standard_Integer_Vectors.Clear(vpath(i));
! 221: vpath(i) := new Standard_Integer_Vectors.Vector'(freq);
! 222: end Update_Frequency;
! 223:
! 224: procedure Update_Frequency_Table
! 225: ( freqv : in out Standard_Floating_VecVecs.VecVec;
! 226: vpath : in out Standard_Integer_VecVecs.VecVec;
! 227: cnt : in out natural; tol : in double_float;
! 228: v : in Standard_Floating_Vectors.Vector;
! 229: pathno : in natural ) is
! 230:
! 231: -- DESCRIPTION :
! 232: -- Scans the current frequency table for the vector v.
! 233:
! 234: done : boolean := false;
! 235:
! 236: begin
! 237: for i in 1..cnt loop
! 238: if Equal(v,freqv(i).all,tol)
! 239: then Update_Frequency(vpath,i,pathno); done := true;
! 240: end if;
! 241: exit when done;
! 242: end loop;
! 243: if not done
! 244: then cnt := cnt + 1;
! 245: freqv(cnt) := new Standard_Floating_Vectors.Vector'(v);
! 246: vpath(cnt) := new Standard_Integer_Vectors.Vector'(1..1 => pathno);
! 247: end if;
! 248: end Update_Frequency_Table;
! 249:
! 250: procedure Frequency_Table
! 251: ( v : in Standard_Floating_VecVecs.VecVec;
! 252: e,r : in Standard_Floating_Vectors.Vector;
! 253: tol : in double_float;
! 254: freqv : out Standard_Floating_VecVecs.Link_to_VecVec;
! 255: vpath : out Standard_Integer_VecVecs.Link_to_VecVec ) is
! 256:
! 257: -- DESCRIPTION :
! 258: -- Computes the frequency table of path directions.
! 259: -- Only those directions for which the error is smaller than tol
! 260: -- and the corresponding residual is higher than tol are considered.
! 261:
! 262: -- ON ENTRY :
! 263: -- v scaled computed path directions;
! 264: -- e errors on the computed path directions;
! 265: -- r residuals of the solutions at the end of the paths;
! 266: -- tol tolerance for precision.
! 267:
! 268: -- ON RETURN :
! 269: -- freqv different path directions;
! 270: -- vpath for each direction, vector of paths with same direction.
! 271:
! 272: cnt : natural := 0;
! 273: freqdirs : Standard_Floating_VecVecs.VecVec(v'range);
! 274: freqpath : Standard_Integer_VecVecs.VecVec(v'range);
! 275:
! 276: begin
! 277: for i in v'range loop
! 278: if r(i) > tol and then e(i) < tol
! 279: then Update_Frequency_Table(freqdirs,freqpath,cnt,tol,v(i).all,i);
! 280: end if;
! 281: end loop;
! 282: freqv := new Standard_Floating_VecVecs.VecVec'(freqdirs(1..cnt));
! 283: vpath := new Standard_Integer_VecVecs.VecVec'(freqpath(1..cnt));
! 284: end Frequency_Table;
! 285:
! 286: function Average_Error ( pathdir : Standard_Integer_Vectors.Vector;
! 287: e : Standard_Floating_Vectors.Vector )
! 288: return double_float is
! 289:
! 290: -- DESCRIPTION :
! 291: -- Returns the average error over all paths in pathdir.
! 292:
! 293: cnt : constant natural := pathdir'length;
! 294: res : double_float := 0.0;
! 295:
! 296: begin
! 297: for i in pathdir'range loop
! 298: res := res + e(pathdir(i));
! 299: end loop;
! 300: res := res/double_float(cnt);
! 301: return res;
! 302: end Average_Error;
! 303:
! 304: procedure Write_Frequency_Table
! 305: ( file : in file_type;
! 306: freqv : in Standard_Floating_VecVecs.VecVec;
! 307: vpath : in Standard_Integer_VecVecs.VecVec;
! 308: m : in Standard_Integer_Vectors.Vector;
! 309: e : in Standard_Floating_Vectors.Vector ) is
! 310:
! 311: -- DESCRIPTION :
! 312: -- Writes the frequency table on file, together with estimated m and error.
! 313:
! 314: begin
! 315: put_line(file,"FREQUENCY TABLE OF PATH DIRECTIONS :");
! 316: for i in freqv'range loop
! 317: put(file,"v("); put(file,i,1); put(file,") : ");
! 318: put(file,freqv(i).all,3,3,3);
! 319: put(file," m : "); put(file,m(vpath(i)(1)),1);
! 320: put(file," err : "); put(file,e(vpath(i)(1)),3,3,3);
! 321: put(file," avgerr : "); put(file,Average_Error(vpath(i).all,e),3,3,3);
! 322: put(file," laterr : "); put(file,Lattice_Error(freqv(i).all),3,3,3);
! 323: new_line(file);
! 324: put(file," with "); put(file,vpath(i)'last,1);
! 325: put(file," paths : "); put(file,vpath(i)); new_line(file);
! 326: end loop;
! 327: end Write_Frequency_Table;
! 328:
! 329: -- COMPUTING FACES OF POLYNOMIAL SYSTEMS :
! 330:
! 331: function Product ( d : Degrees; v : Standard_Floating_Vectors.Vector )
! 332: return double_float is
! 333:
! 334: -- DESCRIPTION : Returns <d,v>.
! 335:
! 336: res : double_float := 0.0;
! 337:
! 338: begin
! 339: for i in d'range loop
! 340: res := res + double_float(d(i))*v(i);
! 341: end loop;
! 342: return res;
! 343: end Product;
! 344:
! 345: function Minimal_Support ( p : Poly; v : Standard_Floating_Vectors.Vector )
! 346: return double_float is
! 347:
! 348: min : double_float := 0.0;
! 349: first : boolean := true;
! 350:
! 351: procedure Scan_Term ( t : in Term; cont : out boolean ) is
! 352:
! 353: ip : double_float := Product(t.dg,v);
! 354:
! 355: begin
! 356: if first
! 357: then min := ip; first := false;
! 358: elsif ip < min
! 359: then min := ip;
! 360: end if;
! 361: cont := true;
! 362: end Scan_Term;
! 363: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
! 364:
! 365: begin
! 366: Scan_Terms(p);
! 367: return min;
! 368: end Minimal_Support;
! 369:
! 370: function Face ( p : Poly; tol,ip : double_float;
! 371: v : Standard_Floating_Vectors.Vector ) return Poly is
! 372:
! 373: -- DESCRIPTION :
! 374: -- Returns those terms t for which abs(<t.dg,v> - ip) <= tol.
! 375:
! 376: res : Poly := Null_Poly;
! 377:
! 378: procedure Visit_Term ( t : in Term; cont : out boolean ) is
! 379: begin
! 380: if abs(Product(t.dg,v) - ip) <= tol
! 381: then Add(res,t);
! 382: end if;
! 383: cont := true;
! 384: end Visit_Term;
! 385: procedure Visit_Terms is new Visiting_Iterator(Visit_Term);
! 386:
! 387: begin
! 388: Visit_Terms(p);
! 389: return res;
! 390: end Face;
! 391:
! 392: function Face ( p : Poly; tol : double_float;
! 393: v : Standard_Floating_Vectors.Vector )
! 394: return Poly is
! 395:
! 396: -- DESCRIPTION :
! 397: -- Returns the face of the polynomial with outer normal v.
! 398:
! 399: ip : double_float := Minimal_Support(p,v);
! 400:
! 401: begin
! 402: return Face(p,tol,ip,v);
! 403: end Face;
! 404:
! 405: function Face ( p : Poly_Sys; tol : double_float;
! 406: v : Standard_Floating_Vectors.Vector )
! 407: return Poly_Sys is
! 408:
! 409: -- DESCRIPTION :
! 410: -- Returns the face of the system with outer normal v.
! 411:
! 412: res : Poly_Sys(p'range);
! 413:
! 414: begin
! 415: for i in p'range loop
! 416: res(i) := Face(p(i),tol,v);
! 417: end loop;
! 418: return res;
! 419: end Face;
! 420:
! 421: procedure Face_Systems
! 422: ( file : in file_type; p : in Poly_Sys; tol : in double_float;
! 423: v : in Standard_Floating_VecVecs.VecVec ) is
! 424:
! 425: -- DESCRIPTION :
! 426: -- Computes the faces and the number of paths that diverged to it.
! 427:
! 428: begin
! 429: put_line(file,"FACE SYSTEMS :");
! 430: for i in v'range loop
! 431: put(file,"path direction ");
! 432: put(file,i,1); put_line(file," :");
! 433: put(file,Face(p,tol,v(i).all));
! 434: end loop;
! 435: end Face_Systems;
! 436:
! 437: -- MAIN PROCEDURE :
! 438:
! 439: procedure Main ( infile,outfile : in file_type ) is
! 440:
! 441: tol : constant double_float := 10.0**(-2);
! 442: lp : Link_to_Poly_Sys;
! 443: n,npaths : natural;
! 444:
! 445: begin
! 446: Scan_System(infile,lp);
! 447: put_line(outfile,"TARGET SYSTEM : "); put(outfile,lp.all);
! 448: Scan_Dimensions(infile,n,npaths);
! 449: put(outfile," n : "); put(outfile,n,1);
! 450: put(outfile," and #paths : "); put(outfile,npaths,1); new_line(outfile);
! 451: declare
! 452: v : Standard_Floating_VecVecs.VecVec(1..npaths);
! 453: m : Standard_Integer_Vectors.Vector(1..npaths);
! 454: e,r,le : Standard_Floating_Vectors.Vector(1..npaths);
! 455: freqv : Standard_Floating_VecVecs.Link_to_VecVec;
! 456: vpath : Standard_Integer_VecVecs.Link_to_VecVec;
! 457: begin
! 458: Scan_Data(infile,outfile,n,npaths,v,m,e,r);
! 459: le := Lattice_Errors(v);
! 460: -- Scale(v,e); -- avoid because, then fractions may arise
! 461: Frequency_Table(v,e,r,tol,freqv,vpath);
! 462: Write_Frequency_Table(outfile,freqv.all,vpath.all,m,e);
! 463: Face_Systems(outfile,lp.all,tol,freqv.all);
! 464: end;
! 465: end Main;
! 466:
! 467: begin
! 468: tstart(timer);
! 469: Main(pocofile,resultfile);
! 470: tstop(timer);
! 471: new_line(resultfile);
! 472: print_times(resultfile,timer,"Validation stage of polyhedral end game");
! 473: new_line(resultfile);
! 474: end valipoco;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to valipoco.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Continuation