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