Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/floating_pruning_methods.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Natural_Vectors;
! 2: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
! 3: with Standard_Floating_Linear_Solvers; use Standard_Floating_Linear_Solvers;
! 4: with Floating_Linear_Inequalities; use Floating_Linear_Inequalities;
! 5: with Lists_of_Floating_Vectors; use Lists_of_Floating_Vectors;
! 6:
! 7: package body Floating_Pruning_Methods is
! 8:
! 9: -- GENERAL AUXILIARIES :
! 10:
! 11: procedure Normalize ( v : in out Standard_Floating_Vectors.Vector ) is
! 12:
! 13: -- DESCRIPTION : Divides every entry by v(v'last).
! 14:
! 15: begin
! 16: for i in v'range loop
! 17: v(i) := v(i)/v(v'last);
! 18: end loop;
! 19: end Normalize;
! 20:
! 21: function Convert ( fa : Face_Array ) return Array_of_Lists is
! 22:
! 23: -- DESCRIPTION :
! 24: -- Converts the array of faces into an array of lists by
! 25: -- converting the first element of each list of faces.
! 26:
! 27: res : Array_of_Lists(fa'range);
! 28:
! 29: begin
! 30: for k in fa'range loop
! 31: res(k) := Shallow_Create(fa(k).all);
! 32: end loop;
! 33: return res;
! 34: end Convert;
! 35:
! 36: -- AUXILIARIES FOR THE PRUNING ALGORITHMS :
! 37:
! 38: procedure Initialize ( n : in natural; mat : out Matrix;
! 39: ipvt : out Standard_Natural_Vectors.Vector ) is
! 40:
! 41: -- DESCRIPTION :
! 42: -- Initializes an n*(n+1) matrix with zeroes and the pivoting vector
! 43: -- with the entries 1..n.
! 44:
! 45: begin
! 46: for i in 1..n loop
! 47: for j in 1..n+1 loop
! 48: mat(i,j) := 0.0;
! 49: end loop;
! 50: end loop;
! 51: for i in 1..n loop
! 52: ipvt(i) := i;
! 53: end loop;
! 54: end Initialize;
! 55:
! 56: function Number_of_Inequalities
! 57: ( mix : Vector; lifted : Array_of_Lists ) return natural is
! 58:
! 59: -- DESCRIPTION :
! 60: -- Returns the maximal number of inequalities for pruning.
! 61:
! 62: res : natural := 0;
! 63:
! 64: begin
! 65: for k in lifted'range loop
! 66: res := res + Length_Of(lifted(k)) - mix(k) - 1;
! 67: end loop;
! 68: return res;
! 69: end Number_of_Inequalities;
! 70:
! 71: procedure Ordered_Inequalities ( k : in natural; mat : in out Matrix ) is
! 72:
! 73: -- DESCRIPTION :
! 74: -- Defines k inequalities mat(k,k) - mat(k+1,k) >= 0.
! 75:
! 76: begin
! 77: for i in mat'first(1)..k loop
! 78: for j in mat'range(2) loop
! 79: mat(i,j) := 0.0;
! 80: end loop;
! 81: mat(i,i) := 1.0; mat(i,i+1) := -1.0;
! 82: end loop;
! 83: end Ordered_Inequalities;
! 84:
! 85: procedure Check_and_Update
! 86: ( mic : in Face_Array; lifted : in Array_of_Lists;
! 87: m : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
! 88: tol : in double_float;
! 89: mixsub,mixsub_last : in out Mixed_Subdivision ) is
! 90:
! 91: -- DESCRIPTION :
! 92: -- Computes the normal to the points in pts, by solving the
! 93: -- linear system defined by m and ipvt.
! 94: -- If the computed normal is an inner normal w.r.t. the lifted points,
! 95: -- then the mixed subdivision will be updated with a new cell.
! 96:
! 97: use Standard_Floating_Vectors;
! 98: v : Standard_Floating_Vectors.Vector(m'range(2)) := Solve(m,tol,ipvt);
! 99: pts : Array_of_Lists(mic'range);
! 100:
! 101: begin
! 102: if abs(v(v'last)) > tol
! 103: then Normalize(v);
! 104: if v(v'last) < 0.0
! 105: then Min(v);
! 106: end if;
! 107: -- pts := Convert(mic);
! 108: -- Update(pts,v,mixsub,mixsub_last);
! 109: declare
! 110: cell : Mixed_Cell;
! 111: begin
! 112: cell.nor := new Standard_Floating_Vectors.Vector'(v);
! 113: cell.pts := new Array_of_Lists'(Convert(mic));
! 114: cell.sub := null;
! 115: Append(mixsub,mixsub_last,cell);
! 116: end;
! 117: end if;
! 118: end Check_and_Update;
! 119:
! 120: procedure Create_Equalities
! 121: ( n : in natural; f : in Face; mat,ineq : in Matrix;
! 122: tol : in double_float;
! 123: ipvt : in Standard_Natural_Vectors.Vector;
! 124: row,rowineq : in natural;
! 125: newmat,newineq : in out Matrix;
! 126: newipvt : in out Standard_Natural_Vectors.Vector;
! 127: newrow,newrowineq : in out natural; fail : out boolean ) is
! 128:
! 129: -- DESCRIPTION :
! 130: -- Creates new equalities and uses them to eliminate unknowns in the
! 131: -- matrices of equalities and inequalities. Failure is reported when
! 132: -- a zero row is encountered. On entry, all new* variables must be
! 133: -- initialized with the corresponding *-ones.
! 134:
! 135: shi : Standard_Floating_Vectors.Vector(1..n+1) := f(f'first).all;
! 136: fl : boolean := false;
! 137: pivot : natural;
! 138:
! 139: begin
! 140: for i in f'first+1..f'last loop
! 141: newrow := newrow + 1;
! 142: for j in f(i)'range loop
! 143: newmat(newrow,j) := f(i)(j) - shi(j);
! 144: end loop;
! 145: Switch(newipvt,newrow,newmat);
! 146: Upper_Triangulate(newrow,newmat,tol,newipvt,pivot);
! 147: fl := (pivot = 0);
! 148: exit when fl;
! 149: Switch(newrow,pivot,ineq'first,rowineq,newineq);
! 150: end loop;
! 151: fail := fl;
! 152: end Create_Equalities;
! 153:
! 154: procedure Complementary_Slackness
! 155: ( tableau : in out matrix;
! 156: tol : in double_float; feasible : out boolean ) is
! 157:
! 158: lastcol : constant integer := tableau'last(2)-1;
! 159: rhs,sol : Standard_Floating_Vectors.Vector(tableau'range(1));
! 160: columns : Vector(sol'range);
! 161:
! 162: begin
! 163: for i in rhs'range loop
! 164: rhs(i) := double_float(tableau(i,tableau'last(2)));
! 165: end loop;
! 166: Complementary_Slackness(tableau,lastcol,rhs,tol,sol,columns,feasible);
! 167: end Complementary_Slackness;
! 168:
! 169: function Check_Feasibility ( k,m,n : natural; ineq : Matrix;
! 170: tol : double_float ) return boolean is
! 171:
! 172: -- DESCRIPTION :
! 173: -- Returns true if -v(n+1) > 0 can be derived, false otherwise.
! 174:
! 175: -- ON ENTRY :
! 176: -- k current unknown that has been eliminated,
! 177: -- for all i <= k : ineq(l,i) = 0, for l in ineq'first..m;
! 178: -- m number of inequalities;
! 179: -- n dimension;
! 180: -- ineq matrix of inequalities.
! 181:
! 182: tableau : Matrix(1..n-k+1,1..m+1);
! 183: feasi : boolean;
! 184:
! 185: begin
! 186: if m = 0
! 187: then feasi := false;
! 188: else for i in k+1..n+1 loop
! 189: for j in 1..m loop
! 190: tableau(i-k,j) := ineq(j,i);
! 191: end loop;
! 192: tableau(i-k,m+1) := 0.0;
! 193: end loop;
! 194: tableau(n-k+1,m+1) := -1.0;
! 195: Complementary_Slackness(tableau,tol,feasi);
! 196: end if;
! 197: return feasi;
! 198: end Check_Feasibility;
! 199:
! 200: procedure Update_Inequalities
! 201: ( k,rowmat1,rowmat2,n : in natural;
! 202: mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
! 203: tol : in double_float;
! 204: rowineq : in out natural; ineq : in out Matrix;
! 205: lifted : in Array_of_Lists; mic : in out Face_Array ) is
! 206:
! 207: -- DESCRIPTION :
! 208: -- The inequalities will be updated w.r.t. the equality
! 209: -- constraints on the inner normal.
! 210:
! 211: tmp : List;
! 212: pt,shi : Standard_Floating_Vectors.Link_to_Vector;
! 213:
! 214: begin
! 215: for i in ineq'first..rowineq loop -- update the old inequalities
! 216: for j in rowmat1..rowmat2 loop
! 217: Upper_Triangulate(j,mat,tol,i,ineq);
! 218: end loop;
! 219: end loop;
! 220: shi := mic(k)(mic(k)'first);
! 221: tmp := lifted(k); -- make new inequalities
! 222: while not Is_Null(tmp) loop
! 223: pt := Head_Of(tmp);
! 224: if not Is_In(mic(k),pt.all)
! 225: then rowineq := rowineq + 1;
! 226: for j in pt'range loop
! 227: ineq(rowineq,j) := pt(j) - shi(j);
! 228: end loop;
! 229: Switch(ipvt,rowineq,ineq);
! 230: for i in 1..rowmat2 loop
! 231: Upper_Triangulate(i,mat,tol,rowineq,ineq);
! 232: end loop;
! 233: end if;
! 234: tmp := Tail_Of(tmp);
! 235: end loop;
! 236: end Update_Inequalities;
! 237:
! 238: -- CONSTRUCTION WITH PRUNING :
! 239:
! 240: procedure Gen1_Create
! 241: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
! 242: lifted : in Array_of_Lists; tol : in double_float;
! 243: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
! 244: mixsub : out Mixed_Subdivision ) is
! 245:
! 246: res,res_last : Mixed_Subdivision;
! 247: accu : Face_Array(fa'range);
! 248: ma : Matrix(1..n,1..n+1);
! 249: ipvt : Standard_Natural_Vectors.Vector(1..n);
! 250: ineqrows : natural;
! 251:
! 252: procedure Compute_Mixed_Cells
! 253: ( k,row : in natural; mat : in Matrix;
! 254: ipvt : in Standard_Natural_Vectors.Vector;
! 255: rowineq : in natural; ineq : in Matrix;
! 256: mic : in out Face_Array; continue : out boolean );
! 257:
! 258: -- DESCRIPTION :
! 259: -- Backtrack recursive procedure to enumerate face-face combinations.
! 260:
! 261: -- ON ENTRY :
! 262: -- k index for current point set;
! 263: -- row number of rows already in the matrix;
! 264: -- mat matrix which determines the inner normal;
! 265: -- ipvt contains the pivoting information;
! 266: -- rowineq number of inequality constraints already in ineq;
! 267: -- ineq matrix for the inequality constraints on the
! 268: -- inner normal v, of type <.,v> >= 0;
! 269: -- mic contains the current selected faces, up to k-1.
! 270:
! 271: -- ON RETURN :
! 272: -- mic updated selected faces.
! 273: -- continue indicates whether to continue the creation or not.
! 274:
! 275: procedure Process_Inequalities
! 276: ( k,rowmat1,rowmat2 : in natural;
! 277: mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
! 278: rowineq : in out natural; ineq : in out Matrix;
! 279: mic : in out Face_Array; cont : out boolean ) is
! 280:
! 281: -- DESCRIPTION :
! 282: -- Updates inequalities and checks feasibility before proceeding.
! 283:
! 284: fl : boolean := false;
! 285: tmp : List;
! 286: pt,shi : Link_to_Vector;
! 287:
! 288: begin
! 289: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,tol,rowineq,ineq,
! 290: lifted,mic);
! 291: if Check_Feasibility(rowmat2,rowineq,n,ineq,tol)
! 292: then nbfail(k) := nbfail(k) + 1.0;
! 293: cont := true;
! 294: else nbsucc(k) := nbsucc(k) + 1.0;
! 295: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic,cont);
! 296: end if;
! 297: end Process_Inequalities;
! 298:
! 299: procedure Compute_Mixed_Cells
! 300: ( k,row : in natural; mat : in Matrix;
! 301: ipvt : in Standard_Natural_Vectors.Vector;
! 302: rowineq : in natural; ineq : in Matrix;
! 303: mic : in out Face_Array; continue : out boolean ) is
! 304:
! 305: old : Mixed_Subdivision := res_last;
! 306: cont : boolean := true;
! 307: tmpfa : Faces;
! 308:
! 309: begin
! 310: if k > mic'last
! 311: then Check_and_Update(mic,lifted,mat,ipvt,tol,res,res_last);
! 312: if old /= res_last
! 313: then Process(Head_Of(res_last),continue);
! 314: else continue := true;
! 315: end if;
! 316: else tmpfa := fa(k);
! 317: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
! 318: mic(k) := Head_Of(tmpfa);
! 319: declare -- update the matrices
! 320: fl : boolean;
! 321: newipvt : Standard_Natural_Vectors.Vector(ipvt'range) := ipvt;
! 322: newmat : Matrix(mat'range(1),mat'range(2)) := mat;
! 323: newineq : Matrix(ineq'range(1),ineq'range(2)) := ineq;
! 324: newrow : natural := row;
! 325: newrowineq : natural := rowineq;
! 326: begin
! 327: Create_Equalities
! 328: (n,mic(k),mat,ineq,tol,ipvt,row,rowineq,newmat,newineq,
! 329: newipvt,newrow,newrowineq,fl);
! 330: if fl
! 331: then nbfail(k) := nbfail(k) + 1.0;
! 332: else Process_Inequalities(k,row+1,newrow,newmat,newipvt,
! 333: newrowineq,newineq,mic,cont);
! 334: end if;
! 335: end;
! 336: tmpfa := Tail_Of(tmpfa);
! 337: exit when not cont;
! 338: end loop;
! 339: continue := cont;
! 340: end if;
! 341: end Compute_Mixed_Cells;
! 342:
! 343: begin
! 344: Initialize(n,ma,ipvt);
! 345: ineqrows := Number_of_Inequalities(mix,lifted);
! 346: declare
! 347: ineq : matrix(1..ineqrows,1..n+1);
! 348: cont : boolean;
! 349: begin
! 350: ineq(1,1) := 0.0;
! 351: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu,cont);
! 352: end;
! 353: mixsub := res;
! 354: end Gen1_Create;
! 355:
! 356: procedure Create
! 357: ( n : in natural; mix : in Vector; fa : in Array_of_Faces;
! 358: lifted : in Array_of_Lists; tol : in double_float;
! 359: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
! 360: mixsub : out Mixed_Subdivision ) is
! 361:
! 362: res,res_last : Mixed_Subdivision;
! 363: accu : Face_Array(fa'range);
! 364: ma : Matrix(1..n,1..n+1);
! 365: ipvt : Standard_Natural_Vectors.Vector(1..n);
! 366: ineqrows : natural;
! 367:
! 368: procedure Compute_Mixed_Cells
! 369: ( k,row : in natural; mat : in Matrix;
! 370: ipvt : in Standard_Natural_Vectors.Vector;
! 371: rowineq : in natural; ineq : in Matrix;
! 372: mic : in out Face_Array );
! 373:
! 374: -- DESCRIPTION :
! 375: -- Backtrack recursive procedure to enumerate face-face combinations.
! 376:
! 377: -- ON ENTRY :
! 378: -- k index for current point set;
! 379: -- row number of rows already in the matrix;
! 380: -- mat matrix which determines the inner normal;
! 381: -- ipvt contains the pivoting information;
! 382: -- rowineq number of inequality constraints already in ineq;
! 383: -- ineq matrix for the inequality constraints on the
! 384: -- inner normal v, of type <.,v> >= 0;
! 385: -- mic contains the current selected faces, up to k-1.
! 386:
! 387: -- ON RETURN :
! 388: -- mic updated selected faces.
! 389:
! 390: procedure Process_Inequalities
! 391: ( k,rowmat1,rowmat2 : in natural;
! 392: mat : in Matrix; ipvt : in Standard_Natural_Vectors.Vector;
! 393: rowineq : in out natural; ineq : in out matrix;
! 394: mic : in out Face_Array) is
! 395:
! 396: -- DESCRIPTION :
! 397: -- Updates inequalities and checks feasibility before proceeding.
! 398:
! 399: tmp : List;
! 400: pt,shi : Link_to_Vector;
! 401:
! 402: begin
! 403: Update_Inequalities(k,rowmat1,rowmat2,n,mat,ipvt,tol,rowineq,ineq,
! 404: lifted,mic);
! 405: if Check_Feasibility(rowmat2,rowineq,n,ineq,tol)
! 406: then nbfail(k) := nbfail(k) + 1.0;
! 407: else nbsucc(k) := nbsucc(k) + 1.0;
! 408: Compute_Mixed_Cells(k+1,rowmat2,mat,ipvt,rowineq,ineq,mic);
! 409: end if;
! 410: end Process_Inequalities;
! 411:
! 412: procedure Compute_Mixed_Cells
! 413: ( k,row : in natural; mat : in Matrix;
! 414: ipvt : in Standard_Natural_Vectors.Vector;
! 415: rowineq : in natural; ineq : in Matrix;
! 416: mic : in out Face_Array ) is
! 417:
! 418: tmpfa : Faces;
! 419:
! 420: begin
! 421: if k > mic'last
! 422: then Check_and_Update(mic,lifted,mat,ipvt,tol,res,res_last);
! 423: else tmpfa := fa(k);
! 424: while not Is_Null(tmpfa) loop -- enumerate faces of kth polytope
! 425: mic(k) := Head_Of(tmpfa);
! 426: declare -- update matrices
! 427: fl : boolean;
! 428: newipvt : Standard_Natural_Vectors.Vector(ipvt'range) := ipvt;
! 429: newmat : Matrix(mat'range(1),mat'range(2)) := mat;
! 430: newineq : Matrix(ineq'range(1),ineq'range(2)) := ineq;
! 431: newrow : natural := row;
! 432: newrowineq : natural := rowineq;
! 433: begin
! 434: Create_Equalities
! 435: (n,mic(k),mat,ineq,tol,ipvt,row,rowineq,newmat,newineq,
! 436: newipvt,newrow,newrowineq,fl);
! 437: if fl
! 438: then nbfail(k) := nbfail(k) + 1.0;
! 439: else Process_Inequalities
! 440: (k,row+1,newrow,newmat,newipvt,newrowineq,newineq,mic);
! 441: end if;
! 442: end;
! 443: tmpfa := Tail_Of(tmpfa);
! 444: end loop;
! 445: end if;
! 446: end Compute_Mixed_Cells;
! 447:
! 448: begin
! 449: Initialize(n,ma,ipvt);
! 450: ineqrows := Number_of_Inequalities(mix,lifted);
! 451: declare
! 452: ineq : Matrix(1..ineqrows,1..n+1);
! 453: begin
! 454: ineq(1,1) := 0.0;
! 455: Compute_Mixed_Cells(accu'first,0,ma,ipvt,0,ineq,accu);
! 456: end;
! 457: mixsub := res;
! 458: end Create;
! 459:
! 460: end Floating_Pruning_Methods;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to floating_pruning_methods.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Stalift