Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/integer_face_enumerators.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
2: with Standard_Common_Divisors; use Standard_Common_Divisors;
3: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
4: with Face_Enumerators_Utilities; use Face_Enumerators_Utilities;
5: with Integer_Linear_Inequalities; use Integer_Linear_Inequalities;
6:
7: package body Integer_Face_Enumerators is
8:
9: -- ELIMINATE TO MAKE UPPER TRIANGULAR :
10:
11: procedure Eliminate ( pivot : in integer; elim : in Vector;
12: target : in out Vector ) is
13:
14: -- DESCRIPTION :
15: -- Makes target(pivot) = 0 by means of making a linear
16: -- combination of the vectors target and elim.
17:
18: -- REQUIRED ON ENTRY :
19: -- target(pivot) /= 0 and elim(pivot) /= 0
20:
21: a,b,lcmab,faca,facb : integer;
22:
23: begin
24: a := elim(pivot); b := target(pivot);
25: lcmab := lcm(a,b);
26: if lcmab < 0 then lcmab := -lcmab; end if;
27: faca := lcmab/a; facb := lcmab/b;
28: if facb < 0
29: then facb := -facb; faca := -faca;
30: end if;
31: for j in target'range loop
32: target(j) := facb*target(j) - faca*elim(j);
33: end loop;
34: Scale(target);
35: end Eliminate;
36:
37: procedure Eliminate ( l : in natural; pivots : in Vector;
38: elim : in VecVec; target : in out Vector ) is
39:
40: -- DESCRIPTION :
41: -- Makes target(pivots(i)) = 0 by means of making a linear
42: -- combination of the vectors target and elim(i), for i in 1..l.
43:
44: -- REQUIRED ON ENTRY :
45: -- elim(i)(pivots(i)) /= 0
46:
47: begin
48: for i in 1..l loop
49: if target(pivots(i)) /= 0
50: then Eliminate(pivots(i),elim(i).all,target);
51: end if;
52: end loop;
53: end Eliminate;
54:
55: function Pivot_after_Elimination
56: ( l,k : in natural; point,face,pivots : in Vector;
57: elim : in VecVec ) return natural is
58:
59: -- DESCRIPTION :
60: -- Returns the first nonzero element of the given point after elimination
61: -- w.r.t. the entries in the face with lower index.
62:
63: work : Vector(point'range) := point - elim(1-k).all;
64: pivot : integer;
65:
66: begin
67: for i in (face'first+1)..face'last loop
68: if (face(i) < l) and then (work(pivots(i)) /= 0)
69: then Eliminate(pivots(i),elim(i).all,work);
70: end if;
71: exit when (face(i) > l);
72: end loop;
73: pivot := 0;
74: for i in work'range loop
75: if work(pivots(i)) /= 0
76: then pivot := i;
77: end if;
78: exit when (pivot /= 0);
79: end loop;
80: return pivot;
81: end Pivot_after_Elimination;
82:
83: procedure Update_Pivots ( point : in Vector; l : in natural;
84: pivots : in out Vector; pivot : out natural ) is
85:
86: -- DESCRIPTION :
87: -- Searches in the point(l..point'last) for the first nonzero entry
88: -- and updates the pivoting information.
89:
90: temp,piv : integer;
91:
92: begin
93: piv := 0;
94: for i in l..point'last loop
95: if point(pivots(i)) /= 0
96: then piv := i;
97: end if;
98: exit when (piv /= 0);
99: end loop;
100: if piv /= 0 and then (piv /= l)
101: then temp := pivots(l);
102: pivots(l) := pivots(piv);
103: pivots(piv) := temp;
104: end if;
105: pivot := piv;
106: end Update_Pivots;
107:
108: procedure Update_Eliminator ( elim : in out VecVec; l : in natural;
109: pivots : in out Vector;
110: point : in Vector; pivot : out natural ) is
111:
112: -- DESCRIPTION :
113: -- Updates the vector of eliminators, by adding the lth elimination
114: -- equation. This procedure ensures the invariant condition on the
115: -- eliminator, which has to be upper triangular. If this cannot be
116: -- achieved, degeneracy is indicated by pivot = 0.
117:
118: -- ON ENTRY :
119: -- elim vector of elimination equations: elim(i)(pivots(i)) /= 0
120: -- and for j in 1..(i-1) : elim(i)(pivots(j)) = 0,
121: -- for i in 1..(l-1), elim(0) contains the basis point;
122: -- l index of current elimination equation to be made;
123: -- pivots contains the pivoting information;
124: -- point new point to make the equation `point - elim(0) = 0'.
125:
126: -- ON RETURN :
127: -- elim if not degen, then elim(l)(pivots(l)) /= 0 and
128: -- for i in 1..(l-1): elim(l)(pivots(i)) = 0;
129: -- pivots updated pivot information;
130: -- pivot the pivot that has been used for elim(l)(pivots(l)) /= 0;
131: -- piv = 0 when the new elimination equation elim(l)
132: -- became entirely zero after ensuring the invariant condition.
133:
134: begin
135: elim(l) := new Vector'(point - elim(0).all);
136: Eliminate(l-1,pivots,elim,elim(l).all);
137: Update_Pivots(elim(l).all,l,pivots,pivot);
138: end Update_Eliminator;
139:
140: procedure Update_Eliminator_for_Sum
141: ( elim : in out VecVec; l : in natural; pivots : in out Vector;
142: point : in Vector; index : in natural; pivot : out natural ) is
143:
144: -- DESCRIPTION :
145: -- Updates the vector of eliminators, by adding the lth elimination
146: -- equation. This procedure ensures the invariant condition on the
147: -- eliminator, which has to be upper triangular. If this cannot be
148: -- achieved, degeneracy is indicated by pivot = 0.
149:
150: -- ON ENTRY :
151: -- elim vector of elimination equations: elim(i)(pivots(i)) /= 0
152: -- and for j in 1..(i-1) : elim(i)(pivots(j)) = 0,
153: -- for i in 1..(l-1), elim(1-index) contains the basis point;
154: -- l index of current elimination equation to be made;
155: -- pivots contains the pivoting information;
156: -- point new point to make the equation `point - elim(1-index) = 0';
157: -- index indicates the current polytope.
158:
159: -- ON RETURN :
160: -- elim if not degen, then elim(l)(pivots(l)) /= 0 and
161: -- for i in 1..(l-1): elim(l)(pivots(i)) = 0;
162: -- pivots updated pivot information;
163: -- piv the pivot that has been used for elim(l)(pivots(l)) /= 0;
164: -- piv = 0 when the new elimination equation elim(l)
165: -- became entirely zero after ensuring the invariant condition.
166:
167: begin
168: elim(l) := new Vector'(point - elim(1-index).all);
169: Eliminate(l-1,pivots,elim,elim(l).all);
170: Update_Pivots(elim(l).all,l,pivots,pivot);
171: end Update_Eliminator_for_Sum;
172:
173: -- CREATE TABLEAU OF INEQUALITIES :
174:
175: procedure Create_Tableau_for_Vertices
176: ( i,n : in integer; pts : in VecVec; tab : out matrix ) is
177:
178: column : integer := tab'first(2);
179:
180: begin
181: for j in pts'range loop
182: if j /= i
183: then
184: for k in pts(j)'range loop
185: tab(k,column) := pts(j)(k);
186: end loop;
187: tab(tab'last(1),column) := 1;
188: column := column + 1;
189: end if;
190: end loop;
191: for k in pts(i)'range loop
192: tab(k,tab'last(2)) := pts(i)(k);
193: end loop;
194: tab(tab'last(1),tab'last(2)) := 1;
195: end Create_Tableau_for_Vertices;
196:
197: procedure Create_Tableau_for_Edges
198: ( i,j,n,pivot : in integer; elim : in Vector;
199: pts : in VecVec; tab : out matrix;
200: lastcol : out integer; degenerate : out boolean ) is
201:
202: column : integer := 1;
203: ineq : Vector(1..n);
204:
205: begin
206: degenerate := false;
207: for k in pts'range loop
208: if (k/=i) and then (k/=j)
209: then
210: ineq := pts(k).all - pts(i).all; -- compute inequality
211: -- put("Inequality before eliminate : "); put(ineq); new_line;
212: if ineq(pivot) /= 0 -- make ineq(pivot) = 0
213: then Eliminate(pivot,elim,ineq);
214: end if;
215: -- put("Inequality after eliminate : "); put(ineq); new_line;
216: if Is_Zero(ineq) -- check if degenerate
217: then
218: if not In_Edge(pts(k).all,pts(i).all,pts(j).all)
219: then degenerate := true; return;
220: end if;
221: else
222: for l in 1..n loop -- fill in the column
223: if l < pivot
224: then tab(l,column) := ineq(l);
225: elsif l > pivot
226: then tab(l-1,column) := ineq(l);
227: end if;
228: end loop;
229: tab(tab'last(1),column) := 1;
230: column := column + 1;
231: end if;
232: end if;
233: end loop;
234: for k in tab'first(1)..(tab'last(1)-1) loop -- right hand side
235: tab(k,tab'last(2)) := 0;
236: end loop;
237: tab(tab'last(1),tab'last(2)) := 1;
238: lastcol := column-1;
239: end Create_Tableau_for_Edges;
240:
241: procedure Create_Tableau_for_Faces
242: ( k,n : in natural; face,pivots : in Vector;
243: pts,elim : in VecVec; tab : out matrix; lastcol : out integer;
244: degenerate : out boolean ) is
245:
246: column : integer := 1;
247: ineq : Vector(1..n);
248:
249: begin
250: degenerate := false;
251: for l in pts'range loop
252: if not Is_In(l,face)
253: then
254: ineq := pts(l).all - elim(0).all; -- new inequality
255: -- put("Inequality before : "); put(ineq); new_line;
256: Eliminate(k,pivots,elim,ineq); -- ineq(pivots(i)) = 0, i=1,2,..k
257: -- put("and after elimination : "); put(ineq); new_line;
258: if Is_Zero(ineq)
259: then
260: if --not In_Face(k,face,pts(l).all,pts)
261: --and then
262: l < face(face'last) -- lexicographic enumeration
263: and then (Pivot_after_Elimination
264: (l,1,pts(l).all,face,pivots,elim) /= 0)
265: then -- put("Degenerate for l = "); put(l,1); new_line;
266: degenerate := true; return;
267: end if;
268: else
269: for ll in (k+1)..n loop -- fill in the column
270: tab(ll-k,column) := ineq(pivots(ll));
271: end loop;
272: tab(tab'last(1),column) := 1;
273: column := column + 1;
274: end if;
275: end if;
276: end loop;
277: for l in tab'first(1)..(tab'last(1)-1) loop -- right hand side
278: tab(l,tab'last(2)) := 0;
279: end loop;
280: tab(tab'last(1),tab'last(2)) := 1;
281: lastcol := column-1;
282: end Create_Tableau_for_Faces;
283:
284: procedure Create_Tableau_for_Faces_of_Sum
285: ( k,n,i,r : in natural; ind,pivots : in Vector;
286: pts,elim,face : in VecVec; tab : out matrix;
287: lastcol : out integer; degenerate : out boolean ) is
288:
289: -- DESCRIPTION :
290: -- Creates the table of inequalities for determining whether the given
291: -- candidate face spans a face of the sum polytope.
292:
293: -- ON ENTRY :
294: -- k dimension of the face on the sum;
295: -- n dimension of the points;
296: -- i current polytope;
297: -- r number of polytopes;
298: -- ind indicate the beginning vector in pts of each polytope;
299: -- etc...
300:
301: column : integer := 1;
302: ineq : Vector(1..n);
303: last : integer;
304:
305: begin
306: degenerate := false;
307: for l1 in face'first..i loop
308: if l1 = r
309: then last := pts'last;
310: else last := ind(l1+1)-1;
311: end if;
312: for l2 in ind(l1)..last loop
313: if not Is_In(l2,face(l1).all)
314: then
315: ineq := pts(l2).all - elim(1-l1).all; -- new inequality
316: Eliminate(k,pivots,elim,ineq); -- ineq(pivots(i)) = 0, i=1,2,..k
317: if Is_Zero(ineq)
318: then
319: if --not In_Face(face(l1)'length-1,face(l1).all,pts(l2).all,pts)
320: --and then
321: face(l1)(face(l1)'first) <= l2
322: and then l2 < face(l1)(face(l1)'last)
323: -- lexicographic enumeration
324: and then (Pivot_after_Elimination
325: (l2,l1,pts(l2).all,face(l1).all,pivots,elim) /= 0)
326: then degenerate := true; return;
327: end if;
328: else
329: for ll in (k+1)..n loop -- fill in the column
330: tab(ll-k,column) := ineq(pivots(ll));
331: end loop;
332: tab(tab'last(1),column) := 1;
333: column := column + 1;
334: end if;
335: end if;
336: end loop;
337: end loop;
338: for l in tab'first(1)..(tab'last(1)-1) loop -- right hand side
339: tab(l,tab'last(2)) := 0;
340: end loop;
341: tab(tab'last(1),tab'last(2)) := 1;
342: lastcol := column-1;
343: end Create_Tableau_for_Faces_of_Sum;
344:
345: procedure Create_Tableau_for_Lower_Edges
346: ( i,j,n,pivot : in integer; elim : in Vector;
347: pts : in VecVec; tab : out matrix;
348: lastcol : out integer; degenerate : out boolean ) is
349:
350: column : integer := 1;
351: ineq : Vector(1..n);
352:
353: begin
354: -- put("The elimination vector : "); put(elim); new_line;
355: degenerate := false;
356: for k in pts'range loop
357: if (k/=i) and then (k/=j)
358: then
359: ineq := pts(k).all - pts(i).all; -- compute inequality
360: -- put("Inequality before eliminate : "); put(ineq); new_line;
361: if ineq(pivot) /= 0 -- make ineq(pivot) = 0
362: then Eliminate(pivot,elim,ineq);
363: end if;
364: -- put("Inequality after eliminate : "); put(ineq); new_line;
365: if Is_Zero(ineq) -- check if degenerate
366: then
367: if not In_Edge(pts(k).all,pts(i).all,pts(j).all)
368: then degenerate := true; return;
369: end if;
370: else
371: for l in 1..n loop -- fill in the column
372: if l < pivot
373: then tab(l,column) := ineq(l);
374: elsif l > pivot
375: then tab(l-1,column) := ineq(l);
376: end if;
377: end loop;
378: column := column + 1;
379: end if;
380: end if;
381: end loop;
382: for k in tab'first(1)..(tab'last(1)-1) loop -- right hand side
383: tab(k,tab'last(2)) := 0;
384: end loop;
385: tab(tab'last(1),tab'last(2)) := -1;
386: lastcol := column-1;
387: end Create_Tableau_for_Lower_Edges;
388:
389: procedure Create_Tableau_for_Lower_Faces
390: ( k,n : in natural; face,pivots : in Vector;
391: pts,elim : in VecVec; tab : out matrix; lastcol : out integer;
392: degenerate : out boolean ) is
393:
394: column : integer := 1;
395: ineq : Vector(1..n);
396:
397: begin
398: degenerate := false;
399: for l in pts'range loop
400: if not Is_In(l,face)
401: then
402: ineq := pts(l).all - elim(0).all; -- new inequality
403: -- put("Inequality before : "); put(ineq); new_line;
404: Eliminate(k,pivots,elim,ineq); -- ineq(pivots(i)) = 0, i=1,2,..k
405: -- put("and after elimination : "); put(ineq); new_line;
406: if Is_Zero(ineq)
407: then
408: if --not In_Face(k,face,pts(l).all,pts)
409: --and then
410: l < face(face'last) -- lexicographic enumeration
411: and then (Pivot_after_Elimination
412: (l,1,pts(l).all,face,pivots,elim) /= 0)
413: then --put_line("Degenerate choice");
414: degenerate := true; return;
415: end if;
416: else
417: for ll in (k+1)..n loop -- fill in the column
418: tab(ll-k,column) := ineq(pivots(ll));
419: end loop;
420: column := column + 1;
421: end if;
422: end if;
423: end loop;
424: for l in tab'first(1)..(tab'last(1)-1) loop -- right hand side
425: tab(l,tab'last(2)) := 0;
426: end loop;
427: tab(tab'last(1),tab'last(2)) := -1;
428: lastcol := column-1;
429: end Create_Tableau_for_Lower_Faces;
430:
431: -- DETERMINE WHETHER THE CANDIDATE IS VERTEX, SPANS AN EDGE OR A K-FACE :
432:
433: function Is_Vertex ( i,m,n : integer; pts : VecVec ) return boolean is
434:
435: tableau : matrix(1..n+1,1..m);
436: feasible : boolean;
437:
438: begin
439: Create_Tableau_for_Vertices(i,n,pts,tableau);
440: -- put_line("The tableau :"); put(tableau);
441: Integer_Complementary_Slackness(tableau,feasible);
442: return not feasible;
443: end Is_Vertex;
444:
445: function Is_Edge ( i,j,m,n : integer; pts : VecVec ) return boolean is
446:
447: elim : Vector(1..n) := pts(i).all - pts(j).all;
448: pivot : integer;
449:
450: begin
451: pivot := elim'first - 1;
452: for k in elim'range loop
453: if elim(k) /= 0
454: then pivot := k;
455: end if;
456: exit when pivot >= elim'first;
457: end loop;
458: if pivot < elim'first
459: then return false;
460: else Scale(elim);
461: declare
462: tab : matrix(1..n,1..m-1);
463: deg,feasible : boolean;
464: lst : integer;
465: begin
466: -- put("The elimination vector : "); put(elim); new_line;
467: Create_Tableau_for_Edges(i,j,n,pivot,elim,pts,tab,lst,deg);
468: -- put_line("The tableau :"); put(tab);
469: if deg
470: then return false;
471: elsif lst = 0
472: then return true;
473: else Integer_Complementary_Slackness(tab,lst,feasible);
474: return not feasible;
475: end if;
476: end;
477: end if;
478: end Is_Edge;
479:
480: function Is_Face ( k,n,m : natural; elim,pts : VecVec;
481: pivots,face : Vector ) return boolean is
482:
483: -- DESCRIPTION :
484: -- Applies Complementary Slackness to determine whether the given
485: -- candidate face is a face of the polytope.
486:
487: -- ON ENTRY :
488: -- k dimension of the candidate face;
489: -- n dimension of the vector space;
490: -- m number of points which span the polytope;
491: -- elim elimination equations in upper triangular form;
492: -- pts the points which span the polytope;
493: -- pivots pivoting information for the elimination equations;
494: -- face entries of the points which span the candidate face.
495:
496: nn : constant natural := n-k+1;
497: tab : matrix(1..nn,1..(m-k));
498: deg,feasible : boolean;
499: lst : integer;
500:
501: begin
502: -- put("The face being tested : "); put(face); new_line;
503: -- put("The pivots : "); put(pivots); new_line;
504: -- put_line("The elimination equations : ");
505: -- for i in 1..elim'last loop
506: -- put(elim(i)); put(" = 0 ");
507: -- end loop;
508: -- new_line;
509: if m - k <= 1
510: then return true;
511: else
512: Create_Tableau_for_Faces(k,n,face,pivots,pts,elim,tab,lst,deg);
513: -- put_line("The tableau of inequalities : "); put(tab);
514: if deg
515: then -- put_line("Tableau is degenerate: no solution :");
516: return false;
517: elsif lst = 0
518: then return true;
519: else Integer_Complementary_Slackness(tab,lst,feasible);
520: -- if feasible
521: -- then put_line(" is feasible");
522: -- else put_line(" is not feasible");
523: -- end if;
524: return not feasible;
525: end if;
526: end if;
527: end Is_Face;
528:
529: function Is_Face_of_Sum
530: ( k,n,m,i,r : natural; elim,pts,face : VecVec;
531: ind,pivots : Vector ) return boolean is
532:
533: -- DESCRIPTION :
534: -- Applies Complementary Slackness to determine whether the given
535: -- candidate face is a face of the polytope.
536:
537: -- ON ENTRY :
538: -- k dimension of the candidate face;
539: -- n dimension of the vector space;
540: -- m number of total points which span the polytope;
541: -- i current polytope;
542: -- r number of different polytopes;
543: -- elim elimination equations in upper triangular form;
544: -- pts the points which span the polytope;
545: -- face entries of the points which span the candidate face;
546: -- ind indicates the starting vector in pts of each polytope;
547: -- pivots pivoting information for the elimination equations.
548:
549: nn : constant natural := n-k+1;
550: tab : matrix(1..nn,1..(m-k));
551: deg,feasible : boolean;
552: lst : integer;
553:
554: begin
555: if m - k <= 1
556: then return true;
557: else
558: Create_Tableau_for_Faces_of_Sum
559: (k,n,i,r,ind,pivots,pts,elim,face,tab,lst,deg);
560: -- put_line("The tableau of inequalities : "); put(tab);
561: if deg
562: then --put_line("Tableau is degenerate: no solution :");
563: return false;
564: elsif lst = 0
565: then return true;
566: else Integer_Complementary_Slackness(tab,lst,feasible);
567: --if feasible
568: -- then put_line(" is feasible");
569: -- else put_line(" is not feasible");
570: --end if;
571: return not feasible;
572: end if;
573: end if;
574: end Is_Face_of_Sum;
575:
576: function Is_Lower_Edge ( i,j,m,n : integer; pts : VecVec )
577: return boolean is
578:
579: elim : Vector(1..n) := pts(i).all - pts(j).all;
580: pivot : integer;
581:
582: begin
583: pivot := elim'first - 1;
584: for k in elim'range loop
585: if elim(k) /= 0
586: then pivot := k;
587: end if;
588: exit when pivot >= elim'first;
589: end loop;
590: if pivot < elim'first or else (pivot = elim'last)
591: then return false;
592: else Scale(elim);
593: declare
594: tab : matrix(1..n-1,1..m-1);
595: deg,feasible : boolean;
596: lst : integer;
597: begin
598: Create_Tableau_for_Lower_Edges(i,j,n,pivot,elim,pts,tab,lst,deg);
599: -- put_line("The tableau :"); put(tab);
600: if deg
601: then return false;
602: elsif lst = 0
603: then return true;
604: else Integer_Complementary_Slackness(tab,lst,feasible);
605: return not feasible;
606: end if;
607: end;
608: end if;
609: end Is_Lower_Edge;
610:
611: function Is_Lower_Face
612: ( k,n,m : in natural; elim,pts : VecVec;
613: pivots,face : Vector ) return boolean is
614:
615: -- DESCRIPTION :
616: -- Applies Complementary Slackness to determine whether the given
617: -- candidate face is a face of the lower hull of the polytope.
618:
619: -- ON ENTRY :
620: -- k dimension of the candidate face;
621: -- n dimension of the vector space;
622: -- m number of points which span the polytope;
623: -- elim elimination equations in upper triangular form;
624: -- pts the points which span the polytope;
625: -- pivots pivoting information for the elimination equations;
626: -- face entries of the points which span the candidate face.
627:
628: nn : constant natural := n-k;
629: tab : matrix(1..nn,1..(m-k));
630: deg,feasible : boolean;
631: lst : integer;
632:
633: begin
634: -- put("The pivots : "); put(pivots); new_line;
635: -- put_line("The elimination equations : ");
636: -- for i in 1..elim'last loop
637: -- put(elim(i)); put(" = 0 ");
638: -- end loop;
639: -- new_line;
640: if m - k <= 1
641: then return true;
642: else
643: Create_Tableau_for_Lower_Faces(k,n,face,pivots,pts,elim,tab,lst,deg);
644: -- put_line("The tableau of inequalities : "); put(tab);
645: if deg
646: then --put_line("Degenerate tableau");
647: return false;
648: elsif lst = 0
649: then --put_line("lst = 0");
650: return true;
651: else Integer_Complementary_Slackness(tab,lst,feasible);
652: --if feasible
653: -- then put_line(" is feasible");
654: -- else put_line(" is not feasible");
655: --end if;
656: return not feasible;
657: end if;
658: end if;
659: end Is_Lower_Face;
660:
661: -- TARGET ROUTINES :
662:
663: procedure Enumerate_Vertices ( pts : in VecVec ) is
664:
665: continue : boolean := true;
666: n : constant natural := pts(pts'first).all'length;
667: m : constant natural := pts'length;
668:
669: begin
670: for i in pts'range loop
671: if Is_Vertex(i,m,n,pts)
672: then Process(i,continue);
673: end if;
674: exit when not continue;
675: end loop;
676: end Enumerate_Vertices;
677:
678: procedure Enumerate_Edges ( pts : in VecVec ) is
679:
680: continue : boolean := true;
681: n : constant natural := pts(pts'first).all'length;
682: m : constant natural := pts'length;
683:
684: procedure Candidate_Edges ( i,n : in integer ) is
685: begin
686: for j in (i+1)..pts'last loop -- enumerate all candidates
687: -- put("Verifying "); put(pts(i).all); put(" and");
688: -- put(pts(j).all); put(" :"); new_line;
689: if Is_Edge(i,j,m,n,pts)
690: then Process(i,j,continue);
691: end if;
692: exit when not continue;
693: end loop;
694: end Candidate_Edges;
695:
696: begin
697: for i in pts'first..(pts'last-1) loop
698: Candidate_Edges(i,n);
699: end loop;
700: end Enumerate_Edges;
701:
702: procedure Enumerate_Lower_Edges ( pts : in VecVec ) is
703:
704: continue : boolean := true;
705: n : constant natural := pts(pts'first).all'length;
706: m : constant natural := pts'length;
707:
708: procedure Candidate_Lower_Edges ( i,n : in integer ) is
709: begin
710: for j in (i+1)..pts'last loop -- enumerate all candidates
711: -- put("Verifying "); put(pts(i).all); put(" and");
712: -- put(pts(j).all); put(" :"); new_line;
713: if Is_Lower_Edge(i,j,m,n,pts)
714: then Process(i,j,continue);
715: end if;
716: exit when not continue;
717: end loop;
718: end Candidate_Lower_Edges;
719:
720: begin
721: for i in pts'first..(pts'last-1) loop
722: Candidate_Lower_Edges(i,n);
723: end loop;
724: end Enumerate_Lower_Edges;
725:
726: procedure Enumerate_Faces ( k : in natural; pts : in VecVec ) is
727:
728: m : constant natural := pts'length;
729: n : constant natural := pts(pts'first).all'length;
730: candidate : Vector(0..k);
731: elim : VecVec(0..k);
732: continue : boolean := true;
733:
734: procedure Candidate_Faces ( ipvt : in Vector; i,l : in integer ) is
735:
736: -- DESCRIPTION :
737: -- Enumerates all candidate k-faces in lexicographic order.
738:
739: piv : natural;
740: pivots : Vector(1..n);
741:
742: begin
743: if l > k
744: then if (k = m)
745: or else Is_Face(k,n,m,elim,pts,ipvt,candidate)
746: then Process(candidate,continue);
747: end if;
748: else for j in (i+1)..pts'last loop
749: candidate(l) := j;
750: pivots := ipvt;
751: Update_Eliminator(elim,l,pivots,pts(j).all,piv);
752: if piv /= 0
753: then Candidate_Faces(pivots,j,l+1);
754: end if;
755: Clear(elim(l));
756: exit when not continue;
757: end loop;
758: end if;
759: end Candidate_Faces;
760:
761: begin
762: if k <= m
763: then
764: declare
765: ipvt : Vector(1..n);
766: begin
767: for i in ipvt'range loop
768: ipvt(i) := i;
769: end loop;
770: for i in pts'first..(pts'last-k) loop
771: candidate(0) := i;
772: elim(0) := pts(i); -- basis point
773: Candidate_Faces(ipvt,i,1);
774: end loop;
775: end;
776: end if;
777: end Enumerate_Faces;
778:
779: procedure Enumerate_Lower_Faces ( k : in natural; pts : VecVec ) is
780:
781: m : constant natural := pts'length;
782: n : constant natural := pts(pts'first).all'length;
783: candidate : Vector(0..k);
784: elim : VecVec(0..k);
785: continue : boolean := true;
786:
787: procedure Candidate_Lower_Faces ( ipvt : in Vector; i,l : in integer ) is
788:
789: -- DESCRIPTION :
790: -- Enumerates all candidate k-faces in lexicographic order.
791:
792: piv : natural;
793: pivots : Vector(1..n);
794:
795: begin
796: if l > k
797: then --put_line("Testing the following candidate face :");
798: --for ii in candidate'first..candidate'last-1 loop
799: -- put(pts(candidate(ii))); put(" & ");
800: --end loop;
801: --put(pts(candidate(candidate'last))); new_line;
802: if (k = m) or else Is_Lower_Face(k,n,m,elim,pts,ipvt,candidate)
803: then Process(candidate,continue);
804: end if;
805: else for j in (i+1)..pts'last loop
806: candidate(l) := j;
807: pivots := ipvt;
808: -- put("Picking "); put(pts(j));
809: -- put(" Pivots : "); put(pivots); new_line;
810: Update_Eliminator(elim,l,pivots,pts(j).all,piv);
811: -- put(" update of eliminator piv = "); put(piv,1);
812: -- put(" Pivots : "); put(pivots); new_line;
813: if (piv /= 0) and (piv /= n)
814: then Candidate_Lower_Faces(pivots,j,l+1);
815: end if;
816: Clear(elim(l));
817: exit when not continue;
818: end loop;
819: end if;
820: end Candidate_Lower_Faces;
821:
822: begin
823: if k <= m
824: then
825: declare
826: ipvt : Vector(1..n);
827: begin
828: for i in ipvt'range loop
829: ipvt(i) := i;
830: end loop;
831: for i in pts'first..(pts'last-k) loop
832: candidate(0) := i;
833: elim(0) := pts(i); -- basis point
834: Candidate_Lower_Faces(ipvt,i,1);
835: end loop;
836: end;
837: end if;
838: end Enumerate_Lower_Faces;
839:
840: procedure Enumerate_Faces_of_Sum
841: ( ind,typ : in Vector; k : in natural; pts : in VecVec ) is
842:
843: m : constant natural := pts'length; -- number of points
844: n : constant natural := pts(pts'first)'length; -- dimension of vectors
845: r : constant natural := ind'length; -- number of polytopes
846: candidates : VecVec(1..r);
847: elim : VecVec(1-r..k);
848: pivots : Vector(1..n);
849: continue : boolean := true;
850: lasti,sum : natural;
851:
852: procedure Candidate_Faces_of_Sum ( ipvt : in Vector; i : in integer ) is
853:
854: -- DESCRIPTION :
855: -- Enumerates all faces of the given type on the sum,
856: -- i indicates the current polytope.
857:
858: procedure Candidate_Faces ( ipvt : in Vector;
859: start,l : in integer ) is
860:
861: -- DESCRIPTION :
862: -- Enumerates all candidate k-faces, with k = typ(i).
863: -- The parameter l indicates the current element of pts to be chosen.
864: -- The previously chosen point is given by start.
865:
866: piv,last : natural;
867:
868: begin
869: if i = r
870: then last := m;
871: else last := ind(i+1)-1;
872: end if;
873: if l > typ(i)
874: then
875: if (typ(i) = last-ind(i)+1)
876: or else Is_Face_of_Sum
877: (sum,n,last-i+1,i,r,elim,pts(pts'first..last),
878: candidates,ind,ipvt)
879: then
880: if i = r
881: then Process(candidates,continue);
882: else Candidate_Faces_of_Sum(ipvt,i+1);
883: end if;
884: end if;
885: else
886: for j in (start+1)..(last-typ(i)+l) loop
887: candidates(i)(l) := j;
888: if l = 0
889: then
890: Candidate_Faces(ipvt,j,l+1);
891: else
892: pivots := ipvt;
893: Update_Eliminator_for_Sum
894: (elim,sum-typ(i)+l,pivots,pts(j).all,i,piv);
895: if piv /= 0
896: then Candidate_Faces(pivots,j,l+1);
897: end if;
898: Clear(elim(sum-typ(i)+l));
899: end if;
900: exit when not continue;
901: end loop;
902: end if;
903: end Candidate_Faces;
904:
905: begin
906: candidates(i) := new Vector(0..typ(i));
907: if i = r
908: then lasti := pts'last;
909: else lasti := ind(i+1)-1;
910: end if;
911: sum := sum + typ(i);
912: for j in ind(i)..(lasti-typ(i)) loop
913: candidates(i)(0) := j;
914: elim(1-i) := pts(j);
915: Candidate_Faces(ipvt,j,1);
916: end loop;
917: sum := sum - typ(i);
918: -- for j in (sum+1)..(sum+typ(i)) loop
919: -- Clear(elim(j));
920: -- end loop;
921: Clear(candidates(i));
922: end Candidate_Faces_of_Sum;
923:
924: begin
925: declare
926: ipvt : Vector(1..n);
927: begin
928: for i in ipvt'range loop
929: ipvt(i) := i;
930: end loop;
931: sum := 0;
932: Candidate_Faces_of_Sum(ipvt,1);
933: end;
934: end Enumerate_Faces_of_Sum;
935:
936: end Integer_Face_Enumerators;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to integer_face_enumerators.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports