Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_shapiro.adb, Revision 1.1.1.1
1.1 maekawa 1: with text_io,integer_io; use text_io,integer_io;
2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
3: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
4: with Standard_Random_Numbers; use Standard_Random_Numbers;
5: with Standard_Floating_Vectors; use Standard_Floating_Vectors;
6: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
7: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
8: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
9: with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
10: with Standard_Floating_Linear_Solvers; use Standard_Floating_Linear_Solvers;
11: with Osculating_Planes; use Osculating_Planes;
12:
13: procedure ts_shapiro is
14:
15: -- DESCRIPTION :
16: -- Tests the generation of special matrices to test Shapiro's conjecture.
17:
18: function Determinant
19: ( mat : Matrix; rows : Standard_Natural_Vectors.Vector )
20: return double_float is
21:
22: -- DESCRIPTION :
23: -- Computes the determinant of the matrix obtained by selecting rows.
24:
25: res : double_float := 1.0;
26: sqm : Matrix(rows'range,rows'range);
27: piv : Standard_Natural_Vectors.Vector(rows'range);
28: inf : natural;
29:
30: begin
31: for i in rows'range loop
32: piv(i) := i;
33: for j in rows'range loop
34: sqm(i,j) := mat(rows(i),j);
35: end loop;
36: end loop;
37: lufac(sqm,rows'last,piv,inf);
38: for i in rows'range loop
39: res := res*sqm(i,i);
40: end loop;
41: for i in piv'range loop
42: if piv(i) > i
43: then res := -res;
44: end if;
45: end loop;
46: return res;
47: end Determinant;
48:
49: procedure Maximal_Minors ( n,d : in natural; mat : in Matrix;
50: min,max : out double_float ) is
51:
52: -- DESCRIPTION :
53: -- Computes all maximal minors of a (nxd)-matrix mat, d < n.
54:
55: rows : Standard_Natural_Vectors.Vector(1..d);
56: first : boolean := true;
57: mindet,maxdet : double_float;
58:
59: procedure Select_Rows ( k,start : in natural ) is
60:
61: det : double_float;
62:
63: begin
64: if k > d
65: then det := Determinant(mat,rows);
66: -- put("Minor "); put(rows); put(" equals "); put(det); new_line;
67: det := abs(det);
68: if first
69: then mindet := det; maxdet := det; first := false;
70: else if det > maxdet
71: then maxdet := det;
72: elsif det < mindet
73: then mindet := det;
74: end if;
75: end if;
76: else for j in start..n loop
77: rows(k) := j;
78: Select_Rows(k+1,j+1);
79: end loop;
80: end if;
81: end Select_Rows;
82:
83: begin
84: Select_Rows(1,1);
85: put("Min : "); put(mindet,3,3,3);
86: put(" Max : "); put(maxdet,3,3,3);
87: put(" Max/Min : "); put(maxdet/mindet,3,3,3); new_line;
88: min := mindet; max := maxdet;
89: end Maximal_Minors;
90:
91: procedure Test_Sample ( n,d : natural; s : double_float ) is
92:
93: cheb_mat : Matrix(1..n,1..d) := Chebychev_Basis(n,d,s);
94: orto_mat : Matrix(1..n,1..d) := Orthogonal_Basis(n,d,s);
95: min,max : double_float;
96:
97: begin
98: put("Osculating "); put(d,1); put("-plane in ");
99: put(n,1); put_line("-space : "); put(cheb_mat);
100: put_line("All maximal minors : ");
101: Maximal_Minors(n,d,cheb_mat,min,max);
102: put("Orthogonal respresentation of osculating ");
103: put(d,1); put("-plane in "); put(n,1); put_line("-space : ");
104: put(orto_mat);
105: put_line("All maximal minors : ");
106: Maximal_Minors(n,d,orto_mat,min,max);
107: end Test_Sample;
108:
109: procedure Sampled_Generation is
110:
111: n,d : natural;
112: s : double_float;
113: ans : character;
114:
115: begin
116: new_line;
117: put("Give n : "); get(n);
118: put("Give d : "); get(d);
119: loop
120: s := Random;
121: put("The s-value : "); put(s); new_line;
122: Test_Sample(n,d,s);
123: put("Do you want to test another sample ? (y/n) "); get(ans);
124: exit when (ans /= 'y');
125: end loop;
126: end Sampled_Generation;
127:
128: procedure Best_Sampled_Generation is
129:
130: n,d,m : natural;
131: ans : character;
132: s,bestratio,bests : double_float;
133: first : boolean;
134:
135: begin
136: new_line;
137: put("Give n : "); get(n);
138: put("Give d : "); get(d);
139: loop
140: put("Give the number of samples : "); get(m);
141: first := true;
142: for i in 1..m loop
143: s := Random;
144: put("The s-value : "); put(s); new_line;
145: declare
146: mat : Matrix(1..n,1..d) := Chebychev_Basis(n,d,s);
147: min,max,ratio : double_float;
148: begin
149: Maximal_Minors(n,d,mat,min,max);
150: ratio := max/min;
151: if first
152: then bestratio := ratio; bests := s; first := false;
153: else if ratio < bestratio
154: then bestratio := ratio; bests := s;
155: end if;
156: end if;
157: end;
158: end loop;
159: put("Best ratio : "); put(bestratio);
160: put(" for s : "); put(bests); new_line;
161: put("Do you want more ratio's to test ? (y/n) "); get(ans);
162: exit when (ans /= 'y');
163: end loop;
164: end Best_Sampled_Generation;
165:
166: procedure Update ( s : in out double_float; inc : in double_float ) is
167: begin
168: s := s + inc;
169: if s >= 1.0
170: then s := s - 2.0;
171: end if;
172: end Update;
173:
174: procedure Equidistant_Sampling
175: ( n,d,nb : in natural; inits : in double_float;
176: rat : out double_float ) is
177:
178: -- DESCRIPTION :
179: -- Generates nb equidistant s-values and computes the maximum
180: -- of all ratios max/min minors.
181:
182: inc : constant double_float := 2.0/double_float(nb);
183: mat : Matrix(1..n,1..d);
184: s : double_float := inits;
185: first : boolean := true;
186: min,max,ratio,maxratio : double_float;
187:
188: begin
189: for i in 1..nb loop
190: mat := Chebychev_Basis(n,d,s);
191: Maximal_Minors(n,d,mat,min,max);
192: ratio := max/min;
193: if first
194: then maxratio := ratio; first := false;
195: else if ratio > maxratio
196: then maxratio := ratio;
197: end if;
198: end if;
199: Update(s,inc);
200: end loop;
201: rat := maxratio;
202: end Equidistant_Sampling;
203:
204: procedure Init ( mat : in out Matrix ) is
205: begin
206: for i in mat'range(1) loop
207: for j in mat'range(2) loop
208: mat(i,j) := 0.0;
209: end loop;
210: end loop;
211: end Init;
212:
213: procedure Div ( mat : in out Matrix; d : in double_float ) is
214: begin
215: for i in mat'range(1) loop
216: for j in mat'range(2) loop
217: mat(i,j) := mat(i,j)/d;
218: end loop;
219: end loop;
220: end Div;
221:
222: procedure Averaged_Equidistant_Sampling
223: ( n,d,nb,nav : in natural; inits : in double_float;
224: rat : out double_float ) is
225:
226: -- DESCRIPTION :
227: -- Generates nb equidistant s-values and computes the maximum
228: -- of all ratios max/min minors.
229: -- Averages over every interval.
230:
231: inc : constant double_float := 2.0/double_float(nb*nav);
232: mat : Matrix(1..n,1..d);
233: s : double_float := inits;
234: first : boolean := true;
235: min,max,ratio,maxratio : double_float;
236:
237: begin
238: for i in 1..nb loop
239: Init(mat);
240: for j in 1..nav loop
241: mat := mat + Chebychev_Basis(n,d,s);
242: Update(s,inc);
243: end loop;
244: Div(mat,double_float(nav));
245: Maximal_Minors(n,d,mat,min,max);
246: ratio := max/min;
247: if first
248: then maxratio := ratio; first := false;
249: else if ratio > maxratio
250: then maxratio := ratio;
251: end if;
252: end if;
253: Update(s,inc);
254: end loop;
255: rat := maxratio;
256: end Averaged_Equidistant_Sampling;
257:
258: procedure Best_Equidistant_Sampling is
259:
260: n,d,m,nb : natural;
261: ans : character;
262: s,ratio,bestratio,bests : double_float;
263: first : boolean;
264:
265: begin
266: new_line;
267: put("Give n : "); get(n);
268: put("Give d : "); get(d);
269: loop
270: put("Give #samples : "); get(m);
271: put("Give #equidistant points : "); get(nb);
272: first := true;
273: for i in 1..m loop
274: s := Random;
275: Equidistant_Sampling(n,d,nb,s,ratio);
276: if first
277: then bestratio := ratio; bests := s; first := false;
278: else if ratio < bestratio
279: then bestratio := ratio; bests := s;
280: end if;
281: end if;
282: end loop;
283: put("Best ratio : "); put(bestratio);
284: put(" for s : "); put(bests); new_line;
285: put("Do you want more ratio's to test ? (y/n) "); get(ans);
286: exit when (ans /= 'y');
287: end loop;
288: end Best_Equidistant_Sampling;
289:
290: procedure Best_Averaged_Equidistant_Sampling is
291:
292: n,d,m,nb,nav : natural;
293: ans : character;
294: s,ratio,bestratio,bests : double_float;
295: first : boolean;
296:
297: begin
298: new_line;
299: put("Give n : "); get(n);
300: put("Give d : "); get(d);
301: put("Give avg #times : "); get(nav);
302: loop
303: put("Give #samples : "); get(m);
304: put("Give #equidistant points : "); get(nb);
305: first := true;
306: for i in 1..m loop
307: s := Random;
308: Averaged_Equidistant_Sampling(n,d,nb,nav,s,ratio);
309: if first
310: then bestratio := ratio; bests := s; first := false;
311: else if ratio < bestratio
312: then bestratio := ratio; bests := s;
313: end if;
314: end if;
315: end loop;
316: put("Best ratio : "); put(bestratio);
317: put(" for s : "); put(bests); new_line;
318: put("Do you want more ratio's to test ? (y/n) "); get(ans);
319: exit when (ans /= 'y');
320: end loop;
321: end Best_Averaged_Equidistant_Sampling;
322:
323: function Distance ( v : Standard_Floating_Vectors.Vector;
324: k : natural ) return double_float is
325:
326: min : double_float := 2.0;
327:
328: begin
329: for i in v'first..(k-1) loop
330: if abs(v(i)-v(k)) < min
331: then min := abs(v(i)-v(k));
332: end if;
333: end loop;
334: return min;
335: end Distance;
336:
337: procedure Spaced_Sampling
338: ( n,d,nb : in natural; mindist : in double_float;
339: rat : out double_float ) is
340:
341: -- DESCRIPTION :
342: -- Generates nb distinct s-values, not closer to each other than
343: -- mindist and computes the maximum of all ratios max/min minors.
344:
345: mat : Matrix(1..n,1..d);
346: sva : Standard_Floating_Vectors.Vector(1..nb);
347: first : boolean := true;
348: min,max,ratio,maxratio : double_float;
349:
350: begin
351: for i in 1..nb loop
352: loop
353: sva(i) := Random;
354: exit when (Distance(sva,i) > mindist);
355: end loop;
356: mat := Chebychev_Basis(n,d,sva(i));
357: Maximal_Minors(n,d,mat,min,max);
358: ratio := max/min;
359: if first
360: then maxratio := ratio; first := false;
361: else if ratio > maxratio
362: then maxratio := ratio;
363: end if;
364: end if;
365: end loop;
366: rat := maxratio;
367: end Spaced_Sampling;
368:
369: procedure Best_Spaced_Sampling is
370:
371: n,d,m,nb : natural;
372: ans : character;
373: inc,mindist,ratio,bestratio : double_float;
374: first : boolean;
375:
376: begin
377: new_line;
378: put("Give n : "); get(n);
379: put("Give d : "); get(d);
380: loop
381: put("Give #samples : "); get(m);
382: put("Give #equidistant points : "); get(nb);
383: inc := 2.0/double_float(nb);
384: put("The increment : "); put(inc); new_line;
385: put("Give Minimal distance : "); get(mindist);
386: first := true;
387: for i in 1..m loop
388: Spaced_Sampling(n,d,nb,mindist,ratio);
389: if first
390: then bestratio := ratio; first := false;
391: else if ratio < bestratio
392: then bestratio := ratio;
393: end if;
394: end if;
395: end loop;
396: put("Best ratio : "); put(bestratio); new_line;
397: put("Do you want more ratio's to test ? (y/n) "); get(ans);
398: exit when (ans /= 'y');
399: end loop;
400: end Best_Spaced_Sampling;
401:
402: procedure Interactive_Generation is
403:
404: n,d : natural;
405: s : double_float;
406: ans : character;
407:
408: begin
409: new_line;
410: put("Give n : "); get(n);
411: put("Give d : "); get(d);
412: loop
413: put("Give s-value : "); get(s);
414: Test_Sample(n,d,s);
415: put("Do you want to test another s-value ? (y/n) "); get(ans);
416: exit when (ans /= 'y');
417: end loop;
418: end Interactive_Generation;
419:
420: procedure Main is
421:
422: ans : character;
423:
424: begin
425: new_line;
426: put_line("Generation of osculating d-planes in n-space.");
427: loop
428: new_line;
429: put_line("Choose one of the following : ");
430: put_line(" 0. Exit this program.");
431: put_line(" 1. Interactively input of s-values.");
432: put_line(" 2. Sampling s-values one after the other.");
433: put_line(" 3. Sample many times for best ratio.");
434: put_line(" 4. Equidistant Sampling s-values and select min ratio.");
435: put_line(" 5. Averaged Equidistant Sampling s-values + min ratio.");
436: put_line(" 6. Spaced Sampling s-values for min ratio.");
437: put("Make your choice (0,1,2,3,4,5 or 6) : "); get(ans);
438: exit when (ans = '0');
439: case ans is
440: when '1' => Interactive_Generation;
441: when '2' => Sampled_Generation;
442: when '3' => Best_Sampled_Generation;
443: when '4' => Best_Equidistant_Sampling;
444: when '5' => Best_Averaged_Equidistant_Sampling;
445: when '6' => Best_Spaced_Sampling;
446: when others => null;
447: end case;
448: end loop;
449: end Main;
450:
451: begin
452: Main;
453: end ts_shapiro;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ts_shapiro.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Schubert