Annotation of OpenXM_contrib/PHC/Ada/Continuation/valipoco.adb, Revision 1.1.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