# Ticket #289: question.sing

File question.sing, 4.7 KB (added by , 12 years ago) |
---|

Line | |
---|---|

1 | // Clade distribution parameterizations |

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 | |

11 | LIB "primdec.lib"; |

12 | LIB "qhmoduli.lib"; |

13 | LIB "solve.lib"; |

14 | LIB "poly.lib"; |

15 | |

16 | system("--browser", "mac"); |

17 | option("redSB"); |

18 | |

19 | // 1 = wordy; 0 = quiet |

20 | int wordy = 0; |

21 | |

22 | // common ring for elimination ideals with parameters |

23 | ring 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)) |

28 | poly bal1 = 1-1/18*x*y^3*z-2/3*x-1/9*x*y^3; |

29 | poly bal2 = -1/18*x*y^3*z+1/3*x-1/9*x*y^3; |

30 | poly bal3 = 1/9*y*x*z+1/18*x*y^3*z; |

31 | poly bal4 = 1/9*y*x*z+1/18*x*y^3*z; |

32 | poly bal5 = -1/18*x*y^3*z+1/3*x-1/9*x*y^3; |

33 | poly bal6 = 1/9*y*x*z+1/18*x*y^3*z; |

34 | poly bal7 = 1/9*y*x*z+1/18*x*y^3*z; |

35 | poly bal8 = -1/9*y*x*z+1/18*x*y^3*z+2/9*y*z; |

36 | poly bal9 = -1/9*y*x*z+1/18*x*y^3*z+2/9*y*z; |

37 | poly bal10 = 1-2/3*z-1/18*y*x*z-1/9*y*z; |

38 | poly 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; |

39 | poly bal12 = -1/9*y*x*z-1/36*x*y^3*z+2/9*y*z; |

40 | poly bal13 = -1/9*y*x*z-1/36*x*y^3*z+2/9*y*z; |

41 | poly bal14 = 1/9*y*x*z-1/36*x*y^3*z; |

42 | poly bal15 = 1/9*y*x*z-1/36*x*y^3*z; |

43 | poly bal16 = -1/12*y*x*z-1/18*x*y^3*z+1/6*y*x+1/18*x*y^3; |

44 | poly bal17 = 1/9*y*x*z-1/36*x*y^3*z; |

45 | poly bal18 = 1/9*y*x*z-1/36*x*y^3*z; |

46 | poly bal19 = -1/12*y*x*z-1/18*x*y^3*z+1/6*y*x+1/18*x*y^3; |

47 | poly 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; |

48 | poly bal21 = 1/3*z-1/12*y*x*z+1/60*x*y^3*z-1/6*y*z; |

49 | poly bal22 = 1/3*z-1/12*y*x*z+1/60*x*y^3*z-1/6*y*z; |

50 | poly bal23 = 1/12*y*x*z+1/60*x*y^3*z+1/3*y-1/6*y*z-1/6*y*x; |

51 | poly bal24 = -1/12*y*x*z+1/60*x*y^3*z+1/6*y*x; |

52 | poly bal25 = -1/12*y*x*z+1/60*x*y^3*z+1/6*y*x; |

53 | |

54 | ideal 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 | |

62 | ideal 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 |

67 | list 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))"; |

75 | ideal Hbal = B; |

76 | Hbal = homog(Hbal,H); |

77 | int sze = size(Hbal); |

78 | |

79 | int m = Max(balDeg); |

80 | if (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];} |

87 | if (wordy==1) |

88 | { |

89 | "Those needing powers of H"; |

90 | for (int i=1; i<=sze; i=i+1) |

91 | Hcat { |

92 | if (m-deg(Hbal[i])) |

93 | { |

94 | i, " ", Hbal[i]; |

95 | } |

96 | } |

97 | } |

98 | |

99 | // homogenize polynomials |

100 | for (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 | |

108 | if (wordy==1){ |

109 | "After making all the same degree.", deg(Hps[1..sze]); |

110 | } |

111 | |

112 | poly gg1 = c1-Hbal[1]; |

113 | poly gg2 = c2-Hbal[2]; |

114 | poly gg3 = c3-Hbal[3]; |

115 | poly gg8 = c8-Hbal[8]; |

116 | poly gg10 = c10-Hbal[10]; |

117 | poly gg11 = c11-Hbal[11]; |

118 | poly gg12 = c12-Hbal[12]; |

119 | poly gg14 = c14-Hbal[14]; |

120 | poly gg16 = c16-Hbal[16]; |

121 | poly gg20 = c20-Hbal[20]; |

122 | poly gg21 = c21-Hbal[21]; |

123 | poly gg23 = c23-Hbal[23]; |

124 | poly gg24 = c24-Hbal[24]; |

125 | |

126 | ideal JHbal = gg1, gg2, gg3, gg8, gg10, gg11, gg12, gg14, gg16, gg20, gg21, gg23, gg24; |

127 | |

128 | // option(prot); |

129 | degBound=20; |

130 | " For the balanced tree: The degree bound is ", degBound; |

131 | |

132 | // JHbal; |

133 | ideal 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); |

146 | degBound=10; |

147 | "**** Resetting degBound ****"; |

148 | " For the balanced tree: The degree bound is ", degBound; |

149 | |

150 | // JHbal; |

151 | ideal 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); |

164 | degBound=5; |

165 | "**** Resetting degBound ****"; |

166 | " For the balanced tree: The degree bound is ", degBound; |

167 | |

168 | // JHbal; |

169 | ideal 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 | |

183 | string cmd; |

184 | cmd = read("question.sing"); |