1 /*
2
3 PhyML: a program that computes maximum likelihood phylogenies from
4 DNA or AA homologous sequences.
5
6 Copyright (C) Stephane Guindon. Oct 2003 onward.
7
8 All parts of the source except where indicated are distributed under
9 the GNU public licence. See http://www.opensource.org for details.
10
11 */
12
13 #include "rw.h"
14
15 //////////////////////////////////////////////////////////////
16 //////////////////////////////////////////////////////////////
17
RW_Lk(t_tree * tree)18 phydbl RW_Lk(t_tree *tree)
19 {
20 phydbl fwd_lk,ic_lk,coal_lk;
21
22 assert(tree->mmod->id == RW);
23
24 #ifdef PHYREX
25 PHYREX_Update_Lindisk_List(tree);
26 PHYREX_Ldsk_To_Tree(tree);
27 if(PHYREX_Total_Number_Of_Intervals(tree) > tree->mmod->max_num_of_intervals) return UNLIKELY;
28 #endif
29
30 fwd_lk = RW_Forward_Lk(tree);
31 coal_lk = TIMES_Lk_Coalescent(tree);
32 ic_lk = RW_Independent_Contrasts(tree);
33
34 // fwd_lk may be larger than ic_lk since the marginal density of lineage location (on an
35 // habitat of infinite area) is 1/infinity
36 /* if(fwd_lk > ic_lk) */
37 /* { */
38 /* PhyML_Fprintf(stderr,"\n. fwd_lk = %f ic_lk = %f sigsq = %f",fwd_lk,ic_lk,tree->mmod->sigsq); */
39 /* PhyML_Fprintf(stderr,"\n. move: %s",tree->mcmc->move_name[tree->mcmc->move_idx]); */
40 /* assert(fwd_lk < ic_lk); */
41 /* } */
42
43 tree->mmod->c_lnL = fwd_lk - ic_lk + coal_lk;
44 /* tree->mmod->c_lnL = fwd_lk + coal_lk; */
45 /* tree->mmod->c_lnL = coal_lk; */
46
47 return(tree->mmod->c_lnL);
48 }
49
50 //////////////////////////////////////////////////////////////
51 //////////////////////////////////////////////////////////////
52 /* Returns probability density of observed lineage location as
53 defined in Felsenstein's (1973) independent contrasts article.
54 */
RW_Independent_Contrasts(t_tree * tree)55 phydbl RW_Independent_Contrasts(t_tree *tree)
56 {
57 phydbl lnP;
58 t_node *n1,*n2;
59 phydbl v1,v2,u,s;
60
61 lnP = 0.0;
62
63 for(int i=0;i<tree->mmod->n_dim;++i)
64 {
65 RW_Init_Contrasts(i,tree);
66 RW_Independent_Contrasts_Post(tree->n_root,tree->n_root->v[1],&lnP,tree);
67 RW_Independent_Contrasts_Post(tree->n_root,tree->n_root->v[2],&lnP,tree);
68
69 n1 = tree->n_root->v[1];
70 n2 = tree->n_root->v[2];
71
72 // v1 and v2 in Equation (7) of Felsenstein's article
73 v1 = fabs(tree->ctrst->tprime[n1->num]-tree->times->nd_t[tree->n_root->num]);
74 v2 = fabs(tree->ctrst->tprime[n2->num]-tree->times->nd_t[tree->n_root->num]);
75
76 // x6
77 tree->ctrst->x[tree->n_root->num] =
78 (v2/(v2+v1))*tree->ctrst->x[n1->num] +
79 (v1/(v2+v1))*tree->ctrst->x[n2->num] ;
80
81 // u1
82 u = tree->ctrst->x[n1->num] - tree->ctrst->x[n2->num];
83
84 // t6'
85 tree->ctrst->tprime[tree->n_root->num] = tree->times->nd_t[tree->n_root->num] + (v1*v2)/(v1+v2);
86
87 s = tree->mmod->sigsq;
88 lnP -= log(sqrt(2.*PI*s*(v1+v2)));
89 lnP -= (1./(2.*s*(v1+v2)))*u*u;
90 }
91
92 return(lnP);
93 }
94
95 //////////////////////////////////////////////////////////////
96 //////////////////////////////////////////////////////////////
97
RW_Independent_Contrasts_Post(t_node * a,t_node * d,phydbl * lnP,t_tree * tree)98 void RW_Independent_Contrasts_Post(t_node *a, t_node *d, phydbl *lnP, t_tree *tree)
99 {
100 if(d->tax == YES) return;
101 else
102 {
103 t_node *n1,*n2;
104 phydbl v1,v2,u,s;
105
106 n1 = n2 = NULL;
107 for(int i=0;i<3;++i)
108 {
109 if(d->v[i] != a && d->b[i] != tree->e_root)
110 {
111 if(!n1) n1 = d->v[i];
112 else n2 = d->v[i];
113 RW_Independent_Contrasts_Post(d,d->v[i],lnP,tree);
114 }
115 }
116
117 // v1 and v2 in Equation (7) of Felsenstein's article
118 v1 = fabs(tree->ctrst->tprime[n1->num]-tree->times->nd_t[d->num]);
119 v2 = fabs(tree->ctrst->tprime[n2->num]-tree->times->nd_t[d->num]);
120
121 // x6
122 tree->ctrst->x[d->num] =
123 (v2/(v2+v1))*tree->ctrst->x[n1->num] +
124 (v1/(v2+v1))*tree->ctrst->x[n2->num] ;
125
126 // u1
127 u = tree->ctrst->x[n1->num] - tree->ctrst->x[n2->num];
128
129 // t6'
130 tree->ctrst->tprime[d->num] = tree->times->nd_t[d->num] + (v1*v2)/(v1+v2);
131
132 s = tree->mmod->sigsq;
133 *lnP -= log(sqrt(2.*PI*s*(v1+v2)));
134 *lnP -= (1./(2.*s*(v1+v2)))*u*u;
135 }
136 }
137
138 //////////////////////////////////////////////////////////////
139 //////////////////////////////////////////////////////////////
140
RW_Forward_Lk(t_tree * tree)141 phydbl RW_Forward_Lk(t_tree *tree)
142 {
143 phydbl lnP,la,ld,sd;
144 int i,err;
145
146 lnP = 0.0;
147 RW_Forward_Lk_Pre(tree->n_root,tree->n_root->v[1],&lnP,tree);
148 RW_Forward_Lk_Pre(tree->n_root,tree->n_root->v[2],&lnP,tree);
149
150
151 for(i=0;i<tree->mmod->n_dim;++i)
152 {
153 ld = tree->n_root->v[1]->ldsk->coord->lonlat[i];
154 la = tree->n_root->ldsk->coord->lonlat[i];
155 sd = sqrt(tree->mmod->sigsq*fabs(tree->times->nd_t[tree->n_root->v[1]->num]-tree->times->nd_t[tree->n_root->num]));
156
157 /* lnP += Log_Dnorm_Trunc(ld, */
158 /* la, */
159 /* sd, */
160 /* tree->mmod->lim_do->lonlat[i], */
161 /* tree->mmod->lim_up->lonlat[i],&err); */
162 lnP += Log_Dnorm(ld,la,sd,&err);
163
164
165 ld = tree->n_root->v[2]->ldsk->coord->lonlat[i];
166 la = tree->n_root->ldsk->coord->lonlat[i];
167 sd = sqrt(tree->mmod->sigsq*fabs(tree->times->nd_t[tree->n_root->v[2]->num]-tree->times->nd_t[tree->n_root->num]));
168
169 /* lnP += Log_Dnorm_Trunc(ld, */
170 /* la, */
171 /* sd, */
172 /* tree->mmod->lim_do->lonlat[i], */
173 /* tree->mmod->lim_up->lonlat[i],&err); */
174 lnP += Log_Dnorm(ld,la,sd,&err);
175 }
176
177 for(i=0;i<tree->mmod->n_dim;++i) lnP -= log(tree->mmod->lim_up->lonlat[i]-tree->mmod->lim_do->lonlat[i]);
178 return(lnP);
179 }
180
181 //////////////////////////////////////////////////////////////
182 //////////////////////////////////////////////////////////////
183
RW_Forward_Lk_Pre(t_node * a,t_node * d,phydbl * lnP,t_tree * tree)184 void RW_Forward_Lk_Pre(t_node *a, t_node *d, phydbl *lnP, t_tree *tree)
185 {
186 int i,err;
187 phydbl ld,la,sd;
188
189 sd = sqrt(tree->mmod->sigsq*fabs(tree->times->nd_t[d->num]-tree->times->nd_t[a->num]));
190
191 for(i=0;i<tree->mmod->n_dim;++i)
192 {
193 ld = d->ldsk->coord->lonlat[i];
194 la = a->ldsk->coord->lonlat[i];
195
196 /* *lnP += Log_Dnorm_Trunc(ld, */
197 /* la, */
198 /* sd, */
199 /* tree->mmod->lim_do->lonlat[i], */
200 /* tree->mmod->lim_up->lonlat[i],&err); */
201 *lnP += Log_Dnorm(ld,la,sd,&err);
202 }
203
204 if(d->tax) return;
205 else
206 {
207 for(i=0;i<3;++i)
208 {
209 if(d->v[i] != a && d->b[i] != tree->e_root)
210 {
211 RW_Forward_Lk_Pre(d,d->v[i],lnP,tree);
212 }
213 }
214 }
215 }
216
217 //////////////////////////////////////////////////////////////
218 //////////////////////////////////////////////////////////////
219 //////////////////////////////////////////////////////////////
220 //////////////////////////////////////////////////////////////
221 //////////////////////////////////////////////////////////////
222 //////////////////////////////////////////////////////////////
223 //////////////////////////////////////////////////////////////
224 //////////////////////////////////////////////////////////////
225