Skip to content

Commit 9a5ce92

Browse files
author
Alexandra (Tally) Balaban
committed
Changed ForwardProbMIG name to ForwardProbA. Now ForwardProbMIG calculates eq. 30.
in cjumpchain one of the comments for SampleFromIS is wrong. For the next level, (i.e., level 3) the initial character state assignment will be: {1:0, 4:1, 5:1} (since 2 & 3 were in state 1, their parent 5 also will be in state 1). {1:0, 4:1, 5:1} as the initial character state assignment I changed this to For the next level, (i.e., level 3) the initial character state assignment will be: {1:0, 4:1, 5:1} (since 2 & 3 were in state 0, their parent 5 also will be in state 0). {1:0, 4:1, 5:0} as the initial character state assignment Re-wrote def GetInfoForForwardProbMIG(A,level): into a multiple functions compatible with A's tree structure: GetNumberOfLevels(A): GetBranchLength(A,level): GetNumberOfMigrationEvents(A,level): GetQ_bbAndQ_bt(A,level): I need help with this one Different inputs, but same function: ProbKMigrationInL(level,l,k,q_bb,q_bt,sigma) PhiJK(j,k,(r_t,q_bb,q_b_arrow_t,q_bt),sigma) Phi(i,(r_t,q_bb,q_b_arrow_t,q_bt),sigma): Implemented GetH(X,i,A,sigma): Gives the rate of event X in the ith level Implemented GetHBar(i,A,sigma): gives the sum of the rates of all possible speciation events in the level i Implemented GetYi(A,i): returns the speciation event that ends level i Theoretically, ForwardProbA should be working now, (which may mean all of is_likelihood is implemented and working) Added to MakeTransitionMatrixForLevel documentation Didn't like PickNextStateofChain implementation. May work, but not intuitive to me, so re-wrote. I think there is an error in the last lines of code of ChooseLineageAndUpdateDelta. If lineages are already in state 1, they shouldn't be included in the random_list. As code was: for i in range(len(current_delta)): if (i != next_coalescing_lineages[0].index and i != next_coalescing_lineages[1].index): random_list.append(i) random.shuffle(random_list) lineage_to_migrate = random_list[0] if (current_delta[lineage_to_migrate] == 1): #the if statement is really redundant current_delta[lineage_to_migrate] = 0 return random_list[0] So I changed the line to if (i != next_coalescing_lineages[0].index and i != next_coalescing_lineages[1].index and current_delta[i]!=1): I also added an assert statement to ChooseLineageAndUpdateDelta changed delta structure to reflect the fact that delta is appended to event history and should be consistent with the event_history description in is_classes. So I changed current_delta from " a list of of character states (boreal or tropical)" to " a map of lineage indices to character states". Made modifications in the code to reflect the fact that delta is a map, not a list. added all_delta, a map of all lineages at all levels to their states. As it is updated as A is created, it will eventually be a map of all lineages to their states right after they came into existence MoveToNextLevel finished Edited PrepareTree in SampleFromIS if(index_of_next_state==len(transition_matrix_for_the_level): whether_in_the_same_level == False #whether_in_the_same_level = WhetherInTheSameLevel(state_of_cond_jump_chain,next_state_of_cond_jump_chain) Since next_state_of_cond_jump_chain = index_in_transition_matrix_to_state[len(index_in_transition_matrix_to_state)] is undefined. A kappa event has an index of len(index_in_transition_matrix_to_state) which is outside of the range of index_in_transition_matrix_to_state. created a Test function to get the inputs for SampleFromIS SampleFromIS when called from the Test function in cjumpchain is returning a value and is filling out the event_histories for each level I haven't gotten around to debugging is_likelihood yet.
1 parent a83b9c0 commit 9a5ce92

File tree

6 files changed

+515
-335
lines changed

6 files changed

+515
-335
lines changed

‎cjumpchain.py‎

Lines changed: 262 additions & 197 deletions
Large diffs are not rendered by default.

‎is_likelihood.py‎

Lines changed: 186 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -13,39 +13,69 @@
1313
import math
1414

1515

