Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/inner_normal_cones.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
2: with Transformations; use Transformations;
3: with Integer_Vectors_Utilities; use Integer_Vectors_Utilities;
4: with Lists_of_Vectors_Utilities; use Lists_of_Vectors_Utilities;
5: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
6: with Vertices; use Vertices;
7:
8: package body Inner_Normal_Cones is
9:
10: -- NOTE ON THE IMPLEMENTATION :
11: -- There is a safety mode in which after each computation of the generators,
12: -- it is automatically tested whether the generators satisfy the
13: -- inequalities of the normal cone. If not, then the bug is reported and
14: -- a program_error is raised.
15:
16: -- AUXILIAIRIES :
17:
18: function Lower ( a1,b1,a2,b2 : integer ) return boolean is
19:
20: -- DESCRIPTION :
21: -- Returns true if a1/b1 < a2/b2, false otherwise.
22:
23: -- REQUIRED : b1 /= 0 and b2 /= 0.
24:
25: begin
26: if b1*b2 > 0
27: then return (a1*b2 - a2*b1 < 0);
28: else return (a1*b2 - a2*b1 > 0);
29: end if;
30: end Lower;
31:
32: function Higher ( a1,b1,a2,b2 : integer ) return boolean is
33:
34: -- DESCRIPTION :
35: -- Returns true if a1/b1 > a2/b2, false otherwise.
36:
37: -- REQUIRED : b1 /= 0 and b2 /= 0.
38:
39: begin
40: if b1*b2 > 0
41: then return (a1*b2 - a2*b1 > 0);
42: else return (a1*b2 - a2*b1 < 0);
43: end if;
44: end Higher;
45:
46: function Compute_Facets ( l : List; x : Vector ) return Faces is
47:
48: -- DESCRIPTION :
49: -- Returns all facets of conv(l) that all contain x.
50: -- Checks whether x belongs to the list or not.
51:
52: res : Faces;
53: n : constant natural := x'length;
54:
55: begin
56: if Is_In(l,x)
57: then res := Create(n-1,n,l,x);
58: else declare
59: wrk : List := l;
60: lx : Link_to_Vector;
61: begin
62: lx := new Vector'(x);
63: Construct(lx,wrk);
64: res := Create(n-1,n,wrk,x);
65: end;
66: end if;
67: return res;
68: end Compute_Facets;
69:
70: procedure Shifted_Points ( l : in List; x : in Vector;
71: res,res_last : in out List ) is
72:
73: -- DESCRIPTION :
74: -- Appends the list of shifted points w.r.t. x to the list res.
75:
76: tmp : List := l;
77:
78: begin
79: tmp := l;
80: while not Is_Null(tmp) loop
81: Append_Diff(res,res_last,Head_Of(tmp).all-x);
82: tmp := Tail_Of(tmp);
83: end loop;
84: end Shifted_Points;
85:
86: function In_Hull ( l : List; x : Vector ) return boolean is
87:
88: -- DESCRIPTION :
89: -- Returns true if the vector x belongs the convex hull spanned by
90: -- the points in l minus the point x itself.
91:
92: res : boolean;
93: lx : List;
94:
95: begin
96: Copy(l,lx);
97: Remove(lx,x);
98: res := Is_In_Hull(x,lx);
99: Clear(lx);
100: return res;
101: end In_Hull;
102:
103: -- AUXILIAIRIES TO CONSTRUCT THE NORMALS :
104:
105: procedure Normal ( m : in out Matrix; v : out Vector ) is
106:
107: -- DESCRIPTION :
108: -- Computes the normal to the vectors stored in the rows of the matrix m.
109:
110: -- REQUIRED : the matrix m is square!
111:
112: res : Vector(m'range(2));
113:
114: begin
115: Upper_Triangulate(m); Scale(m);
116: res := (res'range => 0);
117: Solve0(m,res);
118: v := res;
119: end Normal;
120:
121: procedure Orientate_Normal ( l : in List; f : in Face;
122: normal : in out Vector ) is
123:
124: -- DESCRIPTION :
125: -- Orientates the normal of the face such that it becomes an inner normal
126: -- w.r.t. the points in the list l.
127:
128: tmp : List := l;
129: ip : constant integer := f(f'first).all*normal;
130: pt : Vector(Head_Of(l)'range);
131: done : boolean := false;
132: ptip : integer;
133:
134: begin
135: while not Is_Null(tmp) loop
136: pt := Head_Of(tmp).all;
137: if not Is_In(f,pt)
138: then ptip := pt*normal;
139: if ptip /= ip
140: then if ptip < ip
141: then normal := -normal;
142: end if;
143: done := true;
144: end if;
145: end if;
146: exit when done;
147: tmp := Tail_Of(tmp);
148: end loop;
149: end Orientate_Normal;
150:
151: function Inner_Normal ( l : List; facet : Face ) return Vector is
152:
153: -- DESCRIPTION :
154: -- Returns the inner normal to the facet.
155:
156: res,fst : Vector(facet(facet'first)'range);
157: m : Matrix(facet'first+1..facet'last,facet(facet'first)'range);
158:
159: begin
160: fst := facet(facet'first).all;
161: for i in m'first(1)..m'last(1) loop
162: for j in m'range(2) loop
163: m(i,j) := facet(i)(j) - fst(j);
164: end loop;
165: end loop;
166: Normal(m,res);
167: Orientate_Normal(l,facet,res);
168: return res;
169: end Inner_Normal;
170:
171: procedure Inner_Normals ( l : in List; facets : in Faces;
172: iv,iv_last : in out List ) is
173:
174: -- DESCRIPTION :
175: -- Returns the list of all inner normals to the facets of conv(l).
176: -- If there is only one facet, then both `inner' and `outer' normal
177: -- will be returned, because there is nothing to orientate with.
178: -- This means that also the negation of the normal will be in the list iv.
179:
180: tmp : Faces := facets;
181:
182: begin
183: while not Is_Null(tmp) loop
184: declare
185: f : Face := Head_Of(tmp);
186: v : constant Vector := Inner_Normal(l,Head_Of(tmp));
187: begin
188: Append(iv,iv_last,v);
189: if Length_Of(facets) = 1
190: then Append(iv,iv_last,-v);
191: end if;
192: end;
193: tmp := Tail_Of(tmp);
194: end loop;
195: end Inner_Normals;
196:
197: function Number_of_Rows ( l : List ) return natural is
198:
199: -- DESCRIPTION :
200: -- Returns Maximum(dimension,Length_Of(l)).
201:
202: -- REQUIRED : not Is_Null(l).
203:
204: n : constant natural := Head_Of(l)'length;
205: m : constant natural := Length_Of(l);
206:
207: begin
208: if n > m
209: then return n;
210: else return m;
211: end if;
212: end Number_of_Rows;
213:
214: function Normal ( l : List ) return Vector is
215:
216: -- DESCRIPTION :
217: -- Return a normal to the vectors in the list.
218:
219: -- REQUIRED : Length_Of(l) <= n = Head_Of(l)'length
220: -- or the points in l all lie on the same facet.
221:
222: fst : Link_to_Vector := Head_Of(l);
223: m : Matrix(1..Number_of_Rows(l),fst'range);
224: cnt : natural := 0;
225: res : Vector(fst'range);
226: tmp : List := Tail_Of(l);
227:
228: begin
229: while not Is_Null(tmp) loop
230: res := Head_Of(tmp).all;
231: cnt := cnt+1;
232: for i in res'range loop
233: m(cnt,i) := res(i) - fst(i);
234: end loop;
235: tmp := Tail_Of(tmp);
236: end loop;
237: for i in cnt+1..m'last(1) loop
238: for j in m'range(2) loop
239: m(i,j) := 0;
240: end loop;
241: end loop;
242: Normal(m,res);
243: return res;
244: end Normal;
245:
246: procedure Normals ( l : in List; x : in Vector; iv,iv_last : in out List ) is
247:
248: -- DESCRIPTION :
249: -- Computes a list of all normals to the vectors in the list.
250:
251: -- REQUIRED : Length_Of(l) <= n = x'length.
252: -- or all points in l lie on the same facet.
253:
254: nor : constant Vector := Normal(l);
255: piv : natural := Pivot(nor);
256: pivnor : Vector(nor'range); -- normal with pivnor(piv) > 0
257:
258: begin
259: if piv <= nor'last
260: then if x'length > 2
261: then if nor(piv) < 0
262: then pivnor := -nor;
263: else pivnor := nor;
264: end if;
265: declare
266: t : Transfo := Build_Transfo(pivnor,piv);
267: tx : constant Vector := Reduce(t*x,piv);
268: tt : Transfo := Transpose(t);
269: tl : List := Transform_and_Reduce(t,piv,l);
270: begin
271: if Length_Of(tl) <= nor'length-1
272: then Normals(tl,tx,iv,iv_last);
273: else declare
274: facets : Faces := Compute_Facets(tl,tx);
275: begin
276: if Length_Of(facets) > 1
277: then Inner_Normals(tl,facets,iv,iv_last);
278: else Normals(tl,tx,iv,iv_last);
279: end if;
280: end;
281: end if;
282: Insert_and_Transform(iv,piv,0,tt);
283: iv_last := Pointer_to_Last(iv);
284: Clear(t); Deep_Clear(tl); Clear(tt);
285: end;
286: end if;
287: Append(iv,iv_last,nor); Append(iv,iv_last,-nor);
288: end if;
289: end Normals;
290:
291: -- CONSTRUCTORS FOR PRIMAL REPRESENTATION :
292:
293: function Generators ( l : List; facets : Faces; x : Vector ) return List is
294:
295: res,res_last : List;
296:
297: begin
298: if Length_Of(facets) > 1
299: then Inner_Normals(l,facets,res,res_last); -- compute inner normals
300: elsif In_Hull(l,x)
301: then return res; -- empty normal cone
302: else Normals(l,x,res,res_last); -- compute normals
303: if Length_Of(l) = 2
304: then Shifted_Points(l,x,res,res_last);
305: end if;
306: end if;
307: -- SAFETY MODE :
308: -- if not Contained_in_Cone(l,x,res)
309: -- then put_line("Bug raised for computing generators for list"); put(l);
310: -- put(" and the point : "); put(x); new_line;
311: -- raise PROGRAM_ERROR;
312: -- end if;
313: return res;
314: end Generators;
315:
316: function Generators ( l : List; x : Vector ) return List is
317:
318: res,res_last : List;
319: n : constant natural := x'length;
320: facets : Faces;
321:
322: begin
323: if In_Hull(l,x)
324: then return res; -- empty normal cone
325: else if Length_Of(l) > n
326: then facets := Compute_Facets(l,x);
327: end if;
328: if Length_Of(facets) > 1
329: then res := Generators(l,facets,x);
330: else Normals(l,x,res,res_last);
331: if Length_Of(l) = 2
332: then Shifted_Points(l,x,res,res_last);
333: end if;
334: end if;
335: -- SAFETY MODE :
336: -- if not Contained_in_Cone(l,x,res)
337: -- then put_line("Bug raised for computing generators for list");
338: -- put(l); put(" and the point : "); put(x); new_line;
339: -- raise PROGRAM_ERROR;
340: -- end if;
341: return res;
342: end if;
343: end Generators;
344:
345: -- CONSTRUCTOR FOR DUAL REPRESENTATION :
346:
347: function Inner_Normal_Cone ( l : List; x : Vector ) return Matrix is
348:
349: res : Matrix(x'range,1..Length_Of(l)-1);
350: -- CAUTION : Is_In(l,x) must hold !!
351: y : Vector(x'range);
352: cnt : natural := 0;
353: tmp : List := l;
354:
355: begin
356: while not Is_Null(tmp) loop
357: y := Head_Of(tmp).all;
358: if not Equal(x,y) -- avoid trivial inequality
359: then cnt := cnt+1;
360: for i in x'range loop
361: res(i,cnt) := y(i) - x(i);
362: end loop;
363: end if;
364: tmp := Tail_Of(tmp);
365: end loop;
366: return res;
367: end Inner_Normal_Cone;
368:
369: function Included_Normal_Cone ( gx : List; x : Vector ) return Matrix is
370:
371: res : Matrix(x'first-1..x'last,1..Length_Of(gx));
372: tmp : List := gx;
373: v : Vector(x'range);
374:
375: begin
376: for j in res'range(2) loop
377: v := Head_Of(tmp).all;
378: res(res'first(1),j) := x*v;
379: for i in res'first(1)+1..res'last(1) loop
380: res(i,j) := v(i);
381: end loop;
382: tmp := Tail_Of(tmp);
383: end loop;
384: return res;
385: end Included_Normal_Cone;
386:
387: -- PRIMITIVE SELECTORS :
388:
389: function Evaluate ( m : Matrix; i : integer; v : Vector ) return integer is
390:
391: sum : integer := 0;
392:
393: begin
394: for j in v'range loop
395: sum := sum + m(j,i)*v(j);
396: end loop;
397: return sum;
398: end Evaluate;
399:
400: function Satisfies ( m : Matrix; i : integer; v : Vector ) return boolean is
401: begin
402: if m'first(1) = v'first-1
403: then return (Evaluate(m,i,v) >= m(0,i));
404: else return (Evaluate(m,i,v) >= 0);
405: end if;
406: end Satisfies;
407:
408: function Strictly_Satisfies ( m : Matrix; i : integer; v : Vector )
409: return boolean is
410: begin
411: if m'first(1) = v'first-1
412: then return (Evaluate(m,i,v) > m(0,i));
413: else return (Evaluate(m,i,v) > 0);
414: end if;
415: end Strictly_Satisfies;
416:
417: function Satisfies ( m : Matrix; v : Vector ) return boolean is
418: begin
419: for i in m'range(2) loop
420: if not Satisfies(m,i,v)
421: then return false;
422: end if;
423: end loop;
424: return true;
425: end Satisfies;
426:
427: function Strictly_Satisfies ( m : Matrix; v : Vector ) return boolean is
428: begin
429: for i in m'range(2) loop
430: if not Strictly_Satisfies(m,i,v)
431: then return false;
432: end if;
433: end loop;
434: return true;
435: end Strictly_Satisfies;
436:
437: function Satisfies ( m : Matrix; v : List ) return boolean is
438:
439: tmp : List := v;
440:
441: begin
442: while not Is_Null(tmp) loop
443: if not Satisfies(m,Head_Of(tmp).all)
444: then return false;
445: else tmp := Tail_Of(tmp);
446: end if;
447: end loop;
448: return true;
449: end Satisfies;
450:
451: function Strictly_Satisfies ( m : Matrix; v : List ) return boolean is
452:
453: tmp : List := v;
454:
455: begin
456: while not Is_Null(tmp) loop
457: if not Strictly_Satisfies(m,Head_Of(tmp).all)
458: then return false;
459: else tmp := Tail_Of(tmp);
460: end if;
461: end loop;
462: return true;
463: end Strictly_Satisfies;
464:
465: -- SECONDARY SELECTORS :
466:
467: function Satisfy ( m : Matrix; l : List ) return List is
468:
469: res,res_last,tmp : List;
470:
471: begin
472: tmp := l;
473: while not Is_Null(tmp) loop
474: if Satisfies(m,Head_Of(tmp).all)
475: then Append(res,res_last,Head_Of(tmp).all);
476: end if;
477: tmp := Tail_Of(tmp);
478: end loop;
479: return res;
480: end Satisfy;
481:
482: function Strictly_Satisfy ( m : Matrix; l : List ) return List is
483:
484: res,res_last,tmp : List;
485:
486: begin
487: tmp := l;
488: while not Is_Null(tmp) loop
489: if Strictly_Satisfies(m,Head_Of(tmp).all)
490: then Append(res,res_last,Head_Of(tmp).all);
491: end if;
492: tmp := Tail_Of(tmp);
493: end loop;
494: return res;
495: end Strictly_Satisfy;
496:
497: function Contained_in_Cone ( l : List; x : Vector; v : List )
498: return boolean is
499:
500: ineq : Matrix(x'range,1..Length_Of(l)-1) := Inner_Normal_Cone(l,x);
501:
502: begin
503: return Satisfies(ineq,v);
504: end Contained_in_Cone;
505:
506: function Strictly_Contained_in_Cone ( l : List; x : Vector; v : List )
507: return boolean is
508:
509: ineq : Matrix(x'range,1..Length_Of(l)-1) := Inner_Normal_Cone(l,x);
510:
511: begin
512: return Strictly_Satisfies(ineq,v);
513: end Strictly_Contained_in_Cone;
514:
515: function Contained_in_Cone ( l,v : List ) return boolean is
516:
517: tmp : List := l;
518:
519: begin
520: while not Is_Null(tmp) loop
521: if Contained_in_Cone(l,Head_Of(tmp).all,v)
522: then -- put("true for vector : "); put(Head_Of(tmp).all); new_line;
523: return true;
524: else tmp := Tail_Of(tmp);
525: end if;
526: end loop;
527: return false;
528: end Contained_in_Cone;
529:
530: function Strictly_Contained_in_Cone ( l,v : List ) return boolean is
531:
532: tmp : List := l;
533:
534: begin
535: while not Is_Null(tmp) loop
536: if Strictly_Contained_in_Cone(l,Head_Of(tmp).all,v)
537: then return true;
538: else tmp := Tail_Of(tmp);
539: end if;
540: end loop;
541: return false;
542: end Strictly_Contained_in_Cone;
543:
544: -- CONCERNING THE UNION OF CONES :
545:
546: function In_Union ( v1,v2 : Vector; k1,k2 : Matrix ) return boolean is
547:
548: -- ALGORITHM : consider x(t) = (1-t)*v1 + t*v2, for t in [0,1]
549: -- and compute t1: smallest t such that Satisfies(k1,x(t))
550: -- and t2: largest t such that Satisfies(k2,x(t)).
551: -- Then t1 >= t2 makes the test returning true, false otherwise.
552:
553: diff : constant Vector := v1-v2;
554: t1a,t1b,t2a,t2b,a,b : integer; -- t1 = t1a/t1b and t2 = t2a/t2b
555:
556: begin
557: t1a := 1; t1b := 1; -- initially t1 = 1 = 1/1
558: for i in k1'range(2) loop
559: b := Evaluate(k1,i,diff);
560: if b /= 0
561: then a := Evaluate(k1,i,v1);
562: if Lower(a,b,t1a,t1b)
563: then t1a := a; t1b := b;
564: end if;
565: end if;
566: end loop;
567: t2a := 0; t2b := 1; -- initially t2 = 0 = 0/1
568: for i in k2'range(2) loop
569: b := Evaluate(k2,i,diff);
570: if b /= 0
571: then a := Evaluate(k2,i,v1);
572: if Higher(a,b,t2a,t2b)
573: then t2a := a; t2b := b;
574: end if;
575: end if;
576: end loop;
577: -- put("In_Union for "); put(v1); put(" and "); put(v2);
578: -- put(" : t1 = "); put(t1a,1); put("/"); put(t1b,1);
579: -- put(" : t2 = "); put(t2a,1); put("/"); put(t2b,1); new_line;
580: -- put_line("k1 : "); put(k1); put_line("k2 : "); put(k2);
581: return not Lower(t1a,t1b,t2a,t2b);
582: end In_Union;
583:
584: function In_Union ( v1,v2 : List; k1,k2 : Matrix ) return boolean is
585:
586: tmpv1,tmpv2 : List;
587:
588: begin
589: tmpv1 := v1;
590: while not Is_Null(tmpv1) loop
591: tmpv2 := v2;
592: while not Is_Null(tmpv2) loop
593: if not In_Union(Head_Of(tmpv1).all,Head_Of(tmpv2).all,k1,k2)
594: then return false;
595: else tmpv2 := Tail_Of(tmpv2);
596: end if;
597: end loop;
598: tmpv1 := Tail_Of(tmpv1);
599: end loop;
600: return true;
601: end In_Union;
602:
603: end Inner_Normal_Cones;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to inner_normal_cones.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