Annotation of OpenXM_contrib/PHC/Ada/Schubert/driver_for_sagbi_homotopies.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io; use integer_io;
2: with Timing_Package; use Timing_Package;
3: with Communications_with_User; use Communications_with_User;
4: with Numbers_io; use Numbers_io;
5: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
6: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
7: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
8: with Standard_Random_Numbers; use Standard_Random_Numbers;
9: with Standard_Integer_Vectors;
10: with Standard_Natural_Matrices;
11: with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
12: with Standard_Complex_Vectors;
13: with Standard_Complex_VecVecs;
14: with Standard_Random_Vectors; use Standard_Random_Vectors;
15: with Standard_Floating_Vectors;
16: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
17: with Standard_Floating_Matrices;
18: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
19: with Standard_Complex_Matrices;
20: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
21: with Standard_Random_Matrices; use Standard_Random_Matrices;
22: with Symbol_Table; use Symbol_Table;
23: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
24: with Standard_Complex_Polynomials_io; use Standard_Complex_Polynomials_io;
25: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
26: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
27: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
28: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
29: with Standard_Complex_Jaco_Matrices;
30: with Standard_Complex_Laur_Systems; use Standard_Complex_Laur_Systems;
31: with Standard_Complex_Laur_Functions; use Standard_Complex_Laur_Functions;
32: with Standard_Complex_Laur_Systems; use Standard_Complex_Laur_Systems;
33: with Standard_Complex_Laur_SysFun; use Standard_Complex_Laur_SysFun;
34: with Standard_Complex_Laur_Jacomats;
35: with Exponent_Vectors; use Exponent_Vectors;
36: with Standard_Poly_Laur_Convertors; use Standard_Poly_Laur_Convertors;
37: with Matrix_Indeterminates;
38: with Homotopy;
39: with Standard_Complex_Solutions; use Standard_Complex_Solutions;
40: with Standard_Complex_Solutions_io; use Standard_Complex_Solutions_io;
41: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
42: with Lists_of_Integer_Vectors_io; use Lists_of_Integer_Vectors_io;
43: with Arrays_of_Integer_Vector_Lists; use Arrays_of_Integer_Vector_Lists;
44: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
45: with Standard_Root_Refiners; use Standard_Root_Refiners;
46: with Drivers_for_Poly_Continuation; use Drivers_for_Poly_Continuation;
47: with Power_Lists; use Power_Lists;
48: with Integer_Lifting_Utilities; use Integer_Lifting_Utilities;
49: with Integer_Mixed_Subdivisions; use Integer_Mixed_Subdivisions;
50: with Integer_Polyhedral_Continuation; use Integer_Polyhedral_Continuation;
51: with Triangulations; use Triangulations;
52: with Dynamic_Triangulations; use Dynamic_Triangulations;
53: with Triangulations_and_Subdivisions; use Triangulations_and_Subdivisions;
54: with Bracket_Expansions; use Bracket_Expansions;
55: with SAGBI_Homotopies; use SAGBI_Homotopies;
56: with Matrix_Homotopies;
57: with Matrix_Homotopies_io;
58: with Osculating_Planes; use Osculating_Planes;
59:
60: procedure Driver_for_SAGBI_Homotopies
61: ( file : in file_type; n,d : in natural ) is
62:
63: dim : constant natural := (n-d)*d;
64:
65: function Convert ( mat : Standard_Floating_Matrices.Matrix )
66: return Standard_Complex_Matrices.Matrix is
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_Osculating_SAGBI_Homotopy
80: ( file : file_type; locmap : Standard_Natural_Matrices.Matrix )
81: return Poly_Sys is
82:
83: -- DESCRIPTION :
84: -- Generates planes that are osculating a rational normal curve.
85: -- The system on return is the SAGBI homotopy.
86:
87: p : Poly;
88: l : Poly_Sys(1..dim);
89: s : double_float;
90: inc : double_float := 2.0/(double_float(dim));
91: mat : Standard_Floating_Matrices.Matrix(1..n,1..n-d);
92:
93: begin
94: -- p := Lifted_Localized_Laplace_Expansion(n,d);
95: p := Lifted_Localized_Laplace_Expansion(locmap);
96: new_line(file);
97: put_line(file,"The selected s-values : ");
98: s := Random;
99: for i in l'range loop
100: s := s + inc;
101: if s >= 1.0
102: then s := s - 2.0;
103: end if; -- s lies in [-1.0,+1.0]
104: put(file,s); new_line(file);
105: mat := Orthogonal_Basis(n,n-d,s);
106: Matrix_Homotopies.Add_Target(i,Convert(mat));
107: l(i) := Intersection_Condition(mat,p);
108: end loop;
109: return l;
110: end Random_Osculating_SAGBI_Homotopy;
111:
112: function Given_Osculating_SAGBI_Homotopy
113: ( file : file_type; locmap : Standard_Natural_Matrices.Matrix )
114: return Poly_Sys is
115:
116: -- DESCRIPTION :
117: -- Reads s-values and generates planes that are osculating a
118: -- rational normal curve.
119: -- The system on return is the SAGBI homotopy.
120:
121: p : Poly;
122: l : Poly_Sys(1..dim);
123: s : Standard_Floating_Vectors.Vector(1..dim);
124: mat : Standard_Floating_Matrices.Matrix(1..n,1..n-d);
125:
126: begin
127: -- p := Lifted_Localized_Laplace_Expansion(n,d);
128: p := Lifted_Localized_Laplace_Expansion(locmap);
129: new_line;
130: put("Give "); put(dim,1); put_line(" distinct real values for s : ");
131: for i in s'range loop
132: Read_Double_Float(s(i));
133: end loop;
134: new_line(file);
135: put_line(file,"The selected s-values : ");
136: put_line(file,s);
137: for i in l'range loop
138: mat := Orthogonal_Basis(n,n-d,s(i));
139: Matrix_Homotopies.Add_Target(i,Convert(mat));
140: l(i) := Intersection_Condition(mat,p);
141: end loop;
142: return l;
143: end Given_Osculating_SAGBI_Homotopy;
144:
145: function Read_Input_SAGBI_Homotopy
146: ( file : file_type; locmap : Standard_Natural_Matrices.Matrix )
147: return Poly_Sys is
148:
149: -- DESCRIPTION :
150: -- Reads the input planes from file.
151: -- The system on return is the SAGBI homotopy.
152:
153: planesfile : file_type;
154: p : Poly;
155: l : Poly_Sys(1..dim);
156: mat : Standard_Floating_Matrices.Matrix(1..n,1..n-d);
157:
158: begin
159: -- p := Lifted_Localized_Laplace_Expansion(n,d);
160: p := Lifted_Localized_Laplace_Expansion(locmap);
161: new_line;
162: put_line("Reading the name of the file with the input planes.");
163: Read_Name_and_Open_File(planesfile);
164: new_line(file);
165: put_line(file,"The input planes : ");
166: for i in l'range loop
167: get(planesfile,mat);
168: new_line(file); put(file,mat);
169: mat := Orthogonalize(mat);
170: Matrix_Homotopies.Add_Target(i,Convert(mat));
171: l(i) := Intersection_Condition(mat,p);
172: end loop;
173: Close(planesfile);
174: return l;
175: end Read_Input_SAGBI_Homotopy;
176:
177: function Complex_Random_SAGBI_Homotopy
178: ( locmap : Standard_Natural_Matrices.Matrix ) return Poly_Sys is
179:
180: -- DESCRIPTION :
181: -- Generates a SAGBI homotopy by taking random complex (n-d)-planes.
182: -- This system will be the start system in the Cheater's homotopy.
183:
184: p : Poly;
185: l : Poly_Sys(1..dim);
186: mat : Standard_Complex_Matrices.Matrix(1..n,1..n-d);
187:
188: begin
189: -- p := Lifted_Localized_Laplace_Expansion(n,d);
190: p := Lifted_Localized_Laplace_Expansion(locmap);
191: for i in l'range loop
192: mat := Random_Orthogonal_Matrix(n,n-d);
193: Matrix_Homotopies.Add_Start(i,mat);
194: l(i) := Intersection_Condition(mat,p);
195: end loop;
196: return l;
197: end Complex_Random_SAGBI_Homotopy;
198:
199: function Start_System ( l : Poly_Sys ) return Poly_Sys is
200:
201: -- DESCRIPTION :
202: -- Returns the system l after evaluation for t = 0.
203: -- This will be the start system in the SAGBI homotopy, flat deformation.
204:
205: r : Poly_Sys(l'range);
206: m : constant natural := Number_of_Unknowns(l(l'first));
207:
208: begin
209: for i in l'range loop
210: r(i) := Eval(l(i),Create(0.0),m);
211: end loop;
212: return r;
213: end Start_System;
214:
215: function Target_System ( l : Poly_Sys ) return Poly_Sys is
216:
217: -- DESCRIPTION :
218: -- Returns the system l after evaluation for t = 1.
219: -- This will be the target system in the SAGBI homotopy, flat deformation.
220:
221: r : Poly_Sys(l'range);
222: m : constant natural := Number_of_Unknowns(l(l'first));
223:
224: begin
225: for i in l'range loop
226: r(i) := Eval(l(i),Create(1.0),m);
227: end loop;
228: return r;
229: end Target_System;
230:
231: procedure Polyhedral_Homotopy_Continuation
232: ( file : in file_type;
233: p : in Poly_Sys; sols : in out Solution_List;
234: t : in Triangulation; lifted : in List; rep : in boolean ) is
235:
236: -- DESCRIPTION :
237: -- This is the first continuation stage, the resolution of the start system
238: -- in the SAGBI homotopy by means of polyhedral homotopy continuation.
239:
240: use Standard_Complex_Laur_JacoMats;
241:
242: timer : Timing_Widget;
243: lifted_lq,lq : Laur_Sys(p'range);
244: mix : Standard_Integer_Vectors.Vector(1..1) := (1..1 => p'last);
245: lif : Array_of_Lists(1..1) := (1..1 => lifted);
246: h : Eval_Coeff_Laur_Sys(p'range);
247: c : Standard_Complex_VecVecs.VecVec(h'range);
248: e : Exponent_Vectors_Array(h'range);
249: j : Eval_Coeff_Jaco_Mat(h'range,h'first..h'last+1);
250: m : Mult_Factors(j'range(1),j'range(2));
251: mixsub : Mixed_Subdivision;
252:
253: begin
254: tstart(timer);
255: mixsub := Shallow_Create(p'last,t);
256: lq := Polynomial_to_Laurent_System(p);
257: lifted_lq := Perform_Lifting(p'last,mix,lif,lq);
258: h := Create(lq);
259: for i in c'range loop
260: declare
261: coeff_lq : constant Standard_Complex_Vectors.Vector := Coeff(lq(i));
262: begin
263: c(i) := new Standard_Complex_Vectors.Vector(coeff_lq'range);
264: for k in coeff_lq'range loop
265: c(i)(k) := coeff_lq(k);
266: end loop;
267: end;
268: end loop;
269: e := Create(lq);
270: Create(lq,j,m);
271: if rep
272: then Mixed_Solve(file,lifted_lq,lif,h,c,e,j,m,mix,mixsub,sols);
273: else Mixed_Solve(lifted_lq,lif,h,c,e,j,m,mix,mixsub,sols);
274: end if;
275: tstop(timer);
276: new_line(file);
277: print_times(file,timer,"Polyhedral Homotopy Continuation");
278: end Polyhedral_Homotopy_Continuation;
279:
280: procedure Solve_Start_System
281: ( file : in file_type;
282: p : in Poly_Sys; sols : in out Solution_List;
283: report : in boolean ) is
284:
285: -- DESCRIPTION :
286: -- Computes the volume of the support of p. This is the
287: -- combinatorial-geometric set up for the first continuation stage.
288:
289: timer : Timing_Widget;
290: support : List := Create(p(p'first));
291: lifted,lifted_last : List;
292: t : Triangulation;
293: vol : natural;
294:
295: begin
296: new_line(file);
297: put_line(file,"The support of the start system : ");
298: new_line(file);
299: put(file,support);
300: tstart(timer);
301: Dynamic_Lifting(support,false,false,0,lifted,lifted_last,t);
302: new_line(file);
303: put_line(file,"The lifted support : ");
304: new_line(file);
305: put(file,lifted);
306: vol := Volume(t);
307: tstop(timer);
308: new_line(file);
309: put(file,"The volume : "); put(file,vol,1); new_line(file);
310: new_line(file);
311: print_times(file,timer,"Triangulation and Volume Computation");
312: Polyhedral_Homotopy_Continuation(file,p,sols,t,lifted,report);
313: Clear(t);
314: end Solve_Start_System;
315:
316: procedure Flat_Deformation
317: ( file : in file_type; sh : in Poly_Sys;
318: sols : in out Solution_List; report : in boolean ) is
319:
320: -- DESCRIPTION :
321: -- Performs flat deformation as defined by the SAGBI homotopy.
322: -- This is the second stage in the continuation.
323:
324: use Standard_Complex_Jaco_Matrices;
325:
326: timer : Timing_Widget;
327: sh_eval : Eval_Poly_Sys(sh'range);
328: jac_mat : Jaco_Mat(sh'range,sh'first..sh'last+1);
329: eva_jac : Eval_Jaco_Mat(jac_mat'range(1),jac_mat'range(2));
330:
331: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
332: return Standard_Complex_Vectors.Vector is
333:
334: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
335:
336: begin
337: xt(x'range) := x;
338: xt(xt'last) := t;
339: return Eval(sh_eval,xt);
340: end Eval;
341:
342: function Diff ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
343: return Standard_Complex_Matrices.Matrix is
344:
345: res : Standard_Complex_Matrices.Matrix(x'range,x'range);
346: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
347:
348: begin
349: xt(x'range) := x;
350: xt(xt'last) := t;
351: for i in res'range(1) loop
352: for j in res'range(2) loop
353: res(i,j) := Eval(eva_jac(i,j),xt);
354: end loop;
355: end loop;
356: return res;
357: end Diff;
358:
359: function Diff ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
360: return Standard_Complex_Vectors.Vector is
361:
362: res : Standard_Complex_Vectors.Vector(x'range);
363: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
364:
365: begin
366: xt(x'range) := x;
367: xt(xt'last) := t;
368: for i in res'range loop
369: res(i) := Eval(eva_jac(i,eva_jac'last(2)),xt);
370: end loop;
371: return res;
372: end Diff;
373:
374: procedure Sil_Cont is new Silent_Continue(Max_Norm,Eval,Diff,Diff);
375: procedure Rep_Cont is new Reporting_Continue(Max_Norm,Eval,Diff,Diff);
376:
377: begin
378: tstart(timer);
379: sh_eval := Create(sh);
380: jac_mat := Create(sh);
381: eva_jac := Create(jac_mat);
382: Set_Continuation_Parameter(sols,Create(0.0));
383: if report
384: then Rep_Cont(file,sols,false,Create(1.0));
385: else Sil_Cont(sols,false,Create(1.0));
386: end if;
387: tstop(timer);
388: new_line(file); print_times(file,timer,"Flat Deformation");
389: Clear(sh_eval);
390: Clear(jac_mat);
391: Clear(eva_jac);
392: end Flat_Deformation;
393:
394: procedure Solve_Target_System
395: ( file : in file_type; start,target : in Poly_Sys;
396: sols : in out Solution_List; report : in boolean ) is
397:
398: -- DESCRIPTION :
399: -- This is the third and last continuation stage: Cheater's homotopy.
400: -- It is implemented in the space of polynomials.
401:
402: timer : Timing_Widget;
403: n : constant natural := target'last;
404: a : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
405: b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
406:
407: begin
408: tstart(timer);
409: Homotopy.Create(target,start,2,a,b,true); -- linear cheater
410: Set_Continuation_Parameter(sols,Create(0.0));
411: declare
412: procedure Sil_Cont is
413: new Silent_Continue(Max_Norm,
414: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
415: procedure Rep_Cont is
416: new Reporting_Continue(Max_Norm,
417: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
418: begin
419: if report
420: then Rep_Cont(file,sols,false,Create(1.0));
421: else Sil_Cont(sols,false,Create(1.0));
422: end if;
423: end;
424: tstop(timer);
425: new_line(file);
426: print_times(file,timer,"Cheater's homotopy to target system");
427: end Solve_Target_System;
428:
429: procedure Solve_Target_System
430: ( file : in file_type; symtarget : in Poly;
431: sols : in out Solution_List; report : in boolean ) is
432:
433: -- DESCRIPTION :
434: -- This is the third and last continuation stage: Cheater's homotopy.
435: -- It uses a coefficient-matrix homotopy and takes as input the target
436: -- system with coefficients as brackets.
437: -- This is far more expensive than the linear Cheater's homotopy.
438:
439: use Standard_Complex_Jaco_Matrices;
440:
441: timer : Timing_Widget;
442: h : Standard_Complex_Poly_Functions.Eval_Coeff_Poly := Create(symtarget);
443: homsys : Poly_Sys(1..dim);
444: jacmat : Eval_Coeff_Jaco_Mat(1..dim,1..dim);
445: mulfac : Mult_Factors(1..dim,1..dim);
446: brkcff : Standard_Complex_Vectors.Vector := Coeff(symtarget);
447: syscff : Standard_Complex_VecVecs.VecVec(1..dim);
448:
449: begin
450: tstart(timer);
451: Set_Continuation_Parameter(sols,Create(0.0));
452: for i in homsys'range loop
453: Copy(symtarget,homsys(i));
454: end loop;
455: Create(homsys,jacmat,mulfac);
456: for i in syscff'range loop
457: syscff(i) := new Standard_Complex_Vectors.Vector(brkcff'range);
458: end loop;
459: declare
460:
461: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
462: return Standard_Complex_Vectors.Vector is
463:
464: y : Standard_Complex_Vectors.Vector(x'range);
465: c : Standard_Complex_Vectors.Vector(brkcff'range);
466:
467: begin
468: for i in y'range loop
469: c := Intersection_Coefficients(Matrix_Homotopies.Eval(i,t),brkcff);
470: y(i) := Eval(h,c,x);
471: end loop;
472: return y;
473: end Eval;
474:
475: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
476: return Standard_Complex_Vectors.Vector is
477:
478: y : Standard_Complex_Vectors.Vector(x'range) := x;
479:
480: begin
481: return y;
482: end dHt;
483:
484: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
485: return Standard_Complex_Matrices.Matrix is
486: begin
487: for i in syscff'range loop
488: syscff(i).all
489: := Intersection_Coefficients(Matrix_Homotopies.Eval(i,t),brkcff);
490: end loop;
491: return Eval(jacmat,mulfac,syscff,x);
492: end dHx;
493:
494: procedure Sil_Cont is new Silent_Continue(Max_Norm,Eval,dHt,dHx);
495: procedure Rep_Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);
496:
497: begin
498: if report
499: then Rep_Cont(file,sols,false,Create(1.0));
500: else Sil_Cont(sols,false,Create(1.0));
501: end if;
502: end;
503: Standard_Complex_VecVecs.Clear(syscff);
504: tstop(timer);
505: new_line(file);
506: print_times(file,timer,"Cheater's homotopy to target system");
507: end Solve_Target_System;
508:
509: procedure Refine_Roots ( file : in file_type;
510: p : in Poly_Sys; sols : in out Solution_List ) is
511:
512: epsxa : constant double_float := 10.0**(-8);
513: epsfa : constant double_float := 10.0**(-8);
514: tolsing : constant double_float := 10.0**(-8);
515: max : constant natural := 3;
516: numit : natural := 0;
517:
518: begin
519: Reporting_Root_Refiner(file,p,sols,epsxa,epsfa,tolsing,numit,max,false);
520: end Refine_Roots;
521:
522: procedure Set_Parameters ( file : in file_type; report : out boolean ) is
523:
524: -- DESCRIPTION :
525: -- Interactive determination of the continuation and output parameters.
526:
527: oc : natural;
528:
529: begin
530: new_line;
531: Driver_for_Continuation_Parameters(file);
532: new_line;
533: Driver_for_Process_io(file,oc);
534: report := not (oc = 0);
535: new_line;
536: put_line("See the output file for results.");
537: new_line;
538: end Set_Parameters;
539:
540: function t_Symbol return Symbol is
541:
542: -- DESCRIPTION :
543: -- Returns the symbol to represent the continuation parameter t.
544:
545: res : Symbol;
546:
547: begin
548: res(1) := 't';
549: for i in 2..res'last loop
550: res(i) := ' ';
551: end loop;
552: return res;
553: end t_Symbol;
554:
555: function Lowest_Degree_Localization_Map
556: ( n,d : natural ) return Standard_Natural_Matrices.Matrix is
557:
558: -- DESCRIPTION :
559: -- Returns the lowest-degree localization map.
560:
561: res : Standard_Natural_Matrices.Matrix(1..n,1..d);
562:
563: begin
564: for i in 1..n-d loop
565: for j in 1..d loop
566: res(i,j) := 2;
567: end loop;
568: end loop;
569: for i in n-d+1..n loop
570: for j in 1..d loop
571: if i = j+n-d
572: then res(i,j) := 1;
573: else res(i,j) := 0;
574: end if;
575: end loop;
576: end loop;
577: return res;
578: end Lowest_Degree_Localization_Map;
579:
580: function Read_Localization_Map
581: ( n,d : natural ) return Standard_Natural_Matrices.Matrix is
582:
583: -- DESCRIPTION :
584: -- Reads a localization map from standard input.
585:
586: res : Standard_Natural_Matrices.Matrix(1..n,1..d);
587:
588: begin
589: put_line("Reading localization map: 0,1 for I and 2 = free element.");
590: put_line("The lowest-degree localization map looks like :");
591: put(Lowest_Degree_Localization_Map(n,d));
592: put_line("and the general localization map is :");
593: put(Localization_Map(n,d));
594: put("Give a "); put(n,1); put("-by-"); put(d,1);
595: put_line(" matrix to represent your own map :");
596: get(res); skip_line;
597: return res;
598: end Read_Localization_Map;
599:
600: procedure Main is
601:
602: -- DESCRIPTION :
603: -- There are four parts in the elaboration of the SAGBI homotopy :
604: -- 0. Set up of the polynomials in the SAGBI homotopy;
605: -- 1. Solve the start system by polyhedral continuation;
606: -- 2. Apply the flat deformations to solve complex instance;
607: -- 3. Cheater's homotopy from complex to real instance.
608:
609: timer,totaltimer : Timing_Widget;
610: realsagbih,compsagbih,realtarget,comptarget,starts : Poly_Sys(1..dim);
611: sols : Solution_List;
612: report : boolean;
613: inputchoice,localchoice : character;
614: locmap : Standard_Natural_Matrices.Matrix(1..n,1..d);
615:
616: begin
617: new_line;
618: put_line("MENU for a localization for the output planes : ");
619: put_line(" 1. Lowest-degree map, with I in lower-right corner.");
620: put_line(" 2. General map, able to represent any intersection.");
621: put_line(" 3. Give in your own localization map.");
622: put("Type 1, 2, or 3 to select : "); Ask_Alternative(localchoice,"123");
623: case localchoice is
624: when '1' => locmap := Lowest_Degree_Localization_Map(n,d);
625: when '2' => locmap := Localization_Map(n,d);
626: when '3' => locmap := Read_Localization_Map(n,d);
627: when others => null;
628: end case;
629: new_line;
630: put_line("MENU for constructing target system based on input planes : ");
631: put_line(" 1. Generate input planes osculating a rational normal curve.");
632: put_line(" 2. Interactively give s-values for the "
633: & "osculating input planes.");
634: put_line(" 3. Give the name of the file with input planes.");
635: put("Type 1, 2, or 3 to select : "); Ask_Alternative(inputchoice,"123");
636: tstart(totaltimer);
637: tstart(timer);
638: Matrix_Indeterminates.Initialize_Symbols(n,d);
639: Matrix_Indeterminates.Reduce_Symbols(locmap);
640: Matrix_Homotopies.Init(dim);
641: case inputchoice is
642: when '1' => realsagbih := Random_Osculating_SAGBI_Homotopy(file,locmap);
643: when '2' => realsagbih := Given_Osculating_SAGBI_Homotopy(file,locmap);
644: when '3' => realsagbih := Read_Input_SAGBI_Homotopy(file,locmap);
645: when others => null;
646: end case;
647: realtarget := Target_System(realsagbih);
648: compsagbih := Complex_Random_SAGBI_Homotopy(locmap);
649: comptarget := Target_System(compsagbih);
650: starts := Start_System(compsagbih);
651: Symbol_Table.Add(t_Symbol);
652: new_line(file);
653: put_line(file,"The target polynomial system with real coefficients : ");
654: new_line(file);
655: put(file,realtarget'last,realtarget);
656: new_line(file);
657: put_line(file,"The localization map : "); put(file,locmap);
658: new_line(file);
659: Matrix_Homotopies_io.Write(file);
660: put_line(file,"The target polynomial system with complex coefficients : ");
661: put_line(file,comptarget);
662: new_line(file);
663: put_line(file,"The SAGBI Homotopy as complex polynomial system : ");
664: new_line(file);
665: put_line(file,compsagbih);
666: new_line(file);
667: put_line(file,"The SAGBI Homotopy at t=0, the start system : ");
668: new_line(file);
669: put_line(file,starts);
670: tstop(timer);
671: new_line(file);
672: print_times(file,timer,"Setting up the polynomials in the SAGBI homotopy");
673: Set_Parameters(file,report);
674: Solve_Start_System(file,starts,sols,report);
675: Flat_Deformation(file,compsagbih,sols,report);
676: Solve_Target_System(file,comptarget,realtarget,sols,report);
677: -- Solve_Target_System(file,symtarget,sols,report);
678: -- this is the determinantal cheater's homotopy, very expensive
679: Refine_Roots(file,realtarget,sols);
680: Matrix_Indeterminates.Clear_Symbols;
681: tstop(totaltimer);
682: new_line(file);
683: print_times(file,totaltimer,"Total time for Solving the System");
684: end Main;
685:
686: begin
687: Main;
688: end Driver_for_SAGBI_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to driver_for_sagbi_homotopies.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