Annotation of OpenXM_contrib/PHC/Ada/Schubert/drivers_for_input_planes.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io,Numbers_io; use integer_io,Numbers_io;
2: with Communications_with_User; use Communications_with_User;
3: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
4: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
5: with Standard_Random_Numbers; use Standard_Random_Numbers;
6: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
7: with Standard_Floating_Matrices;
8: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
9: with Standard_Complex_Matrices;
10: with Standard_Random_Matrices; use Standard_Random_Matrices;
11: with Standard_Complex_VecMats_io; use Standard_Complex_VecMats_io;
12: with Osculating_Planes; use Osculating_Planes;
13:
14: package body Drivers_for_Input_Planes is
15:
16: function Finite ( dim : Bracket; fin_sum : natural ) return boolean is
17:
18: -- DESCRIPTION :
19: -- Returns true if the co-dimensions in dim will lead to a finite
20: -- number of solutions.
21:
22: sum : natural := 0;
23:
24: begin
25: for i in dim'range loop
26: sum := sum + dim(i);
27: end loop;
28: if sum = fin_sum
29: then return true;
30: else return false;
31: end if;
32: end Finite;
33:
34: function Read_Codimensions ( m,p,q : natural ) return Bracket is
35:
36: mpq : constant natural := m*p + q*(m+p);
37: codim : Bracket(1..mpq);
38: n : natural;
39:
40: begin
41: loop
42: put("Give number of intersection conditions : "); get(n);
43: put("Give "); put(n,1); put(" co-dimensions (sum = ");
44: put(mpq,1); put(") : ");
45: for i in 1..n loop
46: get(codim(i));
47: end loop;
48: for i in 1..n-1 loop
49: put(codim(i),1); put(" + ");
50: end loop;
51: put(codim(n),1);
52: if Finite(codim(1..n),mpq)
53: then put(" = "); put(mpq,1); put_line(" Finite #sols.");
54: exit;
55: else put(" /= "); put(mpq,1);
56: put_line(" Please try again.");
57: end if;
58: end loop;
59: return codim(1..n);
60: end Read_Codimensions;
61:
62: function Convert ( mat : Standard_Floating_Matrices.Matrix )
63: return Standard_Complex_Matrices.Matrix is
64:
65: -- DESCRIPTION :
66: -- Converts a real matrix into a complex one.
67:
68: res : Standard_Complex_Matrices.Matrix(mat'range(1),mat'range(2));
69:
70: begin
71: for i in mat'range(1) loop
72: for j in mat'range(2) loop
73: res(i,j) := Create(mat(i,j));
74: end loop;
75: end loop;
76: return res;
77: end Convert;
78:
79: function Random_Complex_Planes ( m,p : natural ) return VecMat is
80:
81: res : VecMat(1..m*p);
82: n : constant natural := m+p;
83: use Standard_Complex_Matrices;
84:
85: begin
86: for i in res'range loop
87: res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
88: end loop;
89: return res;
90: end Random_Complex_Planes;
91:
92: function Random_Complex_Planes ( m,p : natural; k : Bracket ) return VecMat is
93:
94: res : VecMat(k'range);
95: n : constant natural := m+p;
96: use Standard_Complex_Matrices;
97:
98: begin
99: for i in res'range loop
100: res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m+1-k(i)));
101: end loop;
102: return res;
103: end Random_Complex_Planes;
104:
105: function Random_Complex_Planes ( m,p,q : natural ) return VecMat is
106:
107: res : VecMat(1..(m*p+q*(m+p)));
108: n : constant natural := m+p;
109: use Standard_Complex_Matrices;
110:
111: begin
112: for i in res'range loop
113: res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
114: end loop;
115: return res;
116: end Random_Complex_Planes;
117:
118: function Random_Real_Planes ( m,p : natural ) return VecMat is
119:
120: res : VecMat(1..m*p);
121: n : constant natural := m+p;
122: fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
123: cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);
124:
125: begin
126: for i in res'range loop
127: fltmat := Random_Orthogonal_Matrix(n,m);
128: cmpmat := Convert(fltmat);
129: res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
130: end loop;
131: return res;
132: end Random_Real_Planes;
133:
134: function Random_Real_Planes ( m,p : natural; k : Bracket ) return VecMat is
135:
136: res : VecMat(k'range);
137: n : constant natural := m+p;
138:
139: begin
140: for i in res'range loop
141: declare
142: fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
143: cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m+1-k(i));
144: begin
145: fltmat := Random_Orthogonal_Matrix(n,m+1-k(i));
146: cmpmat := Convert(fltmat);
147: res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
148: end;
149: end loop;
150: return res;
151: end Random_Real_Planes;
152:
153: function Random_Real_Planes ( m,p,q : natural ) return VecMat is
154:
155: res : VecMat(1..(m*p+q*(m+p)));
156: n : constant natural := m+p;
157: fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
158: cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);
159:
160: begin
161: for i in res'range loop
162: fltmat := Random_Orthogonal_Matrix(n,m);
163: cmpmat := Convert(fltmat);
164: res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
165: end loop;
166: return res;
167: end Random_Real_Planes;
168:
169: function Equidistant_Interpolation_Points ( n : natural ) return Vector is
170:
171: res : Vector(1..n);
172: s : double_float := Random;
173: inc : constant double_float := 2.0/double_float(n);
174:
175: begin
176: for i in res'range loop
177: res(i) := s;
178: s := s+inc;
179: if s >= 1.0
180: then s := s - 2.0;
181: end if;
182: end loop;
183: return res;
184: end Equidistant_Interpolation_Points;
185:
186: function Read_Interpolation_Points ( n : natural ) return Vector is
187:
188: res : Vector(1..n);
189:
190: begin
191: new_line;
192: put("Give "); put(n,1); put_line(" distinct real interpolation points : ");
193: for i in res'range loop
194: Read_Double_Float(res(i));
195: end loop;
196: return res;
197: end Read_Interpolation_Points;
198:
199: function Osculating_Input_Planes ( m,p : natural ) return VecMat is
200:
201: s : constant Vector := Equidistant_Interpolation_Points(m*p);
202:
203: begin
204: return Osculating_Input_Planes(m,p,s);
205: end Osculating_Input_Planes;
206:
207: function Osculating_Input_Planes
208: ( m,p : natural; s : Vector ) return VecMat is
209:
210: res : VecMat(1..m*p);
211: n : constant natural := m+p;
212: realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
213:
214: begin
215: for i in res'range loop
216: realmat := Orthogonal_Basis(n,m,s(i));
217: res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
218: end loop;
219: return res;
220: end Osculating_Input_Planes;
221:
222: function Osculating_Input_Planes
223: ( m,p : natural; k : bracket ) return VecMat is
224:
225: s : constant Vector := Equidistant_Interpolation_Points(k'length);
226:
227: begin
228: return Osculating_Input_Planes(m,p,k,s);
229: end Osculating_Input_Planes;
230:
231: function Osculating_Input_Planes
232: ( m,p : natural; k : bracket; s : Vector ) return VecMat is
233:
234: res : VecMat(k'range);
235: n : constant natural := m+p;
236:
237: begin
238: for i in res'range loop
239: declare
240: realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
241: begin
242: realmat := Orthogonal_Basis(n,m+1-k(i),s(i));
243: res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
244: end;
245: end loop;
246: return res;
247: end Osculating_Input_Planes;
248:
249: function Osculating_Input_Planes ( m,p,q : natural ) return VecMat is
250:
251: dim : constant natural := m*p + q*(m+p);
252: s : constant Vector := Equidistant_Interpolation_Points(dim);
253:
254: begin
255: return Osculating_Input_Planes(m,p,q,s);
256: end Osculating_Input_Planes;
257:
258: function Osculating_Input_Planes
259: ( m,p,q : natural; s : Vector ) return VecMat is
260:
261: res : VecMat(1..m*p+q*(m+p));
262: n : constant natural := m+p;
263: realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
264:
265: begin
266: for i in res'range loop
267: realmat := Orthogonal_Basis(n,m,s(i));
268: res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
269: end loop;
270: return res;
271: end Osculating_Input_Planes;
272:
273: function Read_Input_Planes ( m,p : natural ) return VecMat is
274:
275: res : VecMat(1..m*p);
276: planesfile : file_type;
277: n : constant natural := m+p;
278: realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
279:
280: begin
281: new_line;
282: put_line("Reading the name of the file with the input planes.");
283: Read_Name_and_Open_File(planesfile);
284: for i in res'range loop
285: get(planesfile,realmat);
286: realmat := Orthogonalize(realmat);
287: res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
288: end loop;
289: Close(planesfile);
290: return res;
291: end Read_Input_Planes;
292:
293: function Read_Input_Planes ( m,p : natural; k : in Bracket ) return VecMat is
294:
295: res : VecMat(k'range);
296: planesfile : file_type;
297: n : constant natural := m+p;
298:
299: begin
300: new_line;
301: put_line("Reading the name of the file with the input planes.");
302: Read_Name_and_Open_File(planesfile);
303: for i in res'range loop
304: declare
305: realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
306: begin
307: get(planesfile,realmat);
308: realmat := Orthogonalize(realmat);
309: res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
310: end;
311: end loop;
312: Close(planesfile);
313: return res;
314: end Read_Input_Planes;
315:
316: function Read_Input_Planes ( m,p,q : natural ) return VecMat is
317:
318: res : VecMat(1..m*p+q*(m+p));
319: planesfile : file_type;
320: n : constant natural := m+p;
321: realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
322:
323: begin
324: new_line;
325: put_line("Reading the name of the file with the input planes.");
326: Read_Name_and_Open_File(planesfile);
327: for i in res'range loop
328: get(planesfile,realmat);
329: realmat := Orthogonalize(realmat);
330: res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
331: end loop;
332: Close(planesfile);
333: return res;
334: end Read_Input_Planes;
335:
336: -- MAIN INTERACTIVE DRIVERS :
337:
338: function Select_Input_Choice return character is
339:
340: -- DESCRIPTION :
341: -- Displays the menu to obtain a choice for the kind of input.
342:
343: res : character;
344:
345: begin
346: new_line;
347: put_line("MENU for generating real input planes : ");
348: put_line(" 0. Random real planes at equidistant interpolation points.");
349: put_line(" 1. Generate input planes osculating a rational normal curve.");
350: put_line(" 2. Interactively give s-values for the "
351: & "osculating input planes.");
352: put_line(" 3. Give the name of the file with input planes.");
353: put("Type 0, 1, 2, or 3 to select : "); Ask_Alternative(res,"0123");
354: return res;
355: end Select_Input_Choice;
356:
357: procedure Driver_for_Input_Planes
358: ( file : in file_type; m,p : in natural; planes : out VecMat ) is
359:
360: input_choice : character := Select_Input_Choice;
361: dim : constant natural := m*p;
362: svals : Vector(1..dim);
363:
364: begin
365: case input_choice is
366: when '0' => planes := Random_Real_Planes(m,p);
367: when '1' => planes := Osculating_Input_Planes(m,p);
368: when '2' => svals := Read_Interpolation_Points(dim);
369: planes := Osculating_Input_Planes(m,p,svals);
370: when '3' => planes := Read_Input_Planes(m,p);
371: when others => null;
372: end case;
373: end Driver_for_Input_Planes;
374:
375: procedure Driver_for_Input_Planes
376: ( file : in file_type; m,p : in natural; k : in Bracket;
377: planes : out VecMat ) is
378:
379: input_choice : character := Select_Input_Choice;
380: svals : Vector(k'range);
381:
382: begin
383: case input_choice is
384: when '0' => planes := Random_Real_Planes(m,p,k);
385: when '1' => planes := Osculating_Input_Planes(m,p,k);
386: when '2' => svals := Read_Interpolation_Points(k'length);
387: planes := Osculating_Input_Planes(m,p,k,svals);
388: when '3' => planes := Read_Input_Planes(m,p,k);
389: when others => null;
390: end case;
391: end Driver_for_Input_Planes;
392:
393: procedure Driver_for_Input_Planes
394: ( file : in file_type; m,p,q : in natural;
395: s : out Vector; planes : out VecMat ) is
396:
397: input_choice : character := Select_Input_Choice;
398: dim : constant natural := m*p + q*(m+p);
399: svals : Vector(1..dim);
400:
401: begin
402: case input_choice is
403: when '0' => svals := Equidistant_Interpolation_Points(dim);
404: s := svals;
405: planes := Random_Real_Planes(m,p,q);
406: when '1' => svals := Equidistant_Interpolation_Points(dim);
407: s := svals;
408: planes := Osculating_Input_Planes(m,p,q,svals);
409: when '2' => svals := Read_Interpolation_Points(dim);
410: s := svals;
411: planes := Osculating_Input_Planes(m,p,q,svals);
412: when '3' => svals := Read_Interpolation_Points(dim);
413: s := svals;
414: planes := Read_Input_Planes(m,p,q);
415: when others => null;
416: end case;
417: put_line(file,"The interpolation points : ");
418: put_line(file,svals);
419: put_line(file,"The target planes : ");
420: put(file,planes,3);
421: end Driver_for_Input_Planes;
422:
423: end Drivers_for_Input_Planes;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>