16-
def GetInfoForForwardProbMIG(A,level):
17-
"""
18-
A helper function to get inputs for ProbMIG
19-
20-
Input Parameters
21-
----------------
22-
@param A the tuple (theta,BRL(theta),MIG(theta))
23-
@param level the level in the tree
16+
def GetNumberOfLevels(A):
17+
return len(A.levels)
18+
19+
def GetLevelLength(A,level):
20+
current_level=A.levels[level]
21+
branch_length=current_level.begin_time-current_level.end_time
22+
return branch_length
23+
24+
def GetNumberOfMigrationEvents(A,level):
25+
current_level=A.levels[level]
26+
current_event_history=current_level.event_history
27+
return len(event_history)-1
28+
29+
def GetQ_bbAndQ_bt(A,level):
30+
current_level=A.levels[level]
31+
current_event_history=current_level.event_history
32+
current_lineages=current_level.lineages
33+
q_bb=0
34+
q_bt=0
35+
for current_lineage in level_lineages:
36+
current_lineage_index=current_lineage.index
37+
#see if the lineages is a boreal lineage and not a migrating one
38+
if(current_event_history[0][current_lineage_index]==0 and (current_lineage_index not in current_event_history)):
39+
children=current_lineage.children
40+
child_1=children[0]
41+
child_2=children[1]
42+
child_1_index=child_1.index
43+
child_2_index=child_2.index
44+
A_levels=A.levels
45+
child_1_level_index=-1
46+
child_2_level_index=-1
47+
found="FALSE"
48+
#find the levels at which the two children lineages/nodes first appear
49+
for i in range(1,GetNumberOfLevels(A)+1):
50+
current_event_history=A_levels[i].event_history
51+
if(current_event_history[0].has_key(child_1_index)):
52+
child_1_level_index=i
53+
found="TRUE"
54+
if(found=="TRUE"):
55+
break
56+
57+
for i in range(1,GetNumberOfLevels(A)+1):
58+
current_event_history=A_levels[i].event_history
59+
if(current_event_history[0].has_key(child_2_index)):
60+
child_2_level_index=i
61+
found="TRUE"
62+
if(found=="TRUE"):
63+
break
64+
65+
child_1_level=A.levels[child_1_level_index]
66+
child_2_level=A.levels[child_2_level_index]
67+
child_1_state=child_1_level.event_history[0][child_1_index]
68+
child_2_state=child_2_level.event_history[0][child_2_index]
69+
if(child_1_state==1 or child_2_state=1):
70+
q_bt+=1
71+
72+
else:
73+
q_bb+=1
74+
75+
#end if statement making sure lineages are both boreal and non-migrating
76+
#end for-loop looping through all lineages in the specified level
77+
return(q_bt, q_bb)
2478

25-
Return Values
26-
-------------
27-
k the number of migration (SBArrowT) events that will occur within the level
28-
states_within_level a k+1 length array of all the (r_t,q_bb,q_b_arrow_t,q_bt) states within the level,
29-
where states_within_level[i] is the state right after the ith migration
30-
l the time span of level "level"
31-
n the number of levels in A
32-
33-
Details
34-
------
35-
There are k+1 states within a level since there are k migration events within a level, and at each migration
36-
the values (r_t,q_bb,q_b_arrow_t,q_bt) change.
37-
states_within_level[i]= the tuple (r_t,q_bb,q_b_arrow_t,q_bt) representing the state after the i-th migration
38-
where i ranges from 0 to k.
39-
"""
40-
##right now I am returning one value in order to test out other functions
41-
##here the level is 6+5+4+3=18
42-
k=2
43-
l=100
44-
states_within_level=range(3)
45-
states_within_level[0]=(0,5,4,3)
46-
states_within_level[1]=(1,5,3,3)
47-
states_within_level[2]=(2,5,2,3)
48-
return (l,k,states_within_level,1)
4979

5080
def ForwardProbMIG(A,sigma):
5181
"""
@@ -56,43 +86,140 @@ def ForwardProbMIG(A,sigma):
5686
@param A the tuple (theta,BRL(theta),MIG(theta))
5787
@param level the level in the tree
5888
"""
59-
# level is 1 here for no particular reason, just need n
60-
(l,k,states_within_level,n)= GetInfoForForwardProbMIG(A,1)
61-
mult=1
89+
# level is 1 here for no particular reason, just need n's value which is independent of level
90+
n=GetNumberOfLevels(A)
6291
for level in range(1,n+1):
63-
(l,k,states_within_level,n)= GetInfoForForwardProbMIG(A,level)
64-
mult=mult*ProbKMigrationInL(l,k,states_within_level,sigma)
92+
l=GetLevelLength(A,level)
93+
k=GetNumberOfMigrationEvents(A,level)
94+
(q_bb,q_bt)=GetQ_bbAndQ_bt(A,level)
95+
mult=mult*ProbKMigrationInL(level,l,k,q_bb,q_bt,sigma)
6596

