Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/mixed_volume_computation.adb, Revision 1.1.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>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mixed_volume_computation.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