Annotation of OpenXM_contrib/PHC/Ada/Schubert/osculating_planes.adb, Revision 1.1
1.1 ! maekawa 1: with Chebychev_Polynomials; use Chebychev_Polynomials;
! 2: with Standard_Natural_Vectors;
! 3: with Standard_Floating_Vectors; use Standard_Floating_Vectors;
! 4: with Standard_Random_Numbers; use Standard_Random_Numbers;
! 5: with Standard_Floating_Linear_Solvers; use Standard_Floating_Linear_Solvers;
! 6: with Standard_Floating_QR_Decomposition; use Standard_Floating_QR_Decomposition;
! 7:
! 8: package body Osculating_Planes is
! 9:
! 10: function Standard_Basis
! 11: ( n,d : natural; s : double_float ) return Matrix is
! 12:
! 13: res : Matrix(1..n,1..d);
! 14: acc : double_float;
! 15: fac,lim : integer;
! 16:
! 17: begin
! 18: for i in 1..d loop -- set the zeros and ones
! 19: res(i,i) := 1.0;
! 20: for j in (i+1)..d loop
! 21: res(i,j) := 0.0;
! 22: end loop;
! 23: end loop;
! 24: for j in 1..d loop -- set the powers of the s-values
! 25: acc := s;
! 26: for i in (j+1)..n loop
! 27: res(i,j) := acc;
! 28: acc := acc*s;
! 29: end loop;
! 30: end loop;
! 31: for i in 3..n loop -- compute the factors from derivation
! 32: fac := 1;
! 33: if i-1 > d
! 34: then lim := d;
! 35: else lim := i-1;
! 36: end if;
! 37: for j in 2..lim loop
! 38: fac := fac*(i+1-j);
! 39: res(i,j) := double_float(fac)*res(i,j);
! 40: end loop;
! 41: if i <= d
! 42: then res(i,i) := double_float(fac);
! 43: end if;
! 44: end loop;
! 45: for j in 3..d loop -- scale the columns
! 46: for i in (j+1)..n loop
! 47: res(i,j) := res(i,j)/res(j,j);
! 48: end loop;
! 49: res(j,j) := 1.0;
! 50: end loop;
! 51: return res;
! 52: end Standard_Basis;
! 53:
! 54: function Chebychev_Basis
! 55: ( n,d : natural; s : double_float ) return Matrix is
! 56:
! 57: res : Matrix(1..n,1..d);
! 58: lim : integer;
! 59:
! 60: begin
! 61: for i in 1..d loop -- set the zeros and ones
! 62: res(i,i) := 1.0;
! 63: for j in (i+1)..d loop
! 64: res(i,j) := 0.0;
! 65: end loop;
! 66: end loop;
! 67: for i in 2..n loop -- evaluate and differentiate
! 68: declare
! 69: p : constant Vector := Chebychev_Polynomials.Create(i-1);
! 70: begin
! 71: res(i,1) := Eval(p,s);
! 72: if i > d
! 73: then lim := d;
! 74: else lim := i;
! 75: end if;
! 76: for j in 2..lim loop
! 77: declare
! 78: dp : constant Vector := Chebychev_Polynomials.Diff(p,j-1);
! 79: begin
! 80: res(i,j) := Eval(dp,s);
! 81: end;
! 82: end loop;
! 83: end;
! 84: end loop;
! 85: for j in 3..d loop -- scale the columns
! 86: for i in (j+1)..n loop
! 87: res(i,j) := res(i,j)/res(j,j);
! 88: end loop;
! 89: res(j,j) := 1.0;
! 90: end loop;
! 91: return res;
! 92: end Chebychev_Basis;
! 93:
! 94: function Orthogonal_Basis
! 95: ( n,d : natural; s : double_float ) return Matrix is
! 96:
! 97: res : Matrix(1..n,1..d) := Chebychev_Basis(n,d,s);
! 98: wrk : Matrix(1..n,1..d);
! 99: bas : Matrix(1..n,1..n);
! 100: qra : Standard_Floating_Vectors.Vector(1..n) := (1..n => 0.0);
! 101: pvt : Standard_Natural_Vectors.Vector(1..n) := (1..n => 0);
! 102:
! 103: begin
! 104: wrk := res;
! 105: QRD(wrk,qra,pvt,false);
! 106: for i in wrk'range(1) loop
! 107: for j in wrk'range(2) loop
! 108: bas(i,j) := wrk(i,j);
! 109: end loop;
! 110: for j in d+1..n loop
! 111: bas(i,j) := 0.0;
! 112: end loop;
! 113: end loop;
! 114: Basis(bas,res);
! 115: for i in res'range(1) loop
! 116: for j in res'range(2) loop
! 117: res(i,j) := bas(i,j);
! 118: end loop;
! 119: end loop;
! 120: return res;
! 121: end Orthogonal_Basis;
! 122:
! 123: function Determinant
! 124: ( mat : Matrix; rows : Standard_Natural_Vectors.Vector )
! 125: return double_float is
! 126:
! 127:
! 128: -- DESCRIPTION :
! 129: -- Computes the determinant of the matrix obtained by selecting rows.
! 130:
! 131: res : double_float := 1.0;
! 132: sqm : Matrix(rows'range,rows'range);
! 133: piv : Standard_Natural_Vectors.Vector(rows'range);
! 134: inf : natural;
! 135:
! 136: begin
! 137: for i in rows'range loop
! 138: piv(i) := i;
! 139: for j in rows'range loop
! 140: sqm(i,j) := mat(rows(i),j);
! 141: end loop;
! 142: end loop;
! 143: lufac(sqm,rows'last,piv,inf);
! 144: for i in rows'range loop
! 145: res := res*sqm(i,i);
! 146: end loop;
! 147: for i in piv'range loop
! 148: if piv(i) > i
! 149: then res := -res;
! 150: end if;
! 151: end loop;
! 152: return res;
! 153: end Determinant;
! 154:
! 155: procedure Maximal_Minors ( n,d : in natural; mat : in Matrix;
! 156: min,max : out double_float ) is
! 157:
! 158: -- DESCRIPTION :
! 159: -- Computes all maximal minors of a (nxd)-matrix mat, d < n.
! 160:
! 161: rows : Standard_Natural_Vectors.Vector(1..d);
! 162: first : boolean := true;
! 163: mindet,maxdet : double_float;
! 164:
! 165: procedure Select_Rows ( k,start : in natural ) is
! 166:
! 167: det : double_float;
! 168:
! 169: begin
! 170: if k > d
! 171: then det := Determinant(mat,rows);
! 172: det := abs(det);
! 173: if first
! 174: then mindet := det; maxdet := det; first := false;
! 175: else if det > maxdet
! 176: then maxdet := det;
! 177: elsif det < mindet
! 178: then mindet := det;
! 179: end if;
! 180: end if;
! 181: else for j in start..n loop
! 182: rows(k) := j;
! 183: Select_Rows(k+1,j+1);
! 184: end loop;
! 185: end if;
! 186: end Select_Rows;
! 187:
! 188: begin
! 189: Select_Rows(1,1);
! 190: min := mindet; max := maxdet;
! 191: end Maximal_Minors;
! 192:
! 193: procedure Sampled_Chebychev_Basis
! 194: ( n,d,m : in natural; mat : out Matrix;
! 195: s,ratio : out double_float ) is
! 196:
! 197: rs,min,max,rat,bestrat,bests : double_float;
! 198: first : boolean;
! 199: themat : Matrix(1..n,1..d);
! 200:
! 201: begin
! 202: first := true;
! 203: for i in 1..m loop
! 204: rs := Random;
! 205: themat := Chebychev_Basis(n,d,rs);
! 206: Maximal_Minors(n,d,themat,min,max);
! 207: rat := max/min;
! 208: if first
! 209: then bestrat := rat; bests := rs; first := false;
! 210: else if rat < bestrat
! 211: then bestrat := rat; bests := rs;
! 212: mat := themat;
! 213: end if;
! 214: end if;
! 215: end loop;
! 216: ratio := bestrat;
! 217: s := bests;
! 218: end Sampled_Chebychev_Basis;
! 219:
! 220: end Osculating_Planes;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to osculating_planes.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