Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_org_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_Natural_Vectors; use Standard_Natural_Vectors;
6: with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
7: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
8: with Standard_Complex_Vectors;
9: with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
10: with Standard_Complex_Matrices;
11: with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
12: with Standard_Complex_VecMats; use Standard_Complex_VecMats;
13: with Standard_Complex_VecMats_io; use Standard_Complex_VecMats_io;
14: with Standard_Random_Matrices; use Standard_Random_Matrices;
15: with Standard_Complex_Poly_Matrices;
16: with Standard_Complex_Poly_Matrices_io; use Standard_Complex_Poly_Matrices_io;
17: with Symbol_Table; use Symbol_Table;
18: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
19: with Standard_Complex_Polynomials_io; use Standard_Complex_Polynomials_io;
20: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
21: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
22: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
23: with Matrix_Indeterminates; use Matrix_Indeterminates;
24: with Drivers_for_Poly_Continuation; use Drivers_for_Poly_Continuation;
25: with Brackets,Brackets_io; use Brackets,Brackets_io;
26: with Bracket_Monomials; use Bracket_Monomials;
27: with Bracket_Monomials_io; use Bracket_Monomials_io;
28: with Bracket_Expansions; use Bracket_Expansions;
29: with Bracket_Polynomials; use Bracket_Polynomials;
30: with Bracket_Polynomials_io; use Bracket_Polynomials_io;
31: with Bracket_Systems; use Bracket_Systems;
32: with Pieri_Trees,Pieri_Trees_io; use Pieri_Trees,Pieri_Trees_io;
33: with Pieri_Root_Counts; use Pieri_Root_Counts;
34: with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
35: with Numeric_Minor_Equations; use Numeric_Minor_Equations;
36: with Solve_Pieri_Leaves;
37: with Specialization_of_Planes; use Specialization_of_Planes;
38: with Pieri_Deformations; use Pieri_Deformations;
39:
40: procedure ts_org_pieri is
41:
42: -- DESCRIPTION :
43: -- This program executes the original implementation of the Pieri homotopy
44: -- algorithm, implemented strictly along the lines of the description in
45: -- the paper "Numerical Schubert Calculus" of Birkett Huber, Frank Sottile
46: -- and Bernd Sturmfels.
47:
48: function Maximum ( n1,n2 : natural ) return natural is
49: begin
50: if n1 >= n2
51: then return n1;
52: else return n2;
53: end if;
54: end Maximum;
55:
56: procedure Add_t_Symbol is
57:
58: -- DESCRIPTION :
59: -- Adds the symbol for the continuation parameter t to the symbol table.
60:
61: tsb : Symbol;
62:
63: begin
64: Symbol_Table.Enlarge(1);
65: tsb(1) := 't';
66: for i in 2..tsb'last loop
67: tsb(i) := ' ';
68: end loop;
69: Symbol_Table.Add(tsb);
70: end Add_t_Symbol;
71:
72: -- DISPLAYING THE REPRESENTATIONS OF THE PLANES :
73:
74: procedure Display_Polynomial_Pattern
75: ( file : in file_type; n : in natural; b1,b2 : in Bracket ) is
76:
77: -- DESCRIPTION :
78: -- Displays the pattern by a polynomial matrix.
79:
80: pm : Standard_Complex_Poly_Matrices.Matrix(1..n,b1'range)
81: := Schubert_Pattern(n,b1,b2);
82:
83: begin
84: put(file,pm);
85: Standard_Complex_Poly_Matrices.Clear(pm);
86: end Display_Polynomial_Pattern;
87:
88: -- AUXILIARIES TO SET UP EQUATIONS :
89:
90: procedure Expand_Minors ( file : in file_type;
91: mat : in Standard_Complex_Poly_Matrices.Matrix;
92: bm : in Bracket_Monomial ) is
93:
94: -- DESCRIPTION :
95: -- Expands the quadratic bracket monomial.
96:
97: first : boolean := true;
98: lb : Link_to_Bracket;
99:
100: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
101:
102: p : Poly := Null_Poly;
103:
104: begin
105: if first
106: then lb := new Bracket'(b);
107: first := false;
108: else p := Expanded_Minor(mat,b);
109: if p /= Null_Poly
110: then put(file,lb.all);
111: put(file,"*("); put(file,p); put(file,")");
112: end if;
113: Clear(lb); Clear(p);
114: end if;
115: continue := true;
116: end Visit_Bracket;
117: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
118:
119: begin
120: Visit_Brackets(bm);
121: end Expand_Minors;
122:
123: procedure Expand_Minors ( file : in file_type;
124: mat : in Standard_Complex_Poly_Matrices.Matrix;
125: bp : Bracket_Polynomial ) is
126:
127: -- DESCRIPTION :
128: -- Writes the expansion of the matrix, using the bracket polynomial which
129: -- is a list of quadratic monomials that represent the Laplace expansion.
130:
131: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
132: begin
133: if REAL_PART(t.coeff) > 0.0
134: then put("+");
135: else put("-");
136: end if;
137: Expand_Minors(file,mat,t.monom);
138: continue := true;
139: end Visit_Term;
140: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
141:
142: begin
143: Visit_Terms(bp);
144: end Expand_Minors;
145:
146: procedure Pieri_Equations
147: ( file : in file_type; n,d : in natural; bs : in Bracket_System;
148: b1,b2 : in Bracket ) is
149:
150: -- DESCRIPTION :
151: -- Writes the Pieri equations corresponding to the pair of brackets.
152:
153: cffmat : Standard_Complex_Matrices.Matrix(1..n,1..(n-d))
154: := Random_Matrix(n,n-d);
155: polmat : Standard_Complex_Poly_Matrices.Matrix(1..n,1..d)
156: := Schubert_Pattern(n,b1,b2);
157: sys : Poly_Sys(bs'first+1..bs'last);
158: sol : Standard_Complex_Matrices.Matrix(1..n,1..d);
159:
160: begin
161: put(file,"Plane X for (");
162: put(file,b1); put(file,","); put(file,b2); put_line(file,") :");
163: put(file,polmat);
164: put_line(file,"The system with expanded minors of X : ");
165: for i in 1..bs'last loop
166: Expand_Minors(file,polmat,bs(i)); new_line(file);
167: end loop;
168: -- put("Give a "); put(n,1); put("x"); put(n-d,1);
169: -- put_line("-matrix of complex numbers : ");
170: -- get(cffmat);
171: -- put("Your "); put(n-d,1); put_line("-plane : "); put(cffmat);
172: -- put("A random "); put(n-d,1); put_line("-plane : "); put(cffmat);
173: put("The minors evaluated at a random ");
174: put(n-d,1); put_line("-plane :");
175: sys := Expanded_Minors(cffmat,polmat,bs);
176: put_line(sys);
177: Standard_Complex_Poly_Matrices.Clear(polmat);
178: if Pieri_Condition(n,b1,b2)
179: then put("The "); put(d,1); put_line("-plane at the leaves :");
180: sol := Solve_Pieri_Leaves(Standard_Output,b1,b2,cffmat); put(sol);
181: put_line("The solution evaluated at the system : ");
182: put(Evaluate(sys,sol)); new_line;
183: else put_line("Pair of leaves does not satisfy Pieri's condition.");
184: end if;
185: end Pieri_Equations;
186:
187: procedure Expand_Minors ( n,d : in natural; bs : in Bracket_System ) is
188:
189: -- DESCRIPTION :
190: -- Expands the minors to obtain a symbolic formulation of the equations.
191:
192: b1,b2 : Bracket(1..d);
193:
194: begin
195: put("Give 1st bracket : "); get(b1);
196: put("Give 2nd bracket : "); get(b2);
197: Pieri_Equations(Standard_Output,n,d,bs,b1,b2);
198: end Expand_Minors;
199:
200: procedure Pieri_Equations_for_Paired_Chains
201: ( file : in file_type; n,d : in natural; bs : in Bracket_System;
202: b1,b2 : in bracket_Array ) is
203:
204: -- DESCRIPTION :
205: -- Writes the equations for each node along a pair of chains.
206:
207: maxlen : constant natural := Maximum(b1'last,b2'last);
208: lb1,lb2 : Link_to_Bracket;
209:
210: begin
211: for i in b1'first..maxlen loop
212: if i <= b1'last
213: then lb1 := b1(i);
214: else lb1 := b1(b1'last);
215: end if;
216: if i <= b2'last
217: then lb2 := b2(i);
218: else lb2 := b2(b2'last);
219: end if;
220: Pieri_Equations(file,n,d,bs,lb1.all,lb2.all);
221: end loop;
222: end Pieri_Equations_for_Paired_Chains;
223:
224: procedure Write_Pieri_Equations
225: ( file : in file_type; n,d : in natural; t1,t2 : in Pieri_Tree;
226: bs : in Bracket_System ) is
227:
228: -- DESCRIPTION :
229: -- Writes the Pieri Equations for all pairs of chains.
230:
231: procedure Visit_Chain ( b1,b2 : in Bracket_Array; cont : out boolean ) is
232: begin
233: Pieri_Equations_for_Paired_Chains(file,n,d,bs,b1,b2);
234: cont := true;
235: end Visit_Chain;
236: procedure Visit_Chains is new Enumerate_Paired_Chains(Visit_Chain);
237:
238: begin
239: Visit_Chains(t1,t2);
240: end Write_Pieri_Equations;
241:
242: procedure Write_Pieri_Equations
243: ( n,d : in natural; t1,t2 : in Pieri_Tree ) is
244:
245: k,kd : natural;
246: bm : Bracket_Monomial;
247: file : file_type;
248:
249: begin
250: new_line; skip_line;
251: put_line("Reading the name of the output file.");
252: Read_Name_and_Create_File(file);
253: put("Give k to determine (m-k+1)-plane : "); get(k);
254: kd := n-k+1;
255: bm := Maximal_Minors(n,kd); -- because n = m+p
256: put(file,"All maximal minors : "); put(file,bm); new_line(file);
257: declare
258: bs : constant Bracket_System := Minor_Equations(kd,kd-d,bm);
259: begin
260: put_line(file,"The generic equation in the Laplace expansion : ");
261: put(file,bs(0));
262: put_line(file,"The specific equations in the system : ");
263: for i in 1..bs'last loop
264: put(file,bs(i));
265: end loop;
266: Write_Pieri_Equations(file,n,d,t1,t2,bs);
267: end;
268: Close(file);
269: Clear(bm);
270: end Write_Pieri_Equations;
271:
272: -- AUXILIARIES TO REPRESENT PIERI TREES :
273:
274: procedure Write_a_Chain ( file : in file_type; b : in Bracket_Array ) is
275: begin
276: put(b(b'first).all);
277: for i in b'first+1..b'last loop
278: put(" < "); put(b(i).all);
279: end loop;
280: new_line;
281: end Write_a_Chain;
282:
283: procedure Write_Chains ( t : in Pieri_Tree ) is
284:
285: procedure Write_Chain ( b : in Bracket_Array; cont : out boolean ) is
286: begin
287: Write_a_Chain(Standard_Output,b);
288: cont := true;
289: end Write_Chain;
290: procedure Write is new Enumerate_Chains(Write_Chain);
291:
292: begin
293: Write(t);
294: end Write_Chains;
295:
296: -- AUXILIARIES TO COUNT THE ROOTS :
297:
298: procedure Write_Polynomial_Patterns
299: ( file : in file_type;
300: n : in natural; b1,b2 : in Bracket_Array ) is
301:
302: -- DESCRIPTION :
303: -- Writes the two chains in a paired fashion. If they have unequal
304: -- length, then the last element of the shortest chain appears repeated.
305:
306: maxlen : constant natural := Maximum(b1'last,b2'last);
307: lb1,lb2 : Link_to_Bracket;
308:
309: begin
310: for i in b1'first..maxlen loop
311: put(file,"(");
312: if i <= b1'last
313: then lb1 := b1(i);
314: else lb1 := b1(b1'last);
315: end if;
316: put(file,lb1.all); put(file,",");
317: if i <= b2'last
318: then lb2 := b2(i);
319: else lb2 := b2(b2'last);
320: end if;
321: put(file,lb2.all); put_line(file,") has pattern : ");
322: Display_Polynomial_Pattern(file,n,lb1.all,lb2.all);
323: end loop;
324: end Write_Polynomial_Patterns;
325:
326: procedure Write_Paired_Chain
327: ( file : in file_type;
328: n : in natural; b1,b2 : in Bracket_Array ) is
329:
330: -- DESCRIPTION :
331: -- Writes the two chains in a paired fashion. If they have unequal
332: -- length, then the last element of the shortest chain appears repeated.
333:
334: maxlen : constant natural := Maximum(b1'last,b2'last);
335:
336: begin
337: put(file,"("); put(file,b1(b1'first).all); put(file,",");
338: put(file,b2(b2'first).all); put(file,")");
339: for i in b1'first+1..maxlen loop
340: put(file," < "); put(file,"(");
341: if i <= b1'last
342: then put(file,b1(i).all);
343: else put(file,b1(b1'last).all);
344: end if;
345: put(file,",");
346: if i <= b2'last
347: then put(file,b2(i).all);
348: else put(file,b2(b2'last).all);
349: end if;
350: put(file,")");
351: end loop;
352: new_line(file);
353: Write_Polynomial_Patterns(Standard_Output,n,b1,b2);
354: if Pieri_Condition(n,b1(b1'last).all,b2(b2'last).all)
355: then put_line("Leaves satisfy Pieri's condition.");
356: else put_line("Leaves do not satisfy Pieri's condition.");
357: end if;
358: end Write_Paired_Chain;
359:
360: procedure Write_Pieri_Chains ( n : in natural; t1,t2 : in Pieri_Tree ) is
361:
362: -- DESCRIPTION :
363: -- Enumerates all pairs of chains and checks Pieri's condition
364: -- at the leaves.
365:
366: procedure Visit_Chain ( b1,b2 : in Bracket_Array; cont : out boolean ) is
367: begin
368: Write_Paired_Chain(Standard_Output,n,b1,b2);
369: cont := true;
370: end Visit_Chain;
371: procedure Visit_Chains is new Enumerate_Paired_Chains(Visit_Chain);
372:
373: begin
374: Visit_Chains(t1,t2);
375: end Write_Pieri_Chains;
376:
377: function First_Standard_Plane
378: ( n,m,r : natural ) return Standard_Complex_Matrices.Matrix is
379:
380: -- DESCRIPTION :
381: -- Returns the plane spanned by the first m+1-r standard basis vectors.
382:
383: res : Standard_Complex_Matrices.Matrix := Random_Matrix(n,m+1-r);
384:
385: begin
386: for i in res'range(1) loop
387: for j in res'range(2) loop
388: if i = j
389: then res(i,j) := Create(1.0);
390: else res(i,j) := Create(0.0);
391: end if;
392: end loop;
393: end loop;
394: return res;
395: end First_Standard_Plane;
396:
397: function Last_Standard_Plane
398: ( n,m,r : natural ) return Standard_Complex_Matrices.Matrix is
399:
400: -- DESCRIPTION :
401: -- Returns the plane spanned by the first m+1-r standard basis vectors.
402:
403: res : Standard_Complex_Matrices.Matrix := Random_Matrix(n,m+1-r);
404:
405: begin
406: for i in res'range(1) loop
407: for j in res'range(2) loop
408: if i = res'last(2) + 1 - j
409: then res(i,j) := Create(1.0);
410: else res(i,j) := Create(0.0);
411: end if;
412: end loop;
413: end loop;
414: return res;
415: end Last_Standard_Plane;
416:
417: function First_Random_Input_Sequence
418: ( n,m,a : natural; kp : Vector ) return VecMat is
419:
420: -- DESCRIPTION :
421: -- Returns the first sequence of random input planes. The first plane
422: -- in the sequence is spanned by the first standard basis vectors.
423: -- The dimensions of the planes are m+1-kp(i), for the appropriate i.
424:
425: res : VecMat(0..a-1);
426:
427: begin
428: res(0) := new Standard_Complex_Matrices.Matrix'
429: (First_Standard_Plane(n,m,kp(1)));
430: for i in 1..res'last loop
431: res(i) := new Standard_Complex_Matrices.Matrix'
432: (Random_Matrix(n,m+1-kp(i+1)));
433: end loop;
434: return res;
435: end First_Random_Input_Sequence;
436:
437: function Second_Random_Input_Sequence
438: ( n,m,a : natural; kp : Vector ) return VecMat is
439:
440: -- DESCRIPTION :
441: -- Returns the second sequence of random input planes. The first plane
442: -- in the sequence is spanned by the last standard basis vectors.
443:
444: res : VecMat(0..kp'last-a-2);
445:
446: begin
447: res(0) := new Standard_Complex_Matrices.Matrix'
448: (Last_Standard_Plane(n,m,kp(1)));
449: for i in 1..res'last loop
450: res(i) := new Standard_Complex_Matrices.Matrix'
451: (Random_Matrix(n,m+1-kp(a+1+i)));
452: end loop;
453: return res;
454: end Second_Random_Input_Sequence;
455:
456: procedure Reallify ( c : in out Complex_Number ) is
457:
458: -- DESCRIPTION :
459: -- Sets the imaginary part of c to zero.
460:
461: begin
462: c := Create(REAL_PART(c),0.0);
463: end Reallify;
464:
465: procedure Reallify ( m : in out Standard_Complex_Matrices.Matrix ) is
466:
467: -- DESCRIPTION :
468: -- Sets the imaginary part of every entry in the matrix to zero.
469:
470: begin
471: for i in m'range(1) loop
472: for j in m'range(2) loop
473: Reallify(m(i,j));
474: end loop;
475: end loop;
476: end Reallify;
477:
478: procedure Reallify ( v : in out VecMat ) is
479:
480: -- DESCRIPTION :
481: -- Sets the imaginary part of every entry of every matrix in v to zero
482: -- and makes the matrices orthonormal.
483:
484: begin
485: for i in v'range loop
486: Reallify(v(i).all);
487: v(i).all := Orthogonalize(v(i).all);
488: end loop;
489: end Reallify;
490:
491: procedure Orthogonalize ( v : in out VecMat ) is
492:
493: -- DESCRIPTION :
494: -- Orthonormalizes every matrix in the array.
495:
496: begin
497: for i in v'range loop
498: v(i).all := Orthogonalize(v(i).all);
499: end loop;
500: end Orthogonalize;
501:
502: procedure Set_Parameters ( file : in file_type; report : out boolean ) is
503:
504: -- DESCRIPTION :
505: -- Interactive determination of the continuation and output parameters.
506:
507: oc : natural;
508:
509: begin
510: new_line;
511: Driver_for_Continuation_Parameters(file);
512: new_line;
513: Driver_for_Process_io(file,oc);
514: report := not (oc = 0);
515: end Set_Parameters;
516:
517: function Select_Pairs ( lps : List_of_Paired_Nodes )
518: return List_of_Paired_Nodes is
519:
520: -- DESCRIPTION :
521: -- Returns a selection of the list of pairs.
522:
523: res,res_last : List_of_Paired_Nodes;
524: k : natural;
525: tmp : List_of_Paired_Nodes := lps;
526:
527: begin
528: put("Give the number of pairs : "); get(k);
529: declare
530: sel : Standard_Natural_Vectors.Vector(1..k);
531: ind : natural := 1;
532: begin
533: put("Give an increasing sequence of "); put(k,1); put(" numbers : ");
534: get(sel);
535: for i in 1..Length_Of(lps) loop
536: if i = sel(ind)
537: then ind := ind+1;
538: Append(res,res_last,Head_Of(tmp));
539: end if;
540: tmp := Tail_Of(tmp);
541: end loop;
542: end;
543: skip_line;
544: return res;
545: end Select_Pairs;
546:
547: procedure put ( file : in file_type; lp : in List_of_Paired_Nodes ) is
548:
549: tmp : List_of_Paired_Nodes := lp;
550: pnd : Paired_Nodes;
551:
552: begin
553: put_line(file,"The pairs of leaves that satisfy Pieri's condition :");
554: for i in 1..Length_Of(lp) loop
555: pnd := Head_Of(tmp);
556: put(file,"Leaf "); put(file,i,3); put(file," : ");
557: put(file,"("); put(file,pnd.left.node); put(file,",");
558: put(file,pnd.right.node); put_line(file,")");
559: tmp := Tail_Of(tmp);
560: end loop;
561: end put;
562:
563: procedure Root_Count ( n,d,a : in natural; kp : in Vector ) is
564:
565: -- DESCRIPTION :
566: -- Set up of Pieri trees from a partition of the planes.
567:
568: -- ON ENTRY :
569: -- n the dimension of the space we are working in;
570: -- d dimension of the brackets, of the output planes;
571: -- a number of planes in the first set of the partition;
572: -- kp the dimensions of the input planes are m+1-kp(i),
573: -- where m = n-d.
574:
575: file : file_type;
576: timer : Timing_Widget;
577: ans : character;
578: m : constant natural := n-d;
579: v1 : Vector(0..a-1) := kp(1..a);
580: t1 : Pieri_Tree(d,a-1) := Create(n,d,v1);
581: a2 : constant natural := kp'last-a-1;
582: v2 : Vector(0..a2-1) := kp(a+1..kp'last-1);
583: t2 : Pieri_Tree(d,a2) := Create(n,d,v2);
584: lp : List_of_Paired_Nodes := Create(n,d,t1,t2);
585: np : Nodal_Pair(d) := Create(d,lp);
586: sel_lp : List_of_Paired_Nodes;
587: rc : constant natural := Length_Of(lp);
588: nb : constant natural := Number_of_Paths(np);
589: l1 : VecMat(0..a-1) := First_Random_Input_Sequence(n,m,a,kp);
590: l2 : VecMat(0..kp'last-a-2) := Second_Random_Input_Sequence(n,m,a,kp);
591: ln : Standard_Complex_Matrices.Matrix(1..n,1..m+1-kp(kp'last))
592: := Random_Matrix(n,m+1-kp(kp'last));
593: report,outlog : boolean;
594: sols : VecMat(1..rc);
595:
596: begin
597: skip_line;
598: new_line;
599: put_line("Reading the name of the output file.");
600: Read_Name_and_Create_File(file);
601: new_line;
602: put_line("See the output file for results.");
603: new_line;
604: put(file,"(m = "); put(file,m,1); put(file,",p = "); put(file,d,1);
605: put(file,")");
606: put(file," k = {");
607: for i in 1..a-1 loop
608: put(file,kp(i),1); put(file,",");
609: end loop;
610: put(file,kp(a),1);
611: put(file,"}");
612: put(file,"{");
613: for i in a+1..kp'last-2 loop
614: put(file,kp(i),1); put(file,",");
615: end loop;
616: put(file,kp(kp'last-1),1);
617: put(file,"} ");
618: put(file,kp(kp'last),1);
619: new_line(file);
620: put(file,"The first Pieri tree has height "); put(file,Height(t1),1);
621: put_line(file," :"); Write_Tree(file,t1); -- put(file,t1);
622: put(file,"The second Pieri tree has height "); put(file,Height(t2),1);
623: put_line(file," :"); Write_Tree(file,t2); -- put(file,t2);
624: put(file,"The root count equals : "); put(file,rc,1); new_line(file);
625: put("The root count equals : "); put(rc,1); new_line;
626: put(file,lp);
627: put_line(file,"The tree of paired nodes :");
628: Write(file,np);
629: put("The number of paths : "); put(nb,1); new_line;
630: put(file,"The number of paths : "); put(file,nb,1); new_line(file);
631: new_line;
632: put("Do you want real random input planes ? (y/n) "); get(ans);
633: if ans = 'y'
634: then Reallify(l1); Reallify(l2); Reallify(ln);
635: else Orthogonalize(l1); Orthogonalize(l2);
636: ln := Orthogonalize(ln);
637: end if;
638: put("Do you want moving cycles/polynomial systems on file ? (y/n) ");
639: Ask_Yes_or_No(ans);
640: outlog := (ans = 'y');
641: put("Do you want to select only certain pairs ? (y/n) ? ");
642: Ask_Yes_or_No(ans);
643: if ans = 'y'
644: then sel_lp := Select_Pairs(lp);
645: else sel_lp := lp;
646: end if;
647: Set_Parameters(file,report);
648: new_line(file);
649: put_line(file,"The first sequence of random input planes :");
650: put(file,l1,2);
651: put_line(file,"The second sequence of random input planes :");
652: put(file,l2,2);
653: put_line(file,"The last random input plane :"); put(file,ln,2);
654: put_line(file,"Starting the deformations at the paired leaves.");
655: Matrix_Indeterminates.Initialize_Symbols(n,d);
656: Add_t_Symbol;
657: tstart(timer);
658: Deform_Pairs(file,n,d,sel_lp,l1,l2,ln,report,outlog,sols);
659: tstop(timer);
660: new_line(file);
661: print_times(file,timer,"Pieri Deformations");
662: Matrix_Indeterminates.Clear_Symbols;
663: Clear(t1); Clear(t2); Clear(lp);
664: Clear(l1); Clear(l2); Clear(sols);
665: Close(file);
666: end Root_Count;
667:
668: -- FIVE MAJOR TEST PROGRAMS :
669:
670: procedure Test_Pieri_Condition ( n,d : in natural ) is
671:
672: -- DESCRIPTION :
673: -- Reads two brackets and tests Pieri's condition.
674:
675: b1,b2 : Bracket(1..d);
676: ans : character;
677: cnt : natural := 0;
678:
679: begin
680: new_line;
681: put_line("Interactive test on Pieri's condition for pair of brackets.");
682: Matrix_Indeterminates.Initialize_Symbols(n,d);
683: loop
684: new_line;
685: put("Give first bracket : "); get(b1);
686: put("Give second bracket : "); get(b2);
687: put("("); put(b1); put(","); put(b2); put(")");
688: if Pieri_Condition(n,b1,b2)
689: then put_line(" satisfies Pieri's condition.");
690: else put_line(" does not satisfy Pieri's condition.");
691: end if;
692: put("Do you want more tests ? (y/n) "); get(ans);
693: exit when ans /= 'y';
694: end loop;
695: Matrix_Indeterminates.Clear_Symbols;
696: end Test_Pieri_Condition;
697:
698: procedure Test_Pieri_Tree ( n,d : in natural ) is
699:
700: -- DESCRIPTION : Constructs T(r_0,r_1,..,r_a).
701:
702: a : natural;
703: ans : character;
704:
705: begin
706: new_line;
707: put_line("Interactive test on the construction of Pieri trees");
708: loop
709: new_line;
710: put("Give number of jumped-branching levels : "); get(a);
711: declare
712: v : Vector(0..a);
713: t : Pieri_Tree(d,a);
714: begin
715: put("Give "); put(a+1,1); put(" natural numbers : "); get(v);
716: t := Create(n,d,v);
717: put_line("The Pieri tree : "); Write_Tree(t); --put(t);
718: put_line("The chains in the Pieri tree :"); Write_Chains(t);
719: Clear(t);
720: end;
721: put("Do you want more tests ? (y/n) "); get(ans);
722: exit when ans /= 'y';
723: end loop;
724: end Test_Pieri_Tree;
725:
726: procedure Test_Root_Count ( n,d : in natural ) is
727:
728: -- DESCRIPTION :
729: -- Reads the input planes and sets up a pair of trees to perform
730: -- the combinatorial root count.
731:
732: m : constant natural := n-d;
733: np,sum,a : natural;
734: ans : character;
735:
736: begin
737: new_line;
738: put_line("Counting roots combinatorially by Pieri trees");
739: loop
740: new_line;
741: put("Give the number of planes to intersect : "); get(np);
742: declare
743: kp : Vector(1..np);
744: begin
745: put("Give "); put(np,1); put(" co-dimensions of the planes : ");
746: get(kp);
747: put("m = "); put(m,1); put(" p = "); put(d,1);
748: put(" Sum : "); put(kp(kp'first),1);
749: sum := kp(kp'first);
750: for i in kp'first+1..kp'last loop
751: put(" + "); put(kp(i),1);
752: sum := sum + kp(i);
753: end loop;
754: put(" = "); put(sum,1);
755: if sum = m*d
756: then put_line(" = m*p");
757: loop
758: put("Give number of elements in first set : "); get(a);
759: Root_Count(n,d,a,kp);
760: new_line;
761: put("Do you want to test other partitions ? (y/n) "); get(ans);
762: exit when ans /= 'y';
763: end loop;
764: else put_line(" /= m*p");
765: end if;
766: end;
767: put("Do you want more tests ? (y/n) "); get(ans);
768: exit when ans /= 'y';
769: end loop;
770: end Test_Root_Count;
771:
772: procedure Test_Pieri_Equations ( n,d : in natural ) is
773:
774: -- DESCRIPTION :
775: -- Set up of the expansions of the maximal minors.
776:
777: k,kd,nb : natural;
778: bm : Bracket_Monomial;
779: ans : character;
780: b1,b2 : Bracket(1..d);
781:
782: begin
783: new_line;
784: put_line("Set up equations for certain Schubert conditions.");
785: Matrix_Indeterminates.Initialize_Symbols(n,d);
786: loop
787: new_line;
788: put("Give k to determine (m-k+1)-plane : "); get(k);
789: kd := n-k+1;
790: bm := Maximal_Minors(n,kd); -- because n = m+p
791: put("All maximal minors : "); put(bm); new_line;
792: declare
793: bs : constant Bracket_System := Minor_Equations(kd,kd-d,bm);
794: begin
795: put_line("The generic equation in the Laplace expansion : ");
796: put(bs(0));
797: put_line("The specific equations in the system : ");
798: for i in 1..bs'last loop
799: put(bs(i));
800: end loop;
801: Expand_Minors(n,d,bs);
802: end;
803: put("Do you want more tests ? (y/n) "); get(ans);
804: Clear(bm);
805: exit when (ans /= 'y');
806: end loop;
807: Matrix_Indeterminates.Clear_Symbols;
808: new_line;
809: put("Do you want to test memory consumption ? (y/n) "); get(ans);
810: if ans = 'y'
811: then put("Give number of times : "); get(nb);
812: k := 1;
813: for i in 1..d loop
814: b1(i) := i;
815: b2(i) := i;
816: end loop;
817: for i in 1..nb loop
818: kd := n-k+1;
819: bm := Maximal_Minors(n,kd);
820: declare
821: nva : constant natural := n*d+1;
822: bs : Bracket_System(0..Number_of_Brackets(bm))
823: := Minor_Equations(kd,kd-d,bm);
824: cffmat : Standard_Complex_Matrices.Matrix(1..n,1..(n-d))
825: := Random_Matrix(n,n-d);
826: stamat : Standard_Complex_Matrices.Matrix(1..n,1..(n-d))
827: := Random_Matrix(n,n-d);
828: movmat : Standard_Complex_Poly_Matrices.Matrix
829: (cffmat'range(1),cffmat'range(2))
830: := Moving_U_Matrix(nva,stamat,cffmat);
831: polmat : Standard_Complex_Poly_Matrices.Matrix(1..n,1..d)
832: := Schubert_Pattern(n,b1,b2);
833: sys : Poly_Sys(1..bs'last)
834: := Expanded_Minors(cffmat,polmat,bs);
835: movcyc : Poly_Sys(1..bs'last)
836: := Lifted_Expanded_Minors(movmat,polmat,bs);
837: begin
838: Clear(bm); Clear(bs);
839: Standard_Complex_Poly_Matrices.Clear(movmat);
840: Standard_Complex_Poly_Matrices.Clear(polmat);
841: Clear(sys); Clear(movcyc);
842: end;
843: end loop;
844: end if;
845: end Test_Pieri_Equations;
846:
847: procedure Test_Pieri_Homotopies ( n,d : in natural ) is
848:
849: -- DESCRIPTION :
850: -- Set up of the parametric families in the Pieri homotopy algorithm.
851:
852: b1,b2 : Bracket(1..d);
853: m : constant natural := n-d;
854: nva : constant natural := n*d+1;
855: k,kd : natural;
856: bm : Bracket_Monomial;
857: F1 : constant Standard_Complex_Matrices.Matrix
858: := Random_Upper_Triangular(n);
859: F2 : constant Standard_Complex_Matrices.Matrix
860: := Random_Lower_Triangular(n);
861:
862: begin
863: new_line;
864: put_line("Set up the moving cycles in the Pieri homotopy algorithm.");
865: new_line;
866: Matrix_Indeterminates.Initialize_Symbols(n,d);
867: Add_t_Symbol;
868: put_line("The general upper triangular matrix F : "); put(F1,2);
869: put_line("The general lower triangular matrix F' : "); put(F2,2);
870: put("Give first "); put(d,1); put("-bracket : "); get(b1);
871: put("Give second "); put(d,1); put("-bracket : "); get(b2);
872: put("The pair of brackets : ");
873: put("("); put(b1); put(","); put(b2); put_line(")");
874: put("Give k to determine (m-k+1)-plane : "); get(k);
875: kd := n-k+1;
876: bm := Maximal_Minors(n,kd);
877: declare
878: L : constant Standard_Complex_Matrices.Matrix := Random_Matrix(n,m+1-k);
879: X : Standard_Complex_Poly_Matrices.Matrix(1..n,b1'range);
880: U1 : constant Standard_Complex_Matrices.Matrix := U_Matrix(F1,b1);
881: movU1 : Standard_Complex_Poly_Matrices.Matrix(L'range(1),L'range(2));
882: U2 : constant Standard_Complex_Matrices.Matrix := U_Matrix(F2,b2);
883: movU2 : Standard_Complex_Poly_Matrices.Matrix(L'range(1),L'range(2));
884: bs : constant Bracket_System := Minor_Equations(kd,kd-d,bm);
885: movcyc1,movcyc2 : Poly_Sys(1..bs'last);
886: begin
887: put("The U-matrix for F and "); put(b1); put_line(" :"); put(U1,2);
888: put("The U-matrix for F' and "); put(b2); put_line(" :"); put(U2,2);
889: put("A random "); put(m+1-k,1); put_line("-plane :"); put(L,2);
890: movU1 := Moving_U_Matrix(nva,U1,L);
891: put_line("The first moving U-matrix :"); put(movU1);
892: movU2 := Moving_U_Matrix(nva,U2,L);
893: put_line("The second moving U-matrix :"); put(movU2);
894: X := Schubert_Pattern(n,b1,b2);
895: put("The unknown "); put(d,1); put_line("-plane X :"); put(X);
896: movcyc1 := Lifted_Expanded_Minors(movU1,X,bs);
897: put_line("The polynomial system for the first moving U-matrix :");
898: put_line(movcyc1);
899: movcyc2 := Lifted_Expanded_Minors(movU2,X,bs);
900: put_line("The polynomial system for the second moving U-matrix :");
901: put_line(movcyc2);
902: put_line("Target system at t = 1 :");
903: put_line(Eval(movcyc2,Create(1.0),Number_of_Unknowns(movcyc2(1))));
904: put_line("System that must be solved :");
905: put_line(Expanded_Minors(L,X,bs));
906: end;
907: Matrix_Indeterminates.Clear_Symbols;
908: end Test_Pieri_Homotopies;
909:
910: procedure Main is
911:
912: m,p,n : natural;
913: ans : character;
914:
915: begin
916: new_line;
917: put_line("Interactive test for setting up of Pieri homotopies.");
918: new_line;
919: put("Give m, m+p = n, dimension of space : "); get(m);
920: put("Give p, number of entries in brackets : "); get(p);
921: n := m+p;
922: new_line;
923: put_line("MENU for testing the Pieri Homotopy Algorithm : ");
924: put_line(" 1. Construction of a Pieri tree. ");
925: put_line(" 2. Interactively test Pieri's condition on pairs of brackets.");
926: put_line(" 3. Set up equations for certain Schubert conditions.");
927: put_line(" 4. Set up the moving cycles in the Pieri homotopy algorithm.");
928: put_line(" 5. Test combinatorial root count and start deforming.");
929: put("Make your choice (1, 2, 3, 4 or 5) : "); get(ans);
930: case ans is
931: when '1' => Test_Pieri_Tree(n,p);
932: when '2' => Test_Pieri_Condition(n,p);
933: when '3' => Test_Pieri_Equations(n,p);
934: when '4' => Test_Pieri_Homotopies(n,p);
935: when '5' => Test_Root_Count(n,p);
936: when others => null;
937: end case;
938: end Main;
939:
940: begin
941: Main;
942: end ts_org_pieri;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ts_org_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