Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/contributions_to_mixed_volume.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
2: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
3: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
4: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
5: with Vertices; use Vertices;
6: with Inner_Normal_Cones; use Inner_Normal_Cones;
7: with Normal_Cone_Intersections; use Normal_Cone_Intersections;
8:
9: package body Contributions_to_Mixed_Volume is
10:
11: -- AUXILIAIRIES TO CONSTRUCT THE FACETS :
12:
13: procedure Copy_Remove ( l : in out List; x : in Vector ) is
14:
15: -- DESCRIPTION :
16: -- Replaces the list by a copy of it, without the point x.
17:
18: tmp,res,res_last : List;
19:
20: begin
21: tmp := l;
22: while not Is_Null(tmp) loop
23: declare
24: pt : Link_to_Vector := Head_Of(tmp);
25: begin
26: if not Equal(pt.all,x)
27: then Append(res,res_last,pt.all);
28: end if;
29: end;
30: tmp := Tail_Of(tmp);
31: end loop;
32: Copy(res,l); Deep_Clear(l);
33: end Copy_Remove;
34:
35: function Vertex_Points ( l : Array_of_Lists ) return Array_of_Lists is
36:
37: -- DESCRIPTION :
38: -- returns for each list the list of the vertex points.
39:
40: res : Array_of_Lists(l'range);
41:
42: begin
43: for i in l'range loop
44: res(i) := Vertex_Points(l(i));
45: end loop;
46: return res;
47: end Vertex_Points;
48:
49: procedure Copy ( f1 : in Array_of_Faces; f2 : in out Array_of_Faces ) is
50:
51: -- DESCRIPTION :
52: -- Copies the array f1 into the array f2.
53:
54: begin
55: for i in f1'range loop
56: Deep_Copy(f1(i),f2(i));
57: end loop;
58: end Copy;
59:
60: function Create_Facets ( n : natural; l : List; x : Vector ) return Faces is
61:
62: -- DESCRIPTION :
63: -- Returns a list of all facets of conv(l), that all contain x.
64: -- First it will be checked whether x belongs to l or not.
65:
66: res : Faces;
67: wrk : List;
68: lx : Link_to_Vector;
69:
70: begin
71: if Is_In(l,x)
72: then res := Create(n-1,n,l,x);
73: else wrk := l;
74: lx := new Vector'(x);
75: Construct(lx,wrk);
76: res := Create(n-1,n,wrk,x);
77: end if;
78: return res;
79: end Create_Facets;
80:
81: function All_Facets ( n : natural; l : Array_of_Lists )
82: return Array_of_Faces is
83:
84: -- DESCRIPTION :
85: -- Returns all facets of all sets in l.
86:
87: res : Array_of_Faces(l'range);
88:
89: begin
90: for i in l'range loop
91: res(i) := Create(n-1,n,l(i));
92: end loop;
93: return res;
94: end All_Facets;
95:
96: -- DETERMINE ZERO CONTRIBUTION BASED ON INTERSECTION MATRIX :
97:
98: function Exhaustive_Zero_Contribution
99: ( pts : Array_of_Lists; g : List; i : natural ) return boolean is
100:
101: -- DESCRIPTION :
102: -- Creates an intersection matrix, based on the list of generators of
103: -- the normal cone of a points in the ith component of pts.
104:
105: res : boolean;
106: n1 : constant natural := pts'length - 1;
107: mg : constant natural := Length_Of(g);
108: nc : constant natural := Number_of_Cones(pts,i);
109: ima : Intersection_Matrix(n1,mg,nc);
110:
111: begin
112: ima := Create(pts,g,i);
113: return Contained_in_Union(pts,i,g,ima);
114: end Exhaustive_Zero_Contribution;
115:
116: -- THE CRITERION :
117:
118: function Simple_Zero_Contribution
119: ( pts : Array_of_Lists; x : Vector; i : natural )
120: return boolean is
121:
122: res : boolean := false;
123: g : List := Generators(pts(i),x);
124:
125: begin
126: for j in pts'range loop
127: if j /= i
128: then res := Contained_in_Cone(pts(j),g);
129: end if;
130: exit when res;
131: end loop;
132: Deep_Clear(g);
133: return res;
134: end Simple_Zero_Contribution;
135:
136: function Simple_Zero_Contribution
137: ( pts : Array_of_Lists; ifacets : Faces;
138: x : Vector; i : natural ) return boolean is
139:
140: g : List := Generators(pts(i),ifacets,x);
141: res : boolean := false;
142:
143: begin
144: for j in pts'range loop
145: if j /= i
146: then res := Contained_in_Cone(pts(j),g);
147: end if;
148: exit when res;
149: end loop;
150: Deep_Clear(g);
151: return res;
152: end Simple_Zero_Contribution;
153:
154: function Exhaustive_Zero_Contribution
155: ( pts : Array_of_Lists;
156: x : Vector; i : natural ) return boolean is
157:
158: n : constant natural := x'length;
159: res : boolean := false;
160:
161: begin
162: if Length_Of(pts(i)) > n
163: then declare
164: f : Faces := Create_Facets(n,pts(i),x);
165: begin
166: res := Exhaustive_Zero_Contribution(pts,f,x,i);
167: Clear(f);
168: end;
169: else declare
170: g : List := Generators(pts(i),x);
171: begin
172: res := Exhaustive_Zero_Contribution(pts,g,i);
173: end;
174: end if;
175: return res;
176: end Exhaustive_Zero_Contribution;
177:
178: function Exhaustive_Zero_Contribution
179: ( pts : Array_of_Lists; ifacets : Faces;
180: x : Vector; i : natural ) return boolean is
181:
182: g : List;
183:
184: begin
185: if not Is_Null(ifacets)
186: then g := Generators(pts(i),ifacets,x);
187: else g := Generators(pts(i),x);
188: end if;
189: return Exhaustive_Zero_Contribution(pts,g,i);
190: end Exhaustive_Zero_Contribution;
191:
192: -- SWEEPING THROUGH THE POINT LISTS :
193:
194: function Simple_Sweep ( pts : Array_of_Lists ) return Array_of_Lists is
195:
196: n : constant natural := Head_Of(pts(pts'first))'length;
197: afa : Array_of_Faces(pts'range) := All_Facets(n,pts);
198:
199: begin
200: return Simple_Sweep(pts,afa);
201: end Simple_Sweep;
202:
203: function Simple_Sweep ( pts : Array_of_Lists; facets : Array_of_Faces )
204: return Array_of_Lists is
205:
206: res,res_last,points : Array_of_Lists(pts'range);
207: -- wrkfacets : Array_of_Faces(facets'range);
208:
209: -- SAFETY MODE : checks whether mixed volume does not decrease
210: -- n : constant natural := pts'last;
211: -- mix : constant Vector := (1..n => 1);
212: -- mv : constant natural := Mixed_Volume(n,mix,pts);
213:
214: begin
215: -- Copy(facets,wrkfacets);
216: points := Vertex_Points(pts); -- instead of: Copy(pts,points);
217: for i in points'range loop
218: declare
219: tmp : constant VecVec := Shallow_Create(points(i));
220: begin
221: for j in tmp'range loop
222: declare
223: x : constant Vector := tmp(j).all;
224: -- f : Faces := Extract_Faces(wrkfacets(i),x);
225: begin
226: -- if not Simple_Zero_Contribution(points,f,x,i)
227: if not Simple_Zero_Contribution(points,x,i)
228: then Append(res(i),res_last(i),x);
229: else Remove(points(i),x);
230: -- SAFETY MODE :
231: -- if mv > Mixed_Volume(n,mix,points)
232: -- then put_line("BUG at points : "); put(points);
233: -- put("for the vector : "); put(x); new_line;
234: -- put(" at component "); put(i,1); new_line;
235: -- raise PROGRAM_ERROR;
236: -- end if;
237: -- Clear(wrkfacets(i));
238: -- wrkfacets(i) := Create(x'length-1,x'length,points(i));
239: end if;
240: end;
241: end loop;
242: end;
243: Copy(res(i),points(i));
244: end loop;
245: Deep_Clear(points);
246: return res;
247: end Simple_Sweep;
248:
249: function Exhaustive_Sweep ( pts : Array_of_Lists ) return Array_of_Lists is
250:
251: n : constant natural := Head_Of(pts(pts'first))'length;
252: afa : Array_of_Faces(pts'range) := All_Facets(n,pts);
253:
254: begin
255: return Exhaustive_Sweep(pts,afa);
256: end Exhaustive_Sweep;
257:
258: function Exhaustive_Sweep ( pts : Array_of_Lists; facets : Array_of_Faces )
259: return Array_of_Lists is
260:
261: res,res_last,points : Array_of_Lists(pts'range);
262: -- wrkfacets : Array_of_Faces(facets'range);
263:
264: -- SAFETY MODE : checks whether mixed volume does not decrease
265: -- n : constant natural := pts'last;
266: -- mix : constant Vector := (1..n => 1);
267: -- mv : constant natural := Mixed_Volume(n,mix,pts);
268:
269: begin
270: -- Copy(facets,wrkfacets);
271: points := Vertex_Points(pts); -- instead of: Copy(pts,points);
272: for i in points'range loop
273: declare
274: tmp : constant VecVec := Shallow_Create(points(i));
275: begin
276: for j in tmp'range loop
277: declare
278: x : constant Vector := tmp(j).all;
279: -- f : Faces := Extract_Faces(wrkfacets(i),x);
280: begin
281: -- if not Exhaustive_Zero_Contribution(points,f,x,i)
282: if not Exhaustive_Zero_Contribution(points,x,i)
283: then Append(res(i),res_last(i),x);
284: else Remove(points(i),x);
285: -- SAFETY MODE :
286: -- if mv > Mixed_Volume(n,mix,points)
287: -- then put_line("BUG at points : "); put(points);
288: -- put("for the vector : "); put(x); new_line;
289: -- put(" at component "); put(i,1); new_line;
290: -- raise PROGRAM_ERROR;
291: -- end if;
292: -- Clear(wrkfacets(i));
293: -- wrkfacets(i) := Create(x'length-1,x'length,points(i));
294: end if;
295: end;
296: end loop;
297: end;
298: Copy(res(i),points(i));
299: end loop;
300: Deep_Clear(points);
301: return res;
302: end Exhaustive_Sweep;
303:
304: end Contributions_to_Mixed_Volume;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>