Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/mixed_volume_computation.adb, Revision 1.1
1.1 ! maekawa 1: with integer_io; use integer_io;
! 2: with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
! 3: with Arrays_of_Integer_Vector_Lists_io; use Arrays_of_Integer_Vector_Lists_io;
! 4:
! 5: with Standard_Floating_Vectors;
! 6: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
! 7: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
! 8: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
! 9: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
! 10: with Integer_Lifting_Utilities; use Integer_Lifting_Utilities;
! 11: with Integer_Mixed_Subdivisions_io; use Integer_Mixed_Subdivisions_io;
! 12: with Mixed_Coherent_Subdivisions; use Mixed_Coherent_Subdivisions;
! 13:
! 14: package body Mixed_Volume_Computation is
! 15:
! 16: -- AUXILIAIRY OUTPUT ROUTINES :
! 17:
! 18: procedure put ( file : in file_type; points : in Array_of_Lists;
! 19: n : in natural; mix : in Vector;
! 20: mixsub : in out Mixed_Subdivision; mv : out natural ) is
! 21: begin
! 22: new_line(file);
! 23: put_line(file,"THE LIFTED SUPPORTS :");
! 24: new_line(file);
! 25: put(file,points);
! 26: new_line(file);
! 27: put_line(file,"THE MIXED SUBDIVISION :");
! 28: new_line(file);
! 29: put(file,n,mix,mixsub,mv);
! 30: end put;
! 31:
! 32: procedure Sort ( supports : in out Array_of_Lists; k,nb,n : in natural;
! 33: mxt,perm : in out Vector ) is
! 34:
! 35: -- DESCRIPTION :
! 36: -- Auxiliary operation for Compute_Mixture.
! 37: -- Compares the kth support with the following supports.
! 38: -- Already nb different supports have been found.
! 39:
! 40: begin
! 41: for l in (k+1)..n loop
! 42: if Equal(Supports(k),Supports(l))
! 43: then if l /= k + mxt(nb)
! 44: then declare
! 45: pos : natural := k + mxt(nb);
! 46: tmpdl : List := supports(l);
! 47: tmppos : natural;
! 48: begin
! 49: supports(l) := supports(pos);
! 50: supports(pos) := tmpdl;
! 51: tmppos := perm(l);
! 52: perm(l) := perm(pos);
! 53: perm(pos) := tmppos;
! 54: end;
! 55: end if;
! 56: mxt(nb) := mxt(nb) + 1;
! 57: end if;
! 58: end loop;
! 59: end Sort;
! 60:
! 61: -- TARGET ROUTINES :
! 62:
! 63: procedure Compute_Mixture ( supports : in out Array_of_Lists;
! 64: mix,perms : out Link_to_Vector ) is
! 65:
! 66: n : constant natural := supports'last;
! 67: cnt : natural := 0; -- counts the number of different supports
! 68: mxt : Vector(supports'range) -- counts the number of occurrencies
! 69: := (supports'range => 1);
! 70: perm : Link_to_Vector -- keeps track of the permutations
! 71: := new Standard_Integer_Vectors.Vector(supports'range);
! 72: index : natural := supports'first;
! 73:
! 74: begin
! 75: for k in perm'range loop
! 76: perm(k) := k;
! 77: end loop;
! 78: while index <= supports'last loop
! 79: cnt := cnt + 1;
! 80: Sort(supports,index,cnt,n,mxt,perm.all);
! 81: index := index + mxt(cnt);
! 82: end loop;
! 83: mix := new Standard_Integer_Vectors.Vector'(mxt(mxt'first..cnt));
! 84: perms := perm;
! 85: end Compute_Mixture;
! 86:
! 87: function Compute_Index ( k : natural; mix : Vector ) return natural is
! 88:
! 89: -- DESCRIPTION :
! 90: -- Returns the index of k w.r.t. to the type of mixture.
! 91:
! 92: index : natural := mix(mix'first);
! 93:
! 94: begin
! 95: if k <= index
! 96: then return mix'first;
! 97: else for l in (mix'first+1)..mix'last loop
! 98: index := index + mix(l);
! 99: if k <= index
! 100: then return l;
! 101: end if;
! 102: end loop;
! 103: return mix'last;
! 104: end if;
! 105: end Compute_Index;
! 106:
! 107: function Compute_Permutation ( n : natural; mix : Vector;
! 108: supports : Array_of_Lists )
! 109: return Link_to_Vector is
! 110:
! 111: perms : Link_to_Vector := new Vector(1..n);
! 112:
! 113: begin
! 114: for k in perms'range loop
! 115: perms(k) := k;
! 116: end loop;
! 117: return perms;
! 118: end Compute_Permutation;
! 119:
! 120: function Permute ( p : Poly_Sys; perm : Link_to_Vector ) return Poly_Sys is
! 121:
! 122: res : Poly_Sys(p'range);
! 123:
! 124: begin
! 125: for k in p'range loop
! 126: res(k) := p(perm(k));
! 127: end loop;
! 128: return res;
! 129: end Permute;
! 130:
! 131: function Permute ( supports : Array_of_Lists; perm : Link_to_Vector )
! 132: return Array_of_Lists is
! 133:
! 134: res : Array_of_Lists(supports'range);
! 135:
! 136: begin
! 137: for k in supports'range loop
! 138: res(k) := supports(perm(k));
! 139: end loop;
! 140: return res;
! 141: end Permute;
! 142:
! 143: function Typed_Lists ( mix : Vector; points : Array_of_Lists )
! 144: return Array_of_Lists is
! 145:
! 146: res : Array_of_Lists(mix'range);
! 147: ind : natural := res'first;
! 148:
! 149: begin
! 150: for i in mix'range loop
! 151: res(i) := points(ind);
! 152: ind := ind + mix(i);
! 153: end loop;
! 154: return res;
! 155: end Typed_Lists;
! 156:
! 157: -- MIXED VOLUME COMPUTATIONS BASED ON SUBDIVISIONS :
! 158:
! 159: -- AUXILIARIES :
! 160:
! 161: function Is_Fine ( mix : Vector; mic : Mixed_Cell ) return boolean is
! 162:
! 163: -- DESCRIPTION :
! 164: -- Returns true if the mixed volume can be computed by a determinant.
! 165:
! 166: fine : boolean := true;
! 167:
! 168: begin
! 169: for k in mic.pts'range loop
! 170: fine := (Length_Of(mic.pts(k)) = mix(k) + 1);
! 171: exit when not fine;
! 172: end loop;
! 173: return fine;
! 174: end Is_Fine;
! 175:
! 176: function Reduced_Supports ( n : natural; mix : Vector; mic : Mixed_Cell )
! 177: return Array_of_Lists is
! 178:
! 179: -- DESCRIPTION :
! 180: -- Returns the supports of the cell without the lifting values.
! 181:
! 182: res : Array_of_Lists(1..n);
! 183: cnt : natural := 1;
! 184:
! 185: begin
! 186: for k in mic.pts'range loop
! 187: res(cnt) := Reduce(mic.pts(k),n+1);
! 188: for l in 1..mix(k)-1 loop
! 189: Copy(res(cnt),res(cnt+l));
! 190: end loop;
! 191: cnt := cnt + mix(k);
! 192: end loop;
! 193: return res;
! 194: end Reduced_Supports;
! 195:
! 196: function Fine_Mixed_Volume ( n : natural; mix : Vector; mic : Mixed_Cell )
! 197: return natural is
! 198:
! 199: -- DESCRIPTION :
! 200: -- Computes the mixed volume for a cell that is fine mixed.
! 201:
! 202: -- REQUIRED : Fine(mix,mic).
! 203:
! 204: res,count : natural;
! 205: mat : matrix(1..n,1..n);
! 206: detmat : integer;
! 207: tmp : List;
! 208: sh,pt : Link_to_Vector;
! 209:
! 210: begin
! 211: count := 1;
! 212: for k in mic.pts'range loop
! 213: sh := Head_Of(mic.pts(k));
! 214: tmp := Tail_Of(mic.pts(k));
! 215: while not Is_Null(tmp) loop
! 216: pt := Head_Of(tmp);
! 217: for j in 1..n loop
! 218: mat(count,j) := pt(j) - sh(j);
! 219: end loop;
! 220: tmp := Tail_Of(tmp);
! 221: count := count + 1;
! 222: end loop;
! 223: end loop;
! 224: detmat := Det(mat);
! 225: if detmat >= 0
! 226: then res := detmat;
! 227: else res := -detmat;
! 228: end if;
! 229: return res;
! 230: end Fine_Mixed_Volume;
! 231:
! 232: function Mixed_Volume ( n : natural; mix : Vector;
! 233: mic : Mixed_Cell ) return natural is
! 234:
! 235: -- ALGORITHM :
! 236: -- First check if the cell has a refinement, if so, then use it,
! 237: -- if not, then check if the cell is fine mixed.
! 238: -- If the cell is fine mixed, only a determinant needs to be computed,
! 239: -- otherwise the cell will be refined.
! 240:
! 241: res : natural;
! 242:
! 243: begin
! 244: if (mic.sub /= null) and then not Is_Null(mic.sub.all)
! 245: then res := Mixed_Volume_Computation.Mixed_Volume(n,mix,mic.sub.all);
! 246: elsif Is_Fine(mix,mic)
! 247: then res := Fine_Mixed_Volume(n,mix,mic);
! 248: else declare
! 249: rcell : Array_of_Lists(1..n) := Reduced_Supports(n,mix,mic);
! 250: begin
! 251: res := Mixed_Volume_Computation.Mixed_Volume(n,rcell);
! 252: Deep_Clear(rcell);
! 253: end;
! 254: end if;
! 255: return res;
! 256: end Mixed_Volume;
! 257:
! 258: function Mixed_Volume ( n : natural; mix : Vector;
! 259: mixsub : Mixed_Subdivision ) return natural is
! 260:
! 261: res : natural := 0;
! 262: tmp : Mixed_Subdivision := mixsub;
! 263:
! 264: begin
! 265: while not Is_Null(tmp) loop
! 266: res := res + Mixed_Volume(n,mix,Head_Of(tmp));
! 267: tmp := Tail_Of(tmp);
! 268: end loop;
! 269: return res;
! 270: end Mixed_Volume;
! 271:
! 272: procedure Mixed_Volume ( n : in natural; mix : in Vector;
! 273: mic : in out Mixed_Cell; mv : out natural ) is
! 274: begin
! 275: if (mic.sub /= null) and then not Is_Null(mic.sub.all)
! 276: then mv := Mixed_Volume_Computation.Mixed_Volume(n,mix,mic.sub.all);
! 277: elsif Is_Fine(mix,mic)
! 278: then mv := Fine_Mixed_Volume(n,mix,mic);
! 279: else -- NOTE : keep the same type of mixture!
! 280: declare
! 281: rcell : Array_of_Lists(1..n) := Reduced_Supports(n,mix,mic);
! 282: lifted : Array_of_Lists(mix'range);
! 283: mixsub : Mixed_Subdivision;
! 284: begin
! 285: Mixed_Volume_Computation.Mixed_Volume
! 286: (n,mix,rcell,lifted,mixsub,mv);
! 287: mic.sub := new Mixed_Subdivision'(mixsub);
! 288: Deep_Clear(rcell); Deep_Clear(lifted);
! 289: end;
! 290: end if;
! 291: end Mixed_Volume;
! 292:
! 293: procedure Mixed_Volume ( n : in natural; mix : in Vector;
! 294: mixsub : in out Mixed_Subdivision;
! 295: mv : out natural ) is
! 296:
! 297: res : natural := 0;
! 298: tmp : Mixed_Subdivision := mixsub;
! 299:
! 300: begin
! 301: while not Is_Null(tmp) loop
! 302: declare
! 303: mic : Mixed_Cell := Head_Of(tmp);
! 304: mmv : natural;
! 305: begin
! 306: Mixed_Volume(n,mix,mic,mmv);
! 307: Set_Head(tmp,mic);
! 308: res := res + mmv;
! 309: end;
! 310: tmp := Tail_Of(tmp);
! 311: end loop;
! 312: mv := res;
! 313: end Mixed_Volume;
! 314:
! 315: -- MIXED VOLUME COMPUTATIONS BASED ON SUPPORTS :
! 316:
! 317: function Mixed_Volume ( n : natural; supports : Array_of_Lists )
! 318: return natural is
! 319: mv : natural;
! 320: mix,perm : Link_to_Vector;
! 321: permsupp : Array_of_Lists(supports'range);
! 322:
! 323: begin
! 324: Copy(supports,permsupp);
! 325: Compute_Mixture(permsupp,mix,perm);
! 326: mv := Mixed_Volume(n,mix.all,permsupp);
! 327: Clear(mix); Clear(perm); Deep_Clear(permsupp);
! 328: return mv;
! 329: end Mixed_Volume;
! 330:
! 331: function Mixed_Volume ( file : file_type; n : natural;
! 332: supports : Array_of_Lists ) return natural is
! 333: mv : natural;
! 334: mix,perm : Link_to_Vector;
! 335: permsupp : Array_of_Lists(supports'range);
! 336:
! 337: begin
! 338: Copy(supports,permsupp);
! 339: Compute_Mixture(permsupp,mix,perm);
! 340: mv := Mixed_Volume(file,n,mix.all,permsupp);
! 341: Clear(mix); Clear(perm); Deep_Clear(permsupp);
! 342: return mv;
! 343: end Mixed_Volume;
! 344:
! 345: function Mixed_Volume ( n : natural; mix : Vector;
! 346: supports : Array_of_Lists ) return natural is
! 347:
! 348: res : natural;
! 349: mixsub : Mixed_Subdivision;
! 350: lifted : Array_of_Lists(mix'range);
! 351:
! 352: begin
! 353: Mixed_Volume_Computation.Mixed_Volume(n,mix,supports,lifted,mixsub,res);
! 354: Deep_Clear(lifted);
! 355: Shallow_Clear(mixsub);
! 356: return res;
! 357: end Mixed_Volume;
! 358:
! 359: procedure Mixed_Volume
! 360: ( n : in natural; mix : in Vector; supports : in Array_of_Lists;
! 361: lifted : out Array_of_Lists;
! 362: mixsub : out Mixed_Subdivision; mv : out natural ) is
! 363:
! 364: low : constant Vector := (mix'range => 0);
! 365: upp : constant Vector := Adaptive_Lifting(supports);
! 366: nbsucc,nbfail : Standard_Floating_Vectors.Vector(mix'range)
! 367: := (mix'range => 0.0);
! 368: liftsupp : Array_of_Lists(mix'range);
! 369: sub : Mixed_Subdivision;
! 370:
! 371: begin
! 372: Mixed_Coherent_Subdivision(n,mix,supports,false,low,upp,liftsupp,
! 373: nbsucc,nbfail,sub);
! 374: Mixed_Volume(n,mix,sub,mv);
! 375: lifted := liftsupp;
! 376: mixsub := sub;
! 377: end Mixed_Volume;
! 378:
! 379: function Mixed_Volume ( file : file_type; n : natural; mix : Vector;
! 380: supports : Array_of_Lists ) return natural is
! 381:
! 382: res : natural;
! 383: mixsub : Mixed_Subdivision;
! 384: lifted : Array_of_Lists(mix'range);
! 385:
! 386: begin
! 387: Mixed_Volume_Computation.Mixed_Volume
! 388: (file,n,mix,supports,lifted,mixsub,res);
! 389: Deep_Clear(lifted);
! 390: Shallow_Clear(mixsub);
! 391: return res;
! 392: end Mixed_Volume;
! 393:
! 394: procedure Mixed_Volume
! 395: ( file : in file_type; n : in natural;
! 396: mix : in Vector; supports : in Array_of_Lists;
! 397: lifted : out Array_of_Lists;
! 398: mixsub : out Mixed_Subdivision; mv : out natural ) is
! 399:
! 400: low : constant Vector := (mix'range => 0);
! 401: upp : constant Vector := Adaptive_Lifting(supports);
! 402: sub : Mixed_Subdivision;
! 403: nbsucc,nbfail : Standard_Floating_Vectors.Vector(mix'range)
! 404: := (mix'range => 0.0);
! 405: liftsupp : Array_of_Lists(mix'range);
! 406:
! 407: begin
! 408: Mixed_Coherent_Subdivision
! 409: (n,mix,supports,false,low,upp,liftsupp,nbsucc,nbfail,sub);
! 410: lifted := liftsupp;
! 411: mixsub := sub;
! 412: put(file,liftsupp,n,mix,sub,mv);
! 413: end Mixed_Volume;
! 414:
! 415: end Mixed_Volume_Computation;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>