← Back to branch summary

~ubuntu-branches/ubuntu/vivid/paml/vivid

« back to all changes in this revision

Viewing changes to src/mcmctree.c

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2013-08-08 22:07:27 UTC
  • mfrom: (1.1.3)
  • Revision ID: package-import@ubuntu.com-20130808220727-deh08eeugdwjff1h
Tags: 4.7-1
* New upstream version
* debian/control:
   - cme fix dpkg-control
   - use anonscm in Vcs fields
* debian/rules: provide get-orig-source target to make use of new uscan
  to remove windows binaries

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/mcmctree.c
4
4
   
5
5
                   Ziheng YANG, since December 2002
6
6
 
7
 
   cc -o mcmctree -m64 -march=opteron -mtune=opteron -ansi -funsigned-char -O3 -funroll-loops -fomit-frame-pointer -finline-functions mcmctree.c tools.c -lm
8
 
   cc -o mcmctree -march=pentiumpro -mcpu=pentiumpro -O4 -funroll-loops -fomit-frame-pointer -finline-functions mcmctree.c tools.c -lm
9
 
   cc -o mcmctree -march=athlon -mcpu=athlon -O4 -funroll-loops -fomit-frame-pointer -finline-functions mcmctree.c tools.c -lm
10
 
 
 
7
   cc -o mcmctree -m64 -ansi -funsigned-char -O3 -funroll-loops -fomit-frame-pointer -finline-functions mcmctree.c tools.c -lm
 
8
   cc -o mcmctree -O4 -funroll-loops -fomit-frame-pointer -finline-functions mcmctree.c tools.c -lm
 
9
   
11
10
   cc -o infinitesites -DINFINITESITES -march=athlon -mcpu=athlon -Wall -O2 -funroll-loops -fomit-frame-pointer -finline-functions mcmctree.c tools.c -lm
12
11
 
13
12
   cl -O2 -W3 mcmctree.c tools.c
14
13
 
 
14
   cl -FeInfiniteSites.exe -DINFINITESITES -W3 -D_CRT_SECURE_NO_DEPRECATE -O2 mcmctree.c tools.c
 
15
 
15
16
   mcmctree <ControlFileName>
 
17
 
 
18
   InfiniteSites <ControlFileName>
 
19
   FixedDsClock1.txt or FixedDsClock23.txt are hard-coded file names
16
20
*/
17
21
 
18
22
/*
60
64
double lnptC(void);
61
65
double lnptCalibrationDensity(double t, int fossil, double p[7]);
62
66
int SetupPriorTimesFossilErrors(void);
63
 
double lnpriorTimesBDS_Approach2(void);
 
67
double lnpriorTimesBDS_Approach1(void);
64
68
double lnpriorTimesTipDate(void);
65
69
double lnpriorTimes(void);
66
70
double lnpriorRates(void);
67
71
double logPriorRatioGamma(double xnew, double xold, double a, double b);
68
72
void copySptree(void);
69
73
void printSptree(void);
70
 
double Infinitesites (void);
71
 
double InfinitesitesClock (double *FixedDs);
 
74
double Infinitesites(void);
 
75
double InfinitesitesClock(double *FixedDs);
72
76
int collectx (FILE* fout, double x[]);
73
77
int MCMC(FILE* fout);
74
78
int LabelOldCondP(int spnode);
81
85
double mixing(double *lnL, double finetune);
82
86
double UpdatePFossilErrors(double finetune);
83
87
int getPfossilerr (double postEFossil[], double nround);
84
 
int DescriptiveStatisticsSimpleMCMCTREE (FILE *fout, char infile[], int nbin, int nrho);
 
88
int DescriptiveStatisticsSimpleMCMCTREE (FILE *fout, char infile[], int nbin);
85
89
 
86
90
struct CommonInfo {
87
 
   char *z[NS], *spname[NS], seqf[256],outf[256],treef[256],daafile[96],mcmcf[96],inBVf[96];
 
91
   unsigned char *z[NS];
 
92
   char *spname[NS], seqf[256],outf[256],treef[256],daafile[96],mcmcf[96],inBVf[96];
88
93
   char oldconP[NNODE];       /* update conP for node? (0 yes; 1 no) */
89
94
   int seqtype, ns, ls, ngene, posG[2],lgene[1], *pose, npatt, readpattern;
90
95
   int np, ncode, ntime, nrate, nrgene, nalpha, npi, ncatG, print;
125
130
 
126
131
struct SPECIESTREE {
127
132
   int nbranch, nnode, root, nspecies, nfossil;
128
 
   double RootAge[4], Pfossilerr, *CcomFossilErr;
 
133
   double RootAge[4];
129
134
   struct TREESPN {
130
135
      char name[LSPNAME+1], fossil, usefossil;  /* fossil: 0, 1(L), 2(U), 3(B), 4(G) */
131
136
      int father, nson, sons[2];
146
151
   double rgene[NGENE], kappa[NGENE], omega[NGENE], alpha[NGENE];
147
152
   double BDS[4];  /* parameters in the birth-death-sampling model */
148
153
   double rgenegamma[2], kappagamma[2], omegagamma[2], alphagamma[2], sigma2gamma[2];
149
 
   double pfossilerror[3]; /* (p_beta, q_beta, NminCorrect) */
 
154
   double pfossilerror[3], /* (p_beta, q_beta, NminCorrect) */ Pfossilerr, *CcomFossilErr;
150
155
   double sigma2[NGENE];  /* sigma2[g] are the variances */
151
156
   double *blMLE[NGENE], *Gradient[NGENE], *Hessian[NGENE];
152
157
   int transform;
164
169
int nR=4;
165
170
double PMat[16], Cijk[64], Root[4];
166
171
double _rateSite=1, OldAge=999;
167
 
int LASTROUND=0, BayesEB, debug=0, testlnL=0, NPMat=0; /* no use for this */
 
172
int debug=0, LASTROUND=0, BayesEB, testlnL=0, NPMat=0; /* no use for this */
168
173
 
169
174
/* for sptree.nodes[].fossil: lower, upper, bounds, gamma, inverse-gamma */
170
175
enum {LOWER_F=1, UPPER_F, BOUND_F, GAMMA_F, SKEWN_F, SKEWT_F, S2N_F} FOSSIL_FLAGS;
177
182
#define MCMCTREE  1
178
183
#include "treesub.c"
179
184
 
 
185
 
180
186
int main (int argc, char *argv[])
181
187
{
182
 
   char ctlf[96] = "mcmctree.ctl";
 
188
   char ctlf[128] = "mcmctree.ctl";
183
189
   int i,j,k;
184
190
   FILE  *fout;
185
191
 
191
197
   com.ncode=4;      com.clock=1;
192
198
 
193
199
   printf("MCMCTREE in %s\n", pamlVerStr);
194
 
   if (argc>1) 
195
 
      { strncpy(ctlf, argv[1], 95); printf ("\nctlfile reset to %s.\n", ctlf); }
 
200
   if(argc>1)
 
201
      strncpy(ctlf, argv[1], 127);
196
202
 
197
203
   data.BDS[0]=1;    data.BDS[1]=1;  data.BDS[2]=0; 
198
 
   com.cleandata=0; mcmc.usedata=2;
199
204
   strcpy(com.outf, "out");
200
205
   strcpy(com.mcmcf, "mcmc.out");
201
206
 
206
211
   fprintf(fout, "MCMCTREE (%s) %s\n", pamlVerStr, com.seqf);
207
212
 
208
213
   ReadTreeSeqs(fout);
 
214
   if(data.pfossilerror && (data.pfossilerror[2]<0 || data.pfossilerror[2]>sptree.nfossil))
 
215
      error2("nMinCorrect for fossil errors is out of range.");
 
216
 
209
217
   if(mcmc.usedata==1) {
210
218
      if(com.seqtype!=0) error2("usedata = 1 for nucleotides only");
211
219
      if(com.alpha==0)
224
232
      if(com.seqtype==2) com.ncode = 20;
225
233
   }
226
234
   else if(mcmc.usedata==3) {
227
 
      GenerateBlengthGH("out.BV");
 
235
      GenerateBlengthGH("out.BV");  /* this is used so that the in.BV is not overwritten */
228
236
      exit(0);
229
237
   }
230
238
 
239
247
      for(i=0; i<4; i++) 
240
248
         sptree.nodes[sptree.root].pfossil[i] = sptree.RootAge[i];
241
249
   }
242
 
   printf("\nFossil calibration information used.\n");
 
250
   if( com.TipDate_TimeUnit==0) 
 
251
      printf("\nFossil calibration information used.\n");
243
252
   for(i=sptree.nspecies; i<sptree.nspecies*2-1; i++) {
244
253
      if((k=sptree.nodes[i].fossil) == 0) continue;
245
254
      printf("Node %3d: %3s ( ", i+1, fossils[k]);
409
418
 
410
419
      xtoy(data.pi[locus], com.pi, com.ncode);
411
420
      if(com.model<=TN93)
412
 
         EigenTN93(com.model, com.kappa, com.kappa, com.pi, &nR, Root, Cijk);
 
421
         eigenTN93(com.model, com.kappa, com.kappa, com.pi, &nR, Root, Cijk);
413
422
 
414
423
      if(com.alpha)
415
424
         DiscreteGamma (com.freqK,com.rK,com.alpha,com.alpha,com.ncatG,DGammaUseMedian);
591
600
   }
592
601
   else {              /* independent & correlated rates */
593
602
      for(i=0; i<tree.nnode; i++) {
594
 
         if(i==tree.root) continue;
 
603
         if(i==tree.root) 
 
604
            continue;
595
605
         for(j=nodes[i].ipop,b=0; j!=nodes[nodes[i].father].ipop; j=dad) {
596
606
            dad = sptree.nodes[j].father;
597
 
            t = sptree.nodes[dad].age-sptree.nodes[j].age;
 
607
            t = sptree.nodes[dad].age - sptree.nodes[j].age;
598
608
            b += t * sptree.nodes[j].rates[locus];
599
609
         }
600
610
         nodes[i].branch = b;
649
659
*/
650
660
   int debug=0, i,j, ib, nb=tree.nnode-1-1;  /* don't use tree.nbranch, which is not up to date. */
651
661
   int son1=nodes[tree.root].sons[0], son2=nodes[tree.root].sons[1];
652
 
   double lnL=0, z[NS*2-1], *g=data.Gradient[locus], *H=data.Hessian[locus], JCc=(com.ncode-1.0)/com.ncode;
 
662
   double lnL=0, z[NS*2-1], *g=data.Gradient[locus], *H=data.Hessian[locus], cJC=(com.ncode-1.0)/com.ncode;
653
663
   double bTlog=1e-5, elog=0.1, e;   /* change the same line in ReadBlengthGH() as well  */
654
664
 
655
665
   /* construct branch lengths */
669
679
         if(data.transform==SQRT_B)       /* sqrt transform, MLE transformed */
670
680
            z[ib] = sqrt(z[ib]);
671
681
         else if(data.transform==ARCSIN_B)    /* JC transform, MLE transformed */
672
 
            z[ib] = 2*asin(sqrt( JCc - JCc*exp(-z[ib]/JCc) ));
 
682
            z[ib] = 2*asin(sqrt( cJC - cJC*exp(-z[ib]/cJC) ));
673
683
         z[ib] -= data.blMLE[locus][ib];
674
684
      }
675
685
   }
709
719
   char line[100];
710
720
   int locus, i, j, nb, son1, son2, leftsingle;
711
721
   double small = 1e-20;
712
 
   double dbu[NBRANCH], dbu2[NBRANCH], u, JCc=(com.ncode-1.0)/com.ncode, sin2u, cos2u;
 
722
   double dbu[NBRANCH], dbu2[NBRANCH], u, cJC=(com.ncode-1.0)/com.ncode, sin2u, cos2u;
713
723
   double bTlog=1e-5, elog=0.1, e;   /* change the line in lnpD_locus_Approx() as well */
714
724
   int debug = 0;
715
725
 
801
811
            if(data.blMLE[locus][i] < 1e-20) 
802
812
               error2("blength should be > 0 for JC transform");
803
813
            */
804
 
            u = 2*asin(sqrt(JCc - JCc*exp(- data.blMLE[locus][i]/JCc )));
 
814
            u = 2*asin(sqrt(cJC - cJC*exp(- data.blMLE[locus][i]/cJC )));
805
815
            sin2u = sin(u/2);  cos2u = cos(u/2);
806
 
            dbu[i]  = sin2u*cos2u/(1-sin2u*sin2u/JCc);
807
 
            dbu2[i] = (cos2u*cos2u-sin2u*sin2u)/2/(1-sin2u*sin2u/JCc) + dbu[i]*dbu[i]/JCc;
 
816
            dbu[i]  = sin2u*cos2u/(1-sin2u*sin2u/cJC);
 
817
            dbu2[i] = (cos2u*cos2u-sin2u*sin2u)/2/(1-sin2u*sin2u/cJC) + dbu[i]*dbu[i]/cJC;
808
818
         }
809
819
      }