6697
return mult
98+
99+
def GetHBar(i,A,sigma):
100+
"""
101+
gives the sum of the rates of all possible speciation events in the level i
102+
103+
Input Parameters
104+
----------------
105+
i the level
106+
A the tree G plus character state assignments, of class Tree
107+
sigma see ForwardRateSTT
108+
Return Value
109+
-----------
110+
Sum of GetH(X,i,A,sigma) for X in {'s_tt,s_bb,s_bt'}
111+
"""
112+
return GetH('s_tt',i,A,sigma)+GetH('s_bb',i,A,sigma)+GetH('s_bt',i,A,sigma)
113+
114+
def GetYi(A,i):
115+
"""
116+
returns the speciation event that ends level i
117+
118+
Input
119+
-----
120+
@param i the level
121+
@param A The tree G plus character assignments, of class Tree
122+
123+
Return value
124+
------------
125+
the speciation event that ends level i in Tree A, either s_tt, s_bb, or s_bt
126+
"""
127+
level_i_event_history=A.levels[i].eventhistory
128+
next_level_event_history=A.levels[i+1].eventhistory
129+
set_lineage_indices_level_i=set(level_i_event_history[0].keys())
130+
set_lineage_indices_next_level=set(next_level_event_history[0].keys())
131+
set_new_lineages_indices=set_lineage_indices_next_level.difference(set_lineage_indices_level_i)
132+
set_bifurcated_lineage_index=set_lineage_indices_level_i.difference(set_lineage_indices_next_level)
133+
new_lineages_indices=[set_new_lineages.pop(),set_new_lineages.pop()]
134+
bifurcated_lineage_index=[set_bifurcated_lineage_index.pop()]
135+
136+
character_bifurcated=level_i_event_history[0][bifurcated_lineage_index]
137+
character_new_1=next_level_event_history[0][new_lineages_indices[0]]
138+
character_new_2=next_level_event_history[0][new_lineages_indices[1]]
139+
140+
if(character_bifurcated==1):
141+
return 's_tt'
142+
if(character_new_1==1 or character_new_2==1):
143+
return 's_bt'
144+
else:
145+
return 's_bb'
146+
return ()
147+
148+
return ()
149+
def GetH(X,i,A,sigma):
150+
"""
151+
Gives the rate of event X in the ith level
67152
153+
Input Parameters
154+
----------------
155+
@param X either 's_tt', 's_bb' or 's_bt'
156+
@param i the level
157+
@param A the tree G plus character state assignments, of class Tree
158+
@param sigma see ForwardRateSTT
68159
69-
def ProbKMigrationInL(l,k,states_within_level,sigma):
160+
Return Value
161+
------------
162+
the rate of event X in the ith level, where X is either STT, SBB or SBT
70163
"""
71-
Returns the forward probability of the k migrations that occur in the level lev. Defined in equation 29.
164+
(q_bb,q_bt)=GetQ_bbAndQ_bt(A,i)
165+
q_b_arrow_t=len(A.levels[i].event_history)-1
166+
r_t=i-q_b_arrow_t-q_bb-q_bt
167+
if(X=='s_tt'):
168+
return ForwardRateSTT((r_t,q_bb,q_b_arrow_t,q_bt),sigma)
169+
if(X=='s_bb'):
170+
return ForwardRateSBB((r_t,q_bb,q_b_arrow_t,q_bt),sigma)
171+
if(X=='s_bt'):
172+
return ForwardRateSBT((r_t,q_bb,q_b_arrow_t,q_bt),sigma)
173+
return ()
174+
175+
def ForwardProbThetaAndBRLGivenMIG(A,sigma)
176+
"""
177+
"""
178+
n=GetNumberOfLevels(A)
179+
for i in range(1,n):
180+
l=GetLevelLength(A,i)
181+
hbar_i=GetHBar(i,A,sigma)
182+
Yi=GetYi(A,i)
183+
mult=mult*hbar_i*math.exp(-hbar_i*l_i)*GetH(Y_i,i,A,sigma)/hbar_i
184+
185+
l_n=GetLevelLength(A,n)
186+
hbar_n=GetHBar(n,A,sigma)
187+
mult=mult*hbar_n*math.exp(-hbar_n*l_n)
188+
return mult
189+
190+
def ForwardProbA(A,sigma)
191+
"""
192+
"""
193+
return ForwardProbThetaAndBRLGivenMIG(A,sigma)*ForwardProbMIG(A,sigma)
194+
195+
196+
197+
def ProbKMigrationInL(level,l,k,q_bb,q_bt,sigma):
198+
"""
199+
Returns the forward probability of the k migrations that occur in the level. Defined in equation 29.
72200
Input Parameters
73201
---------------
202+
@param level the level number
74203
@param l the duration of the level
75204
76205
@param k the number of migrations in the level
77206
78-
@param states_within_level a k+1 length array of all the (r_t,q_bb,q_b_arrow_t,q_bt) states within the level,
79-
where states_within_level[i] is the state right after the ith migration
207+
@param q_bb the number of lineages in the level whose next event is a SBB event
208+
@param q_bt the number of lineages in the level whose next event is a SBT event
80209
81210
@param sigma sigma =(lambda, b, B, alpha, mu), see ForwarRateSTT for more details
82211
"""
83212
coef=1
213+
q_b_arrow_t=len(A.levels[i].event_history)-1
214+
r_t=level-q_b_arrow_t-q_bb-q_bt
84215
for i in range(1,k+1):
85-
#states_within_level[i-1] is the state right before the ith migration within the level
86-
coef=coef*Phi(states_within_level[i-1],sigma)
87-
216+
coef=coef*Phi(i,(r_t,q_bb,q_b_arrow_t,q_bt),sigma)
88217
sum=0
89-
#if the input to PhiJK is (j,k,...) the input to Phi, a term in PhiJK, is j-1
90218
for j in range(1,k+1):
91-
sum=sum+PhiJK(j,k,states_within_level,sigma)*math.exp(-Phi(states_within_level[j-1],sigma)*l)
92-
219+
sum=sum+PhiJK(j,k,(r_t,q_bb,q_b_arrow_t,q_bt),sigma)*math.exp(-Phi(j,(r_t,q_bb,q_b_arrow_t,q_bt),sigma)*l)
93220
return coef*sum
94221

