Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/vertices.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 2: with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
! 3: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
! 4: with Dictionaries;
! 5: with Linear_Programming; use Linear_Programming;
! 6: with Integer_Face_Enumerators; use Integer_Face_Enumerators;
! 7:
! 8: with Integer_Vectors_Utilities; use Integer_Vectors_Utilities;
! 9: with Lists_of_Vectors_Utilities; use Lists_of_Vectors_Utilities;
! 10: with Transformations; use Transformations;
! 11: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
! 12:
! 13: package body Vertices is
! 14:
! 15: function Is_In_Hull ( point : Standard_Integer_Vectors.Vector; l : List )
! 16: return boolean is
! 17:
! 18: fpt : Standard_Floating_Vectors.Vector(point'range);
! 19:
! 20: begin
! 21: for i in point'range loop
! 22: fpt(i) := double_float(point(i));
! 23: end loop;
! 24: return Is_In_Hull(fpt,l);
! 25: end Is_In_Hull;
! 26:
! 27: function Is_In_Hull ( point : Standard_Floating_Vectors.Vector; l : List )
! 28: return boolean is
! 29:
! 30: -- ALGORITHM:
! 31: -- The following linear program will be solved:
! 32: --
! 33: -- min u_0 + u_1 + .. + u_n
! 34: --
! 35: -- l_1*p_1i + l_2*p_2i + .. + l_m*p_mi + u_i*q_i = q_i i=1,2,..,n
! 36: --
! 37: -- l_1 + l_2 + .. + l_m + u_0 = 1
! 38: --
! 39: -- to determine whether q belongs to the convex hull spanned by the
! 40: -- vectors p_1,p_2,..,p_m.
! 41: -- If all u_i are zero and all constraints are satisfied,
! 42: -- then q belongs to the convex hull.
! 43:
! 44: -- CONSTANTS :
! 45:
! 46: n : constant natural := point'last;
! 47: m : constant natural := 2*n+2; -- number of constraints
! 48: nb : constant natural := Length_Of(l)+n+1; -- number of unknowns
! 49: eps : constant double_float := 10.0**(-10);
! 50:
! 51: -- VARIABLES :
! 52:
! 53: dic : matrix(0..m,0..nb);
! 54: sol : Standard_Floating_Vectors.vector(1..nb);
! 55: inbas : Standard_Integer_Vectors.Vector(1..m) := (1..m => 0);
! 56: outbas : Standard_Integer_Vectors.Vector(1..nb) := (1..nb => 0);
! 57: nit : natural := 0;
! 58: feasi : boolean;
! 59:
! 60: tmp : List;
! 61: pt : Link_to_Vector;
! 62: s : double_float;
! 63:
! 64: begin
! 65:
! 66: -- INITIALIZATION OF target :
! 67:
! 68: for i in 0..(nb-n-1) loop
! 69: dic(0,i) := 0.0; -- sum of the lambda's
! 70: end loop;
! 71: for i in (nb-n)..nb loop
! 72: dic(0,i) := -1.0; -- sum of the mu's
! 73: end loop;
! 74:
! 75: -- INITIALIZATION OF dic :
! 76:
! 77: for i in 0..(nb-n) loop
! 78: dic(m-1,i) := 1.0; -- sum of the lambda's + mu_0
! 79: dic(m,i) := -1.0;
! 80: end loop;
! 81: for i in (nb-n+1)..nb loop
! 82: dic(m-1,i) := 0.0;
! 83: dic(m,i) := 0.0;
! 84: end loop;
! 85: for i in 1..n loop
! 86: for j in (nb-n)..nb loop
! 87: if i /= j-nb+n
! 88: then dic(i,j) := 0.0;
! 89: dic(i+n,j) := 0.0;
! 90: end if;
! 91: end loop;
! 92: end loop;
! 93:
! 94: tmp := l;
! 95: for j in 1..(nb-n-1) loop
! 96: pt := Head_Of(tmp);
! 97: for i in pt'range loop
! 98: dic(i,j) := double_float(pt(i));
! 99: dic(i+n,j) := -double_float(pt(i));
! 100: end loop;
! 101: tmp := Tail_Of(tmp);
! 102: end loop;
! 103:
! 104: for i in point'range loop
! 105: dic(i,0) := point(i);
! 106: dic(i+n,0) := -point(i);
! 107: dic(i,i+nb-n) := point(i);
! 108: dic(i+n,i+nb-n) := -point(i);
! 109: end loop;
! 110:
! 111: -- SOLVE THE LINEAR PROGRAM :
! 112:
! 113: Dictionaries.Init_Basis(inbas,outbas);
! 114: Dual_Simplex(dic,eps,inbas,outbas,nit,feasi);
! 115: if not feasi
! 116: then return false;
! 117: else
! 118: sol := Dictionaries.Primal_Solution(dic,inbas,outbas);
! 119:
! 120: -- CHECK THE SOLUTION :
! 121:
! 122: s := 0.0;
! 123: for i in 1..(nb-n-1) loop
! 124: s := s + sol(i);
! 125: end loop;
! 126: if abs(s - 1.0) > eps
! 127: then return false;
! 128: end if;
! 129: s := 0.0;
! 130: for i in (nb-n)..nb loop
! 131: s := s + sol(i);
! 132: end loop;
! 133: if abs(s) > eps
! 134: then return false;
! 135: end if;
! 136:
! 137: return true;
! 138:
! 139: end if;
! 140:
! 141: end Is_In_Hull;
! 142:
! 143: function Vertex_Points ( l : List ) return List is
! 144:
! 145: result,result_last : List;
! 146: tmp,rest,origrest,rest_last : List;
! 147: pt : Link_to_Vector;
! 148:
! 149: begin
! 150: if Is_Null(l) or else Is_Null(Tail_Of(l))
! 151: then return l;
! 152: else Copy(Tail_Of(l),rest);
! 153: origrest := rest;
! 154: rest_last := rest;
! 155: while not Is_Null(Tail_Of(rest_last)) loop
! 156: rest_last := Tail_Of(rest_last);
! 157: end loop;
! 158: tmp := l;
! 159: for i in 1..Length_Of(l) loop
! 160: pt := Head_Of(tmp);
! 161: if not Is_In_Hull(pt.all,rest)
! 162: then Append(result,result_last,pt.all);
! 163: Append(rest,rest_last,pt.all);
! 164: end if;
! 165: rest := Tail_Of(rest);
! 166: tmp := Tail_Of(tmp);
! 167: end loop;
! 168: Clear(origrest);
! 169: return result;
! 170: end if;
! 171: end Vertex_Points;
! 172:
! 173: function Vertex_Points1 ( l : List ) return List is
! 174:
! 175: len : constant natural := Length_Of(l);
! 176: pts : VecVec(1..len) := Shallow_Create(l);
! 177: result,result_last : List;
! 178:
! 179: procedure Collect_Vertex ( i : in integer; cont : out boolean ) is
! 180: begin
! 181: Append(result,result_last,pts(i).all);
! 182: cont := true;
! 183: end Collect_Vertex;
! 184: procedure Enum_Vertices is new Enumerate_Vertices(Collect_Vertex);
! 185:
! 186: begin
! 187: Enum_Vertices(pts);
! 188: return result;
! 189: end Vertex_Points1;
! 190:
! 191: procedure Add ( pt : Link_to_Vector; l : in out List ) is
! 192:
! 193: -- DESCRIPTION :
! 194: -- Constructs the point to the list, but without sharing.
! 195:
! 196: newpt : Link_to_Vector := new Vector'(pt.all);
! 197:
! 198: begin
! 199: Construct(newpt,l);
! 200: end Add;
! 201:
! 202: function Extremal_Points ( l : List; v : Link_to_Vector ) return List is
! 203:
! 204: min,max,sp : integer;
! 205: tmp,res : List;
! 206: minpt,maxpt,pt : Link_to_Vector;
! 207:
! 208: begin
! 209: if Length_Of(l) <= 1
! 210: then Copy(l,res);
! 211: else pt := Head_Of(l);
! 212: min := pt*v; max := min;
! 213: minpt := pt; maxpt := pt;
! 214: tmp := Tail_Of(l);
! 215: while not Is_Null(tmp) loop
! 216: pt := Head_Of(tmp);
! 217: sp := pt*v;
! 218: if sp > max
! 219: then max := sp; maxpt := pt;
! 220: elsif sp < min
! 221: then min := sp; minpt := pt;
! 222: end if;
! 223: tmp := Tail_Of(tmp);
! 224: end loop;
! 225: Add(minpt,res);
! 226: if min /= max
! 227: then Add(maxpt,res);
! 228: end if;
! 229: end if;
! 230: return res;
! 231: end Extremal_Points;
! 232:
! 233: function Extremal_Points ( k,n : natural; exl,l : List ) return List is
! 234:
! 235: res,tmp,nres : List;
! 236: iv : VecVec(1..k);
! 237: v : Link_to_Vector := new Vector(1..n);
! 238: sp : integer;
! 239: pt : Link_to_Vector;
! 240: done : boolean;
! 241:
! 242: begin
! 243: tmp := exl;
! 244: for j in iv'range loop
! 245: iv(j) := Head_Of(tmp);
! 246: tmp := Tail_Of(tmp);
! 247: end loop;
! 248: v := Compute_Normal(iv);
! 249: sp := Head_Of(exl)*v;
! 250: nres := Extremal_Points(l,v);
! 251: tmp := nres;
! 252: res := exl;
! 253: done := false;
! 254: while not done and not Is_Null(tmp) loop
! 255: pt := Head_Of(tmp);
! 256: if pt*v /= sp
! 257: then Add(pt,res);
! 258: done := true;
! 259: else tmp := Tail_Of(tmp);
! 260: end if;
! 261: end loop;
! 262: Clear(v);
! 263: return res;
! 264: end Extremal_Points;
! 265:
! 266: function Max_Extremal_Points ( k,n : natural; exl,l : List ) return List is
! 267:
! 268: res,tmp,nres : List;
! 269: iv : VecVec(1..k);
! 270: v : Link_to_Vector := new Vector(1..n);
! 271: sp : integer;
! 272: pt : Link_to_Vector;
! 273: done : boolean;
! 274:
! 275: begin
! 276: tmp := exl;
! 277: for j in iv'range loop
! 278: iv(j) := Head_Of(tmp);
! 279: tmp := Tail_Of(tmp);
! 280: end loop;
! 281: v := Compute_Normal(iv);
! 282: sp := Head_Of(exl)*v;
! 283: nres := Extremal_Points(l,v);
! 284: --put("The computed normal : "); put(v); new_line;
! 285: --put("The inner product : "); put(sp,1); new_line;
! 286: --put_line("The list of extremal points :"); put(nres);
! 287: tmp := nres;
! 288: Copy(exl,res);
! 289: done := false;
! 290: while not done and not Is_Null(tmp) loop
! 291: pt := Head_Of(tmp);
! 292: if pt*v /= sp
! 293: then Add(pt,res);
! 294: done := true;
! 295: else tmp := Tail_Of(tmp);
! 296: end if;
! 297: end loop;
! 298: Clear(nres);
! 299: if not done
! 300: then -- for all points x in l : <x,v> = sp
! 301: if n >= 2
! 302: then declare
! 303: i : integer := Pivot(v);
! 304: t,invt : Transfo;
! 305: texl,tl,tres : List;
! 306: begin
! 307: if i <= v'last
! 308: then t := Build_Transfo(v,i);
! 309: invt := Invert(t);
! 310: texl := Transform_and_Reduce(t,i,exl);
! 311: tl := Transform_and_Reduce(t,i,l);
! 312: tres := Max_Extremal_Points(k,n-1,texl,tl);
! 313: --put("The normal : "); put(v); new_line;
! 314: --put_line("The list of points : "); put(l);
! 315: --put_line("The list of extremal points :"); put(exl);
! 316: --put_line("The transformed list of points : "); put(tl);
! 317: --put_line("The transformed list of extremal points : ");
! 318: --put(texl);
! 319: --put_line("The computed transformed extremal points : ");
! 320: --put(tres);
! 321: --put("k : "); put(k,1); new_line;
! 322: Clear(t); Clear(texl); Clear(tl);
! 323: if Length_Of(tres) = k+1
! 324: then res := Insert_and_Transform(tres,i,sp,invt);
! 325: else Copy(exl,res);
! 326: end if;
! 327: --put_line("The computed extremal points : "); put(res);
! 328: Clear(tres); --responsible for segmentation violation...
! 329: Clear(invt);
! 330: else Copy(exl,res);
! 331: end if;
! 332: end;
! 333: else Copy(exl,res);
! 334: end if;
! 335: end if;
! 336: Clear(v);
! 337: --put_line("The new list of extremal points : "); put(res);
! 338: return res;
! 339: end Max_Extremal_Points;
! 340:
! 341: function Extremal_Points ( n : natural; l : List ) return List is
! 342:
! 343: res : List;
! 344: nor : Link_to_Vector := new Vector'(1..n => 1);
! 345: k : natural;
! 346:
! 347: begin
! 348: res := Extremal_Points(l,nor);
! 349: Clear(nor);
! 350: if Length_Of(res) < 2
! 351: then return res;
! 352: else k := 2;
! 353: while k < n+1 loop
! 354: res := Extremal_Points(k,n,res,l);
! 355: exit when Length_Of(res) = k;
! 356: k := k+1;
! 357: end loop;
! 358: return res;
! 359: end if;
! 360: end Extremal_Points;
! 361:
! 362: function Two_Extremal_Points ( n : natural; l : List ) return List is
! 363:
! 364: -- DESCRIPTION :
! 365: -- Computes two extremal points of the list l.
! 366:
! 367: res,tmp : List;
! 368: pt,minpt,maxpt : Link_to_Vector;
! 369: min,max : integer;
! 370:
! 371: begin
! 372: if Length_Of(l) <= 2
! 373: then Copy(l,res);
! 374: else for i in 1..n loop
! 375: pt := Head_Of(l);
! 376: max := pt(i); min := max;
! 377: maxpt := pt; minpt := pt;
! 378: tmp := Tail_Of(l);
! 379: while not Is_Null(tmp) loop
! 380: pt := Head_Of(tmp);
! 381: if pt(i) < min
! 382: then min := pt(i); minpt := pt;
! 383: elsif pt(i) > max
! 384: then max := pt(i); maxpt := pt;
! 385: end if;
! 386: tmp := Tail_Of(tmp);
! 387: end loop;
! 388: if Is_Null(res)
! 389: then Add(minpt,res);
! 390: if min /= max
! 391: then Add(maxpt,res);
! 392: end if;
! 393: else pt := Head_Of(res);
! 394: if pt.all /= minpt.all
! 395: then Add(minpt,res);
! 396: elsif pt.all /= maxpt.all
! 397: then Add(maxpt,res);
! 398: end if;
! 399: end if;
! 400: exit when (Length_Of(res) = 2);
! 401: end loop;
! 402: end if;
! 403: return res;
! 404: end Two_Extremal_Points;
! 405:
! 406: function Max_Extremal_Points ( n : natural; l : List ) return List is
! 407:
! 408: res,wl,wres,newres : List;
! 409: k : natural;
! 410: nullvec,shiftvec : Link_to_Vector;
! 411: -- shifted : boolean;
! 412:
! 413: begin
! 414: if Length_Of(l) <= 2
! 415: then Copy(l,res);
! 416: else res := Two_Extremal_Points(n,l);
! 417: k := Length_Of(res);
! 418: if k = 2
! 419: then --nullvec := new Vector'(1..n => 0);
! 420: --if Is_In(nullvec,res)
! 421: -- then Copy(l,wl);
! 422: -- Copy(res,wres);
! 423: -- shifted := false;
! 424: -- else shiftvec := Head_Of(res);
! 425: -- wl := Shift(shiftvec,l);
! 426: -- wres := Shift(shiftvec,res);
! 427: -- shifted := true;
! 428: --end if;
! 429: --Clear(nullvec);
! 430: Copy(res,wres);
! 431: while k < n+1 loop
! 432: newres := Max_Extremal_Points(k,n,wres,l);
! 433: Copy(newres,wres);
! 434: exit when Length_Of(wres) = k;
! 435: k := k+1;
! 436: end loop;
! 437: res := wres;
! 438: -- Clear(res);
! 439: -- if shifted
! 440: -- then Min_Vector(shiftvec);
! 441: -- res := Shift(shiftvec,wres);
! 442: -- Clear(wres);
! 443: -- else Copy(wres,res);
! 444: -- end if;
! 445: -- Clear(wl);
! 446: end if;
! 447: end if;
! 448: return res;
! 449: end Max_Extremal_Points;
! 450:
! 451: end Vertices;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to vertices.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 / Implift