/* A general model combining Tuffley-Steel, Huelsenbeck and Galtier models */ /* This code was based on comte_ml.c in NHML (Galtier & Gouy, 1998; Galtier, 2001) */ #include "mynhmlg.h" #include #include #define PointGamma(prob,alpha,beta) PointChi2(prob,2.0*(alpha))/(2.0*(beta)) time_t tdeb, tfin, ttotdeb, ttotfin; double probatime=0., dliketime=0., d2liketime=0., tottime; FILE* probaoutcov; extern tree* s_tree; extern compute_option compute_o; extern noeud root; extern char infileName[100]; extern char intreeName[100]; extern double freqaa[20]; extern double jttQmatrix[20][20]; extern double eigmat[]; extern double probmat[20][20]; extern double JTTeigenvect[20][20]; #define pos(i,j,n) ((i)*(n)+(j)) upgmanoeud upgmaroot; double min_value[4]={frmin, lmin, covmin, pimin}; double max_value[4]={frmax, lmax, covmax, pimax}; double gamma_range[GAMMA_RANGE_SIZE]={0.1, 0.2, 0.4, 0.6, 0.8, 1., 2., 5., 10.}; int paramcase[33][17]; double *gamma_clmean; int gamma_nbcl; int gamma_nbcl1; // gamma_nbcl1 = gamma_nbcl*2 double covar; // not used in current version double covar1, covar2, covar3; // corresponding to s01, s10, s11 double pi; // Galtier's pi for proportion of covarions sites double **P; int noprint; proba **p1prov, **p2prov; dproba **dp1prov, **dp2prov; /* utilities */ dproba **d2p1prov, **d2p2prov; FILE* dbg, *dbgprov; char aachar[26]="ARNDCQEGHILKMFPSTWYVBZX?*-"; void fprinttime(double rtime, FILE *out){ int d, h, m, s; d=(int)rtime/86400; rtime-=d*86400; if(d) fprintf(out,"%d day", d); if(d>1) fprintf(out,"s"); if(d) fprintf(out," "); h=(int)rtime/3600; rtime-=h*3600; if(h) fprintf(out,"%d hour", h); if(h>1) fprintf(out,"s"); if(h) fprintf(out," "); m=(int)rtime/60; rtime-=m*60; if(m) fprintf(out,"%d minute", m); if(m>1) fprintf(out,"s"); if(m) fprintf(out," "); s=(int)rtime; if(s) fprintf(out,"%d second", s); if(s>1) fprintf(out,"s"); if(d==0 && h==0 && m==0 && s==0) fprintf(out,"less than 1 second"); fprintf(out,"\n"); } int eigen(int job, double A[], int n, double rr[], double ri[], double vr[], double vi[], double w[]); /* init_cas */ /* Initialise le tableau paramcase qui sert au calcul des derivees secondes */ /* et croisees de la vraisemblance. */ /* premier indice (1->32): type de cas */ /* Un cas est un type de "relations" entre deux parametres vis-a-vis */ /* du processus de derivation. Il en existe 32, selon la nature des 2 */ /* parametres (branch length, Ts/Tv ration, GC equilibre, ...) et leur */ /* position dans l'arbre par rapport au noeud ou l'on calcule la derivee. */ /* deuxieme indice (1->16): rang des 16 termes de l'expression de la derivee */ /* seconde d'un produit de 4 facteurs: */ /* */ /* d2(p1.L1.p2.L2)/dxdy = d2p1/dxdy . L1 . p2 . L2 (terme 1) */ /* + dp1/dx . dL1/dy . p2 . L2 (terme 2) */ /* + dp1/dx . L1 . dp2/dy . L2 (terme 3) */ /* + ... */ /* ... */ /* + p1 . L1 . p2 . d2L2/dxdy (terme 16) */ /* */ /* param_case explicite les termes de la somme qui valent 0 (valeur a 0) et ceux */ /* qu'il faut calculer (valeur a 1), pour chaque cas. */ void init_cas(){ int i, j; for(i=0;i<33;i++) for(j=0;j<17;j++) paramcase[i][j]=0; for(i=1;i<=16;i++) paramcase[i][i]=1; for(i=1;i<=16;i++) paramcase[17][i]=1; for(i=1;i<=8;i++) paramcase[18][2*i-1]=1; for(i=1;i<=4;i++) paramcase[19][4*i-3]=1; for(i=1;i<=4;i++) paramcase[20][4*i-1]=1; for(i=1;i<=4;i++) paramcase[21][4*i-2]=1; for(i=1;i<=4;i++) paramcase[22][4*i]=1; paramcase[23][1]=1; paramcase[23][3]=1; paramcase[23][9]=1; paramcase[23][11]=1; paramcase[24][2]=1; paramcase[24][10]=1; paramcase[25][4]=1; paramcase[25][12]=1; paramcase[26][1]=1; paramcase[26][9]=1; paramcase[27][3]=1; paramcase[27][11]=1; paramcase[28][5]=1; paramcase[28][7]=1; paramcase[29][13]=1; paramcase[29][15]=1; paramcase[30][1]=1; paramcase[30][3]=1; paramcase[31][9]=1; paramcase[31][11]=1; } /* which_cas */ /* Determine dans quel cas (vis-a-vis du processus de derivation) */ /* se trouve une paire de parametres, decrits par leur nature */ /* np1 et np2 (0=Ts/Tv ratio, 1=root fraction, 2=GCanc, 3=branch length, */ /* 4=GC eq) et par leur position dans l'arbre par rapport au noeud */ /* dont on calcule la derivee, und1 et und2 (0 = au dessus, 1=en dessous */ /* a gauche mais pas branche fille, 2=en dessous a droite mais pas branche */ /* fille, 10=branche fille gauche, 20=branche fille droite, 12=branche */ /* racine). Si le noeud ou l'on calcule la derivee est la racine, rootnode */ /* est passe a 1 (cas particulier). */ int which_cas(int np1, int np2, int und1, int und2, int rootnode){ // if(np1==2 || np2==2) { printf("GCanc exception\n"); return -1; } if(np2b) return(a); else return(b); } */ void printmemory_init(){ char *line, *ret; int maxlline=1000; dbg=fopen("debugNHML", "w"); system("ps -aux > debugNHML.prov"); dbgprov=fopen("debugNHML.prov", "r"); if(dbgprov==NULL){printf("Warning: problem when debugging\n"); return;} line=check_alloc(maxlline+1, sizeof(char)); ret=fgets(line, maxlline, dbgprov); if(ret!=NULL) fprintf(dbg, " %s", line); fclose(dbgprov); fclose(dbg); free(line); } void printmemory(char title_arg[]){ char* line, *ret, *syscall, *title, *prov; int maxlline=1000; int thispid, lgtit=10, lgarg, i; thispid=(int)getpid(); lgarg=(int)strlen(title_arg); title=check_alloc(lgtit+3, sizeof(char)); sprintf(title, "%.10s", title_arg); prov=title; while(*prov) prov++; for(i=0;i debugNHML.prov", thispid); system(syscall); dbgprov=fopen("debugNHML.prov", "r"); dbg=fopen("debugNHML", "a"); line=check_alloc(maxlline+1, sizeof(char)); while(1){ ret=fgets(line, maxlline, dbgprov); if(ret==NULL) break; fprintf(dbg, "%s: %s", title, line); } fclose(dbgprov); fclose(dbg); system("rm -f debugNHML.prov"); free(line); free(syscall); free(title); } /* ludcmp */ /* from Numerical Recipes in C */ /* Replace matrix a by a rowwise permutation of its LU decomposition. */ int ludcmp(double **a, int n, int *indx, double *d){ int i,imax,j,k; double big,dum,sum,temp; double *vv; vv=vector(1,n); *d=1.0; for (i=1;i<=n;i++) { big=0.0; for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; if (big == 0.0) return 0; vv[i]=1.0/big; } for (j=1;j<=n;j++) { for (i=1;i= big) { big=dum; imax=i; } } if (j != imax) { for (k=1;k<=n;k++) { dum=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=dum; } *d = -(*d); vv[imax]=vv[j]; } indx[j]=imax; if (a[j][j] == 0.0) a[j][j]=1.0e-20; if (j != n) { dum=1.0/(a[j][j]); for (i=j+1;i<=n;i++) a[i][j] *= dum; } } free_vector(vv,1,n); return 1; } /* lubksb */ /* from Numerical Recipes in C */ void lubksb(double **a, int n, int *indx, double b[]){ int i,ii=0,ip,j; double sum; for (i=1;i<=n;i++) { ip=indx[i]; sum=b[ip]; b[ip]=b[i]; if (ii) for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; else if (sum) ii=i; b[i]=sum; } for (i=n;i>=1;i--) { sum=b[i]; for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; b[i]=sum/a[i][i]; } } /* matvect */ /* Return the product vector of matrix m1(nbl1, nbc1) by column vector v2. */ int matvect(double** m1, int nbl1, int nbc1, double* v2, int nbl2, double* vprod){ int i, k; if(nbc1!=nbl2) return 0; for(i=0;i0, accurate to 10 decimal places. Stirling's formula is used for the central polynomial part of the procedure. Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function. Communications of the Association for Computing Machinery, 9:684 */ double x=alpha, f=0, z; if (x<7) { f=1; z=x-1; while (++z<7) f*=z; x=z; f=-log(f); } z = 1/(x*x); return f + (x-0.5)*log(x) - x + .918938533204673 + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z +.083333333333333)/x; } double IncompleteGamma (double x, double alpha, double ln_gamma_alpha) { /* returns the incomplete gamma ratio I(x,alpha) where x is the upper limit of the integration and alpha is the shape parameter. returns (-1) if in error ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant. (1) series expansion if (alpha>x || x<=1) (2) continued fraction otherwise RATNEST FORTRAN by Bhattacharjee GP (1970) The incomplete gamma integral. Applied Statistics, 19: 285-287 (AS32) */ int i; double p=alpha, g=ln_gamma_alpha; double accurate=1e-8, overflow=1e30; double factor, gin=0, rn=0, a=0,b=0,an=0,dif=0, term=0, *pn; if (x==0) return (0); if (x<0 || p<=0) return (-1); pn=check_alloc(6, sizeof(double)); factor=exp(p*log(x)-x-g); if (x>1 && x>=p) goto l30; /* (1) series expansion */ gin=1; term=1; rn=p; l20: rn++; term*=x/rn; gin+=term; if (term > accurate) goto l20; gin*=factor/p; goto l50; l30: /* (2) continued fraction */ a=1-p; b=a+x+1; term=0; pn[0]=1; pn[1]=x; pn[2]=x+1; pn[3]=x*b; gin=pn[2]/pn[3]; l32: a++; b+=2; term++; an=a*term; for (i=0; i<2; i++) pn[i+4]=b*pn[i+2]-an*pn[i]; if (pn[5] == 0) goto l35; rn=pn[4]/pn[5]; dif=fabs(gin-rn); if (dif>accurate) goto l34; if (dif<=accurate*rn) goto l42; l34: gin=rn; l35: for (i=0; i<4; i++) pn[i]=pn[i+2]; if (fabs(pn[4]) < overflow) goto l32; for (i=0; i<4; i++) pn[i]/=overflow; goto l32; l42: gin=1-factor*gin; l50: free(pn); return (gin); } double PointChi2 (double prob, double v) { /* returns z so that Prob{x.999998 || v<=0) return (-1); g = LnGamma (v/2); xx=v/2; c=xx-1; if (v >= -1.24*log(p)) goto l1; ch=pow((p*xx*exp(g+xx*aa)), 1/xx); if (ch-e<0) return (ch); goto l4; l1: if (v>.32) goto l3; ch=0.4; a=log(1-p); l2: q=ch; p1=1+ch*(4.67+ch); p2=ch*(6.73+ch*(6.66+ch)); t=-0.5+(4.67+2*ch)/p1 - (6.73+ch*(13.32+3*ch))/p2; ch-=(1-exp(a+g+.5*ch+c*aa)*p2/p1)/t; if (fabs(q/ch-1)-.01 <= 0) goto l4; else goto l2; l3: x=PointNormal (p); p1=0.222222/v; ch=v*pow((x*sqrt(p1)+1-p1), 3.0); if (ch>2.2*v+6) ch=-2*(log(1-p)-c*log(.5*ch)+g); l4: q=ch; p1=.5*ch; if ((t=IncompleteGamma (p1, xx, g))<0) { printf ("\nerr IncompleteGamma"); return (-1); } p2=p-t; t=p2*exp(xx*aa+g+p1-c*log(ch)); b=t/ch; a=0.5*t-b*c; s1=(210+a*(140+a*(105+a*(84+a*(70+60*a))))) / 420; s2=(420+a*(735+a*(966+a*(1141+1278*a))))/2520; s3=(210+a*(462+a*(707+932*a)))/2520; s4=(252+a*(672+1182*a)+c*(294+a*(889+1740*a)))/5040; s5=(84+264*a+c*(175+606*a))/2520; s6=(120+c*(346+127*c))/5040; ch+=t*(1+0.5*t*s1-b*c*(s1-b*(s2-b*(s3-b*(s4-b*(s5-b*s6)))))); if (fabs(q/ch-1) > e) goto l4; return (ch); } int DiscreteGamma (double freqK[], double rK[], double alfa, double beta, int K, int median) { /* discretization of gamma distribution with equal proportions in each category */ int i; double gap05=1.0/(2.0*K), t, factor=alfa/beta*K, lnga1; if (median) { for (i=0; inbseq;j++) if(s_tree->node[j]->seq[i]!=site[j]) return 0; return 1; } char* readtree(FILE* treef, int nbsp){ char *line; char *ctree, *ctreefin, *prov; int i; line=(char*)check_alloc(MAXLLINE+1, sizeof(char)); ctree=check_alloc(3*MAXLNAME*nbsp, sizeof(char)); ctreefin=ctree; line[0]='\n'; while(line[0]=='\n') fgets(line, MAXLLINE, treef); i=0; while(line[i]==' ' || line[i]=='\t') i++; if(line[i]!='(' && line[i]!='['){ printf("Tree first character should be '[' or '('\n"); return NULL; } if(line[i]=='['){ while((prov=strchr(line, ']'))==NULL){ if(!fgets(line, MAXLLINE, treef)){ printf("Unterminated comment in tree ('[' missing)\n"); return NULL; } } prov++; while(*prov==' ' || *prov=='\t') prov++; if(*prov=='\n'){ line[0]='\n'; while(line[0]=='\n' && fgets(line, MAXLLINE, treef)); prov=line; while(*prov==' ' || *prov=='\t') prov++; } if(*prov!='('){ printf("Comments should be followed by '('\n"); return NULL; } } else prov=line+i; /* prov = premiere parenthese */ sprintf(ctree, "%s", prov); while(*ctreefin!='\n' && *ctreefin!=';' && *ctreefin) ctreefin++; if(*ctreefin==';') return ctree; while(fgets(line, MAXLLINE, treef)){ sprintf(ctreefin, "%s", line); if( (prov=strchr(line, ';')) ){ *(prov+1)=0; return ctree; } while(*ctreefin!='\n' && *ctreefin) ctreefin++; } *ctreefin=0; free(line); return ctree; } /* upgma_ctreewrite */ char* upgma_ctreewrite(upgmanoeud comesfrom, upgmanoeud node, char* ctree, int nodegc){ upgmanoeud vd, vg; double l, b; if(comesfrom==node->v1) { vg=node->v2; vd=node->v3; l=node->l1; b=node->b1;} else if(comesfrom==node->v2) { vg=node->v1; vd=node->v3; l=node->l2; b=node->b2;} else { vg=node->v1; vd=node->v2; l=node->l3; b=node->b3;} if (vg==NULL){ sprintf(ctree,"%s", node->nom); while(*ctree) ctree++; if (l>-0.5) sprintf(ctree,":%.4f",l); while(*ctree) ctree++; return ctree; } else{ *(ctree++)='('; ctree=upgma_ctreewrite(node, vg, ctree, nodegc); } sprintf(ctree,", "); ctree+=2; ctree=upgma_ctreewrite(node, vd, ctree, nodegc); *(ctree++)=')'; if (b>-0.1) sprintf(ctree,"%.2f",b); while(*ctree) ctree++; if (l>-0.9) sprintf(ctree,":%.6f", l); while(*ctree) ctree++; return ctree; } /* upgma_stoc */ void upgma_stoc(upgmanoeud* arbre_s, int racine, int notu, char* ctree, int nodegc){ char* tring, *prov; tring=(char*)check_alloc((50+MAXLNAME)*notu,sizeof(char)); if (racine) { tring=upgma_ctreewrite(NULL, upgmaroot, ctree, nodegc); prov=ctree+strlen(ctree)-1; while(*prov!=')') *(prov--)='\0'; strcat(ctree,";"); } else{ tring=upgma_ctreewrite(arbre_s[2*notu-3]->v3, arbre_s[2*notu-3], ctree, nodegc); prov=ctree+strlen(ctree)-1; while(*prov!=')') *(prov--)='\0'; *(prov++)=','; *(prov++)=' '; tring=upgma_ctreewrite(arbre_s[2*notu-3], arbre_s[2*notu-3]->v3, prov, nodegc); strcat(ctree,");"); } } /* ctreewrite */ /* Write string describing the s_tree underlying to noeud node at address ctree. */ char* ctreewrite(noeud comesfrom, noeud node, char* ctree, int nodegc){ noeud vd, vg; double l, b; if(comesfrom==node->v1) { vg=node->v2; vd=node->v3; l=node->l1;} else if(comesfrom==node->v2) { vg=node->v1; vd=node->v3; l=node->l2;} else { vg=node->v1; vd=node->v2; l=node->l3;} if(nodegc==1) b=node->c+node->g; else b=node->ceq+node->geq; if(nodegc<0) b=-1.; b*=100.; if (vg==NULL){ sprintf(ctree,"%s", node->nom); while(*ctree) ctree++; // if (b>-0.1) sprintf(ctree," %.2f",b); while(*ctree) ctree++; if (l>-0.5) sprintf(ctree,":%.4f",l); while(*ctree) ctree++; return ctree; } else{ *(ctree++)='('; ctree=ctreewrite(node, vg, ctree, nodegc); } sprintf(ctree,", "); ctree+=2; ctree=ctreewrite(node, vd, ctree, nodegc); *(ctree++)=')'; // if (b>-0.1) sprintf(ctree,"%.2f",b); while(*ctree) ctree++; if (l>-0.9) sprintf(ctree,":%.6f", l); while(*ctree) ctree++; return ctree; } /* stoc */ /* StructTOChaine. Write string describing s_tree arbre_s at address ctree */ /* racine = unrooted/rooted (0/other) , notu = number of leaves. */ void stoc(noeud* arbre_s, int racine, int notu, char* ctree, int nodegc){ char* tring, *prov; prov=ctree; while(*prov) {*prov=0; prov++;} tring=(char*)check_alloc(50*notu,sizeof(char)); if (racine) { tring=ctreewrite(NULL, root, ctree, nodegc); prov=ctree+strlen(ctree)-1; while(*prov!=')') *(prov--)='\0'; strcat(ctree,";"); } else{ tring=ctreewrite(arbre_s[2*notu-3]->v3, arbre_s[2*notu-3], ctree, nodegc); prov=ctree+strlen(ctree)-1; while(*prov!=')') *(prov--)='\0'; *(prov++)=','; *(prov++)=' '; tring=ctreewrite(arbre_s[2*notu-3], arbre_s[2*notu-3]->v3, prov, nodegc); strcat(ctree,");"); } } int retder(int *liste){ int i=0, j; while (liste[i] != 0) i++; j = *(liste + i - 1); *(liste + i - 1) = 0; return j; } void aj(int *liste, int nb) { int i=0; while (liste[i] != 0) i++; *(liste + i) = nb; return; } /* ctob */ /* ChaineTOBits. Read c_tree carbre and write b_tree barbre. */ int ctob(char *carbre, char *nom[], int *barbre[]){ int i=0, j, k, fin=0, nbpo=0, nbpf = 0, cptv1 = 0, nbotu, br_ouverte, *listecour, otu = -1, pomoinspf, cpttree=0, t=1; char c, cc, dejalu = '\0', readname[30]; sscanf(carbre+(cpttree++), "%c", &c); if (c == '[') { while ((c != ']') && c) sscanf(carbre+(cpttree++), "%c", &c); if (c != ']'){ printf("Unmatched '[' ']'\n"); return -1; } } else if (c == '(') cpttree=0; else{ printf("Tree file 1st character must be '(' or '['\n"); return -1; } while ((c != ';') && c ) { sscanf(carbre+(cpttree++), "%c", &c); if (c == '(') nbpo++; if (c == ')') nbpf++; if ((nbpo == nbpf + 1) && (c == ',')) cptv1++; } if (c != ';'){ printf("';' missing at end of tree\n"); return -1; } if (nbpo != nbpf){ printf("Unmatched parenthesis\n"); return -1; } if (cptv1 == 2) /* unrooted : ok */; else if (cptv1 == 1){ /* rooted : problem */ printf("Unexpected rooted tree.\n"); return -1; } else{ /* bad tree string */ printf("Bad number of ',' in tree\n"); return -1; } nbotu=nbpo+2; listecour=(int*) check_alloc(nbotu, sizeof(int)); cpttree=0; sscanf(carbre+(cpttree++), "%c", &c); if (c == '['){ while (c != ']') sscanf(carbre+(cpttree++), "%c", &c); while((c==']') || (c==' ') || (c=='\n') || (c=='\t')) sscanf(carbre+(cpttree++), "%c", &c); if (c!='(') return -1; } else while(c!='(') sscanf(carbre+(cpttree++), "%c", &c); pomoinspf=1; for (i = 0; i < nbotu; i++) listecour[i] = 0; br_ouverte = 0; for (k = 0; t==1; k++) { if (dejalu == '\0') sscanf(carbre+(cpttree++), "%c", &c); else { c = dejalu; dejalu = '\0'; } switch (c) { case ';': fin = 1; break; case ',': case '\n': case '\t': case '\'': case ' ': break; case '(': pomoinspf ++; br_ouverte++; aj(listecour, br_ouverte); break; case ')': pomoinspf--; if (*(carbre+cpttree+1) == ';' || pomoinspf==0) { fin = 1; break; } j = retder(listecour); while(carbre[cpttree]!=',' && carbre[cpttree]!=')') cpttree++; break; default: otu=-1; cc = c; i = 0; while ((cc != ':') && (cc != ',') && (cc != ')') && (cc != '\n') && (cc != ' ')) { if (cc != '\'') { readname[i++] = cc; readname[i]='\0'; } sscanf(carbre+(cpttree++), "%c", &cc); } for(j=0;j r (rooted) or n (not). */ int ctot(char *input, int *arbre[], double *lgbi, double *lgbp, double* bootstrap, char *nom[], char *racine, int* nbbi){ int i=0, j, k, fdf=0, fin=0, nbpo=0, nbpf = 0, cptv1 = 0, nbotu, br_ouverte, *listecour, otu = -1, pomoinspf, cpttree=0, t=1, nbnom=0; char c, cc, cas, dejalu = '\0'; double f; sscanf(input+(cpttree++), "%c", &c); if (c == '[') { while ((c != ']') && (fdf != EOF)) fdf = sscanf(input+(cpttree++), "%c", &c); if (c != ']'){ printf("Unmatched '[' ']'\n"); return -1; } } else if (c == '(') cpttree=0; else{ printf("Tree file 1st character must be '(' or '['\n"); return -1; } while ((c != ';') && (fdf != EOF)) { fdf = sscanf(input+(cpttree++), "%c", &c); if (c == '(') nbpo++; if (c == ')') nbpf++; if ((nbpo == nbpf + 1) && (c == ',')) cptv1++; } if (c != ';'){ printf("';' missing at end of tree file\n"); return -1; } if (nbpo != nbpf){ printf("Unmatched parenthesis\n"); return -1; } if (cptv1 == 1) cas = 'c'; if (cptv1 == 2) cas = 'a'; if ((cptv1!=1) && (cptv1!=2)){ printf("Bad number of ',' in tree file\n"); return -1; } nbotu = nbpo + 2; listecour=(int*) check_alloc(nbotu, sizeof(int)); if (cas == 'a') *racine='n'; else *racine='r'; if (cas=='c'){ if (lgbp) lgbp[0] = 0.; strcpy(nom[0],"ROOT"); for(i=0;iv1 = v1; nd->v2 = v2; nd->v3 = v3; nd->l1 = l1; nd->l2 = l2; nd->l3 = l3; nd->alive1 = (int)b1; nd->alive2 = (int)b2; nd->alive3 = (int)b3; if (nom!=NULL) {strncpy(nd->nom, nom, MAXLNAME); nd->nom[MAXLNAME]='\0'; } return nd; } /* bottomnode */ /* Return the root node of rooted s_tree including node nd. */ /* !!BEWARE!! */ /* If the s_tree including node nd is unrooted, program will bug */ /* (infinite loop). */ noeud bottomnode(noeud nd){ if(nd->v3 == NULL) return nd; return (bottomnode(nd->v3)); } /* ttos */ /* TableTOStruct. Create s_tree arbre_s from : t_tree arbre_t, leaves number notu, */ /* branch lengths lgbi (internal branches) , lgbp (terminal branches). */ /* If tree has no branch length, arguments 3 and 4 must be NULL. */ int ttos(int** arbre_t, int notu, double* lgbi, double* lgbp, double* alive, char** nom, noeud* arbre_s){ noeud p1, p2; int notuv, i, j, k, sommebi, tax1, tax2, n1, n2, n3, cpt1=-1; int *kill_tax, *kill_bi; double arg6, arg9; kill_tax = (int *)check_alloc(notu, sizeof(int)); kill_bi = (int *)check_alloc(notu, sizeof(int)); notuv = notu; for (i = 0; i l3, p2->l3, arg6, (double)p1->alive3, (double)p2->alive3, arg9, NULL); p1->v3 = arbre_s[notu + cpt1]; p2->v3 = arbre_s[notu + cpt1]; kill_tax[tax1] = kill_bi[j] = 1; notuv--; } /* last node */ for (i = 0; i < 2 * notu - 3; i++) if (arbre_s[i]->v3 == NULL) { n1 = i; break; } for (i = n1 + 1; i < 2 * notu - 3; i++) if (arbre_s[i]->v3 == NULL) { n2 = i; break; } for (i = n2 + 1; i < 2 * notu - 3; i++) if (arbre_s[i]->v3 == NULL) { n3 = i; break; } arbre_s[2 * notu - 3] = create_node(arbre_s[n1], arbre_s[n2], arbre_s[n3], arbre_s[n1]->l3, arbre_s[n2]->l3, arbre_s[n3]->l3, arbre_s[n1]->alive3, arbre_s[n2]->alive3, arbre_s[n3]->alive3, NULL); arbre_s[n1]->v3 = arbre_s[2 * notu - 3]; arbre_s[n2]->v3 = arbre_s[2 * notu - 3]; arbre_s[n3]->v3 = arbre_s[2 * notu - 3]; free(kill_tax); free(kill_bi); return 0; } /* mindepth */ /* Set node's depth in a rooted tree. */ /* mindepth(terminal node) = 0 */ /* mindepth(internal node) = number of branches between node and closest terminal node */ /* node comesfrom (one of node's neighbors) is root. */ int mindepth (noeud comesfrom, noeud node){ noeud droit, gauche; if (node->v1==NULL) { node->depth=0; return 0; } if (node->v2==NULL) { node->depth=0; return 0; } if (comesfrom==node->v1) { gauche=node->v2; droit=node->v3; } else if (comesfrom==node->v2) { gauche=node->v1; droit=node->v3; } else { gauche=node->v1; droit=node->v2; } node->depth=mini(mindepth(node, gauche), mindepth(node, droit))+1; return node->depth; } void makelistbr_unrooted(noeud from, noeud nd, branche* br, int* alivebr,int nbbranch){ noeud to1, to2; double lto1, lto2; int i=0, alive1, alive2; if (nd==root){ if(from==root->v1) makelistbr_unrooted(root, nd->v2, br, alivebr, nbbranch); if(from==root->v2) makelistbr_unrooted(root, nd->v1, br, alivebr, nbbranch); return; } if (from==nd->v1) {to1=nd->v2; to2=nd->v3; lto1=nd->l2; lto2=nd->l3; alive1=nd->alive2; alive2=nd->alive3;} else if (from==nd->v2) {to1=nd->v1; to2=nd->v3; lto1=nd->l1; lto2=nd->l3; alive1=nd->alive1; alive2=nd->alive3;} else if (from==nd->v3) {to1=nd->v1; to2=nd->v2; lto1=nd->l1; lto2=nd->l2; alive1=nd->alive1; alive2=nd->alive2;} else { printf("erreur making branch list.\n"); } while(ibout1) i++; if (to1 && to1==root) { if(nd==root->v1){ lto1+=root->l2; br[i]->bout1=nd; br[i]->bout2=root->v2; br[i]->length=lto1; i++; } else{ lto1+=root->l1; br[i]->bout1=nd; br[i]->bout2=root->v1; br[i]->length=lto1; i++; } if(alivebr) alivebr[i-1]=1; } if (to2 && to2==root) { if(nd==root->v1){ lto2+=root->l2; br[i]->bout1=nd; br[i]->bout2=root->v2; br[i]->length=lto2; i++; } else{ lto2+=root->l1; br[i]->bout1=nd; br[i]->bout2=root->v1; br[i]->length=lto2; i++; } if(alivebr) alivebr[i-1]=1; } if (to1 && to1!=root) { br[i]->bout1=nd; br[i]->bout2=to1; br[i]->length=lto1; if(alivebr) alivebr[i]=alive1; i++; } if (to2 && to2!=root) { br[i]->bout1=nd; br[i]->bout2=to2; br[i]->length=lto2; if(alivebr) alivebr[i]=alive2; } if (to1) makelistbr_unrooted(nd, to1, br, alivebr, nbbranch); if (to2) makelistbr_unrooted(nd, to2, br, alivebr, nbbranch); } void organize_tree(noeud from, noeud nd){ noeud prov; double lprov; int alprov; static int numappel=0; if(!nd->v1) return; if(nd!=root){ numappel++; sprintf(nd->nom, "int%d", numappel); } if(from==nd->v1) { prov=nd->v3; nd->v3=nd->v1; nd->v1=prov; lprov=nd->l3; nd->l3=nd->l1; nd->l1=lprov; alprov=nd->alive3; nd->alive3=nd->alive1; nd->alive1=alprov; } if(from==nd->v2) { prov=nd->v3; nd->v3=nd->v2; nd->v2=prov; lprov=nd->l3; nd->l3=nd->l2; nd->l2=lprov; alprov=nd->alive3; nd->alive3=nd->alive2; nd->alive2=alprov; } organize_tree(nd, nd->v1); organize_tree(nd, nd->v2); if(nd==root) numappel=0; } void invert_branch(branche br){ noeud prov; prov=br->bout1; br->bout1=br->bout2; br->bout2=prov; } void organize_listbr_node(noeud nd, branche* listbr, int nbseq){ int i, break1=0, break2=0; if(nd->v1==NULL) return; for(i=0;ibout1==nd->v1 && listbr[i]->bout2==nd)){ invert_branch(listbr[i]); break1=1; } if(!break1 && (listbr[i]->bout1==nd && listbr[i]->bout2==nd->v1)){ break1=1; } if(!break2 && (listbr[i]->bout1==nd->v2 && listbr[i]->bout2==nd)){ invert_branch(listbr[i]); break2=1; } if(!break2 && (listbr[i]->bout1==nd && listbr[i]->bout2==nd->v2)){ break2=1; } if(break1 && break2) break; } organize_listbr_node(nd->v1, listbr, nbseq); organize_listbr_node(nd->v2, listbr, nbseq); } void organize_listbr(tree* s_tree){ int i, nbseq; branche* listbr; listbr=s_tree->listbr; nbseq=s_tree->nbseq; for(i=0;i<2*nbseq-3;i++){ if(listbr[i]->bout1==root->v2 && listbr[i]->bout2==root->v1){ invert_branch(listbr[i]); break; } } organize_listbr_node(root->v1, listbr, nbseq); organize_listbr_node(root->v2, listbr, nbseq); } int setpattern_site(noeud nd, int site) { int i, son1, son2; if(nd->v1==NULL){ if (nd->seq[site]=='A') return 0; else if (nd->seq[site]=='R') return 1; else if (nd->seq[site]=='N') return 2; else if (nd->seq[site]=='D') return 3; else if (nd->seq[site]=='C') return 4; else if (nd->seq[site]=='Q') return 5; else if (nd->seq[site]=='E') return 6; else if (nd->seq[site]=='G') return 7; else if (nd->seq[site]=='H') return 8; else if (nd->seq[site]=='I') return 9; else if (nd->seq[site]=='L') return 10; else if (nd->seq[site]=='K') return 11; else if (nd->seq[site]=='M') return 12; else if (nd->seq[site]=='F') return 13; else if (nd->seq[site]=='P') return 14; else if (nd->seq[site]=='S') return 15; else if (nd->seq[site]=='T') return 16; else if (nd->seq[site]=='W') return 17; else if (nd->seq[site]=='Y') return 18; else if (nd->seq[site]=='V') return 19; else if (nd->seq[site]=='B') return 20; else if (nd->seq[site]=='Z') return 21; return 22; // this is for the cases of X?*- } son1=setpattern_site(nd->v1, site); son2=setpattern_site(nd->v2, site); for(i=0;inbpattern;i++){ if(nd->patternson1[i]==son1 && nd->patternson2[i]==son2) { return i; } } nd->patternson1[nd->nbpattern]=son1; nd->patternson2[nd->nbpattern]=son2; nd->nbpattern++; return nd->nbpattern-1; } void setpattern(tree* s_tree){ int i, nbnd; static int nbappel=0; nbappel++; nbnd=2*s_tree->nbseq-2; for(i=0;inode[i]->nbpattern=0; root->nbpattern=0; for(i=0;ilgseq;i++) setpattern_site(root, i); } int setunderlyingbr(noeud nd, branche br, int i){ int value=0; if(nd->v1==NULL){ nd->underlyingbr[i]=0; return 0; } if(setunderlyingbr(nd->v1, br, i)) value=1; else if(setunderlyingbr(nd->v2, br, i)) value=2; else if(nd==br->bout1 && nd->v1==br->bout2) value=10; else if(nd==br->bout2 && nd->v1==br->bout1) value=10; else if(nd==br->bout1 && nd->v2==br->bout2) value=20; else if(nd==br->bout2 && nd->v2==br->bout1) value=20; else if(nd->v1==br->bout1 && nd->v2==br->bout2) value=12; /* nd=root, br=root branch */ else if(nd->v2==br->bout1 && nd->v1==br->bout2) value=12; /* nd=root, br=root branch */ nd->underlyingbr[i]=value; return value; } void setunderlyinggc(noeud nd, int nb){ int i; for(i=0;i<2*nb-3;i++){ nd->underlyinggc[i]=nd->underlyingbr[i]; if(nd->underlyinggc[i]==12) nd->underlyinggc[i]=10; } if(nd==root) nd->underlyinggc[2*nb-3]=20; else nd->underlyinggc[2*nb-3]=0; if(nd->v1==NULL) return; setunderlyinggc(nd->v1, nb); setunderlyinggc(nd->v2, nb); } void set_numerobr(noeud nd, branche* listbr, int nbbr){ int i; noeud br1, br2; if(nd==root) nd->numerobr=-1; else{ for(i=0;ibout1; br2=listbr[i]->bout2; if(nd->v3==root){ if((root->v1==br1 && root->v2==br2)||(root->v2==br1 && root->v1==br2)){ nd->numerobr=i; break; } } else if((nd->v3==br1 && nd==br2)||(nd==br1 && nd->v3==br2)){ nd->numerobr=i; break; } } } if(!nd->v1) return; set_numerobr(nd->v1, listbr, nbbr); set_numerobr(nd->v2, listbr, nbbr); } void root_s(branche br, double frac1){ noeud nd1, nd2; double l; nd1=br->bout1; nd2=br->bout2; l=br->length; root->v1=nd1; root->v2=nd2; root->l1=l*frac1; root->l2=l*(1-frac1); root->alive1=root->alive2=1; root->alive3=-1; nd1->remain=nd2; nd2->remain=nd1; nd1->lremain=l; nd2->lremain=l; if (nd1->v1==nd2) { nd1->v1=root; nd1->l1=l*frac1; nd1->alremain=nd1->alive1; nd1->alive1=1;} else if (nd1->v2==nd2) { nd1->v2=root; nd1->l2=l*frac1; nd1->alremain=nd1->alive2; nd1->alive2=1; } else if (nd1->v3==nd2) { nd1->v3=root; nd1->l3=l*frac1; nd1->alremain=nd1->alive3; nd1->alive3=1; } else printf("Error rooting.\n"); if (nd2->v1==nd1) { nd2->v1=root; nd2->l1=l*(1-frac1); nd2->alremain=nd2->alive1; nd2->alive1=1; } else if (nd2->v2==nd1) { nd2->v2=root; nd2->l2=l*(1-frac1); nd2->alremain=nd2->alive2; nd2->alive2=1; } else if (nd2->v3==nd1) { nd2->v3=root; nd2->l3=l*(1-frac1); nd2->alremain=nd2->alive3; nd2->alive3=1; } else printf("Error rooting.\n"); } /* unroot */ /* Unroot t_tree treer. All parameters are modified (branch lengths, names, ...). */ /* Root info may be kept in list1-list2 l1-l2 fracroot1 if l1 is non NULL and list1 and list2 */ /* are allocated as nbtaxa-sized char* arrays. */ int unroot(int** treer, int notu, double* lgbi, double* lgbp, double* bootvals, char** nom, char** list1, char** list2, int* l1, int* l2, double* fracroot1){ int i, j, j1, j2, br_a_virer=-1, bp_a_modifier=-1, bi_a_modifier=-1, somme1, somme2, drapeau; double l; for(j=0;jv1) { n1=nd->v2; n2=nd->v3; } else if(from==nd->v2) { n1=nd->v1; n2=nd->v3; } else if(from==nd->v3) { n1=nd->v1; n2=nd->v2; } else { printf("Erreur\n"); } if(!n1){ for(i=0;inom, list[i])) {flag=1; break; } } return flag; } return(equallist(nd, n1, nb, list) && equallist(nd, n2, nb, list)); } /* isbranch */ /* Check wether lists of taxa name list1 and list2 (respectively l1 */ /* and l2 taxa) exactly match the lists of taxa laying left-side */ /* and right-side branch br. */ int isbranch(branche br, int l1, char** list1, int l2, char** list2){ if( equallist(br->bout1, br->bout2, l1, list1) && equallist(br->bout2, br->bout1, l2, list2)) return 2; if( equallist(br->bout2, br->bout1, l1, list1) && equallist(br->bout1, br->bout2, l2, list2)) return 1; return 0; } /* taxa_count */ /* Return the number of taxa in c_tree ctree. */ int taxa_count(char* ctree){ char* prov; int nbpo=0, nbpf=0, nbv=0; prov=ctree; while(*prov!='[' && *prov!='(' && *prov) prov++; if(*prov=='\0'){ printf("Trees should begin with '[' or '('\n"); return -1; } if(*prov=='['){ while(*prov!=']') prov++; while(*prov!='(') prov++; } /* debut arbre */ while(*prov && *prov!=';'){ if(*prov=='(') nbpo++; if(*prov==')') nbpf++; if(nbpo-nbpf==1 && *prov==',') nbv++; prov++; } if(nbv==1){ /* rooted */ return (nbpo+1); } if(nbv==2){ /* unrooted */ return (nbpo+2); } printf("Bad number of commas ',' in tree\n"); return -1; } noeud* ctos(char* ctree, int* notu, int *tree_is_rooted){ int i, nb, **ttree, l1, l2, nbbi; double *lgbp, *lgbi, *bootvals, fracroot1; char **nom, racine, *list1[MAXNSP], *list2[MAXNSP]; noeud* stree; branche* br_list; nb=taxa_count(ctree); lgbp=(double*)check_alloc(nb+1, sizeof(double)); lgbi=(double*)check_alloc(nb+1, sizeof(double)); bootvals=(double*)check_alloc(nb+1, sizeof(double)); /*list1=check_alloc(2*nb+2, sizeof(char*)); list2=check_alloc(2*nb+2, sizeof(char*));*/ nom=check_alloc(nb+1, sizeof(char*)); ttree=check_alloc(nb+1, sizeof(int**)); for(i=0;i0.) *sce+=(pdist[i]-cdist[i])*(pdist[i]-cdist[i]); else *sce+=cdist[i]*cdist[i]; } if(totallength){ *totallength=0.; for(i=0;i0.) *totallength+=lgbr[i]; } free(cdist); free(ind1); free(ind2); for(i=0;iv1) return congruent(root->v2, root, btree, name, nbseq, j); else if(from==root->v2) return congruent(root->v1, root, btree, name, nbseq, j); } if(nd->v1==NULL){ for(i=0;inom, name[i])){ if(testbit(btree[j], i+1)) return 1; else return 0; } } if(i==nbseq){ printf("erreur cong\n"); exit(0); } } if(from==nd->v1) {to1=nd->v2; to2=nd->v3;} else if(from==nd->v2) {to1=nd->v1; to2=nd->v3;} else if(from==nd->v3) {to1=nd->v1; to2=nd->v2;} else if(from==root->v1 && nd==root->v2) {to1=nd->v1; to2=nd->v2;} else if(from==root->v2 && nd==root->v1) {to1=nd->v1; to2=nd->v2;} else { printf("erreur cong2 : nd=%s , from=%s\n", nd->nom, from->nom); exit(0);} c1=congruent(to1, nd, btree, name, nbseq, j); c2=congruent(to2, nd, btree, name, nbseq, j); if(c1==c2) return c1; return 10; } void reset_bl(tree* s_tree, noeud nd){ double l; if(nd==NULL) return; if(nd!=root){ l=s_tree->listbr[nd->numerobr]->length; if(nd==root->v1) nd->l3=root->l1=l*s_tree->froot; else if(nd==root->v2) nd->l3=root->l2=l*(1.-s_tree->froot); else{ nd->l3=l; if(nd==nd->v3->v1) nd->v3->l1=l; if(nd==nd->v3->v2) nd->v3->l2=l; } } reset_bl(s_tree, nd->v1); reset_bl(s_tree, nd->v2); } void set_parambl(tree* s_tree, noeud nd){ double l; if(nd==NULL) return; if(nd!=root){ if(nd->v3==root) l=root->l1+root->l2; else l=nd->l3; s_tree->listbr[nd->numerobr]->length=l; } set_parambl(s_tree, nd->v1); set_parambl(s_tree, nd->v2); } void btltostl(tree* stree, int** btree, char** name, double* treelg){ int i, j, nbseq, left, right; branche br; noeud boutterm; nbseq=stree->nbseq; for(i=0;i<2*nbseq-3;i++){ br=stree->listbr[i]; if(br->bout1->v1==NULL || br->bout2->v1==NULL){ /* terminal */ boutterm=(br->bout1->v1)?br->bout2:br->bout1; for(j=0;jnom, name[j])){ br->length=treelg[j]; break; } } if(j==nbseq) printf("erreur1\n"); } else{ /* internal */ for(j=0;jbout1, br->bout2, btree, name, nbseq, j); right=congruent(br->bout2, br->bout1, btree, name, nbseq, j); if((left==0 && right==1) || (left==1 && right==0)){ br->length=treelg[nbseq+j]; break; } } if(j==nbseq-3) printf("erreur2\n"); } } reset_bl(stree, root); } /* unroot_c */ /* Remove root from rooted c_tree carbre. */ int unroot_c(char* carbre){ int i=0, diff; while((carbre[i]==' ' || carbre[i]=='\n' || carbre[i]=='\t') && carbre[i]) i++; if(carbre[i]!='[' && carbre[i]!='(') { printf("Tree first char must be ( or [\n"); return -1; } if(carbre[i]=='[') while(carbre[i]!=']' && carbre[i]) i++; if(!carbre[i]){ printf("Unmatched '[' ']'\n"); return -1; } while(carbre[i]!='(' && carbre[i]) i++; if(!carbre[i]) {printf("No initial parenthesis\n"); return -1; } i++; while(carbre[i]!='(' && carbre[i] && carbre[i]!=';') i++; if(!carbre[i] || carbre[i]==';') return 0; carbre[i]=' '; diff=0; while(diff!=-1 && carbre[i] && carbre[i]!=';'){ if(carbre[i]=='(') diff++; if(carbre[i]==')') diff--; i++; } if(!carbre[i] || carbre[i]==';') return 0; carbre[i-1]=' '; while(carbre[i]!=',' && carbre[i]!=')') { carbre[i]=' '; i++; } return 1; } void init_bl(tree* stree){ int i, nbseq, lliste, **btree; char* ctree; double* newbl, sumalive=0., sum=0., rap; nbseq=stree->nbseq; ctree=(char*)check_alloc(50*nbseq, sizeof(char)); lliste=(nbseq+lmot-1)/lmot; btree=(int**)check_alloc(nbseq-3, sizeof(int*)); for(i=0;inode, 1, nbseq, ctree, 0); unroot_c(ctree); ctob(ctree, stree->names, btree); newbl=lslgbr(btree, nbseq, stree->dist, NULL, NULL); for(i=0;i<2*nbseq-3;i++){ if(newbl[i]lmax) newbl[i]=lmax; } for(i=0;i<2*nbseq-3;i++){ sum+=newbl[i]; if(stree->alivebr[i]) sumalive+=newbl[i]; } rap=sum/sumalive; for(i=0;i<2*nbseq-3;i++) if(stree->alivebr[i]) newbl[i]*=rap; btltostl(stree, btree, stree->names, newbl); for(i=0;il1 , nd->l2 */ /* root location : nd->l1/(nd->l1+nd->l2) */ double rr1, rr2; /* relative rate */ double p0, p1, p2, q; long i,j,k,l,m,kk; int parameters; double elambdat[20],delambdat[20],ddelambdat[20]; if(nd->v1==NULL) return; tdeb=time(NULL); /* calculs des rr (=rij) */ effective_length(nd->l1, cl1, cl2, covar3, &(nd->rr1[cl1][cl2]), &(nd->drr1[cl1][cl2]), &(nd->d2rr1[cl1][cl2])); effective_length(nd->l2, cl1, cl2, covar3, &(nd->rr2[cl1][cl2]), &(nd->drr2[cl1][cl2]), &(nd->d2rr2[cl1][cl2])); rr1=nd->rr1[cl1][cl2]; rr2=nd->rr2[cl1][cl2]; /* proba. to branch 1 */ for (k=0;k<20;k++) { elambdat[k] = exp(nd->l1*rr1*eigmat[k]); delambdat[k] = (elambdat[k]*rr1*eigmat[k]); ddelambdat[k] = (delambdat[k]*rr1*eigmat[k]); } for (m=0;m<20;m++) { for (l=0;l<20;l++) { p0=0.0; p1=0.0; p2=0.0; for (k=0;k<20;k++) { q=probmat[k][m]*probmat[k][l]; p0 += (q*elambdat[k]); p1 += (q*delambdat[k]); p2 += (q*ddelambdat[k]); } if(nd->alive1) { if(p0<0.) p0=0.; // why negative happens ??? nd->p1[cl1][cl2][m][l]=p0/freqaa[m]; for (parameters=0;parameters<3;parameters++){ nd->dp1[cl1][cl2][parameters][m][l]=p1/freqaa[m]; nd->d2p1[cl1][cl2][parameters][m][l]=p2/freqaa[m]; } } } } if(nd->p1[cl1][cl2][0][1]<0.){ printf("l1=%f, rr1=%f\n", nd->l1, rr1); printf("cl1=%d, cl2=%d\n", cl1, cl2); printf("alive1: %d, alive2: %d\n", nd->alive1, nd->alive2); } /* proba. to branch 2 */ for (k=0;k<20;k++) { elambdat[k] = exp(nd->l2*rr2*eigmat[k]); delambdat[k] = (elambdat[k]*rr2*eigmat[k]); ddelambdat[k] = (delambdat[k]*rr2*eigmat[k]); } for (m=0;m<20;m++) { for (l=0;l<20;l++) { p0=0.0; p1=0.0; p2=0.0; for (k=0;k<20;k++) { q=probmat[k][m]*probmat[k][l]; p0 += (q*elambdat[k]); p1 += q*delambdat[k]; p2 += q*ddelambdat[k]; } if(nd->alive2){ if (p0<0.) p0=0.; nd->p2[cl1][cl2][m][l]=p0/freqaa[m]; for (parameters=0;parameters<3;parameters++){ nd->dp2[cl1][cl2][parameters][m][l]=p1/freqaa[m]; nd->d2p2[cl1][cl2][parameters][m][l]=p2/freqaa[m]; } } } } /* IF UNRESOLVED NODE */ if(nd->alive1==0){ for(i=0;i<20;i++) for(j=0;j<20;j++) if(i==j) nd->p1[cl1][cl2][i][j]=1.; else nd->p1[cl1][cl2][i][j]=0.; for(i=0;i<20;i++) for(j=0;j<20;j++) for(kk=0;kk<20;kk++) nd->dp1[cl1][cl2][i][j][kk]=0; for(i=0;i<20;i++) for(j=0;j<20;j++) for(kk=0;kk<20;kk++) nd->d2p1[cl1][cl2][i][j][kk]=0; } if(nd->alive2==0){ for(i=0;i<20;i++) for(j=0;j<20;j++) if(i==j) nd->p2[cl1][cl2][i][j]=1.; else nd->p2[cl1][cl2][i][j]=0.; for(i=0;i<20;i++) for(j=0;j<20;j++) for(kk=0;kk<20;kk++) nd->dp2[cl1][cl2][i][j][kk]=0; for(i=0;i<20;i++) for(j=0;j<20;j++) for(kk=0;kk<20;kk++) nd->d2p2[cl1][cl2][i][j][kk]=0; } tfin=time(NULL); probatime+=difftime(tfin, tdeb); compute_proba_t3_cov(nd->v1, OPTIMIZE_LENGTH, OPTIMIZE_ROOT, OPTIMIZE_COV, cl1, cl2); compute_proba_t3_cov(nd->v2, OPTIMIZE_LENGTH, OPTIMIZE_ROOT, OPTIMIZE_COV, cl1, cl2); } int callMyEigen(int n, double **R, double **vr, double *rr) { int i, ret; double *ri, **vi; vi=check_alloc(n, sizeof(double*)); for(i=0;iv1==NULL) return 1; /* bl1: compute P */ rtop(dim, nd->l1, P, vr, invvr, rr); /* bl1: change to nhml proba */ for(m1=0;m1p1[ri][rf][i][j]=P[m1][m2]; if(P[m1][m2]<0.){ if(P[m1][m2]<-0.0000001) return 0; else nd->p1[ri][rf][i][j]=0.; } if(P[m1][m2]>1.){ if(P[m1][m2]>1.0000001) return 0; else nd->p1[ri][rf][i][j]=1.; } } } /* bl2: compute P */ rtop(dim, nd->l2, P, vr, invvr, rr); /* bl2: change to nhml proba */ for(m1=0;m1p2[ri][rf][i][j]=P[m1][m2]; if(P[m1][m2]<0.){ if(P[m1][m2]<-0.0000001) return 0; // prob. cannot be computed else nd->p2[ri][rf][i][j]=0.; } if(P[m1][m2]>1.){ if(P[m1][m2]>1.0000001) return 0; else nd->p2[ri][rf][i][j]=1.; } } } ret1=compute_exact_probacov(nd->v1, vr, invvr, rr); ret2=compute_exact_probacov(nd->v2, vr, invvr, rr); return ret1*ret2; /* return 1 is good */ } void reset_tree(tree* s_tree){ if(s_tree->which_var_param[1]>=0 || s_tree->which_var_param[0]>=0) reset_bl(s_tree, root); // if(s_tree->which_var_param[4]>=0) reset_gc(s_tree, root); // if(s_tree->which_var_param[2]>=0) set_node_comp(root, s_tree->GCanc, s_tree); } void boundaries(tree* s_tree, int* bounded) { int i, limit; limit=s_tree->nbalivebr; /* limit=2*s_tree->nbseq-3; */ for(i=0;inb_var_param;i++) bounded[i]=0; if(s_tree->which_var_param[0]>=0){ if(*(s_tree->var_param[s_tree->which_var_param[0]])var_param[s_tree->which_var_param[0]])=frmin; bounded[s_tree->which_var_param[0]]=1; } if(*(s_tree->var_param[s_tree->which_var_param[0]])>frmax){ *(s_tree->var_param[s_tree->which_var_param[0]])=frmax; bounded[s_tree->which_var_param[0]]=1; } } if(s_tree->which_var_param[1]>=0){ for(i=0;ivar_param[s_tree->which_var_param[1]+i])var_param[s_tree->which_var_param[1]+i])=lmin; bounded[s_tree->which_var_param[1]+i]=1; } } if(s_tree->which_var_param[2]>=0){ if(*(s_tree->var_param[s_tree->which_var_param[2]])var_param[s_tree->which_var_param[2]])=covmin; bounded[s_tree->which_var_param[2]]=1; } if(*(s_tree->var_param[s_tree->which_var_param[2]])>covmax){ *(s_tree->var_param[s_tree->which_var_param[2]])=covmax; bounded[s_tree->which_var_param[2]]=1; } } if(s_tree->which_var_param[3]>=0){ if(*(s_tree->var_param[s_tree->which_var_param[3]])var_param[s_tree->which_var_param[3]])=covmin; bounded[s_tree->which_var_param[3]]=1; } if(*(s_tree->var_param[s_tree->which_var_param[3]])>covmax){ *(s_tree->var_param[s_tree->which_var_param[3]])=covmax; bounded[s_tree->which_var_param[3]]=1; } } if(s_tree->which_var_param[4]>=0){ if(*(s_tree->var_param[s_tree->which_var_param[4]])var_param[s_tree->which_var_param[4]])=covmin; bounded[s_tree->which_var_param[4]]=1; } if(*(s_tree->var_param[s_tree->which_var_param[4]])>covmax){ *(s_tree->var_param[s_tree->which_var_param[4]])=covmax; bounded[s_tree->which_var_param[4]]=1; } } if(s_tree->which_var_param[5]>=0){ if(*(s_tree->var_param[s_tree->which_var_param[5]])var_param[s_tree->which_var_param[5]])=pimin; bounded[s_tree->which_var_param[5]]=1; } if(*(s_tree->var_param[s_tree->which_var_param[5]])>pimax){ *(s_tree->var_param[s_tree->which_var_param[5]])=pimax; bounded[s_tree->which_var_param[5]]=1; } } } void resetdx(noeud nd, int param){ int k, i, cl; if(nd->v1==NULL) return; for(cl=0;clnbpattern;i++){ for(k=0;k<20;k++){ nd->dx[cl][i][k]=0.; } } } resetdx(nd->v1, param); resetdx(nd->v2, param); } void resetd2x(noeud nd){ int k, i, cl; if(nd->v1==NULL) return; for(cl=0;clnbpattern;i++){ for(k=0;k<20;k++){ nd->d2x[cl][i][k]=0.; } } } resetd2x(nd->v1); resetd2x(nd->v2); } void reset_like(tree* s_tree){ int i, j, k, l, nbpattern; for(i=s_tree->nbseq;i<2*s_tree->nbseq-2;i++){ nbpattern=s_tree->node[i]->nbpattern; for(j=0;jnode[i]->x[l][j][k]=-1.; } } } } } void dlike_dpi(tree* s_tree, int param, double* dlike_s, int cross){ int i; for(i=0;ilgseq;i++) dlike_s[i]=s_tree->like_s[i]-s_tree->likeasrv_s[i]; } void printMatrix(double **matrx, int nbparam1, int offset, FILE *out) { int i, j; for(i=offset;inb_var_param;i++){ switch(s_tree->param_nature[i]){ case COVA: fprintf (out, "cov1\t"); break; case COVB: fprintf (out, "cov2\t"); break; case COVC: fprintf (out, "cov3\t"); break; case GaltierPI: fprintf (out, "GaltierPi\t"); break; } } if (compute_o->OPTIMIZE_GAMMA) fprintf(out, "alpha"); fprintf (out, "\n"); for (i=0; inb_var_param;i++){ switch(s_tree->param_nature[i]){ case ROOT: case LENGTH: offset ++; break; case COVA: case COVB: case COVC: case GaltierPI: nbparam++; break; } } if (compute_o->OPTIMIZE_GAMMA) nbparam1=nbparam+1; else nbparam1=nbparam; remparam = (double *) check_alloc(s_tree->nb_var_param, sizeof(double)); useless=(double*)check_alloc(gamma_nbcl, sizeof(double)); lkh1 = (double *) check_alloc(nbparam1, sizeof(double)); covarianceM = (double **)check_alloc(nbparam1, sizeof(double *)); fisherI = (double **)check_alloc(nbparam1, sizeof(double *)); for (i=0; inb_var_param;i++) remparam[i]=*(s_tree->var_param[i]); /* store current log likelihoods */ l0 = s_tree->lkh; for (i=0; ivar_param[i+offset])= remparam[i+offset]; for(i=0;inb_var_param;i++) remparam[i]=*(s_tree->var_param[i]); if (compute_o->OPTIMIZE_GAMMA) { for (i=0; ivar_param[i+offset])= remparam[i+offset]; } /* get diagonal entry for alpha in the fisher information matrix */ gamma_shp += h2; DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); l1 = -minus_log_likelihood(remparam); gamma_shp -= 2*h2; DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); l4 = -minus_log_likelihood(remparam); fisherI[nbparam][nbparam] = -1*(l1+l4-2*l0)/h2square; gamma_shp += h2; DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); } s_tree->lkh = l0; fprintf(out, "fisher information matrix\n"); printMatrix(fisherI, nbparam1, offset, out); /* inverse fisher information matrix to get covariance matrix if ret = 0 if the fisher information matrix is singular */ ret = invmat(fisherI, nbparam1, covarianceM); if (ret == 0) fprintf (out, "Singular fisher information matrix! Covariance matrix cannot be computed.\n"); else { fprintf(out, "covariance matrix\n"); printMatrix (covarianceM, nbparam1, offset, out); } for (i=0; inb_var_param; lgseq=s_tree->lgseq; resetdx(root, 0); h = 0.00001; hsquare=h*h; // store likelihood and site likelihood for the tree of current parameters l0 = s_tree->lkh; like_s = (double *) check_alloc(lgseq, sizeof(double)); likeasrv_s = (double *) check_alloc(lgseq, sizeof(double)); for (i=0; ilike_s[i]; likeasrv_s[i] = s_tree->likeasrv_s[i]; } remparam = (double *) check_alloc(nbparam, sizeof(double)); for(i=0;ivar_param[i]); for (i=0; iparam_nature[i] != GaltierPI) { // we do analytical deriv for PI_galtier remparam[i] += h; // x+h l1 = -minus_log_likelihood(remparam); // f(x+h) remparam[i] += h; // x+2h l2 = -minus_log_likelihood(remparam); // f(x+2h) remparam[i] -= 2*h; // restore the original parameter *(s_tree->var_param[i])= remparam[i]; s_tree->deriv[i] = (l1-l0)/h; s_tree->deriv2[i] = (l2+l0-2*l1)/hsquare; } else { dlike_dpi(s_tree, i, s_tree->dlike_s[i],0); s_tree->deriv[i]=s_tree->deriv2[i]=0.; for(j=0;jdlike_s[i][j]; d2like = 0; s_tree->deriv[i]+=s_tree->weight[j]*dlike/like; s_tree->deriv2[i]+=s_tree->weight[j]*(-(dlike*dlike/(like*like))); } } } // restore tree likelihood and site likelihoods to the original value before calling the function s_tree->lkh = l0; for (i=0; ilike_s[i] = like_s[i]; s_tree->likeasrv_s[i] = likeasrv_s[i]; } free(remparam); free(like_s);free(likeasrv_s); } /* like_node */ /* compute likelihood for node nd */ void ASRVlike_node(noeud nd){ int i, j, k, cl; double s1, s2; if(nd->v1==NULL) return; ASRVlike_node(nd->v1); ASRVlike_node(nd->v2); for(i=0;inbpattern;i++){ for(cl=0;clp1[cl][cl][j][k]*nd->v1->x[cl][nd->patternson1[i]][k]; s2+=nd->p2[cl][cl][j][k]*nd->v2->x[cl][nd->patternson2[i]][k]; } nd->x[cl][i][j]=s1*s2; if(nd->x[cl][i][j]<0.){ if(nd->x[cl][i][j]x[cl][i][j]:%.20f\n", i, j, nd->x[cl][i][j]); exit(0); } else nd->x[cl][i][j]=0.; } } } } } void Covarionlike_node(noeud nd){ int i, j, k, cl, cl2, deb, fin; double s1, s2; if(nd->v1==NULL) return; Covarionlike_node(nd->v1); Covarionlike_node(nd->v2); for(i=0;inbpattern;i++){ for(cl=0;clp1[cl][cl2][j][k]*nd->v1->x[cl2][nd->patternson1[i]][k]; s2+=nd->p2[cl][cl2][j][k]*nd->v2->x[cl2][nd->patternson2[i]][k]; } } nd->x[cl][i][j]=s1*s2; if(nd->x[cl][i][j]<0.){ if(nd->x[cl][i][j]x[cl][i][j]:%.20f\n", i, j, nd->x[cl][i][j]); exit(0); } else nd->x[cl][i][j]=0.; } } } } } /* like_node_2sites */ /* Compute joint likelihood of the current two sites for node nd */ void like_node_2sites(noeud nd){ int i, j, i1, j1, cl, cl1, deb, fin; double s1, s2, t1, t2; if(nd->v1==NULL){ return; } like_node_2sites(nd->v1); like_node_2sites(nd->v2); for(cl=0;cl0.){deb=0; fin=gamma_nbcl;} else {deb=cl; fin=cl+1;} for(i=0;i<20;i++){ for(j=0;j<20;j++){ s1=s2=0.; for(cl1=deb;cl1v1->x2sites[cl1][i1][j1]*nd->p1[cl][cl1][i][i1]*nd->p1[cl][cl1][j][j1]; t2=nd->v2->x2sites[cl1][i1][j1]*nd->p2[cl][cl1][i][i1]*nd->p2[cl][cl1][j][j1]; s1+=t1; s2+=t2; } } } nd->x2sites[cl][i][j]=s1*s2; if(nd->x2sites[cl][i][j]<0.){ printf("i:%d , j:%d , nd->x2sites[cl][i][j]:%.20f\n", i, j, nd->x2sites[cl][i][j]); exit(0); } } } } } /* init_x2sites */ /* Initialize x2sites at tips for sites s1 and s2 */ void init_x2sites(noeud nd, int s1, int s2){ int j, k, l, p1, p2; int aa, aa1; if(nd->v1){ init_x2sites(nd->v1, s1, s2); init_x2sites(nd->v2, s1 ,s2); return;} p1=setpattern_site(nd, s1); p2=setpattern_site(nd, s2); printf("p1=%d\tp2=%d\n",p1,p2); for(j=0;jx2sites[j][k][l]=1.; else nd->x2sites[j][k][l]=0.; } } } // special cases if (p1==20) { // B if (l==p2) {nd->x2sites[j][2][l]=1.; nd->x2sites[j][3][l]=1.;} if (p2==20) { nd->x2sites[j][2][2]=1.; nd->x2sites[j][2][3]=1.; nd->x2sites[j][3][2]=1.; nd->x2sites[j][3][3]=1.; } if (p2==21) { nd->x2sites[j][2][5]=1.; nd->x2sites[j][2][6]=1.; nd->x2sites[j][3][5]=1.; nd->x2sites[j][3][6]=1.; } if (p2==22) { for (aa=0;aa<20;aa++) {nd->x2sites[j][2][aa]=1.; nd->x2sites[j][3][aa]=1.;} } } if (p1==21) { // Z if (l==p2) {nd->x2sites[j][5][l]=1.; nd->x2sites[j][6][l]=1.;} if (p2==20) { nd->x2sites[j][5][2]=1.; nd->x2sites[j][5][3]=1.; nd->x2sites[j][6][2]=1.; nd->x2sites[j][6][3]=1.; } if (p2==21) { nd->x2sites[j][5][5]=1.; nd->x2sites[j][5][6]=1.; nd->x2sites[j][6][5]=1.; nd->x2sites[j][6][6]=1.; } if (p2==22) { for (aa=0;aa<20;aa++) {nd->x2sites[j][5][aa]=1.; nd->x2sites[j][6][aa]=1.;} } } if (p1==22) { // gaps, unknnwn residues if (l==p2) { for (aa=0;aa<20;aa++) nd->x2sites[j][aa][l]=1.; } if (p2==20) { for (aa=0;aa<20;aa++) { nd->x2sites[j][aa][2]=1.; nd->x2sites[j][aa][3]=1.; } } if (p2==21) { for (aa=0;aa<20;aa++) { nd->x2sites[j][aa][5]=1.; nd->x2sites[j][aa][6]=1.; } } if (p2==22) { for (aa=0;aa<20;aa++) for (aa1=0;aa1<20;aa1++) nd->x2sites[j][aa][aa1]=1.; } } if (p2==20) { // B if (k==p1) {nd->x2sites[j][k][2]=1.; nd->x2sites[j][k][3]=1.;} if (p1==21) { nd->x2sites[j][5][2]=1.; nd->x2sites[j][5][3]=1.; nd->x2sites[j][6][2]=1.; nd->x2sites[j][6][3]=1.; } if (p1==22) { for (aa=0;aa<20;aa++) {nd->x2sites[j][aa][2]=1.; nd->x2sites[j][aa][3]=1.;} } } if (p2==21) { // Z if (k==p1) {nd->x2sites[j][k][5]=1.; nd->x2sites[j][k][6]=1.;} if (p1==20) { nd->x2sites[j][2][5]=1.; nd->x2sites[j][2][6]=1.; nd->x2sites[j][3][5]=1.; nd->x2sites[j][3][6]=1.; } if (p1==22) { for (aa=0;aa<20;aa++) {nd->x2sites[j][aa][5]=1.; nd->x2sites[j][aa][6]=1.;} } } if (p2==22) { // gaps, unknnwn residues if (k==p1) { for (aa=0;aa<20;aa++) nd->x2sites[j][k][aa]=1.; } if (p1==20) { for (aa=0;aa<20;aa++) { nd->x2sites[j][2][aa]=1.; nd->x2sites[j][3][aa]=1.; } } if (p1==21) { for (aa=0;aa<20;aa++) { nd->x2sites[j][5][aa]=1.; nd->x2sites[j][6][aa]=1.; } } } } /* like_2sites */ /* Compute the joint likelihood of sites s1 and s2 assuming */ /* common changes of rate class */ /* note: patterns and patternsons are not used for this calculation */ double like_2sites(tree* s_tree, int s1, int s2){ int i, j, k; double like; init_x2sites(root, s1, s2); like_node_2sites(root); like=0.; for(k=0;kx2sites[k][i][j]*freqaa[i]*freqaa[j]; } if(like<0.){ printf("x2sites %d %d: like <0\n", s1, s2); return 1.; } like/=(gamma_nbcl*gamma_nbcl); return like; } double ASRVsiteLikelihood (tree* s_tree) { int i, j, k; double memcovar; // if(covar3>=0. && pi>=0.){ /* remove the condition 23/03/2006*/ if (1) { memcovar=covar3; covar3 = -1; for(k=0;kOPTIMIZE_LENGTH, compute_o->OPTIMIZE_ROOT, compute_o->OPTIMIZE_COV3, k, k); ASRVlike_node(root); for(k=0;klgseq;i++){ s_tree->classlike_s[k][i]=0; for(j=0;j<20;j++) s_tree->classlike_s[k][i] +=root->x[k][i][j]*freqaa[j]; if(s_tree->classlike_s[k][i]<0.){ printf("site : %d , argument=%e\n", i, s_tree->classlike_s[k][i]); printf("%d : %e %e %e %e\n", i, root->x[k][i][0], root->x[k][i][1], root->x[k][i][2], root->x[k][i][3]); return 1.; // site likehood cannot be computed } } } for(i=0;ilgseq;i++){ s_tree->likeasrv_s[i]=0.; for(k=0;klikeasrv_s[i]+=s_tree->classlike_s[k][i]; s_tree->likeasrv_s[i]/=gamma_nbcl; } covar3=memcovar; } return 0.; // site likehood computed } double CovarionSiteLikelihood (tree* s_tree) { int m, dim, ret, i, j, k; double pi_on, **vr, **invvr, *rr; m = 20; // row size of amino acid subst. matrix dim = m*gamma_nbcl1; rr = check_alloc(dim, sizeof(double)); vr = check_alloc(dim, sizeof(double*)); invvr = check_alloc(dim, sizeof(double*)); for (i=0; ilgseq;i++){ s_tree->classlike_s[k][i]=0; for(j=0;j<20;j++) s_tree->classlike_s[k][i] +=root->x[k][i][j]*freqaa[j]; if(s_tree->classlike_s[k][i]<0.){ printf("site : %d , argument=%e\n", i, s_tree->classlike_s[k][i]); printf("%d : %e %e %e %e\n", i, root->x[k][i][0], root->x[k][i][1], root->x[k][i][2], root->x[k][i][3]); return 1.; // site likehood cannot be computed } } } for(i=0;ilgseq;i++){ s_tree->like_s[i]=0.; for(k=0;klike_s[i]+=s_tree->classlike_s[k][i]; pi_on = covar1/(covar1+covar2); s_tree->like_s[i] *= pi_on; for(k=gamma_nbcl;klike_s[i]+=s_tree->classlike_s[k][i]*(1-pi_on); s_tree->like_s[i]/=gamma_nbcl; } return 0.; // site likehood computed } /* minus_log_likelihood Compute minus the log-likelihood of tree s_tree (global) */ /* for parameters paramvprov and Gamma distribution given by */ /* gamma_nbcl and gamma_clmean (global) */ double minus_log_likelihood(double *paramprov){ int i, nbparam, ret; double l=0; int *bounded; nbparam=s_tree->nb_var_param; bounded = (int *) check_alloc(nbparam,sizeof(int)); for(i=0;ivar_param[i])=paramprov[i]; bounded[i] = 0; } boundaries(s_tree, bounded); /* let the parameters bounded in the preset lower and upper bounds: this only works for those parameters set to be optimized, as defined in the option file */ reset_tree(s_tree); if (!(pi == 1. && !compute_o->OPTIMIZE_PI)) ASRVsiteLikelihood (s_tree); /* compute ASRV for non Huelsenbeck models */ if (pi>0.) { /* add a condition for doing covarion 23/03/2006*/ ret = CovarionSiteLikelihood (s_tree); if (ret == 1.) return 1.; // transition prob. and log likelihood cannot be computed } if (pi>0.) { if (pi==1. && !compute_o->OPTIMIZE_PI) { /* Huelsenbeck's model*/ for(i=0;ilgseq;i++) l+=log(s_tree->like_s[i])*s_tree->weight[i]; } else { /* Galtier or General models*/ for(i=0;ilgseq;i++) l+=log(pi*s_tree->like_s[i]+(1-pi)*s_tree->likeasrv_s[i])*s_tree->weight[i]; } } else { /* ASRV model only */ for(i=0;ilgseq;i++) l+=log(s_tree->likeasrv_s[i])*s_tree->weight[i]; /* change like_s to likeasrv_s. 23/03/2006 */ } if(l>0.) { printf("Warning: log likelihood greater than 0! parameters :\n"); for(i=0;ivar_param[i])); exit(0); } return -l; } void set_proba_zero(noeud nd){ int i, j, cl1, cl2, p; if(nd==NULL) return; for(cl1=0;cl1p1[cl1][cl2][i][j]=0.; nd->p2[cl1][cl2][i][j]=0.; for(p=0;p<5;p++){ nd->dp1[cl1][cl2][p][i][j]=0.; nd->d2p1[cl1][cl2][p][i][j]=0.; nd->dp2[cl1][cl2][p][i][j]=0.; nd->d2p2[cl1][cl2][p][i][j]=0.; } } } } } set_proba_zero(nd->v1); set_proba_zero(nd->v2); } double true_length(noeud nd){ double mylg, teta, lat, lcg, glat, glcg, coeff; noeud v3; int i; double alpha=2; // temporary setting if(nd==root) return -1.; coeff=0.; v3=nd->v3; mylg=nd->l3; teta=nd->ceq+nd->geq; lat=(1+alpha*teta)/2.; lcg=(1+alpha*(1-teta))/2.; for(i=0;ia*glat + v3->c*glcg + v3->g*glcg + v3->t*glat; } return mylg*coeff/gamma_nbcl; } double my_length(noeud nd, double alpha){ double truelg, coeff, teta, lat, lcg, glat, glcg; noeud v3; int i; if(nd==root) return -1.; coeff=0.; v3=nd->v3; truelg=nd->l3; teta=nd->ceq+nd->geq; lat=(1+alpha*teta)/2.; lcg=(1+alpha*(1-teta))/2.; for(i=0;iaeq*glat + nd->ceq*glcg + nd->geq*glcg + nd->teq*glat; /* note: equil base comp are used here instead of node base comp */ /* because this (approximate) calculation occurs at the beginning */ /* of the prog, i.e. when node base comp are unknown */ } return truelg*gamma_nbcl/coeff; } void mylg_to_truelg(noeud nd){ double tl; noeud v3; if(nd==NULL) return; tl=true_length(nd); if(tl>=0.){ nd->l3=tl; v3=nd->v3; if(nd==v3->v1) v3->l1=tl; else v3->l2=tl; } mylg_to_truelg(nd->v1); mylg_to_truelg(nd->v2); } void truelg_to_mylg(noeud nd, double alpha){ double ml; noeud v3; if(nd==NULL) return; ml=my_length(nd, alpha); if(ml>=0.){ nd->l3=ml; v3=nd->v3; if(nd==v3->v1) v3->l1=ml; else v3->l2=ml; } truelg_to_mylg(nd->v1, alpha); truelg_to_mylg(nd->v2, alpha); } void output_tree_comp(noeud nd, FILE* out, tree* s_tree){ static int numappel=0; int indbl; float m1=-1.; double truelg; if(nd==root){ numappel++; if(out) fprintf(out, ";;data set %d\n", numappel); else return; // set_node_deduced_comp(root); } /* trouver le numero de parametre de la branche et du GC du noeud */ indbl=s_tree->which_var_param[1]+nd->numerobr; /* which_var_param[4] in the following 2 lines of code cannot be computed: subscript out of range. param[4] is for GC content which is not defined in the protein model */ /* indgc=s_tree->which_var_param[4]+nd->numerobr; if(nd==root->v2) indgc=s_tree->which_var_param[4]+2*s_tree->nbseq-3; */ truelg=true_length(nd); if(nd->v1==NULL){ fprintf(out, "noeud %s (depth 0): bl=%f , aeq=%f , ceq=%f , geq=%f , teq=%f, a=%f, c=%f, g=%f, t=%f, truel=%f, dlkh/dbl=%f, d2lkh/dbl2=%f\n" /* dlkh/dgc=%f, d2lkh/dgc2=%f, gc=%f\n" */ , nd->nom, nd->l3, nd->aeq, nd->ceq, nd->geq, nd->teq, nd->a, nd->c, nd->g, nd->t, truelg, (nd==root || indbl<0)?m1:s_tree->deriv[indbl], (nd==root|| indbl<0)?m1:s_tree->deriv2[indbl] /*, (nd==root || indgc<0)?m1:s_tree->deriv[indgc], (nd==root || indgc<0)?m1:s_tree->deriv2[indgc], nd->c+nd->g */); return; } fprintf(out, "noeud %s (v1=%s , v2=%s) (depth %d) : bl=%f , aeq=%f , ceq=%f , geq=%f , teq=%f, a=%f, c=%f, g=%f, t=%f, truel=%f, dlkh/dbl=%f, d2lkh/dbl2=%f\n" /* , dlkh/dgc=%f, d2lkh/dgc2=%f, gc=%f\n" */, nd->nom, nd->v1->nom, nd->v2->nom, nd->depth, nd->l3, nd->aeq, nd->ceq, nd->geq, nd->teq, nd->a, nd->c, nd->g, nd->t, truelg, (nd==root || indbl<0)?m1:s_tree->deriv[indbl], (nd==root || indbl<0)?m1:s_tree->deriv2[indbl] /*, (nd==root || indgc<0)?m1:s_tree->deriv[indgc], (nd==root || indgc<0)?m1:s_tree->deriv2[indgc], nd->c+nd->g */ ); output_tree_comp(nd->v1, out, s_tree); output_tree_comp(nd->v2, out, s_tree); } double simplex(double* initpoint, double bound, double* finalvect){ int i, j, nbparam, nbiter, min; double ** point, *tempoint, *pointlike, ran, **contraintes, ret; nbparam=s_tree->nb_var_param; tempoint = (double *) check_alloc(nbparam, sizeof(double)); /* set contraintes */ contraintes=(double**)check_alloc(nbparam+1, sizeof(double*)); for(i=1;i<=nbparam;i++) contraintes[i]=(double*)check_alloc(2, sizeof(double)); for(i=1;i<=nbparam;i++){ contraintes[i][0]=min_value[s_tree->param_nature[i-1]]; contraintes[i][1]=max_value[s_tree->param_nature[i-1]]; } /* choose N+1 initial points around initial values given as argument */ /* WARNING : arrays start from 1 */ point=(double**)check_alloc(nbparam+2, sizeof(double*)); for(i=0;i<=nbparam+1;i++) point[i]=(double*)check_alloc(nbparam+1, sizeof(double)); pointlike=(double*)check_alloc(nbparam+2, sizeof(double)); for(i=2;i<=nbparam+1;i++){ for(j=1;j<=nbparam;j++){ ran=drand48(); /* 0contraintes[j][1]) point[i][j]=contraintes[j][1]; } } for(j=1;j<=nbparam;j++) point[1][j]=initpoint[j-1]; /* compute the likelihoods of each initial point */ for(i=1;i<=nbparam+1;i++){ for (j=0; jv1==NULL){ nd->alive=0; return; } set_alive_0(nd->v1); set_alive_0(nd->v2); } void set_tiplg_min(upgmanoeud nd, double tiplgmin){ if(nd->v1==NULL){ nd->l3-=tiplgmin; nd->v3->l1=nd->v3->l2=nd->l3; return; } set_tiplg_min(nd->v1, tiplgmin); set_tiplg_min(nd->v2, tiplgmin); } char* labelled_upgma(int nb, char** nom, double** dist, double *label){ int i, j, j1, j2, k1, k2, iter, nbn, nbm1, ndnum; double dmin, d1, d2, tiplgmin; char* ctree; upgmanoeud* stree; /* alloc nodes */ nbn=2*nb-1; stree=check_alloc(nbn, sizeof(upgmanoeud)); /* set tip nodes */ for(i=0;iv1=stree[i]->v2=NULL; stree[i]->l1=stree[i]->l2=-1; stree[i]->num=i; stree[i]->alive=1; stree[i]->depth=0.; stree[i]->nom=nom[i]; stree[i]->b3=label[i]; } /* main loop */ nbm1=nb-1; ndnum=nb; for(iter=0;iteralive==0) continue; for(j=i+1;jalive==0) continue; if(dist[i][j]num==j1){k1=i; break;} for(i=0;inum==j2){k2=i; break;} /* cluster */ stree[ndnum]=check_alloc(1, sizeof(struct upgmanoeud)); stree[ndnum]->v1=stree[k1]; stree[ndnum]->v2=stree[k2]; stree[k1]->v3=stree[ndnum]; stree[k2]->v3=stree[ndnum]; stree[ndnum]->l1=stree[k1]->l3=dmin/2.-stree[k1]->depth; stree[ndnum]->l2=stree[k2]->l3=dmin/2.-stree[k2]->depth; if(iter==0) tiplgmin=stree[ndnum]->l1; stree[ndnum]->depth=dmin/2.; stree[ndnum]->num=stree[k1]->num; stree[ndnum]->b3=(stree[k1]->b3+stree[k2]->b3)/2.; stree[k1]->num=-stree[k1]->num-1; stree[k2]->num=-stree[k2]->num-1; set_alive_0(stree[k2]); ndnum++; /* update distances */ for(i=0;ij1) dist[j1][i]=(d1+d2)/2.; } } /* struct to string */ ctree=check_alloc(nb*100, sizeof(char)); upgmaroot=stree[ndnum-1]; upgmaroot->v3=NULL; set_tiplg_min(upgmaroot, tiplgmin); upgma_stoc(stree, 1, nb, ctree, 0); return ctree; } /* output_covar_pattern */ /* Compute and print into file outcovfile : */ /* - the "variation index" of each site (contribution to the total likelihood */ /* of scenarii involving at least one change of rate class ) */ /* - the "covariation index" of each pair of site (likelihood conditional */ /* on simultaneous chanegs of rate class divided by unconstrained likelihood) */ void output_covar_pattern(tree* s_tree, double lcovar, double* bestvect, int optcov, FILE *outcov){ int i, j, k, ii, jj, lgseq, truelgseq, *pattnumber; int maxsitepername=5; double lnoch, lcov, *varind, **covarind; char **name, **longname, **finlongname, *prov, *covartree; // FILE* outcov; printf("computing variation and covariation indices...\n"); // outcov=fopen("covarfile", "w"); if(outcov==NULL) return; lgseq=s_tree->lgseq; truelgseq=s_tree->weightsum; varind=check_alloc(lgseq, sizeof(double)); covarind=check_alloc(lgseq, sizeof(double*)); for(k=0;klike_s[i]; for(j=i;jlike_s[i]*s_tree->like_s[j]; } /* variation indices */ covar=-1.; if(optcov) bestvect[s_tree->which_var_param[2]]=-1.; // let covarion rate = -1. // log_like_total(s_tree); // temporary disable it for(i=0;ilike_s[i]; varind[i]=1.-(lnoch/varind[i]); } covar=lcovar; if(optcov) bestvect[s_tree->which_var_param[2]]=lcovar; // change covar to the original value /* covariation indices */ for(i=0;isitetopatt[i]; k=pattnumber[j]; pattnumber[j]++; if(kmaxsitepername){ prov=name[j]; sscanf(prov, "set%d", &jj); jj--; sprintf(finlongname[jj], "_%d", i+1); while(*(finlongname[jj])) finlongname[jj]++; } } /* compute distances */ for(i=0;i=truelgseq) break; fprintf(outcov, "%d\t%.4e\n", i+1, varind[s_tree->sitetopatt[i]]); } if(i>=truelgseq) break; } if(i>=truelgseq) break; ii++; } // fclose(outcov); } /* my_eigen */ /* Diagonalize square matrix A (dim=n) */ /* job=0 -> compute eigen values only */ /* job=1 -> compute eigen values and eigen vectors */ /* eigen values returned in rr (real part) and ri (complex part) */ /* right eigen vectors returned as vr and vi */ /* A = v.diag(r).v-1 , where diag(r) is the diagonal */ /* matrix of eigenvalues, r is the matrix of right */ /* eigen vectors, and r-1 is the inverse of r */ /* rr, ri vr and vi must be passed allocated */ int my_eigen(int job, double** A, int n, double *rr, double *ri , double** vr, double** vi){ int i,j, ret; double* newA; double* newvr, *newvi, *w; newA=check_alloc(n*n, sizeof(double)); newvr=check_alloc(n*n, sizeof(double)); newvi=check_alloc(n*n, sizeof(double)); w=check_alloc(2*n, sizeof(double)); for(i=0;iv1) tl1=compute_treeLength(nd->v1)+nd->l1; else tl1=0.; if (nd->v2) tl2=compute_treeLength(nd->v2)+nd->l2; else tl2=0.; return (tl1+tl2); } void getTreelength(double *oriTreeLen, double *treeLen) { int i;double sum0; sum0 = 0.; for(i=0;i<2*s_tree->nbseq-3;i++) sum0 += s_tree->lgbr_init[i]; *oriTreeLen = sum0; if(compute_o->OPTIMIZE_LENGTH) *treeLen = compute_treeLength(root); else *treeLen = *oriTreeLen; } double compute(tree* s_tree, init_option init_o, compute_option compute_opt, converge_option converge_o, print_option print_o) { int i, j, ii, nbseq, lgseq, nbparam, firstcomp, nbreduction, best_gamma; int OPTIMIZE, nbiter, *bounded, do_simplex=0, do_simplex_now=0; double rem, bound, *toadd, totallength, newtotallength, gamma_l, best_gamma_l; double llike, maxllike, **invmatparam, *paramvect, *remvect, exact_llike; double prec, gamma_shp, *useless, lll, oriTreeLen, treeLen; double* bestvect, *finalvect; char type[20]; static FILE* out=NULL; int den; double l1, l_1, delta, d_gamma, d2_gamma, toadd_gamma, old_gamma_shp; char *ctree; // originally defined in eval_nh.c in Galtier's code char *resultFileName; ttotdeb=time(NULL); /* printmemory("deb compute");*/ if(print_o->PRINT3){ resultFileName = check_alloc(strlen(infileName)+12, sizeof(char)); strcpy(resultFileName, infileName); strcat(resultFileName,".result.dat"); if(!out) out=fopen(resultFileName, "w"); if(!out) { perror("Cannot create file: "); exit(1); } } if(print_o->PRINT1==0) noprint=1; else noprint=0; nbseq=s_tree->nbseq; lgseq=s_tree->lgseq; nbparam=s_tree->nb_var_param; if(!noprint) printf("%d species, %d sites, %d patterns\n", nbseq, (int)s_tree->weightsum, lgseq); if (print_o->PRINT3) { fprintf(out,"input sequences: %s\n", infileName); fprintf(out,"input tree: %s\n\n", intreeName); fprintf(out,"%d species, %d sites, %d patterns\n\n", nbseq, (int)s_tree->weightsum, lgseq); } prec=converge_o->PRECISION; compute_o=compute_opt; OPTIMIZE = compute_o->OPTIMIZE_LENGTH || compute_o->OPTIMIZE_ROOT || compute_o->OPTIMIZE_COV1 || compute_o->OPTIMIZE_COV2 || compute_o->OPTIMIZE_COV3 || compute_o->OPTIMIZE_GAMMA || compute_o->OPTIMIZE_PI; bound=pow(10., -1.*(prec)); /* INIT PARAMETERS */ /* root */ s_tree->froot=init_o->INIT_ROOT; /* Gamma distribution */ gamma_shp=init_o->INIT_GAMMA; gamma_nbcl=init_o->GAMMA_NCL; // original number of gamma categories in the input file gamma_nbcl1=gamma_nbcl*2; if(gamma_nbcl>1 && !(gamma_shp<=0. && compute_o->OPTIMIZE_GAMMA==0)){ if (gamma_shp<=0.) gamma_shp=gammamax; gamma_clmean=(double*)check_alloc(gamma_nbcl, sizeof(double)); useless=(double*)check_alloc(gamma_nbcl, sizeof(double)); DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); } else{ if(gamma_shp<=0.) gamma_shp=gammamax; gamma_nbcl=1; gamma_nbcl1=2; gamma_clmean=(double*)check_alloc(gamma_nbcl, sizeof(double)); gamma_clmean[0]=1.0; compute_o->OPTIMIZE_GAMMA=0; useless=(double*)check_alloc(gamma_nbcl, sizeof(double)); } /* lengths */ if(strcmp(init_o->INIT_LENGTH, "REDO")==0) init_bl(s_tree); /* covarions */ covar1=init_o->INIT_COV1; covar2=init_o->INIT_COV2; covar3=init_o->INIT_COV3; pi=init_o->INIT_PI; if(covar1OPTIMIZE_COV1==1) covar1=covmin; if(covar1>covmax && compute_o->OPTIMIZE_COV1==1) covar1=covmax; if(covar2OPTIMIZE_COV2==1) covar2=covmin; if(covar2>covmax && compute_o->OPTIMIZE_COV2==1) covar2=covmax; if(covar3OPTIMIZE_COV3==1) covar3=covmin; if(covar3>covmax && compute_o->OPTIMIZE_COV3==1) covar3=covmax; p1prov=check_alloc(gamma_nbcl, sizeof(proba*)); p2prov=check_alloc(gamma_nbcl, sizeof(proba*)); for(i=0;iNBRANDOM>0){ if(init_o->NOISE<0){ totallength=newtotallength=0.; for(i=0;i<2*nbseq-3;i++) totallength+=s_tree->listbr[i]->length; s_tree->froot=drand48(); for(i=0;i<2*nbseq-3;i++){ s_tree->listbr[i]->length=drand48(); newtotallength+=s_tree->listbr[i]->length; } for(i=0;i<2*nbseq-3;i++){ s_tree->listbr[i]->length*=(totallength/newtotallength); } } else{ for(i=0;ivar_param[i])*= 1+ (drand48()-0.5)*2.*init_o->NOISE/100.; } } } toadd=(double*)check_alloc(nbparam, sizeof(double)); invmatparam=(double**)check_alloc(nbparam, sizeof(double*)); for(i=0;ivar_param[i]); while(TRUE){ /* start of Newton-Raphson loop */ nbiter++; printf ("\nnbiter=%d\n", nbiter); firstcomp=1; nbreduction=0; for(i=0;ivar_param[i]); /* compute likelihood, compare to previous, */ /* redo (up to 10 times) until it increases */ while(firstcomp || llikevar_param[i])=remvect[i]+toadd[i]; for(i=0;ivar_param[i]); llike=-minus_log_likelihood(paramvect); firstcomp=0; // if(nbiterPRINT2) printf("likelihood decrease (%f->%f): moving backward\n", rem, llike); for(i=0;ivar_param[i])=remvect[i]; /* if likelihood decrease 10 times no more modification of the parameters. restore the parameters to the values before the first time likelihood decreases. */ reset_tree(s_tree); do_simplex_now=1; break; } } } if(do_simplex_now){ do_simplex=1; break; } if(nbiter==1) maxllike=llike; else{ if(maxllike0.) { printf("log negatif\n"); if(print_o->PRINT3){ fprintf(out, "Negative log\n"); output_tree_comp(root, NULL, NULL); } return 1.; } if(print_o->PRINT2) printf("pi: %.8f\tgamma classes: %d\tgamma_shp: %.8f\t\tlog likelihood: %.12f\n", pi, gamma_nbcl,gamma_shp,llike); if(print_o->PRINT2) printf("cov1: %.12f\tcov2: %.12f\tcov3: %.12f\n", covar1,covar2,covar3); if(!OPTIMIZE) break; /* leave if converged */ /* nbiter>10 allows always do more than 10 iterations*/ if(nbiter!=1 && llike-rem-bound && nbiter>10){ for(i=0;ilkh=llike; /* first steps: try several gamma distributions */ old_gamma_shp=gamma_shp; if(compute_o->OPTIMIZE_GAMMA && nbiterbest_gamma_l){ best_gamma=i; best_gamma_l=gamma_l; } } if(best_gamma!=-1){ llike=best_gamma_l; rem=llike; s_tree->lkh=llike; gamma_shp=gamma_range[best_gamma]; } else gamma_shp=old_gamma_shp; DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); minus_log_likelihood(paramvect); } /* analytical derivatives for ASRV and numerical derivative for covarion (1st + 2nd) */ derivatives(s_tree, paramvect); /* not first steps: */ /* numerical derivative with respect to shape param of gamma distribution, */ /* and modification of gamma shape parameter. */ if(compute_o->OPTIMIZE_GAMMA && nbiter>GAMMA_RANGE_ITER) { delta=0.00001; lll=s_tree->lkh; gamma_shp+=delta; DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); l1=-minus_log_likelihood(paramvect); gamma_shp-=2*delta; DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); l_1=-minus_log_likelihood(paramvect); gamma_shp+=delta; d_gamma=(l1-l_1)/(2*delta); d2_gamma=(l1+l_1-2*lll)/(delta*delta); /* we changed Galtier's equation for 2nd deriv.*/ den=1; old_gamma_shp=gamma_shp; if(d2_gamma<0.){ toadd_gamma=-d_gamma/d2_gamma; while(den<1000){ gamma_shp=old_gamma_shp+toadd_gamma/den; if(gamma_shp1000 || d2_gamma>=0.){ printf("den=%d\td2_gamma=%.10f\n",den,d2_gamma); gamma_shp=old_gamma_shp; DiscreteGamma (useless, gamma_clmean, gamma_shp, gamma_shp, gamma_nbcl, 0); } else llike=lll; } /* modify parameters */ for(ii=0;iinb_var_param;ii++){ if (s_tree->deriv2[ii]==0.0) { if (s_tree->deriv[ii]==0.0) toadd[ii] = 0.0; else if (s_tree->deriv[ii]>0.0) switch(s_tree->param_nature[ii]){ case LENGTH: toadd[ii] = lmin-remvect[ii];break; case COVA:case COVB:case COVC: toadd[ii] = covmin-remvect[ii];break; case GaltierPI: toadd[ii] = pimin-remvect[ii]; break; } else switch(s_tree->param_nature[ii]){ case LENGTH: toadd[ii] = lmax-remvect[ii];break; case COVA:case COVB:case COVC: toadd[ii] = covmax-remvect[ii];break; case GaltierPI: toadd[ii] = pimax-remvect[ii]; break; } } else toadd[ii]=-s_tree->deriv[ii]/s_tree->deriv2[ii]; } for(i=0;ideriv2[i]>0.){ toadd[i]=0.; if(!print_o->PRINT2) continue; switch(s_tree->param_nature[i]){ case ROOT: sprintf(type, "root"); break; case LENGTH: sprintf(type, "length"); break; case COVA: sprintf(type, "covar1"); break; case COVB: sprintf(type, "covar2"); break; case COVC: sprintf(type, "covar3"); break; case GaltierPI: sprintf(type, "pi"); break; } printf("param %d (%s) diverging\n", i+1, type); } } } /* end of Newton-Raphson loop */ covarianceMatrix(s_tree, paramvect, gamma_shp, out); exact_llike=maxllike; if(covar>0. && print_o->PRINT2) printf("Exact log likelihood: %f\n", exact_llike); fprintf (out, "log likelihood for exact covarion before simplex: %.6f\n", exact_llike); if(covar>0.1) do_simplex=1; // if(covar<=0.1) exact_cov=0; do_simplex = 1; // simplex does not work properly /* SIMPLEX */ if((converge_o->SIMPLEX || covar>0.1) && do_simplex){ if(print_o->PRINT2) printf("Simplex...\nWarning: Simplex may take long time...\n"); llike=simplex(bestvect, 0.01, finalvect); if(llike>maxllike){ for(i=0;ivar_param[i])=finalvect[i]; reset_tree(s_tree); maxllike=llike; } else{ for(i=0;ivar_param[i])=bestvect[i]; reset_tree(s_tree); llike=maxllike; } if(print_o->PRINT3) fprintf (out, "log likelihood for exact covarion after simplex: %.6f\n\n", llike); } getTreelength(&oriTreeLen, &treeLen); /* OUTPUT */ s_tree->gamma_shp=gamma_shp; s_tree->covar=covar; s_tree->pi=pi; if (print_o->PRINT1){ printf("\nln(L) = %.6f\tnbiter = %d\n\n", maxllike, nbiter); if(print_o->PRINT3) fprintf (out, "maximum log likelihood is: %.10f\n\n", maxllike); fprintf(out,"Input number of gamma categories: %d\n",init_o->GAMMA_NCL); if(compute_o->OPTIMIZE_GAMMA) fprintf(out,"Initial gamma shape parameter: %.4f\tEstimated Gamma shape parameter: %.4f\n", init_o->INIT_GAMMA, gamma_shp); else if(init_o->INIT_GAMMA>0. && gamma_nbcl>1) fprintf(out,"Fixed Gamma shape parameter: %.4f\n", gamma_shp); else fprintf(out,"Constant rates accross sites\n"); if(compute_o->OPTIMIZE_COV1) fprintf(out,"Initial covarion1: %.4f\tEstimated covarion1 parameter: %.4f\n", init_o->INIT_COV1, covar1); else if(init_o->INIT_COV1>0.) fprintf(out,"Fixed covarion1 parameter: %.4f\n", covar1); else fprintf(out,"No covarion1\n"); if(compute_o->OPTIMIZE_COV2) fprintf(out,"Initial covarion2 parameter: %.4f\tEstimated covarion2 parameter: %.4f\n", init_o->INIT_COV2, covar2); else if(init_o->INIT_COV2>0.) fprintf(out,"Fixed covarion2 parameter: %.4f\n", covar2); else fprintf(out,"No covarion2\n"); if(compute_o->OPTIMIZE_COV3) fprintf(out,"Initial covarion3 parameter: %.4f\tEstimated covarion3 parameter: %.4f\n", init_o->INIT_COV3, covar3); else if(init_o->INIT_COV3>0.) fprintf(out,"Fixed covarion3 parameter: %.4f\n", covar3); else fprintf(out,"No covarion3\n"); if(compute_o->OPTIMIZE_PI) fprintf(out,"Initial proportion of covarions: %.4f\tEstimated proportion of covarions: %.4f\n\n", init_o->INIT_PI, pi); else if(init_o->INIT_PI>=0.) fprintf(out,"Fixed proportion of covarions: %.4f\n", pi); if(compute_o->OPTIMIZE_LENGTH) fprintf(out,"Original tree length: %.4f\tEstimated tree length: %.4f\n\n", oriTreeLen, treeLen); else fprintf(out,"Fixed tree length: %.4f\n\n", oriTreeLen); } /* if(print_o->PRINT1 && print_o->PRINT3){ mindepth(NULL, root); output_tree_comp(root, out, s_tree); printf("Detailed results are written into file : resfile\n"); } */ /* write log site likelihoods so they can be used by consel */ for(i=0;iweightsum;i++) { j = s_tree->sitetopatt[i]; fprintf(out, "%.6f\t", log(pi*s_tree->like_s[j]+(1-pi)*s_tree->likeasrv_s[j])); } fprintf (out, "\n"); /* if((compute_o->OPTIMIZE_COV || init_o->INIT_COV>0.) && print_o->COV_SITE_PATTERN){ if(covar>covmin){ output_covar_pattern(s_tree, covar, bestvect, compute_o->OPTIMIZE_COV, out); printf("Site-specific covariation patterns are written into file : result.dat\n"); } } */ /* ESTIMATE ACTUAL BASE COMPOSITIONS AND BRANCH LENGTHS */ /* disable the following because we cannot compute base comp. and "true" branch length for the protein alignment - we don't know GC content and ti/tv for(k=0;knode, nbseq, 1, ctree, 0); if(print_o->PRINT3) fprintf(out, "%s\n", ctree); /* FREE */ free(resultFileName); free(toadd); free(bestvect); free(paramvect); free(finalvect); free(remvect); for(i=0;inb_var_param;i++) free(invmatparam[i]); free(invmatparam); free(bounded); free(gamma_clmean); free(useless); for(i=0;i