Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/inner_normal_cones.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
! 2: with Transformations; use Transformations;
! 3: with Integer_Vectors_Utilities; use Integer_Vectors_Utilities;
! 4: with Lists_of_Vectors_Utilities; use Lists_of_Vectors_Utilities;
! 5: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
! 6: with Vertices; use Vertices;
! 7:
! 8: package body Inner_Normal_Cones is
! 9:
! 10: -- NOTE ON THE IMPLEMENTATION :
! 11: -- There is a safety mode in which after each computation of the generators,
! 12: -- it is automatically tested whether the generators satisfy the
! 13: -- inequalities of the normal cone. If not, then the bug is reported and
! 14: -- a program_error is raised.
! 15:
! 16: -- AUXILIAIRIES :
! 17:
! 18: function Lower ( a1,b1,a2,b2 : integer ) return boolean is
! 19:
! 20: -- DESCRIPTION :
! 21: -- Returns true if a1/b1 < a2/b2, false otherwise.
! 22:
! 23: -- REQUIRED : b1 /= 0 and b2 /= 0.
! 24:
! 25: begin
! 26: if b1*b2 > 0
! 27: then return (a1*b2 - a2*b1 < 0);
! 28: else return (a1*b2 - a2*b1 > 0);
! 29: end if;
! 30: end Lower;
! 31:
! 32: function Higher ( a1,b1,a2,b2 : integer ) return boolean is
! 33:
! 34: -- DESCRIPTION :
! 35: -- Returns true if a1/b1 > a2/b2, false otherwise.
! 36:
! 37: -- REQUIRED : b1 /= 0 and b2 /= 0.
! 38:
! 39: begin
! 40: if b1*b2 > 0
! 41: then return (a1*b2 - a2*b1 > 0);
! 42: else return (a1*b2 - a2*b1 < 0);
! 43: end if;
! 44: end Higher;
! 45:
! 46: function Compute_Facets ( l : List; x : Vector ) return Faces is
! 47:
! 48: -- DESCRIPTION :
! 49: -- Returns all facets of conv(l) that all contain x.
! 50: -- Checks whether x belongs to the list or not.
! 51:
! 52: res : Faces;
! 53: n : constant natural := x'length;
! 54:
! 55: begin
! 56: if Is_In(l,x)
! 57: then res := Create(n-1,n,l,x);
! 58: else declare
! 59: wrk : List := l;
! 60: lx : Link_to_Vector;
! 61: begin
! 62: lx := new Vector'(x);
! 63: Construct(lx,wrk);
! 64: res := Create(n-1,n,wrk,x);
! 65: end;
! 66: end if;
! 67: return res;
! 68: end Compute_Facets;
! 69:
! 70: procedure Shifted_Points ( l : in List; x : in Vector;
! 71: res,res_last : in out List ) is
! 72:
! 73: -- DESCRIPTION :
! 74: -- Appends the list of shifted points w.r.t. x to the list res.
! 75:
! 76: tmp : List := l;
! 77:
! 78: begin
! 79: tmp := l;
! 80: while not Is_Null(tmp) loop
! 81: Append_Diff(res,res_last,Head_Of(tmp).all-x);
! 82: tmp := Tail_Of(tmp);
! 83: end loop;
! 84: end Shifted_Points;
! 85:
! 86: function In_Hull ( l : List; x : Vector ) return boolean is
! 87:
! 88: -- DESCRIPTION :
! 89: -- Returns true if the vector x belongs the convex hull spanned by
! 90: -- the points in l minus the point x itself.
! 91:
! 92: res : boolean;
! 93: lx : List;
! 94:
! 95: begin
! 96: Copy(l,lx);
! 97: Remove(lx,x);
! 98: res := Is_In_Hull(x,lx);
! 99: Clear(lx);
! 100: return res;
! 101: end In_Hull;
! 102:
! 103: -- AUXILIAIRIES TO CONSTRUCT THE NORMALS :
! 104:
! 105: procedure Normal ( m : in out Matrix; v : out Vector ) is
! 106:
! 107: -- DESCRIPTION :
! 108: -- Computes the normal to the vectors stored in the rows of the matrix m.
! 109:
! 110: -- REQUIRED : the matrix m is square!
! 111:
! 112: res : Vector(m'range(2));
! 113:
! 114: begin
! 115: Upper_Triangulate(m); Scale(m);
! 116: res := (res'range => 0);
! 117: Solve0(m,res);
! 118: v := res;
! 119: end Normal;
! 120:
! 121: procedure Orientate_Normal ( l : in List; f : in Face;
! 122: normal : in out Vector ) is
! 123:
! 124: -- DESCRIPTION :
! 125: -- Orientates the normal of the face such that it becomes an inner normal
! 126: -- w.r.t. the points in the list l.
! 127:
! 128: tmp : List := l;
! 129: ip : constant integer := f(f'first).all*normal;
! 130: pt : Vector(Head_Of(l)'range);
! 131: done : boolean := false;
! 132: ptip : integer;
! 133:
! 134: begin
! 135: while not Is_Null(tmp) loop
! 136: pt := Head_Of(tmp).all;
! 137: if not Is_In(f,pt)
! 138: then ptip := pt*normal;
! 139: if ptip /= ip
! 140: then if ptip < ip
! 141: then normal := -normal;
! 142: end if;
! 143: done := true;
! 144: end if;
! 145: end if;
! 146: exit when done;
! 147: tmp := Tail_Of(tmp);
! 148: end loop;
! 149: end Orientate_Normal;
! 150:
! 151: function Inner_Normal ( l : List; facet : Face ) return Vector is
! 152:
! 153: -- DESCRIPTION :
! 154: -- Returns the inner normal to the facet.
! 155:
! 156: res,fst : Vector(facet(facet'first)'range);
! 157: m : Matrix(facet'first+1..facet'last,facet(facet'first)'range);
! 158:
! 159: begin
! 160: fst := facet(facet'first).all;
! 161: for i in m'first(1)..m'last(1) loop
! 162: for j in m'range(2) loop
! 163: m(i,j) := facet(i)(j) - fst(j);
! 164: end loop;
! 165: end loop;
! 166: Normal(m,res);
! 167: Orientate_Normal(l,facet,res);
! 168: return res;
! 169: end Inner_Normal;
! 170:
! 171: procedure Inner_Normals ( l : in List; facets : in Faces;
! 172: iv,iv_last : in out List ) is
! 173:
! 174: -- DESCRIPTION :
! 175: -- Returns the list of all inner normals to the facets of conv(l).
! 176: -- If there is only one facet, then both `inner' and `outer' normal
! 177: -- will be returned, because there is nothing to orientate with.
! 178: -- This means that also the negation of the normal will be in the list iv.
! 179:
! 180: tmp : Faces := facets;
! 181:
! 182: begin
! 183: while not Is_Null(tmp) loop
! 184: declare
! 185: f : Face := Head_Of(tmp);
! 186: v : constant Vector := Inner_Normal(l,Head_Of(tmp));
! 187: begin
! 188: Append(iv,iv_last,v);
! 189: if Length_Of(facets) = 1
! 190: then Append(iv,iv_last,-v);
! 191: end if;
! 192: end;
! 193: tmp := Tail_Of(tmp);
! 194: end loop;
! 195: end Inner_Normals;
! 196:
! 197: function Number_of_Rows ( l : List ) return natural is
! 198:
! 199: -- DESCRIPTION :
! 200: -- Returns Maximum(dimension,Length_Of(l)).
! 201:
! 202: -- REQUIRED : not Is_Null(l).
! 203:
! 204: n : constant natural := Head_Of(l)'length;
! 205: m : constant natural := Length_Of(l);
! 206:
! 207: begin
! 208: if n > m
! 209: then return n;
! 210: else return m;
! 211: end if;
! 212: end Number_of_Rows;
! 213:
! 214: function Normal ( l : List ) return Vector is
! 215:
! 216: -- DESCRIPTION :
! 217: -- Return a normal to the vectors in the list.
! 218:
! 219: -- REQUIRED : Length_Of(l) <= n = Head_Of(l)'length
! 220: -- or the points in l all lie on the same facet.
! 221:
! 222: fst : Link_to_Vector := Head_Of(l);
! 223: m : Matrix(1..Number_of_Rows(l),fst'range);
! 224: cnt : natural := 0;
! 225: res : Vector(fst'range);
! 226: tmp : List := Tail_Of(l);
! 227:
! 228: begin
! 229: while not Is_Null(tmp) loop
! 230: res := Head_Of(tmp).all;
! 231: cnt := cnt+1;
! 232: for i in res'range loop
! 233: m(cnt,i) := res(i) - fst(i);
! 234: end loop;
! 235: tmp := Tail_Of(tmp);
! 236: end loop;
! 237: for i in cnt+1..m'last(1) loop
! 238: for j in m'range(2) loop
! 239: m(i,j) := 0;
! 240: end loop;
! 241: end loop;
! 242: Normal(m,res);
! 243: return res;
! 244: end Normal;
! 245:
! 246: procedure Normals ( l : in List; x : in Vector; iv,iv_last : in out List ) is
! 247:
! 248: -- DESCRIPTION :
! 249: -- Computes a list of all normals to the vectors in the list.
! 250:
! 251: -- REQUIRED : Length_Of(l) <= n = x'length.
! 252: -- or all points in l lie on the same facet.
! 253:
! 254: nor : constant Vector := Normal(l);
! 255: piv : natural := Pivot(nor);
! 256: pivnor : Vector(nor'range); -- normal with pivnor(piv) > 0
! 257:
! 258: begin
! 259: if piv <= nor'last
! 260: then if x'length > 2
! 261: then if nor(piv) < 0
! 262: then pivnor := -nor;
! 263: else pivnor := nor;
! 264: end if;
! 265: declare
! 266: t : Transfo := Build_Transfo(pivnor,piv);
! 267: tx : constant Vector := Reduce(t*x,piv);
! 268: tt : Transfo := Transpose(t);
! 269: tl : List := Transform_and_Reduce(t,piv,l);
! 270: begin
! 271: if Length_Of(tl) <= nor'length-1
! 272: then Normals(tl,tx,iv,iv_last);
! 273: else declare
! 274: facets : Faces := Compute_Facets(tl,tx);
! 275: begin
! 276: if Length_Of(facets) > 1
! 277: then Inner_Normals(tl,facets,iv,iv_last);
! 278: else Normals(tl,tx,iv,iv_last);
! 279: end if;
! 280: end;
! 281: end if;
! 282: Insert_and_Transform(iv,piv,0,tt);
! 283: iv_last := Pointer_to_Last(iv);
! 284: Clear(t); Deep_Clear(tl); Clear(tt);
! 285: end;
! 286: end if;
! 287: Append(iv,iv_last,nor); Append(iv,iv_last,-nor);
! 288: end if;
! 289: end Normals;
! 290:
! 291: -- CONSTRUCTORS FOR PRIMAL REPRESENTATION :
! 292:
! 293: function Generators ( l : List; facets : Faces; x : Vector ) return List is
! 294:
! 295: res,res_last : List;
! 296:
! 297: begin
! 298: if Length_Of(facets) > 1
! 299: then Inner_Normals(l,facets,res,res_last); -- compute inner normals
! 300: elsif In_Hull(l,x)
! 301: then return res; -- empty normal cone
! 302: else Normals(l,x,res,res_last); -- compute normals
! 303: if Length_Of(l) = 2
! 304: then Shifted_Points(l,x,res,res_last);
! 305: end if;
! 306: end if;
! 307: -- SAFETY MODE :
! 308: -- if not Contained_in_Cone(l,x,res)
! 309: -- then put_line("Bug raised for computing generators for list"); put(l);
! 310: -- put(" and the point : "); put(x); new_line;
! 311: -- raise PROGRAM_ERROR;
! 312: -- end if;
! 313: return res;
! 314: end Generators;
! 315:
! 316: function Generators ( l : List; x : Vector ) return List is
! 317:
! 318: res,res_last : List;
! 319: n : constant natural := x'length;
! 320: facets : Faces;
! 321:
! 322: begin
! 323: if In_Hull(l,x)
! 324: then return res; -- empty normal cone
! 325: else if Length_Of(l) > n
! 326: then facets := Compute_Facets(l,x);
! 327: end if;
! 328: if Length_Of(facets) > 1
! 329: then res := Generators(l,facets,x);
! 330: else Normals(l,x,res,res_last);
! 331: if Length_Of(l) = 2
! 332: then Shifted_Points(l,x,res,res_last);
! 333: end if;
! 334: end if;
! 335: -- SAFETY MODE :
! 336: -- if not Contained_in_Cone(l,x,res)
! 337: -- then put_line("Bug raised for computing generators for list");
! 338: -- put(l); put(" and the point : "); put(x); new_line;
! 339: -- raise PROGRAM_ERROR;
! 340: -- end if;
! 341: return res;
! 342: end if;
! 343: end Generators;
! 344:
! 345: -- CONSTRUCTOR FOR DUAL REPRESENTATION :
! 346:
! 347: function Inner_Normal_Cone ( l : List; x : Vector ) return Matrix is
! 348:
! 349: res : Matrix(x'range,1..Length_Of(l)-1);
! 350: -- CAUTION : Is_In(l,x) must hold !!
! 351: y : Vector(x'range);
! 352: cnt : natural := 0;
! 353: tmp : List := l;
! 354:
! 355: begin
! 356: while not Is_Null(tmp) loop
! 357: y := Head_Of(tmp).all;
! 358: if not Equal(x,y) -- avoid trivial inequality
! 359: then cnt := cnt+1;
! 360: for i in x'range loop
! 361: res(i,cnt) := y(i) - x(i);
! 362: end loop;
! 363: end if;
! 364: tmp := Tail_Of(tmp);
! 365: end loop;
! 366: return res;
! 367: end Inner_Normal_Cone;
! 368:
! 369: function Included_Normal_Cone ( gx : List; x : Vector ) return Matrix is
! 370:
! 371: res : Matrix(x'first-1..x'last,1..Length_Of(gx));
! 372: tmp : List := gx;
! 373: v : Vector(x'range);
! 374:
! 375: begin
! 376: for j in res'range(2) loop
! 377: v := Head_Of(tmp).all;
! 378: res(res'first(1),j) := x*v;
! 379: for i in res'first(1)+1..res'last(1) loop
! 380: res(i,j) := v(i);
! 381: end loop;
! 382: tmp := Tail_Of(tmp);
! 383: end loop;
! 384: return res;
! 385: end Included_Normal_Cone;
! 386:
! 387: -- PRIMITIVE SELECTORS :
! 388:
! 389: function Evaluate ( m : Matrix; i : integer; v : Vector ) return integer is
! 390:
! 391: sum : integer := 0;
! 392:
! 393: begin
! 394: for j in v'range loop
! 395: sum := sum + m(j,i)*v(j);
! 396: end loop;
! 397: return sum;
! 398: end Evaluate;
! 399:
! 400: function Satisfies ( m : Matrix; i : integer; v : Vector ) return boolean is
! 401: begin
! 402: if m'first(1) = v'first-1
! 403: then return (Evaluate(m,i,v) >= m(0,i));
! 404: else return (Evaluate(m,i,v) >= 0);
! 405: end if;
! 406: end Satisfies;
! 407:
! 408: function Strictly_Satisfies ( m : Matrix; i : integer; v : Vector )
! 409: return boolean is
! 410: begin
! 411: if m'first(1) = v'first-1
! 412: then return (Evaluate(m,i,v) > m(0,i));
! 413: else return (Evaluate(m,i,v) > 0);
! 414: end if;
! 415: end Strictly_Satisfies;
! 416:
! 417: function Satisfies ( m : Matrix; v : Vector ) return boolean is
! 418: begin
! 419: for i in m'range(2) loop
! 420: if not Satisfies(m,i,v)
! 421: then return false;
! 422: end if;
! 423: end loop;
! 424: return true;
! 425: end Satisfies;
! 426:
! 427: function Strictly_Satisfies ( m : Matrix; v : Vector ) return boolean is
! 428: begin
! 429: for i in m'range(2) loop
! 430: if not Strictly_Satisfies(m,i,v)
! 431: then return false;
! 432: end if;
! 433: end loop;
! 434: return true;
! 435: end Strictly_Satisfies;
! 436:
! 437: function Satisfies ( m : Matrix; v : List ) return boolean is
! 438:
! 439: tmp : List := v;
! 440:
! 441: begin
! 442: while not Is_Null(tmp) loop
! 443: if not Satisfies(m,Head_Of(tmp).all)
! 444: then return false;
! 445: else tmp := Tail_Of(tmp);
! 446: end if;
! 447: end loop;
! 448: return true;
! 449: end Satisfies;
! 450:
! 451: function Strictly_Satisfies ( m : Matrix; v : List ) return boolean is
! 452:
! 453: tmp : List := v;
! 454:
! 455: begin
! 456: while not Is_Null(tmp) loop
! 457: if not Strictly_Satisfies(m,Head_Of(tmp).all)
! 458: then return false;
! 459: else tmp := Tail_Of(tmp);
! 460: end if;
! 461: end loop;
! 462: return true;
! 463: end Strictly_Satisfies;
! 464:
! 465: -- SECONDARY SELECTORS :
! 466:
! 467: function Satisfy ( m : Matrix; l : List ) return List is
! 468:
! 469: res,res_last,tmp : List;
! 470:
! 471: begin
! 472: tmp := l;
! 473: while not Is_Null(tmp) loop
! 474: if Satisfies(m,Head_Of(tmp).all)
! 475: then Append(res,res_last,Head_Of(tmp).all);
! 476: end if;
! 477: tmp := Tail_Of(tmp);
! 478: end loop;
! 479: return res;
! 480: end Satisfy;
! 481:
! 482: function Strictly_Satisfy ( m : Matrix; l : List ) return List is
! 483:
! 484: res,res_last,tmp : List;
! 485:
! 486: begin
! 487: tmp := l;
! 488: while not Is_Null(tmp) loop
! 489: if Strictly_Satisfies(m,Head_Of(tmp).all)
! 490: then Append(res,res_last,Head_Of(tmp).all);
! 491: end if;
! 492: tmp := Tail_Of(tmp);
! 493: end loop;
! 494: return res;
! 495: end Strictly_Satisfy;
! 496:
! 497: function Contained_in_Cone ( l : List; x : Vector; v : List )
! 498: return boolean is
! 499:
! 500: ineq : Matrix(x'range,1..Length_Of(l)-1) := Inner_Normal_Cone(l,x);
! 501:
! 502: begin
! 503: return Satisfies(ineq,v);
! 504: end Contained_in_Cone;
! 505:
! 506: function Strictly_Contained_in_Cone ( l : List; x : Vector; v : List )
! 507: return boolean is
! 508:
! 509: ineq : Matrix(x'range,1..Length_Of(l)-1) := Inner_Normal_Cone(l,x);
! 510:
! 511: begin
! 512: return Strictly_Satisfies(ineq,v);
! 513: end Strictly_Contained_in_Cone;
! 514:
! 515: function Contained_in_Cone ( l,v : List ) return boolean is
! 516:
! 517: tmp : List := l;
! 518:
! 519: begin
! 520: while not Is_Null(tmp) loop
! 521: if Contained_in_Cone(l,Head_Of(tmp).all,v)
! 522: then -- put("true for vector : "); put(Head_Of(tmp).all); new_line;
! 523: return true;
! 524: else tmp := Tail_Of(tmp);
! 525: end if;
! 526: end loop;
! 527: return false;
! 528: end Contained_in_Cone;
! 529:
! 530: function Strictly_Contained_in_Cone ( l,v : List ) return boolean is
! 531:
! 532: tmp : List := l;
! 533:
! 534: begin
! 535: while not Is_Null(tmp) loop
! 536: if Strictly_Contained_in_Cone(l,Head_Of(tmp).all,v)
! 537: then return true;
! 538: else tmp := Tail_Of(tmp);
! 539: end if;
! 540: end loop;
! 541: return false;
! 542: end Strictly_Contained_in_Cone;
! 543:
! 544: -- CONCERNING THE UNION OF CONES :
! 545:
! 546: function In_Union ( v1,v2 : Vector; k1,k2 : Matrix ) return boolean is
! 547:
! 548: -- ALGORITHM : consider x(t) = (1-t)*v1 + t*v2, for t in [0,1]
! 549: -- and compute t1: smallest t such that Satisfies(k1,x(t))
! 550: -- and t2: largest t such that Satisfies(k2,x(t)).
! 551: -- Then t1 >= t2 makes the test returning true, false otherwise.
! 552:
! 553: diff : constant Vector := v1-v2;
! 554: t1a,t1b,t2a,t2b,a,b : integer; -- t1 = t1a/t1b and t2 = t2a/t2b
! 555:
! 556: begin
! 557: t1a := 1; t1b := 1; -- initially t1 = 1 = 1/1
! 558: for i in k1'range(2) loop
! 559: b := Evaluate(k1,i,diff);
! 560: if b /= 0
! 561: then a := Evaluate(k1,i,v1);
! 562: if Lower(a,b,t1a,t1b)
! 563: then t1a := a; t1b := b;
! 564: end if;
! 565: end if;
! 566: end loop;
! 567: t2a := 0; t2b := 1; -- initially t2 = 0 = 0/1
! 568: for i in k2'range(2) loop
! 569: b := Evaluate(k2,i,diff);
! 570: if b /= 0
! 571: then a := Evaluate(k2,i,v1);
! 572: if Higher(a,b,t2a,t2b)
! 573: then t2a := a; t2b := b;
! 574: end if;
! 575: end if;
! 576: end loop;
! 577: -- put("In_Union for "); put(v1); put(" and "); put(v2);
! 578: -- put(" : t1 = "); put(t1a,1); put("/"); put(t1b,1);
! 579: -- put(" : t2 = "); put(t2a,1); put("/"); put(t2b,1); new_line;
! 580: -- put_line("k1 : "); put(k1); put_line("k2 : "); put(k2);
! 581: return not Lower(t1a,t1b,t2a,t2b);
! 582: end In_Union;
! 583:
! 584: function In_Union ( v1,v2 : List; k1,k2 : Matrix ) return boolean is
! 585:
! 586: tmpv1,tmpv2 : List;
! 587:
! 588: begin
! 589: tmpv1 := v1;
! 590: while not Is_Null(tmpv1) loop
! 591: tmpv2 := v2;
! 592: while not Is_Null(tmpv2) loop
! 593: if not In_Union(Head_Of(tmpv1).all,Head_Of(tmpv2).all,k1,k2)
! 594: then return false;
! 595: else tmpv2 := Tail_Of(tmpv2);
! 596: end if;
! 597: end loop;
! 598: tmpv1 := Tail_Of(tmpv1);
! 599: end loop;
! 600: return true;
! 601: end In_Union;
! 602:
! 603: end Inner_Normal_Cones;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to inner_normal_cones.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