810
820
 
827
837
      }
828
838
      else if(data.transform==ARCSIN_B)    /* JC transform on branch lengths */
829
839
         for(i=0; i<nb; i++)
830
 
            data.blMLE[locus][i] = 2*asin(sqrt(JCc - JCc*exp(-data.blMLE[locus][i]/JCc )));
 
840
            data.blMLE[locus][i] = 2*asin(sqrt(cJC - cJC*exp(-data.blMLE[locus][i]/cJC )));
831
841
   }    /* for(locus) */
832
842
 
833
843
   fclose(fBGH);
896
906
      if(com.alpha) 
897
907
         fprintf(fctl, "fix_alpha = 0\nalpha = 0.5\nncatG = %d\n", com.ncatG);
898
908
      fprintf(fctl, "Small_Diff = 0.1e-6\ngetSE = 2\n");
899
 
      fprintf(fctl, "method = %d\n", (com.alpha||com.ns<10 ? 0 : 1));
 
909
      fprintf(fctl, "method = %d\n", (com.alpha==0||com.ns>100 ? 1 : 0));
900
910
 
901
911
      fclose(fseq);  fclose(ftree);  fclose(fctl);
902
912
 
914
924
 
915
925
int GetOptions (char *ctlf)
916
926
{
917
 
   int iopt,i, nopt=28, lline=4096;
 
927
   int iopt,i, nopt=29, lline=4096;
918
928
   char line[4096],*pline, opt[32], *comment="*#";
919
929
   char *optstr[] = {"seed", "seqfile","treefile", "outfile", "mcmcfile", 
920
 
        "seqtype", "aaRatefile", "icode", "usedata", "ndata", "model", "clock", 
 
930
        "seqtype", "aaRatefile", "icode", "noisy", "usedata", "ndata", "model", "clock", 
921
931
        "TipDate", "RootAge", "fossilerror", "alpha", "ncatG", "cleandata", 
922
932
        "BDparas", "kappa_gamma", "alpha_gamma", 
923
933
        "rgene_gamma", "sigma2_gamma", "print", "burnin", "sampfreq", 
926
936
   FILE  *fctl=gfopen (ctlf, "r");
927
937
 
928
938
   strcpy(com.inBVf, "in.BV");
929
 
   data.transform = ARCSIN_B;
 
939
   data.transform = ARCSIN_B;  /* SQRT_B, LOG_B, ARCSIN_B */
930
940
   if (fctl) {
931
941
      if (noisy) printf ("\nReading options from %s..\n", ctlf);
932
942
      for (;;) {
951
961
                  case ( 5): com.seqtype=(int)t;    break;
952
962
                  case ( 6): sscanf(pline+2,"%s", com.daafile);  break;
953
963
                  case ( 7): com.icode=(int)t;      break;
954
 
                  case ( 8): sscanf(pline+1, "%d %s%d", &mcmc.usedata, com.inBVf, &data.transform);
 
964
                  case ( 8): noisy=(int)t;          break;
 
965
                  case ( 9): sscanf(pline+1, "%d %s%d", &mcmc.usedata, com.inBVf, &data.transform);
955
966
                     if(mcmc.usedata==2 && strchr(com.inBVf, '*'))
956
967
                        strcpy(com.inBVf, "in.BV");
957
968
                     break;
958
 
                  case ( 9): com.ndata=(int)t;      break;
959
 
                  case (10): com.model=(int)t;      break;
960
 
                  case (11): com.clock=(int)t;      break;
961
 
                  case (12): 
 
969
                  case (10): com.ndata=(int)t;      break;
 
970
                  case (11): com.model=(int)t;      break;
 
971
                  case (12): com.clock=(int)t;      break;
 
972
                  case (13): 
962
973
                     sscanf(pline+2,"%lf%lf", &com.TipDate, &com.TipDate_TimeUnit);
963
974
                     if(com.TipDate && com.TipDate_TimeUnit==0) error2("should set com.TipDate_TimeUnit");
 
975
                     data.transform = SQRT_B;  /* SQRT_B, LOG_B, ARCSIN_B */
964
976
                     break;
965
 
                  case (13):
 
977
                  case (14):
966
978
                     sptree.RootAge[2] = sptree.RootAge[3] = 0.025;  /* default tail probs */
967
979
                     if((strchr(line, '>') || strchr(line, '<')) && (strstr(line, "U(") || strstr(line, "B(")))
968
980
                        error2("don't mix < U B on the RootAge line");
977
989
                     else if((pline=strstr(line, "B(")))
978
990
                        sscanf(pline+2, "%lf,%lf,%lf,%lf", &sptree.RootAge[0], &sptree.RootAge[1], &sptree.RootAge[2], &sptree.RootAge[3]);
979
991
                     break;
980
 
                  case (14):
 
992
                  case (15):
981
993
                     data.pfossilerror[0] = 0.0;
982
 
                     data.pfossilerror[2] = 2;  /* default: minimum 2 good fossils */
 
994
                     data.pfossilerror[2] = 1;  /* default: minimum 2 good fossils */
983
995
                     sscanf(pline+1, "%lf%lf%lf", data.pfossilerror, data.pfossilerror+1, data.pfossilerror+2);
984
996
                     break;
985
 
                  case (15): com.alpha=t;           break;
986
 
                  case (16): com.ncatG=(int)t;      break;
987
 
                  case (17): com.cleandata=(int)t;  break;
988
 
                  case (18): 
 
997
                  case (16): com.alpha=t;           break;
 
998
                  case (17): com.ncatG=(int)t;      break;
 
999
                  case (18): com.cleandata=(int)t;  break;
 
1000
                  case (19): 
989
1001
                     sscanf(pline+1,"%lf%lf%lf%lf", &data.BDS[0],&data.BDS[1],&data.BDS[2],&data.BDS[3]);
990
1002
                     break;
991
 
                  case (19): 
 
1003
                  case (20): 
992
1004
                     sscanf(pline+1,"%lf%lf", data.kappagamma, data.kappagamma+1); break;
993
 
                  case (20): 
 
1005
                  case (21): 
994
1006
                     sscanf(pline+1,"%lf%lf", data.alphagamma, data.alphagamma+1); break;
995
 
                  case (21): 
 
1007
                  case (22): 
996
1008
                     sscanf(pline+1,"%lf%lf", data.rgenegamma, data.rgenegamma+1); break;
997
 
                  case (22): 
 
1009
                  case (23): 
998
1010
                     sscanf(pline+1,"%lf%lf", data.sigma2gamma, data.sigma2gamma+1); break;
999
 
                  case (23): mcmc.print=(int)t;     break;
1000
 
                  case (24): mcmc.burnin=(int)t;    break;
1001
 
                  case (25): mcmc.sampfreq=(int)t;  break;
1002
 
                  case (26): mcmc.nsample=(int)t;   break;
1003
 
                  case (27):
 
1011
                  case (24): mcmc.print=(int)t;     break;
 
1012
                  case (25): mcmc.burnin=(int)t;    break;
 
1013
                  case (26): mcmc.sampfreq=(int)t;  break;
 
1014
                  case (27): mcmc.nsample=(int)t;   break;
 
1015
                  case (28):
1004
1016
                     sscanf(pline+1,"%d:%lf%lf%lf%lf%lf%lf", &mcmc.resetFinetune, eps,eps+1,eps+2,eps+3,eps+4,eps+5);
1005
1017
                     break;
1006
1018
               }
1066
1078
   int i,j, ir, nround=0, nsaved=0;
1067
1079
   int s=sptree.nspecies, g=data.ngene;
1068
1080
   double t, tnew, naccept=0;
1069
 
   double e=mcmc.finetune[0], lnp, lnpnew, lnacceptance, c, *x;
1070
 
   double tmean, t025,t975;
 
1081
   double e=mcmc.finetune[0], lnp, lnpnew, lnacceptance, c, *x, Pjump=0;
 
1082
   double tmean, t025,t975, tL, tU;
1071
1083
   char timestr[32];
1072
1084
 
1073
1085
   matout2(F0, FixedDs, 1, s-1+g-1, 8, 4);
1078
1090
   if(x==NULL) error2("oom x");
1079
1091
 
1080
1092
   for(ir=-mcmc.burnin,tmean=0; ir<mcmc.nsample*mcmc.sampfreq; ir++) {
1081
 
      if(ir==0) { nround=0; naccept=0; tmean=0; }
1082
 
      lnacceptance = e*rnd2NormalSym(m2NormalSym);
 
1093
      if(ir==0 || (ir<0 && ir%(mcmc.burnin/2)==0)) {
 
1094
         nround=0; naccept=0; tmean=0; 
 
1095
         ResetFinetuneSteps(NULL, &Pjump, &e, 1);
 
1096
      }
 
1097
      lnacceptance = e*rndBactrian(mBactrian);
1083
1098
      c = exp(lnacceptance);
1084
1099
      tnew = t*c;
1085
1100
      lnpnew = lnpInfinitesitesClock(tnew, FixedDs);
1094
1109
      if(ir>=0 && (ir+1)%mcmc.sampfreq==0)
1095
1110
         x[nsaved++]=t;
1096
1111
 
 
1112
      Pjump = naccept/nround;
1097
1113
      if((ir+1)%max2(mcmc.sampfreq, mcmc.sampfreq*mcmc.nsample/10000)==0)
1098
 
         printf("\r%3.0f%%  %7.2f mean t0 = %9.6f", (ir+1.)/(mcmc.nsample*mcmc.sampfreq)*100, naccept/nround,tmean);
1099
 
      if(mcmc.sampfreq*mcmc.nsample>20 && (ir+1)%(mcmc.sampfreq*mcmc.nsample/20)==0)
 
1114
         printf("\r%3.0f%%  %7.2f mean t0 = %9.6f", (ir+1.)/(mcmc.nsample*mcmc.sampfreq)*100, Pjump,tmean);
 
1115
      if(mcmc.sampfreq*mcmc.nsample>20 && (ir+1)%(mcmc.sampfreq*mcmc.nsample/10)==0)
1100
1116
         printf(" %s\n", printtime(timestr));
1101
1117
   }
1102
1118
 
1103
1119
   qsort(x, (size_t)mcmc.nsample, sizeof(double), comparedouble);
1104
 
   t025=x[(int)(mcmc.nsample*.025+.5)]; t975=x[(int)(mcmc.nsample*.975+.5)];
 
1120
   t025 = x[(int)(mcmc.nsample*.025+.5)];
 
1121
   t975 = x[(int)(mcmc.nsample*.975+.5)];
1105
1122
 
1106
1123
   /* Below x[] is used to collect the posterior means and 95% CIs */
1107
1124
   for(i=0; i<3; i++) {
1110
1127
      for(j=s; j<sptree.nnode; j++)
1111
1128
         x[i*(s+g)+(j-s)] = sptree.nodes[j].age;
1112
1129
      for(j=0; j<g; j++) 
1113
 
         x[i*(s+g)+s-1+j]=data.rgene[j];
1114
 
   }
1115
 
   printf("\nmean & 95%% CI for times\n\n");
1116
 
   for(j=s; j<sptree.nnode; j++)
1117
 
      printf("Node %2d: %9.5f (%9.5f, %9.5f)\n", j+1,x[j-s],x[(s+g)+j-s],x[2*(s+g)+j-s]);
 
1130
         x[i*(s+g)+s-1+j] = data.rgene[j];
 
1131
   }
 
1132
   printf("\nmean (95%% CI) CI-width for times\n\n");
 
1133
   for(j=s; j<sptree.nnode; j++) {
 
1134
      tL = x[(s+g)+j-s]; 
 
1135
      tU = x[2*(s+g)+j-s];
 
1136
      printf("Node %2d: %9.6f (%9.6f, %9.6f) %9.6f\n", j+1, x[j-s], tL, tU, tU-tL);
 
1137
   }
1118
1138
   printf("\nmean & 95%% CI for rates\n\n");
1119
1139
   for(j=0; j<g; j++)
1120
 
      printf("gene %2d: %9.5f (%9.5f, %9.5f)\n", j+1,x[s-1+j], x[2*(s+g)+s-1+j], x[(s+g)+s-1+j]);
 
1140
      printf("gene %2d: %9.6f (%9.6f, %9.6f)\n", j+1,x[s-1+j], x[2*(s+g)+s-1+j], x[(s+g)+s-1+j]);
1121
1141
 
1122
1142
   printf("\nNote: the posterior has only one dimension.\n");
1123
1143
   free(x);
1199
1219
   int nxpr[2]={12,3};
1200
1220
   int MixingStep=1;
1201
1221
   char line[10000], timestr[36];
1202
 
   double *e=mcmc.finetune, y, ynew, yb[2], naccept[4]={0};
 
1222
   double *e=mcmc.finetune, y, ynew, yb[2], naccept[4]={0}, Pjump[4]={0};
1203
1223
   double lnL, lnLnew, lnacceptance, c,lnc;
1204
1224
   double *x,*mx, *FixedDs,*b, maxt0,tson[2];
1205
1225
   char *FidedDf[2]={"FixedDsClock1.txt", "FixedDsClock23.txt"};
1211
1231
   GetMem();
1212
1232
   GetInitials();
1213
1233
 
1214
 
   printf("\nInfinite sites\n");
1215
 
   printf("Fixed branch lengths from %s (s = %d  g = %d):\n\n", FidedDf[com.clock>1], s,g);
 
1234
   printf("\nInfiniteSites, reading fixed distance data from %s (s = %d  g = %d):\n\n", FidedDf[com.clock>1], s,g);
1216
1235
   np = s-1 + g + g + g;   /* times, rates, mu & sigma2 */
1217
1236
   FixedDs = (double*)malloc((g*sptree.nnode+np*2)*sizeof(double));
1218
1237
   if(FixedDs==NULL) error2("oom");
1220
1239
   mx = x+np;
1221
1240
   fscanf(fdin, "%d", &i);
1222
1241
   if(i!=s) error2("wrong number of species in FixedDs.txt");
1223
 
   fgets(line, lline, fdin);
1224
1242
 
1225
1243
   if(data.pfossilerror[0]) {
1226
1244
      puts("model of fossil errors for infinite data not tested yet.");
1229
1247
   }
1230
1248
 
1231
1249
   if(com.clock==1) { /* global clock: read FixedDs[] and close file. */
1232
 
      for(i=0; i<s-1+g-1; i++) fscanf(fdin, "%lf", &FixedDs[i]);
 
1250
      for(i=0; i<s-1+g-1; i++) 
 
1251
         fscanf(fdin, "%lf", &FixedDs[i]);
1233
1252
      fclose(fdin);
1234
1253
      InfinitesitesClock(FixedDs);
1235
 
     return(0);
 
1254
      return(0);
1236
1255
   }
1237
1256
  
1238
1257
   for(locus=0,b=FixedDs,nodes=nodes_t; locus<g; locus++,b+=sptree.nnode) {
1276
1295
   printf("\nUsing finetune parameters from the control file\n");
1277
1296
 
1278
1297
   for(ir=-mcmc.burnin,nround=0; ir<mcmc.sampfreq*mcmc.nsample; ir++) {
1279
 
      if(ir==0) { /* reset after burnin */
 
1298
      if(ir==0 || (ir<0 && ir%(mcmc.burnin/2)==0)) {
1280
1299
         nround=0; naccept[0]=naccept[1]=naccept[2]=naccept[3]=0;  zero(mx,com.np);
 
1300
         ResetFinetuneSteps(NULL, Pjump, mcmc.finetune, 4);
1281
1301
      }
1282
1302
      for(ip=0; ip<np; ip++) {
1283
1303
         lnacceptance = 0;
1298
1318
                  yb[0] = max2(yb[0], sptree.nodes[root].age-maxt0);
1299
1319
               }
1300
1320
            }
1301
 
            ynew = y + e[0]*rnd2NormalSym(m2NormalSym);
 
1321
            ynew = y + e[0]*rndBactrian(mBactrian);
1302
1322
            ynew = sptree.nodes[s+ip].age = reflect(ynew,yb[0],yb[1]);
1303
1323
         }
1304
1324
         else if(ip-(s-1)<g) {    /* rate r0 for root son0 for loci (r0*t0<b0) */
1308
1328
            yb[0] = 0;
1309
1329
            yb[1] = FixedDs[(ip-(s-1))*sptree.nnode+sons[0]]/y;
1310
1330
            y = sptree.nodes[sons[0]].rates[ip-(s-1)];
1311
 
            ynew = y + e[1]*rnd2NormalSym(m2NormalSym);
 
1331
            ynew = y + e[1]*rndBactrian(mBactrian);
1312
1332
            ynew = sptree.nodes[sons[0]].rates[ip-(s-1)] = reflect(ynew,yb[0],yb[1]);
1313
1333
         }
1314
1334
         else if (ip-(s-1+g)<g) { /* mu for loci */
1315
 
            lnacceptance = e[3]*rnd2NormalSym(m2NormalSym);
 
1335
            lnacceptance = e[3]*rndBactrian(mBactrian);
1316
1336
            c=exp(lnacceptance);
1317
1337
            y = data.rgene[ip-(s-1+g)];
1318
1338
            ynew = data.rgene[ip-(s-1+g)] *= c;
1319
 
            lnacceptance+=logPriorRatioGamma(ynew, y, data.rgenegamma[0], data.rgenegamma[1]);
 
1339
            lnacceptance += logPriorRatioGamma(ynew, y, data.rgenegamma[0], data.rgenegamma[1]);
1320
1340
         }
1321
1341
         else {                   /* sigma2 for loci */
1322
 
            lnacceptance = e[3]*rnd2NormalSym(m2NormalSym);
 
1342
            lnacceptance = e[3]*rndBactrian(mBactrian);
1323
1343
            c=exp(lnacceptance);
1324
1344
            y = data.sigma2[ip-(s-1+g*2)];
1325
1345
            ynew = data.sigma2[ip-(s-1+g*2)] *= c;
1326
 
            lnacceptance+=logPriorRatioGamma(ynew, y, data.sigma2gamma[0], data.sigma2gamma[1]);
 
1346
            lnacceptance += logPriorRatioGamma(ynew, y, data.sigma2gamma[0], data.sigma2gamma[1]);
1327
1347
         }
1328
1348
 
1329
1349
         lnLnew = lnpInfinitesites(FixedDs);
1348
1368
      }
1349
1369
 
1350
1370
      if(MixingStep) {  /* this multiples times by c and divides r and mu by c. */
1351
 
         lnc = e[2]*rnd2NormalSym(m2NormalSym);
 
1371
         lnc = e[2]*rndBactrian(mBactrian);
1352
1372
         c = exp(lnc);
1353
1373
         lnacceptance = (s - 1 - g - g)*(lnc);
1354
1374
         for(j=s; j<sptree.nnode; j++)  sptree.nodes[j].age *= c;
1356
1376
         for(i=0; i<g; i++) {
1357
1377
            y = data.rgene[i]; 
1358
1378
            ynew = data.rgene[i] /= c;
1359
 
            lnacceptance+=logPriorRatioGamma(ynew, y, data.rgenegamma[0], data.rgenegamma[1]);
 
1379
            lnacceptance += logPriorRatioGamma(ynew, y, data.rgenegamma[0], data.rgenegamma[1]);
1360
1380
         }
1361
1381
         lnLnew = lnpInfinitesites(FixedDs);
1362
1382
         lnacceptance += lnLnew-lnL;
1374
1394
         }
1375
1395
      }
1376
1396
      nround++;
1377
 
      for(i=0; i<np; i++) mx[i]=(mx[i]*(nround-1)+x[i])/nround;
 
1397
      for(i=0; i<np; i++) mx[i] = (mx[i]*(nround-1)+x[i])/nround;
 
1398
      Pjump[0] = naccept[0]/((s-1)*nround);
 
1399
      Pjump[1] = naccept[1]/(g*nround);
 
1400
      Pjump[2] = naccept[2]/nround;
 
1401
      Pjump[3] = naccept[3]/(g*2*nround);
1378
1402
 
1379
1403
      if((ir+1)%max2(mcmc.sampfreq, mcmc.sampfreq*mcmc.nsample/1000)==0) {
1380
1404
         printf("\r%3.0f%%", (ir+1.)/(mcmc.nsample*mcmc.sampfreq)*100.);
1381
 
         printf(" %4.2f %4.2f %4.2f %4.2f  ", naccept[0]/((s-1)*nround),naccept[1]/(g*nround),naccept[2]/nround,naccept[3]/(g*2*nround));
1382
 
 
1383
 
 
 
1405
         printf(" %4.2f %4.2f %4.2f %4.2f  ", Pjump[0], Pjump[1], Pjump[2], Pjump[3]);
1384
1406
         if(np<nxpr[0]+nxpr[1]) { nxpr[0]=com.np; nxpr[1]=0; }
1385
1407
         for(j=0; j<nxpr[0]; j++) printf(" %5.3f", mx[j]);
1386
1408
         if(np>nxpr[0]+nxpr[1] && nxpr[1]) printf(" -");
1399
1421
   free(FixedDs);
1400
1422
 
1401
1423
   printf("\nSummarizing MCMC results, time reset.\n");
1402
 
   DescriptiveStatisticsSimpleMCMCTREE(NULL, com.mcmcf, 1, 1);
 
1424
   DescriptiveStatisticsSimpleMCMCTREE(NULL, com.mcmcf, 1);
1403
1425
 
1404
1426
   exit(0);
1405
1427
}
1443
1465
   possible, even though this means that the chain might start from a poor 
1444
1466
   place.
1445
1467
*/
1446
 
   int np=0, i,j,k, dad, nchanges, g=data.ngene;
 
1468
   int np=0, i,j,k, jj, dad, nchanges, g=data.ngene;
1447
1469
   double maxlower=0; /* rough age for root */
1448
 
   double *p=sptree.nodes[sptree.root].pfossil, smallt=(p[1]-p[0])/sptree.nspecies*2;
 
1470
   double *p=sptree.nodes[sptree.root].pfossil, smallt=(p[1]-p[0])/(40*sqrt(sptree.nspecies));
1449
1471
   double a_r=data.rgenegamma[0], b_r=data.rgenegamma[1], a,b, smallr=1e-3, d;
 
1472
   double AgeLow[NS]={0}, tz;
1450
1473
 
1451
1474
   com.rgene[0]=-1;  /* com.rgene[] is not used.  -1 to force crash */
1452
1475
   puts("\ngetting initial values to start MCMC.");
1453
1476
 
1454
 
   /* set up rough time unit by looking at the fossil info */
1455
1477
   if(com.TipDate) {  /* TipDate model */
 
1478
      /* set up initial node ages by looking at the minimum ages at each node */
 
1479
      if(sptree.nodes[sptree.root].fossil == BOUND_F)
 
1480
         sptree.nodes[sptree.root].age = p[0] + (p[1]-p[0])*rndu();
 
1481
      else
 
1482
         error2("\nthere should be something on the root age for the TipDate model\n");
 
1483
 
 
1484
      for(j=0; j<sptree.nspecies; j++) {
 
1485
         tz = sptree.nodes[j].age;
 
1486
         for(k=sptree.nodes[j].father; k!=-1; k=sptree.nodes[k].father)
 
1487
            if(tz < AgeLow[k-sptree.nspecies]) 
 
1488
               break;
 
1489
            else 
 
1490
               AgeLow[k-sptree.nspecies] = tz;
 
1491
      }
 
1492
      for(j=sptree.nspecies+1; j<sptree.nnode; j++) {
 
1493
         jj = j-sptree.nspecies;
 
1494
         sptree.nodes[j].age = AgeLow[jj] + (sptree.nodes[sptree.nodes[j].father].age - AgeLow[jj])*(0.5+0.5*rndu());
 
1495
      }
 
1496
 
1456
1497
      for(i=0; i<1000; i++) {
1457
1498
         for(j=0,nchanges=0; j<sptree.nspecies; j++)  {
1458
1499
            for(k=j; k!=sptree.root; k=dad) {
1459
1500
               dad = sptree.nodes[k].father;
1460
1501
               if(sptree.nodes[dad].age <= sptree.nodes[k].age) {
1461
 
                  sptree.nodes[dad].age = max2(sptree.nodes[k].age, smallt) * (1 + smallt*2*rndu());
 
1502
                  sptree.nodes[dad].age = max2(sptree.nodes[k].age, smallt) * (1 + smallt*0.1*rndu());
1462
1503
                  nchanges ++;
1463
1504
               }
1464
1505
            }
1465
1506
         }
 
1507
         if(nchanges) puts("nchanges should be 0 here??");
1466
1508
         if(!nchanges) break;
1467
1509
      }
1468
1510
      if(sptree.nodes[sptree.root].fossil == BOUND_F) {
1469
1511
         if(sptree.nodes[sptree.root].age<p[0]) {
1470
 
            sptree.nodes[sptree.root].age = p[0] + (p[1]-p[0])*(0.2+rndu()*1.6);
 
1512
            sptree.nodes[sptree.root].age = p[0] + (p[1]-p[0])*rndu();
1471
1513
         }
1472
1514
      }
 
1515
 
1473
1516
   }
1474
1517
   else {             /* not TipDate model */
 
1518
      /* set up initial node ages by looking at the fossil info */
1475
1519
      for(j=sptree.nspecies; j<sptree.nnode; j++)  {
1476
1520
         sptree.nodes[j].age = -1;
1477
1521
         if(sptree.nodes[j].fossil == 0) continue;
1561
1605
   if(data.pfossilerror[0]) {
1562
1606
      a = data.pfossilerror[0];
1563
1607
      b = data.pfossilerror[1];
1564
 
      sptree.Pfossilerr = a/(a+b)*(0.4+0.6*rndu());
 
1608
      data.Pfossilerr = a/(a+b)*(0.4+0.6*rndu());
1565
1609
      np ++;
1566
1610
   }
1567
1611
 
1631
1675
   }
1632
1676
   if(data.pfossilerror[0]) {
1633
1677
      if(firsttime && fout)  fprintf(fout, "\tpFossilErr");
1634
 
      x[np++] = sptree.Pfossilerr;
 
1678
      x[np++] = data.Pfossilerr;
1635
1679
   }
1636
1680
 
1637
1681
   if(np!=com.np) {
1647
1691
}
1648
1692
 
1649
1693
 
1650
 
double lnpriorTimesBDS_Approach2 (void)
 
1694
double lnpriorTimesBDS_Approach1 (void)
1651
1695
{
1652
 
/* Approach 2 in Tanja's notes.  27 July 2011.  
 
1696
/* Approach 1 or Theorem #.# in Stadler & Yang (2013 Syst Biol).  Based on Tanja's notes of 27 July 2011.
1653
1697
   This calculates g(z*), which is constant throughout the MCMC, and thus wastes time. 
1654
1698
*/
1655
 
   int i,j, k;
 
1699
   int j, k;
1656
1700
   double lambda=data.BDS[0], mu=data.BDS[1], rho=data.BDS[2], psi=data.BDS[3];
1657
 
   double lnp=0, t, t1=sptree.nodes[sptree.root].age, p0, p1, gt, gt1, e, c1, c2;
1658
 
   double a, b, tailL, tailR, thetaL, thetaR;
 
1701
   double lnp=0, t, t1=sptree.nodes[sptree.root].age, gt, gt1, a, e, c1, c2;
1659
1702
   double zstar, z0, z1;
1660
1703
 
1661
1704
   if(com.ndata>1 && com.TipDate) 
1667
1710
   if(psi==0 && fabs(lambda-mu)<1e-20) {
1668
1711
      c1 = 1/t1 + rho*lambda;
1669
1712
      for(j=sptree.nspecies; j<sptree.nnode; j++)
1670
 
         if(j!=sptree.root) {
 
1713
         if(j != sptree.root) {
1671
1714
            c2 = 1 + rho*lambda*sptree.nodes[j].age;
1672
1715
            lnp += log(c1/(c2*c2));
1673
1716
         }
1676
1719
      a = lambda - rho*lambda - mu;
1677
1720
      e = exp((mu - lambda)*t1);
1678
1721
      c1 = (rho*lambda + a*e)/(1 - e);
 
1722
      if(fabs(1-e) < 1E-100)
 
1723
         printf("e too close to 1..");
1679
1724
      /* a & c1 are constant over loop */
1680
1725
      for(j=sptree.nspecies; j<sptree.nnode; j++)
1681
1726
         if(j!=sptree.root) {
1682
1727
            e = exp((mu - lambda)*sptree.nodes[j].age);
1683
1728
            c2 = (lambda-mu)/(rho*lambda + a*e);
1684
1729
            c2 *= c2*e*c1;
 
1730
            if(c2 < 1E-300 || c2>1E300)
 
1731
               printf("c2 not numeric.");
1685
1732
            lnp += log(c2);
1686
1733
         }
1687
1734
   }
1689
1736
      c1 = sqrt(square(lambda - mu - psi) + 4*lambda*psi);
1690
1737
      c2 = -(lambda - mu - 2*lambda*rho - psi)/c1;
1691
1738
      gt1 = 1/( exp(-c1*t1)*(1-c2) + (1+c2) );
 
1739
      if(gt1<-1E-300 || gt1>1E300)
 
1740
         printf("gt1 not numeric.");
1692
1741
 
1693
1742
      for(j=sptree.nspecies; j<sptree.nnode; j++)  {
1694
1743
         if(j==sptree.root) continue;
1717
1766
   return(lnp);
1718
1767
}
1719
1768
 
1720
 
double p0qt_BDS(double t, double lambda, double mu, double rho, double psi, double *p0t)
 
1769
double logqtp0_BDS_Approach2 (double t, double lambda, double mu, double rho, double psi, double *p0t)
1721
1770
{
1722
 
/* This calculates p0t and qt.
 
1771
/* This calculates p0t and log(qt) for Approach 2 of Stadler & Yang (2013 Syst Biol).
1723
1772
*/
1724
 
   double c1, c2, e, r, qt;
 
1773
   double c1, c2, e=0, r, logqt;
1725
1774
 
1726
 
   if(psi==0 && fabs(lambda-mu)<1e-20) {
 
1775
   if(psi==0 && fabs(lambda-mu)<=1e-20) {
1727
1776
      r = 1 + rho*lambda*t;
1728
1777
      *p0t = 1 - rho/(1 + rho*lambda*t);
1729
 
      qt = 4*r*r;
 
1778
      logqt = log(4*r*r);
1730
1779
   }
1731
1780
   else if(psi==0 && fabs(lambda-mu)>1e-20) {
1732
1781
      e = exp((mu - lambda)*t);
1733
1782
      r = (rho*lambda + (lambda-rho*lambda-mu)*e) / (lambda-mu);
1734
1783
      *p0t = 1 - rho*(lambda - mu)/(rho*lambda + (lambda-rho*lambda-mu)*e);
1735
 
      qt = 4*r*r/e;
 
1784
      logqt = log(4*r*r) - (mu - lambda)*t;
1736
1785
   }
1737
1786
   else {
1738
1787
      c1 = sqrt(square(lambda - mu - psi) + 4*lambda*psi);
1742
1791
      r = (e*(1-c2) - (1+c2)) / (e*(1-c2) + (1+c2));
1743
1792
 
1744
1793
      *p0t = (lambda + mu + psi + c1*r) / (2*lambda);
1745
 
      qt = 2*(1-c2*c2) + e*square(1-c2) + square(1+c2)/e;
 
1794
      logqt = c1*t;
 
1795
      logqt += log( e*e*square(1-c2) + e*2*(1-c2*c2) + square(1+c2) );
1746
1796
   }
1747
 
   /* printf("lmrp %8.6f %8.6f %7.4f %7.4f t=%7.4f p0qt= %7.4f %7.4f\n", lambda, mu, rho, psi, t, *p0t, qt); */
1748
 
 
1749
 
   return(qt);
 
1797
   /*
 
1798
   if(e<=1e-300)
 
1799
      printf("lmrp %8.4f %8.4f %7.4f %7.4f t=%7.4f p0 logqt= %7.4f %7.4e\n", lambda, mu, rho, psi, t, *p0t, logqt);
 
1800
   */
 
1801
   return(logqt);
1750
1802
}
1751
1803
 
1752
 
 
1753
 
double lnpriorTimesTipDate (void)
 
1804
double lnpriorTimesTipDate_Approach2 (void)
1754
1805
{
1755
 
/* Equation 6 in Stadler T. 2010. Sampling-through-time in birth-death trees. 
1756
 
   J Theor Biol 267:396-404.
 
1806
/* Approach 2 of Stadler & Yang (2013 Syst Biol).  The BDS parameters are assumed to be fixed.
 
1807
   This ignores terms involving y, which are constants when the BDSS parameters are fixed.
1757
1808
   t1 is root age.
1758
1809
*/
1759
1810
   int i, m=0, n=sptree.nspecies, k=0;
1760
1811
   double lambda=data.BDS[0], mu=data.BDS[1], rho=data.BDS[2], psi=data.BDS[3];
1761
 
   double lnp=0, t, t1=sptree.nodes[sptree.root].age, p0, p1, qt, qt1, z, e;
1762
 
   double a, b, tailL, tailR, thetaL, thetaR;
 
1812
   double lnp=0, t1=sptree.nodes[sptree.root].age, p0, p1, logqt, logqt1, p0t1, z, e;
1763
1813
 
1764
1814
   if(com.ndata>1 && com.TipDate) 
1765
1815
      error2("don't know how ndata works with TipDate.");
1778
1828
   if(rho>0 && n<2)
1779
1829
      error2("we want more than 2 extant sampled species");
1780
1830
 
1781
 
   /* the numerator */
1782
 
   qt1 = p0qt_BDS(t1, lambda, mu, rho, psi, &p0);
1783
 
   lnp -= 2*log(qt1);
 
1831
   /* the numerator f(x, z | L(tmrca) = 2).   */
 
1832
   logqt1 = logqtp0_BDS_Approach2(t1, lambda, mu, rho, psi, &p0t1);
 
1833
   lnp -= 2*logqt1;
1784
1834
   for(i=sptree.nspecies; i<sptree.nnode; i++) { /* loop over n+m-1 internal node ages except t1  */
1785
1835
      if(i == sptree.root)  continue;
1786
 
      qt = p0qt_BDS(sptree.nodes[i].age, lambda, mu, rho, psi, &p0);
1787
 
      lnp -= log(qt);
 
1836
      logqt = logqtp0_BDS_Approach2(sptree.nodes[i].age, lambda, mu, rho, psi, &p0);
 
1837
      lnp -= logqt;
1788
1838
   }
1789
1839
 
1790
 
   /* the denominator */
 
1840
   /* the denominator, h(tmrca, n) */
1791
1841
   if(rho) {
1792
1842
      if(fabs(lambda-mu-psi) > 1e-20) {
1793
 
         e = exp(-(lambda - mu - psi)*t1);
 
1843
         e = exp(-(lambda - mu - psi)*t1);   /* this causes overflow for large t */
1794
1844
         z = lambda*rho + (lambda - lambda*rho - mu - psi)*e;
1795
1845
         p1 = rho*square(lambda - mu - psi)*e/(z*z);
1796
1846
         if(p1<=0)
1797
 
            printf("p1 = %.6f <=0\n", p1);
 
1847
            printf("lnpriorTimesTipDate: p1 = %.6f <=0 (t1 = %9.5g)\n", p1, t1);
1798
1848
         lnp -= log(p1*p1);
1799
1849
         if(n>2) {
1800
1850
            z = rho*lambda*(1 - e)/z;
1810
1860
      }
1811
1861
   }
1812
1862
   else {  /* rho = 0 for viruses */
1813
 
      p0qt_BDS(t1, lambda, mu, rho, psi, &p0);
1814
 
      z = square(1 - p0);
1815
 
      lnp -= log(z);
1816
 
   }
1817
 
 
1818
 
   /* prior on t1 */
1819
 
   if(sptree.nodes[sptree.nspecies].fossil != BOUND_F)
1820
 
      error2("Only bounds for the root age are implemented.");
1821
 
   lnp += lnptCalibrationDensity(t1, sptree.nodes[sptree.nspecies].fossil, sptree.nodes[sptree.nspecies].pfossil);
1822
 
 
 
1863
      lnp -= 2*log(1 - p0t1);
 
1864
   }
 
1865
 
 
1866
   /* prior on t1 */
 
1867
   if(sptree.nodes[sptree.nspecies].fossil != BOUND_F)
 
1868
      error2("Only bounds for the root age are implemented.");
 
1869
   lnp += lnptCalibrationDensity(t1, sptree.nodes[sptree.nspecies].fossil, sptree.nodes[sptree.nspecies].pfossil);
 
1870
 
 
1871
   return(lnp);
 
1872
}
 
1873
 
 
1874
double p01t_BDSEquation6Stadler2010 (double t, double lambda, double mu, double rho, double psi, double *p0t)
 
1875
{
 
1876
/* This calculates p0t and p1t, defined in equations 1 & 2 in Stadler (2010).
 
1877
*/
 
1878
   double c1, c2, e, r, p1t;
 
1879
 
 
1880
   if(fabs(lambda-mu)<1e-20 && psi==0) {
 
1881
      r = 1/(1/rho + lambda*t);
 
1882
      *p0t = 1 - r;
 
1883
      p1t = r*r/rho;
 
1884
   }
 
1885
   else if(fabs(lambda-mu)>1e-20 && psi==0) {
 
1886
      e = exp((mu - lambda)*t);
 
1887
      r = rho*(lambda-mu) / (rho*lambda + (lambda-mu-rho*lambda)*e);
 
1888
      *p0t = 1 - r;
 
1889
      p1t = r*r/rho*e;
 
1890
   }
 
1891
   else {
 
1892
      c1 = sqrt(square(lambda - mu - psi) + 4*lambda*psi);
 
1893
      if(c1==0) error2("c1==0");
 
1894
      c2 = -(lambda - mu - 2*lambda*rho - psi)/c1;
 
1895
      e = exp(-c1*t);
 
1896
      r = (e*(1-c2) - (1+c2)) / (e*(1-c2) + (1+c2));
 
1897
 
 
1898
      *p0t = (lambda + mu + psi + c1*r) / (2*lambda);
 
1899
      p1t = 4*rho/( 2*(1-c2*c2) + e*square(1-c2) + square(1+c2)/e );
 
1900
   }
 
1901
 
 
1902
   /* printf("lmrp %8.6f %8.6f %7.4f %7.4f t=%7.4f p01= %7.4f %7.4f\n", lambda, mu, rho, psi, t, *p0t, p1t); */
 
1903
 
 
1904
   return(p1t);
 
1905
}
 
1906
 
 
1907
double lnpriorTimesTipDateEquation6Stadler2010 (void)
 
1908
{
 
1909
/* Equation 6 in Stadler T. 2010. Sampling-through-time in birth-death trees. 
 
1910
   J Theor Biol 267:396-404.
 
1911
   t1 is root age.
 
1912
*/
 
1913
   int i, m, n, k=0;
 
1914
   double lambda=data.BDS[0], mu=data.BDS[1], rho=data.BDS[2], psi=data.BDS[3];
 
1915
   double lnp=0, t, t1=sptree.nodes[sptree.root].age, p0, p1, z, e;
 
1916
   double a, b, tailL, tailR, thetaL, thetaR;
 
1917
 
 
1918
   if(com.ndata>1 && com.TipDate) 
 
1919
      error2("don't know how ndata works with TipDate.");
 
1920
 
 
1921
   if(lambda<=0 || mu<=0 || (rho<=0 && psi<=0))
 
1922
      error2("wrong B-D-S parameters..");
 
1923
   for(i=0,m=0; i<sptree.nspecies; i++)
 
1924
      m += (sptree.nodes[i].age > 0);
 
1925
   n = sptree.nspecies - m;
 
1926
 
 
1927
   if(n>2) {
 
1928
      if(fabs(lambda-mu)<1e-20)
 
1929
         z = 1 + 1/(rho*lambda*t1);
 
1930
      else {
 
1931
         e = exp((mu-lambda)*t1);
 
1932
         z = (lambda*rho + (lambda-mu-lambda*rho)*e) / (rho*lambda*(1-e));
 
1933
      }
 
1934
      lnp = (n-2)*log(z*z);
 
1935
   }
 
1936
   if(n>1) {
 
1937
      p1 = p01t_BDSEquation6Stadler2010(t1, lambda, mu, rho, 0, &p0);
 
1938
      lnp -= log((n-1)*p1*p1);
 
1939
   }
 
1940
   /*  this term is constant when the BDS parameters are fixed.
 
1941
   lnp += (n+m-2)*log(lambda) + (k+m)*log(psi);
 
1942
   */
 
1943
 
 
1944
   p1 = p01t_BDSEquation6Stadler2010(t1, lambda, mu, rho, psi, &p0);
 
1945
   lnp += 2*log(p1);
 
1946
   for(i=sptree.nspecies; i<sptree.nnode; i++)  /* loop over n+m-1 internal node ages except t1  */
 
1947
      if(i != sptree.root) 
 
1948
         lnp += log( p01t_BDSEquation6Stadler2010(sptree.nodes[i].age, lambda, mu, rho, psi, &p0) );
 
1949
 
 
1950
   /* the y terms do not change when the BDS parameters are fixed.
 
1951
   for(i=0; i<sptree.nspecies; i++) {
 
1952
      if( (t=sptree.nodes[i].age) ) {
 
1953
         p1 = p01t_BDSEquation6Stadler2010(t, lambda, mu, rho, psi, &p0);
 
1954
         lnp += log(p0/p1);
 
1955
      }
 
1956
   }
 
1957
   */
 
1958
 
 
1959
   /* prior on t1 */
 
1960
   if(sptree.nodes[sptree.nspecies].fossil != BOUND_F)
 
1961
      error2("Only bounds for the root age are implemented.");
 
1962
   lnp += lnptCalibrationDensity(t1, sptree.nodes[sptree.nspecies].fossil, sptree.nodes[sptree.nspecies].pfossil);
1823
1963
   return(lnp);
1824
1964
}
1825
1965
 
1837
1977
   else {
1838
1978
      expmlt = exp(mlt);
1839
1979
      lnp = P0t_YR07(expmlt);
 
1980
      if(lnp<0) puts("this should not be negative, in lnPDFkernelBDS_YR07");
1840
1981
      lnp = log( lnp*lnp * lambda/(vt1*rho) * expmlt );
1841
1982
   }
1842
1983
   return(lnp);
1995
2136
   }
1996
2137
   if(debug==1) printf("\npdf = %.12f\n", exp(lnpt));
1997
2138
 
 
2139
 
 
2140
/*
 
2141
{
 
2142
double t1=sptree.nodes[4].age, t2=sptree.nodes[5].age, t3=sptree.nodes[6].age, Gt2, Gt3;
 
2143
Gt2 = CDFkernelBDS_YR07(t2, t1, vt1, lambda, mu, rho);
 
2144
Gt3 = CDFkernelBDS_YR07(t3, t1, vt1, lambda, mu, rho);
 
2145
if(sptree.nodes[5].usefossil==0 && sptree.nodes[6].usefossil==0)
 
2146
   lnpt += log(1.0/2);
 
2147
if(sptree.nodes[5].usefossil==1 && sptree.nodes[6].usefossil==0)
 
2148
   if(t3<t2) lnpt += log(Gt2);
 
2149
   else      lnpt += log(1-Gt2);
 
2150
if(sptree.nodes[5].usefossil==0 && sptree.nodes[6].usefossil==1)
 
2151
   if(t2<t3) lnpt += log(Gt3);
 
2152
   else      lnpt += log(1-Gt3);
 
2153
 
 
2154
}
 
2155
*/
1998
2156
   return (lnpt);
1999
2157
}
2000
2158
 
2293
2451
/* This prepares for lnpriorTimes() under models of fossil errors.  It calculates 
2294
2452
   the scaling factor for the probability density of times for a given combination
2295
2453
   of the indicator variables, indicating which fossil is in error and not used.
2296
 
   The combination in which all fossils are wrong is excluded.
 
2454
   nMinCorrect is the minimum number of correct fossils.
2297
2455
*/
2298
 
   int  is,i, icom, nused=0, it;
2299
 
   int nMinCorrect = (int)data.pfossilerror[2];
2300
 
   int ncomFerr = (1<<sptree.nfossil) - 1 - (nMinCorrect==2?sptree.nfossil:0);
 
2456
   int nMinCorrect = (int)data.pfossilerror[2], ncomFerr, is,i, icom, nused=0, it;
 
2457
   double t;
2301
2458
   int debug=1;
2302
2459
 
2303
 
   if(data.pfossilerror[0]==0 || sptree.nfossil>MaxNFossils || sptree.nfossil<2) 
 
2460
   if(data.pfossilerror[0]==0 || sptree.nfossil>MaxNFossils || sptree.nfossil<nMinCorrect)
2304
2461
      error2("Fossil-error model is inappropriate."); 
2305
2462
 
2306
 
   sptree.CcomFossilErr = (double*)malloc(ncomFerr*sizeof(double));
2307
 
   if(sptree.CcomFossilErr == NULL) error2("oom for CcomFossilErr");
 
2463
   for (i=nMinCorrect,ncomFerr=0; i<=sptree.nfossil; i++) {
 
2464
      ncomFerr += (int)( Binomial((double)sptree.nfossil, i, &t) + .1 );
 
2465
      if(t!=0) error2("too many fossils to deal with? ");
 
2466
   }
 
2467
 
 
2468
   data.CcomFossilErr = (double*)malloc(ncomFerr*sizeof(double));
 
2469
   if(data.CcomFossilErr == NULL) error2("oom for CcomFossilErr");
2308
2470
 
2309
2471
   /* cycle through combinations of indicators. */
2310
 
   for (i=0,icom=0; i < (1<<sptree.nfossil)-1; i++) {
 
2472
   for (i=0,icom=0; i < (1<<sptree.nfossil); i++) {
2311
2473
      it = i;
2312
2474
      for (is=sptree.nspecies, nused=0; is<sptree.nnode; is++) {
2313
2475
         if(sptree.nodes[is].fossil) {
2317
2479
            if(debug) printf("%2d (%2d)", sptree.nodes[is].usefossil, is+1);
2318
2480
         }
2319
2481
      }
2320
 
      if(nMinCorrect==2 && nused<2) continue;
2321
 
      icom++;
2322
 
      if(debug) printf("\n\n********* Com %2d/%2d:  %2d fossils used", icom, ncomFerr, nused);
2323
 
      sptree.CcomFossilErr[icom-1] = getScaleFossilCombination();
 
2482
      if(nused<nMinCorrect) continue;
 
2483
      if(debug) printf("\n\n********* Combination %2d/%2d:  %2d fossils used", icom+1, ncomFerr, nused);
 
2484
      data.CcomFossilErr[icom++] = log( getScaleFossilCombination() );
2324
2485
   }
2325
2486
   return(0);
2326
2487
}
2327
2488
 
2328
2489
 
 
2490
 
2329
2491
double lnpriorTimes (void)
2330
2492
{
2331
 
/* This calculates the prior density of node times in the master species tree:
2332
 
   sptree.nodes[].age. 
 
2493
/* This calculates the prior density of node times in the master species tree.
 
2494
      Node ages are in sptree.nodes[].age. 
2333
2495
*/
2334
2496
   int nMinCorrect = (int)data.pfossilerror[2];
2335
 
   int ncomFerr = (1<<sptree.nfossil) - 1 - (nMinCorrect==2?sptree.nfossil:0);
2336
2497
   int  is, i,icom, nused, it;
2337
 
   double pE = sptree.Pfossilerr, ln1pE=log(1-pow(pE,(double)sptree.nfossil));
 
2498
   double pE = data.Pfossilerr, lnpE, ln1pE, pEconst=0, lnC, lnNC;
2338
2499
   double lnpt=0, scaleF=-1e300, y;
 
2500
   int debug=0;
2339
2501
 
2340
2502
   if(sptree.nfossil==0 || com.TipDate) {
2341
 
      if(1)
2342
 
         lnpt = lnpriorTimesBDS_Approach2();
2343
 
      else
2344
 
         lnpt = lnpriorTimesTipDate();
2345
 
      
2346
 
   }
2347
 
   else if(data.pfossilerror[0]==0)  /* no fossil errors in model */
2348
 
      lnpt = lnptC() + lnptNCgiventC();
 
2503
      if(1)        lnpt = lnpriorTimesBDS_Approach1();
 
2504
      else if(0)   lnpt = lnpriorTimesTipDate_Approach2();
 
2505
      else if(0)   lnpt = lnpriorTimesTipDateEquation6Stadler2010();
 
2506
   }
 
2507
   else if(data.pfossilerror[0]==0) { /* no fossil errors in model */
 
2508
      lnC = lnptC();
 
2509
      lnNC = lnptNCgiventC();
 
2510
      lnpt = lnC + lnNC;
 
2511
 
 
2512
      if(debug) printf("\nftC ftNC ft = %9.5f%9.5f%9.5f", exp(lnC), exp(lnNC), exp(lnC)*exp(lnNC));
 
2513
   }
2349
2514
   else {   /* fossil errors: cycle through combinations of indicators, using nMinCorrect. */
2350
 
      for(i=0,icom=0; i < (1<<sptree.nfossil) - 1; i++) {
 
2515
      if(nMinCorrect) {
 
2516
         pEconst = y = pow(pE, (double)sptree.nfossil);
 
2517
         if(y<1e-300) error2("many fossils & small pE.  Rewrite the code?");
 
2518
         for(i=1; i<nMinCorrect; i++) /* i is the number of correct fossils used */
 
2519
            pEconst += y *= (sptree.nfossil-i+1)*(1-pE)/(i*pE);
 
2520
         pEconst = log(1 - pEconst);
 
2521
      }
 
2522
      lnpE=log(pE); ln1pE=log(1-pE);
 
2523
      for(i=0,icom=0; i < (1<<sptree.nfossil); i++) { /* sum over U */
2351
2524
         it = i;
2352
2525
         for(is=sptree.nspecies, nused=0; is<sptree.nnode; is++) {
2353
2526
            if(sptree.nodes[is].fossil) {
2354
2527
               sptree.nodes[is].usefossil = 1 - it%2;
2355
2528
               nused += sptree.nodes[is].usefossil;
2356
 
               it/=2;
 
2529
               it /= 2;
2357
2530
            }
2358
2531
         }
2359
 
         if(nMinCorrect==2 && nused<2) continue;
2360
 
         icom ++;
2361
 
 
2362
 
         y = nused*log(1-pE)+(sptree.nfossil-nused)*log(pE)-ln1pE
2363
 
           - log(sptree.CcomFossilErr[icom-1]) + lnptC() + lnptNCgiventC();
 
2532
         if(nused<nMinCorrect) continue;
 
2533
         lnC = lnptC();
 
2534
         lnNC = lnptNCgiventC();
 
2535
 
 
2536
         y = nused*ln1pE + (sptree.nfossil-nused)*lnpE - pEconst - data.CcomFossilErr[icom++]
 
2537
           + lnC + lnNC;
 
2538
 
 
2539
 
 
2540
         if(debug) 
 
2541
            printf("\nU%d%d ftC ftNC ft = %9.5f%9.5f%9.5f", sptree.nodes[5].usefossil, sptree.nodes[6].usefossil, 
 
2542
                      exp(lnC), exp(lnNC), exp(lnC)*exp(lnNC));
2364
2543
 
2365
2544
         if(y < scaleF + 100)
2366
2545
            lnpt += exp(y-scaleF);
2367
2546
         else {         
2368
 
            lnpt = lnpt*exp(scaleF-y) + 1;
 
2547
            lnpt = 1;
2369
2548
            scaleF = y;
2370
2549
         }
2371
2550
         lnpt += exp(y-scaleF);
2372
2551
      }
2373
2552
      lnpt = scaleF + log(lnpt);
2374
2553
   }
 
2554
   if(debug) exit(0);
2375
2555
 
2376
2556
   return (lnpt);
2377
2557
}
2382
2562
   for each fossil given the times and pE.
2383
2563
*/
2384
2564
   int nMinCorrect = (int)data.pfossilerror[2];
2385
 
   int ncomFerr = (1<<sptree.nfossil) - 1 - (nMinCorrect==2?sptree.nfossil:0);
2386
2565
   int  is,i,icom, k, nf=sptree.nfossil, nused, it;
2387
 
   double pE = sptree.Pfossilerr, ln1pE=log(1-pow(pE,(double)nf));
2388
 
   double pt, pEf[MaxNFossils]={0}, scaleF=-1e300, y;
2389
 
   char Ef[MaxNFossils];  /* indicators of fossil errors */
 
2566
   double pE = data.Pfossilerr, lnpE, ln1pE, pEconst=0;
 
2567
   double pt, pU[MaxNFossils]={0}, scaleF=-1e300, y;
 
2568
   char U[MaxNFossils];  /* indicators of fossil errors */
2390
2569
 
2391
2570
   if(data.priortime==1) 
2392
2571
      error2("beta kernel for prior time not yet implemented for model of fossil errors.");
 
2572
   if(nMinCorrect) {
 
2573
      pEconst = y = pow(pE, (double)sptree.nfossil);
 
2574
      for(i=1; i<nMinCorrect; i++) /* i is the number of correct fossils used */
 
2575
         pEconst += y *= (sptree.nfossil-i+1)*(1-pE)/(i*pE);
 
2576
      pEconst = log(1 - pEconst);
 
2577
   }
2393
2578
 
2394
 
   for(i=0,icom=0,pt=0; i < (1<<sptree.nfossil) - 1; i++) {
 
2579
   lnpE=log(pE); ln1pE=log(1-pE);
 
2580
   for(i=0,icom=0,pt=0; i < (1<<sptree.nfossil); i++) { /* sum over U */
2395
2581
      it = i;
2396
2582
      for(is=sptree.nspecies, k=0, nused=0; is<sptree.nnode; is++) {
2397
2583
         if(sptree.nodes[is].fossil) {
2398
2584
            sptree.nodes[is].usefossil = 1 - it%2;
2399
2585
            nused += sptree.nodes[is].usefossil;
2400
 
            Ef[k++] = !sptree.nodes[is].usefossil;
 
2586
            U[k++] = !sptree.nodes[is].usefossil;
2401
2587
            it /= 2;
2402
2588
         }
2403
2589
      }
2404
 
      if(nMinCorrect==2 && nused<2) continue;
2405
 
      icom ++;
2406
 
      if(k != nf) error2("k == nf?");
2407
 
 
2408
 
      y = nused*log(1-pE)+(nf-nused)*log(pE)-ln1pE
2409
 
        - log(sptree.CcomFossilErr[icom-1]) + lnptC() + lnptNCgiventC();
2410
 
 
2411
 
      if(y < scaleF + 200) {
 
2590
      if(nused<nMinCorrect) continue;
 
2591
      if(k != nf) error2("k != nf in getPfossilerr()?");
 
2592
 
 
2593
      y = nused*ln1pE+(nf-nused)*lnpE - pEconst - data.CcomFossilErr[icom++] + lnptC() + lnptNCgiventC();
 
2594
 
 
2595
      if(y < scaleF + 100) {
2412
2596
         pt += y = exp(y-scaleF);
 
2597
         for(k=0; k<nf; k++)
 
2598
            if(U[k]) pU[k] += y;
2413
2599
      }
2414
2600
      else {         /* change scaleF */
2415
 
         pt = pt*exp(scaleF-y) + 1;
2416
 
         for(k=0; k<nf; k++)
2417
 
            pEf[k] *= exp(scaleF-y);
2418
2601
         scaleF = y;
2419
 
         y = 1;
 
2602
         pt = 1;
 
2603
         for(k=0; k<nf; k++) pU[k] = U[k];  /*  1 if U[k]=1 or 0 if U[k]=0 */
2420
2604
      }
2421
 
      for(k=0; k<nf; k++)
2422
 
         if(Ef[k]) pEf[k] += y;
2423
2605
   }
2424
 
   for(k=0; k<nf; k++)
2425
 
      pEf[k] /= pt;
2426
 
 
2427
 
   for(k=0; k<nf; k++)
2428
 
      postEFossil[k] = (postEFossil[k]*(nround-1) + pEf[k])/nround;
 
2606
   for(k=0; k<nf; k++) pU[k] /= pt;
 
2607
   for(k=0; k<nf; k++)
 
2608
      postEFossil[k] = (postEFossil[k]*(nround-1) + pU[k])/nround;
2429
2609
   
2430
2610
   return (0);
2431
2611
}
2500
2680
      else                   yb[1] = 99;
2501
2681
 
2502
2682
      y = log(t);
2503
 
      ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
2683
      ynew = y + finetune*rndBactrian(mBactrian);
2504
2684
      ynew = reflect(ynew, yb[0], yb[1]);
2505
2685
      sptree.nodes[is].age = tnew = exp(ynew);
 
2686
 
2506
2687
      lnacceptance = ynew - y;
2507
2688
      lnpTnew = lnpriorTimes();
2508
2689
      lnacceptance += lnpTnew - data.lnpT;
2578
2759
      else                   yb[1] = 99;
2579
2760
 
2580
2761
      y = log(t);
2581
 
      ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
2762
      ynew = y + finetune*rndBactrian(mBactrian);
2582
2763
      ynew = reflect(ynew, yb[0], yb[1]);
2583
2764
      sptree.nodes[is].age = tnew = exp(ynew);
2584
2765
      lnacceptance = ynew - y;
2678
2859
   clock=3: the root rate is mu or data.rgene[].  The algorithm cycles through 
2679
2860
   the ancestral nodes and deals with the two daughter branches.
2680
2861
*/
2681
 
   int i, inode, dad=-1, g=data.ngene, s=sptree.nspecies, sons[2], root=sptree.root;
 
2862
   int i, inode, locus, dad=-1, g=data.ngene, s=sptree.nspecies, sons[2];
2682
2863
   double lnpR=-log(2*Pi)/2.0*(2*s-2)*g, t,tA,t1,t2,Tinv[4], detT;
2683
2864
   double zz, r=-1, rA,r1,r2, y1,y2;
2684
 
   double alpha, beta;
 
2865
   double a, b;
2685
2866
 
2686
 
   if(data.priorrate==1 && com.clock==3)
 
2867
   if(com.clock==3 && data.priorrate==1)
2687
2868
      error2("gamma prior for rates for clock3 not implemented yet.");
2688
 
   else if(data.priorrate==1 && com.clock==2) {   /* gamma rate prior, clock2 */
 
2869
   else if(com.clock==2 && data.priorrate==1) {   /* clock2, gamma rate prior */
2689
2870
      lnpR = 0;
2690
 
      for(i=0; i<g; i++) {
2691
 
         alpha = data.rgene[i]*data.rgene[i]/data.sigma2[i];
2692
 
         beta  = data.rgene[i]/data.sigma2[i];
2693
 
         lnpR += (alpha*log(beta) - LnGamma(alpha)) * (2*s-2);
 
2871
      for(locus=0; locus<g; locus++) {
 
2872
         a = data.rgene[locus]*data.rgene[locus]/data.sigma2[locus];
 
2873
         b = data.rgene[locus]/data.sigma2[locus];
 
2874
         lnpR += (a*log(b) - LnGamma(a)) * (2*s-2);
2694
2875
         for(inode=0; inode<sptree.nnode; inode++) {
2695
 
            if(inode==root) continue;
2696
 
            r = sptree.nodes[inode].rates[i];
2697
 
            lnpR += -beta*r + (alpha-1)*log(r);
 
2876
            if(inode==sptree.root) continue;
 
2877
            r = sptree.nodes[inode].rates[locus];
 
2878
            lnpR += -b*r + (a-1)*log(r);
2698
2879
         }
2699
2880
      }
2700
2881
   }
2701
 
   else if(data.priorrate ==0 && com.clock==2) {  /* LN rate prior, clock2 */
2702
 
      for(i=0; i<g; i++)
2703
 
         lnpR -= log(data.sigma2[i])/2.*(2*s-2);
 
2882
   else if(com.clock==2 && data.priorrate ==0) {  /* clock2, LN rate prior */
 
2883
      for(locus=0; locus<g; locus++)
 
2884
         lnpR -= log(data.sigma2[locus])/2.*(2*s-2);
2704
2885
      for(inode=0; inode<sptree.nnode; inode++) {
2705
 
         if(inode==root) continue;
2706
 
         for(i=0; i<g; i++) {
2707
 
            r = sptree.nodes[inode].rates[i];
2708
 
            zz = log(r/data.rgene[i]) + data.sigma2[i]/2;
2709
 
            lnpR += -zz*zz/(2*data.sigma2[i]) - log(r);
 
2886
         if(inode==sptree.root) continue;
 
2887
         for(locus=0; locus<g; locus++) {
 
2888
            r = sptree.nodes[inode].rates[locus];
 
2889
            zz = log(r/data.rgene[locus]) + data.sigma2[locus]/2;
 
2890
            lnpR += -zz*zz/(2*data.sigma2[locus]) - log(r);
2710
2891
         }
2711
2892
      }
2712
2893
   }
2713
 
   else if(data.priorrate ==0 && com.clock==3) {  /* LN rate prior, clock3 */
 
2894
   else if(com.clock==3 && data.priorrate ==0) {  /* clock3, LN rate prior */
2714
2895
      for(inode=0; inode<sptree.nnode; inode++) {
2715
2896
         if(sptree.nodes[inode].nson==0) continue; /* skip the tips */
2716
2897
         dad = sptree.nodes[inode].father;
2717
2898
         for(i=0; i<2; i++) sons[i] = sptree.nodes[inode].sons[i];
2718
2899
         t = sptree.nodes[inode].age;
2719
 
         if(inode==root) tA = 0;
2720
 
         else            tA = (sptree.nodes[dad].age - t)/2;
 
2900
         if(inode==sptree.root) tA = 0;
 
2901
         else                   tA = (sptree.nodes[dad].age - t)/2;
2721
2902
         t1 = (t-sptree.nodes[sons[0]].age)/2;
2722
2903
         t2 = (t-sptree.nodes[sons[1]].age)/2;
2723
2904
         detT = t1*t2+tA*(t1+t2);
2724
2905
         Tinv[0] = (tA+t2)/detT; 
2725
2906
         Tinv[1] = Tinv[2] = -tA/detT; 
2726
2907
         Tinv[3] = (tA+t1)/detT;
2727
 
         for(i=0; i<g; i++) {
2728
 
            rA = (inode==root||dad==root ? data.rgene[i] : sptree.nodes[dad].rates[i]);
2729
 
            r1 = sptree.nodes[sons[0]].rates[i];
2730
 
            r2 = sptree.nodes[sons[1]].rates[i];
2731
 
            y1 = log(r1/rA)+(tA+t1)*data.sigma2[i]/2;
2732
 
            y2 = log(r2/rA)+(tA+t2)*data.sigma2[i]/2;
 
2908
         for(locus=0; locus<g; locus++) {
 
2909
            rA = (inode==sptree.root ? data.rgene[locus] : sptree.nodes[inode].rates[locus]);
 
2910
            r1 = sptree.nodes[sons[0]].rates[locus];
 
2911
            r2 = sptree.nodes[sons[1]].rates[locus];
 
2912
            y1 = log(r1/rA)+(tA+t1)*data.sigma2[locus]/2;
 
2913
            y2 = log(r2/rA)+(tA+t2)*data.sigma2[locus]/2;
2733
2914
            zz = (y1*y1*Tinv[0]+2*y1*y2*Tinv[1]+y2*y2*Tinv[3]);
2734
 
            lnpR -= zz/(2*data.sigma2[i]) + log(detT*square(data.sigma2[i]))/2 + log(r1*r2);
 
2915
            lnpR -= zz/(2*data.sigma2[locus]) + log(detT*square(data.sigma2[locus]))/2 + log(r1*r2);
2735
2916
         }
2736
2917
      }
2737
2918
   }
2738
2919
   return lnpR;
2739
2920
}
2740
2921
 
 
2922
double lnpriorRatioRates (int locus, int inodeChanged, double rold)
 
2923
{
 
2924
/* This calculates the lnpriorRatio when one rate is changed.
 
2925
   If inodeChanged is tip, we sum over 1 term (equation 7 in RY2007).  
 
2926
   If inodeChanged is not tip, we sum over 2 terms.
 
2927
*/
 
2928
   double rnew=sptree.nodes[inodeChanged].rates[locus], lnpRd=0, a, b, z, znew;
 
2929
   double zz, r=-1, rA,r1,r2, y1,y2, t,tA,t1,t2,Tinv[4], detT;
 
2930
   int i, inode, ir, dad=-1, sons[2], OldNew;
 
2931
   
 
2932
   if(com.clock==2 && data.priorrate==0) {         /* clock2, LN rate prior */
 
2933
      z    = log(rold/data.rgene[locus]) + data.sigma2[locus]/2;;
 
2934
      znew = log(rnew/data.rgene[locus]) + data.sigma2[locus]/2;
 
2935
      lnpRd = -log(rnew/rold) - (znew*znew - z*z)/(2*data.sigma2[locus]);
 
2936
   }
 
2937
   else if (com.clock==2 && data.priorrate==1) {   /* clock2, gamma rate prior */
 
2938
      a = data.rgene[locus]*data.rgene[locus]/data.sigma2[locus];
 
2939
      b = data.rgene[locus]/data.sigma2[locus];
 
2940
      lnpRd = -b*(rnew-rold) + (a-1)*log(rnew/rold);
 
2941
   }
 
2942
   else {                                          /* clock3 */
 
2943
      for(ir=0; ir<(sptree.nodes[inodeChanged].nson==0 ? 1 : 2); ir++) {
 
2944
         inode = (ir==0 ? sptree.nodes[inodeChanged].father : inodeChanged);
 
2945
         dad = sptree.nodes[inode].father;
 
2946
         for(i=0; i<2; i++) sons[i] = sptree.nodes[inode].sons[i];
 
2947
         t = sptree.nodes[inode].age;
 
2948
         if(inode==sptree.root) tA = 0;
 
2949
         else                   tA = (sptree.nodes[dad].age - t)/2;
 
2950
         t1 = (t-sptree.nodes[sons[0]].age)/2;
 
2951
         t2 = (t-sptree.nodes[sons[1]].age)/2;
 
2952
         detT = t1*t2+tA*(t1+t2);
 
2953
         Tinv[0] = (tA+t2)/detT; 
 
2954
         Tinv[1] = Tinv[2] = -tA/detT; 
 
2955
         Tinv[3] = (tA+t1)/detT;
 
2956
         for(OldNew=0; OldNew<2; OldNew++) {  /* old rate & new rate */
 
2957
            sptree.nodes[inodeChanged].rates[locus] = (OldNew==0 ? rold : rnew);
 
2958
            rA = (inode==sptree.root ? data.rgene[locus] : sptree.nodes[inode].rates[locus]);
 
2959
            r1 = sptree.nodes[sons[0]].rates[locus];
 
2960
            r2 = sptree.nodes[sons[1]].rates[locus];
 
2961
            y1 = log(r1/rA)+(tA+t1)*data.sigma2[locus]/2;
 
2962
            y2 = log(r2/rA)+(tA+t2)*data.sigma2[locus]/2;
 
2963
            zz = (y1*y1*Tinv[0]+2*y1*y2*Tinv[1]+y2*y2*Tinv[3]);
 
2964
            zz = zz/(2*data.sigma2[locus]) + log(r1*r2);
 
2965
            lnpRd -= (OldNew==0 ? -1 : 1) * zz;
 
2966
         }
 
2967
      }
 
2968
   }
 
2969
   return(lnpRd);
 
2970
}
 
2971
 
 
2972
 
2741
2973
double UpdateRates (double *lnL, double finetune)
2742
2974
{
2743
2975
/* This updates rates under the rate-drift models (clock=2 or 3).
2748
2980
   tree, thus wasting computation, if rates are not correlated across loci.
2749
2981
*/
2750
2982
   int locus, inode, j, g=data.ngene;
2751
 
   double naccept=0, lnpRnew, lnpDinew, lnacceptance, lnLd, rold;
 
2983
   double naccept=0, lnpDinew, lnacceptance, lnLd, r, rnew, lnpRd;
2752
2984
   double yb[2]={-99,99},y,ynew;
2753
2985
 
2754
2986
   if(finetune<=0) error2("steplength = 0 in UpdateRates");
2758
2990
   for(locus=0; locus<g; locus++) {
2759
2991
      for(inode=0; inode<sptree.nnode; inode++) {
2760
2992
         if(inode == sptree.root) continue;
2761
 
         rold = sptree.nodes[inode].rates[locus];
2762
 
         y = log(rold);
2763
 
         ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
2993
         r = sptree.nodes[inode].rates[locus];
 
2994
         y = log(r);
 
2995
         ynew = y + finetune*rndBactrian(mBactrian);
2764
2996
         ynew = reflect(ynew, yb[0], yb[1]);
2765
 
         sptree.nodes[inode].rates[locus] = exp(ynew);
2766
 
         lnacceptance = ynew - y;
 
2997
         sptree.nodes[inode].rates[locus] = rnew = exp(ynew);
 
2998
         lnacceptance = ynew - y;  /* proposal ratio */
2767
2999
 
2768
3000
         UseLocus(locus, 1, mcmc.usedata, 0);  /* copyconP=1 */
2769
3001
 
2771
3003
            FOR(j,sptree.nspecies*2-1) com.oldconP[j]=1;
2772
3004
            LabelOldCondP(inode);
2773
3005
         }
2774
 
 
2775
 
         lnpRnew = lnpriorRates();
2776
 
         lnacceptance += lnpRnew - data.lnpR;
 
3006
         lnpRd = lnpriorRatioRates(locus, inode, r);
 
3007
         lnacceptance += lnpRd;
2777
3008
         lnpDinew = lnpD_locus(locus);
2778
3009
         lnLd = lnpDinew - data.lnpDi[locus];
2779
3010
         lnacceptance +=  lnLd;
2780
 
 
2781
3011
         if(lnacceptance>0 || rndu()<exp(lnacceptance)) {
2782
3012
            naccept++;
2783
3013
            if(mcmc.usedata==1) AcceptRejectLocus(locus,1);
2784
3014
   
2785
 
            data.lnpR = lnpRnew;
 
3015
            data.lnpR += lnpRd;
2786
3016
            data.lnpDi[locus] = lnpDinew;
2787
3017
            *lnL += lnLd;
2788
3018
         }
2789
3019
         else {
2790
3020
            if(mcmc.usedata==1) AcceptRejectLocus(locus,0);
2791
 
            sptree.nodes[inode].rates[locus] = rold;
 
3021
            sptree.nodes[inode].rates[locus] = r;
2792
3022
         }
2793
3023
      }
2794
3024
   }
2821
3051
   for(locus=0; locus<g; locus++) {
2822
3052
      rold = data.rgene[locus];
2823
3053
      y = log(rold);
2824
 
      ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
3054
      ynew = y + finetune*rndBactrian(mBactrian);
2825
3055
      ynew = reflect(ynew, yb[0], yb[1]);
2826
3056
      data.rgene[locus] = exp(ynew);
2827
3057
      lnacceptance = ynew - y;
2875
3105
         if(ip==0 && !com.fix_kappa) {  /* kappa */
2876
3106
            pold = data.kappa[locus];
2877
3107
            y = log(pold);
2878
 
            ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
3108
            ynew = y + finetune*rndBactrian(mBactrian);
2879
3109
            ynew = reflect(ynew, yb[0], yb[1]);
2880
3110
            data.kappa[locus] = pnew = exp(ynew);
2881
3111
            gammaprior = data.kappagamma;
2883
3113
         else {  /* alpha */
2884
3114
            pold = data.alpha[locus];
2885
3115
            y = log(pold);
2886
 
            ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
3116
            ynew = y + finetune*rndBactrian(mBactrian);
2887
3117
            ynew = reflect(ynew, yb[0], yb[1]);
2888
3118
            data.alpha[locus] = pnew = exp(ynew);
2889
3119
            gammaprior = data.alphagamma;
2943
3173
         if(ip==0) {  /* rgene (mu) */
2944
3174
            pold = data.rgene[i];
2945
3175
            y = log(pold);
2946
 
            ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
3176
            ynew = y + finetune*rndBactrian(mBactrian);
2947
3177
            ynew = reflect(ynew, yb[0], yb[1]);
2948
3178
            data.rgene[i] = pnew = exp(ynew);
2949
3179
            gammaprior = data.rgenegamma;
2951
3181
         else {         /* sigma2 */
2952
3182
            pold = data.sigma2[i];
2953
3183
            y = log(pold);
2954
 
            ynew = y + finetune*rnd2NormalSym(m2NormalSym);
 
3184
            ynew = y + finetune*rndBactrian(mBactrian);
2955
3185
            ynew = reflect(ynew, yb[0], yb[1]);
2956
3186
            data.sigma2[i] = pnew = exp(ynew);
2957
3187
            gammaprior = data.sigma2gamma;
3009
3239
          jj = j-sptree.nspecies;
3010
3240
          x[jj] = (sptree.nodes[j].age - AgeLow[jj])/(sptree.nodes[sptree.nodes[j].father].age - AgeLow[jj]);
3011
3241
   }
3012
 
   lnacceptance = lnc = finetune*rnd2NormalSym(m2NormalSym);
 
3242
   lnacceptance = lnc = finetune*rndBactrian(mBactrian);
3013
3243
   c = exp(lnacceptance);
3014
3244
   sptree.nodes[sptree.root].age = AgeLow[0] + (sptree.nodes[sptree.root].age - AgeLow[0]) * c;
3015
3245
   for(j=sptree.nspecies; j<sptree.nnode; j++) {
3016
3246
      if(j==sptree.root) continue;
3017
 
          jj = j-sptree.nspecies;
3018
 
          tz = sptree.nodes[j].age;
3019
 
          sptree.nodes[j].age = AgeLow[jj] + x[jj] * (sptree.nodes[sptree.nodes[j].father].age - AgeLow[jj]);
3020
 
          lnacceptance += log( (sptree.nodes[j].age - AgeLow[jj])/(tz - AgeLow[jj]) );
 
3247
           jj = j-sptree.nspecies;
 
3248
           tz = sptree.nodes[j].age;
 
3249
           sptree.nodes[j].age = AgeLow[jj] + x[jj] * (sptree.nodes[sptree.nodes[j].father].age - AgeLow[jj]);
 
3250
           lnacceptance += log( (sptree.nodes[j].age - AgeLow[jj])/(tz - AgeLow[jj]) );
3021
3251
   }
3022
3252
 
3023
3253
   lnpTnew = lnpriorTimes();
3053
3283
   if(lnacceptance>0 || rndu()<exp(lnacceptance)) { /* accept */
3054
3284
      naccept = 1;
3055
3285
      data.lnpT = lnpTnew;
3056
 
      if(com.clock==3) data.lnpR = lnpRnew;
 
3286
      data.lnpR = lnpRnew;
3057
3287
      for(locus=0; locus<data.ngene; locus++)
3058
3288
             data.lnpDi[locus] = lnpDinew[locus];
3059
3289
      if(mcmc.usedata==1) switchconPin();
3103
3333
   if(com.TipDate)
3104
3334
      return mixingTipDate(lnL, finetune);
3105
3335
 
3106
 
   lnc = finetune*rnd2NormalSym(m2NormalSym);
 
3336
   lnc = finetune*rndBactrian(mBactrian);
3107
3337
   c = exp(lnc);
3108
3338
   for(j=sptree.nspecies; j<sptree.nnode; j++) 
3109
3339
      sptree.nodes[j].age *= c;
3154
3384
   double p = data.pfossilerror[0], q = data.pfossilerror[1];
3155
3385
 
3156
3386
   if(finetune<=0) error2("steplength = 0 in UpdatePFossilErrors");
3157
 
   pold = sptree.Pfossilerr;
3158
 
   pnew = pold + finetune*rnd2NormalSym(m2NormalSym);
3159
 
   sptree.Pfossilerr = pnew = reflect(pnew,0,1);
 
3387
   pold = data.Pfossilerr;
 
3388
   pnew = pold + finetune*rndBactrian(mBactrian);
 
3389
   data.Pfossilerr = pnew = reflect(pnew,0,1);
3160
3390
   lnacceptance = (p-1)*log(pnew/pold) + (q-1)*log((1-pnew)/(1-pold));
3161
3391
   lnpTnew = lnpriorTimes();
3162
 
   lnacceptance += lnpTnew-data.lnpT;
 
3392
   lnacceptance += lnpTnew - data.lnpT;
3163
3393
 
3164
3394
   if(lnacceptance>=0 || rndu()<exp(lnacceptance)) {
3165
3395
      naccept++;
3166
3396
      data.lnpT = lnpTnew;
3167
3397
   }
3168
3398
   else {
3169
 
      sptree.Pfossilerr = pold;
 
3399
      data.Pfossilerr = pold;
3170
3400
   }
3171
3401
   return(naccept);
3172
3402
}
3173
3403
 
3174
3404
 
3175
 
int DescriptiveStatisticsSimpleMCMCTREE (FILE *fout, char infile[], int nbin, int nrho)
 
3405
int DescriptiveStatisticsSimpleMCMCTREE (FILE *fout, char infile[], int nbin)
3176
3406
{
3177
3407
   FILE *fin=gfopen(infile,"r"), *fFigTree;
3178
3408
   char *fmt=" %9.6f", *fmt1=" %9.1f", timestr[32], FigTreef[96]="FigTree.tre";
3179
 
   double *x, *mean, *median, *minx, *maxx, *x005,*x995,*x025,*x975,*x25,*x75;
3180
 
   double *y;
 
3409
   double *x, *mean, *median, *minx, *maxx, *x005,*x995,*x025,*x975,*xHPD025,*xHPD975;
 
3410
   double *y, tmp[2];
3181
3411
   int  n, p, i,j, jj;
3182
3412
   int  lline=640000, ifields[MAXNFIELDS], Ignore1stColumn=1, ReadHeader=0;
3183
3413
   char *line;
3188
3418
   if(ReadHeader)
3189
3419
      for(i=0; i<p; i++) sscanf(line+ifields[i], "%s", varstr[i]);
3190
3420
 
3191
 
   x = (double*)malloc((p*13+p*p + p*nrho)*sizeof(double));
 
3421
   x = (double*)malloc((p*13+p*p)*sizeof(double));
3192
3422
   if (x==NULL) { puts("did not get enough space."); exit(-1); }
3193
 
   for(j=0;j<p*13+p*p + p*nrho; j++) x[j]=0;
 
3423
   for(j=0;j<p*13+p*p; j++) x[j]=0;
3194
3424
   mean=x+p; median=mean+p; minx=median+p; maxx=minx+p; 
3195
 
   x005=maxx+p; x995=x005+p; x025=x995+p; x975=x025+p; x25=x975+p; x75=x25+p;
 
3425
   x005=maxx+p; x995=x005+p; x025=x995+p; x975=x025+p; xHPD025=x975+p; xHPD975=xHPD025+p;
3196
3426
 
3197
3427
   for(i=0; i<n; i++) {
3198
3428
      for(j=0;j<p;j++) fscanf(fin,"%lf", &x[j]);
3214
3444
      else        median[jj]=y[(n+1)/2];
3215
3445
      x005[jj] = y[(int)(n*.005)];  x995[jj] = y[(int)(n*.995)];
3216
3446
      x025[jj] = y[(int)(n*.025)];  x975[jj] = y[(int)(n*.975)];
3217
 
      x25[jj]  = y[(int)(n*.25)];    x75[jj] = y[(int)(n*.75)];
 
3447
      HPDinterval(y, n, tmp, 0.05);
 
3448
      xHPD025[jj] = tmp[0];         xHPD975[jj] = tmp[1];
3218
3449
      printf("Summarizing... %d/%d%10s\r", jj+1, p, printtime(timestr));
3219
3450
   }
3220
3451
   printf("%40s\n", "");
3221
3452
 
3222
 
 
3223
3453
   for(j=sptree.nspecies; j<sptree.nnode; j++) {
3224
3454
      nodes[j].nodeStr = (char*)malloc(32*sizeof(char));
3225
3455
      jj = Ignore1stColumn + j - sptree.nspecies;
3263
3493
   printf("\n\nPosterior mean (95%% credibility interval)\n\n");
3264
3494
   for (j=Ignore1stColumn; j<p; j++) {
3265
3495
      printf("%-10s", varstr[j]);
3266
 
      printf("%7.4f (%6.4f, %6.4f)", mean[j], x025[j], x975[j]);
 
3496
      printf("%7.4f (%6.4f, %6.4f) (%6.4f, %6.4f)", mean[j], x025[j], x975[j], xHPD025[j], xHPD975[j]);
3267
3497
      if(j<Ignore1stColumn + sptree.nspecies-1) {
3268
3498
         printf("  (Jeffnode %2d)", 2*sptree.nspecies-1-1-j+Ignore1stColumn);
3269
3499
      
3270
 
         if(com.TipDate) printf(" time: %7.1f (%5.1f, %5.1f)", com.TipDate-mean[j]*com.TipDate_TimeUnit, 
3271
 
            com.TipDate-x975[j]*com.TipDate_TimeUnit, com.TipDate-x025[j]*com.TipDate_TimeUnit);
 
3500
         if(com.TipDate) 
 
3501
            printf(" time: %7.1f (%5.1f, %5.1f)", com.TipDate-mean[j]*com.TipDate_TimeUnit, 
 
3502
                  com.TipDate-x975[j]*com.TipDate_TimeUnit, com.TipDate-x025[j]*com.TipDate_TimeUnit);
3272
3503
      }
3273
3504
      printf("\n");
3274
3505
   }
3303
3534
   double Pjump[6]={0}, lnL=0, nround=0, *x, *mx, *vx, *px, postEFossil[MaxNFossils]={0};
3304
3535
   char timestr[36];
3305
3536
 
3306
 
   noisy = 2;
3307
3537
   mcmc.saveconP = 1;
3308
3538
   if(mcmc.usedata!=1) mcmc.saveconP = 0;
3309
3539
   for(j=0; j<sptree.nspecies*2-1; j++) 
3352
3582
 
3353
3583
/*
3354
3584
sptree.nodes[5].age=0.99;  sptree.nodes[6].age=0.8;  sptree.nodes[7].age=0.6;  sptree.nodes[8].age=0.55;     for(i=sptree.nspecies; i<sptree.nnode; i++) printf("%9.5f", sptree.nodes[i].age);
3355
 
printf("  lnpt = %9.6f\n", (lnp1=lnpriorTimesBDS_Approach2()));
 
3585
printf("  lnpt = %9.6f\n", (lnp1=lnpriorTimesBDS_Approach1()));
3356
3586
sptree.nodes[5].age=1.2;  sptree.nodes[6].age=0.7;  sptree.nodes[7].age=0.6;  sptree.nodes[8].age=0.55;     for(i=sptree.nspecies; i<sptree.nnode; i++) printf("%9.5f", sptree.nodes[i].age);
3357
 
printf("  lnpt = %9.6f", (lnp2=lnpriorTimesBDS_Approach2()));  
 
3587
printf("  lnpt = %9.6f", (lnp2=lnpriorTimesBDS_Approach1()));  
3358
3588
printf("  diff = %9.6f\n", lnp2-lnp1);
3359
3589
*/
3360
 
/* lnpriorTimesBDS_Approach2
3361
 
*/
3362
3590
 
3363
3591
exit(0);
3364
3592
}
3395
3623
   if(com.clock>1) printf("\tG(%.4f, %.4f) for sigma2\n", data.sigma2gamma[0], data.sigma2gamma[1]);
3396
3624
 
3397
3625
   /* calculates prior for times and likelihood for each locus */
 
3626
/*
 
3627
data.Pfossilerr=0.4;
 
3628
sptree.nodes[4].age = 1;
 
3629
sptree.nodes[5].age = 0.3;
 
3630
sptree.nodes[6].age = 0.7;
 
3631
printSptree();
 
3632
*/
3398
3633
   if(data.pfossilerror[0]) 
3399
3634
      SetupPriorTimesFossilErrors();
3400
3635
   data.lnpT = lnpriorTimes();
3411
3646
   printf("finetune steps (time rate mixing para RatePara ...):");
3412
3647
   for(j=0; j<nsteps; j++) 
3413
3648
      printf(" %7.4f",mcmc.finetune[j]);
3414
 
   if(com.clock==1) 
 
3649
   if(com.clock==1)
3415
3650
      printf("\n  paras: %d times, %d mu, (& kappa, alpha)\n", 
3416
3651
             sptree.nspecies-1, data.ngene);
3417
3652
   else
3433
3668
         zero(mx, com.np); 
3434
3669
         zero(vx, com.np); 
3435
3670
      }
 
3671
      nround++;
3436
3672
 
3437
 
      nround++;
3438
3673
      Pjump[0] = (Pjump[0]*(nround-1) + UpdateTimes(&lnL, mcmc.finetune[0]))/nround;
3439
3674
      Pjump[1] = (Pjump[1]*(nround-1) + UpdateRates(&lnL, mcmc.finetune[1]))/nround;
3440
3675
      Pjump[2] = (Pjump[2]*(nround-1) + mixing(&lnL, mcmc.finetune[2]))/nround;
3442
3677
         Pjump[3] = (Pjump[3]*(nround-1) + UpdateParameters(&lnL, mcmc.finetune[3]))/nround;
3443
3678
      if(com.clock>1)
3444
3679
         Pjump[4] = (Pjump[4]*(nround-1) + UpdateParaRates(mcmc.finetune[4],com.space))/nround;
3445
 
 
3446
3680
      if(data.pfossilerror[0])
3447
3681
         Pjump[5] = (Pjump[5]*(nround-1) + UpdatePFossilErrors(mcmc.finetune[5]))/nround;
3448
3682
 
 
3683
 
 
3684
      /* printf("root age = %.4f\n", sptree.nodes[sptree.root].age); */
 
3685
 
3449
3686
      collectx(fmcmc, x);
3450
3687
 
3451
3688
      for(j=0; j<com.np; j++) vx[j] += square(x[j] - mx[j]) * (nround-1.)/nround;
3460
3697
         if(mcmc.usedata) fprintf(fmcmc,"\t%.3f",lnL);
3461
3698
         FPN(fmcmc);
3462
3699
      }
3463
 
      if((ir+1)%max2(mcmc.sampfreq, mcmc.sampfreq*mcmc.nsample/10000)==0) {
 
3700
      if((ir+1)%max2(mcmc.sampfreq, mcmc.sampfreq*mcmc.nsample/100)==0) {
3464
3701
         printf("\r%3.0f%%", (ir+1.)/(mcmc.nsample*mcmc.sampfreq)*100.);
3465
3702
 
3466
3703
         for(j=0; j<nsteps; j++) 
3488
3725
 
3489
3726
         fflush(fmcmc);
3490
3727
         if(mcmc.print) fflush(fmcmc);
 
3728
 
3491
3729
      }
3492
3730
   }
3493
3731
   for(i=0; i<com.np; i++)
3517
3755
   FPN(fout); OutTreeN(fout,1,1); FPN(fout);
3518
3756
 
3519
3757
   if(mcmc.print)
3520
 
      DescriptiveStatisticsSimpleMCMCTREE(fout, com.mcmcf, 1, 1);
 
3758
      DescriptiveStatisticsSimpleMCMCTREE(fout, com.mcmcf, 1);
3521
3759
   if(data.pfossilerror[0]) {
3522
 
      free(sptree.CcomFossilErr);
 
3760
      free(data.CcomFossilErr);
3523
3761
      printf("\nPosterior probabilities that each fossil is in error.\n");
3524
3762
      for(i=k=0; i<sptree.nspecies*2-1; i++) {
3525
3763
         if( (j=sptree.nodes[i].fossil) )
3526
 
            printf("Node %2d: %s (%9.4f, %9.4f)%9.3f\n", i+1, fossils[j], 
 
3764
            printf("Node %2d: %s (%9.4f, %9.4f) %8.4f\n", i+1, fossils[j], 
3527
3765
                    sptree.nodes[i].pfossil[0], sptree.nodes[i].pfossil[1], postEFossil[k++]);
3528
3766
      }
3529
3767
      fprintf(fout, "\nPosterior probabilities that each fossil is in error.\n");
3530
3768
      for(i=k=0; i<sptree.nspecies*2-1; i++) {
3531
3769
         if( (j=sptree.nodes[i].fossil) )
3532
 
            fprintf(fout, "Node %2d: %s (%9.4f, %9.4f)%9.3f\n", i+1, fossils[j], 
 
3770
            fprintf(fout, "Node %2d: %s (%9.4f, %9.4f) %8.4f\n", i+1, fossils[j], 
3533
3771
                        sptree.nodes[i].pfossil[0], sptree.nodes[i].pfossil[1], postEFossil[k++]);
3534
3772
      }
3535
3773
   }