/* PHYML : a program that computes maximum likelihood phylogenies from DNA or AA homologous sequences Copyright (C) Stephane Guindon. Oct 2003 onward All parts of the source except where indicated are distributed under the GNU public licence. See http://www.opensource.org for details. */ #include "utilities.h" #include "lk.h" #include "optimiz.h" #include "models.h" #include "free.h" /* int LIM_SCALE; */ /* double LIM_SCALE_VAL; */ /* double MDBL_MAX; */ /* double MDBL_MIN; */ /*********************************************************/ void Init_Tips_At_One_Site_Nucleotides(char state, double **p_lk) { switch(state){ case 'A' : (*p_lk)[0]=1.; (*p_lk)[1]=(*p_lk)[2]=(*p_lk)[3]=.0; break; case 'C' : (*p_lk)[1]=1.; (*p_lk)[0]=(*p_lk)[2]=(*p_lk)[3]=.0; break; case 'G' : (*p_lk)[2]=1.; (*p_lk)[1]=(*p_lk)[0]=(*p_lk)[3]=.0; break; case 'T' : (*p_lk)[3]=1.; (*p_lk)[1]=(*p_lk)[2]=(*p_lk)[0]=.0; break; case 'U' : (*p_lk)[3]=1.; (*p_lk)[1]=(*p_lk)[2]=(*p_lk)[0]=.0; break; case 'M' : (*p_lk)[0]=(*p_lk)[1]=1.; (*p_lk)[2]=(*p_lk)[3]=.0; break; case 'R' : (*p_lk)[0]=(*p_lk)[2]=1.; (*p_lk)[1]=(*p_lk)[3]=.0; break; case 'W' : (*p_lk)[0]=(*p_lk)[3]=1.; (*p_lk)[1]=(*p_lk)[2]=.0; break; case 'S' : (*p_lk)[1]=(*p_lk)[2]=1.; (*p_lk)[0]=(*p_lk)[3]=.0; break; case 'Y' : (*p_lk)[1]=(*p_lk)[3]=1.; (*p_lk)[0]=(*p_lk)[2]=.0; break; case 'K' : (*p_lk)[2]=(*p_lk)[3]=1.; (*p_lk)[0]=(*p_lk)[1]=.0; break; case 'B' : (*p_lk)[1]=(*p_lk)[2]=(*p_lk)[3]=1.; (*p_lk)[0]=.0; break; case 'D' : (*p_lk)[0]=(*p_lk)[2]=(*p_lk)[3]=1.; (*p_lk)[1]=.0; break; case 'H' : (*p_lk)[0]=(*p_lk)[1]=(*p_lk)[3]=1.; (*p_lk)[2]=.0; break; case 'V' : (*p_lk)[0]=(*p_lk)[1]=(*p_lk)[2]=1.; (*p_lk)[3]=.0; break; case 'N' : case 'X' : case '?' : case 'O' : case '-' : (*p_lk)[0]=(*p_lk)[1]=(*p_lk)[2]=(*p_lk)[3]=1.;break; default : { printf("\n. Unknown character state : %c\n",state); Exit("\n. Init failed (check the data type)\n"); break; } } } /*********************************************************/ void Init_Tips_At_One_Site_AA(char aa, double **p_lk) { int i; For(i,20) (*p_lk)[i] = .0; switch(aa){ case 'A' : (*p_lk)[0]= 1.; break;/* Alanine */ case 'R' : (*p_lk)[1]= 1.; break;/* Arginine */ case 'N' : (*p_lk)[2]= 1.; break;/* Asparagine */ case 'D' : (*p_lk)[3]= 1.; break;/* Aspartic acid */ case 'C' : (*p_lk)[4]= 1.; break;/* Cysteine */ case 'Q' : (*p_lk)[5]= 1.; break;/* Glutamine */ case 'E' : (*p_lk)[6]= 1.; break;/* Glutamic acid */ case 'G' : (*p_lk)[7]= 1.; break;/* Glycine */ case 'H' : (*p_lk)[8]= 1.; break;/* Histidine */ case 'I' : (*p_lk)[9]= 1.; break;/* Isoleucine */ case 'L' : (*p_lk)[10]=1.; break;/* Leucine */ case 'K' : (*p_lk)[11]=1.; break;/* Lysine */ case 'M' : (*p_lk)[12]=1.; break;/* Methionine */ case 'F' : (*p_lk)[13]=1.; break;/* Phenylalanin */ case 'P' : (*p_lk)[14]=1.; break;/* Proline */ case 'S' : (*p_lk)[15]=1.; break;/* Serine */ case 'T' : (*p_lk)[16]=1.; break;/* Threonine */ case 'W' : (*p_lk)[17]=1.; break;/* Tryptophan */ case 'Y' : (*p_lk)[18]=1.; break;/* Tyrosine */ case 'V' : (*p_lk)[19]=1.; break;/* Valine */ case 'B' : (*p_lk)[2]= 1.; break;/* Asparagine */ case 'Z' : (*p_lk)[5]= 1.; break;/* Glutamine */ case 'X' : case '?' : case '-' : For(i,20) (*p_lk)[i] = 1.; break; default : { printf("\n. Unknown character state : %c\n",aa); Exit("\n. Init failed (check the data type)\n"); break; } } } /*********************************************************/ void Get_All_Partial_Lk(arbre *tree, edge *b_fcus, node *a, node *d) { int i,j; double p1_lk1,p2_lk2; double ***p_lk,***p_lk_v1,***p_lk_v2; int catg,site; double ***Pij1,***Pij2; if(d->tax) return; else { int dir1,dir2; dir1=dir2=-1; For(i,3) if(d->v[i] != a) (dir1<0)?(dir1=i):(dir2=i); if(b_fcus->l < BL_MIN) b_fcus->l = BL_MIN; p_lk = (d == b_fcus->left)? (b_fcus->p_lk_left): (b_fcus->p_lk_rght); p_lk_v1 = (d == d->b[dir1]->left)? (d->b[dir1]->p_lk_rght): (d->b[dir1]->p_lk_left); p_lk_v2 = (d == d->b[dir2]->left)? (d->b[dir2]->p_lk_rght): (d->b[dir2]->p_lk_left); Pij1 = d->b[dir1]->Pij_rr; Pij2 = d->b[dir2]->Pij_rr; For(catg,tree->mod->n_catg) { For(site,tree->n_pattern) { For(i,tree->mod->ns) /*sort sum terms ? No global effect*/ { p1_lk1 = p2_lk2 = .0; For(j,tree->mod->ns) { p1_lk1 += (Pij1[catg][i][j] * p_lk_v1[site][catg][j]); p2_lk2 += (Pij2[catg][i][j] * p_lk_v2[site][catg][j]); } p_lk[site][catg][i] = p1_lk1*p2_lk2; if(p_lk[site][catg][i] < MDBL_MIN) { printf("\nWARNING : scaling is required at site %d\n",site); /* printf("Alpha = %f\n",tree->mod->alpha); */ /* Exit(""); */ } } } } } } /*********************************************************/ void Get_All_Partial_Lk_Scale(arbre *tree, edge *b_fcus, node *a, node *d) { int i,j; double p1_lk1,p2_lk2; double ***p_lk,***p_lk_v1,***p_lk_v2; int catg,site; double ***Pij1,***Pij2; double max_p_lk; /* double min_p_lk; */ double sum_scale_d1,sum_scale_d2; double try; p1_lk1 = p2_lk2 = .0; if(d->tax) return; else { int dir1,dir2; dir1=dir2=-1; For(i,3) if(d->v[i] != a) (dir1<0)?(dir1=i):(dir2=i); if(b_fcus->l < BL_MIN) b_fcus->l = BL_MIN; p_lk = (d == b_fcus->left)? (b_fcus->p_lk_left): (b_fcus->p_lk_rght); p_lk_v1 = (d == d->b[dir1]->left)? (d->b[dir1]->p_lk_rght): (d->b[dir1]->p_lk_left); p_lk_v2 = (d == d->b[dir2]->left)? (d->b[dir2]->p_lk_rght): (d->b[dir2]->p_lk_left); Pij1 = d->b[dir1]->Pij_rr; Pij2 = d->b[dir2]->Pij_rr; For(site,tree->n_pattern) { sum_scale_d1 = sum_scale_d2 = .0; (d == d->b[dir1]->left)? (sum_scale_d1 = d->b[dir1]->sum_scale_f_rght[site]): (sum_scale_d1 = d->b[dir1]->sum_scale_f_left[site]); (d == d->b[dir2]->left)? (sum_scale_d2 = d->b[dir2]->sum_scale_f_rght[site]): (sum_scale_d2 = d->b[dir2]->sum_scale_f_left[site]); (d == b_fcus->left)? (b_fcus->sum_scale_f_left[site] = sum_scale_d1 + sum_scale_d2): (b_fcus->sum_scale_f_rght[site] = sum_scale_d1 + sum_scale_d2); /* min_p_lk = MDBL_MAX; */ max_p_lk = MDBL_MIN; For(catg,tree->mod->n_catg) { For(i,tree->mod->ns) /*sort sum terms ? No global effect*/ { p1_lk1 = p2_lk2 = .0; For(j,tree->mod->ns) { p1_lk1 += (Pij1[catg][i][j] * p_lk_v1[site][catg][j]); p2_lk2 += (Pij2[catg][i][j] * p_lk_v2[site][catg][j]); } try = p1_lk1*p2_lk2; p_lk[site][catg][i]=try; if((p_lk[site][catg][i] > max_p_lk)) max_p_lk = p_lk[site][catg][i]; /* if((p_lk[site][catg][i] < min_p_lk)) min_p_lk = p_lk[site][catg][i]; */ } } if(max_p_lk < LIM_SCALE_VAL) { For(catg,tree->mod->n_catg) { For(i,tree->mod->ns) { p_lk[site][catg][i] /= max_p_lk; if(p_lk[site][catg][i] > MDBL_MAX) { Exit("\n. Numerical underflow ! (send me an e-mail : s.guindon@auckland.ac.nz)\n"); /* p_lk[site][catg][i] = p_lk[site][catg-1][i] ; */ } } } (d == b_fcus->left)? (b_fcus->sum_scale_f_left[site] += log(max_p_lk)): (b_fcus->sum_scale_f_rght[site] += log(max_p_lk)); } if(max_p_lk > (1./LIM_SCALE_VAL)) { For(catg,tree->mod->n_catg) { For(i,tree->mod->ns) { p_lk[site][catg][i] /= max_p_lk; if(p_lk[site][catg][i] < MDBL_MIN) { Exit("\n. Numerical overflow ! (send me an e-mail : s.guindon@auckland.ac.nz)\n"); /* p_lk[site][catg][i] = p_lk[site][catg-1][i] ; */ } } } (d == b_fcus->left)? (b_fcus->sum_scale_f_left[site] += log(max_p_lk)): (b_fcus->sum_scale_f_rght[site] += log(max_p_lk)); } } } } /*********************************************************/ void Post_Order_Lk(node *pere, node *fils, arbre *tree) { int i,dir1,dir2,dir3; dir1 = dir2 = dir3 = -1; if(fils->tax) return; else { For(i,3) { if(fils->v[i] != pere) { Post_Order_Lk(fils,fils->v[i],tree); if(dir1 < 0) dir1 = i; else dir2 = i; } else dir3 = i; } (tree->n_otu > LIM_SCALE)? (Get_All_Partial_Lk_Scale(tree,fils->b[dir3],pere,fils)): (Get_All_Partial_Lk(tree,fils->b[dir3],pere,fils)); } } /*********************************************************/ void Pre_Order_Lk(node *pere, node *fils, arbre *tree) { int i,j,dir1,dir2; dir1 = dir2 = -1; if(fils->tax) return; else { For(i,3) { if(fils->v[i] != pere) { For(j,3) { if(j != i) { if(dir1 < 0) dir1 = j; else dir2 = j; } } (tree->n_otu > LIM_SCALE)? (Get_All_Partial_Lk_Scale(tree,fils->b[i],fils->v[i],fils)): (Get_All_Partial_Lk(tree,fils->b[i],fils->v[i],fils)); dir1 = dir2 = -1; Pre_Order_Lk(fils,fils->v[i],tree); } } } } /*********************************************************/ void Lk(arbre *tree,allseq *alldata) { int br,site,catg; double len; Set_Model_Parameters(tree); For(br,2*tree->n_otu-3) { For(site,tree->n_pattern) { tree->t_edges[br]->sum_scale_f_rght[site] = .0; tree->t_edges[br]->sum_scale_f_left[site] = .0; } if(tree->t_edges[br]->l < BL_MIN) tree->t_edges[br]->l = BL_MIN; if(tree->t_edges[br]->l > BL_MAX) tree->t_edges[br]->l = BL_MAX; For(catg,tree->mod->n_catg) { len = tree->t_edges[br]->l*tree->mod->rr[catg]; if(len < BL_MIN) len = BL_MIN; PMat(len, tree->mod, &tree->t_edges[br]->Pij_rr[catg]); } } Post_Order_Lk(tree->noeud[0],tree->noeud[0]->v[0],tree); if(tree->both_sides) Pre_Order_Lk(tree->noeud[0], tree->noeud[0]->v[0], tree); tree->tot_loglk = .0; tree->curr_catg = 0; tree->curr_site = 0; For(site,tree->n_pattern) { tree->tot_loglk_sorted[site] = .0; tree->site_lk[site] = .0; tree->curr_site = site; Site_Lk(tree,alldata); } qksort(tree->tot_loglk_sorted, 0, tree->n_pattern-1); /* qsort(tree->tot_loglk_sorted,tree->n_pattern,sizeof(double),Sort_Double_Decrease); */ tree->tot_loglk = .0; For(site, tree->n_pattern) { if(tree->tot_loglk_sorted[site] < .0) /* WARNING : change cautiously */ tree->tot_loglk += tree->tot_loglk_sorted[site]; } For(br,2*tree->n_otu-3) { if(tree->t_edges[br]->get_p_lk_left) tree->t_edges[br]->ud_p_lk_left = 1; if(tree->t_edges[br]->get_p_lk_rght) tree->t_edges[br]->ud_p_lk_rght = 1; } } /*********************************************************/ void Site_Lk(arbre *tree, allseq *alldata) { int j,k,l,m; double *site_dlk=NULL, *site_d2lk=NULL; double log_site_lk,site_lk,site_lk_previous,aux; int left; edge *eroot; int is_ambigu; int state_root, state_elsewhere, site = tree->curr_site; if(alldata->wght[site] < MDBL_MIN) { tree->tot_loglk_sorted[tree->curr_site] = 1.; /* WARNING : change cautiously */ return; } if(tree->mod->s_opt->opt_bl) { site_dlk = (double *)mCalloc(2*tree->n_otu-3,sizeof(double)); site_d2lk = (double *)mCalloc(2*tree->n_otu-3,sizeof(double)); } eroot = tree->noeud[0]->b[0]; (eroot->rght->tax)?(left=1):(left=0); /* state = tree->data->invar[site]; */ is_ambigu = alldata->ambigu[site]; state_root = -1; state_root = Assign_State(alldata->c_seq[eroot->rght->num]->state + site, tree->mod->datatype, tree->mod->stepsize); state_elsewhere = -1; state_elsewhere = tree->data->invar[site]; /**/ /* is_ambigu = 1; */ /**/ if(tree->mod->s_opt->opt_bl) For(j,2*tree->n_otu-3) site_dlk[j] = site_d2lk[j] = 0.0; log_site_lk = site_lk = .0; site_lk_previous = .0; For(j,tree->mod->n_catg) { if(is_ambigu) { For(k,tree->mod->ns) /*sort sum terms ? No global effect*/ { For(l,tree->mod->ns) { site_lk += tree->mod->r_proba[j] * tree->mod->pi[k] * eroot->p_lk_left[site][j][k] * eroot->Pij_rr[j][k][l] * eroot->p_lk_rght[site][j][l]; } } } else { For(k,tree->mod->ns) /*sort sum terms ? No global effect*/ { site_lk += tree->mod->r_proba[j] * tree->mod->pi[k] * eroot->p_lk_left[site][j][k] * eroot->Pij_rr[j][k][state_root]; } } tree->log_site_lk_cat[j][tree->curr_site] = (site_lk/tree->mod->r_proba[j] - site_lk_previous); site_lk_previous = site_lk/tree->mod->r_proba[j]; if(tree->mod->s_opt->opt_bl) { For(k,2*tree->n_otu-3) { tree->t_edges[k]->site_dlk_rr[j] = .0; tree->t_edges[k]->site_d2lk_rr[j] = .0; For(l,tree->mod->ns) /*sort sum terms ? No global effect*/ { For(m,tree->mod->ns) { tree->t_edges[k]->site_dlk_rr[j] += tree->mod->pi[l] * tree->t_edges[k]->p_lk_left[site][j][l]* tree->t_edges[k]->p_lk_rght[site][j][m]* tree->t_edges[k]->dPij_rr[j][l][m]; tree->t_edges[k]->site_d2lk_rr[j] += tree->mod->pi[l] * tree->t_edges[k]->p_lk_left[site][j][l]* tree->t_edges[k]->p_lk_rght[site][j][m]* tree->t_edges[k]->d2Pij_rr[j][l][m]; } } site_dlk[k] += tree->t_edges[k]->site_dlk_rr[j]* tree->mod->r_proba[j]; site_d2lk[k] += tree->t_edges[k]->site_d2lk_rr[j]* tree->mod->r_proba[j]; } } } /* code 2.3 begin*/ if (!tree->mod->invar) { log_site_lk += log(site_lk) + eroot->sum_scale_f_left[site] + eroot->sum_scale_f_rght[site]; } else { if ((double)tree->data->invar[site] > -0.5) { if (!(eroot->sum_scale_f_left[site] + eroot->sum_scale_f_rght[site]==0.0)) site_lk *= exp(eroot->sum_scale_f_left[site] + eroot->sum_scale_f_rght[site]); log_site_lk = log(site_lk*(1.0-tree->mod->pinvar) + tree->mod->pinvar*tree->mod->pi[state_elsewhere]); } else { log_site_lk = log(site_lk*(1.0-tree->mod->pinvar)) + eroot->sum_scale_f_left[site] + eroot->sum_scale_f_rght[site]; } } For(j,tree->mod->n_catg) tree->log_site_lk_cat[j][tree->curr_site] = log(tree->log_site_lk_cat[j][tree->curr_site]) + + eroot->sum_scale_f_left[site] + eroot->sum_scale_f_rght[site]; /* code 2.3 end*/ if(log_site_lk < -MDBL_MAX) { printf("%d %E %f %f %f %f\n", site, log_site_lk, tree->mod->alpha, eroot->sum_scale_f_left[site], eroot->sum_scale_f_rght[site], tree->mod->pinvar); Exit("\nlog_site_lk < -MDBL_MAX\n"); } tree->site_lk[site] = log_site_lk; if(tree->mod->s_opt->opt_bl) { For(k,2*tree->n_otu-3) { aux = exp(tree->t_edges[k]->sum_scale_f_rght[site]+ tree->t_edges[k]->sum_scale_f_left[site]); site_dlk[k] *= aux; site_d2lk[k] *= aux; tree->t_edges[k]->site_dlk = site_dlk[k]; tree->t_edges[k]->site_d2lk = site_d2lk[k]; } } /* tree->tot_loglk += alldata->wght[site]*log_site_lk; */ tree->tot_loglk_sorted[site] = alldata->wght[site]*log_site_lk; if((tree->mod->s_opt->opt_bl) && (fabs(site_lk = exp(log_site_lk)) > sqrt(MDBL_MIN))) { For(k,2*tree->n_otu-3) { aux = site_dlk[k] / site_lk; tree->tot_dloglk[k] += alldata->wght[site] * aux; tree->tot_d2loglk[k] += alldata->wght[site] * (site_d2lk[k]/site_lk - aux*aux); } } if(tree->mod->s_opt->opt_bl) { Free(site_dlk); Free(site_d2lk); } } /*********************************************************/ double Lk_At_Given_Edge(arbre *tree, edge *b_fcus) { int site,catg,k,l,edge_num, ns = tree->mod->ns; double site_lk,log_site_lk,site_dlk,site_d2lk,aux; edge_num = b_fcus->num; tree->tot_loglk = log_site_lk = .0; tree->tot_dloglk[edge_num] = tree->tot_d2loglk[edge_num] = .0; tree->n_pattern = tree->data->crunch_len/tree->mod->stepsize; if(b_fcus->l < BL_MIN) b_fcus->l = BL_MIN; if(b_fcus->l > BL_MAX) b_fcus->l = BL_MAX; For(catg,tree->mod->n_catg) { aux = b_fcus->l*tree->mod->rr[catg]; if(aux < BL_MIN) aux = BL_MIN; PMat(aux, tree->mod,&b_fcus->Pij_rr[catg]); } if((tree->mod->s_opt->opt_bl) && (tree->mod->datatype == NT)) { For(catg,tree->mod->n_catg) { dPMat(tree->t_edges[edge_num]->l, tree->mod->rr[catg], tree->mod, &b_fcus->dPij_rr[catg]); d2PMat(tree->t_edges[edge_num]->l, tree->mod->rr[catg], tree->mod, &b_fcus->d2Pij_rr[catg]); } } For(site,tree->data->crunch_len) { if(tree->data->wght[site]) { log_site_lk = site_lk = .0; /* see equation (2) in phyml_tech_doc.pdf */ For(catg,tree->mod->n_catg) { For(k,ns) /*sort sum terms ? No global effect*/ { For(l,ns) { site_lk += tree->mod->r_proba[catg] * tree->mod->pi[k] * b_fcus->p_lk_left[site][catg][k] * b_fcus->p_lk_rght[site][catg][l] * b_fcus->Pij_rr[catg][k][l]; } } } /* code 2.3 begin*//* see equations in phyml_tech_doc.pdf */ /* compute log_site_lk */ if (!tree->mod->invar) { log_site_lk += log(site_lk) + b_fcus->sum_scale_f_left[site] + b_fcus->sum_scale_f_rght[site]; } else { if ((double)tree->data->invar[site] > -0.5) { if (!(b_fcus->sum_scale_f_left[site] + b_fcus->sum_scale_f_rght[site]==0.0)) site_lk *= exp(b_fcus->sum_scale_f_left[site] + b_fcus->sum_scale_f_rght[site]); log_site_lk = log(site_lk*(1.0-tree->mod->pinvar) + tree->mod->pinvar*tree->mod->pi[tree->data->invar[site]]); } else { log_site_lk = log(site_lk*(1.0-tree->mod->pinvar)) + b_fcus->sum_scale_f_left[site] + b_fcus->sum_scale_f_rght[site]; } } /* code 2.3 end*/ if(log_site_lk < -MDBL_MAX) Exit("\nlog_site_lk < -MDBL_MAX\n"); tree->site_lk[site] = log_site_lk; /*tree->tot_loglk += *//* old code */ tree->tot_loglk_sorted[site] = /* code 2.3 */ tree->data->wght[site]*log_site_lk; } else tree->tot_loglk_sorted[site] = 1.; /* WARNING : change cautiously */ } /* code 2.3 begin*/ /* sort and add numbers from smallest to biggest */ qksort(tree->tot_loglk_sorted, 0, tree->n_pattern-1); /* qsort(tree->tot_loglk_sorted,tree->n_pattern,sizeof(double),Sort_Double_Decrease); */ tree->tot_loglk = .0; For(k, tree->data->crunch_len) if(tree->tot_loglk_sorted[k] < .0) /* WARNING : change cautiously */ tree->tot_loglk += tree->tot_loglk_sorted[k]; /* code 2.3 end*/ if((tree->mod->s_opt->opt_bl) && (tree->mod->datatype == NT)) { For(site,tree->n_pattern) { if(tree->data->wght[site]) { site_dlk = site_d2lk = .0; b_fcus->site_dlk = b_fcus->site_d2lk = .0; For(catg,tree->mod->n_catg) { For(k,ns) /*sort sum terms ? No global effect*/ { For(l,ns) { aux = tree->mod->r_proba[catg] * tree->mod->pi[k] * b_fcus->p_lk_left[site][catg][k]; site_dlk += aux * b_fcus->p_lk_rght[site][catg][l] * b_fcus->dPij_rr[catg][k][l]; site_d2lk += aux * b_fcus->p_lk_rght[site][catg][l] * b_fcus->d2Pij_rr[catg][k][l]; } } } if(tree->n_otu > LIM_SCALE) { aux = exp(b_fcus->sum_scale_f_rght[site]+ b_fcus->sum_scale_f_left[site]); site_dlk *= aux; site_d2lk *= aux; } b_fcus->site_dlk = site_dlk; b_fcus->site_d2lk = site_d2lk; if(fabs(exp(log_site_lk)) > sqrt(MDBL_MIN)) { aux = site_dlk / tree->site_lk[site]; tree->tot_dloglk[edge_num] += tree->data->wght[site] * aux; tree->tot_d2loglk[edge_num] += tree->data->wght[site] * (site_d2lk/tree->site_lk[site] - aux*aux); } } } } return tree->tot_loglk; } /*********************************************************/ void Update_P(arbre *tree, int t_edge_num) { int i; double len; len = -1.0; For(i,tree->mod->n_catg) { tree->curr_catg = i; len = tree->t_edges[t_edge_num]->l*tree->mod->rr[i]; if(len < BL_MIN) len = BL_MIN; tree->mod->update_eigen = 0; PMat(len,tree->mod,&tree->t_edges[t_edge_num]->Pij_rr[i]); /* Derivatives(tree->t_edges[t_edge_num],tree); */ } } /*********************************************************/ double Return_Lk(arbre *tree) { Lk(tree,tree->data); return tree->tot_loglk; } /*********************************************************/ double Return_Abs_Lk(arbre *tree) { Lk(tree,tree->data); return fabs(tree->tot_loglk); } /*********************************************************/ matrix *ML_Dist(allseq *data, model *mod) { int j,k,l; double init; int n_catg; double d_max; matrix *mat; allseq *twodata,*tmpdata; tmpdata = (allseq *)mCalloc(1,sizeof(allseq)); tmpdata->n_otu=2; tmpdata->c_seq = (seq **)mCalloc(2,sizeof(seq *)); tmpdata->b_frq = (double *)mCalloc(mod->ns,sizeof(double)); tmpdata->ambigu = (int *)mCalloc(data->crunch_len,sizeof(int)); tmpdata->crunch_len = tmpdata->init_len = data->crunch_len; mat = (mod->datatype == NT) ? ((mod->whichmodel < 10)?(K2P_dist(data,2000)):(JC69_Dist(data,mod))): (JC69_Dist(data,mod)); n_catg = -1; For(j,data->n_otu-1) { tmpdata->c_seq[0]=data->c_seq[j]; tmpdata->c_seq[0]->name=data->c_seq[j]->name; tmpdata->wght = data->wght; for(k=j+1;kn_otu;k++) { tmpdata->c_seq[1]=data->c_seq[k]; tmpdata->c_seq[1]->name=data->c_seq[k]->name; twodata = Compact_CSeq(tmpdata,mod); For(l,mod->ns) twodata->b_frq[l]=data->b_frq[l]; Check_Ambiguities(twodata,mod->datatype,1); Hide_Ambiguities(twodata); init = mat->dist[j][k]; if((init == DIST_MAX) || (init < .0)) init = 0.1; n_catg = mod->n_catg; mod->n_catg = 1; /* BRENT */ d_max = Optimize_Dist(mod,init,twodata); /* d_max = init; */ /* NEWTON-RAPHSON */ /* if(d_max < .0) */ /* { */ /* d_max = Optimize_One_Dist(twodata,0,1,init,mod); */ /* d_max = init; */ /* } */ mod->n_catg = n_catg; if(d_max >= DIST_MAX) { printf("\n. Large distance encountered between %s and %s sequences\n", tmpdata->c_seq[1]->name, tmpdata->c_seq[0]->name); d_max = DIST_MAX; } mat->dist[j][k] = d_max; mat->dist[k][j] = mat->dist[j][k]; Free_Cseq(twodata); } } Free(tmpdata->ambigu); Free(tmpdata->b_frq); Free(tmpdata->c_seq); free(tmpdata); return mat; } /*********************************************************/ double Lk_Given_Two_Seq(allseq *data, int numseq1, int numseq2, double dist, model *mod, double *loglk, double *dloglk, double *d2loglk) { seq *seq1,*seq2; double site_lk,site_dlk,site_d2lk,log_site_lk; int i,j,k,l; double **p_lk_l,**p_lk_r; double len; DiscreteGamma (mod->r_proba, mod->rr, mod->alpha, mod->alpha,mod->n_catg,0); seq1 = data->c_seq[numseq1]; seq2 = data->c_seq[numseq2]; p_lk_l = (double **)mCalloc(data->c_seq[0]->len,sizeof(double *)); p_lk_r = (double **)mCalloc(data->c_seq[0]->len,sizeof(double *)); For(i,data->c_seq[0]->len) { p_lk_l[i] = (double *)mCalloc(mod->ns,sizeof(double)); p_lk_r[i] = (double *)mCalloc(mod->ns,sizeof(double)); } if(dist < BL_MIN) dist = BL_START; For(i,mod->n_catg) { len = dist*mod->rr[i]; if(len < BL_MIN) len = BL_MIN; PMat(len,mod,&(mod->Pij_rr[i])); } if(mod->datatype == NT) { For(i,mod->n_catg) { dPMat(dist,mod->rr[i],mod,&(mod->dPij_rr[i])); d2PMat(dist,mod->rr[i],mod,&(mod->d2Pij_rr[i])); } For(i,data->c_seq[0]->len) { Init_Tips_At_One_Site_Nucleotides(seq1->state[i], &p_lk_l[i]); Init_Tips_At_One_Site_Nucleotides(seq2->state[i], &p_lk_r[i]); } } else { For(i,data->c_seq[0]->len) { Init_Tips_At_One_Site_AA(seq1->state[i], &p_lk_l[i]); Init_Tips_At_One_Site_AA(seq2->state[i], &p_lk_r[i]); } } site_lk = site_dlk = site_d2lk = .0; *loglk = *dloglk = *d2loglk = 0; For(i,data->c_seq[0]->len) { if(data->wght[i]) { site_lk = log_site_lk = .0; if(!data->ambigu[i]) { For(k,mod->ns) {if(p_lk_l[i][k] > .0001) break;} For(l,mod->ns) {if(p_lk_r[i][l] > .0001) break;} For(j,mod->n_catg) { site_lk += mod->r_proba[j] * mod->pi[k] * p_lk_l[i][k] * mod->Pij_rr[j][k][l] * p_lk_r[i][l]; } } else { For(j,mod->n_catg) { For(k,mod->ns) /*sort sum terms ? No global effect*/ { For(l,mod->ns) { site_lk += mod->r_proba[j] * mod->pi[k] * p_lk_l[i][k] * mod->Pij_rr[j][k][l] * p_lk_r[i][l]; } } } } /* printf("'%c' '%c' -> %f\n",seq1->state[i],seq2->state[i],site_lk); */ if(site_lk <= .0) { printf("'%c' '%c'\n",seq1->state[i],seq2->state[i]); Exit("\n. Err: site lk <= 0\n"); } log_site_lk += log(site_lk); *loglk += data->wght[i] * log_site_lk;/* sort sum terms ? No global effect*/ } } For(i,data->c_seq[0]->len) { Free(p_lk_l[i]); Free(p_lk_r[i]); } Free(p_lk_l); Free(p_lk_r); return *loglk; } /*********************************************************/ double ***Get_Partial_Lk_Struct(arbre *tree, int len, int n_catg) { double ***p_lk; int j,k; p_lk = (double ***)mCalloc(len,sizeof(double **)); For(j,len) { p_lk[j] = (double **)mCalloc((int)n_catg,sizeof(double *)); For(k,n_catg) p_lk[j][k] = (double *)mCalloc(tree->mod->ns,sizeof(double )); } return p_lk; } /*********************************************************/ void Unconstraint_Lk(arbre *tree) { int i; tree->unconstraint_lk = .0; For(i,tree->data->crunch_len) { tree->unconstraint_lk += tree->data->wght[i]*log(tree->data->wght[i]); } tree->unconstraint_lk -= tree->data->init_len*log(tree->data->init_len); } /*********************************************************/ void Update_P_Lk(arbre *tree, edge *b_fcus, node *n) { /* | |<- b_cus | n / \ / \ / \ */ int k,l; int site, catg; double ***p_lk, ***p_lk_v1, ***p_lk_v2; double **Pij1, **Pij2; double *n_scale_f, *d1_scale_f, *d2_scale_f; double p1_lk1,p2_lk2; double max_p_lk; edge *b1, *b2; b1 = b2 = NULL; p_lk = p_lk_v1 = p_lk_v2 = NULL; max_p_lk = MDBL_MIN; if(n == b_fcus->left) { /* if(b_fcus->ud_p_lk_left) { printf("This p_lk is up to date\n"); return;} */ p_lk = b_fcus->p_lk_left; p_lk_v1 = (n == n->b[b_fcus->l_v1]->left)? (n->b[b_fcus->l_v1]->p_lk_rght): (n->b[b_fcus->l_v1]->p_lk_left); p_lk_v2 = (n == n->b[b_fcus->l_v2]->left)? (n->b[b_fcus->l_v2]->p_lk_rght): (n->b[b_fcus->l_v2]->p_lk_left); n_scale_f = b_fcus->sum_scale_f_left; d1_scale_f = (n == n->b[b_fcus->l_v1]->left)? (n->b[b_fcus->l_v1]->sum_scale_f_rght): (n->b[b_fcus->l_v1]->sum_scale_f_left); d2_scale_f = (n == n->b[b_fcus->l_v2]->left)? (n->b[b_fcus->l_v2]->sum_scale_f_rght): (n->b[b_fcus->l_v2]->sum_scale_f_left); b_fcus->get_p_lk_left = 1; b_fcus->ud_p_lk_left = 1; } else { p_lk = b_fcus->p_lk_rght; p_lk_v1 = (n == n->b[b_fcus->r_v1]->left)? (n->b[b_fcus->r_v1]->p_lk_rght): (n->b[b_fcus->r_v1]->p_lk_left); p_lk_v2 = (n == n->b[b_fcus->r_v2]->left)? (n->b[b_fcus->r_v2]->p_lk_rght): (n->b[b_fcus->r_v2]->p_lk_left); n_scale_f = b_fcus->sum_scale_f_rght; d1_scale_f = (n == n->b[b_fcus->r_v1]->left)? (n->b[b_fcus->r_v1]->sum_scale_f_rght): (n->b[b_fcus->r_v1]->sum_scale_f_left); d2_scale_f = (n == n->b[b_fcus->r_v2]->left)? (n->b[b_fcus->r_v2]->sum_scale_f_rght): (n->b[b_fcus->r_v2]->sum_scale_f_left); b_fcus->get_p_lk_rght = 1; b_fcus->ud_p_lk_rght = 1; } if(b_fcus->l < BL_MIN) b_fcus->l = BL_MIN; if(n == b_fcus->left) { b1 = n->b[b_fcus->l_v1]; b2 = n->b[b_fcus->l_v2]; } else { b1 = n->b[b_fcus->r_v1]; b2 = n->b[b_fcus->r_v2]; } if(tree->n_otu <= LIM_SCALE) { /* NO SCALE */ For(site,tree->n_pattern) { if(tree->data->wght[site]) { For(catg,tree->mod->n_catg) { Pij1 = b1->Pij_rr[catg]; Pij2 = b2->Pij_rr[catg]; For(k,tree->mod->ns) /*sort sum terms ? No global effect*/ { p1_lk1 = p2_lk2 = .0; For(l,tree->mod->ns) { p1_lk1 += Pij1[k][l] * p_lk_v1[site][catg][l]; p2_lk2 += Pij2[k][l] * p_lk_v2[site][catg][l]; } p_lk[site][catg][k] = p1_lk1 * p2_lk2; } } } } } else { /* SCALE */ For(site,tree->n_pattern) { if(tree->data->wght[site]) { For(catg,tree->mod->n_catg) { Pij1 = b1->Pij_rr[catg]; Pij2 = b2->Pij_rr[catg]; if(!catg) { n_scale_f[site] = d1_scale_f[site] + d2_scale_f[site]; max_p_lk = -MDBL_MAX; } For(k,tree->mod->ns) /*sort sum terms ? No global effect*/ { p_lk[site][catg][k] = .0; p1_lk1 = p2_lk2 = .0; For(l,tree->mod->ns) { p1_lk1 += Pij1[k][l] * p_lk_v1[site][catg][l]; p2_lk2 += Pij2[k][l] * p_lk_v2[site][catg][l]; } p_lk[site][catg][k] = p1_lk1 * p2_lk2; if((p_lk[site][catg][k] > max_p_lk)) max_p_lk = p_lk[site][catg][k]; } } if(max_p_lk < LIM_SCALE_VAL) { For(catg,tree->mod->n_catg) { For(k,tree->mod->ns) { p_lk[site][catg][k] /= max_p_lk; if(p_lk[site][catg][k] > MDBL_MAX) { Exit("\n. Numerical overflow ! (send me an e-mail : s.guindon@auckland.ac.nz)\n"); } } } n_scale_f[site] += log(max_p_lk); } if(max_p_lk > (1./LIM_SCALE_VAL)) { For(catg,tree->mod->n_catg) { For(k,tree->mod->ns) { p_lk[site][catg][k] /= max_p_lk; if(p_lk[site][catg][k] < MDBL_MIN) { printf("\n. Numerical underflow ! (send me an e-mail : s.guindon@auckland.ac.nz)\n"); } } } n_scale_f[site] += log(max_p_lk); } } } } } /*********************************************************/ void Make_Tree_4_Lk(arbre *tree, allseq *alldata) { int i; if(!tree->tot_loglk_sorted) tree->tot_loglk_sorted = (double *)mCalloc(tree->n_pattern, sizeof(double)); if(!tree->tot_dloglk) tree->tot_dloglk = (double *)mCalloc(2*tree->n_otu-3,sizeof(double)); if(!tree->tot_d2loglk) tree->tot_d2loglk = (double *)mCalloc(2*tree->n_otu-3,sizeof(double)); if(!tree->site_lk) tree->site_lk = (double *)mCalloc(alldata->crunch_len,sizeof(double)); if(!tree->log_site_lk_cat) { tree->log_site_lk_cat = (double **)mCalloc(tree->mod->n_catg,sizeof(double *)); For(i,tree->mod->n_catg) tree->log_site_lk_cat[i] = (double *)mCalloc(alldata->crunch_len,sizeof(double)); } tree->root = tree->noeud[0]; For(i,2*tree->n_otu-3) { Make_Edge_Lk(tree->t_edges[i]->left, tree->t_edges[i]->rght, tree); } Alloc_All_P_Lk(tree); /* Make_P_Lk_Struct(tree); */ Init_P_Lk_Tips(tree); } /*********************************************************/ void Init_P_Lk_Tips(arbre *tree) { int curr_site,i,j,k; Fors(curr_site,tree->data->crunch_len,tree->mod->stepsize) { For(i,tree->n_otu) { if (tree->mod->datatype == NT) { Init_Tips_At_One_Site_Nucleotides(tree->data->c_seq[i]->state[curr_site], &tree->noeud[i]->b[0]->p_lk_rght[curr_site][0]); } else Init_Tips_At_One_Site_AA(tree->data->c_seq[i]->state[curr_site], &tree->noeud[i]->b[0]->p_lk_rght[curr_site][0]); if((tree->noeud[i]->b[0]->p_lk_rght) && (tree->noeud[i]->b[0]->get_p_lk_rght)) { for(j=1;jmod->n_catg;j++) { For(k,tree->mod->ns) { tree->noeud[i]->b[0]->p_lk_rght[curr_site][j][k]= tree->noeud[i]->b[0]->p_lk_rght[curr_site][0][k]; } } } else { for(j=1;jmod->n_catg;j++) For(k,tree->mod->ns) tree->noeud[i]->b[0]->p_lk_rght[curr_site][j][k]= tree->noeud[i]->b[0]->p_lk_rght[curr_site][0][k]; } } } } /*********************************************************/