95-
def PhiJK(j,k,states_within_level,sigma):
222+
def PhiJK(j,k,(r_t,q_bb,q_b_arrow_t,q_bt),sigma):
96223
"""
97224
PhiJK is defined in Equation 29
98225
This is a helper function to ProbKMigrationInL
@@ -101,56 +228,37 @@ def PhiJK(j,k,states_within_level,sigma):
101228
----------------
102229
@param j the input to Phi will be j-1
103230
@param k The number of migrations within the level
104-
@param states_within_level a k+1 length array of all the (r_t,q_bb,q_b_arrow_t,q_bt) states within the level,
105-
where states_within_level[i] is the state right after the ith migration
231+
@param (r_t,q_bb,q_b_arrow_t,q_bt) the (r_t,q_bb,q_b_arrow_t,q_bt) tuple before the first migration within the level
232+
See ForwardRateSTT for more details
106233
@param sigma sigma =(lambda, b, B, alpha, mu), where
107234
lambda is the turnover rate in state 0 (the boreal region),
108235
b = the effective migration rate from state 0 to 1 (boreal to tropical),
109236
B = the total number of species in state 0 (the boreal region)
110237
alpha = the per lineage speciation rate in state 1 (the neotropical region), and
111238
mu = the per lineage extinction rate in state 1 (the neotropical region).
112-
Details
113-
-------
114-
#states_within_level[i-1] is the state right before the ith migration within the level
115239
"""
116240
product=1
117-
#print("j is" + str(j))
118241
if(j==1):
119242
product=1
120-
#print("j equals 1")
121-
122243
else:
123244
for i in range(1,j):
124-
current_state_i=states_within_level[i-1]
125-
current_state_j=states_within_level[j-1]
126-
product=product*(Phi(current_state_i,sigma)-Phi(current_state_j,sigma))
127-
#print("in first")
128-
#print("range "+str(1)+" "+str(k))
129-
#print("i is " +str(i))
130-
#print(product)
131-
132-
#if(j==k):
133-
#print("j equals k")
245+
product=product*(Phi(i,(r_t,q_bb,q_b_arrow_t,q_bt),sigma)-Phi(j,(r_t,q_bb,q_b_arrow_t,q_bt),sigma))
246+
134247
if(j!=k):
135248
for i in range(j+1,k+1):
136-
current_state_i=states_within_level[i-1]
137-
current_state_j=states_within_level[j-1]
138-
product=product*(Phi(current_state_i,sigma)-Phi(current_state_j,sigma))
139-
#print("in second")
140-
#print("range "+str(j+1)+" "+str(k))
141-
#print("i is " +str(i))
142-
#print(product)
249+
product=product*(Phi(i,(r_t,q_bb,q_b_arrow_t,q_bt),sigma)-Phi(j,(r_t,q_bb,q_b_arrow_t,q_bt),sigma))
143250

