/* tools.c */ #include "paml.h" /************************ sequences *************************/ char BASEs[]="TCAGUYRMKSWHBVDN?-"; char nBASEs[]={1,1,1,1, 1,2,2,2,2,2,2, 3,3,3,3, 4,4,4}; char *EquateNUC[]={"T","C","A","G", "T", "TC","AG","CA","TG","CG","TA", "TCA","TCG","CAG","TAG", "TCAG","TCAG","TCAG"}; char AAs[] = "ARNDCQEGHILKMFPSTWYV*-?x"; char AA3Str[]= {"AlaArgAsnAspCysGlnGluGlyHisIleLeuLysMetPheProSerThrTrpTyrVal***"}; char BINs[] ="TC"; int GeneticCode[][64] = {{13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,-1,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 0:universal */ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15,-1,-1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 1:vertebrate mt.*/ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17, 16,16,16,16,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 2:yeast mt. */ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 3:mold mt. */ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15,15,15, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 4:invertebrate mt. */ {13,13,10,10,15,15,15,15,18,18, 5, 5, 4, 4,-1,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 5:ciliate nuclear*/ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9, 9,12,16,16,16,16, 2, 2, 2,11,15,15,15,15, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 6:echinoderm mt.*/ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4, 4,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 7:euplotid mt. */ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,-1,17, 10,10,10,15,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 8:alternative yeast nu.*/ {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15, 7, 7, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 9:ascidian mt. */ {13,13,10,10,15,15,15,15,18,18,-1, 5, 4, 4,-1,17, 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1, 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1, 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 10:blepharisma nu.*/ { 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9,10,10,10,10,11,11,11,11,12,12,12,12, 13,13,13,13,14,14,14,14,15,15,15,15,16,16,16,16} /* 11:Ziheng's regular code */ }; /* GeneticCode[icode][#codon] */ int noisy=0, Iround=0, NFunCall=0, NEigenQ, NPMatUVRoot; double SIZEp=0; int blankline (char *str) { char *p=str; while (*p) if (isalnum(*p++)) return(0); return(1); } int PopEmptyLines (FILE* fseq, int lline, char line[]) { /* pop out empty lines in the sequence data file. returns -1 if EOF. */ char *eqdel=".-?", *p; int i; for (i=0; ;i++) { p=fgets (line, lline, fseq); if (p==NULL) return(-1); while (*p) if (*p==eqdel[0] || *p==eqdel[1] || *p==eqdel[2] || isalpha(*p)) return(0); else p++; } } int picksite (char *z, int l, int begin, int gap, char *result) { /* pick every gap-th site, e.g., the third codon position for example. */ int il=begin; for (il=0, z+=begin; il=9) printf ("\nwarning: strange character '%c' ", b); return (-1); } int dnamaker (char z[], int ls, double pi[]) { /* sequences z[] are coded 0,1,2,3 */ int i, j; double p[4], r, small=1e-5; xtoy(pi, p, 4); for (i=1; i<4; i++) p[i]+=p[i-1]; if(fabs(p[3]-1)>small) error2("sum pi != 1.."); for(i=0; iualpha /* && fabs(u2)>ualpha */ ) { fprintf (fout,"\n"); FOR (j, lfrag) fprintf (fout,"%1c", BASEs[ib[j]]); fprintf (fout,"%6d %8.1f%7.2f %8.1f%7.2f ",ffrag[i],ne1,u1,ne2,u2); if (u1<-ualpha && u2<-ualpha) fprintf (fout, " %c", '-'); else if (u1>ualpha && u2>ualpha) fprintf (fout, " %c", '+'); else if (u1*u2<0 && fabs(u1) > ualpha && fabs(u2) > ualpha) fprintf (fout, " %c", '?'); else fprintf (fout, " %c", ' '); } } return (0); } int zztox ( int n31, int l, char *z1, char *z2, double *x ) { /* x[n31][4][4] */ double t = 1./(double) (l / n31); int i, ib[2]; int il; zero (x, n31*16); for (i=0; i1) it=-1; t += x[i*4+j]; } if (fabs(t-1) > 1e-4) it =-1; return(it); } int difcodonNG (char codon1[], char codon2[], double *SynSite,double *AsynSite, double *SynDif, double *AsynDif, int transfed, int icode) { /* # of synonymous and non-synonymous sites and differences. Nei, M. and T. Gojobori (1986) returns the number of differences between two codons. The two codons (codon1 & codon2) do not contain ambiguity characters. dmark[i] (=0,1,2) is the i_th different codon position, with i=0,1,ndiff step[j] (=0,1,2) is the codon position to be changed at step j (j=0,1,ndiff) b[i][j] (=0,1,2,3) is the nucleotide at position j (0,1,2) in codon i (0,1) I made some arbitrary decisions when the two codons have ambiguity characters 20 September 2002. */ int i,j,k, i1,i2, iy[2]={0}, iaa[2],ic[2]; int ndiff,npath,nstop,sdpath,ndpath,dmark[3],step[3],b[2][3],bt1[3],bt2[3]; int by[3] = {16, 4, 1}; char str[4]=""; for (i=0,*SynSite=0,nstop=0; i<2; i++) { for (j=0; j<3; j++) { if (transfed) b[i][j]=(i?codon1[j]:codon2[j]); else b[i][j]=(int)CodeChara((char)(i?codon1[j]:codon2[j]),0); iy[i]+=by[j]*b[i][j]; if(b[i][j]<0||b[i][j]>3) { if(noisy>=9) printf("\nwarning ambiguity in difcodonNG: %s %s", codon1,codon2); *SynSite = 0.5; *AsynSite = 2.5; *SynDif = (codon1[2]!=codon2[2])/2; *AsynDif = *SynDif + (codon1[0]!=codon2[0])+(codon1[1]!=codon2[1]); return((int)(*SynDif + *AsynDif)); } } iaa[i]=GeneticCode[icode][iy[i]]; if(iaa[i]==-1) { printf("\nNG86: stop codon %s.\n",getcodon(str,iy[i])); exit(-1); } for(j=0; j<3; j++) FOR (k,4) { if (k==b[i][j]) continue; i1=GeneticCode[icode][iy[i] + (k-b[i][j])*by[j]]; if (i1==-1) nstop++; else if (i1==iaa[i]) (*SynSite)++; } } *SynSite*=3/18.; /* 2 codons, 2*9 possibilities. */ *AsynSite=3*(1-nstop/18.) - *SynSite; #if 0 /* MEGA 1.1 */ *AsynSite = 3 - *SynSite; #endif ndiff=0; *SynDif=*AsynDif=0; FOR (k,3) dmark[k]=-1; FOR (k,3) if (b[0][k]-b[1][k]) dmark[ndiff++]=k; if (ndiff==0) return(0); npath=1; nstop=0; if(ndiff>1) npath=(ndiff==2)?2:6 ; if (ndiff==1) { if (iaa[0]==iaa[1]) (*SynDif)++; else (*AsynDif)++; } else { /* ndiff=2 or 3 */ FOR (k, npath) { FOR (i1,3) step[i1]=-1; if (ndiff==2) { step[0]=dmark[k]; step[1]=dmark[1-k]; } else { step[0]=k/2; step[1]=k%2; if (step[0]<=step[1]) step[1]++; step[2]=3-step[0]-step[1]; } FOR (i1,3) bt1[i1]=bt2[i1]=b[0][i1]; sdpath=ndpath=0; /* mutations for each path */ FOR (i1, ndiff) { /* mutation steps for each path */ bt2[step[i1]] = b[1][step[i1]]; for (i2=0,ic[0]=ic[1]=0; i2<3; i2++) { ic[0]+=bt1[i2]*by[i2]; ic[1]+=bt2[i2]*by[i2]; } FOR (i2,2) iaa[i2]=GeneticCode[icode][ic[i2]]; if (iaa[1]==-1) { nstop++; sdpath=ndpath=0; break; } if (iaa[0]==iaa[1]) sdpath++; else ndpath++; FOR (i2,3) bt1[i2]=bt2[i2]; } *SynDif+=(double)sdpath; *AsynDif+=(double)ndpath; } } if (npath==nstop) { puts ("NG86: All paths are through stop codons.."); if (ndiff==2) { *SynDif=0; *AsynDif=2; } else { *SynDif=1; *AsynDif=2; } } else { *SynDif/=(double)(npath-nstop); *AsynDif/=(double)(npath-nstop); } return (ndiff); } int difcodonLWL85 (char codon1[], char codon2[], double sites[3], double sdiff[3], double vdiff[3], int transfed, int icode) { /* This partitions codon sites according to degeneracy, that is sites[3] has L0, L2, L4, averaged between the two codons. It also compares the two codons and add the differences to the transition and transversion differences for use in LWL85 and similar methods. The two codons (codon1 & codon2) should not contain ambiguity characters. */ int b[2][3], by[3] = {16, 4, 1}, i,j, ifold[2], c[2], ct, aa[2], ibase,nsame; char str[4]=""; for(i=0; i<3; i++) sites[i]=sdiff[i]=vdiff[i]=0; /* check the two codons and code them */ for (i=0; i<2; i++) { for (j=0,c[i]=0; j<3; j++) { if (transfed) b[i][j]=(i?codon1[j]:codon2[j]); else b[i][j]=(int)CodeChara((char)(i?codon1[j]:codon2[j]),0); c[i]+=b[i][j]*by[j]; if(b[i][j]<0||b[i][j]>3) { if(noisy>=9) printf("\nwarning ambiguity in difcodonLWL85: %s %s", codon1,codon2); return(0); } } aa[i]=GeneticCode[icode][c[i]]; if(aa[i]==-1) { printf("\nLWL85: stop codon %s.\n",getcodon(str,c[i])); exit(-1); } } for (j=0; j<3; j++) { /* three codon positions */ for (i=0; i<2; i++) { /* two codons */ for(ibase=0,nsame=0; ibase<4; ibase++) { ct=c[i]+(ibase-b[i][j])*by[j]; if(ibase!=b[i][j] && aa[i]==GeneticCode[icode][ct]) nsame++; } if(nsame==0) ifold[i]=0; else if (nsame==1 || nsame==2) ifold[i]=1; else ifold[i]=2; sites[ifold[i]]+=.5; } if(b[0][j]==b[1][j]) continue; if(b[0][j]+b[1][j]==1 || b[0][j]+b[1][j]==5) { /* transition */ sdiff[ifold[0]]+=.5; sdiff[ifold[1]]+=.5; } else { /* transversion */ vdiff[ifold[0]]+=.5; vdiff[ifold[1]]+=.5; } } return (0); } int testTransP (double P[], int n) { int i,j, status=0; double sum, small=1e-10; for (i=0; ismall && status==0) { printf ("\nrow sum (#%2d) = 1 = %10.6f", i+1, sum); status=-1; } } return (status); } int PMatUVRoot (double P[], double t, int n, double U[], double V[], double Root[]) { /* P(t) = U * exp{Root*t} * V */ int i,j,k; double expt, uexpt, *pP; double smallp = 0; if (t < 0.0) { printf("PMatUVRoot is called with negative branchlength! %f\n", t); } NPMatUVRoot++; if (t<-0.1) printf ("\nt = %.5f in PMatUVRoot", t); if (t<1e-100) { identity (P, n); return(0); } for (k=0,zero(P,n*n); k=5) if (testTransP(P,n)) { printf("\nP(%.6f) err in PMatUVRoot.\n", t); exit(-1); } #endif return (0); } int PMatQRev(double Q[], double pi[], double t, int n, double space[]) { /* This calculates P(t) = exp(Q*t), where Q is the rate matrix for a time-reversible Markov process. Q[] or P[] has the rate matrix as input, and P(t) in return. space[n*n*2+n*2] */ double *U=space, *V=U+n*n, *Root=V+n*n, *spacesqrtpi=Root+n; eigenQREV(Q, pi, n, Root, U, V, spacesqrtpi); PMatUVRoot(Q, t, n, U, V, Root); return(0); } void pijJC69 (double pij[2], double t) { if (t<-0.0001) printf ("\nt = %.5f in pijJC69", t); if (t<1e-100) { pij[0]=1; pij[1]=0; } else { pij[0] = (1.+3*exp(-4*t/3.))/4; pij[1] = (1-pij[0])/3; } } int PMatK80 (double P[], double t, double kappa) { /* PMat for JC69 and K80 */ int i,j; double e1, e2; if (t<-0.01) printf ("\nt = %.5f in PMatK80", t); if (t<1e-100) { identity (P, 4); return(0); } e1=exp(-4*t/(kappa+2)); if (fabs(kappa-1)<1e-5) { FOR (i,4) FOR (j,4) if (i==j) P[i*4+j]=(1+3*e1)/4; else P[i*4+j]=(1-e1)/4; } else { e2=exp(-2*t*(kappa+1)/(kappa+2)); FOR (i,4) P[i*4+i]=(1+e1+2*e2)/4; P[0*4+1]=P[1*4+0]=P[2*4+3]=P[3*4+2]=(1+e1-2*e2)/4; P[0*4+2]=P[0*4+3]=P[2*4+0]=P[3*4+0]= P[1*4+2]=P[1*4+3]=P[2*4+1]=P[3*4+1]=(1-e1)/4; } return (0); } int PMatT92 (double P[], double t, double kappa, double pGC) { /* PMat for Tamura'92 t is branch lnegth, number of changes per site. */ double e1, e2; t/=(pGC*(1-pGC)*kappa + .5); if (t<-0.1) printf ("\nt = %.5f in PMatT92", t); if (t<1e-100) { identity (P, 4); return(0); } e1=exp(-t); e2=exp(-(kappa+1)*t/2); P[0*4+0]=P[2*4+2] = (1-pGC)/2*(1+e1)+pGC*e2; P[1*4+1]=P[3*4+3] = pGC/2*(1+e1)+(1-pGC)*e2; P[1*4+0]=P[3*4+2] = (1-pGC)/2*(1+e1)-(1-pGC)*e2; P[0*4+1]=P[2*4+3] = pGC/2*(1+e1)-pGC*e2; P[0*4+2]=P[2*4+0]=P[3*4+0]=P[1*4+2] = (1-pGC)/2*(1-e1); P[1*4+3]=P[3*4+1]=P[0*4+3]=P[2*4+1] = pGC/2*(1-e1); return (0); } int PMatTN93 (double P[], double a1t, double a2t, double bt, double pi[]) { double T=pi[0],C=pi[1],A=pi[2],G=pi[3], Y=T+C, R=A+G; double e1, e2, e3, small=-1e-3; if(noisy && (a1t1e-5) e3=exp(-(Y*a1t+R*bt)); P[0*4+0] = T+R*T/Y*e1+C/Y*e3; P[0*4+1] = C+R*C/Y*e1-C/Y*e3; P[0*4+2] = A*(1-e1); P[0*4+3] = G*(1-e1); P[1*4+0] = T+R*T/Y*e1-T/Y*e3; P[1*4+1] = C+R*C/Y*e1+T/Y*e3; P[1*4+2] = A*(1-e1); P[1*4+3] = G*(1-e1); P[2*4+0] = T*(1-e1); P[2*4+1] = C*(1-e1); P[2*4+2] = A+Y*A/R*e1+G/R*e2; P[2*4+3] = G+Y*G/R*e1-G/R*e2; P[3*4+0] = T*(1-e1); P[3*4+1] = C*(1-e1); P[3*4+2] = A+Y*A/R*e1-A/R*e2; P[3*4+3] = G+Y*G/R*e1+A/R*e2; return(0); } int EvolveHKY85 (char source[], char target[], int ls, double t, double rates[], double pi[4], double kappa, int isHKY85) { /* isHKY85=1 if HKY85, =0 if F84 Use NULL for rates if rates are identical among sites. */ int i,j,h,n=4; double TransP[16],a1t,a2t,bt,r, Y=pi[0]+pi[1],R=pi[2]+pi[3]; if (isHKY85) a1t=a2t=kappa; else { a1t=1+kappa/Y; a2t=1+kappa/R; } bt=t/(2*(pi[0]*pi[1]*a1t+pi[2]*pi[3]*a2t)+2*Y*R); a1t*=bt; a2t*=bt; FOR (h, ls) { if (h==0 || (rates && rates[h]!=rates[h-1])) { r=(rates?rates[h]:1); PMatTN93 (TransP, a1t*r, a2t*r, bt*r, pi); for (i=0;i1e-5) error2("TransP err"); } } for (j=0,i=source[h],r=rndu();j1). Rates are converted into the c.d.f. if cdf=1, which is useful for simulation under JC69-like models. space[ncatG*5] */ int h, ir, j, K=ncatG, *Lalias=(int*)(space+3*K), *counts=(int*)(space+4*K); double *rK=space, *freqK=space+K, *Falias=space+2*K; if (alpha==0) { if(rates) FOR(h,ls) rates[h]=1; } else { if (K>1) { DiscreteGamma (freqK, rK, alpha, alpha, K, DGammaMean); MultiNomialAliasSetTable(K, freqK, Falias, Lalias, space+5*K); MultiNomialAlias(ls, K, Falias, Lalias, counts); for (ir=0,h=0; ir63) { printf("\ncodon %d\n", icodon); error2("getcodon."); } codon[0]=BASEs[icodon/16]; codon[1]=BASEs[(icodon%16)/4]; codon[2]=BASEs[icodon%4]; codon[3]=0; return (codon); } char *getAAstr(char *AAstr, int iaa) { /* iaa (0,20) with 20 meaning termination */ if (iaa<0 || iaa>20) error2("getAAstr: iaa err. \n"); strncpy (AAstr, AA3Str+iaa*3, 3); return (AAstr); } int NucListall(char b, int *nb, int ib[4]) { /* Resolve an ambiguity nucleotide b into all possibilities. nb is number of bases and ib (0,1,2,3) list all of them. Data are complete if (nb==1). */ int j,k; k=strchr(BASEs,b)-BASEs; if(k<0) { printf("NucListall: strange character %c\n",b); return(-1);} if(k<4) { *nb=1; ib[0]=k; } else { *nb=nBASEs[k]; FOR(j,*nb) ib[j]=strchr(BASEs,EquateNUC[k][j])-BASEs; } return(0); } int Codon2AA(char codon[3], char aa[3], int icode, int *iaa) { /* translate a triplet codon[] into amino acid (aa[] and iaa), using genetic code icode. This deals with ambiguity nucleotides. *iaa=(0,...,19), 20 for stop or missing data. Distinquish between stop codon and missing data? naa=0: only stop codons; 1: one AA; 2: more than 1 AA. Returns 0: if one amino acid 1: if multiple amino acids (ambiguity data) -1: if stop codon */ int nb[3],ib[3][4], ic, i, i0,i1,i2, iaa0=-1,naa=0; for(i=0;i<3;i++) NucListall(codon[i], &nb[i], ib[i]); FOR(i0,nb[0]) FOR(i1,nb[1]) FOR(i2,nb[2]) { ic=ib[0][i0]*16+ib[1][i1]*4+ib[2][i2]; *iaa=GeneticCode[icode][ic]; if(*iaa==-1) continue; if(naa==0) { iaa0=*iaa; naa++; } else if (*iaa!=iaa0) naa=2; } if(naa==0) { printf("stop codon %c%c%c\n",codon[0],codon[1],codon[2]); *iaa=20; } else if(naa==2) *iaa=20; else *iaa=iaa0; strncpy(aa, AA3Str+*iaa*3, 3); return(naa==1?0: (naa==0?-1:1)); } int DNA2protein(char dna[], char protein[], int lc, int icode) { /* translate a DNA into a protein, using genetic code icode, with lc codons. dna[] and protein[] can be the same string. */ int h, iaa, k; char aa3[4]; for(h=0; h0) fprintf(fout, " "); else { fprintf(fout, "%s %c", aa3,(iaa<20?AAs[iaa]:'*')); strcpy(ss3[k], aa3); } fprintf(fout, " %s", codon); if (fcodon) fprintf(fout, "%*.*f", wc,wd, fcodon[it] ); if (k<3) fprintf(fout, " %c ", word[0]); } FPN (fout); } fputs (noodle, fout); } return(0); } int printcums (FILE *fout, int ns, double fcodons[], int icode) { int neach0=6, neach=neach0, wc=3,wd=0; /* wc: for codon, wd: decimal */ int iaa,it, i,j,k, i1, ngroup, igroup; char *word="|-", aa3[4]=" ",codon[4]=" ", ss3[4][4], *noodle; ngroup=(ns-1)/neach+1; for(igroup=0; igroup0) fprintf(fout, " "); else { fprintf(fout, "%s", aa3); strcpy(ss3[k], aa3); } fprintf(fout, " %s", codon); FOR (i1,neach) fprintf(fout, " %*.*f", wc-1, wd, fcodons[(igroup*neach0+i1)*64+it] ); if (k<3) fprintf(fout, " %c ", word[0]); } FPN (fout); } fputs (noodle, fout); } } return(0); } int QtoPi (double Q[], double pi[], int n, double space[]) { /* from rate matrix Q[] to pi, the stationary frequencies: Q' * pi = 0 pi * 1 = 1 space[] is of size n*(n+1). */ int i,j; double *T = space; /* T[n*(n+1)] */ for(i=0;i255) error2("line >255 in strc"); FOR (i,n) s[i]=(char)c; s[n]=0; return (s); } int putdouble(FILE*fout, double a) { double aa=fabs(a); return fprintf(fout, (aa<1e-5||aa>1e6 ? " %11.4e" : " %11.6f"), a); } void strcase (char *str, int direction) { /* direction = 0: to lower; 1: to upper */ char *p=str; if(direction) while(*p) { *p=(char)toupper(*p); p++; } else while(*p) { *p=(char)tolower(*p); p++; } } FILE *gfopen(char *filename, char *mode) { FILE *fp; if(filename==NULL || filename[0]==0) error2("file name empty."); fp=(FILE*)fopen(filename, mode); if(fp==NULL) { printf("\nerror when opening file %s\n", filename); if(!strchr(mode,'r')) exit(-1); printf("tell me the full path-name of the file? "); scanf("%s", filename); if((fp=(FILE*)fopen(filename, mode))!=NULL) return(fp); puts("Can't find the file. I give up."); exit(-1); } return(fp); } int appendfile(FILE*fout, char*filename) { FILE *fin=fopen(filename,"r"); int ch; if(fin) { while((ch=fgetc(fin))!=EOF) fputc(ch,fout); fclose(fin); } return(0); } void error2 (char * message) { printf("\nError: %s.\n", message); exit(-1); } int zero (double x[], int n) { int i; FOR (i,n) x[i]=0; return (0);} double sum (double x[], int n) { int i; double t=0; for(i=0; i x[ibest[0]]) { nbest=1; ibest[0]=i; } else if(x[i] == x[ibest[0]]) ibest[nbest++]=i; } for(i=0; it) { t=x[j]; it=j; } } else { for ( ; jx else x->f. The iterative variable x[] and frequency f[0,1,n-2] are related as: freq[k] = exp(x[k])/(1+SUM(exp(x[k]))), k=0,1,...,n-2, x[] and freq[] may be the same vector. The last element (f[n-1] or x[n-1]=1) is updated only if(LastItem). */ int i; double tot; if (fromf) { /* f => x */ if((tot=1-sum(f,n-1))<1e-80) error2("f[n-1]==1, not dealt with."); tot=1/tot; FOR(i,n-1) x[i]=log(f[i]*tot); if(LastItem) x[n-1]=0; } else { /* x => f */ for(i=0,tot=1; ia && x=b) x = b - (x - b) - small; } if(rounds>maxround) printf("\nx=%.6g (%.6g, %.6g), reflected %d rounds.\n", x0,a,b, maxround); return(x); } void randorder(int order[], int n, int space[]) { /* This orders 0,1,2,...,n-1 at random space[n] */ int i,k, *item=space; for(i=0; i=0 && sumsquare<=1) break; } return (u1*sqrt(-2*log(sumsquare)/sumsquare)); } int rndpoisson (double m) { /* m is the rate parameter of the poisson Numerical Recipes in C, 2nd ed. pp. 293-295 */ static double sq, alm, g, oldm=-1; double em, t, y; /* search from the origin if (m<5) { if (m!=oldm) { oldm=m; g=exp(-m); } y=rndu(); sq=alm=g; for (em=0; ; ) { if (yt); } return ((int) em); } double rndgamma1 (double s); double rndgamma2 (double s); double rndgamma (double s) { /* random standard gamma (Mean=Var=s, with shape parameter=s, scale para=1) r^(s-1)*exp(-r) J. Dagpunar (1988) Principles of random variate generation, Clarendon Press, Oxford calling rndgamma1() if s<1 or rndgamma2() if s>1 or exponential if s=1 This is unsafe, and is found to return 0 when s is very small. */ double r=0; if (s<=0) puts ("jgl gamma.."); else if (s<1) r=rndgamma1 (s); else if (s>1) r=rndgamma2 (s); else r=-log(rndu()); return (r); } double rndgamma1 (double s) { /* random standard gamma for s<1 switching method */ double r, x=0,small=1e-37,w; static double a,p,uf,ss=10,d; if (s!=ss) { a=1-s; p=a/(a+s*exp(-a)); uf=p*pow(small/a,s); d=a*log(a); ss=s; } for (;;) { r=rndu(); if (r>p) x=a-log((1-r)/(1-p)), w=a*log(x)-d; else if (r>uf) x=a*pow(r/p,1/s), w=x; else return (0); r=rndu (); if (1-r<=w && r>0) if (r*(w+1)>=1 || -log(r)<=w) continue; break; } return (x); } double rndgamma2 (double s) { /* random standard gamma for s>1 Best's (1978) t distribution method */ double r,d,f,g,x; static double b,h,ss=0; if (s!=ss) { b=s-1; h=sqrt(3*s-0.75); ss=s; } for (;;) { r=rndu (); g=r-r*r; f=(r-0.5)*h/sqrt(g); x=b+f; if (x <= 0) continue; r=rndu(); d=64*r*r*g*g*g; if (d*x < x-2*f*f || log(d) < 2*(b*log(x/b)-f)) break; } return (x); } double rndbeta (double p, double q) { /* this generates a random beta(p,q) variate */ double gamma1, gamma2; gamma1=rndgamma(p); gamma2=rndgamma(q); return gamma1/(gamma1+gamma2); } int rndNegBinomial (double shape, double mean) { /* mean=mean, var=mean^2/shape+m */ return (rndpoisson(rndgamma(shape)/shape*mean)); } int MultiNomialAliasSetTable (int ncat, double prob[], double F[], int L[], double space[]) { /* This sets up the tables F and L for the alias algorithm for generating samples from the multinomial distribution MN(ncat, p) (Walker 1974; Kronmal & Peterson 1979). F[i] has cutoff probabilities, L[i] has aliases. I[i] is an indicator: -1 for F[i]<1; +1 for F[i]>=1; 0 if the cell is now empty. Should perhaps check whether prob[] sums to 1. */ signed char *I = (signed char *)space; int i,j,k, nsmall; for(i=0; i=1) I[i]=1; else { I[i]=-1; nsmall++; } } for(i=0; nsmall>0; i++) { for(j=0; j=1; 0 if the cell is now empty. */ int i,k; double r; for(i=0; i20), ncrude, lcrude[200]; double r, *pcdf=(space==NULL?prob:space), small=1e-5; ncrude=max2(5,ncat/20); ncrude=min2(200,ncrude); FOR(i,ncat) nobs[i]=0; if (space) { xtoy(prob, pcdf, ncat); for(i=1; i small) error2("sum P!=1 in MultiNomial2"); if (crude) { for(j=1,lcrude[0]=i=0; j 38, log(F(x)) can't be distinguished from 0. F(5) = 1 - 1.89E-8, and when x>5, double numbers lose digits, so the smaller tail area F(-5) is used for the calculation instead, using the expansion log(1-z) = -z(1 + z/2 + z*z/3), where z=F(-5) is small. For 3 < x < 7, both approaches are close, but when x = 8, Mathematica and log(CDFNormal) give the incorrect answer -6.66133815E-16, when the truth is log(F(8)) ~= -F(-8) = -6.22096057E-16. F(x) when x<-10 is reliably calculatd using the series expansion, even though log(CDFNormal) works until F(-38) = 2.88542835E-316. Regarding calculation of the logarithm of Pr(a < X < b), note that F(-9) - F(-10) = F(10) - F(9), but the left side is better computationally. */ double lnF, z=fabs(x), C, low=-10, high=5; /* calculate the log of the smaller area */ if(x>=low && x<=high) return log(CDFNormal(x)); if(x>high && x<-low) lnF = log(CDFNormal(-z)); else { C = 1 - 1/(z*z) + 3/(z*z*z*z) - 15/(z*z*z*z*z*z) + 105/(z*z*z*z*z*z*z*z); lnF = -z*z/2 - log(sqrt(2*Pi)*z) + log(C); } if(x>0) { z = exp(lnF); lnF = -z*(1 + z/2 + z*z/3 + z*z*z/4 + z*z*z*z/5); } return(lnF); } double LnGamma (double x) { /* returns ln(gamma(x)) for x>0, 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 f=0, fneg=0, z, lng; int nx=(int)x-1; if((double)nx==x && nx>0 && nx<10) lng=log((double)factorial(nx)); else { if(x<=0) { error2("lnGamma not implemented for x<0"); if((int)x-x==0) { puts("lnGamma undefined"); return(-1); } for (fneg=1; x<0; x++) fneg/=x; if(fneg<0) error2("strange!! check lngamma"); fneg=log(fneg); } if (x<7) { f=1; z=x-1; while (++z<7) f*=z; x=z; f=-log(f); } z = 1/(x*x); lng = fneg+ f + (x-0.5)*log(x) - x + .918938533204673 + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z +.083333333333333)/x; } return lng; } double PDFGamma (double x, double alpha, double beta) { /* gamma density: mean=alpha/beta; var=alpha/beta^2 */ if (x<=0 || alpha<=0 || beta<=0) { printf("x=%.6f a=%.6f b=%.6f", x, alpha, beta); error2("x a b outside range in PDFGamma()"); } if (alpha>100) error2("large alpha in PDFGamma()"); return pow(beta*x,alpha)/x * exp(-beta*x - LnGamma(alpha)); } double PDF_IGamma (double x, double alpha, double beta) { /* inverse-gamma density: mean=beta/(alpha-1); var=beta^2/[(alpha-1)^2*(alpha-2)] */ if (x<=0 || alpha<=0 || beta<=0) { printf("x=%.6f a=%.6f b=%.6f", x, alpha, beta); error2("x a b outside range in PDF_IGamma()"); } if (alpha>100) error2("large alpha in PDF_IGamma()"); return pow(beta/x,alpha)/x * exp(-beta/x - LnGamma(alpha)); } 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 accurate=1e-10, overflow=1e60; double factor, gin=0, rn=0, a=0,b=0,an=0,dif=0, term=0, pn[6]; if (x==0) return (0); if (x<0 || p<=0) return (-1); 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: return (gin); } double PointChi2 (double prob, double v) { /* returns z so that Prob{x1-small) return(9999); if (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) error2 ("\nIncompleteGamma"); 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 InverseGamma=0; /* 0 for gamma; 1 for inverse gamma */ double *gamma_modetiles, gamma_pexp[2]; /* gamma_modetiles[3]: mode, 2.5%, 97.5% */ double fun_ab_gamma (double x) { int i; double diff=0, a=x, b, p[2]={0.025,0.975}, *pe=gamma_pexp; if(a<=1) error2("a should be > 1"); if(!InverseGamma) b=(a-1)/gamma_modetiles[0]; else b=(a+1)/gamma_modetiles[0]; for(i=0; i<2; i++) pe[i] = CDFGamma(gamma_modetiles[i+1],a,b); diff = fabs(p[0]+1-p[1]- (pe[0]+1-pe[1])) * 2 + fabs(p[0]-pe[0]); /* diff = fabs(p[0]-pe[0]) + fabs(p[1]-pe[1]); */ return(diff); } int getab_gamma (double *a, double *b, double xmodetiles[3]) { /* This gets parameters a (alpha) and b (beta) for the gamma distribution, given the mode, and the 2.5% and 97.5% percentiles. This is used to construct the gamma prior for fossil date using soft bounds. The mode is matched perfectly so that b = (a-1)/mode, and a is solved to minimize the differences in CDFs at the soft bounds. a > 1 is assumed. */ int noisy0=noisy; double x=InverseGamma+1+1+2*rndu(), xb[2]={1.01,999}, S, e=1e-15; if(xmodetiles[0]<=xmodetiles[1] || xmodetiles[0]>xmodetiles[2]) error2("modetiles for gamma out of range."); noisy=0; gamma_modetiles = xmodetiles; LineSearch(fun_ab_gamma, &S, &x, xb, 0.2, e); *a=x; if(!InverseGamma) *b=(*a-1)/gamma_modetiles[0]; else *b=(*a+1)/gamma_modetiles[0]; noisy=noisy0; return(0); } int DiscreteGamma (double freqK[], double rK[], double alpha, double beta, int K, int mean) { /* discretization of gamma distribution with equal proportions in each category. */ int i; double t, factor=alpha/beta*K, lnga1; if (mean) { lnga1=LnGamma(alpha+1); for (i=0; i1-1e-4) error2("rho out of range"); */ FOR (i,K-1) point[i]=PointNormal((i+1.0)/K); for (i=0; i0) y-= M[(i-1)*K+j]; if (j>0) y-= M[i*K+(j-1)]; if (i>0 && j>0) y += M[(i-1)*K+(j-1)]; M[i*K+j] = (M[i*K+j]+y)*K; if (M[i*K+j]<0) printf("M(%d,%d) =%12.8f<0\n", i+1, j+1, M[i*K+j]); } } DiscreteGamma (freqK, rK, alpha, alpha, K, DGammaMean); for (i=0,v1=*rho1=0; ih, y>k), where x and y are standard binormal, with r = corr(x,y). (1) Drezner Z., and G.O. Wesolowsky (1990) On the computation of the bivariate normal integral. J. Statist. Comput. Simul. 35:101-107. (2) Genz, A.C., Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statist. Comput., 2004. 14: p. 1573-1375. This uses the algorithm of Genz (2004) for r>0, and L(h,k,r) = F(-k) - L(-h, k, -r), if r<0. Here hk, a swapping between h and k is performed. Gauss-Legendre quadrature points used. |r| Genz nGL <0.3 (eq. 3) 6 10 <0.75 (eq. 3) 12 20 <0.925 (eq. 3) 20 20 <1 (eq. 6) 20 20 */ int nGL = (fabs(r)<0.3 ? 10 : 20), i,j, signr = (r>=0 ? 1 : -1); double *x, *w; /* Gauss-Legendre quadrature points */ double hk=h*k, h0=h,k0=k, L=0, t[2], hk2, y, a=0,b,c,d, bs, as, rs, smallr=1e-10; h=min2(h0,k0); k=max2(h0,k0); if(signr==-1) { r=-r; h=-h; hk=h*k; } if(r>1) error2("|r| > 1 in LBinormal"); GaussLegendreRule(&x, &w, nGL); if(r < 0.925) { /* equation 3 */ if(r>smallr) { hk2 = (h*h + k*k)/2; a = asin(r)/2; for(i=0,L=0; i-1+smallr && r<1-smallr) { /* first term in equation (6), analytical */ as = 1-r*r; a = sqrt(as); b = fabs(h - k); bs = b*b; c = (4 - hk)/8 ; d = (12 - hk)/16; y = -(bs/as + hk)/2; if(y > -500) L = a*exp(y)*(1 - c*(bs-as)*(1-d*bs/5)/3 + c*d*as*as/5); if(hk > -500) { L -= exp(-hk/2)*sqrt(2*Pi) * CDFNormal(-b/a) * b * (1 - c*bs*(1 - d*bs/5)/3); } /* second term in equation (6), numerical */ a /= 2; for(i=0; i -500) L += a*w[i]*exp(y)*(exp(-hk*(1-rs)/(2*(1+rs)))/rs - (1+c*t[j]*(1+d*t[j]))); } } } L = CDFNormal(-max2(h, k)) - L/(2*Pi); } if(signr<0) L = CDFNormal(-k) - L; if(L<-1e-12) printf("L = %.9g very negative. Let me know please.\n", L); if(L<0) L=0; return (L); } double logLBinormal (double h, double k, double r) { /* this is to calculate the logarithm of tail probabilities of bivariate normal distribution, modified from LBinormal(). When r < 0, we use L(h,k,r) = F(-k) - L(-h, k, -r), where h < k. In other words, L(-10, 9, -1) = F(-9) - F(-10) is better than L(-10, 9, -1) = F(10) - F(9). So we use L(-10, 9, -0.3) = L(-9) - L(10, 9, 0.3). not L(-10, 9, -0.3) = L(10) - L(-10, -9, 0.3). Results for the left tail, the very negative log(L), are reliable. Results for the right tail are not reliable, that is, if log(L) is close to 0 and L ~= 1. Perhaps try to use the following to reflect: L(-5,-9,r) = 1 - [ F(-5) + F(-9) - L(5,9,r) ] Check the routine logCDFNormal() to see the idea. */ int nGL = (fabs(r)<0.3 ? 10 : 20), i,j, signr = (r>=0 ? 1 : -1); double *x, *w; /* Gauss-Legendre quadrature points */ double hk=h*k, h0=h,k0=k, L, t[2], hk2, a,b,c,d, bs, as, rs; double S1,S2=-1e300,S3=-1e300, y,L2=0,L3=0, largeneg=-1e300, smallr=1e-10; h=min2(h0,k0); k=max2(h0,k0); if(signr==-1) { r=-r; h=-h; hk=-hk; } if(fabs(r)>1+smallr) error2("|r| > 1 in LBinormal"); GaussLegendreRule(&x, &w, nGL); if(r < 0.925) { /* equation 3 */ S1 = L = logCDFNormal(-h) + logCDFNormal(-k); if(r>smallr) { hk2 = (h*h + k*k)/2; a = asin(r)/2; for(i=0,L=0,S2=-hk2; iS2+30) { L *= exp(S2-y); S2 = y; } L += w[i]*exp(y-S2); } } y = max2(S1,S2); a = exp(S1-y) + L*a/(2*Pi)*exp(S2-y); L = (a>0 ? y + log(a) : largeneg); } } else { /* equation 6, using equation 7 instead of equation 5. */ /* log(L) = exp(S1) + L2*exp(S2) + L3*exp(S3) */ if(r>-1+smallr && r<1-smallr) { /* first term in equation (6), analytical: L2 & S2 */ as = 1-r*r; a = sqrt(as); b = fabs(h - k); bs = b*b; c = (4 - hk)/8; d = (12 - hk)/16; S2 = -(bs/as + hk)/2; /* is this too large? */ L2 = a*(1 - c*(bs-as)*(1-d*bs/5)/3 + c*d*as*as/5); y = -hk/2 + logCDFNormal(-b/a); if(y>S2+30) { L2 *= exp(S2-y); S2 = y; } L2 -= sqrt(2*Pi) * exp(y-S2) * b * (1 - c*bs*(1 - d*bs/5)/3); /* second term in equation (6), numerical: L3 & S3 */ a /= 2; for(i=0,L3=0,S3=-1e300; i S3+30) { L3 *= exp(S3-y); S3 = y; } L3 += a*w[i]*exp(y-S3) * (exp(-hk*(1-rs)/(2*(1+rs)))/rs - (1+c*t[j]*(1+d*t[j]))); } } } S1 = logCDFNormal(-max2(h, k)); y = max2(S1,S2); y = max2(y,S3); a = exp(S1-y) - (L2*exp(S2-y) + L3*exp(S3-y))/(2*Pi); L = (a>0 ? y + log(a) : largeneg); } if(signr==-1) { /* r<0: L(h,k,-r) = F(-k) - L(-h,k,r) */ S1 = logCDFNormal(-k); y = max2(S1,L); a = exp(S1-y) - exp(L-y); L = (a > 0 ? y + log(a) : largeneg); } if(L>1e-12) printf("ln L(%2g, %.2g, %.2g) = %.6g is very large.\n", h0,k0,r*signr, L); if(L>0) L=0; return(L); } #if (0) void testLBinormal (void) { double x,y,r, L, lnL; for(r=-1; r<1.01; r+=0.05) { if(fabs(r-1)<1e-5) r=1; printf("\nr = %.2f\n", r); for(x=-10; x<=10.01; x+=0.5) { for(y=-10; y<=10.01; y+=0.5) { printf("x y r? "); scanf("%lf%lf%lf", &x, &y, &r); L = LBinormal(x,y,r); lnL = logLBinormal(x,y,r); if(fabs(L-exp(lnL))>1e-10) printf("L - exp(lnL) = %.10g very different.\n", L - exp(lnL)); // if(L<0 || L>1) printf("%6.2f %6.2f %6.2f L %15.8g = %15.8g %9.5g\n", x,y,r, L, exp(lnL), lnL); if(lnL>0) exit(-1); } } } } #endif double probBinomial (int n, int k, double p) { /* calculates {n\choose k} * p^k * (1-p)^(n-k) */ double C, up, down; if (n<40 || (n<1000&&k<10)) { for (down=min2(k,n-k),up=n,C=1; down>0; down--,up--) C*=up/down; if (fabs(p-.5)<1e-6) C *= pow(p,(double)n); else C *= pow(p,(double)k)*pow((1-p),(double)(n-k)); } else { C = exp((LnGamma(n+1.)-LnGamma(k+1.)-LnGamma(n-k+1.))/n); C = pow(p*C,(double)k) * pow((1-p)*C,(double)(n-k)); } return C; } double probBetaBinomial (int n, int k, double p, double q) { /* This calculates beta-binomial probability of k succeses out of n trials, The binomial probability parameter has distribution beta(p, q) prob(x) = C1(-a,k) * C2(-b,n-k)/C3(-a-b,n) */ double a=p,b=q, C1,C2,C3,scale1,scale2,scale3; if(a<=0 || b<=0) return(0); C1=Binomial(-a, k, &scale1); C2=Binomial(-b, n-k, &scale2); C3=Binomial(-a-b, n, &scale3); C1*=C2/C3; if(C1<0) error2("error in probBetaBinomial"); return C1*exp(scale1+scale2-scale3); } double PDFBeta(double x, double p, double q) { /* Returns pdf of beta(p,q) */ double y, small=1e-20; if(x1-small) error2("bad x in PDFbeta"); y = (p-1)*log(x) + (q-1)*log(1-x); y-= LnGamma(p)+LnGamma(q)-LnGamma(p+q); return(exp(y)); } double CDFBeta(double x, double pin, double qin, double lnbeta) { /* Returns distribution function of the standard form of the beta distribution, that is, the incomplete beta ratio I_x(p,q). lnbeta is log of the complete beta function; provide it if known, and otherwise use 0. This is called from InverseCDFBeta() in a root-finding loop. This routine is a translation into C of a Fortran subroutine by W. Fullerton of Los Alamos Scientific Laboratory. Bosten and Battiste (1974). Remark on Algorithm 179, CACM 17, p153, (1974). */ double ans, c, finsum, p, ps, p1, q, term, xb, xi, y, small=1e-15; int n, i, ib; static double eps = 0, alneps = 0, sml = 0, alnsml = 0; if(x1-small) return 1; if(pin<=0 || qin<=0) { printf("p=%.4f q=%.4f: parameter outside range in CDFBeta",pin,qin); return (-1); } if (eps == 0) {/* initialize machine constants ONCE */ eps = pow((double)FLT_RADIX, -(double)DBL_MANT_DIG); alneps = log(eps); sml = DBL_MIN; alnsml = log(sml); } y = x; p = pin; q = qin; /* swap tails if x is greater than the mean */ if (p / (p + q) < x) { y = 1 - y; p = qin; q = pin; } if(lnbeta==0) lnbeta=LnGamma(p)+LnGamma(q)-LnGamma(p+q); if ((p + q) * y / (p + 1) < eps) { /* tail approximation */ ans = 0; xb = p * log(max2(y, sml)) - log(p) - lnbeta; if (xb > alnsml && y != 0) ans = exp(xb); if (y != x || p != pin) ans = 1 - ans; } else { /* evaluate the infinite sum first. term will equal */ /* y^p / beta(ps, p) * (1 - ps)-sub-i * y^i / fac(i) */ ps = q - floor(q); if (ps == 0) ps = 1; xb=LnGamma(ps)+LnGamma(p)-LnGamma(ps+p); xb = p * log(y) - xb - log(p); ans = 0; if (xb >= alnsml) { ans = exp(xb); term = ans * p; if (ps != 1) { n = (int)max2(alneps/log(y), 4.0); for(i=1 ; i<= n ; i++) { xi = i; term = term * (xi - ps) * y / xi; ans = ans + term / (p + xi); } } } /* evaluate the finite sum. */ if (q > 1) { xb = p * log(y) + q * log(1 - y) - lnbeta - log(q); ib = (int) (xb/alnsml); if(ib<0) ib=0; term = exp(xb - ib * alnsml); c = 1 / (1 - y); p1 = q * c / (p + q - 1); finsum = 0; n = (int) q; if (q == (double)n) n = n - 1; for(i=1 ; i<=n ; i++) { if (p1 <= 1 && term / eps <= finsum) break; xi = i; term = (q - xi + 1) * c * term / (p + q - xi); if (term > 1) { ib = ib - 1; term = term * sml; } if (ib == 0) finsum = finsum + term; } ans = ans + finsum; } if (y != x || p != pin) ans = 1 - ans; if(ans>1) ans=1; if(ans<0) ans=0; } return ans; } double InverseCDFBeta(double prob, double p, double q, double lnbeta) { /* This calculates the inverseCDF of the beta distribution Cran, G. W., K. J. Martin and G. E. Thomas (1977). Remark AS R19 and Algorithm AS 109, Applied Statistics, 26(1), 111-114. Remark AS R83 (v.39, 309-310) and correction (v.40(1) p.236). My own implementation of the algorithm did not bracket the variable well. This version is Adpated from the pbeta and qbeta routines from "R : A Computer Language for Statistical Data Analysis". It fails for extreme values of p and q as well, although it seems better than my previous version. Ziheng Yang, May 2001 */ double fpu=3e-308, acu_min=1e-300, lower=fpu, upper=1-2.22e-16; /* acu_min>= fpu: Minimal value for accuracy 'acu' which will depend on (a,p); */ int swap_tail, i_pb, i_inn, niterations=2000; double a, adj, g, h, pp, prev=0, qq, r, s, t, tx=0, w, y, yprev; double acu, xinbta; if(prob<0 || prob>1 || p<0 || q<0) error2("out of range in InverseCDFBeta"); /* define accuracy and initialize */ xinbta = prob; /* test for admissibility of parameters */ if(p<0 || q<0 || prob<0 || prob>1) error2("beta par err"); if (prob == 0 || prob == 1) return prob; if(lnbeta==0) lnbeta=LnGamma(p)+LnGamma(q)-LnGamma(p+q); /* change tail if necessary; afterwards 0 < a <= 1/2 */ if (prob <= 0.5) { a = prob; pp = p; qq = q; swap_tail = 0; } else { a = 1. - prob; pp = q; qq = p; swap_tail = 1; } /* calculate the initial approximation */ r = sqrt(-log(a * a)); y = r - (2.30753+0.27061*r)/(1.+ (0.99229+0.04481*r) * r); if (pp > 1. && qq > 1.) { r = (y * y - 3.) / 6.; s = 1. / (pp*2. - 1.); t = 1. / (qq*2. - 1.); h = 2. / (s + t); w = y * sqrt(h + r) / h - (t - s) * (r + 5./6. - 2./(3.*h)); xinbta = pp / (pp + qq * exp(w + w)); } else { r = qq*2.; t = 1. / (9. * qq); t = r * pow(1. - t + y * sqrt(t), 3.); if (t <= 0.) xinbta = 1. - exp((log((1. - a) * qq) + lnbeta) / qq); else { t = (4.*pp + r - 2.) / t; if (t <= 1.) xinbta = exp((log(a * pp) + lnbeta) / pp); else xinbta = 1. - 2./(t+1.); } } /* solve for x by a modified newton-raphson method, using CDFBeta */ r = 1. - pp; t = 1. - qq; yprev = 0.; adj = 1.; /* Changes made by Ziheng to fix a bug in qbeta() qbeta(0.25, 0.143891, 0.05) = 3e-308 wrong (correct value is 0.457227) */ if(xinbta<=lower || xinbta>=upper) xinbta=(a+.5)/2; /* Desired accuracy should depend on (a,p) * This is from Remark .. on AS 109, adapted. * However, it's not clear if this is "optimal" for IEEE double prec. * acu = fmax2(acu_min, pow(10., -25. - 5./(pp * pp) - 1./(a * a))); * NEW: 'acu' accuracy NOT for squared adjustment, but simple; * ---- i.e., "new acu" = sqrt(old acu) */ acu = pow(10., -13. - 2.5/(pp * pp) - 0.5/(a * a)); acu = max2(acu, acu_min); for (i_pb=0; i_pb= 0. && tx <= 1.) { if (prev <= acu || fabs(y) <= acu) goto L_converged; if (tx != 0. && tx != 1.) break; } } g /= 3.; } if (fabs(tx-xinbta)=xb[1]) x=.5; LineSearch(diff_InverseCDF, &sdiff, &x, xb, step, e); noisy=noisy0; return(x); } int GaussLegendreRule(double **x, double **w, int npoints) { /* this returns the Gauss-Legendre nodes and weights in x[] and w[]. */ int status=0; static double x10[]={0.148874338981631210884826001130, 0.433395394129247190799265943166, 0.679409568299024406234327365115, 0.865063366688984510732096688423, 0.973906528517171720077964012084}; static double w10[]={0.295524224714752870173892994651, 0.269266719309996355091226921569, 0.219086362515982043995534934228, 0.149451349150580593145776339658, 0.066671344308688137593568809893}; static double x20[]={0.076526521133497333754640409399, 0.227785851141645078080496195369, 0.373706088715419560672548177025, 0.510867001950827098004364050955, 0.636053680726515025452836696226, 0.746331906460150792614305070356, 0.839116971822218823394529061702, 0.912234428251325905867752441203, 0.963971927277913791267666131197, 0.993128599185094924786122388471}; static double w20[]={0.152753387130725850698084331955, 0.149172986472603746787828737002, 0.142096109318382051329298325067, 0.131688638449176626898494499748, 0.118194531961518417312377377711, 0.101930119817240435036750135480, 0.083276741576704748724758143222, 0.062672048334109063569506535187, 0.040601429800386941331039952275, 0.017614007139152118311861962352}; static double x32[]={0.048307665687738316234812570441, 0.144471961582796493485186373599, 0.239287362252137074544603209166, 0.331868602282127649779916805730, 0.421351276130635345364119436172, 0.506899908932229390023747474378, 0.587715757240762329040745476402, 0.663044266930215200975115168663, 0.732182118740289680387426665091, 0.794483795967942406963097298970, 0.849367613732569970133693004968, 0.896321155766052123965307243719, 0.934906075937739689170919134835, 0.964762255587506430773811928118, 0.985611511545268335400175044631, 0.997263861849481563544981128665}; static double w32[]={0.0965400885147278005667648300636, 0.0956387200792748594190820022041, 0.0938443990808045656391802376681, 0.0911738786957638847128685771116, 0.0876520930044038111427714627518, 0.0833119242269467552221990746043, 0.0781938957870703064717409188283, 0.0723457941088485062253993564785, 0.0658222227763618468376500637069, 0.0586840934785355471452836373002, 0.0509980592623761761961632446895, 0.0428358980222266806568786466061, 0.0342738629130214331026877322524, 0.0253920653092620594557525897892, 0.0162743947309056706051705622064, 0.0070186100094700966004070637389}; static double x64[]={0.024350292663424432508955842854, 0.072993121787799039449542941940, 0.121462819296120554470376463492, 0.169644420423992818037313629748, 0.217423643740007084149648748989, 0.264687162208767416373964172510, 0.311322871990210956157512698560, 0.357220158337668115950442615046, 0.402270157963991603695766771260, 0.446366017253464087984947714759, 0.489403145707052957478526307022, 0.531279464019894545658013903544, 0.571895646202634034283878116659, 0.611155355172393250248852971019, 0.648965471254657339857761231993, 0.685236313054233242563558371031, 0.719881850171610826848940217832, 0.752819907260531896611863774886, 0.783972358943341407610220525214, 0.813265315122797559741923338086, 0.840629296252580362751691544696, 0.865999398154092819760783385070, 0.889315445995114105853404038273, 0.910522137078502805756380668008, 0.929569172131939575821490154559, 0.946411374858402816062481491347, 0.961008799652053718918614121897, 0.973326827789910963741853507352, 0.983336253884625956931299302157, 0.991013371476744320739382383443, 0.996340116771955279346924500676, 0.999305041735772139456905624346}; static double w64[]={0.0486909570091397203833653907347, 0.0485754674415034269347990667840, 0.0483447622348029571697695271580, 0.0479993885964583077281261798713, 0.0475401657148303086622822069442, 0.0469681828162100173253262857546, 0.0462847965813144172959532492323, 0.0454916279274181444797709969713, 0.0445905581637565630601347100309, 0.0435837245293234533768278609737, 0.0424735151236535890073397679088, 0.0412625632426235286101562974736, 0.0399537411327203413866569261283, 0.0385501531786156291289624969468, 0.0370551285402400460404151018096, 0.0354722132568823838106931467152, 0.0338051618371416093915654821107, 0.0320579283548515535854675043479, 0.0302346570724024788679740598195, 0.0283396726142594832275113052002, 0.0263774697150546586716917926252, 0.0243527025687108733381775504091, 0.0222701738083832541592983303842, 0.0201348231535302093723403167285, 0.0179517157756973430850453020011, 0.0157260304760247193219659952975, 0.0134630478967186425980607666860, 0.0111681394601311288185904930192, 0.0088467598263639477230309146597, 0.0065044579689783628561173604000, 0.0041470332605624676352875357286, 0.0017832807216964329472960791450}; if(npoints==10) { *x=x10; *w=w10; } else if(npoints==20) { *x=x20; *w=w20; } else if(npoints==32) { *x=x32; *w=w32; } else if(npoints==64) { *x=x64; *w=w64; } else { puts("use 10, 20, 32 or 64 nodes for legendre."); status=-1; } return(status); } double NIntegrateGaussLegendre(double(*fun)(double x), double a, double b, int npoints) { /* this approximates the integral Nintegrate[fun[x], {x,a,b}]. npoints is 10, 20, 32 or 64 nodes for legendre."); */ int i,j; double *x, *w, s, t[2]; GaussLegendreRule(&x, &w, npoints); for(j=0,s=0; j1. */ char *chart,ch, *fmt[2]={"%*.*e ", "%*.*f "}, symbol[]="*~^@",overlap='&'; int i,j,is,iy,ny=1, ncolr=ncol+3, irow=0, icol=0, w=10, wd=2; double large=1e32, xmin, xmax, xgap, ymin[2], ymax[2], ygap[2]; for (i=1,xmin=xmax=x[0]; ix[i]) xmin=x[i]; if(xmax1)*nseries; j++) if (yLorR[j]==1) ny=2; else if (yLorR[j]!=0) printf ("err: y axis %d", yLorR[j]); for (j=0; jy[j*n+i]) ymin[iy]=y[j*n+i]; if (ymax[iy]maxT) error2("temperature rescaling needed."); *R = (int)fabs( -5157.3*T*T*T*T + 9681.4*T*T*T - 5491.9*T*T + 1137.7*T + 6.2168 ); *G = (int)fabs( -1181.4*T*T*T + 964.8*T*T + 203.66*T + 1.2028 ); *B = (int)fabs( 92.463*T*T*T - 595.92*T*T + 481.11*T + 21.769 ); if(*R>255) *R=255; if(*G>255) *G=255; if(*B>255) *B=255; } void GetIndexTernary(int *ix, int *iy, double *x, double *y, int itriangle, int K) { /* This gives the indices (ix, iy) and the coordinates (x, y, 1-x-y) for the itriangle-th triangle, with itriangle from 0, 1, ..., KK-1. The ternary graph (0-1 on each axis) is partitioned into K*K equal-sized triangles. In the first row (ix=0), there is one triangle (iy=0); In the second row (ix=1), there are 3 triangles (iy=0,1,2); In the i-th row (ix=i), there are 2*i+1 triangles (iy=0,1,...,2*i). x rises when ix goes up, but y decreases when iy increases. (x,y) is the centroid in the ij-th small triangle. x and y each takes on 2*K-1 possible values. */ *ix = (int)sqrt((double)itriangle); *iy = itriangle - square(*ix); *x = (1 + (*iy/2)*3 + (*iy%2))/(3.*K); *y = (1 + (K-1- *ix)*3 + (*iy%2))/(3.*K); } long factorial(int n) { long f, i; if (n>10) error2("n>10 in factorial"); for (i=2,f=1; i<=(long)n; i++) f*=i; return (f); } double Binomial(double n, int k, double *scale) { /* calculates n choose k. n is any real number, and k is integer. If(*scale!=0) the result should be c+exp(*scale). */ double c=1,i,large=1e99; *scale=0; if((int)k!=k) error2("k is not a whole number in Binomial."); if(n<0 && k%2==1) c=-1; if(k==0) return(1); if(n>0 && (k<0 || k>n)) return (0); if(n>0 && (int)n==n) k=min2(k,(int)n-k); for (i=1; i<=k; i++) { c*=(n-k+i)/i; if(c>large) { *scale+=log(c); c=1; } } return(c); } /**************************** Vectors and matrices *****************************/ double Det3x3 (double x[3*3]) { return x[0*3+0]*x[1*3+1]*x[2*3+2] + x[0*3+1]*x[1*3+2]*x[2*3+0] + x[0*3+2]*x[1*3+0]*x[2*3+1] - x[0*3+0]*x[1*3+2]*x[2*3+1] - x[0*3+1]*x[1*3+0]*x[2*3+2] - x[0*3+2]*x[1*3+1]*x[2*3+0] ; } int matby (double a[], double b[], double c[], int n,int m,int k) /* a[n*m], b[m*k], c[n*k] ...... c = a*b */ { int i,j,i1; double t; FOR (i,n) FOR(j,k) { for (i1=0,t=0; i1 y[m][n] */ int i,j; FOR (i,n) FOR (j,m) y[j*n+i]=x[i*m+j]; return (0); } int matinv (double x[], int n, int m, double space[]) { /* x[n*m] ... m>=n space[n]. This puts the fabs(|x|) into space[0]. Check and calculate |x|. Det may have the wrong sign. Check and fix. */ int i,j,k; int *irow=(int*) space; double ee=1e-100, t,t1,xmax, det=1; for(i=0; i=0; i--) { if (irow[i] == i) continue; FOR(j,n) { t = x[j*m+i]; x[j*m+i] = x[j*m + irow[i]]; x[j*m + irow[i]] = t; } } space[0]=det; return(0); } int matexp (double Q[], double t, int n, int TimeSquare, double space[]) { /* This calculates the matrix exponential P(t) = exp(t*Q). Input: Q[] has the rate matrix, and t is the time or branch length. TimeSquare is the number of times the matrix is squared and should be from 5 to 31. Output: Q[] has the transition probability matrix, that is P(Qt). space[n*n]: required working space. P(t) = (I + Qt/m + (Qt/m)^2/2)^m, with m = 2^TimeSquare. T[it=0] is the current matrix, and T[it=1] is the squared result matrix, used to avoid copying matrices. Use an even TimeSquare to avoid one round of matrix copying. */ int it, i; double *T[2]; if(TimeSquare<2 || TimeSquare>31) error2("TimeSquare not good"); T[0]=Q; T[1]=space; for(i=0; i=i) are used. */ int i,j,k; double t; for (i=0; i=0) L[i*n+i]=sqrt(t); else return (-1); for (j=i+1; j=j) are used */ int i,j; double t; for (i=0; i=0; i--) { /* solve L'x=y, and store results in x */ for (j=i+1,t=x[i]; jsmall) pi_sqrt[nnew++]=sqrt(pi[j]); /* store in U the symmetrical matrix S = sqrt(D) * Q * sqrt(-D) */ if(nnew==n) { for(i=0; ismall) { for(j=0,jnew=0; jsmall) { U[inew*nnew+jnew] = U[jnew*nnew+inew] = Q[i*n+j] * pi_sqrt[inew]/pi_sqrt[jnew]; jnew++; } U[inew*nnew+inew] = Q[i*n+i]; inew++; } } status=eigenRealSym(U, nnew, Root, V); for(i=n-1,inew=nnew-1; i>=0; i--) /* construct Root */ Root[i] = (pi[i]>small ? Root[inew--] : 0); for(i=n-1,inew=nnew-1; i>=0; i--) { /* construct V */ if(pi[i]>small) { for(j=n-1,jnew=nnew-1; j>=0; j--) if(pi[j]>small) { V[i*n+j] = U[jnew*nnew+inew]*pi_sqrt[jnew]; jnew--; } else V[i*n+j] = (i==j); inew--; } else for(j=0; j=0; i--) { /* construct U */ if(pi[i]>small) { for(j=n-1,jnew=nnew-1;j>=0;j--) if(pi[j]>small) { U[i*n+j] = U[inew*nnew+jnew]/pi_sqrt[inew]; jnew--; } else U[i*n+j] = (i==j); inew--; } else for(j=0;j1e-10 && noisy) printf("Root[0] = %.5e\n",Root[0]); Root[0]=0; */ return(status); } /* eigen solution for real symmetric matrix */ void HouseholderRealSym(double a[], int n, double d[], double e[]); int EigenTridagQLImplicit(double d[], double e[], int n, double z[]); void EigenSort(double d[], double U[], int n); int eigenRealSym(double A[], int n, double Root[], double work[]) { /* This finds the eigen solution of a real symmetrical matrix A[n*n]. In return, A has the right vectors and Root has the eigenvalues. work[n] is the working space. The matrix is first reduced to a tridiagonal matrix using HouseholderRealSym(), and then using the QL algorithm with implicit shifts. Adapted from routine tqli in Numerical Recipes in C, with reference to LAPACK Ziheng Yang, 23 May 2001 */ int status=0; HouseholderRealSym(A, n, Root, work); status=EigenTridagQLImplicit(Root, work, n, A); EigenSort(Root, A, n); return(status); } void EigenSort(double d[], double U[], int n) { /* this sorts the eigen values d[] and rearrange the (right) eigen vectors U[] */ int k,j,i; double p; for (i=0;i= p) p=d[k=j]; if (k != i) { d[k]=d[i]; d[i]=p; for (j=0;j=1;i--) { m=i-1; h=scale=0; if (m > 0) { for (k=0;k<=m;k++) scale += fabs(a[i*n+k]); if (scale == 0) e[i]=a[i*n+m]; else { for (k=0;k<=m;k++) { a[i*n+k] /= scale; h += a[i*n+k]*a[i*n+k]; } f=a[i*n+m]; g=(f >= 0 ? -sqrt(h) : sqrt(h)); e[i]=scale*g; h -= f*g; a[i*n+m]=f-g; f=0; for (j=0;j<=m;j++) { a[j*n+i]=a[i*n+j]/h; g=0; for (k=0;k<=j;k++) g += a[j*n+k]*a[i*n+k]; for (k=j+1;k<=m;k++) g += a[k*n+j]*a[i*n+k]; e[j]=g/h; f += e[j]*a[i*n+j]; } hh=f/(h*2); for (j=0;j<=m;j++) { f=a[i*n+j]; e[j]=g=e[j]-hh*f; for (k=0;k<=j;k++) a[j*n+k] -= (f*e[k]+g*a[i*n+k]); } } } else e[i]=a[i*n+m]; d[i]=h; } d[0]=e[0]=0; /* Get eigenvectors */ for (i=0;i= 0.0 ? fabs(a) : -fabs(a)) int EigenTridagQLImplicit(double d[], double e[], int n, double z[]) { /* This finds the eigen solution of a tridiagonal matrix represented by d and e. d[] is the diagonal (eigenvalues), e[] is the off-diagonal z[n*n]: as input should have the identity matrix to get the eigen solution of the tridiagonal matrix, or the output from HouseholderRealSym() to get the eigen solution to the original real symmetric matrix. z[n*n]: has the orthogonal matrix as output Adapted from routine tqli in Numerical Recipes in C, with reference to LAPACK fortran code. Ziheng Yang, May 2001 */ int m,j,iter,niter=30, status=0, i,k; double s,r,p,g,f,dd,c,b, aa,bb; for (i=1;i1) r=aa*sqrt(1+1/(g*g)); else r=sqrt(1+g*g); g=d[m]-d[j]+e[j]/(g+SIGN(r,g)); s=c=1; p=0; for (i=m-1;i>=j;i--) { f=s*e[i]; b=c*e[i]; /* r=pythag(f,g); */ aa=fabs(f); bb=fabs(g); if(aa>bb) { bb/=aa; r=aa*sqrt(1+bb*bb); } else if(bb==0) r=0; else { aa/=bb; r=bb*sqrt(1+aa*aa); } e[i+1]=r; if (r == 0) { d[i+1] -= p; e[m]=0; break; } s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2*c*b; d[i+1]=g+(p=s*r); g=c*r-b; for (k=0;k= j) continue; d[j]-=p; e[j]=g; e[m]=0; } } while (m != j); } return(status); } #undef SIGN int MeanVar (double x[], int n, double *m, double *v) { int i; for (i=0,*m=0; i1) *v/=(n-1.); return(0); } int variance (double x[], int n, int nx, double mx[], double vx[]) { /* x[nx][n], mx[nx], vx[nx][nx] */ int i, j, k; FOR(i,nx) mx[i]=0; FOR(i,nx) FOR (k,n) mx[i]=(mx[i]*k+x[i*n+k])/(k+1.); FOR (i, nx*nx) vx[i]=0; for (i=0; i0.0 && *v22>0.0) *r=*v12/sqrt(*v11 * *v22); else *r=-9; return(0); } int bubblesort (float x[], int n) { /* inefficient bubble sort */ int i,j; float t=0; for(i=0;i bb ? 1 : -1); } int comparedouble (const void *a, const void *b) { double aa = *(double*)a, bb= *(double*)b; return (aa==bb ? 0 : aa > bb ? 1 : -1); } int splitline (char line[], int fields[]) { /* This finds out how many fields there are in the line, and marks the starting positions of the fields. Fields are separated by spaces, and texts are allowed as well. */ int lline=64000, i, nfields=0, InSpace=1; char *p=line; for(i=0; iMAXNFIELDS) puts("raise MAXNFIELDS?"); */ } } } return(nfields); } int scanfile(FILE*fin, int *nrecords, int *nx) { int lline=64000, nxline, fields[1000]; char line[64000]; for (*nrecords=0; ; ) { if (!fgets(line,lline,fin)) break; nxline=splitline(line,fields); if(nxline==0) continue; if(*nrecords==0) *nx=nxline; else if (*nx!=nxline) break; (*nrecords)++; /* printf("line # %3d: %3d variables\n", *nrecords+1, nxline); */ } rewind(fin); return(0); } #define MAXNF2D 5 #define SQRT5 2.2360679774997896964 #define Epanechnikov(t) ((0.75-0.15*(t)*(t))/SQRT5) int splitline (char line[], int fields[]); int scanfile(FILE*fin, int *nrecords, int *nx); /* density1d and density2d need to be reworked to account for edge effects. October 2006 */ int density1d (FILE* fout, double y[], int n, int nbin, double minx, double gap, double h, double space[], int propternary) { /* This collects the histogram and uses kernel smoothing and adaptive kernel smoothing to estimate the density. The kernel is Epanechnikov. The histogram is collected into fobs, smoothed density into fE, and density from adaptive kernel smoothing into fEA[]. Data y[] are sorted in increasing order before calling this routine. space[bin+n] */ int adaptive=0, i,k, iL, nused; double *fobs=space, *lambda=fobs+nbin, fE, fEA, xt, d, G, alpha=0.5; double maxx=minx+gap*nbin, edge; char timestr[32]; for(i=0; i1000) printf("\r"); } /* smoothing and printing */ for (k=0; kSQRT5) continue; edge = (y[i]-xt > xt-minx || xt-y[i]>maxx-xt ? 2 : 1); fE += edge*Epanechnikov(d)/(n*h); if(adaptive) { d/=lambda[i]; fEA += edge*Epanechnikov(d)/(n*h*lambda[i]); } } fprintf(fout, "%.6f\t%.6f\t%.6f\t%.6f\n", xt, fobs[k], fE, fEA); } return(0); } int density2d (FILE* fout, double y1[], double y2[], int n, int nbin, double minx1, double minx2, double gap1, double gap2, double var[4], double h, double space[], int propternary) { /* This collects the histogram and uses kernel smoothing and adaptive kernel smoothing to estimate the 2-D density. The kernel is Epanechnikov. The histogram is collected into f, smoothed density into fE, and density from adaptive kernel smoothing into fEA[]. Data y1 and y2 are not sorted, unlike the 1-D routine. alpha goes from 0 to 1, with 0 being equivalent to fixed width smoothing. var[] has the 2x2 variance matrix, which is copied into S[4] and inverted. space[nbin*nbin+n] for observed histogram f[nbin*nbin] and for lambda[n]. */ char timestr[32]; int i,j,k; double *fO=space, *fE=fO+nbin*nbin, *fEA=fE+nbin*nbin, *lambda=fEA+nbin*nbin; double h2=h*h, c2d, a,b,c,d, S[4], detS, G, x1,x2, alpha=0.5; /* histogram */ for(i=0; i1000) printf("\r"); /* smoothing and printing */ puts("\t\tSmoothing and printing."); for(j=0; j1) continue; for(i=0;i1) continue; fO[j*nbin+k] = fO[k*nbin+j] = (fO[j*nbin+k] + fO[k*nbin+j])/2; fE[j*nbin+k] = fE[k*nbin+j] = (fE[j*nbin+k] + fE[k*nbin+j])/2; fEA[j*nbin+k] = fEA[k*nbin+j] = (fEA[j*nbin+k] + fEA[k*nbin+j])/2; } } } for(j=0; j1) { printf("Great offer! I can smooth a few 2-D densities for free. How many do you want? "); scanf("%d", &nf2d); } */ if(nf2d>MAXNF2D) error2("I don't want to do that many!"); for(i=0; ix[j]) minx[j]=x[j]; else if (maxx[j]0))+nbin*nbin*3)*sizeof(double)); if(y==NULL) { printf("not enough mem for %d variables\n",n); exit(-1); } space=y+(1+(nf2d>0))*n; /* space[] points to y after the data */ for(jj=0;jj0?var[j*p+k]/t:-9)); } } fprintf(fout, "\n\n(C) Serial correlation coefficients within variables\n"); fprintf(fout, "\nlag variables ...\n\n"); for(i=0; i=r) return(0); r=fabs(f0); if(r=r) return(0); return (1); } int AlwaysCenter=0; double Small_Diff=1e-6; /* reasonable values 1e-5, 1e-7 */ int gradient (int n, double x[], double f0, double g[], double (*fun)(double x[],int n), double space[], int Central) { /* f0 = fun(x) is always given. */ int i,j; double *x0=space, *x1=space+n, eh0=Small_Diff, eh; /* 1e-7 */ if (Central) { for(i=0; i1) t=1; FOR (i,n) x[i] = x0[i] + t * p[i]; (*fun) (x, y, n, ny); for (i=0,s=0; ivmax) { istate=1; break; } } else { v*=smaller; xtoy (x, x0, n); s0=s; } } /* ii, maxround */ if (increase) *sx=s0; else { *sx=s; xtoy(x, x0, n); } if (ii == maxround) istate=-1; free (space); return (istate); } double bound (int nx, double x0[], double p[], double x[], int(*testx)(double x[], int nx)) { /* find largest t so that x[]=x0[]+t*p[] is still acceptable. for bounded minimization, p is possibly changed in this function using testx() */ int i, nd=0; double factor=20, by=1, small=1e-8; /* small=(SIZEp>1?1e-7:1e-8) */ xtoy (x0, x, nx); FOR (i,nx) { x[i]=x0[i]+small*p[i]; if ((*testx) (x, nx)) { p[i]=0.0; nd++; } x[i]=x0[i]; } if (nd==nx) { if (noisy) puts ("bound:no move.."); return (0); } for (by=0.75; ; ) { FOR (i,nx) x[i]=x0[i]+factor*p[i]; if ((*testx)(x,nx)==0) break; factor *= by; } return(factor); } double LineSearch (double(*fun)(double x),double *f,double *x0,double xb[2],double step, double e) { /* linear search using quadratic interpolation From Wolfe M. A. 1978. Numerical methods for unconstrained optimization: An introduction. Van Nostrand Reinhold Company, New York. pp. 62-73. step is used to find the bracket (a1,a2,a3) This is the same routine as LineSearch2(), but I have not got time to test and improve it properly. Ziheng note, September, 2002 */ int ii=0, maxround=100, i; double factor=2, step1, percentUse=0; double a0,a1,a2,a3,a4=-1,a5,a6, f0,f1,f2,f3,f4=-1,f5,f6; /* find a bracket (a1,a2,a3) with function values (f1,f2,f3) so that a1xb[1]) error2("err LineSearch: x0 out of range"); f2=f0=fun(a2=a0=*x0); step1=min2(step,(a0-xb[0])/4); step1=max2(step1,e); for(i=0,a1=a0,f1=f0; ; i++) { a1-=(step1*=factor); if(a1>xb[0]) { f1=fun(a1); if(f1>f2) break; else { a3=a2; f3=f2; a2=a1; f2=f1; } } else { a1=xb[0]; f1=fun(a1); if(f1<=f2) { a2=a1; f2=f1; } break; } if(noisy>2) printf("\ta = %.6f\tf = %.6f %5d\n", a2, f2, NFunCall); } if(i==0) { /* *x0 is the best point during the previous search */ step1=min2(step,(xb[1]-a0)/4); for(i=0,a3=a2,f3=f2; ; i++) { a3+=(step1*=factor); if(a3f2) break; else { a1=a2; f1=f2; a2=a3; f2=f3; } } else { a3=xb[1]; f3=fun(a3); if(f32) printf("\ta = %.6f\tf = %.6f %5d\n", a3, f3, NFunCall); } } /* iteration by quadratic interpolation, fig 2.2.9-10 (pp 71-71) */ for (ii=0; iia2+1e-99 || a3f1+1e-99 || f2>f3+1e-99) /* for linux */ { printf("\npoints out of order (ii=%d)!",ii+1); break; } a4 = (a2-a3)*f1+(a3-a1)*f2+(a1-a2)*f3; if (fabs(a4)>1e-100) a4=((a2*a2-a3*a3)*f1+(a3*a3-a1*a1)*f2+(a1*a1-a2*a2)*f3)/(2*a4); if (a4>a3 || a42) printf("\ta = %.6f\tf = %.6f %5d\n", a4, f4, NFunCall); if (fabs(f2-f4)*(1+fabs(f2))<=e && fabs(a2-a4)*(1+fabs(a2))<=e) break; if (a1<=a4 && a4<=a2) { /* fig 2.2.10 */ if (fabs(a2-a4)>.2*fabs(a1-a2)) { if (f1>=f4 && f4<=f2) { a3=a2; a2=a4; f3=f2; f2=f4; } else { a1=a4; f1=f4; } } else { if (f4>f2) { a5=(a2+a3)/2; f5=fun(a5); if (f5>f2) { a1=a4; a3=a5; f1=f4; f3=f5; } else { a1=a2; a2=a5; f1=f2; f2=f5; } } else { a5=(a1+a4)/2; f5=fun(a5); if (f5>=f4 && f4<=f2) { a3=a2; a2=a4; a1=a5; f3=f2; f2=f4; f1=f5; } else { a6=(a1+a5)/2; f6=fun(a6); if (f6>f5) { a1=a6; a2=a5; a3=a4; f1=f6; f2=f5; f3=f4; } else { a2=a6; a3=a5; f2=f6; f3=f5; } } } } } else { /* fig 2.2.9 */ if (fabs(a2-a4)>.2*fabs(a2-a3)) { if (f2>=f4 && f4<=f3) { a1=a2; a2=a4; f1=f2; f2=f4; } else { a3=a4; f3=f4; } } else { if (f4>f2) { a5=(a1+a2)/2; f5=fun(a5); if (f5>f2) { a1=a5; a3=a4; f1=f5; f3=f4; } else { a3=a2; a2=a5; f3=f2; f2=f5; } } else { a5=(a3+a4)/2; f5=fun(a5); if (f2>=f4 && f4<=f5) { a1=a2; a2=a4; a3=a5; f1=f2; f2=f4; f3=f5; } else { a6=(a3+a5)/2; f6=fun(a6); if (f6>f5) { a1=a4; a2=a5; a3=a6; f1=f4; f2=f5; f3=f6; } else { a1=a5; a2=a6; f1=f5; f2=f6; } } } } } } /* for (ii) */ if (f2<=f4) { *f=f2; a4=a2; } else *f=f4; return (*x0=(a4+a2)/2); } double fun_LineSearch (double t, double (*fun)(double x[],int n), double x0[], double p[], double x[], int n); double fun_LineSearch (double t, double (*fun)(double x[],int n), double x0[], double p[], double x[], int n) { int i; FOR (i,n) x[i]=x0[i] + t*p[i]; return( (*fun)(x, n) ); } double LineSearch2 (double(*fun)(double x[],int n), double *f, double x0[], double p[], double step, double limit, double e, double space[], int n) { /* linear search using quadratic interpolation from x0[] in the direction of p[], x = x0 + a*p a ~(0,limit) returns (a). *f: f(x0) for input and f(x) for output x0[n] x[n] p[n] space[n] adapted from Wolfe M. A. 1978. Numerical methods for unconstrained optimization: An introduction. Van Nostrand Reinhold Company, New York. pp. 62-73. step is used to find the bracket and is increased or reduced as necessary, and is not terribly important. */ int ii=0, maxround=10, status, i, nsymb=0; double *x=space, factor=4, small=1e-10, smallgapa=0.2; double a0,a1,a2,a3,a4=-1,a5,a6, f0,f1,f2,f3,f4=-1,f5,f6; /* look for bracket (a1, a2, a3) with function values (f1, f2, f3) step length step given, and only in the direction a>=0 */ if (noisy>2) printf ("\n%3d h-m-p %7.4f %6.4f %8.4f ",Iround+1,step,limit,norm(p,n)); if (step<=0 || limit=limit) { if (noisy>2) printf ("\nh-m-p:%20.8e%20.8e%20.8e %12.6f\n",step,limit,norm(p,n),*f); return (0); } a0=a1=0; f1=f0=*f; a2=a0+step; f2=fun_LineSearch(a2, fun,x0,p,x,n); if (f2>f1) { /* reduce step length so the algorithm is decreasing */ for (; ;) { step/=factor; if (step2) { printf("-"); nsymb++; } } } else { /* step length is too small? */ for (; ;) { step*=factor; if (step>limit) step=limit; a3=a0+step; f3=fun_LineSearch(a3, fun,x0,p,x,n); if (f3>=f2) break; if(!PAML_RELEASE && noisy>2) { printf("+"); nsymb++; } a1=a2; f1=f2; a2=a3; f2=f3; if (step>=limit) { if(!PAML_RELEASE && noisy>2) for(; nsymb<5; nsymb++) printf(" "); if (noisy>2) printf(" %12.6f%3c %6.4f %5d", *f=f3, 'm', a3, NFunCall); *f=f3; return(a3); } } } /* iteration by quadratic interpolation, fig 2.2.9-10 (pp 71-71) */ for (ii=0; ii1e-100) a4 = ((a2*a2-a3*a3)*f1+(a3*a3-a1*a1)*f2+(a1*a1-a2*a2)*f3)/(2*a4); if (a4>a3 || a4smallgapa*(a2-a1)) || (a4>a2 && a4-a2>smallgapa*(a3-a2))) status='Y'; else status='C'; } f4 = fun_LineSearch(a4, fun,x0,p,x,n); if(!PAML_RELEASE && noisy>2) putchar(status); if (fabs(f2-f4)2) for(nsymb+=ii+1; nsymb<5; nsymb++) printf(" "); break; } /* possible multiple local optima during line search */ if(!PAML_RELEASE && noisy>2 && ((a4f1) || (a4>a2&&f4>f3))) { printf("\n\na %12.6f %12.6f %12.6f %12.6f", a1,a2,a3,a4); printf( "\nf %12.6f %12.6f %12.6f %12.6f\n", f1,f2,f3,f4); for(a5=a1; a5<=a3; a5+=(a3-a1)/20) { printf("\t%.6e ",a5); if(n<5) FOR(i,n) printf("\t%.6f",x0[i] + a5*p[i]); printf("\t%.6f\n", fun_LineSearch(a5, fun,x0,p,x,n)); } puts("Linesearch2 a4: multiple optima?"); } if (a4<=a2) { /* fig 2.2.10 */ if (a2-a4>smallgapa*(a2-a1)) { if (f4<=f2) { a3=a2; a2=a4; f3=f2; f2=f4; } else { a1=a4; f1=f4; } } else { if (f4>f2) { a5=(a2+a3)/2; f5=fun_LineSearch(a5, fun,x0,p,x,n); if (f5>f2) { a1=a4; a3=a5; f1=f4; f3=f5; } else { a1=a2; a2=a5; f1=f2; f2=f5; } } else { a5=(a1+a4)/2; f5=fun_LineSearch(a5, fun,x0,p,x,n); if (f5>=f4) { a3=a2; a2=a4; a1=a5; f3=f2; f2=f4; f1=f5; } else { a6=(a1+a5)/2; f6=fun_LineSearch(a6, fun,x0,p,x,n); if (f6>f5) { a1=a6; a2=a5; a3=a4; f1=f6; f2=f5; f3=f4; } else { a2=a6; a3=a5; f2=f6; f3=f5; } } } } } else { /* fig 2.2.9 */ if (a4-a2>smallgapa*(a3-a2)) { if (f2>=f4) { a1=a2; a2=a4; f1=f2; f2=f4; } else { a3=a4; f3=f4; } } else { if (f4>f2) { a5=(a1+a2)/2; f5=fun_LineSearch(a5, fun,x0,p,x,n); if (f5>f2) { a1=a5; a3=a4; f1=f5; f3=f4; } else { a3=a2; a2=a5; f3=f2; f2=f5; } } else { a5=(a3+a4)/2; f5=fun_LineSearch(a5, fun,x0,p,x,n); if (f5>=f4) { a1=a2; a2=a4; a3=a5; f1=f2; f2=f4; f3=f5; } else { a6=(a3+a5)/2; f6=fun_LineSearch(a6, fun,x0,p,x,n); if (f6>f5) { a1=a4; a2=a5; a3=a6; f1=f4; f2=f5; f3=f6; } else { a1=a5; a2=a6; f1=f5; f2=f6; } } } } } } if (f2>f0 && f4>f0) a4=0; if (f2<=f4) { *f=f2; a4=a2; } else *f=f4; if(noisy>2) printf(" %12.6f%3d %6.4f %5d", *f, ii, a4, NFunCall); return (a4); } #define Safeguard_Newton int Newton (FILE *fout, double *f, double (* fun)(double x[], int n), int (* ddfun) (double x[], double *fx, double dx[], double ddx[], int n), int (*testx) (double x[], int n), double x0[], double space[], double e, int n) { int i,j, maxround=500; double f0=1e40, small=1e-10, h, SIZEp, t, *H, *x, *g, *p, *tv; H=space, x=H+n*n; g=x+n; p=g+n, tv=p+n; printf ("\n\nIterating by Newton\tnp:%6d\nInitial:", n); FOR (i,n) printf ("%8.4f", x0[i]); FPN (F0); if (fout) fprintf (fout, "\n\nNewton\tnp:%6d\n", n); if (testx (x0, n)) error2("Newton..invalid initials."); FOR (Iround, maxround) { if (ddfun) (*ddfun) (x0, f, g, H, n); else { *f = (*fun)(x0, n); Hessian (n, x0, *f, g, H, fun, tv); } matinv(H, n, n, tv); FOR (i,n) for (j=0,p[i]=0; j4) { while (t>small) { FOR (i,n) x[i]=x0[i]+t*p[i]; if ((*f=fun(x,n)) < f0) break; else t/=2; } } if (t2) { printf ("\n%3d h:%7.4f %12.6f x", Iround+1, SIZEp, *f); FOR (i,n) printf ("%7.4f ", x0[i]); } if (fout) { fprintf (fout, "\n%3d h:%7.4f%12.6f x", Iround+1, SIZEp, *f); FOR (i,n) fprintf (fout, "%7.4f ", x0[i]); fflush(fout); } if ((h=norm(x0,n))=xb[i][1]) { x[i]=xb[i][1]; xmark[i]= 1; continue; } ix[nfree++]=i; } if(noisy>2 && nfree=3) printf("warning: ming2, %d paras at boundary.",n-nfree); } f0=*f=(*fun)(x,n); xtoy(x,x0,n); SIZEp=99; if (noisy>2) { printf ("\nIterating by ming2\nInitial: fx= %12.6f\nx=",f0); FOR(i,n) printf(" %8.5f",x[i]); FPN(F0); } if (dfun) (*dfun) (x0, &f0, g0, n); else gradientB (n, x0, f0, g0, fun, tv, xmark); identity (H,nfree); for(Iround=0; Iround1) e_branches=0.1; else if(SIZEp>.5) e_branches=0.05; else if(SIZEp>.1) e_branches=0.01; else if(SIZEp>.05) e_branches=0.005; else e_branches=max2(e,1e-6); */ for (i=0,am=maxstep; i0 && (xb[i][1]-x0[i])/p[i]2) printf(".. "); identity(H,nfree); fail=1; } } else { fail=0; FOR(i,n) x[i]=x0[i]+alpha*p[i]; w=min2(2,e*1000); if(e<1e-4 && e>1e-6) w=0.01; if(Iround==0 || SIZEp=Ngoodtimes) && SIZEp<(e>1e-5?1:.001) && H_end(x0,x,f0,*f,e,e,n)) break; } if (dfun) (*dfun) (x, f, g, n); else gradientB (n, x, *f, g, fun, tv, xmark); /* for(i=0; i0) xmark[i]=1; if (xmark[i]==0) continue; xtoy (H, C, nfree*nfree); FOR (it, nfree) if (ix[it]==i) break; for (i1=it; i1=it))*(nfree+1) + i2+(i2>=it)]; } for (i=0,it=0,w=0; iw) { it=i; w=-g[i]; } else if (xmark[i]==1 && -g[i]<-w) { it=i; w=g[i]; } } if (w>10*SIZEp/nfree) { /* *** */ xtoy (H, C, nfree*nfree); FOR (i1,nfree) FOR (i2,nfree) H[i1*(nfree+1)+i2]=C[i1*nfree+i2]; FOR (i1,nfree+1) H[i1*(nfree+1)+nfree]=H[nfree*(nfree+1)+i1]=0; H[(nfree+1)*(nfree+1)-1]=1; xmark[it]=0; ix[nfree++]=it; } if (noisy>2) { printf (" | %d/%d", n-nfree, n); /* FOR (i,n) if (xmark[i]) printf ("%4d", i+1); */ } for (i=0,f0=*f; i2) FPN(F0); if (Iround==maxround) { if (fout) fprintf (fout,"\ncheck convergence!\n"); return(-1); } if(nfree==n) { xtoy(H, space, n*n); /* H has variance matrix, or inverse of Hessian */ return(1); } return(0); } int ming1 (FILE *fout, double *f, double (* fun)(double x[], int n), int (*dfun) (double x[], double *f, double dx[], int n), int (*testx) (double x[], int n), double x0[], double space[], double e, int n) { /* n-D minimization using quasi-Newton or conjugate gradient algorithms, using function and its gradient. g0[n] g[n] p[n] x[n] y[n] s[n] z[n] H[n*n] tv[2*n] using bound() */ int i,j, maxround=1000, fail=0; double small=1.e-20; /* small value for checking |w|=0 */ double f0, *g0, *g, *p, *x, *y, *s, *z, *H, *tv; double w,v, t, h; if (testx (x0, n)) { printf ("\nInvalid initials..\n"); matout(F0,x0,1,n); return(-1); } f0 = *f = (*fun)(x0, n); if (noisy>2) { printf ("\n\nIterating by ming1\nInitial: fx= %12.6f\nx=", f0); FOR (i,n) printf ("%8.4f", x0[i]); FPN (F0); } if (fout) { fprintf (fout, "\n\nIterating by ming1\nInitial: fx= %12.6f\nx=", f0); FOR (i,n) fprintf (fout, "%10.6f", x0[i]); } g0=space; g=g0+n; p=g+n; x=p+n; y=x+n; s=y+n; z=s+n; H=z+n; tv=H+n*n; if (dfun) (*dfun) (x0, &f0, g0, n); else gradient (n, x0, f0, g0, fun, tv, AlwaysCenter); SIZEp=0; xtoy (x0, x, n); xtoy (g0,g,n); identity (H,n); FOR (Iround, maxround) { FOR (i,n) for (j=0,p[i]=0.; j1e32) { if (fail) { if(SIZEp>.1 && noisy>2) printf("\nSIZEp:%9.4f Iround:%5d", SIZEp, Iround+1); if (AlwaysCenter) { Iround=maxround; break; } else { AlwaysCenter=1; identity(H,n); fail=1; } } else { identity(H, n); fail=1; } } else { fail=0; FOR(i,n) x[i]=x0[i]+t*p[i]; if (fout) { fprintf (fout, "\n%3d %7.4f%14.6f x", Iround+1, SIZEp, *f); FOR (i,n) fprintf (fout, "%8.5f ", x[i]); fflush (fout); } if (SIZEp<0.001 && H_end (x0,x,f0,*f,e,e,n)) { xtoy(x,x0,n); break; } } if (dfun) (*dfun) (x, f, g, n); else gradient (n,x,*f,g,fun,tv, (AlwaysCenter||fail||SIZEp<0.01)); for (i=0,f0=*f; i