Annotation of OpenXM_contrib/PHC/Ada/Schubert/driver_for_quantum_pieri.adb, Revision 1.1.1.1
1.1 maekawa 1: with text_io,integer_io; use text_io,integer_io;
2: with Communications_with_User; use Communications_with_User;
3: with Timing_Package; use Timing_Package;
4: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
5: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
6: with Standard_Floating_Vectors;
7: with Standard_Natural_Vectors;
8: with Standard_Complex_Vectors;
9: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
10: with Standard_Natural_Matrices;
11: with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
12: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
13: with Standard_Random_Vectors; use Standard_Random_Vectors;
14: with Standard_Complex_VecMats; use Standard_Complex_VecMats;
15: with Symbol_Table; use Symbol_Table;
16: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
17: with Standard_Complex_Poly_Matrices;
18: with Standard_Complex_Poly_Matrices_io; use Standard_Complex_Poly_Matrices_io;
19: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
20: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
21: with Standard_Complex_Solutions; use Standard_Complex_Solutions;
22: with Homotopy;
23: with Drivers_for_Poly_Continuation; use Drivers_for_Poly_Continuation;
24: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
25: with Standard_Root_Refiners; use Standard_Root_Refiners;
26: with Brackets,Brackets_io; use Brackets,Brackets_io;
27: with Localization_Posets; use Localization_Posets;
28: with Localization_Posets_io; use Localization_Posets_io;
29: with Curves_into_Grassmannian; use Curves_into_Grassmannian;
30: with Curves_into_Grassmannian_io; use Curves_into_Grassmannian_io;
31: with Deformation_Posets; use Deformation_Posets;
32: with Determinantal_Systems; use Determinantal_Systems;
33: with Plane_Representations; use Plane_Representations;
34: with Drivers_for_Input_Planes; use Drivers_for_Input_planes;
35:
36: procedure Driver_for_Quantum_Pieri
37: ( file : in file_type; n,d,q : in natural ) is
38:
39: function Solution_Plane
40: ( top,bottom : Bracket;
41: locmap : Standard_Natural_Matrices.Matrix;
42: mat : Standard_Complex_Matrices.Matrix )
43: return Solution is
44:
45: -- DESCRIPTION :
46: -- Returns the representation of the solution plane as a vector.
47:
48: solloc : constant Standard_Complex_Vectors.Vector
49: := Column_Vector_Rep(locmap,Localize(locmap,mat));
50: sol : Solution(solloc'length);
51:
52: begin
53: sol.m := 1;
54: sol.t := Create(0.0);
55: sol.res := 0.0;
56: sol.err := 0.0;
57: sol.rco := 0.0;
58: sol.v := solloc;
59: return sol;
60: end Solution_Plane;
61:
62: function Solution_Planes ( top,bottom : Bracket;
63: locmap : Standard_Natural_Matrices.Matrix;
64: vm : VecMat ) return Solution_List is
65:
66: -- DESCRIPTION :
67: -- Returns the representation of the vector of planes as a solution list.
68:
69: res,res_last : Solution_List;
70:
71: begin
72: for i in vm'range loop
73: Append(res,res_last,Solution_Plane(top,bottom,locmap,vm(i).all));
74: end loop;
75: return res;
76: end Solution_Planes;
77:
78: function Create_Polynomial_System
79: ( top,bottom : Bracket;
80: locmap : Standard_Natural_Matrices.Matrix;
81: xpm : Standard_Complex_Poly_Matrices.Matrix;
82: svals : Standard_Complex_Vectors.Vector;
83: planes : VecMat ) return Poly_Sys is
84:
85: -- DESCRIPTION :
86: -- Returns the polynomial system that collects the intersection
87: -- conditions for meeting the given m-planes at the specified s-values.
88: -- The system is localized according to the given localization map.
89:
90: res,wrksys : Poly_Sys(svals'range);
91:
92: begin
93: for i in svals'range loop
94: declare
95: eva : Standard_Complex_Poly_Matrices.Matrix(xpm'range(1),xpm'range(2))
96: := Elim(xpm,svals(i),Create(1.0));
97: wrk : constant Poly_Sys := Polynomial_Equations(planes(i).all,eva);
98: begin
99: wrksys(i) := wrk(wrk'first);
100: Standard_Complex_Poly_Matrices.Clear(eva);
101: end;
102: end loop;
103: res := Column_Localize(top,bottom,locmap,wrksys);
104: Clear(wrksys);
105: return res;
106: end Create_Polynomial_System;
107:
108: procedure Refine_Roots ( file : in file_type;
109: p : in Poly_Sys; sols : in out Solution_List ) is
110:
111: -- DESCRIPTION :
112: -- Calls the root refiner, for square polynomial systems only.
113:
114: epsxa : constant double_float := 10.0**(-8);
115: epsfa : constant double_float := 10.0**(-8);
116: tolsing : constant double_float := 10.0**(-8);
117: max : constant natural := 3;
118: numit : natural := 0;
119:
120: begin
121: Reporting_Root_Refiner(file,p,sols,epsxa,epsfa,tolsing,numit,max,false);
122: end Refine_Roots;
123:
124: procedure Solve_Target_System
125: ( file : in file_type; start,target : in Poly_Sys;
126: sols : in out Solution_List; report : in boolean ) is
127:
128: -- DESCRIPTION :
129: -- Calls the standard continuation routines to solve a specific
130: -- target system, starting at the solutions of the start system.
131:
132: -- REQUIRED : not Is_Null(sols).
133:
134: timer : Timing_Widget;
135: n : constant natural := target'last;
136: a : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
137: b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
138:
139: begin
140: tstart(timer);
141: Homotopy.Create(target,start,1,a,b,true); -- linear cheater
142: Set_Continuation_Parameter(sols,Create(0.0));
143: declare
144: procedure Sil_Cont is
145: new Silent_Continue(Max_Norm,
146: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
147: procedure Rep_Cont is
148: new Reporting_Continue(Max_Norm,
149: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
150: begin
151: if report
152: then Rep_Cont(file,sols,false,Create(1.0));
153: else Sil_Cont(sols,false,Create(1.0));
154: end if;
155: end;
156: tstop(timer);
157: new_line(file);
158: print_times(file,timer,"Cheater's homotopy to target system");
159: Refine_Roots(file,target,sols);
160: end Solve_Target_System;
161:
162: procedure Solve_Hypersurface_Target_System
163: ( file : in file_type; m,p,q : in natural;
164: start_svals,target_svals
165: : in Standard_Complex_Vectors.Vector;
166: start_planes,target_planes : in VecMat;
167: index_poset : in Array_of_Array_of_Nodes;
168: deform_poset : in Array_of_Array_of_VecMats;
169: report : in boolean ) is
170:
171: -- DESCRIPTION :
172: -- This procedure tests the output of the deformation poset,
173: -- creating polynomial representations of the intersection conditions
174: -- and solution lists representing the solution planes.
175:
176: -- ON ENTRY :
177: -- file to write intermediate output on;
178: -- m dimension of the input planes;
179: -- p dimension of the solution planes;
180: -- q degree of the maps;
181: -- start_svals interpolation points at the start;
182: -- target_svals interpolation points at the target;
183: -- start_planes input m-planes in general position;
184: -- target_planes specific input m-planes;
185: -- index_poset indexed localization poset;
186: -- deform_poset poset with the solution p-planes;
187: -- report switch for intermediate output during continuation.
188:
189: dim : constant natural := m*p+q*(m+p);
190: top : constant Bracket := index_poset(dim)(1).top;
191: bot : constant Bracket := index_poset(dim)(1).bottom;
192: xpm : Standard_Complex_Poly_Matrices.Matrix(1..m+p,1..p)
193: := Symbolic_Create(m,p,q,top,bot);
194: solplanes : constant VecMat := deform_poset(dim)(1).all;
195: locmap : Standard_Natural_Matrices.Matrix(1..(m+p)*(q+1),1..p)
196: := Standard_Coordinate_Frame(m,p,q,top,bot,solplanes(1).all);
197: locsys : Poly_Sys(start_planes'range)
198: := Create_Polynomial_System
199: (top,bot,locmap,xpm,start_svals,start_planes);
200: target : Poly_Sys(target_planes'range)
201: := Create_Polynomial_System
202: (top,bot,locmap,xpm,target_svals,target_planes);
203: sols : Solution_List := Solution_Planes(top,bot,locmap,solplanes);
204:
205: begin
206: One_Set_up_Symbol_Table(m,p,q,top,bot);
207: new_line(file);
208: put(file,"The "); put(file,q,1); put(file,"-map of "); put(file,p,1);
209: put_line(file,"-planes representation : ");
210: put(file,xpm);
211: put_line(file,"with as localization map :"); put(file,locmap);
212: new_line(file);
213: Reduce_Symbols(top,bot,locmap);
214: put_line(file,"THE GENERIC SYSTEM : ");
215: put_line(file,locsys);
216: new_line(file);
217: Refine_Roots(file,locsys,sols);
218: new_line(file);
219: put_line(file,"THE TARGET SYSTEM : ");
220: put_line(file,target);
221: new_line(file);
222: Solve_Target_System(file,locsys,target,sols,report);
223: end Solve_Hypersurface_Target_System;
224:
225: procedure Set_Parameters ( file : in file_type; report : out boolean ) is
226:
227: -- DESCRIPTION :
228: -- Interactive determination of the continuation and output parameters.
229:
230: oc : natural;
231:
232: begin
233: new_line;
234: Driver_for_Continuation_Parameters(file);
235: new_line;
236: Driver_for_Process_io(file,oc);
237: report := not (oc = 0);
238: new_line;
239: put_line("No more input expected. See output file for results...");
240: new_line;
241: new_line(file);
242: end Set_Parameters;
243:
244: procedure Write_Poset_Times
245: ( file : in file_type; timer : in Timing_Widget;
246: npaths : in Standard_Natural_Vectors.Vector;
247: timings : in Duration_Array ) is
248:
249: -- DESCRIPTION :
250: -- Writes a overview of #paths and timings spent during deformations
251: -- along the poset structure.
252:
253: begin
254: new_line(file);
255: put_line(file,"--------------------------------------");
256: put_line(file,"| TIMING INFORMATION OVERVIEW |");
257: put_line(file,"--------------------------------------");
258: put_line(file,"| n | #paths | user cpu time |");
259: put_line(file,"--------------------------------------");
260: for i in npaths'range loop
261: if npaths(i) /= 0
262: then put(file,"|"); put(file,i,4); put(file," |");
263: put(file,npaths(i),7); put(file," | ");
264: print_hms(file,timings(i)); put(file," |"); new_line(file);
265: end if;
266: end loop;
267: put_line(file,"--------------------------------------");
268: put(file,"| total |");
269: put(file,Standard_Natural_Vectors.Sum(npaths),7);
270: put(file," | ");
271: print_hms(file,Elapsed_User_Time(timer));
272: put(file," |"); new_line(file);
273: put_line(file,"--------------------------------------");
274: end Write_Poset_Times;
275:
276: procedure Solve_Deformation_Poset
277: ( file : in file_type; m,p,q : in natural;
278: index_poset : in out Array_of_Array_of_Nodes ) is
279:
280: -- DESCRIPTION :
281: -- Writes the symbolic form of the maps.
282:
283: lnd : Link_to_Node;
284: deform_poset : Array_of_Array_of_VecMats(index_poset'range)
285: := Create(index_poset);
286: dim : constant natural := (m*p)+q*(m+p);
287: input : VecMat(1..dim) := Random_Complex_Planes(m,p,q);
288: svals : Standard_Complex_Vectors.Vector(1..dim) := Random_Vector(1,dim);
289: target_planes : VecMat(1..dim);
290: target_svals : Standard_Floating_Vectors.Vector(1..dim);
291: comp_target_svals : Standard_Complex_Vectors.Vector(1..dim);
292: root : Node := index_poset(dim)(1).all;
293: ans : character;
294: report,outlog : boolean;
295: timer : Timing_Widget;
296: npaths : Standard_Natural_Vectors.Vector(1..dim) := (1..dim => 0);
297: timings : Duration_Array(1..dim) := (1..dim => 0.0);
298:
299: begin
300: new_line;
301: put("Do you want to have the homotopies on file ? (y/n) ");
302: Ask_Yes_or_No(ans);
303: outlog := (ans = 'y');
304: Driver_for_Input_Planes(file,m,p,q,target_svals,target_planes);
305: Set_Parameters(file,report);
306: tstart(timer);
307: Solve(file,m+p,q,deform_poset,root,input,svals,report,outlog,
308: npaths,timings);
309: tstop(timer);
310: new_line(file);
311: print_times(file,timer,"Solving along the deformation poset");
312: Write_Poset_Times(file,timer,npaths,timings);
313: for i in target_svals'range loop
314: comp_target_svals(i) := Create(target_svals(i));
315: end loop;
316: Solve_Hypersurface_Target_System
317: (file,m,p,q,svals,comp_target_svals,input,target_planes,
318: index_poset,deform_poset,report);
319: end Solve_Deformation_Poset;
320:
321: procedure Create_Hypersurface_Localization_Poset
322: ( file : in file_type;
323: lnkroot : in Link_to_Node; m,p,q : in natural ) is
324:
325: -- DESCRIPTION :
326: -- Creates the posets and outputs them to the screen and on file.
327: -- Calls the solver afterwards.
328:
329: nq : constant natural := m*p + q*(m+p);
330: level_poset : Array_of_Nodes(0..nq);
331: index_poset : Array_of_Array_of_Nodes(0..nq);
332: nbp : natural;
333:
334: begin
335: level_poset := Create_Leveled_Poset(lnkroot);
336: Count_Roots(level_poset);
337: index_poset := Create_Indexed_Poset(level_poset);
338: put(index_poset);
339: put(file,index_poset);
340: put_line("The size of the poset : "); put_roco(index_poset);
341: put_line(file,"The size of the poset : "); put_roco(file,index_poset);
342: nbp := Root_Count_Sum(level_poset);
343: put("The number of paths : "); put(nbp,1); new_line;
344: put(file,"The number of paths : "); put(file,nbp,1); new_line(file);
345: Solve_Deformation_Poset(file,m,p,q,index_poset);
346: end Create_Hypersurface_Localization_Poset;
347:
348: procedure Create_General_Localization_Poset
349: ( file : in file_type; lnkroot : in Link_to_Node;
350: m,p,q : in natural; codim : in Bracket ) is
351:
352: -- DESCRIPTION :
353: -- Creates the posets and outputs them to the screen and on file.
354:
355: nq : constant natural := m*p + q*(m+p);
356: level_poset : Array_of_Nodes(0..nq);
357: index_poset : Array_of_Array_of_Nodes(0..nq);
358: nbp : natural;
359:
360: begin
361: level_poset := Create_Leveled_Poset(lnkroot);
362: Count_Roots(level_poset);
363: index_poset := Create_Indexed_Poset(level_poset);
364: put(index_poset);
365: put(file,index_poset);
366: put_line("The size of the poset : "); put_roco(index_poset);
367: put_line(file,"The size of the poset : "); put_roco(file,index_poset);
368: nbp := Root_Count_Sum(level_poset);
369: put("The number of paths : "); put(nbp,1); new_line;
370: put(file,"The number of paths : "); put(file,nbp,1); new_line(file);
371: -- Solve_Deformation_Poset(file,m,p,q,index_poset);
372: end Create_General_Localization_Poset;
373:
374: procedure Create_Top_Hypersurface_Poset
375: ( file : in file_type; m,p,q : in natural ) is
376:
377: -- DESCRIPTION :
378: -- Creates the poset by incrementing only top pivots.
379:
380: timer : Timing_Widget;
381: root : Node(p) := Trivial_Root(m,p,q);
382: lnkroot : Link_to_Node := new Node'(root);
383:
384: begin
385: tstart(timer);
386: Q_Top_Create(lnkroot,root.bottom(p),m+p);
387: put_line("The poset created from the top : ");
388: put_line(file,"The poset created from the top : ");
389: Create_Hypersurface_Localization_Poset(file,lnkroot,m,p,q);
390: tstop(timer);
391: new_line(file);
392: print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
393: end Create_Top_Hypersurface_Poset;
394:
395: procedure Create_Bottom_Hypersurface_Poset
396: ( file : in file_type; m,p,q : in natural ) is
397:
398: -- DESCRIPTION :
399: -- Creates the poset by decrementing only bottom pivots.
400:
401: timer : Timing_Widget;
402: root : Node(p) := Trivial_Root(m,p,q);
403: lnkroot : Link_to_Node := new Node'(root);
404:
405: begin
406: tstart(timer);
407: Q_Bottom_Create(lnkroot,m+p);
408: put_line("The poset created from the bottom : ");
409: put_line(file,"The poset created from the bottom : ");
410: Create_Hypersurface_Localization_Poset(file,lnkroot,m,p,q);
411: tstop(timer);
412: new_line(file);
413: print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
414: end Create_Bottom_Hypersurface_Poset;
415:
416: procedure Create_Mixed_Hypersurface_Poset
417: ( file : in file_type; m,p,q : in natural ) is
418:
419: -- DESCRIPTION :
420: -- Creates the poset by incrementing top and decrementing bottom pivots.
421:
422: timer : Timing_Widget;
423: root : Node(p) := Trivial_Root(m,p,q);
424: lnkroot : Link_to_Node := new Node'(root);
425:
426: begin
427: tstart(timer);
428: Q_Top_Bottom_Create(lnkroot,root.bottom(p),m+p);
429: put_line("The poset created in a mixed fashion : ");
430: put_line(file,"The poset created in a mixed fashion :");
431: Create_Hypersurface_Localization_Poset(file,lnkroot,m,p,q);
432: tstop(timer);
433: new_line(file);
434: print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
435: end Create_Mixed_Hypersurface_Poset;
436:
437: procedure Create_Top_General_Poset
438: ( file : in file_type; m,p,q : in natural ) is
439:
440: -- DESCRIPTION :
441: -- Creates the poset by incrementing top pivots.
442:
443: timer : Timing_Widget;
444: root : Node(p) := Trivial_Root(m,p,q);
445: lnkroot : Link_to_Node := new Node'(root);
446: codim : constant Bracket := Read_Codimensions(m,p,q);
447:
448: begin
449: tstart(timer);
450: Q_Top_Create(lnkroot,codim,root.bottom(p),m+p);
451: put_line("The poset created from the top : ");
452: put_line(file,"The poset created from the top :");
453: Create_General_Localization_Poset(file,lnkroot,m,p,q,codim);
454: tstop(timer);
455: new_line(file);
456: print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
457: end Create_Top_General_Poset;
458:
459: procedure Create_Bottom_General_Poset
460: ( file : in file_type; m,p,q : in natural ) is
461:
462: -- DESCRIPTION :
463: -- Creates the poset by decrementing bottom pivots.
464:
465: timer : Timing_Widget;
466: root : Node(p) := Trivial_Root(m,p,q);
467: lnkroot : Link_to_Node := new Node'(root);
468: codim : constant Bracket := Read_Codimensions(m,p,q);
469:
470: begin
471: tstart(timer);
472: Q_Bottom_Create(lnkroot,codim,m+p);
473: put_line("The poset created from the bottom : ");
474: put_line(file,"The poset created from the bottom :");
475: Create_General_Localization_Poset(file,lnkroot,m,p,q,codim);
476: tstop(timer);
477: new_line(file);
478: print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
479: end Create_Bottom_General_Poset;
480:
481: procedure Create_Mixed_General_Poset
482: ( file : in file_type; m,p,q : in natural ) is
483:
484: -- DESCRIPTION :
485: -- Creates the poset by incrementing top and decrementing bottom pivots.
486:
487: timer : Timing_Widget;
488: root : Node(p) := Trivial_Root(m,p,q);
489: lnkroot : Link_to_Node := new Node'(root);
490: codim : constant Bracket := Read_Codimensions(m,p,q);
491:
492: begin
493: tstart(timer);
494: Q_Top_Bottom_Create(lnkroot,codim,root.bottom(p),m+p);
495: put_line("The poset created from the bottom : ");
496: put_line(file,"The poset created from the bottom :");
497: Create_General_Localization_Poset(file,lnkroot,m,p,q,codim);
498: tstop(timer);
499: new_line(file);
500: print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
501: end Create_Mixed_General_Poset;
502:
503: procedure Main is
504:
505: p : constant natural := d;
506: m : constant natural := n-d;
507: ans : character;
508:
509: begin
510: new_line;
511: put_line("MENU for interpolating maps of fixed degree in Grassmannian.");
512: put_line(" 1. k_i = 1 consistently incrementing top pivots.");
513: put_line(" 2. consistently decrementing bottom pivots.");
514: put_line(" 3. mixed top-bottom sequence for poset creation.");
515: put_line(" 4. k_i >= 1 consistently incrementing top pivots.");
516: put_line(" 5. consistently incrementing bottom pivots.");
517: put_line(" 6. mixed top-bottom sequence for poset creation.");
518: put("Type 1, 2, 3, 4, 5, or 6 to choose : ");
519: Ask_Alternative(ans,"123456");
520: new_line;
521: case ans is
522: when '1' => Create_Top_Hypersurface_Poset(file,m,p,q);
523: when '2' => Create_Bottom_Hypersurface_Poset(file,m,p,q);
524: when '3' => Create_Mixed_Hypersurface_Poset(file,m,p,q);
525: when '4' => Create_Top_General_Poset(file,m,p,q);
526: when '5' => Create_Bottom_General_Poset(file,m,p,q);
527: when '6' => Create_Mixed_General_Poset(file,m,p,q);
528: when others => put_line("Option not recognized. Please try again.");
529: end case;
530: end Main;
531:
532: begin
533: Main;
534: end Driver_for_Quantum_Pieri;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to driver_for_quantum_pieri.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