144251
return pow(product,-1)
145252

146-
def Phi(current_state,sigma):
253+
def Phi(i,(r_t,q_bb,q_b_arrow_t,q_bt),sigma):
147254
"""
148255
Phi(i) is the forward rate of SBarrowT right before the i-th migration
149256
150257
Input Parameters
151258
----------------
152-
@param current_state the current (r_t,q_bb,q_b_arrow_t,q_bt) tuple before the ith migration within the level
153-
See ForwardRateSTT for more details
259+
@i represents the ith migration
260+
@param (r_t,q_bb,q_b_arrow_t,q_bt) the (r_t,q_bb,q_b_arrow_t,q_bt) tuple before the first migration within the level
261+
See ForwardRateSTT for more details
154262
155263
@param sigma sigma =(lambda, b, B, alpha, mu), See ForwardRateSTT for more details
156264
@@ -159,7 +267,7 @@ def Phi(current_state,sigma):
159267
ForwardRateSBArrowT(current_state,sigma), the forward rate of SBArrowT right before the ith migration within
160268
the level
161269
"""
162-
return ForwardRateSBArrowT(current_state,sigma)
270+
return ForwardRateSBArrowT((r_t+i-1,q_bb,q_b_arrow_t-i+1,q_bt),sigma)
163271

164272

165273
def ForwardRateSTT(current_state,sigma,num_sum=20):
@@ -364,7 +472,7 @@ def LikelihoodOfParameters(G, delta, sigma):
364472
# used in place of G.
365473
# calculate forward probability of A given sigma, Pr(A | sigma)
366474
# according to the procedure described in Section 5.4.
367-
forward_probability_of_A_given_sigma = ForwardProbMIG(A, sigma)
475+
forward_probability_of_A_given_sigma = ForwardProbA(A,sigma)
368476

369477
# It is guaranteed that density_of_A > 0
370478
importance_sampling_weight = forward_probability_of_A_given_sigma/density_of_A

‎is_main.py‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,8 @@
4545
# standard python imports
4646
import os
4747
import sys
48-
import re
48+
import re
49+
import levels
4950

5051
# Ganesh imports
5152
import parse
@@ -82,7 +83,7 @@
8283
taxa_block = False
8384
states_block = False
8485

85-
TREE_FILE = open('validStrings.nex', "r")
86+
TREE_FILE = open('onetree.nex', "r")
8687

8788
for line in TREE_FILE.readlines():
8889
# skip if line contains just whitespaces
@@ -130,11 +131,11 @@
130131
tree_name = tree_tuple[0].strip()
131132
tree_string = tree_tuple[1].strip()
132133
parsed_tree = parse.Read(tree_name, taxon_table, state_table, tree_string)
133-
#parsed_tree.name = tree_name
134+
parsed_tree.name = tree_name
134135
#tree_traverse.PostOrderTraverse(parsed_tree.root, tree_traverse.CalculateNodeAges)
135136
#tree_traverse.PrintLabel(parsed_tree.root)
136137
#tree_traverse.PostOrderTraverse(parsed_tree.root, tree_traverse.PrintSelfParentLabel)
137-
#levels.DemarcateLevels(parsed_tree)
138+
levels.DemarcateLevels(parsed_tree)
138139
#levels.PrintAllLevels(parsed_tree)
139140
#tree_traverse.TreeLeaves(parsed_tree)
140141

‎levels.py‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -132,4 +132,6 @@ def DemarcateLevels(tree):
132132
for child in curr_node.children:
133133
new_level.lineages.remove(child)
134134
new_level.lineages.append(curr_node)
135-
level = new_level
135+
level = new_level
136+
137+

0 commit comments

Comments
 (0)