# Ticket #289: question.sing

File question.sing, 4.7 KB (added by e.allman@…, 12 years ago)

sample Singular code to illustrate the problem

Line
2//
3// balanced tree ((ab)c,(de))
4// pseudo-cat ((ab)(de),c)
5// caterpillar ((((ab)c)d)e)
6//
7// modified from clade-dist-homog.sing on August 30, 2010
8//
9// In particular, the species trees were changed!
10
11LIB "primdec.lib";
12LIB "qhmoduli.lib";
13LIB "solve.lib";
14LIB "poly.lib";
15
16system("--browser", "mac");
17option("redSB");
18
19// 1 = wordy; 0 = quiet
20int wordy = 0;
21
22// common ring for elimination ideals with parameters
23ring r = 0,(x,y,z,H,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25), dp;
24
25// Parameterization for the three non-caterpillar maps
26
27// balanced ((ab)c,(de))
28poly bal1 = 1-1/18*x*y^3*z-2/3*x-1/9*x*y^3;
29poly bal2 = -1/18*x*y^3*z+1/3*x-1/9*x*y^3;
30poly bal3 = 1/9*y*x*z+1/18*x*y^3*z;
31poly bal4 = 1/9*y*x*z+1/18*x*y^3*z;
32poly bal5 = -1/18*x*y^3*z+1/3*x-1/9*x*y^3;
33poly bal6 = 1/9*y*x*z+1/18*x*y^3*z;
34poly bal7 = 1/9*y*x*z+1/18*x*y^3*z;
35poly bal8 = -1/9*y*x*z+1/18*x*y^3*z+2/9*y*z;
36poly bal9 = -1/9*y*x*z+1/18*x*y^3*z+2/9*y*z;
37poly bal10 = 1-2/3*z-1/18*y*x*z-1/9*y*z;
38poly bal11 = 1-1/18*y*x*z+1/12*x*y^3*z-2/3*y-1/9*y*z-1/3*y*x+1/6*x*y^3;
39poly bal12 = -1/9*y*x*z-1/36*x*y^3*z+2/9*y*z;
40poly bal13 = -1/9*y*x*z-1/36*x*y^3*z+2/9*y*z;
41poly bal14 = 1/9*y*x*z-1/36*x*y^3*z;
42poly bal15 = 1/9*y*x*z-1/36*x*y^3*z;
43poly bal16 = -1/12*y*x*z-1/18*x*y^3*z+1/6*y*x+1/18*x*y^3;
44poly bal17 = 1/9*y*x*z-1/36*x*y^3*z;
45poly bal18 = 1/9*y*x*z-1/36*x*y^3*z;
46poly bal19 = -1/12*y*x*z-1/18*x*y^3*z+1/6*y*x+1/18*x*y^3;
47poly bal20 = 1/12*y*x*z-1/18*x*y^3*z+1/3*y-1/6*y*z-1/6*y*x+1/18*x*y^3;
48poly bal21 = 1/3*z-1/12*y*x*z+1/60*x*y^3*z-1/6*y*z;
49poly bal22 = 1/3*z-1/12*y*x*z+1/60*x*y^3*z-1/6*y*z;
50poly bal23 = 1/12*y*x*z+1/60*x*y^3*z+1/3*y-1/6*y*z-1/6*y*x;
51poly bal24 = -1/12*y*x*z+1/60*x*y^3*z+1/6*y*x;
52poly bal25 = -1/12*y*x*z+1/60*x*y^3*z+1/6*y*x;
53
54ideal Bi = c1-bal1,c2-bal2,c3-bal3,c4-bal4,c5-bal5,c6-bal6,c7-bal7,c8-bal8,c9-bal9,c10-bal10,
55           c11-bal11,c12-bal12,c13-bal13,c14-bal14,c15-bal15,c16-bal16,c17-bal17,c18-bal18,c19-bal19,c20-bal20,
56           c21-bal21,c22-bal22,c23-bal23,c24-bal24,c25-bal25;
57
58" ";
59"Working to homogenize parameterizations ....";
60
61
62ideal B = bal1,bal2,bal3,bal4,bal5,bal6,bal7,bal8,bal9,bal10,
63          bal11,bal12,bal13,bal14,bal15,bal16,bal17,bal18,bal19,bal20,
64          bal21,bal22,bal23,bal24,bal25;
65
66// get homogeneous versions
67list balDeg = deg(bal1), deg(bal2), deg(bal3), deg(bal8), deg(bal10), deg(bal11), deg(bal12), deg(bal14),
68             deg(bal16), deg(bal20), deg(bal21), deg(bal23), deg(bal24);
69"   balDeg";
70"    ",balDeg[1..size(balDeg)];
71
72
73" ";
74"Working on 5-taxon balanced  ((ab)c,(de))";
75ideal Hbal = B;
76Hbal = homog(Hbal,H);
77int sze = size(Hbal);
78
79int m = Max(balDeg);
80if (wordy==1){
81  "The maximum degree before homogenization is", m;
82}
83
84//"Before making all the same degree...";
85//"deg,  poly";
86//for (int i=1; i<=sze; i=i+1){deg(Hbal[i]), "   ", Hbal[i];}
87if (wordy==1)
88  {
89  "Those needing powers of H";
90  for (int i=1; i<=sze; i=i+1)
91Hcat    {
92     if (m-deg(Hbal[i]))
93       {
94        i, "   ", Hbal[i];
95       }
96     }
97  }
98
99// homogenize polynomials
100for (int i=1; i<=sze; i=i+1)
101  {
102   if (m-deg(Hbal[i]))
103     {
104      Hbal[i] = Hbal[i]*H^(m-deg(Hbal[i]));
105     }
106  }
107
108if (wordy==1){
109  "After making all the same degree.", deg(Hps[1..sze]);
110}
111
112poly gg1 = c1-Hbal[1];
113poly gg2 = c2-Hbal[2];
114poly gg3 = c3-Hbal[3];
115poly gg8 = c8-Hbal[8];
116poly gg10 = c10-Hbal[10];
117poly gg11 = c11-Hbal[11];
118poly gg12 = c12-Hbal[12];
119poly gg14 = c14-Hbal[14];
120poly gg16 = c16-Hbal[16];
121poly gg20 = c20-Hbal[20];
122poly gg21 = c21-Hbal[21];
123poly gg23 = c23-Hbal[23];
124poly gg24 = c24-Hbal[24];
125
126ideal JHbal = gg1, gg2, gg3, gg8, gg10, gg11, gg12, gg14, gg16, gg20, gg21, gg23, gg24;
127
128// option(prot);
129degBound=20;
130"  For the balanced tree:  The degree bound is ", degBound;
131
132// JHbal;
133ideal stdJHbal = elim1(JHbal,x*y*z*H);
134"  The size of stdJHbal is", size(stdJHbal);
135
136" ";
137"Summary:  ";
138"  The size of stdJHbal is", size(stdJHbal);
139
140"The degrees of the std's are ";
141"  ", deg(stdJHbal[1..size(stdJHbal)]);
142
143
144" ";
145// option(prot);
146degBound=10;
147"**** Resetting degBound ****";
148"  For the balanced tree:  The degree bound is ", degBound;
149
150// JHbal;
151ideal stdJHbal = elim1(JHbal,x*y*z*H);
152"  The size of stdJHbal is", size(stdJHbal);
153
154" ";
155"Summary:  ";
156"  The size of stdJHbal is", size(stdJHbal);
157
158"The degrees of the std's are ";
159"  ", deg(stdJHbal[1..size(stdJHbal)]);
160
161
162" ";
163// option(prot);
164degBound=5;
165"**** Resetting degBound ****";
166"  For the balanced tree:  The degree bound is ", degBound;
167
168// JHbal;
169ideal stdJHbal = elim1(JHbal,x*y*z*H);
170"  The size of stdJHbal is", size(stdJHbal);
171
172" ";
173"Summary:  ";
174"  The size of stdJHbal is", size(stdJHbal);
175
176"The degrees of the std's are ";
177"  ", deg(stdJHbal[1..size(stdJHbal)]);
178
179
180" ";
181"THE END";
182
183string cmd;