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