Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_shapiro.adb, Revision 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