Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/simplices.adb, Revision 1.1
1.1 ! maekawa 1: with unchecked_deallocation;
! 2: with Standard_Integer_Norms; use Standard_Integer_Norms;
! 3: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
! 4: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
! 5: with Transformations; use Transformations;
! 6:
! 7: package body Simplices is
! 8:
! 9: -- DATA STRUCTURES :
! 10:
! 11: type Point is record
! 12: pt : Link_to_Vector; -- the point
! 13: si : Simplex; -- by pivoting a new simplex is obtained
! 14: end record;
! 15: type Points is array ( integer range <> ) of Point;
! 16:
! 17: type Simplex_Rep ( n : natural ) is record
! 18: nor : vector(1..n); -- normal to the simplex
! 19: tra : Transfo; -- transformation to triangulate the cell
! 20: pts : points(1..n); -- vector of points
! 21: end record;
! 22:
! 23: -- AUXILIAIRIES :
! 24:
! 25: function Diagonalize ( s : simplex ) return matrix is
! 26:
! 27: -- DESCRIPTION :
! 28: -- Places the vertices of the simplex, shifted w.r.t. the first point,
! 29: -- in a matrix. Returns the triangulated matrix.
! 30:
! 31: a : matrix(1..s.n-1,1..s.n-1);
! 32: x,y : Vector(1..s.n-1);
! 33:
! 34: begin
! 35: for k in 2..s.n loop
! 36: x := s.pts(k).pt(x'range) - s.pts(1).pt(x'range);
! 37: y := s.tra*x;
! 38: for i in a'range(1) loop
! 39: a(i,k-1) := y(i);
! 40: end loop;
! 41: end loop;
! 42: return a;
! 43: end Diagonalize;
! 44:
! 45: function Diagonalize ( s : simplex; pt : Vector ) return matrix is
! 46:
! 47: -- DESCRIPTION :
! 48: -- Places the vertices of the simplex, shifted w.r.t. the first point,
! 49: -- in a matrix. Returns the triangulated matrix.
! 50:
! 51: a : matrix(1..s.n-1,1..s.n);
! 52: x,y : Vector(1..s.n-1);
! 53:
! 54: begin
! 55: for k in 2..s.n loop
! 56: x := s.pts(k).pt(x'range) - s.pts(1).pt(x'range); y := s.tra*x;
! 57: for i in a'range(1) loop
! 58: a(i,k-1) := y(i);
! 59: end loop;
! 60: end loop;
! 61: x := pt(x'range) - s.pts(1).pt(x'range); y := s.tra*x;
! 62: for i in a'range(1) loop
! 63: a(i,s.n) := y(i);
! 64: end loop;
! 65: return a;
! 66: end Diagonalize;
! 67:
! 68: function Create ( pts : VecVec ) return Transfo is
! 69:
! 70: a,l : matrix(pts'first..pts'last-1,pts'first..pts'last-1);
! 71: x : vector(pts(pts'first)'range);
! 72:
! 73: begin
! 74: for k in pts'first+1..pts'last loop
! 75: x := pts(k).all - pts(pts'first).all;
! 76: for i in a'range(1) loop
! 77: a(i,k-1) := x(i);
! 78: end loop;
! 79: end loop;
! 80: Upper_Triangulate(l,a);
! 81: return Create(l);
! 82: end Create;
! 83:
! 84: function Create ( pts : VecVec ) return Vector is
! 85:
! 86: a : matrix(pts'first..pts'last - 1,pts'range);
! 87: res : Vector(pts'range);
! 88:
! 89: begin
! 90: for k in a'range(1) loop
! 91: for i in a'range(2) loop
! 92: a(k,i) := pts(k+1)(i) - pts(pts'first)(i);
! 93: end loop;
! 94: end loop;
! 95: Upper_Triangulate(a);
! 96: Scale(a);
! 97: res := (res'range => 0);
! 98: Solve0(a,res);
! 99: Normalize(res);
! 100: if res(res'last) < 0
! 101: then return -res;
! 102: else return res;
! 103: end if;
! 104: end Create;
! 105:
! 106: -- CREATORS :
! 107:
! 108: function Create ( x : VecVec ) return Simplex is
! 109:
! 110: n : natural := x'last - x'first + 1;
! 111: res : Simplex := new Simplex_Rep(n);
! 112: cnt : natural := x'first;
! 113:
! 114: begin
! 115: for k in res.pts'range loop
! 116: res.pts(k).pt := x(cnt);
! 117: cnt := cnt + 1;
! 118: res.pts(k).si := Null_Simplex;
! 119: end loop;
! 120: res.tra := Create(x);
! 121: res.nor := Create(x);
! 122: return res;
! 123: end Create;
! 124:
! 125: procedure Update ( s : in out Simplex; x : in Link_to_Vector;
! 126: k : in natural ) is
! 127:
! 128: pts : VecVec(1..s.n);
! 129: nei : Simplex;
! 130:
! 131: begin
! 132: if s.pts(k).si = Null_Simplex
! 133: then for i in pts'range loop
! 134: if i = k
! 135: then pts(i) := x;
! 136: else pts(i) := s.pts(i).pt;
! 137: end if;
! 138: end loop;
! 139: nei := Create(pts);
! 140: s.pts(k).si := nei;
! 141: nei.pts(k).si := s;
! 142: end if;
! 143: end Update;
! 144:
! 145: procedure Update ( s : in out Simplex; x : in Link_to_Vector;
! 146: pos : in Vector ) is
! 147: begin
! 148: for k in pos'first..pos'last-1 loop
! 149: if pos(k)*pos(pos'last) > 0
! 150: then Update(s,x,k+1);
! 151: end if;
! 152: end loop;
! 153: end Update;
! 154:
! 155: procedure Update_One ( s : in out Simplex; x : in Link_to_Vector;
! 156: pos : in Vector ) is
! 157:
! 158: done : boolean := false;
! 159: nei : Simplex;
! 160:
! 161: begin
! 162: for k in pos'first..pos'last-1 loop
! 163: if pos(k)*pos(pos'last) > 0
! 164: then nei := s.pts(k+1).si;
! 165: if nei /= Null_Simplex
! 166: then Update_One(nei,x,Position(nei,x.all));
! 167: else Update(s,x,k+1);
! 168: end if;
! 169: done := true;
! 170: end if;
! 171: exit when done;
! 172: end loop;
! 173: end Update_One;
! 174:
! 175: procedure Update_One ( s : in out Simplex; x : in Link_to_Vector;
! 176: pos : in Vector; news : out Simplex ) is
! 177:
! 178: done : boolean := false;
! 179: nei,newsnei : Simplex;
! 180:
! 181: begin
! 182: -- LOOK FIRST FOR NULL SIMPLEX IN THE DIRECTION TO x :
! 183: for k in pos'first..pos'last-1 loop
! 184: if pos(k)*pos(pos'last) > 0
! 185: then nei := s.pts(k+1).si;
! 186: if nei = Null_Simplex
! 187: then Update(s,x,k+1);
! 188: news := s.pts(k+1).si;
! 189: done := true;
! 190: end if;
! 191: end if;
! 192: exit when done;
! 193: end loop;
! 194: -- WALK FURTHER IN THE DIRECTION TO x :
! 195: if not done
! 196: then Update_One(nei,x,Position(nei,x.all),newsnei);
! 197: if newsnei /= Null_Simplex
! 198: then news := newsnei;
! 199: s := nei;
! 200: end if;
! 201: end if;
! 202: end Update_One;
! 203:
! 204: procedure Update_All ( s : in out Simplex; x : in Link_to_Vector;
! 205: pos : in Vector; ancestor : in Simplex ) is
! 206:
! 207: nei : Simplex;
! 208: continue : boolean := true;
! 209:
! 210: begin
! 211: for k in pos'first..pos'last-1 loop
! 212: if pos(k)*pos(pos'last) > 0
! 213: then nei := s.pts(k+1).si;
! 214: if nei /= Null_Simplex
! 215: then if not Is_Vertex(nei,x.all) and nei /= ancestor
! 216: then Update_All(nei,x,Position(nei,x.all),s);
! 217: end if;
! 218: else Update(s,x,k+1);
! 219: Process(s.pts(k+1).si,continue);
! 220: end if;
! 221: end if;
! 222: exit when not continue;
! 223: end loop;
! 224: end Update_All;
! 225:
! 226: procedure Connect ( s1,s2 : in out Simplex ) is
! 227:
! 228: neighb : boolean;
! 229: index1,index2 : natural;
! 230:
! 231: begin
! 232: neighb := true; -- assume they are neighbors
! 233: index1 := 0;
! 234: -- SEARCH FOR INDEX OF POINT IN s1 THAT DOES NOT BELONG TO s2 :
! 235: for k in s1.pts'range loop
! 236: if not Is_Vertex(s2,s1.pts(k).pt.all)
! 237: then if (index1 = 0) and (s1.pts(k).si = Null_Simplex)
! 238: then index1 := k; -- kth point is not common
! 239: else neighb := false; -- more than one point not common
! 240: -- or there is already a neighbor
! 241: end if;
! 242: end if;
! 243: exit when not neighb;
! 244: end loop; -- either not neighb or index1 -> point in s1, not in s2
! 245: -- SEARCH FOR INDEX OF POINT IN s2 THAT DOES NOT BELONG TO s1 :
! 246: if neighb
! 247: then index2 := 0;
! 248: for k in s2.pts'range loop
! 249: if not Is_Vertex(s1,s2.pts(k).pt.all)
! 250: then if (index2 = 0) and (s2.pts(k).si = Null_Simplex)
! 251: then index2 := k; -- kth point is not common
! 252: else neighb := false; -- more than one point not common
! 253: -- or there is already a neighbor
! 254: end if;
! 255: end if;
! 256: exit when not neighb;
! 257: end loop; -- either no neighb or index2 -> point in s2, not in s1
! 258: -- CONNECT THE SIMPLICES WITH EACH OTHER :
! 259: if neighb
! 260: then s1.pts(index1).si := s2;
! 261: s2.pts(index2).si := s1;
! 262: end if;
! 263: end if;
! 264: end Connect;
! 265:
! 266: procedure Flatten ( s : in out Simplex ) is
! 267: begin
! 268: s.nor := (s.nor'range => 0);
! 269: s.nor(s.n) := 1;
! 270: for k in s.pts'range loop
! 271: s.pts(k).pt(s.n) := 0;
! 272: end loop;
! 273: end Flatten;
! 274:
! 275: -- SELECTORS :
! 276:
! 277: function Dimension ( s : Simplex ) return natural is
! 278: begin
! 279: return s.n;
! 280: end Dimension;
! 281:
! 282: function Normal ( s : Simplex ) return Vector is
! 283: begin
! 284: return s.nor;
! 285: end Normal;
! 286:
! 287: function Is_Flat ( s : Simplex ) return boolean is
! 288: begin
! 289: for i in s.nor'first..(s.nor'last-1) loop
! 290: if s.nor(i) /= 0
! 291: then return false;
! 292: end if;
! 293: end loop;
! 294: return (s.nor(s.nor'last) = 1);
! 295: end Is_Flat;
! 296:
! 297: function Vertices ( s : Simplex ) return VecVec is
! 298:
! 299: res : VecVec(s.pts'range);
! 300:
! 301: begin
! 302: for k in res'range loop
! 303: res(k) := s.pts(k).pt;
! 304: end loop;
! 305: return res;
! 306: end Vertices;
! 307:
! 308: function Vertex ( s : Simplex; k : natural ) return Vector is
! 309: begin
! 310: return s.pts(k).pt.all;
! 311: end Vertex;
! 312:
! 313: function Is_Vertex ( s : Simplex; x : Vector ) return boolean is
! 314: begin
! 315: for k in s.pts'range loop
! 316: if s.pts(k).pt.all = x
! 317: then return true;
! 318: end if;
! 319: end loop;
! 320: return false;
! 321: end Is_Vertex;
! 322:
! 323: function Equal ( s1,s2 : Simplex ) return boolean is
! 324:
! 325: found : boolean;
! 326:
! 327: begin
! 328: if s1.nor /= s2.nor -- check if normals are the same
! 329: then return false;
! 330: else for k in s1.pts'range loop -- check if vertices are the same
! 331: -- check whether s1.pts(k).pt.all occurs in s2.pts
! 332: found := false;
! 333: for l in s2.pts'range loop
! 334: if s1.pts(k).pt.all = s2.pts(l).pt.all
! 335: then found := true;
! 336: end if;
! 337: exit when found;
! 338: end loop;
! 339: if not found
! 340: then return false;
! 341: end if;
! 342: end loop;
! 343: return true;
! 344: end if;
! 345: end Equal;
! 346:
! 347: function Index ( s : Simplex; x : Vector ) return natural is
! 348: begin
! 349: for k in s.pts'range loop
! 350: if s.pts(k).pt.all = x
! 351: then return k;
! 352: end if;
! 353: end loop;
! 354: return 0;
! 355: end Index;
! 356:
! 357: function Neighbor ( s : Simplex; k : natural ) return Simplex is
! 358: begin
! 359: return s.pts(k).si;
! 360: end Neighbor;
! 361:
! 362: function Neighbor ( s : Simplex; k : natural; pos : Vector )
! 363: return Simplex is
! 364: begin
! 365: if pos(k-1)*pos(pos'last) > 0
! 366: then return s.pts(k).si;
! 367: else return Null_Simplex;
! 368: end if;
! 369: end Neighbor;
! 370:
! 371: function Position ( s : Simplex; x : Vector ) return Vector is
! 372:
! 373: m : matrix(x'first..x'last-1,x'range);
! 374: pos : Vector(x'range);
! 375: res : Vector(0..pos'last);
! 376:
! 377: begin
! 378: -- nbpos := nbpos + 1;
! 379: -- put("# position computations : "); put(nbpos,1); new_line;
! 380: -- transform point and simplex
! 381: m := Diagonalize(s,x);
! 382: -- solve the system
! 383: pos := (pos'range => 0);
! 384: Solve0(m,pos);
! 385: res(pos'first..pos'last) := pos;
! 386: res(0) := 0;
! 387: for k in pos'range loop
! 388: res(0) := res(0) + pos(k);
! 389: end loop;
! 390: res(0) := -res(0);
! 391: return res;
! 392: end Position;
! 393:
! 394: function Is_In ( s : Simplex; x : Vector ) return boolean is
! 395:
! 396: pos : Vector(0..x'last) := Position(s,x);
! 397:
! 398: begin
! 399: return Is_In(s,x,pos);
! 400: end Is_In;
! 401:
! 402: function Is_In ( s : Simplex; x,pos : Vector ) return boolean is
! 403: begin
! 404: for k in pos'first..pos'last-1 loop
! 405: if pos(k)*pos(pos'last) > 0
! 406: then return false;
! 407: end if;
! 408: end loop;
! 409: return true;
! 410: end Is_In;
! 411:
! 412: function Is_In_All ( s : Simplex; x : Vector ) return boolean is
! 413:
! 414: pos : Vector(0..x'last) := Position(s,x);
! 415:
! 416: begin
! 417: return Is_In_All(s,x,pos);
! 418: end Is_In_All;
! 419:
! 420: function Is_In_All ( s : Simplex; x : Vector ) return Simplex is
! 421:
! 422: pos : Vector(0..x'last) := Position(s,x);
! 423:
! 424: begin
! 425: return Is_In_All(s,x,pos);
! 426: end Is_In_All;
! 427:
! 428: function Is_In_All ( s : Simplex; x,pos : Vector ) return boolean is
! 429:
! 430: ins : boolean := true; -- assumes that x belongs to s
! 431:
! 432: begin
! 433: for k in pos'first..pos'last-1 loop
! 434: if pos(k)*pos(pos'last) > 0
! 435: then if s.pts(k+1).si /= Null_Simplex
! 436: then return Is_In_All(s.pts(k+1).si,x);
! 437: else ins := false;
! 438: end if;
! 439: end if;
! 440: end loop;
! 441: return ins;
! 442: end Is_In_All;
! 443:
! 444: function Is_In_All ( s : Simplex; x,pos : Vector ) return Simplex is
! 445:
! 446: ins : boolean := true; -- assumes that x belongs to s
! 447:
! 448: begin
! 449: for k in pos'first..pos'last-1 loop
! 450: if pos(k)*pos(pos'last) > 0
! 451: then if s.pts(k+1).si /= Null_Simplex
! 452: then return Is_In_All(s.pts(k+1).si,x);
! 453: else ins := false;
! 454: end if;
! 455: end if;
! 456: end loop;
! 457: if ins
! 458: then return s;
! 459: else return Null_Simplex;
! 460: end if;
! 461: end Is_In_All;
! 462:
! 463: procedure Neighbors ( s : in out Simplex; x : in Vector ) is
! 464:
! 465: cont : boolean;
! 466:
! 467: procedure Neighbors ( s : in out Simplex; x : in Vector;
! 468: cont : out boolean ) is
! 469:
! 470: pos : Vector(0..x'last) := Position(s,x);
! 471: continue : boolean := true;
! 472:
! 473: begin
! 474: for k in pos'first..pos'last-1 loop
! 475: if pos(k)*pos(pos'last) > 0
! 476: then if s.pts(k+1).si /= Null_Simplex
! 477: then Neighbors(s.pts(k+1).si,x,continue);
! 478: else Process_Neighbor(s,k+1,continue);
! 479: end if;
! 480: end if;
! 481: exit when not continue;
! 482: end loop;
! 483: cont := continue;
! 484: end Neighbors;
! 485:
! 486: begin
! 487: Neighbors(s,x,cont);
! 488: end Neighbors;
! 489:
! 490: function Volume ( s : Simplex ) return natural is
! 491:
! 492: m : constant matrix := Diagonalize(s);
! 493: vol : integer := 1;
! 494:
! 495: begin
! 496: for k in m'range(1) loop
! 497: vol := vol*m(k,k);
! 498: end loop;
! 499: if vol >= 0
! 500: then return vol;
! 501: else return -vol;
! 502: end if;
! 503: end Volume;
! 504:
! 505: -- DESTRUCTORS :
! 506:
! 507: procedure Destroy_Neighbor ( s : in out Simplex; k : in natural ) is
! 508: begin
! 509: s.pts(k).si := Null_Simplex;
! 510: end Destroy_Neighbor;
! 511:
! 512: procedure Destroy_Neighbors ( s : in out Simplex ) is
! 513: begin
! 514: for k in s.pts'range loop
! 515: Destroy_Neighbor(s,k);
! 516: end loop;
! 517: end Destroy_Neighbors;
! 518:
! 519: procedure Clear_Neighbor ( s : in out Simplex; k : in natural ) is
! 520: begin
! 521: Clear(s.pts(k).si);
! 522: end Clear_Neighbor;
! 523:
! 524: procedure Clear_Neighbors ( s : in out Simplex ) is
! 525: begin
! 526: for k in s.pts'range loop
! 527: Clear_Neighbor(s,k);
! 528: end loop;
! 529: end Clear_Neighbors;
! 530:
! 531: procedure Clear ( s : in out Simplex ) is
! 532:
! 533: procedure free is new unchecked_deallocation(Simplex_Rep,Simplex);
! 534:
! 535: begin
! 536: Clear(s.tra);
! 537: free(s);
! 538: end Clear;
! 539:
! 540: end Simplices;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to simplices.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