Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/integer_pruning_methods.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Floating_Matrices;
3: with Standard_Integer_Norms; use Standard_Integer_Norms;
4: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
5: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
6: with Standard_Integer_Linear_Equalities; use Standard_Integer_Linear_Equalities;
7: with Integer_Linear_Inequalities; use Integer_Linear_Inequalities;
8: with Floating_Linear_Inequality_Solvers; use Floating_Linear_Inequality_Solvers;
9:
10: package body Integer_Pruning_Methods is
11:
12: -- GENERAL AUXILIARIES :
13:
14: function Convert ( fa : Face_Array ) return Array_of_Lists is
15:
16: -- DESCRIPTION :
17: -- Converts the array of faces into an array of lists by
18: -- converting the first element of each list of faces.
19:
20: res : Array_of_Lists(fa'range);
21:
22: begin
23: for k in fa'range loop
24: res(k) := Shallow_Create(fa(k).all);
25: end loop;
26: return res;
27: end Convert;
28:
29: -- AUXILIARIES FOR THE PRUNING ALGORITHMS :
30:
31: procedure Initialize ( n : in natural; mat : out Matrix;
32: ipvt : out Vector ) is
33:
34: -- DESCRIPTION :
35: -- Initializes the matrix of the equality constraints on the normals
36: -- and sets the pivoting vector ipvt to 1..n+1.
37:
38: begin
39: for i in 1..n+1 loop
40: ipvt(i) := i;
41: for j in 1..n+1 loop
42: mat(i,j) := 0;
43: end loop;
44: end loop;
45: end Initialize;
46:
47: function Number_of_Inequalities
48: ( mix : Vector; lifted : Array_of_Lists ) return natural is
49:
50: -- DESCRIPTION :
51: -- Returns the maximal number of inequalities on the inner normal.
52:
53: res : natural := 0;
54:
55: begin
56: for k in lifted'range loop
57: res := res + Length_Of(lifted(k)) - mix(k) - 1;
58: end loop;
59: return res;
60: end Number_of_Inequalities;
61:
62: procedure Ordered_Inequalities ( k : in natural; mat : in out matrix ) is
63:
64: -- DESCRIPTION :
65: -- Defines k inequalities mat(k,k) - mat(k+1,k) >= 0.
66:
67: begin
68: for i in mat'first(1)..k loop
69: for j in mat'range(2) loop
70: mat(i,j) := 0;
71: end loop;
72: mat(i,i) := 1; mat(i,i+1) := -1;
73: end loop;
74: end Ordered_Inequalities;
75:
76: procedure Check_and_Update
77: ( mic : in Face_Array; lifted : in Array_of_Lists;
78: m : in matrix; ipvt : in Vector;
79: mixsub,mixsub_last : in out Mixed_Subdivision ) is
80:
81: -- DESCRIPTION :
82: -- Computes the normal to the points in mic, by solving the
83: -- linear system defined by m and ipvt. Updates the mixed subdivision
84: -- if the computed normal is an inner normal.
85:
86: v : Vector(m'range(2)) := Solve(m,ipvt);
87: pts : Array_of_Lists(mic'range);
88:
89: begin
90: if v(v'last) /= 0
91: then Normalize(v);
92: if v(v'last) < 0
93: then Min(v);
94: end if;
95: pts := Convert(mic);
96: Update(pts,v,mixsub,mixsub_last);
97: end if;
98: end Check_and_Update;
99:
100: procedure Create_Equalities
101: ( n : in natural; f : in Face; mat,ineq : in matrix;
102: ipvt : in Vector; row,rowineq : in natural;
103: newmat,newineq : in out matrix; newipvt : in out Vector;
104: newrow,newrowineq : in out natural; fail : out boolean ) is
105:
106: -- DESCRIPTION :
107: -- Creates new equalities and uses them to eliminate unknowns in the
108: -- matrices of equalities and inequalities. Failure is reported when
109: -- a zero row is encountered. On entry, all new* variables must be
110: -- initialized with the corresponding *-ones.
111:
112: shi : vector(1..n+1) := f(f'first).all;
113: fl : boolean := false;
114: pivot : natural;
115:
116: begin
117: for i in f'first+1..f'last loop
118: newrow := newrow + 1;
119: for j in f(i)'range loop
120: newmat(newrow,j) := f(i)(j) - shi(j);
121: end loop;
122: Triangulate(newrow,newmat,newipvt,pivot);
123: fl := (pivot = 0);
124: exit when fl;
125: Switch(newrow,pivot,ineq'first,rowineq,newineq);
126: --Scale(newrow,newmat);
127: end loop;
128: fail := fl;
129: end Create_Equalities;
130:
131: function Check_Feasibility
132: ( k,m,n : natural; ineq : matrix ) return boolean is
133:
134: -- DESCRIPTION :
135: -- Returns true if -v(n+1) > 0 can be derived, false otherwise.
136:
137: -- ON ENTRY :
138: -- k current unknown that has been eliminated,
139: -- for all i <= k : ineq(l,i) = 0, for l in ineq'first..m;
140: -- m number of inequalities;
141: -- n dimension;
142: -- ineq matrix of inequalities.
143:
144: tableau : matrix(1..n-k+1,1..m+1);
145: feasi : boolean;
146:
147: begin
148: if m = 0
149: then feasi := false;
150: else for i in k+1..n+1 loop
151: for j in 1..m loop
152: tableau(i-k,j) := ineq(j,i);
153: end loop;
154: tableau(i-k,m+1) := 0;
155: end loop;
156: tableau(n-k+1,m+1) := -1;
157: Integer_Complementary_Slackness(tableau,feasi);
158: end if;
159: return feasi;
160: end Check_Feasibility;
161:
162: function New_Check_Feasibility
163: ( k,m,n : natural; ineq : matrix; tol : double_float;
164: solu : Standard_Floating_Vectors.Vector ) return boolean is
165:
166: -- DESCRIPTION :
167: -- Returns true if -v(n+1) > 0 can be derived, false otherwise.
168:
169: -- ON ENTRY :
170: -- k current unknown that has been eliminated,
171: -- for all i <= k : ineq(l,i) = 0, for l in ineq'first..m;
172: -- m number of inequalities;
173: -- n dimension;
174: -- ineq matrix of inequalities.
175:
176: tableau : Standard_Floating_Matrices.matrix(1..n-k+1,1..m);
177: tabsolu : Standard_Floating_Vectors.Vector(1..n-k);
178: cols : Vector(tabsolu'range);
179: kfail : integer;
180: max,tmpabs : double_float;
181: -- used for scaling of the last component of the inner normal
182:
183: begin
184: for i in k+1..n loop
185: for j in 1..m loop
186: tableau(i-k,j) := double_float(ineq(j,i));
187: end loop;
188: tabsolu(i-k) := solu(i);
189: end loop;
190: max := 0.0;
191: for j in 1..m loop
192: tableau(n-k+1,j) := -double_float(ineq(j,n+1));
193: tmpabs := abs(tableau(n-k+1,j));
194: if tmpabs > max
195: then max := tmpabs;
196: end if;
197: end loop;
198: if max > 1.0
199: then for j in 1..m loop
200: tableau(n-k+1,j) := tableau(n-k+1,j)/max;
201: end loop;
202: end if;
203: Scale(tableau,tol);
204: Solve(tableau,tol,tabsolu,kfail,cols);
205: return (kfail <= m);
206: end New_Check_Feasibility;
207:
208: procedure Update_Inequalities
209: ( k,rowmat1,rowmat2,n : in natural;
210: mat : in matrix; ipvt : in vector;
211: rowineq : in out natural; ineq : in out matrix;
212: lifted : in Array_of_Lists; mic : in out Face_Array ) is
213:
214: -- DESCRIPTION :
215: -- The inequalities will be updated w.r.t. the equality
216: -- constraints on the inner normal.
217:
218: tmp : List;
219: pt,shi : Link_to_Vector;
220:
221: begin
222: for i in ineq'first..rowineq loop -- triangulate old inequalities
223: for j in rowmat1..rowmat2 loop
224: Triangulate(j,mat,i,ineq);
225: end loop;
226: --Scale(i,ineq);
227: end loop;
228: tmp := lifted(k); -- create and triangulate
229: shi := mic(k)(mic(k)'first); -- new inequalities
230: while not Is_Null(tmp) loop
231: pt := Head_Of(tmp);
232: if not Is_In(mic(k),pt.all)
233: then rowineq := rowineq + 1;
234: for j in pt'range loop
235: ineq(rowineq,j) := pt(j) - shi(j);
236: end loop;
237: Switch(ipvt,rowineq,ineq);
238: for i in 1..rowmat2 loop
239: Triangulate(i,mat,rowineq,ineq);
240: --Scale(rowineq,ineq);
241: end loop;
242: end if;
243: tmp := Tail_Of(tmp);
244: end loop;
245: end Update_Inequalities;
246:
247: procedure Eliminate ( k : in natural; m : in Matrix; tol : in double_float;
248: x : in out Standard_Floating_Vectors.Vector ) is
249:
250: -- DESCRIPTION :
251: -- Eliminates the kth component of x by using the matrix.
252:
253: fac : double_float;
254:
255: begin
256: if abs(x(k)) > tol
257: then fac := x(k)/double_float(m(k,k));
258: for j in k..x'last loop
259: x(j) := x(j) - fac*double_float(m(k,j));
260: end loop;
261: end if;
262: end Eliminate;
263:
264: procedure Switch ( l,pivot : in integer;
265: x : in out Standard_Floating_Vectors.Vector ) is
266:
267: -- DESCRIPTION :
268: -- Applies the pivotation information to the vector x.
269:
270: tmp : double_float;
271:
272: begin
273: if l /= pivot
274: then tmp := x(l); x(l) := x(pivot); x(pivot) := tmp;
275: end if;
276: end Switch;
277:
278: procedure Eliminate ( rowmat1,rowmat2 : in natural; mat : in Matrix;
279: ipvt : in Vector; tol : in double_float;
280: x : in out Standard_Floating_Vectors.Vector ) is
281:
282: -- DESCRIPTION :
283: -- Returns an eliminated solution vector.
284:
285: begin
286: for k in rowmat1..rowmat2 loop
287: Switch(k,ipvt(k),x);
288: Eliminate(k,mat,tol,x);
289: end loop;
290: end Eliminate;
291:
292: -- GENERAL CONSTRUCTOR in several versions :
293:
294: -- procedure Compute_Mixed_Cells ( ); -- general specification.
295:
296: -- DESCRIPTION :
297: -- Backtrack recursive procedure to enumerate face-face combinations.
298:
299: -- ON ENTRY :
300: -- k index for current point set;
301: -- row number of rows already in the matrix;
302: -- mat matrix which determines the inner normal;
303: -- ipvt contains the pivoting information;
304: -- rowineq number of inequality constraints already in ineq;
305: -- ineq matrix for the inequality constraints on the
306: -- inner normal v, of type <.,v> >= 0;
307: -- testsolu test solution to check current set of inequalilities;
308: -- mic contains the current selected faces, up to k-1.
309:
310: -- ON RETURN :
311: -- mic updated selected faces;
312: -- issupp true if current tuple of faces is supported;
313: -- continue indicates whether to continue the creation or not.
314:
315: -- CONSTRUCTION WITH PRUNING :
316:
317: procedure Gen1_Create_CS
318: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
319: lifted : in Array_of_Lists;
320: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
321: mixsub : out Mixed_Subdivision ) is
322:
323: res,res_last : Mixed_Subdivision;
324: accu : Face_Array(fa'range);
325: ma : matrix(1..n+1,1..n+1);
326: ipvt : vector(1..n+1);
327: ineqrows : natural;
328:
329: procedure Compute_Mixed_Cells
330: ( k,row : in natural; mat : in matrix; ipvt : in vector;
331: rowineq : in natural; ineq : in matrix;
332: mic : in out Face_Array; continue : out boolean );
333:
334: -- DESCRIPTION :
335: -- Backtrack recursive procedure to enumerate face-face combinations.
336:
337: procedure Process_Inequalities
338: ( k,rowmat1,rowmat2 : in natural;
339: mat : in matrix; ipvt : in vector;
340: rowineq : in out natural; ineq : in out matrix;
341: mic : in out Face_Array; cont : out boolean ) is
342:
343: -- DESCRIPTION :
344: -- Updates inequalities and checks feasibility before proceeding.
345:
346: begin
347: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,lifted,mic);
348: if Check_Feasibility(rowmat2,rowineq,n,ineq)
349: then nbfail(k) := nbfail(k) + 1.0;
350: cont := true;
351: else nbsucc(k) := nbsucc(k) + 1.0;
352: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic,cont);
353: end if;
354: end Process_Inequalities;
355:
356: procedure Compute_Mixed_Cells
357: ( k,row : in natural; mat : in matrix; ipvt : in vector;
358: rowineq : in natural; ineq : in matrix;
359: mic : in out Face_Array; continue : out boolean ) is
360:
361: old : Mixed_Subdivision := res_last;
362: cont : boolean := true;
363: tmpfa : Faces;
364:
365: begin
366: if k > mic'last
367: then Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
368: if old /= res_last
369: then Process(Head_Of(res_last),continue);
370: else continue := true;
371: end if;
372: else tmpfa := fa(k);
373: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
374: mic(k) := Head_Of(tmpfa);
375: declare -- update matrices
376: fl : boolean;
377: newipvt : vector(ipvt'range) := ipvt;
378: newmat : matrix(mat'range(1),mat'range(2)) := mat;
379: newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
380: newrow : natural := row;
381: newrowineq : natural := rowineq;
382: begin
383: Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,
384: newmat,newineq,newipvt,newrow,newrowineq,fl);
385: if fl
386: then nbfail(k) := nbfail(k) + 1.0;
387: else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
388: newrowineq,newineq,mic,cont);
389: end if;
390: end;
391: exit when not cont;
392: tmpfa := Tail_Of(tmpfa);
393: end loop;
394: continue := cont;
395: end if;
396: end Compute_Mixed_Cells;
397:
398: begin
399: Initialize(n,ma,ipvt);
400: ineqrows := Number_of_Inequalities(mix,lifted);
401: declare
402: ineq : matrix(1..ineqrows,1..n+1);
403: cont : boolean;
404: begin
405: ineq(1,1) := 0;
406: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu,cont);
407: end;
408: mixsub := res;
409: end Gen1_Create_CS;
410:
411: procedure Create_CS
412: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
413: lifted : in Array_of_Lists;
414: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
415: mixsub : out Mixed_Subdivision ) is
416:
417: res,res_last : Mixed_Subdivision;
418: accu : Face_Array(fa'range);
419: ma : matrix(1..n+1,1..n+1);
420: ipvt : vector(1..n+1);
421: ineqrows : natural;
422:
423: procedure Compute_Mixed_Cells
424: ( k,row : in natural; mat : in Matrix; ipvt : in Vector;
425: rowineq : in natural; ineq : in Matrix;
426: mic : in out Face_Array );
427:
428: -- DESCRIPTION :
429: -- Backtrack recursive procedure to enumerate face-face combinations.
430:
431: procedure Process_Inequalities
432: ( k,rowmat1,rowmat2 : in natural;
433: mat : in matrix; ipvt : in vector;
434: rowineq : in out natural; ineq : in out matrix;
435: mic : in out Face_Array ) is
436:
437: -- DESCRIPTION :
438: -- Updates inequalities and checks feasibility before proceeding.
439:
440: begin
441: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,lifted,mic);
442: if Check_Feasibility(rowmat2,rowineq,n,ineq)
443: then nbfail(k) := nbfail(k) + 1.0;
444: else nbsucc(k) := nbsucc(k) + 1.0;
445: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic);
446: end if;
447: end Process_Inequalities;
448:
449: procedure Compute_Mixed_Cells
450: ( k,row : in natural; mat : in Matrix; ipvt : in Vector;
451: rowineq : in natural; ineq : in Matrix;
452: mic : in out Face_Array ) is
453:
454: tmpfa : Faces;
455:
456: begin
457: if k > mic'last
458: then Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
459: else tmpfa := fa(k);
460: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytype
461: mic(k) := Head_Of(tmpfa);
462: declare -- update matrices
463: fl : boolean;
464: newipvt : vector(ipvt'range) := ipvt;
465: newmat : matrix(mat'range(1),mat'range(2)) := mat;
466: newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
467: newrow : natural := row;
468: newrowineq : natural := rowineq;
469: begin
470: Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,newmat,
471: newineq,newipvt,newrow,newrowineq,fl);
472: if fl
473: then nbfail(k) := nbfail(k) + 1.0;
474: else Process_Inequalities
475: (k,row+1,newrow,newmat,newipvt,newrowineq,newineq,mic);
476: end if;
477: end;
478: tmpfa := Tail_Of(tmpfa);
479: end loop;
480: end if;
481: end Compute_Mixed_Cells;
482:
483: begin
484: Initialize(n,ma,ipvt);
485: ineqrows := Number_of_Inequalities(mix,lifted);
486: declare
487: ineq : matrix(1..ineqrows,1..n+1);
488: begin
489: ineq(1,1) := 0;
490: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu);
491: end;
492: mixsub := res;
493: end Create_CS;
494:
495: procedure New_Create_CS
496: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
497: lifted : in Array_of_Lists;
498: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
499: mixsub : out Mixed_Subdivision ) is
500:
501: res,res_last : Mixed_Subdivision;
502: accu : Face_Array(fa'range);
503: ma : matrix(1..n+1,1..n+1);
504: ipvt : vector(1..n+1);
505: tol : constant double_float := double_float(n)*10.0**(-11);
506: ineqrows : natural;
507:
508: procedure Compute_Mixed_Cells
509: ( k,row : in natural; mat : in matrix; ipvt : in vector;
510: rowineq : in natural; ineq : in matrix;
511: testsolu : in Standard_Floating_Vectors.Vector;
512: mic : in out Face_Array );
513:
514: -- DESCRIPTION :
515: -- Backtrack recursive procedure to enumerate face-face combinations.
516:
517: procedure Process_Inequalities
518: ( k,rowmat1,rowmat2 : in natural;
519: mat : in matrix; ipvt : in vector;
520: rowineq : in out natural; ineq : in out matrix;
521: testsolu : in Standard_Floating_Vectors.Vector;
522: mic : in out Face_Array ) is
523:
524: -- DESCRIPTION :
525: -- Updates the inequalities and checks feasibility before proceeding.
526:
527: tmp : List;
528: pt,shi : Link_to_Vector;
529: newsolu : Standard_Floating_Vectors.Vector(testsolu'range) := testsolu;
530:
531: begin
532: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,
533: lifted,mic);
534: Eliminate(rowmat1,rowmat2,mat,ipvt,tol,newsolu);
535: if New_Check_Feasibility(rowmat2,rowineq,n,ineq,tol,newsolu)
536: then nbfail(k) := nbfail(k) + 1.0;
537: else nbsucc(k) := nbsucc(k) + 1.0;
538: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,testsolu,mic);
539: end if;
540: end Process_Inequalities;
541:
542: procedure Compute_Mixed_Cells
543: ( k,row : in natural; mat : in matrix; ipvt : in vector;
544: rowineq : in natural; ineq : in matrix;
545: testsolu : in Standard_Floating_Vectors.Vector;
546: mic : in out Face_Array ) is
547:
548: -- DESCRIPTION :
549: -- Backtrack recursive procedure to enumerate face-face combinations.
550:
551: tmpfa : Faces;
552:
553: begin
554: if k > mic'last
555: then Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
556: else tmpfa := fa(k);
557: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
558: mic(k) := Head_Of(tmpfa);
559: declare -- update matrices
560: fl : boolean;
561: newipvt : vector(ipvt'range) := ipvt;
562: newmat : matrix(mat'range(1),mat'range(2)) := mat;
563: newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
564: newrow : natural := row;
565: newrowineq : natural := rowineq;
566: begin
567: Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,newmat,
568: newineq,newipvt,newrow,newrowineq,fl);
569: if fl
570: then nbfail(k) := nbfail(k) + 1.0;
571: else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
572: newrowineq,newineq,testsolu,mic);
573: end if;
574: end;
575: tmpfa := Tail_Of(tmpfa);
576: end loop;
577: end if;
578: end Compute_Mixed_Cells;
579:
580: begin
581: Initialize(n,ma,ipvt);
582: ineqrows := Number_of_Inequalities(mix,lifted);
583: declare
584: ineq : matrix(1..ineqrows,1..n+1);
585: solu : Standard_Floating_Vectors.Vector(1..n+1);
586: begin
587: ineq(1,1) := 0;
588: solu := (solu'range => 0.0);
589: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,solu,accu);
590: end;
591: mixsub := res;
592: end New_Create_CS;
593:
594: -- AUXILIARIES FOR THE CONSTRAINT PRUNING :
595:
596: function Is_Supported ( f : Face; normal : Vector; supp : integer )
597: return boolean is
598:
599: ip : integer;
600:
601: begin
602: for i in f'range loop
603: ip := f(i).all*normal;
604: if ip /= supp
605: then return false;
606: end if;
607: end loop;
608: return true;
609: end Is_Supported;
610:
611: function Is_Supported ( f : Face; k : natural; normals,suppvals : List )
612: return boolean is
613:
614: -- DESCRIPTION :
615: -- Returns true if there exists a normal for which the inner product
616: -- with each point in the face equals the corresponding supporting value
617: -- of the kth component.
618:
619: tmpnor : List := normals;
620: tmpsup : List := suppvals;
621: support : integer;
622:
623: begin
624: while not Is_Null(tmpnor) loop
625: support := Head_Of(tmpsup)(k);
626: if Is_Supported(f,Head_Of(tmpnor).all,support)
627: then return true;
628: else tmpnor := Tail_Of(tmpnor);
629: tmpsup := Tail_Of(tmpsup);
630: end if;
631: end loop;
632: return false;
633: end Is_Supported;
634:
635: procedure Update ( mic : in Mixed_Cell; normals,suppvals : in out List ) is
636:
637: -- DESCRIPTION :
638: -- Updates the list of normals and supporting values with the information
639: -- of a mixed cell.
640:
641: normal : Link_to_Vector := new Vector'(mic.nor.all);
642: suppval : Link_to_Vector := new Vector(mic.pts'range);
643:
644: begin
645: Construct(normal,normals);
646: for i in suppval'range loop
647: suppval(i) := normal*Head_Of(mic.pts(i));
648: end loop;
649: Construct(suppval,suppvals);
650: end Update;
651:
652: procedure Create_CCS
653: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
654: lifted : in Array_of_Lists;
655: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
656: normals,suppvals : in out List;
657: mixsub : out Mixed_Subdivision ) is
658:
659: res,res_last : Mixed_Subdivision;
660: accu : Face_Array(fa'range);
661: ma : matrix(1..n+1,1..n+1);
662: ipvt : vector(1..n+1);
663: ineqrows : natural;
664:
665: procedure Compute_Mixed_Cells
666: ( k,row : in natural; mat : in matrix; ipvt : in vector;
667: rowineq : in natural; ineq : in matrix;
668: mic : in out Face_Array; issupp : in out boolean );
669:
670: -- DESCRIPTION :
671: -- Backtrack recursive procedure to enumerate face-face combinations.
672:
673: procedure Process_Inequalities
674: ( k,rowmat1,rowmat2 : in natural;
675: mat : in matrix; ipvt : in vector;
676: rowineq : in out natural; ineq : in out matrix;
677: mic : in out Face_Array; issupp : in out boolean ) is
678:
679: -- DESCRIPTION :
680: -- Updates inequalities and checks feasibility before proceeding.
681:
682: begin
683: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,rowineq,ineq,lifted,mic);
684: if issupp
685: then issupp := Is_Supported(mic(k),k,normals,suppvals);
686: end if;
687: if not issupp and then Check_Feasibility(rowmat2,rowineq,n,ineq)
688: then nbfail(k) := nbfail(k) + 1.0;
689: else nbsucc(k) := nbsucc(k) + 1.0;
690: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic,issupp);
691: end if;
692: end Process_Inequalities;
693:
694: procedure Compute_Mixed_Cells
695: ( k,row : in natural; mat : in matrix; ipvt : in vector;
696: rowineq : in natural; ineq : in matrix;
697: mic : in out Face_Array; issupp : in out boolean ) is
698:
699: -- DESCRIPTION :
700: -- Backtrack recursive procedure to enumerate face-face combinations.
701:
702: old : Mixed_Subdivision;
703: tmpfa : Faces;
704:
705: begin
706: if k > mic'last
707: then old := res_last;
708: Check_and_Update(mic,lifted,mat,ipvt,res,res_last);
709: if old /= res_last
710: then Update(Head_Of(res_last),normals,suppvals);
711: end if;
712: else tmpfa := fa(k);
713: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
714: mic(k) := Head_Of(tmpfa);
715: declare -- update matrices
716: fl : boolean;
717: newipvt : vector(ipvt'range) := ipvt;
718: newmat : matrix(mat'range(1),mat'range(2)) := mat;
719: newineq : matrix(ineq'range(1),ineq'range(2)) := ineq;
720: newrow : natural := row;
721: newrowineq : natural := rowineq;
722: begin
723: Create_Equalities(n,mic(k),mat,ineq,ipvt,row,rowineq,newmat,
724: newineq,newipvt,newrow,newrowineq,fl);
725: if fl
726: then nbfail(k) := nbfail(k) + 1.0;
727: else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
728: newrowineq,newineq,mic,issupp);
729: end if;
730: end;
731: tmpfa := Tail_Of(tmpfa);
732: end loop;
733: end if;
734: end Compute_Mixed_Cells;
735:
736: begin
737: Initialize(n,ma,ipvt);
738: ineqrows := Number_of_Inequalities(mix,lifted);
739: declare
740: ineq : matrix(1..ineqrows,1..n+1);
741: supported : boolean := true;
742: begin
743: ineq(1,1) := 0;
744: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu,supported);
745: end;
746: mixsub := res;
747: end Create_CCS;
748:
749: end Integer_Pruning_Methods;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to integer_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