Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/floating_face_enumerators.adb, Revision 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