2121deaths = np .linspace (0 , 5 , num = 100 )
2222params = [(x , y ) for (x , y ) in itertools .product (births , deaths ) if x > y ]
2323
24+
2425def get_bd (r , a ):
2526 """Converts turnover and relative extinction to birth and death rates."""
2627 return - r / (a - 1 ), - a * r / (a - 1 )
2728
29+
2830def get_ra (b , d ):
2931 return (b - d , d / b )
3032
33+
3134def optim_bd_r (ages , sampling ):
3235 """Optimizes birth death using TreePar and R"""
3336 script = """cat(optim(c({birth}, {death}), function(v, ...) TreePar::LikConstant(v[1], v[2], ...), x = c({ages}), sampling = {sampling}, lower=c(.Machine$double.xmin,0), method = "L-BFGS-B")$par, " dum")"""
@@ -36,18 +39,21 @@ def optim_bd_r(ages, sampling):
3639 b , d , _ = output .split (None , 2 )
3740 return float (b ), float (d )
3841
42+
3943def optim_bd_grid (ages , sampling ):
4044 """Optimizes birth death using a grid search"""
4145 res = [lik_constant (x , sampling , ages ) for x in params ]
4246 return params [np .argmin (res )]
4347
48+
4449def update_multiplier_freq (q , d = 1.1 ):
4550 u = np .random .uniform (0 , 1 , 2 )
4651 l = 2 * log (d )
4752 m = np .exp (l * (u - 0.5 ))
4853 new_q = q * m
4954 return new_q
5055
56+
5157def optim_bd_mcmc (ages , sampling ):
5258 """Optimizes birth death using a cheap MCMC-like algorithm"""
5359 new_vec = [0 , 0 ]
@@ -66,6 +72,7 @@ def optim_bd_mcmc(ages, sampling):
6672 vec = new_vec
6773 return vec
6874
75+
6976def wrapped_lik_constant (x , sampling , ages ):
7077 return lik_constant (get_bd (* x ), sampling , ages )
7178
@@ -80,10 +87,12 @@ def optim_bd_scipy(ages, sampling):
8087 bounds = ((1e-6 , None ), (0 , 1 - 1e-6 ))
8188 return get_bd (* minimize (wrapped_lik_constant , (init_r , 0.0 ), args = (sampling , ages ), bounds = bounds , method = "TNC" )["x" ].tolist ())
8289
90+
8391def optim_bd (ages , sampling ):
8492 return optim_bd_scipy (ages , sampling )
8593
86- def optim_yule (ages ,sampling ):
94+
95+ def optim_yule (ages , sampling ):
8796 """Optimizes a Yule model using Scipy"""
8897 if max (ages ) < 0.000001 :
8998 init_r = 1e-3
@@ -93,6 +102,7 @@ def optim_yule(ages,sampling):
93102 bounds = ((1e-6 , None ), (0 , 0 ))
94103 return get_bd (* minimize (wrapped_lik_constant , (init_r , 0.0 ), args = (sampling , ages ), bounds = bounds , method = "TNC" )["x" ].tolist ())
95104
105+
96106def get_lik (vec , rho , x ):
97107 l = vec [0 ]
98108 m = vec [1 ]
@@ -102,19 +112,22 @@ def get_lik(vec, rho, x):
102112 lik3 = - (root + 1 ) * np .log (1 - p0 (x [0 ], l , m , rho ))
103113 return lik1 + lik2 + lik3
104114
115+
105116def p0_exact (t , l , m , rho ):
106117 t = D (t )
107118 l = D (l )
108119 m = D (m )
109120 rho = D (rho )
110121 return D (1 ) - rho * (l - m ) / (rho * l + (l * (D (1 ) - rho ) - m ) * (- (l - m ) * t ).exp ())
111122
123+
112124def p0 (t , l , m , rho ):
113125 try :
114126 return 1 - rho * (l - m ) / (rho * l + (l * (1 - rho ) - m ) * exp (- (l - m ) * t ))
115127 except FloatingPointError :
116128 return float (p0_exact (t , l , m , rho ))
117129
130+
118131def p1_exact (t , l , m , rho ):
119132 """Exact version of p1 using Decimal math."""
120133 t = D (t )
@@ -125,6 +138,7 @@ def p1_exact(t, l, m, rho):
125138 denom = (rho * l + (l * (1 - rho ) - m ) * (- (l - m ) * t ).exp ()) ** D (2 )
126139 return num / denom
127140
141+
128142def p1_orig (t , l , m , rho ):
129143 try :
130144 num = rho * (l - m ) ** 2 * np .exp (- (l - m ) * t )
@@ -133,6 +147,7 @@ def p1_orig(t, l, m, rho):
133147 except (OverflowError , FloatingPointError ):
134148 return float (p1_exact (t , l , m , rho ))
135149
150+
136151def p1 (t , l , m , rho ):
137152 # Optimized version of p1 using common subexpression elimination and strength reduction from
138153 # exponentiation to multiplication.
@@ -144,6 +159,7 @@ def p1(t, l, m, rho):
144159 except (OverflowError , FloatingPointError ):
145160 return float (p1_exact (t , l , m , rho ))
146161
162+
147163def intp1_exact (t , l , m ):
148164 """Exact version of intp1 using Decimal math."""
149165 l = D (l )
@@ -153,12 +169,14 @@ def intp1_exact(t, l, m):
153169 denom = (l - m * (- (l - m ) * t ).exp ())
154170 return num / denom
155171
172+
156173def intp1 (t , l , m ):
157174 try :
158- return (1 - exp (- (l - m ) * t ))/ (l - m * exp (- (l - m ) * t ))
175+ return (1 - exp (- (l - m ) * t )) / (l - m * exp (- (l - m ) * t ))
159176 except OverflowError :
160177 return float (intp1_exact (t , l , m ))
161178
179+
162180def lik_constant (vec , rho , t , root = 1 , survival = 1 , p1 = p1 ):
163181 """
164182 Calculates the likelihood of a constant-rate birth-death process, conditioned
@@ -210,6 +228,7 @@ def crown_capture_probability(n, k):
210228 return 0 # not technically correct but it works for our purposes
211229 return 1 - 2 * (n - k ) / ((n - 1 ) * (k + 1 ))
212230
231+
213232def get_monophyletic_node (tree , species ):
214233 """Returns the node or None that is the MRCA of the `species` in `tree`."""
215234 mrca = tree .mrca (taxon_labels = species )
@@ -218,6 +237,7 @@ def get_monophyletic_node(tree, species):
218237 if mrca and species .issuperset (get_tip_labels (mrca )):
219238 return mrca
220239
240+
221241def get_birth_death_rates (node , sampfrac , yule = False , include_root = False ):
222242 """
223243 Estimates the birth and death rates for the subtree descending from
@@ -229,18 +249,21 @@ def get_birth_death_rates(node, sampfrac, yule=False, include_root=False):
229249 else :
230250 return optim_bd (get_ages (node , include_root ), sampfrac )
231251
252+
232253def get_ages (node , include_root = False ):
233254 ages = [x .age for x in node .ageorder_iter (include_leaves = False , descending = True )]
234255 if include_root :
235256 ages += [node .age ]
236257 return ages
237258
259+
238260def get_tip_labels (tree_or_node ):
239261 try :
240262 return set ([x .taxon .label for x in tree_or_node .leaf_node_iter ()])
241263 except AttributeError :
242264 return set ([x .taxon .label for x in tree_or_node .leaf_iter ()])
243265
266+
244267def edge_iter (node , filter_fn = None ):
245268 """
246269 Iterates over the child edge of `node` and all its descendants.
@@ -253,6 +276,7 @@ def edge_iter(node, filter_fn=None):
253276 yield edge
254277 stack .extend (edge .head_node .child_edge_iter ())
255278
279+
256280def get_tree (path , namespace = None ):
257281 """
258282 Gets a DendroPy tree from a path and precalculate its node ages and bipartition bitmask.
@@ -262,19 +286,23 @@ def get_tree(path, namespace=None):
262286 tree .encode_bipartitions ()
263287 return tree
264288
289+
265290def is_binary (node ):
266291 """Is the subtree under `node` a fully bifurcating tree?"""
267292 for x in node .preorder_internal_node_iter ():
268293 if len (x .child_nodes ()) != 2 :
269294 return False
270295 return True
271296
297+
272298def get_short_branches (node ):
273299 for edge in edge_iter (node ):
274300 if edge .length <= 0.001 :
275301 yield edge
276302
277303# TODO: This could probably be optimized
304+
305+
278306def get_new_times (ages , birth , death , missing , told = None , tyoung = None ):
279307 """
280308 Simulates new speciation events in an incomplete phylogeny assuming a
0 commit comments