![[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|
Return to vertices.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift |
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;