Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/simplices.adb, Revision 1.1.1.1
1.1 maekawa 1: with unchecked_deallocation;
2: with Standard_Integer_Norms; use Standard_Integer_Norms;
3: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
4: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
5: with Transformations; use Transformations;
6:
7: package body Simplices is
8:
9: -- DATA STRUCTURES :
10:
11: type Point is record
12: pt : Link_to_Vector; -- the point
13: si : Simplex; -- by pivoting a new simplex is obtained
14: end record;
15: type Points is array ( integer range <> ) of Point;
16:
17: type Simplex_Rep ( n : natural ) is record
18: nor : vector(1..n); -- normal to the simplex
19: tra : Transfo; -- transformation to triangulate the cell
20: pts : points(1..n); -- vector of points
21: end record;
22:
23: -- AUXILIAIRIES :
24:
25: function Diagonalize ( s : simplex ) return matrix is
26:
27: -- DESCRIPTION :
28: -- Places the vertices of the simplex, shifted w.r.t. the first point,
29: -- in a matrix. Returns the triangulated matrix.
30:
31: a : matrix(1..s.n-1,1..s.n-1);
32: x,y : Vector(1..s.n-1);
33:
34: begin
35: for k in 2..s.n loop
36: x := s.pts(k).pt(x'range) - s.pts(1).pt(x'range);
37: y := s.tra*x;
38: for i in a'range(1) loop
39: a(i,k-1) := y(i);
40: end loop;
41: end loop;
42: return a;
43: end Diagonalize;
44:
45: function Diagonalize ( s : simplex; pt : Vector ) return matrix is
46:
47: -- DESCRIPTION :
48: -- Places the vertices of the simplex, shifted w.r.t. the first point,
49: -- in a matrix. Returns the triangulated matrix.
50:
51: a : matrix(1..s.n-1,1..s.n);
52: x,y : Vector(1..s.n-1);
53:
54: begin
55: for k in 2..s.n loop
56: x := s.pts(k).pt(x'range) - s.pts(1).pt(x'range); y := s.tra*x;
57: for i in a'range(1) loop
58: a(i,k-1) := y(i);
59: end loop;
60: end loop;
61: x := pt(x'range) - s.pts(1).pt(x'range); y := s.tra*x;
62: for i in a'range(1) loop
63: a(i,s.n) := y(i);
64: end loop;
65: return a;
66: end Diagonalize;
67:
68: function Create ( pts : VecVec ) return Transfo is
69:
70: a,l : matrix(pts'first..pts'last-1,pts'first..pts'last-1);
71: x : vector(pts(pts'first)'range);
72:
73: begin
74: for k in pts'first+1..pts'last loop
75: x := pts(k).all - pts(pts'first).all;
76: for i in a'range(1) loop
77: a(i,k-1) := x(i);
78: end loop;
79: end loop;
80: Upper_Triangulate(l,a);
81: return Create(l);
82: end Create;
83:
84: function Create ( pts : VecVec ) return Vector is
85:
86: a : matrix(pts'first..pts'last - 1,pts'range);
87: res : Vector(pts'range);
88:
89: begin
90: for k in a'range(1) loop
91: for i in a'range(2) loop
92: a(k,i) := pts(k+1)(i) - pts(pts'first)(i);
93: end loop;
94: end loop;
95: Upper_Triangulate(a);
96: Scale(a);
97: res := (res'range => 0);
98: Solve0(a,res);
99: Normalize(res);
100: if res(res'last) < 0
101: then return -res;
102: else return res;
103: end if;
104: end Create;
105:
106: -- CREATORS :
107:
108: function Create ( x : VecVec ) return Simplex is
109:
110: n : natural := x'last - x'first + 1;
111: res : Simplex := new Simplex_Rep(n);
112: cnt : natural := x'first;
113:
114: begin
115: for k in res.pts'range loop
116: res.pts(k).pt := x(cnt);
117: cnt := cnt + 1;
118: res.pts(k).si := Null_Simplex;
119: end loop;
120: res.tra := Create(x);
121: res.nor := Create(x);
122: return res;
123: end Create;
124:
125: procedure Update ( s : in out Simplex; x : in Link_to_Vector;
126: k : in natural ) is
127:
128: pts : VecVec(1..s.n);
129: nei : Simplex;
130:
131: begin
132: if s.pts(k).si = Null_Simplex
133: then for i in pts'range loop
134: if i = k
135: then pts(i) := x;
136: else pts(i) := s.pts(i).pt;
137: end if;
138: end loop;
139: nei := Create(pts);
140: s.pts(k).si := nei;
141: nei.pts(k).si := s;
142: end if;
143: end Update;
144:
145: procedure Update ( s : in out Simplex; x : in Link_to_Vector;
146: pos : in Vector ) is
147: begin
148: for k in pos'first..pos'last-1 loop
149: if pos(k)*pos(pos'last) > 0
150: then Update(s,x,k+1);
151: end if;
152: end loop;
153: end Update;
154:
155: procedure Update_One ( s : in out Simplex; x : in Link_to_Vector;
156: pos : in Vector ) is
157:
158: done : boolean := false;
159: nei : Simplex;
160:
161: begin
162: for k in pos'first..pos'last-1 loop
163: if pos(k)*pos(pos'last) > 0
164: then nei := s.pts(k+1).si;
165: if nei /= Null_Simplex
166: then Update_One(nei,x,Position(nei,x.all));
167: else Update(s,x,k+1);
168: end if;
169: done := true;
170: end if;
171: exit when done;
172: end loop;
173: end Update_One;
174:
175: procedure Update_One ( s : in out Simplex; x : in Link_to_Vector;
176: pos : in Vector; news : out Simplex ) is
177:
178: done : boolean := false;
179: nei,newsnei : Simplex;
180:
181: begin
182: -- LOOK FIRST FOR NULL SIMPLEX IN THE DIRECTION TO x :
183: for k in pos'first..pos'last-1 loop
184: if pos(k)*pos(pos'last) > 0
185: then nei := s.pts(k+1).si;
186: if nei = Null_Simplex
187: then Update(s,x,k+1);
188: news := s.pts(k+1).si;
189: done := true;
190: end if;
191: end if;
192: exit when done;
193: end loop;
194: -- WALK FURTHER IN THE DIRECTION TO x :
195: if not done
196: then Update_One(nei,x,Position(nei,x.all),newsnei);
197: if newsnei /= Null_Simplex
198: then news := newsnei;
199: s := nei;
200: end if;
201: end if;
202: end Update_One;
203:
204: procedure Update_All ( s : in out Simplex; x : in Link_to_Vector;
205: pos : in Vector; ancestor : in Simplex ) is
206:
207: nei : Simplex;
208: continue : boolean := true;
209:
210: begin
211: for k in pos'first..pos'last-1 loop
212: if pos(k)*pos(pos'last) > 0
213: then nei := s.pts(k+1).si;
214: if nei /= Null_Simplex
215: then if not Is_Vertex(nei,x.all) and nei /= ancestor
216: then Update_All(nei,x,Position(nei,x.all),s);
217: end if;
218: else Update(s,x,k+1);
219: Process(s.pts(k+1).si,continue);
220: end if;
221: end if;
222: exit when not continue;
223: end loop;
224: end Update_All;
225:
226: procedure Connect ( s1,s2 : in out Simplex ) is
227:
228: neighb : boolean;
229: index1,index2 : natural;
230:
231: begin
232: neighb := true; -- assume they are neighbors
233: index1 := 0;
234: -- SEARCH FOR INDEX OF POINT IN s1 THAT DOES NOT BELONG TO s2 :
235: for k in s1.pts'range loop
236: if not Is_Vertex(s2,s1.pts(k).pt.all)
237: then if (index1 = 0) and (s1.pts(k).si = Null_Simplex)
238: then index1 := k; -- kth point is not common
239: else neighb := false; -- more than one point not common
240: -- or there is already a neighbor
241: end if;
242: end if;
243: exit when not neighb;
244: end loop; -- either not neighb or index1 -> point in s1, not in s2
245: -- SEARCH FOR INDEX OF POINT IN s2 THAT DOES NOT BELONG TO s1 :
246: if neighb
247: then index2 := 0;
248: for k in s2.pts'range loop
249: if not Is_Vertex(s1,s2.pts(k).pt.all)
250: then if (index2 = 0) and (s2.pts(k).si = Null_Simplex)
251: then index2 := k; -- kth point is not common
252: else neighb := false; -- more than one point not common
253: -- or there is already a neighbor
254: end if;
255: end if;
256: exit when not neighb;
257: end loop; -- either no neighb or index2 -> point in s2, not in s1
258: -- CONNECT THE SIMPLICES WITH EACH OTHER :
259: if neighb
260: then s1.pts(index1).si := s2;
261: s2.pts(index2).si := s1;
262: end if;
263: end if;
264: end Connect;
265:
266: procedure Flatten ( s : in out Simplex ) is
267: begin
268: s.nor := (s.nor'range => 0);
269: s.nor(s.n) := 1;
270: for k in s.pts'range loop
271: s.pts(k).pt(s.n) := 0;
272: end loop;
273: end Flatten;
274:
275: -- SELECTORS :
276:
277: function Dimension ( s : Simplex ) return natural is
278: begin
279: return s.n;
280: end Dimension;
281:
282: function Normal ( s : Simplex ) return Vector is
283: begin
284: return s.nor;
285: end Normal;
286:
287: function Is_Flat ( s : Simplex ) return boolean is
288: begin
289: for i in s.nor'first..(s.nor'last-1) loop
290: if s.nor(i) /= 0
291: then return false;
292: end if;
293: end loop;
294: return (s.nor(s.nor'last) = 1);
295: end Is_Flat;
296:
297: function Vertices ( s : Simplex ) return VecVec is
298:
299: res : VecVec(s.pts'range);
300:
301: begin
302: for k in res'range loop
303: res(k) := s.pts(k).pt;
304: end loop;
305: return res;
306: end Vertices;
307:
308: function Vertex ( s : Simplex; k : natural ) return Vector is
309: begin
310: return s.pts(k).pt.all;
311: end Vertex;
312:
313: function Is_Vertex ( s : Simplex; x : Vector ) return boolean is
314: begin
315: for k in s.pts'range loop
316: if s.pts(k).pt.all = x
317: then return true;
318: end if;
319: end loop;
320: return false;
321: end Is_Vertex;
322:
323: function Equal ( s1,s2 : Simplex ) return boolean is
324:
325: found : boolean;
326:
327: begin
328: if s1.nor /= s2.nor -- check if normals are the same
329: then return false;
330: else for k in s1.pts'range loop -- check if vertices are the same
331: -- check whether s1.pts(k).pt.all occurs in s2.pts
332: found := false;
333: for l in s2.pts'range loop
334: if s1.pts(k).pt.all = s2.pts(l).pt.all
335: then found := true;
336: end if;
337: exit when found;
338: end loop;
339: if not found
340: then return false;
341: end if;
342: end loop;
343: return true;
344: end if;
345: end Equal;
346:
347: function Index ( s : Simplex; x : Vector ) return natural is
348: begin
349: for k in s.pts'range loop
350: if s.pts(k).pt.all = x
351: then return k;
352: end if;
353: end loop;
354: return 0;
355: end Index;
356:
357: function Neighbor ( s : Simplex; k : natural ) return Simplex is
358: begin
359: return s.pts(k).si;
360: end Neighbor;
361:
362: function Neighbor ( s : Simplex; k : natural; pos : Vector )
363: return Simplex is
364: begin
365: if pos(k-1)*pos(pos'last) > 0
366: then return s.pts(k).si;
367: else return Null_Simplex;
368: end if;
369: end Neighbor;
370:
371: function Position ( s : Simplex; x : Vector ) return Vector is
372:
373: m : matrix(x'first..x'last-1,x'range);
374: pos : Vector(x'range);
375: res : Vector(0..pos'last);
376:
377: begin
378: -- nbpos := nbpos + 1;
379: -- put("# position computations : "); put(nbpos,1); new_line;
380: -- transform point and simplex
381: m := Diagonalize(s,x);
382: -- solve the system
383: pos := (pos'range => 0);
384: Solve0(m,pos);
385: res(pos'first..pos'last) := pos;
386: res(0) := 0;
387: for k in pos'range loop
388: res(0) := res(0) + pos(k);
389: end loop;
390: res(0) := -res(0);
391: return res;
392: end Position;
393:
394: function Is_In ( s : Simplex; x : Vector ) return boolean is
395:
396: pos : Vector(0..x'last) := Position(s,x);
397:
398: begin
399: return Is_In(s,x,pos);
400: end Is_In;
401:
402: function Is_In ( s : Simplex; x,pos : Vector ) return boolean is
403: begin
404: for k in pos'first..pos'last-1 loop
405: if pos(k)*pos(pos'last) > 0
406: then return false;
407: end if;
408: end loop;
409: return true;
410: end Is_In;
411:
412: function Is_In_All ( s : Simplex; x : Vector ) return boolean is
413:
414: pos : Vector(0..x'last) := Position(s,x);
415:
416: begin
417: return Is_In_All(s,x,pos);
418: end Is_In_All;
419:
420: function Is_In_All ( s : Simplex; x : Vector ) return Simplex is
421:
422: pos : Vector(0..x'last) := Position(s,x);
423:
424: begin
425: return Is_In_All(s,x,pos);
426: end Is_In_All;
427:
428: function Is_In_All ( s : Simplex; x,pos : Vector ) return boolean is
429:
430: ins : boolean := true; -- assumes that x belongs to s
431:
432: begin
433: for k in pos'first..pos'last-1 loop
434: if pos(k)*pos(pos'last) > 0
435: then if s.pts(k+1).si /= Null_Simplex
436: then return Is_In_All(s.pts(k+1).si,x);
437: else ins := false;
438: end if;
439: end if;
440: end loop;
441: return ins;
442: end Is_In_All;
443:
444: function Is_In_All ( s : Simplex; x,pos : Vector ) return Simplex is
445:
446: ins : boolean := true; -- assumes that x belongs to s
447:
448: begin
449: for k in pos'first..pos'last-1 loop
450: if pos(k)*pos(pos'last) > 0
451: then if s.pts(k+1).si /= Null_Simplex
452: then return Is_In_All(s.pts(k+1).si,x);
453: else ins := false;
454: end if;
455: end if;
456: end loop;
457: if ins
458: then return s;
459: else return Null_Simplex;
460: end if;
461: end Is_In_All;
462:
463: procedure Neighbors ( s : in out Simplex; x : in Vector ) is
464:
465: cont : boolean;
466:
467: procedure Neighbors ( s : in out Simplex; x : in Vector;
468: cont : out boolean ) is
469:
470: pos : Vector(0..x'last) := Position(s,x);
471: continue : boolean := true;
472:
473: begin
474: for k in pos'first..pos'last-1 loop
475: if pos(k)*pos(pos'last) > 0
476: then if s.pts(k+1).si /= Null_Simplex
477: then Neighbors(s.pts(k+1).si,x,continue);
478: else Process_Neighbor(s,k+1,continue);
479: end if;
480: end if;
481: exit when not continue;
482: end loop;
483: cont := continue;
484: end Neighbors;
485:
486: begin
487: Neighbors(s,x,cont);
488: end Neighbors;
489:
490: function Volume ( s : Simplex ) return natural is
491:
492: m : constant matrix := Diagonalize(s);
493: vol : integer := 1;
494:
495: begin
496: for k in m'range(1) loop
497: vol := vol*m(k,k);
498: end loop;
499: if vol >= 0
500: then return vol;
501: else return -vol;
502: end if;
503: end Volume;
504:
505: -- DESTRUCTORS :
506:
507: procedure Destroy_Neighbor ( s : in out Simplex; k : in natural ) is
508: begin
509: s.pts(k).si := Null_Simplex;
510: end Destroy_Neighbor;
511:
512: procedure Destroy_Neighbors ( s : in out Simplex ) is
513: begin
514: for k in s.pts'range loop
515: Destroy_Neighbor(s,k);
516: end loop;
517: end Destroy_Neighbors;
518:
519: procedure Clear_Neighbor ( s : in out Simplex; k : in natural ) is
520: begin
521: Clear(s.pts(k).si);
522: end Clear_Neighbor;
523:
524: procedure Clear_Neighbors ( s : in out Simplex ) is
525: begin
526: for k in s.pts'range loop
527: Clear_Neighbor(s,k);
528: end loop;
529: end Clear_Neighbors;
530:
531: procedure Clear ( s : in out Simplex ) is
532:
533: procedure free is new unchecked_deallocation(Simplex_Rep,Simplex);
534:
535: begin
536: Clear(s.tra);
537: free(s);
538: end Clear;
539:
540: end Simplices;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to simplices.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 / Dynlift