![[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|
Return to osculating_planes.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Schubert |
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;