Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/triangulations.adb, Revision 1.1
1.1 ! maekawa 1: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
! 2: with Standard_Random_Numbers; use Standard_Random_Numbers;
! 3:
! 4: --with text_io,Simplices_io; use text_io,Simplices_io;
! 5: --with Integer_Vectors_io; use Integer_Vectors_io;
! 6: --with integer_io; use integer_io;
! 7:
! 8: package body Triangulations is
! 9:
! 10: -- AUXILIAIRY HASH TABLE FOR Update_One :
! 11:
! 12: type Matrix_of_Triangulations
! 13: is array ( integer range <>, integer range <> ) of Triangulation;
! 14:
! 15: type Hash_Table ( n : natural ) is
! 16: record
! 17: weight1,weight2 : Link_to_Vector;
! 18: data : Matrix_of_Triangulations(0..n,0..n);
! 19: end record;
! 20:
! 21: procedure Init ( ht : in out Hash_Table; first,last : in integer ) is
! 22:
! 23: -- Initializes the hash table with Null_Triangulation.
! 24: -- Determines the weights for the hash function.
! 25: -- The numbers first and last determine the length of the weight function.
! 26:
! 27: begin
! 28: ht.weight1 := new Vector(first..last);
! 29: ht.weight2 := new Vector(first..last);
! 30: for i in first..last loop
! 31: ht.weight1(i) := Random(-last,last);
! 32: ht.weight2(i) := Random(-last,last);
! 33: end loop;
! 34: ht.weight1(last) := 1; -- TO AVOID PROBLEMS WITH HIGH LIFTING
! 35: ht.weight2(last) := 1;
! 36: for i in ht.data'range(1) loop
! 37: for j in ht.data'range(2) loop
! 38: ht.data(i,j) := Null_Triangulation;
! 39: end loop;
! 40: end loop;
! 41: end Init;
! 42:
! 43: procedure Hash_Function ( ht : in Hash_Table; s : in Simplex;
! 44: i,j : out integer ) is
! 45:
! 46: -- DESCRIPTION :
! 47: -- Computes the entries in the hash table for the simplex s.
! 48:
! 49: r1 : constant integer := ht.weight1.all*Vertex(s,1);
! 50: r2 : constant integer := ht.weight2.all*Vertex(s,2);
! 51:
! 52: begin
! 53: i := r1 mod (ht.n + 1);
! 54: j := r2 mod (ht.n + 1);
! 55: end Hash_Function;
! 56:
! 57: procedure Update ( ht : in out Hash_Table; s : in Simplex ) is
! 58:
! 59: -- DESCRIPTION :
! 60: -- Updates the hash table with the given simplex s.
! 61:
! 62: i,j : integer;
! 63:
! 64: begin
! 65: Hash_Function(ht,s,i,j);
! 66: Construct(s,ht.data(i,j));
! 67: end Update;
! 68:
! 69: function Is_In ( ht : Hash_Table; s : Simplex ) return boolean is
! 70:
! 71: -- DESCRIPTION :
! 72: -- Checks whether the simplex belongs to the hash table.
! 73:
! 74: i,j : integer;
! 75:
! 76: begin
! 77: Hash_Function(ht,s,i,j);
! 78: return Is_In(ht.data(i,j),s);
! 79: end Is_In;
! 80:
! 81: -- procedure Write_Distribution ( ht : in Hash_Table ) is
! 82: --
! 83: -- -- DESCRIPTION :
! 84: -- -- Writes the lengths of the lists in ht.
! 85: --
! 86: -- begin
! 87: -- for i in ht.data'range(1) loop
! 88: -- for j in ht.data'range(2) loop
! 89: -- put(Length_Of(ht.data(i,j)),1); put(' ');
! 90: -- end loop;
! 91: -- new_line;
! 92: -- end loop;
! 93: -- end Write_Distribution;
! 94:
! 95: procedure Clear ( ht : in out Hash_Table ) is
! 96:
! 97: -- DESCRIPTION :
! 98: -- Clears the allocated memory space for the hash table,
! 99: -- only a shallow copy is performed.
! 100:
! 101: begin
! 102: Clear(ht.weight1);
! 103: Clear(ht.weight2);
! 104: for i in ht.data'range(1) loop
! 105: for j in ht.data'range(2) loop
! 106: Lists_of_Simplices.Clear(Lists_of_Simplices.List(ht.data(i,j)));
! 107: end loop;
! 108: end loop;
! 109: end Clear;
! 110:
! 111: -- CREATORS :
! 112:
! 113: function Create ( s : Simplex ) return Triangulation is
! 114:
! 115: res : Triangulation;
! 116:
! 117: begin
! 118: Construct(s,res);
! 119: return res;
! 120: end Create;
! 121:
! 122: function Is_Inner ( n : natural; s : Simplex ) return boolean is
! 123: begin
! 124: for k in 1..n loop
! 125: if Neighbor(s,k) = Null_Simplex
! 126: then return false;
! 127: end if;
! 128: end loop;
! 129: return true;
! 130: end Is_Inner;
! 131:
! 132: procedure Update ( t : in out Triangulation; s : in out Simplex;
! 133: x : in Link_to_Vector ) is
! 134:
! 135: pos : vector(0..x'last) := Position(s,x.all);
! 136:
! 137: begin
! 138: for k in x'range loop
! 139: if Neighbor(s,k) = Null_Simplex
! 140: then if pos(k-1)*pos(pos'last) > 0
! 141: then Update(s,x,k);
! 142: Construct(Neighbor(s,k),t);
! 143: end if;
! 144: end if;
! 145: end loop;
! 146: end Update;
! 147:
! 148: procedure Connect_and_Update
! 149: ( s : in out Simplex; t : in out Triangulation ) is
! 150:
! 151: -- DESCRIPTION :
! 152: -- Connects the simplex with each other simplex in the triangulation
! 153: -- and adds it to the list t.
! 154:
! 155: tmp : Triangulation := t;
! 156:
! 157: begin
! 158: while not Is_Null(tmp) loop
! 159: declare
! 160: s2 : Simplex := Head_Of(tmp);
! 161: begin
! 162: Connect(s,s2);
! 163: --Set_Head(tmp,s2);
! 164: end;
! 165: tmp := Tail_Of(tmp);
! 166: end loop;
! 167: Construct(s,t);
! 168: end Connect_and_Update;
! 169:
! 170: procedure Connect_and_Update
! 171: ( t1 : in Triangulation; t2 : in out Triangulation ) is
! 172:
! 173: -- DESCRIPTION :
! 174: -- Connects all simplices in t1 properly with each other
! 175: -- and adds them to the list t2.
! 176:
! 177: tmp1 : Triangulation := t1;
! 178: tmp2 : Triangulation;
! 179: s1,s2 : Simplex;
! 180:
! 181: begin
! 182: while not Is_Null(tmp1) loop
! 183: s1 := Head_Of(tmp1);
! 184: tmp2 := Tail_Of(tmp1);
! 185: while not Is_Null(tmp2) loop
! 186: s2 := Head_Of(tmp2);
! 187: Connect(s1,s2);
! 188: tmp2 := Tail_Of(tmp2);
! 189: end loop;
! 190: tmp1 := Tail_Of(tmp1);
! 191: Construct(s1,t2);
! 192: end loop;
! 193: end Connect_and_Update;
! 194:
! 195: procedure Update ( t : in out Triangulation; x : in Link_to_Vector;
! 196: newt : out Triangulation ) is
! 197:
! 198: tmp : Triangulation := t;
! 199: s : Simplex;
! 200: nwt : Triangulation; -- for the new simplices that will contain x
! 201:
! 202: begin
! 203: while not Is_Null(tmp) loop
! 204: s := Head_Of(tmp);
! 205: if not Is_Inner(x'last,s)
! 206: then Update(nwt,s,x);
! 207: end if;
! 208: tmp := Tail_Of(tmp);
! 209: end loop;
! 210: Connect_and_Update(nwt,t);
! 211: newt := nwt;
! 212: end Update;
! 213:
! 214: procedure Update ( t : in out Triangulation; x : in Link_to_Vector ) is
! 215:
! 216: nwt : Triangulation; -- for the new simplices that will contain x
! 217:
! 218: begin
! 219: Update(t,x,nwt);
! 220: -- SHALLOW CLEAR :
! 221: Lists_of_Simplices.Clear(Lists_of_Simplices.List(nwt));
! 222: end Update;
! 223:
! 224: procedure Collect_Vertices
! 225: ( s : in Simplex; l : in out List ) is
! 226:
! 227: -- DESCRIPTION :
! 228: -- Collects the vertices in s and puts them in the list l.
! 229: -- Auxiliary routine for Update_One.
! 230:
! 231: pts : constant VecVec := Vertices(s);
! 232:
! 233: begin
! 234: for k in pts'range loop
! 235: if not Is_In(l,pts(k).all)
! 236: then declare
! 237: newpt : Link_to_Vector := new Vector'(pts(k).all);
! 238: begin
! 239: Construct(newpt,l);
! 240: end;
! 241: end if;
! 242: end loop;
! 243: end Collect_Vertices;
! 244:
! 245: function Has_Vertices_in_List ( s : Simplex; l : List ) return boolean is
! 246:
! 247: -- DESCRIPTION :
! 248: -- Returns true if the simplex s has vertices in the list l.
! 249:
! 250: pts : constant VecVec := Vertices(s);
! 251:
! 252: begin
! 253: for k in pts'range loop
! 254: if Is_In(l,pts(k).all)
! 255: then return true;
! 256: end if;
! 257: end loop;
! 258: return false;
! 259: end Has_Vertices_in_List;
! 260:
! 261: procedure Update_One
! 262: ( t : in out Triangulation; x : in Link_to_Vector ) is
! 263:
! 264: s : Simplex := Head_Of(t);
! 265: news,nei : Simplex;
! 266: nwt,leaves : Triangulation;
! 267: root : Hash_Table(2*x'last-1);
! 268: pos : Vector(0..x'last);
! 269: border_vertices : List;
! 270:
! 271: procedure Process_Tree_of_Updates is
! 272:
! 273: -- DESCRIPTION :
! 274: -- Constructs a tree of adjencencies which consists of
! 275: -- connected simplices on the edge of the triangulation.
! 276:
! 277: tmp,newleaves : Triangulation;
! 278: si,neisi : Simplex;
! 279:
! 280: begin
! 281: -- INITIALIZATION :
! 282: --Construct(s,root);
! 283: Init(root,x'first,x'last);
! 284: Update(root,s);
! 285: for k in x'range loop
! 286: declare
! 287: nei : Simplex := Neighbor(s,k);
! 288: begin
! 289: if (nei /= Null_Simplex) and then not Is_Vertex(nei,x.all)
! 290: then Construct(nei,leaves);
! 291: end if;
! 292: end;
! 293: end loop;
! 294: -- PERFORM THE UPDATES :
! 295: loop
! 296: -- INVARIANT CONDITION :
! 297: -- the simplices considered have vertices on the border close to x
! 298: newleaves := Null_Triangulation;
! 299: tmp := leaves;
! 300: while not Is_Null(tmp) loop
! 301: si := Head_Of(tmp);
! 302: -- UPDATE root TO PREVENT TURNING BACK :
! 303: Update(root,si);
! 304: -- COMPUTE NEW SIMPLICES AND LEAVES :
! 305: if not Is_Inner(x'last,si)
! 306: then pos := Position(si,x.all);
! 307: for k in x'range loop
! 308: -- LOOK IN THE DIRECTION TOWARDS x :
! 309: if pos(k-1)*pos(pos'last) > 0
! 310: then neisi := Neighbor(si,k);
! 311: if neisi = Null_Simplex
! 312: then Update(si,x,k); -- NEW SIMPLEX
! 313: neisi := Neighbor(si,k);
! 314: Construct(neisi,nwt);
! 315: Collect_Vertices(neisi,border_vertices);
! 316: end if;
! 317: end if;
! 318: end loop;
! 319: end if;
! 320: for l in x'range loop -- GO FURTHER
! 321: neisi := Neighbor(si,l);
! 322: if (neisi /= Null_Simplex)
! 323: and then not Is_Vertex(neisi,x.all)
! 324: and then Has_Vertices_in_List(neisi,border_vertices)
! 325: and then not Is_In(leaves,neisi)
! 326: and then not Is_In(newleaves,neisi)
! 327: and then not Is_In(root,neisi)
! 328: then Construct(neisi,newleaves);
! 329: end if;
! 330: end loop;
! 331: tmp := Tail_Of(tmp);
! 332: end loop;
! 333: exit when (newleaves = Null_Triangulation);
! 334: Lists_of_Simplices.Clear(Lists_of_Simplices.List(leaves));
! 335: leaves := newleaves;
! 336: end loop;
! 337: -- SHALLOW CLEAR :
! 338: -- Lists_of_Simplices.Clear(Lists_of_Simplices.List(root));
! 339: -- put_line("Distribution of root : "); Write_Distribution(root);
! 340: Clear(root);
! 341: Lists_of_Simplices.Clear(Lists_of_Simplices.List(leaves));
! 342: end Process_Tree_of_Updates;
! 343:
! 344: begin
! 345: nei := s; -- USED TO SEE WHETHER s WILL CHANGE
! 346: pos := Position(s,x.all);
! 347: Update_One(s,x,pos,news);
! 348: Construct(news,nwt);
! 349: Collect_Vertices(news,border_vertices);
! 350: if s /= nei
! 351: then pos := Position(s,x.all); -- BECAUSE s IS CHANGED !!!
! 352: end if;
! 353: for k in x'range loop -- COMPUTE OTHER NEIGHBORS
! 354: if pos(k-1)*pos(pos'last) > 0
! 355: then nei := Neighbor(s,k);
! 356: if nei = Null_Simplex
! 357: then Update(s,x,k);
! 358: nei := Neighbor(s,k);
! 359: Construct(nei,nwt);
! 360: Collect_Vertices(nei,border_vertices);
! 361: end if;
! 362: end if;
! 363: end loop;
! 364: Process_Tree_of_Updates;
! 365: Clear(border_vertices);
! 366: Connect_and_Update(nwt,t);
! 367: Lists_of_Simplices.Clear(Lists_of_Simplices.List(nwt));
! 368: end Update_One;
! 369:
! 370: procedure Connect ( t : in out Triangulation ) is
! 371:
! 372: tmp1 : Triangulation := t;
! 373: tmp2 : Triangulation;
! 374: s1,s2 : Simplex;
! 375:
! 376: begin
! 377: while not Is_Null(tmp1) loop
! 378: s1 := Head_Of(tmp1);
! 379: tmp2 := Tail_Of(tmp1);
! 380: while not Is_Null(tmp2) loop
! 381: s2 := Head_Of(tmp2);
! 382: Connect(s1,s2);
! 383: --Set_Head(tmp2,s2);
! 384: tmp2 := Tail_Of(tmp2);
! 385: end loop;
! 386: --Set_Head(tmp1,s1);
! 387: tmp1 := Tail_Of(tmp1);
! 388: end loop;
! 389: end Connect;
! 390:
! 391: -- THE OPTIMAL DISCRETE CONSERVATIVE LIFTING FUNCTION :
! 392:
! 393: -- procedure Write_Neighbors ( s : in Simplex ) is
! 394: -- begin
! 395: -- put("The normal of simplex s : "); put(Normal(s)); new_line;
! 396: -- put_line(" with vertices : "); put(s);
! 397: -- for k in Normal(s)'range loop
! 398: -- if Neighbor(s,k) /= Null_Simplex
! 399: -- then put("normal of a not null neighbors of s :");
! 400: -- put(Normal(Neighbor(s,k))); new_line;
! 401: -- end if;
! 402: -- end loop;
! 403: -- end Write_Neighbors;
! 404:
! 405: -- procedure Write ( t : in Triangulation ) is
! 406: --
! 407: -- tmp : Triangulation := t;
! 408: --
! 409: -- begin
! 410: -- while not Is_Null(tmp) loop
! 411: -- Write_Neighbors(Head_Of(tmp));
! 412: -- tmp := Tail_Of(tmp);
! 413: -- end loop;
! 414: -- end Write;
! 415:
! 416: procedure Flatten ( t : in out Triangulation ) is
! 417:
! 418: tmp : Triangulation := t;
! 419: s : Simplex;
! 420:
! 421: -- IMPLEMENTATION REQUIREMENT :
! 422: -- Cells that are already flattened are grouped at the end of the list.
! 423:
! 424: begin
! 425: -- put_line("THE TRIANGULATION BEFORE LIFTING : "); Write(t);
! 426: while not Is_Null(tmp) loop
! 427: s := Head_Of(tmp);
! 428: exit when Is_Flat(s);
! 429: -- put_line("NEIGHBORS BEFORE FLATTENING : "); Write_Neighbors(s);
! 430: Flatten(s);
! 431: -- put_line("NEIGHBORS AFTER FLATTENING : "); Write_Neighbors(s);
! 432: --Set_Head(tmp,s);
! 433: tmp := Tail_Of(tmp);
! 434: end loop;
! 435: -- put_line("THE TRIANGULATION AFTER LIFTING : "); Write(t);
! 436: end Flatten;
! 437:
! 438: -- SELECTORS :
! 439:
! 440: function Is_Vertex ( t : Triangulation; x : Vector ) return boolean is
! 441:
! 442: tmp : Triangulation := t;
! 443:
! 444: begin
! 445: while not Is_Null(tmp) loop
! 446: if Is_Vertex(Head_Of(tmp),x)
! 447: then return true;
! 448: end if;
! 449: tmp := Tail_Of(tmp);
! 450: end loop;
! 451: return false;
! 452: end Is_Vertex;
! 453:
! 454: function Vertices ( t : Triangulation ) return List is
! 455:
! 456: -- DESCRIPTION :
! 457: -- Returns a list with all vertices in the simplices of t.
! 458:
! 459: res,res_last : List;
! 460: tmp : Triangulation := t;
! 461:
! 462: begin
! 463: res_last := res;
! 464: while not Is_Null(tmp) loop
! 465: declare
! 466: s : Simplex := Head_Of(tmp);
! 467: v : constant VecVec := Vertices(s);
! 468: begin
! 469: for i in v'range loop
! 470: if not Is_In(res,v(i))
! 471: then Append(res,res_last,v(i));
! 472: end if;
! 473: end loop;
! 474: end;
! 475: tmp := Tail_Of(tmp);
! 476: end loop;
! 477: return res;
! 478: end Vertices;
! 479:
! 480: function Vertices ( t : Triangulation; x : vector ) return List is
! 481:
! 482: res,res_last : List;
! 483: tmp : Triangulation := t;
! 484:
! 485: begin
! 486: res_last := res;
! 487: while not Is_Null(tmp) loop
! 488: declare
! 489: s : Simplex := Head_Of(tmp);
! 490: begin
! 491: if Is_Vertex(s,x)
! 492: then
! 493: declare
! 494: v : constant VecVec := Vertices(s);
! 495: begin
! 496: for i in v'range loop
! 497: if not Is_In(res,v(i))
! 498: then Append(res,res_last,v(i));
! 499: end if;
! 500: end loop;
! 501: end;
! 502: end if;
! 503: tmp := Tail_Of(tmp);
! 504: end;
! 505: end loop;
! 506: return res;
! 507: end Vertices;
! 508:
! 509: function Vertices ( t : Triangulation ) return VecVec is
! 510:
! 511: vertri : List := Vertices(t);
! 512: res : VecVec(1..Length_Of(vertri));
! 513: tmp : List := vertri;
! 514:
! 515: begin
! 516: for i in res'range loop
! 517: res(i) := Head_Of(tmp);
! 518: tmp := Tail_Of(tmp);
! 519: end loop;
! 520: Shallow_Clear(vertri);
! 521: -- Lists_of_Link_to_Integer_Vectors.Clear
! 522: -- (Lists_of_Link_to_Integer_Vectors.List(vertri));
! 523: return res;
! 524: end Vertices;
! 525:
! 526: function Vertices ( t : Triangulation; x : vector ) return VecVec is
! 527:
! 528: vertri : List := Vertices(t,x);
! 529: res : VecVec(1..Length_Of(vertri));
! 530: tmp : List := vertri;
! 531:
! 532: begin
! 533: for i in res'range loop
! 534: res(i) := Head_Of(tmp);
! 535: tmp := Tail_Of(tmp);
! 536: end loop;
! 537: Shallow_Clear(vertri);
! 538: -- Lists_of_Link_to_Integer_Vectors.Clear
! 539: -- (Lists_of_Link_to_Integer_Vectors.List(vertri));
! 540: return res;
! 541: end Vertices;
! 542:
! 543: function Is_In1 ( t : Triangulation; x : Vector ) return boolean is
! 544:
! 545: tmp : Triangulation := t;
! 546:
! 547: begin
! 548: while not Is_Null(tmp) loop
! 549: if Is_In(Head_Of(tmp),x)
! 550: then return true;
! 551: end if;
! 552: tmp := Tail_Of(tmp);
! 553: end loop;
! 554: return false;
! 555: end Is_In1;
! 556:
! 557: function Is_In ( t : Triangulation; x : Vector ) return boolean is
! 558: begin
! 559: return Is_In_All(Head_Of(t),x);
! 560: end Is_In;
! 561:
! 562: function Is_In ( t : Triangulation; x : Vector ) return Simplex is
! 563: begin
! 564: return Is_In_All(Head_Of(t),x);
! 565: end Is_In;
! 566:
! 567: function Is_In ( t : Triangulation; s : Simplex ) return boolean is
! 568:
! 569: tmp : Triangulation := t;
! 570:
! 571: begin
! 572: -- put("Length of triangulation : "); put(Length_Of(t),1); new_line;
! 573: while not Is_Null(tmp) loop
! 574: if Equal(s,Head_Of(tmp))
! 575: then return true;
! 576: else tmp := Tail_Of(tmp);
! 577: end if;
! 578: end loop;
! 579: return false;
! 580: end Is_In;
! 581:
! 582: function Volume ( t : Triangulation ) return natural is
! 583:
! 584: tmp : Triangulation := t;
! 585: vol : natural := 0;
! 586:
! 587: begin
! 588: while not Is_Null(tmp) loop
! 589: vol := vol + Volume(Head_Of(tmp));
! 590: tmp := Tail_Of(tmp);
! 591: end loop;
! 592: return vol;
! 593: end Volume;
! 594:
! 595: -- DESTRUCTOR :
! 596:
! 597: procedure Clear ( t : in out Triangulation ) is
! 598:
! 599: tmp : Triangulation := t;
! 600:
! 601: begin
! 602: while not Is_Null(tmp) loop
! 603: declare
! 604: s : Simplex := Head_Of(tmp);
! 605: begin
! 606: Clear(s);
! 607: end;
! 608: tmp := Tail_Of(tmp);
! 609: end loop;
! 610: Lists_of_Simplices.Clear(Lists_of_Simplices.List(t));
! 611: end Clear;
! 612:
! 613: end Triangulations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to triangulations.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 / Dynlift