Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/floating_pruning_methods.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Natural_Vectors;
2: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
3: with Standard_Floating_Linear_Solvers; use Standard_Floating_Linear_Solvers;
4: with Floating_Linear_Inequalities; use Floating_Linear_Inequalities;
5: with Lists_of_Floating_Vectors; use Lists_of_Floating_Vectors;
6:
7: package body Floating_Pruning_Methods is
8:
9: -- GENERAL AUXILIARIES :
10:
11: procedure Normalize ( v : in out Standard_Floating_Vectors.Vector ) is
12:
13: -- DESCRIPTION : Divides every entry by v(v'last).
14:
15: begin
16: for i in v'range loop
17: v(i) := v(i)/v(v'last);
18: end loop;
19: end Normalize;
20:
21: function Convert ( fa : Face_Array ) return Array_of_Lists is
22:
23: -- DESCRIPTION :
24: -- Converts the array of faces into an array of lists by
25: -- converting the first element of each list of faces.
26:
27: res : Array_of_Lists(fa'range);
28:
29: begin
30: for k in fa'range loop
31: res(k) := Shallow_Create(fa(k).all);
32: end loop;
33: return res;
34: end Convert;
35:
36: -- AUXILIARIES FOR THE PRUNING ALGORITHMS :
37:
38: procedure Initialize ( n : in natural; mat : out Matrix;
39: ipvt : out Standard_Natural_Vectors.Vector ) is
40:
41: -- DESCRIPTION :
42: -- Initializes an n*(n+1) matrix with zeroes and the pivoting vector
43: -- with the entries 1..n.
44:
45: begin
46: for i in 1..n loop
47: for j in 1..n+1 loop
48: mat(i,j) := 0.0;
49: end loop;
50: end loop;
51: for i in 1..n loop
52: ipvt(i) := i;
53: end loop;
54: end Initialize;
55:
56: function Number_of_Inequalities
57: ( mix : Vector; lifted : Array_of_Lists ) return natural is
58:
59: -- DESCRIPTION :
60: -- Returns the maximal number of inequalities for pruning.
61:
62: res : natural := 0;
63:
64: begin
65: for k in lifted'range loop
66: res := res + Length_Of(lifted(k)) - mix(k) - 1;
67: end loop;
68: return res;
69: end Number_of_Inequalities;
70:
71: procedure Ordered_Inequalities ( k : in natural; mat : in out Matrix ) is
72:
73: -- DESCRIPTION :
74: -- Defines k inequalities mat(k,k) - mat(k+1,k) >= 0.
75:
76: begin
77: for i in mat'first(1)..k loop
78: for j in mat'range(2) loop
79: mat(i,j) := 0.0;
80: end loop;
81: mat(i,i) := 1.0; mat(i,i+1) := -1.0;
82: end loop;
83: end Ordered_Inequalities;
84:
85: procedure Check_and_Update
86: ( mic : in Face_Array; lifted : in Array_of_Lists;
87: m : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
88: tol : in double_float;
89: mixsub,mixsub_last : in out Mixed_Subdivision ) is
90:
91: -- DESCRIPTION :
92: -- Computes the normal to the points in pts, by solving the
93: -- linear system defined by m and ipvt.
94: -- If the computed normal is an inner normal w.r.t. the lifted points,
95: -- then the mixed subdivision will be updated with a new cell.
96:
97: use Standard_Floating_Vectors;
98: v : Standard_Floating_Vectors.Vector(m'range(2)) := Solve(m,tol,ipvt);
99: pts : Array_of_Lists(mic'range);
100:
101: begin
102: if abs(v(v'last)) > tol
103: then Normalize(v);
104: if v(v'last) < 0.0
105: then Min(v);
106: end if;
107: -- pts := Convert(mic);
108: -- Update(pts,v,mixsub,mixsub_last);
109: declare
110: cell : Mixed_Cell;
111: begin
112: cell.nor := new Standard_Floating_Vectors.Vector'(v);
113: cell.pts := new Array_of_Lists'(Convert(mic));
114: cell.sub := null;
115: Append(mixsub,mixsub_last,cell);
116: end;
117: end if;
118: end Check_and_Update;
119:
120: procedure Create_Equalities
121: ( n : in natural; f : in Face; mat,ineq : in Matrix;
122: tol : in double_float;
123: ipvt : in Standard_Natural_Vectors.Vector;
124: row,rowineq : in natural;
125: newmat,newineq : in out Matrix;
126: newipvt : in out Standard_Natural_Vectors.Vector;
127: newrow,newrowineq : in out natural; fail : out boolean ) is
128:
129: -- DESCRIPTION :
130: -- Creates new equalities and uses them to eliminate unknowns in the
131: -- matrices of equalities and inequalities. Failure is reported when
132: -- a zero row is encountered. On entry, all new* variables must be
133: -- initialized with the corresponding *-ones.
134:
135: shi : Standard_Floating_Vectors.Vector(1..n+1) := f(f'first).all;
136: fl : boolean := false;
137: pivot : natural;
138:
139: begin
140: for i in f'first+1..f'last loop
141: newrow := newrow + 1;
142: for j in f(i)'range loop
143: newmat(newrow,j) := f(i)(j) - shi(j);
144: end loop;
145: Switch(newipvt,newrow,newmat);
146: Upper_Triangulate(newrow,newmat,tol,newipvt,pivot);
147: fl := (pivot = 0);
148: exit when fl;
149: Switch(newrow,pivot,ineq'first,rowineq,newineq);
150: end loop;
151: fail := fl;
152: end Create_Equalities;
153:
154: procedure Complementary_Slackness
155: ( tableau : in out matrix;
156: tol : in double_float; feasible : out boolean ) is
157:
158: lastcol : constant integer := tableau'last(2)-1;
159: rhs,sol : Standard_Floating_Vectors.Vector(tableau'range(1));
160: columns : Vector(sol'range);
161:
162: begin
163: for i in rhs'range loop
164: rhs(i) := double_float(tableau(i,tableau'last(2)));
165: end loop;
166: Complementary_Slackness(tableau,lastcol,rhs,tol,sol,columns,feasible);
167: end Complementary_Slackness;
168:
169: function Check_Feasibility ( k,m,n : natural; ineq : Matrix;
170: tol : double_float ) return boolean is
171:
172: -- DESCRIPTION :
173: -- Returns true if -v(n+1) > 0 can be derived, false otherwise.
174:
175: -- ON ENTRY :
176: -- k current unknown that has been eliminated,
177: -- for all i <= k : ineq(l,i) = 0, for l in ineq'first..m;
178: -- m number of inequalities;
179: -- n dimension;
180: -- ineq matrix of inequalities.
181:
182: tableau : Matrix(1..n-k+1,1..m+1);
183: feasi : boolean;
184:
185: begin
186: if m = 0
187: then feasi := false;
188: else for i in k+1..n+1 loop
189: for j in 1..m loop
190: tableau(i-k,j) := ineq(j,i);
191: end loop;
192: tableau(i-k,m+1) := 0.0;
193: end loop;
194: tableau(n-k+1,m+1) := -1.0;
195: Complementary_Slackness(tableau,tol,feasi);
196: end if;
197: return feasi;
198: end Check_Feasibility;
199:
200: procedure Update_Inequalities
201: ( k,rowmat1,rowmat2,n : in natural;
202: mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
203: tol : in double_float;
204: rowineq : in out natural; ineq : in out Matrix;
205: lifted : in Array_of_Lists; mic : in out Face_Array ) is
206:
207: -- DESCRIPTION :
208: -- The inequalities will be updated w.r.t. the equality
209: -- constraints on the inner normal.
210:
211: tmp : List;
212: pt,shi : Standard_Floating_Vectors.Link_to_Vector;
213:
214: begin
215: for i in ineq'first..rowineq loop -- update the old inequalities
216: for j in rowmat1..rowmat2 loop
217: Upper_Triangulate(j,mat,tol,i,ineq);
218: end loop;
219: end loop;
220: shi := mic(k)(mic(k)'first);
221: tmp := lifted(k); -- make new inequalities
222: while not Is_Null(tmp) loop
223: pt := Head_Of(tmp);
224: if not Is_In(mic(k),pt.all)
225: then rowineq := rowineq + 1;
226: for j in pt'range loop
227: ineq(rowineq,j) := pt(j) - shi(j);
228: end loop;
229: Switch(ipvt,rowineq,ineq);
230: for i in 1..rowmat2 loop
231: Upper_Triangulate(i,mat,tol,rowineq,ineq);
232: end loop;
233: end if;
234: tmp := Tail_Of(tmp);
235: end loop;
236: end Update_Inequalities;
237:
238: -- CONSTRUCTION WITH PRUNING :
239:
240: procedure Gen1_Create
241: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
242: lifted : in Array_of_Lists; tol : in double_float;
243: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
244: mixsub : out Mixed_Subdivision ) is
245:
246: res,res_last : Mixed_Subdivision;
247: accu : Face_Array(fa'range);
248: ma : Matrix(1..n,1..n+1);
249: ipvt : Standard_Natural_Vectors.Vector(1..n);
250: ineqrows : natural;
251:
252: procedure Compute_Mixed_Cells
253: ( k,row : in natural; mat : in Matrix;
254: ipvt : in Standard_Natural_Vectors.Vector;
255: rowineq : in natural; ineq : in Matrix;
256: mic : in out Face_Array; continue : out boolean );
257:
258: -- DESCRIPTION :
259: -- Backtrack recursive procedure to enumerate face-face combinations.
260:
261: -- ON ENTRY :
262: -- k index for current point set;
263: -- row number of rows already in the matrix;
264: -- mat matrix which determines the inner normal;
265: -- ipvt contains the pivoting information;
266: -- rowineq number of inequality constraints already in ineq;
267: -- ineq matrix for the inequality constraints on the
268: -- inner normal v, of type <.,v> >= 0;
269: -- mic contains the current selected faces, up to k-1.
270:
271: -- ON RETURN :
272: -- mic updated selected faces.
273: -- continue indicates whether to continue the creation or not.
274:
275: procedure Process_Inequalities
276: ( k,rowmat1,rowmat2 : in natural;
277: mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
278: rowineq : in out natural; ineq : in out Matrix;
279: mic : in out Face_Array; cont : out boolean ) is
280:
281: -- DESCRIPTION :
282: -- Updates inequalities and checks feasibility before proceeding.
283:
284: fl : boolean := false;
285: tmp : List;
286: pt,shi : Link_to_Vector;
287:
288: begin
289: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,tol,rowineq,ineq,
290: lifted,mic);
291: if Check_Feasibility(rowmat2,rowineq,n,ineq,tol)
292: then nbfail(k) := nbfail(k) + 1.0;
293: cont := true;
294: else nbsucc(k) := nbsucc(k) + 1.0;
295: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic,cont);
296: end if;
297: end Process_Inequalities;
298:
299: procedure Compute_Mixed_Cells
300: ( k,row : in natural; mat : in Matrix;
301: ipvt : in Standard_Natural_Vectors.Vector;
302: rowineq : in natural; ineq : in Matrix;
303: mic : in out Face_Array; continue : out boolean ) is
304:
305: old : Mixed_Subdivision := res_last;
306: cont : boolean := true;
307: tmpfa : Faces;
308:
309: begin
310: if k > mic'last
311: then Check_and_Update(mic,lifted,mat,ipvt,tol,res,res_last);
312: if old /= res_last
313: then Process(Head_Of(res_last),continue);
314: else continue := true;
315: end if;
316: else tmpfa := fa(k);
317: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
318: mic(k) := Head_Of(tmpfa);
319: declare -- update the matrices
320: fl : boolean;
321: newipvt : Standard_Natural_Vectors.Vector(ipvt'range) := ipvt;
322: newmat : Matrix(mat'range(1),mat'range(2)) := mat;
323: newineq : Matrix(ineq'range(1),ineq'range(2)) := ineq;
324: newrow : natural := row;
325: newrowineq : natural := rowineq;
326: begin
327: Create_Equalities
328: (n,mic(k),mat,ineq,tol,ipvt,row,rowineq,newmat,newineq,
329: newipvt,newrow,newrowineq,fl);
330: if fl
331: then nbfail(k) := nbfail(k) + 1.0;
332: else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
333: newrowineq,newineq,mic,cont);
334: end if;
335: end;
336: tmpfa := Tail_Of(tmpfa);
337: exit when not cont;
338: end loop;
339: continue := cont;
340: end if;
341: end Compute_Mixed_Cells;
342:
343: begin
344: Initialize(n,ma,ipvt);
345: ineqrows := Number_of_Inequalities(mix,lifted);
346: declare
347: ineq : matrix(1..ineqrows,1..n+1);
348: cont : boolean;
349: begin
350: ineq(1,1) := 0.0;
351: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu,cont);
352: end;
353: mixsub := res;
354: end Gen1_Create;
355:
356: procedure Create
357: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
358: lifted : in Array_of_Lists; tol : in double_float;
359: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
360: mixsub : out Mixed_Subdivision ) is
361:
362: res,res_last : Mixed_Subdivision;
363: accu : Face_Array(fa'range);
364: ma : Matrix(1..n,1..n+1);
365: ipvt : Standard_Natural_Vectors.Vector(1..n);
366: ineqrows : natural;
367:
368: procedure Compute_Mixed_Cells
369: ( k,row : in natural; mat : in Matrix;
370: ipvt : in Standard_Natural_Vectors.Vector;
371: rowineq : in natural; ineq : in Matrix;
372: mic : in out Face_Array );
373:
374: -- DESCRIPTION :
375: -- Backtrack recursive procedure to enumerate face-face combinations.
376:
377: -- ON ENTRY :
378: -- k index for current point set;
379: -- row number of rows already in the matrix;
380: -- mat matrix which determines the inner normal;
381: -- ipvt contains the pivoting information;
382: -- rowineq number of inequality constraints already in ineq;
383: -- ineq matrix for the inequality constraints on the
384: -- inner normal v, of type <.,v> >= 0;
385: -- mic contains the current selected faces, up to k-1.
386:
387: -- ON RETURN :
388: -- mic updated selected faces.
389:
390: procedure Process_Inequalities
391: ( k,rowmat1,rowmat2 : in natural;
392: mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
393: rowineq : in out natural; ineq : in out matrix;
394: mic : in out Face_Array) is
395:
396: -- DESCRIPTION :
397: -- Updates inequalities and checks feasibility before proceeding.
398:
399: tmp : List;
400: pt,shi : Link_to_Vector;
401:
402: begin
403: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,tol,rowineq,ineq,
404: lifted,mic);
405: if Check_Feasibility(rowmat2,rowineq,n,ineq,tol)
406: then nbfail(k) := nbfail(k) + 1.0;
407: else nbsucc(k) := nbsucc(k) + 1.0;
408: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic);
409: end if;
410: end Process_Inequalities;
411:
412: procedure Compute_Mixed_Cells
413: ( k,row : in natural; mat : in Matrix;
414: ipvt : in Standard_Natural_Vectors.Vector;
415: rowineq : in natural; ineq : in Matrix;
416: mic : in out Face_Array ) is
417:
418: tmpfa : Faces;
419:
420: begin
421: if k > mic'last
422: then Check_and_Update(mic,lifted,mat,ipvt,tol,res,res_last);
423: else tmpfa := fa(k);
424: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
425: mic(k) := Head_Of(tmpfa);
426: declare -- update matrices
427: fl : boolean;
428: newipvt : Standard_Natural_Vectors.Vector(ipvt'range) := ipvt;
429: newmat : Matrix(mat'range(1),mat'range(2)) := mat;
430: newineq : Matrix(ineq'range(1),ineq'range(2)) := ineq;
431: newrow : natural := row;
432: newrowineq : natural := rowineq;
433: begin
434: Create_Equalities
435: (n,mic(k),mat,ineq,tol,ipvt,row,rowineq,newmat,newineq,
436: newipvt,newrow,newrowineq,fl);
437: if fl
438: then nbfail(k) := nbfail(k) + 1.0;
439: else Process_Inequalities
440: (k,row+1,newrow,newmat,newipvt,newrowineq,newineq,mic);
441: end if;
442: end;
443: tmpfa := Tail_Of(tmpfa);
444: end loop;
445: end if;
446: end Compute_Mixed_Cells;
447:
448: begin
449: Initialize(n,ma,ipvt);
450: ineqrows := Number_of_Inequalities(mix,lifted);
451: declare
452: ineq : Matrix(1..ineqrows,1..n+1);
453: begin
454: ineq(1,1) := 0.0;
455: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu);
456: end;
457: mixsub := res;
458: end Create;
459:
460: end Floating_Pruning_Methods;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to floating_pruning_methods.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Stalift