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