← Back to branch summary

~ubuntu-branches/ubuntu/utopic/paml/utopic

« back to all changes in this revision

Viewing changes to src/treesub.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/treesub.c
3
3
   such as baseml, basemlg, codeml, and pamp.
4
4
*/
5
5
 
6
 
extern char BASEs[], *EquateBASE[], AAs[], BINs[], CODONs[][4], nChara[], CharaMap[][64];
 
6
extern char BASEs[], *EquateBASE[], BASEs5[], *EquateBASE5[], AAs[], BINs[], CODONs[][4], nChara[], CharaMap[][64];
7
7
 
8
8
extern int noisy;
9
9
 
149
149
   int n31=(com.seqtype==CODONseq||com.seqtype==CODON2AAseq?3:1);
150
150
   int gap=(n31==3?3:10), nchar=(com.seqtype==AAseq?20:4);
151
151
   int h,b[3]={0};
152
 
   char *pch=((com.seqtype<=1||com.seqtype==CODON2AAseq)?BASEs:(com.seqtype==2?AAs:BINs));
 
152
   char *pch=((com.seqtype<=1||com.seqtype==CODON2AAseq) ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5 ? BASEs5 : BINs)));
153
153
   char str[4]="   ";
154
154
   double lst;
155
155
 
177
177
      if(com.spname[j]) free(com.spname[j]);
178
178
      com.spname[j] = (char*)malloc((lspname+1)*sizeof(char));
179
179
      for(i=0; i<lspname+1; i++) com.spname[j][i]=0;
180
 
      if((com.z[j]=(char*)realloc(com.z[j],com.ls*sizeof(char))) == NULL) 
 
180
      if((com.z[j] = (unsigned char*)realloc(com.z[j],com.ls*sizeof(unsigned char))) == NULL)
181
181
         error2("oom z");
182
182
   }
183
183
   com.rgene[0] = 1;   com.ngene = 1;  
259
259
            for (; i<com.ngene; i++)
260
260
               if (fscanf(fseq,"%d", &com.lgene[i])!=1) error2("EOF at lgene");
261
261
 
262
 
            if(!com.readpattern)
263
 
               for(i=0,k=0; i<com.ngene; k+=com.lgene[i],i++)
264
 
                  for(j=0; j<com.lgene[i]; j++)
265
 
                     com.pose[k+j]=i;
266
 
 
267
262
            for(i=0,k=0; i<com.ngene; i++) 
268
263
               k += com.lgene[i];
269
264
            if(k!=com.ls/n31) {
276
271
               puts("Sequence length in number of nucleotides.");
277
272
               exit(-1);
278
273
            }
 
274
            if(!com.readpattern)
 
275
               for(i=0,k=0; i<com.ngene; k+=com.lgene[i],i++)
 
276
                  for(j=0; j<com.lgene[i]; j++)
 
277
                     com.pose[k+j] = i;
 
278
 
279
279
         }
280
280
         else {                   /* site marks on later line(s)  */
281
281
            if(com.readpattern) 
304
304
   if (Sequential)  {    /* sequential */
305
305
      if (noisy) printf ("Reading sequences, sequential format..\n");
306
306
      for (j=0; j<com.ns; j++) {
307
 
         lspname=LSPNAME;
 
307
         lspname = LSPNAME;
308
308
         for (i=0; i<2*lspname; i++) line[i]='\0';
309
309
         if (!fgets (line, lline, fseq)) error2("EOF?");
310
310
         if (blankline(line)) {
322
322
            if (!isgraph(com.spname[j][k]))   com.spname[j][k]=0;
323
323
            else    break;
324
324
 
325
 
         if (noisy>=2) printf ("Reading seq #%2d: %s     \n", j+1, com.spname[j]);
 
325
         if (noisy>=2) printf ("Reading seq #%2d: %s     \r", j+1, com.spname[j]);
326
326
         for (k=0; k<com.ls; p++) {
327
327
            while (*p=='\n' || *p=='\0') {
328
328
               p=fgets(line, lline, fseq);
343
343
               com.z[j][k++] = *p;
344
344
            else if (isalpha(*p)) {
345
345
               printf("\nError in sequence data file: %c at %d seq %d.\n",*p,k+1,j+1); 
346
 
               puts("Perhaps you did not separate the sequence from its name by >2 spaces?");
 
346
               puts("Make sure to separate the sequence from its name by 2 or more spaces.");
347
347
               exit(0); 
348
348
            }
349
349
            else if (*p == (char)EOF) error2("EOF?");
380
380
            }
381
381
            p=line;
382
382
            if (igroup==0) {
383
 
               lspname=LSPNAME;
 
383
               lspname = LSPNAME;
384
384
               while(isspace(*p)) p++;
385
 
               if ((ch=strstr(p,"  ")-p)<lspname && ch>0) lspname=ch;
 
385
               if ((ch=strstr(p,"  ")-p)<lspname && ch>0)
 
386
                  lspname = ch;
386
387
               strncpy (com.spname[j], p, lspname);
387
 
               k=strlen(com.spname[j]);
388
 
               p+=(k<lspname?k:lspname);
 
388
               k = strlen(com.spname[j]);
 
389
               p += (k<lspname?k:lspname);
389
390
 
390
391
               for (; k>0; k--)   /* trim spaces */
391
 
                  if (!isgraph(com.spname[j][k]))  com.spname[j][k]=0;
392
 
                  else   break;
 
392
                  if (!isgraph(com.spname[j][k]))
 
393
                     com.spname[j][k]=0;
 
394
                  else
 
395
                     break;
393
396
               if(noisy>=2) printf("Reading seq #%2d: %s     \r",j+1,com.spname[j]);
394
397
            }
395
398
            for (; *p && *p!='\n'; p++) {
396
399
               if (lt[j]==com.ls) break;
397
 
               *p=(char)toupper(*p);
 
400
               *p = (char)toupper(*p);
398
401
               if((com.seqtype==BASEseq || com.seqtype==CODONseq) && *p=='U') 
399
 
                  *p='T';
400
 
               p1=strchr(pch,*p);
 
402
                  *p = 'T';
 
403
               p1 = strchr(pch, *p);
401
404
               if (p1 && p1-pch>=nchar) 
402
405
                  miss = 1;
403
406
               if (*p == eq) {
411
414
               else if (p1)
412
415
                  com.z[j][lt[j]++]=*p;
413
416
               else if (isalpha(*p)) {
414
 
                  printf("\nerr:%c at %d seq %d block %d.",
 
417
                  printf("\nerr: unrecognised character %c at %d seq %d block %d.",
415
418
                          *p,lt[j]+1,j+1,igroup+1);
416
419
                  exit(-1);
417
420
               }
428
431
   }
429
432
   free(line);
430
433
 
 
434
#ifdef CODEML
 
435
   if(com.seqtype==1) MarkStopCodons();
 
436
#endif
 
437
 
431
438
   if(!miss)
432
439
      com.cleandata = 1;
433
440
   else if (cleandata) {  /* forced removal of ambiguity characters */
474
481
   if(com.ngene>1 && com.Mgene==1 && com.verbose)  printSeqsMgenes ();
475
482
 
476
483
   if(com.bootstrap) { BootstrapSeq("boot.txt");  exit(0); }
 
484
   /* mask stop codons as ???.  */
477
485
#endif
478
486
 
479
487
   if(noisy>=2) printf ("\nSequences read..\n");
533
541
   return (0);
534
542
}
535
543
 
 
544
 
 
545
#if (defined CODEML)
 
546
 
 
547
void MarkStopCodons(void)
 
548
{
 
549
/* this converts the whole column into ??? if there is a stop codon in one sequence.
 
550
   Data in com.z[] are just read in and not encoded yet.
 
551
*/
 
552
   int i,j,h,k, fixed=0;
 
553
   char codon[4]="", stops[10][4]={"","",""}, nstops=0;
 
554
 
 
555
   if(com.seqtype!=1) error2("should not be here");
 
556
 
 
557
   for(i=0; i<64; i++) 
 
558
      if(GeneticCode[com.icode][i]==-1) 
 
559
         getcodon(stops[nstops++], i);
 
560
 
 
561
   for(h=0; h<com.ls/3; h++) {
 
562
      for(i=0; i<com.ns; i++) {
 
563
         codon[0] = com.z[i][h*3+0];
 
564
         codon[1] = com.z[i][h*3+1];
 
565
         codon[2] = com.z[i][h*3+2];
 
566
         for(j=0; j<nstops; j++) 
 
567
            if(strcmp(codon, stops[j])==0) break;
 
568
         if(j<nstops) break;
 
569
      }
 
570
      if(i<com.ns) {
 
571
         for(i=0; i<com.ns; i++) 
 
572
            com.z[i][h*3+0] = com.z[i][h*3+1] = com.z[i][h*3+2] = '?';
 
573
         fixed++;
 
574
      }
 
575
   }
 
576
   if(fixed) {
 
577
      printf("\n%2d columns are converted into ??? because of stop codons\nPress Enter to continue", fixed);
 
578
      getchar();
 
579
   }
 
580
}
 
581
 
 
582
#endif
 
583
 
 
584
 
536
585
void RemoveEmptySequences(void)
537
586
{
538
587
/* this removes empty sequences (? or - only) and adjust com.ns
569
618
 
570
619
int printPatterns(FILE *fout)
571
620
{
572
 
   int j,h, n31=(com.seqtype==CODONseq||com.seqtype==CODON2AAseq?3:1);
573
 
   int gap=(n31==3?3:10), nchar=(com.seqtype==AAseq?20:4);
 
621
   int j,h, n31 = (com.seqtype==CODONseq||com.seqtype==CODON2AAseq ? 3 : 1);
 
622
   int gap=(n31==3?3:10), n=(com.seqtype==AAseq?20:4);
574
623
 
575
624
   fprintf(fout,"\n%10d %10d  P", com.ns, com.npatt*n31);
576
625
   if(com.ngene>1) {
592
641
       fprintf(fout, "\n");
593
642
       for(h=0; h<com.npatt; h++) {
594
643
          fprintf(fout," %4.0f", com.fpatt[h]);
595
 
          if((h+1)%15==0) FPN(fout);
 
644
          if((h+1)%15 == 0) FPN(fout);
596
645
       }
597
646
       fprintf(fout, "\n\n");
598
647
   }
605
654
{
606
655
/* This encodes sequences and set up com.TipMap[][], called after sites are collapsed 
607
656
   into patterns.
608
 
   
609
 
   For codonml, codons are coded into 0, 1, ..., 60 for the universal code.
610
 
   For    yn00, codons are coded into 0, 1, ..., 63 for the universal code.
611
 
   This does not look like a good idea, and perhaps should be changed.
612
657
*/
613
658
   int n=com.ncode, nA, is,h, i, j, k,ic, indel=0, ch, b[3];
614
 
   char *pchar = ((com.seqtype==0||com.seqtype==1) ? BASEs : (com.seqtype==2 ? AAs : BINs));
 
659
   char *pch = ((com.seqtype==0||com.seqtype==1) ? BASEs : (com.seqtype==2 ? AAs : (com.seqtype==5 ? BASEs5: BINs)));
615
660
   unsigned char c[4]="", str[4]="   ";
616
661
 
617
 
   if(com.seqtype==0 || com.seqtype==2) {
 
662
   if(com.seqtype != 1) {
618
663
      for(is=0; is<com.ns; is++) {
619
664
         for (h=0; h<com.npatt; h++) {
620
665
            ch = com.z[is][h];
621
 
            k = strchr(pchar, ch) - pchar;
 
666
            com.z[is][h] = (char)(k = strchr(pch, ch) - pch);
622
667
            if(k<0) {
623
668
               printf("strange character %c in seq %d site %d\n", ch, is+1, h+1);
624
669
               exit(-1);
625
670
            }
626
 
            com.z[is][h] = k;
627
671
         }
628
672
      }
629
673
   }
643
687
               if(b[k]<0) printf("strange nucleotide %c in seq %d\n", c[k], j+1);
644
688
            }
645
689
            if(b[0]<4 && b[1]<4 && b[2]<4) {
646
 
               k = FROM64[b[0]*16+b[1]*4+b[2]];
 
690
               k = FROM64[b[0]*16 + b[1]*4 + b[2]];
647
691
               if(k<0) {
648
 
                  printf("\nstop codon %s in seq #%2d: %s\n", c, j+1,com.spname[j]);
 
692
                  printf("\nstop codon %s in seq #%2d: %s\n", c, j+1, com.spname[j]);
 
693
                  printf("\ncodons in other sequences are\n");
 
694
                  for(i=0; i<com.ns; i++) {
 
695
                     for(k=0; k<3; k++) c[k] = com.z[i][h*3+k]; 
 
696
                     printf("seq #%2d %-30s %s\n", i+1, com.spname[i], c);
 
697
                  }
649
698
                  exit(-1);
650
699
               }
651
700
            }
677
726
/* This sets up CharaMap, the map from the ambiguity characters to resolved characters.
678
727
*/
679
728
   int n=com.ncode, i,j, i0,i1,i2, nb[3], ib[3][4], ic;
680
 
   char *pchar = (com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs : BINs));
 
729
   char *pch = (com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs : (com.seqtype==5 ? BASEs5: BINs)));
 
730
   char *pbases = (com.seqtype==0 ? BASEs : (com.seqtype==5 ? BASEs5: NULL));
 
731
   char **pEquateBASE = (com.seqtype==0 ? EquateBASE : (com.seqtype==5 ? EquateBASE5 : NULL));
 
732
   char debug=0;
681
733
 
682
 
   for(j=0; j<n; j++) {
683
 
      nChara[j] = 1;
684
 
      CharaMap[j][0] = j;
 
734
   for(j=0; j<n; j++) {  /* basic characters, coded according to the definition in pch. */
 
735
      nChara[j] = (char)1;
 
736
      CharaMap[j][0] = (char)j;
685
737
   }
686
738
 
687
 
   if(com.seqtype==0 || com.seqtype==2) {
688
 
      for(j=n,pchar+=n; *pchar; j++,pchar++) {
689
 
         if(com.seqtype==0) {
690
 
            nChara[j] = strlen(EquateBASE[j]);
 
739
   if(com.seqtype != 1) {
 
740
      for(j=n,pch+=n; *pch; j++,pch++) {
 
741
         if(com.seqtype==0 || com.seqtype==5) {  /* ambiguities are allowed for those 2 types */
 
742
            nChara[j] = (char)strlen(pEquateBASE[j]);
691
743
            for(i=0; i<nChara[j]; i++)
692
 
               CharaMap[j][i] = (char)(strchr(BASEs, EquateBASE[j][i]) - BASEs);
 
744
               CharaMap[j][i] = (char)(strchr(pbases, pEquateBASE[j][i]) - pbases);
693
745
         }
694
 
         else {
695
 
            nChara[j] = n;
 
746
         else {  /* for non-nucleotide characters, ambiguity characters must be ? or -. */
 
747
            nChara[j] = (char)n;
696
748
            for(i=0; i<n; i++)
697
 
               CharaMap[j][i] = i;
 
749
               CharaMap[j][i] = (char)i;
 
750
         }
 
751
         if(debug) {
 
752
            printf("character %c (%d): ", pbases[j], nChara[j]);
 
753
            for(i=0; i<nChara[j]; i++)
 
754
               printf("%c", pbases[CharaMap[j][i]]);
 
755
            printf("\n");
698
756
         }
699
757
      }
700
758
   }
701
759
#ifdef CODEML
702
 
   else if(com.seqtype==1) {
 
760
   else {
703
761
      for(j=n; j<256 && CODONs[j][0]; j++) {
704
 
         nChara[j] = 0;
 
762
         nChara[j] = (char)0;
705
763
         for(i=0; i<3; i++)
706
764
            NucListall(CODONs[j][i], &nb[i], ib[i]);
707
765
         for(i0=0; i0<nb[0]; i0++) {
777
835
   analysis.
778
836
*/
779
837
   int j, h, it, ic;
780
 
   char codon[4]="   ", b[3];
 
838
   char codon[4]="   ";
781
839
   int n31=(com.seqtype==CODONseq||com.seqtype==CODON2AAseq?3:1);
782
840
   int gap=(n31==3?3:10);
783
841
 
999
1057
   This routine is called by baseml and aaml.  codonml uses another
1000
1058
   routine InitializeCodon()
1001
1059
*/
1002
 
   char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs)), indel[]="-?";
1003
 
   int wname=30, h,js,k, ig, nconstp, nc=com.ncode;
 
1060
   char *pch = (com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs : (com.seqtype==5 ? BASEs5: BINs)));
 
1061
   char indel[]="-?";
 
1062
   int wname=30, h,js,k, ig, nconstp, n=com.ncode;
1004
1063
   int irf, nrf=20;
1005
1064
   double pi0[20], t,lmax=0, X2G[2], *pisg;  /* freq for species & gene, for X2 & G */
1006
1065
 
1007
1066
   if(noisy) printf("Counting frequencies..");
1008
1067
   if(fout)  fprintf(fout,"\nFrequencies..");
1009
 
   if((pisg=(double*)malloc(com.ns*nc*sizeof(double))) == NULL)
 
1068
   if((pisg=(double*)malloc(com.ns*n*sizeof(double))) == NULL)
1010
1069
      error2("oom pisg");
1011
1070
   for(h=0,nconstp=0; h<com.npatt; h++) {
1012
1071
      for (js=1; js<com.ns; js++)
1014
1073
      if (js==com.ns && com.z[0][h]!=indel[0] && com.z[0][h]!=indel[1])
1015
1074
         nconstp += (int)com.fpatt[h];
1016
1075
   }
1017
 
   for (ig=0,zero(com.pi,nc); ig<com.ngene; ig++) {
 
1076
   for (ig=0,zero(com.pi,n); ig<com.ngene; ig++) {
1018
1077
      if (com.ngene>1)
1019
1078
         fprintf (fout,"\n\nGene %2d (len %4d)", ig+1, com.lgene[ig]-(ig==0?0:com.lgene[ig-1]));
1020
1079
      fprintf(fout,"\n%*s",wname, "");
1021
 
      for(k=0; k<nc; k++) fprintf(fout,"%7c", pch[k]);
 
1080
      for(k=0; k<n; k++) fprintf(fout,"%7c", pch[k]);
1022
1081
 
1023
1082
      /* The following block calculates freqs in each species for each gene.  
1024
1083
         Ambiguities are resolved in each species.  com.pi and com.piG are 
1025
1084
         used for output only, and are not be used later with missing data.
1026
1085
      */
1027
 
      zero(com.piG[ig], nc);
1028
 
      zero(pisg, com.ns*nc);
 
1086
      zero(com.piG[ig], n);
 
1087
      zero(pisg, com.ns*n);
1029
1088
      for(js=0; js<com.ns; js++) {
1030
 
         fillxc(pi0, 1.0/nc, nc);
 
1089
         fillxc(pi0, 1.0/n, n);
1031
1090
         for(irf=0; irf<nrf; irf++) {
1032
 
            zero(com.pi, nc);
 
1091
            zero(com.pi, n);
1033
1092
            AddFreqSeqGene(js, ig, pi0, com.pi);
1034
 
            t = sum(com.pi, nc);
 
1093
            t = sum(com.pi, n);
1035
1094
            if(t<1e-10) {
1036
1095
               printf("Some sequences are empty.\n");
1037
 
               fillxc(com.pi, 1.0/nc, nc);
 
1096
               fillxc(com.pi, 1.0/n, n);
1038
1097
            }
1039
1098
            else 
1040
 
               abyx(1/t, com.pi, nc);
1041
 
            if(com.cleandata || com.cleandata || (t=distance(com.pi,pi0,nc))<1e-8)
 
1099
               abyx(1/t, com.pi, n);
 
1100
            if(com.cleandata || com.cleandata || (t=distance(com.pi,pi0,n))<1e-8)
1042
1101
               break;
1043
 
            xtoy(com.pi, pi0, nc);
 
1102
            xtoy(com.pi, pi0, n);
1044
1103
         }   /* for(irf) */
1045
1104
         fprintf(fout,"\n%-*s", wname, com.spname[js]);
1046
 
         for(k=0; k<nc; k++) fprintf(fout, "%7.4f", com.pi[k]);
1047
 
         for(k=0; k<nc; k++) com.piG[ig][k] += com.pi[k]/com.ns;
1048
 
         xtoy(com.pi, pisg+js*nc, nc);
 
1105
         for(k=0; k<n; k++) fprintf(fout, "%8.5f", com.pi[k]);
 
1106
         for(k=0; k<n; k++) com.piG[ig][k] += com.pi[k]/com.ns;
 
1107
         xtoy(com.pi, pisg+js*n, n);
1049
1108
      }    /* for(js,ns) */
1050
1109
      if(com.ngene>1) {
1051
1110
         fprintf(fout,"\n\n%-*s",wname,"Mean");
1052
 
         for(k=0; k<nc; k++) fprintf(fout,"%7.4f",com.piG[ig][k]);
 
1111
         for(k=0; k<n; k++) fprintf(fout,"%7.4f",com.piG[ig][k]);
1053
1112
      }
1054
1113
 
1055
 
      Chi2FreqHomo(pisg, com.ns, nc, X2G);
 
1114
      Chi2FreqHomo(pisg, com.ns, n, X2G);
1056
1115
 
1057
1116
      fprintf(fout,"\n\nHomogeneity statistic: X2 = %.5f G = %.5f ",X2G[0], X2G[1]);
1058
1117
 
1067
1126
      each gene for com.piG[]).  These are used in ML calculation later.
1068
1127
   */
1069
1128
   if(com.cleandata) {
1070
 
      for (ig=0,zero(com.pi,nc); ig<com.ngene; ig++) {
 
1129
      for (ig=0,zero(com.pi,n); ig<com.ngene; ig++) {
1071
1130
         t = (ig==0 ? com.lgene[0] : com.lgene[ig]-com.lgene[ig-1])/(double)com.ls;
1072
 
         for(k=0; k<nc; k++)  com.pi[k] += com.piG[ig][k]*t;
 
1131
         for(k=0; k<n; k++)  com.pi[k] += com.piG[ig][k]*t;
1073
1132
      }
1074
1133
   }
1075
1134
   else {
1076
1135
      for (ig=0; ig<com.ngene; ig++) { 
1077
 
         xtoy(com.piG[ig], pi0, nc);
 
1136
         xtoy(com.piG[ig], pi0, n);
1078
1137
         for(irf=0; irf<nrf; irf++) {  /* com.piG[] */
1079
 
            zero(com.piG[ig], nc);
 
1138
            zero(com.piG[ig], n);
1080
1139
            for(js=0; js<com.ns; js++)
1081
1140
               AddFreqSeqGene(js, ig, pi0, com.piG[ig]);
1082
 
            t = sum(com.piG[ig], nc);
 
1141
            t = sum(com.piG[ig], n);
1083
1142
            if(t<1e-10) 
1084
1143
               puts("empty sequences?");
1085
 
            abyx(1/t, com.piG[ig], nc);
1086
 
            if(distance(com.piG[ig], pi0, nc)<1e-8) break;
1087
 
            xtoy(com.piG[ig], pi0, nc);
 
1144
            abyx(1/t, com.piG[ig], n);
 
1145
            if(distance(com.piG[ig], pi0, n)<1e-8) break;
 
1146
            xtoy(com.piG[ig], pi0, n);
1088
1147
         }         /* for(irf) */
1089
1148
      }            /* for(ig) */
1090
 
      zero(pi0, nc);
1091
 
      for(k=0; k<nc; k++) for(ig=0; ig<com.ngene; ig++) 
 
1149
      zero(pi0, n);
 
1150
      for(k=0; k<n; k++) for(ig=0; ig<com.ngene; ig++) 
1092
1151
         pi0[k] += com.piG[ig][k]/com.ngene;
1093
1152
      for(irf=0; irf<nrf; irf++) {  /* com.pi[] */
1094
 
         zero(com.pi,nc);
 
1153
         zero(com.pi,n);
1095
1154
         for(ig=0; ig<com.ngene; ig++)  for(js=0; js<com.ns; js++)
1096
1155
            AddFreqSeqGene(js, ig, pi0, com.pi);
1097
 
         abyx(1/sum(com.pi,nc), com.pi, nc);
1098
 
         if(distance(com.pi, pi0, nc)<1e-8) break;
1099
 
         xtoy(com.pi,pi0,nc);
 
1156
         abyx(1/sum(com.pi,n), com.pi, n);
 
1157
         if(distance(com.pi, pi0, n)<1e-8) break;
 
1158
         xtoy(com.pi, pi0, n);
1100
1159
      }            /* for(ig) */
1101
1160
   }
1102
1161
   fprintf (fout, "\n\n%-*s", wname, "Average");
1103
 
   for(k=0; k<nc; k++) fprintf(fout," %7.4f", com.pi[k]);
 
1162
   for(k=0; k<n; k++) fprintf(fout,"%8.5f", com.pi[k]);
1104
1163
   if(!com.cleandata) fputs("\n(Ambiguity characters are used to calculate freqs.)\n",fout);
1105
1164
 
1106
1165
   fprintf (fout,"\n\n# constant sites: %6d (%.2f%%)",
1107
1166
            nconstp, (double)nconstp*100./com.ls);
1108
1167
 
1109
1168
   if (com.model==0 || (com.seqtype==BASEseq && com.model==1)) {
1110
 
      fillxc(com.pi, 1./nc, nc);
1111
 
      FOR(ig,com.ngene) xtoy (com.pi, com.piG[ig], nc);
 
1169
      fillxc(com.pi, 1./n, n);
 
1170
      for(ig=0; ig<com.ngene; ig++)
 
1171
         xtoy(com.pi, com.piG[ig], n);
1112
1172
   }
1113
1173
   if (com.seqtype==BASEseq && com.model==5) { /* T92 model */
1114
 
      com.pi[0]=com.pi[2]=(com.pi[0]+com.pi[2])/2;
1115
 
      com.pi[1]=com.pi[3]=(com.pi[1]+com.pi[3])/2;
 
1174
      com.pi[0] = com.pi[2] = (com.pi[0] + com.pi[2])/2;
 
1175
      com.pi[1] = com.pi[3] = (com.pi[1] + com.pi[3])/2;
1116
1176
      for(ig=0; ig<com.ngene; ig++) {
1117
1177
         com.piG[ig][0] = com.piG[ig][2] = (com.piG[ig][0] + com.piG[ig][2])/2;
1118
1178
         com.piG[ig][1] = com.piG[ig][3] = (com.piG[ig][1] + com.piG[ig][3])/2;
1121
1181
 
1122
1182
   /* this is used only for REV & REVu in baseml and model==3 in aaml */
1123
1183
   if(com.seqtype==AAseq) {
1124
 
      for (k=0,t=0; k<nc; k++) t+=(com.pi[k]>0);
 
1184
      for (k=0,t=0; k<n; k++)  t += (com.pi[k]>0);
1125
1185
      if (t<=4)
1126
1186
         puts("\n\a\t\tAre these a.a. sequences?");
1127
1187
   }
1128
1188
   if(com.cleandata && com.ngene==1) {
1129
1189
      for(h=0,lmax=-(double)com.ls*log((double)com.ls); h<com.npatt; h++)
1130
 
         if(com.fpatt[h]>1) lmax+=com.fpatt[h]*log((double)com.fpatt[h]);
 
1190
         if(com.fpatt[h]>1) lmax += com.fpatt[h]*log((double)com.fpatt[h]);
1131
1191
   }
1132
1192
   if(fout) {
1133
1193
      if(lmax) fprintf(fout, "\nln Lmax (unconstrained) = %.6f\n", lmax);
1145
1205
   using pi0, by resolving ambiguities.  The data are coded.  com.cleandata==1 or 0.
1146
1206
   This is for nucleotide and amino acid sequences only.
1147
1207
*/
1148
 
   char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs));
1149
 
   int k, h, b, nc=com.ncode;
 
1208
   char *pch = (com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs : (com.seqtype==5 ? BASEs5: BINs)));
 
1209
   int k, h, b, n=com.ncode;
1150
1210
   double t;
1151
1211
 
1152
1212
   if(com.cleandata) {
1156
1216
   else {
1157
1217
      for(h=com.posG[ig]; h<com.posG[ig+1]; h++) {
1158
1218
         b = com.z[js][h];
1159
 
         if(b<nc)
 
1219
         if(b<n)
1160
1220
            pi[b] += com.fpatt[h];
1161
1221
         else {
1162
1222
            /*
1174
1234
                  pi[CharaMap[b][k]] += pi0[CharaMap[b][k]]/t * com.fpatt[h];
1175
1235
            }
1176
1236
            else if(com.seqtype==AAseq)  /* unrecognized AAs are treated as "?". */
1177
 
               for(k=0; k<nc; k++) pi[k] += pi0[k]*com.fpatt[h];
 
1237
               for(k=0; k<n; k++) pi[k] += pi0[k]*com.fpatt[h];
1178
1238
         }
1179
1239
      }
1180
1240
   }
1190
1250
   All characters in com.z[][] not found in the character string pch are
1191
1251
   considered ambiguity characters and are removed.
1192
1252
*/
1193
 
   int  h,k, j,js,lnew,nindel, n31,nchar;
1194
 
   char b, *pch, *miss;  /* miss[h]=1 if site (codon) h is missing, 0 otherwise */
 
1253
   int  n=com.ncode, h,k, j,js,lnew,nindel, n31=1;
 
1254
   char b, *miss;  /* miss[h]=1 if site (codon) h is missing, 0 otherwise */
 
1255
   char *pch=((com.seqtype<=1||com.seqtype==CODON2AAseq)?BASEs:(com.seqtype==2?AAs: (com.seqtype==5?BASEs5:BINs)));
1195
1256
 
1196
 
   if(com.seqtype==CODONseq||com.seqtype==CODON2AAseq)
1197
 
      { n31=3; nchar=4; pch=BASEs; }
1198
 
   else {
1199
 
      n31=1;
1200
 
      if(com.seqtype==AAseq)        { nchar=20; pch=AAs; }
1201
 
      else if(com.seqtype==BASEseq) { nchar=4; pch=BASEs; }
1202
 
      else                          { nchar=2; pch=BINs; }
1203
 
    }
 
1257
   if(com.seqtype==CODONseq || com.seqtype==CODON2AAseq) {
 
1258
      n31=3; n=4;
 
1259
   }
1204
1260
 
1205
1261
   if (com.ls%n31) error2("ls in RemoveIndel.");
1206
1262
   if((miss=(char*)malloc(com.ls/n31 *sizeof(char)))==NULL)
1207
1263
      error2("oom miss");
1208
 
   FOR (h,com.ls/n31) miss[h]=0;
 
1264
   for(h=0; h<com.ls/n31; h++)
 
1265
      miss[h] = 0;
1209
1266
   for (js=0; js<com.ns; js++) {
1210
1267
      for (h=0,nindel=0; h<com.ls/n31; h++) {
1211
1268
         for (k=0; k<n31; k++) {
1212
 
            b=(char)toupper(com.z[js][h*n31+k]);
1213
 
            FOR(j,nchar) if(b==pch[j]) break;
1214
 
            if(j==nchar) { miss[h]=1; nindel++; }
 
1269
            b = (char)toupper(com.z[js][h*n31+k]);
 
1270
            for(j=0; j<n; j++) 
 
1271
               if(b==pch[j]) break;
 
1272
            if(j==n) {
 
1273
               miss[h]=1; nindel++; 
 
1274
            }
1215
1275
         }
1216
1276
      }
1217
1277
      if (noisy>2 && nindel) 
1247
1307
   Uses transformed sequences.  
1248
1308
   Not used for a long time.  Does not work if com.pose is NULL.  
1249
1309
*/
1250
 
   char *imark, *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs));
 
1310
   char *imark;
 
1311
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
1251
1312
   int h, i, markb[NS], inf, lsinf;
1252
1313
   FILE *finf, *fninf;
1253
1314
 
1306
1367
   If com.pose is not NULL, the routine also updates com.pose.  This allows 
1307
1368
   the program to work if com.readpattern==1.
1308
1369
*/
1309
 
   char zh[NS], b, gap, *pch=(com.seqtype==0 ? BASEs : AAs);
 
1370
   char zh[NS], b, gap;
 
1371
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
1310
1372
   int npatt0=com.npatt, h, ht, j,k, same=0, ig, recode;
1311
1373
 
1312
1374
   if(com.seqtype==1) 
1313
1375
      error2("PatternWeightJC69like does not work for codon seqs");
1314
1376
   if(noisy) printf("Counting site patterns again, for JC69.\n");
1315
 
   gap = strchr(pch, (int)'-') - pch;
 
1377
   gap = (char) (strchr(pch, (int)'-') - pch);
1316
1378
   for (h=0,com.npatt=0,ig=-1; h<npatt0; h++) {
1317
1379
      if (ig<com.ngene-1 && h==com.posG[ig+1])
1318
1380
         com.posG[++ig] = com.npatt; 
1374
1436
            if(com.pose[k]==h) com.pose[k] = ht;
1375
1437
   }     /* for (h)   */
1376
1438
   com.posG[com.ngene] = com.npatt;
1377
 
   if (noisy) printf ("\nnew no. site patterns:%7d\n", com.npatt);
 
1439
   if (noisy) printf ("new no. site patterns:%7d\n", com.npatt);
1378
1440
 
1379
1441
   if(fout) {
1380
1442
      fprintf(fout, "\n\nPrinting out site pattern counts\n");
1408
1470
   This uses com.seqtype.
1409
1471
*/
1410
1472
   int h, hp, gap=10;
1411
 
   char *pch = (com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs)), str[4]="";
 
1473
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
 
1474
   char str[4]="";
1412
1475
   int nb = (com.seqtype==CODONseq?3:1);
1413
1476
 
1414
1477
   for(h=0; h<ls; h++) {
1600
1663
   double x[16], sumx, larged=9;
1601
1664
 
1602
1665
   zero(x, 16);
1603
 
   if(com.cleandata) {
 
1666
   if(com.cleandata && com.seqtype==0) {
1604
1667
      for (h=0; h<com.npatt; h++)
1605
1668
         x[com.z[is][h]*n+com.z[js][h]] += com.fpatt[h];
1606
1669
   }
1692
1755
}
1693
1756
 
1694
1757
 
1695
 
int EigenTN93 (int model, double kappa1, double kappa2, double pi[],
 
1758
int eigenTN93 (int model, double kappa1, double kappa2, double pi[],
1696
1759
    int *nR, double Root[], double Cijk[])
1697
1760
{
1698
1761
/* initialize Cijk[] & Root[], which are the only part to be changed
1720
1783
   V[2*4+0]=V[2*4+1]=0;  V[2*4+2]=1;   V[2*4+3]=-1;
1721
1784
   V[3*4+0]=1;  V[3*4+1]=-1;   V[3*4+2]=V[3*4+3]=0;
1722
1785
 
1723
 
   FOR (i,4) FOR (j,4) {
 
1786
   for(i=0; i<4; i++) for(j=0; j<4; j++) {
1724
1787
      Cijk[i*4*nr+j*nr+0]=U[i*4+0]*V[0*4+j];
1725
1788
      switch (model) {
1726
1789
      case JC69:
1738
1801
         for (k=1; k<4; k++)  Cijk[i*4*nr+j*nr+k] = U[i*4+k]*V[k*4+j];
1739
1802
         break;
1740
1803
      default:
1741
 
         error2("model in EigenTN93");
 
1804
         error2("model in eigenTN93");
1742
1805
      }
1743
1806
   }
1744
1807
#ifdef BASEMLG
1768
1831
}
1769
1832
 
1770
1833
 
1771
 
int printsmaCodon (FILE *fout,char * z[],int ns,int ls,int lline,int simple)
 
1834
int printsmaCodon (FILE *fout, unsigned char * z[],int ns,int ls,int lline,int simple)
1772
1835
{
1773
1836
/* print, in blocks, multiple aligned and transformed codon sequences.
1774
1837
   indels removed.
1775
1838
   This is needed as codons are coded 0,1, 2, ..., 60, and 
1776
1839
   printsma won't work.
1777
1840
*/
1778
 
   int ig, ngroup, lt, il,is, i,b, lspname=20;
 
1841
   int ig, ngroup, lt, il,is, i,b, lspname=30;
1779
1842
   char equal='.',*pz, c0[4],c[4];
1780
1843
 
1781
1844
   if(ls==0) return(1);
1783
1846
   for (ig=0,FPN(fout); ig<ngroup; ig++)  {
1784
1847
      /* fprintf (fout,"%-8d\n", ig*lline+1); */
1785
1848
      for (is=0; is<ns; is++) {
1786
 
         fprintf(fout,"%-*s  ", lspname,com.spname[is]);
 
1849
         fprintf(fout,"%-*s  ", lspname, com.spname[is]);
1787
1850
         lt=0; 
1788
1851
         for(il=ig*lline,pz=z[is]+il; lt<lline && il<ls; il++,lt++,pz++) {
1789
1852
            b = *pz;  
1925
1988
         if(dN==1) dN = 0;
1926
1989
         else      dN = (dN<=0 ? -1 : 3./4*(alpha==0?-log(dN):alpha*(pow(dN,-1/alpha)-1)));
1927
1990
 
1928
 
         dN_dS = (dS>0 ? dN/dS : -1);
 
1991
         dN_dS = (dS>0 && dN>0 ? dN/dS : -1);
1929
1992
         if(fout) fprintf(fout,"%7.4f (%5.4f %5.4f)",   dN_dS, dN, dS);
1930
1993
 
1931
1994
         if(N>0 && dN<0)  dN = bigD; 
1961
2024
 
1962
2025
#ifdef BASEML
1963
2026
 
1964
 
int EigenQREVbase (FILE* fout, double kappa[], 
1965
 
                   double pi[], int *nR, double Root[], double Cijk[])
 
2027
int eigenQREVbase (FILE* fout, double kappa[], double pi[], int *nR, double Root[], double Cijk[])
1966
2028
{
1967
2029
/* pi[] is constant
1968
2030
*/
1969
 
   int i,j,k, nr=(com.ngene>1&&com.Mgene>=3?com.nrate/com.ngene:com.nrate);
1970
 
   double Q[16], U[16], V[16], mr, space_pisqrt[4];
 
2031
   int n=com.ncode, i,j,k;
 
2032
   int nr = (com.ngene>1 && com.Mgene>=3 ? com.nrate/com.ngene : com.nrate);
 
2033
   double Q[NCODE*NCODE], U[NCODE*NCODE], V[NCODE*NCODE], mr, space_pisqrt[NCODE*NCODE];
1971
2034
 
1972
2035
   NPMatUVRoot=0;
1973
 
   *nR=4;
1974
 
   zero (Q, 16);
 
2036
   *nR=n;
 
2037
   zero(Q, n*n);
1975
2038
   if(com.model==REV) {
1976
 
      for(i=0,k=0,Q[3*4+2]=Q[2*4+3]=1; i<3; i++) for (j=i+1; j<4; j++)
1977
 
         if(i*4+j!=11) Q[i*4+j]=Q[j*4+i]=kappa[k++];
 
2039
      if(n!=4) error2("ncode != 4 for REV");
 
2040
      Q[3*n+2] = Q[2*n+3] = 1;  /* r_AG = r_GA = 1. */
 
2041
      for(i=0,k=0; i<n-1; i++) for (j=i+1; j<n; j++)
 
2042
         if(i*n+j != 2*n+3)
 
2043
            Q[i*n+j] = Q[j*n+i] = kappa[k++];
1978
2044
   }
1979
2045
   else       /* (model==REVu) */
1980
 
      FOR(i,3) for(j=i+1; j<4; j++)
1981
 
         Q[i*4+j]=Q[j*4+i] = (StepMatrix[i*4+j] ? kappa[StepMatrix[i*4+j]-1] : 1);
1982
 
 
1983
 
   FOR(i,4) FOR(j,4) Q[i*4+j] *= pi[j];
1984
 
 
1985
 
   for (i=0,mr=0; i<4; i++) 
1986
 
      { Q[i*4+i]=0; Q[i*4+i]=-sum(Q+i*4, 4); mr-=pi[i]*Q[i*4+i]; }
1987
 
   abyx (1/mr, Q, 16);
 
2046
      for(i=0; i<n-1; i++) for(j=i+1; j<n; j++)
 
2047
         Q[i*n+j]=Q[j*n+i] = (StepMatrix[i*n+j] ? kappa[StepMatrix[i*n+j]-1] : 1);
 
2048
 
 
2049
   for(i=0; i<n; i++) for(j=0; j<n; j++)
 
2050
      Q[i*n+j] *= pi[j];
 
2051
 
 
2052
   for (i=0,mr=0; i<n; i++) {
 
2053
      Q[i*n+i] = 0; 
 
2054
      Q[i*n+i] = -sum(Q+i*n, n);
 
2055
      mr -= pi[i]*Q[i*n+i]; 
 
2056
   }
 
2057
   abyx(1/mr, Q, n*n);
1988
2058
 
1989
2059
   if (fout) {
1990
 
      mr=2*(pi[0]*Q[0*4+1]+pi[2]*Q[2*4+3]);
 
2060
      mr = 2*pi[0]*Q[0*n+1] + 2*pi[2]*Q[2*n+3];
1991
2061
      if(com.nhomo==0) {
1992
2062
         fprintf(fout, "\nRate parameters:  ");
1993
2063
         for(j=0; j<nr; j++) 
1994
2064
            fprintf(fout, " %8.5f", kappa[j]);
1995
2065
         fprintf(fout, "\nBase frequencies: ");
1996
 
         for(j=0; j<4; j++) 
 
2066
         for(j=0; j<n; j++) 
1997
2067
            fprintf(fout," %8.5f", pi[j]);
1998
2068
      }
1999
2069
      fprintf (fout, "\nRate matrix Q, Average Ts/Tv =%9.4f", mr/(1-mr));
2000
 
      matout (fout, Q, 4,4);
 
2070
      matout (fout, Q, n, n);
2001
2071
   }
2002
2072
   else {
2003
 
      eigenQREV(Q, pi, 4, Root, U, V, space_pisqrt);
2004
 
      FOR (i,4) FOR(j,4) FOR(k,4) Cijk[i*4*4+j*4+k] = U[i*4+k]*V[k*4+j];
 
2073
      eigenQREV(Q, pi, n, Root, U, V, space_pisqrt);
 
2074
      for(i=0; i<n; i++) for(j=0; j<n; j++) for(k=0; k<n; k++) 
 
2075
         Cijk[i*n*n+j*n+k] = U[i*n+k]*V[k*n+j];
2005
2076
   }
2006
2077
   return (0);
2007
2078
}
2012
2083
/* This constructs the rate matrix Q for the unrestricted model.
2013
2084
   pi[] is changed in the routine.
2014
2085
*/
2015
 
   int i,j,k;
2016
 
   double mr, space[20];
 
2086
   int n=com.ncode, i,j,k;
 
2087
   double mr, ts, space[20];
2017
2088
 
2018
2089
   if(com.model==UNREST) {
2019
 
      for (i=0,k=0,Q[14]=1; i<4; i++) FOR(j,4) 
2020
 
         if (i!=j && i*4+j != 14)  Q[i*4+j]=rate[k++];
 
2090
      if(n!=4) error2("ncode != 4 for UNREST");
 
2091
      for (i=0,k=0,Q[14]=1; i<n; i++) for(j=0; j<n; j++) 
 
2092
         if (i!=j && i*n+j != 14)  Q[i*n+j] = rate[k++];
2021
2093
   }
2022
2094
   else  /* (model==UNRESTu) */
2023
 
      FOR(i,4) FOR(j,4)
2024
 
         if(i!=j) 
2025
 
            Q[i*4+j] = (StepMatrix[i*4+j] ? rate[StepMatrix[i*4+j]-1] : 1);
 
2095
      for(i=0; i<n; i++) for(j=0; j<n; j++) 
 
2096
         if(i != j) 
 
2097
            Q[i*n+j] = (StepMatrix[i*n+j] ? rate[StepMatrix[i*n+j]-1] : 1);
2026
2098
 
2027
 
   FOR(i,4)  { Q[i*4+i]=0; Q[i*4+i]=-sum(Q+i*4, 4); }
 
2099
   for(i=0; i<n; i++) {
 
2100
      Q[i*n+i] = 0; 
 
2101
      Q[i*n+i] = -sum(Q+i*n, n); 
 
2102
   }
2028
2103
 
2029
2104
   /* get pi */
2030
 
 
2031
 
   QtoPi (Q, com.pi, 4, space);
2032
 
 
2033
 
   for (i=0,mr=0; i<4; i++)  mr -= pi[i]*Q[i*4+i];
2034
 
   for (i=0; i<4*4; i++)  Q[i]/=mr;
 
2105
   QtoPi(Q, com.pi, n, space);
 
2106
 
 
2107
   for (i=0,mr=0; i<n; i++)  mr -= pi[i]*Q[i*n+i];
 
2108
   for (i=0; i<n*n; i++)  Q[i] /= mr;
2035
2109
 
2036
2110
   if (fout) {
2037
 
      mr=pi[0]*Q[0*4+1]+pi[1]*Q[1*4+0]+pi[2]*Q[2*4+3]+pi[3]*Q[3*4+2];
 
2111
      ts = pi[0]*Q[0*n+1] + pi[1]*Q[1*n+0] + pi[2]*Q[2*n+3] + pi[3]*Q[3*n+2];
2038
2112
 
2039
2113
      fprintf(fout, "Rate parameters:  ");
2040
 
      FOR(j,com.nrate) fprintf(fout, " %8.5f", rate[j]);
 
2114
      for(j=0; j<com.nrate; j++)  fprintf(fout, " %8.5f", rate[j]);
2041
2115
      fprintf(fout, "\nBase frequencies: ");
2042
 
      FOR(j,4) fprintf(fout," %8.5f", pi[j]);
2043
 
      fprintf (fout,"\nrate matrix Q, Average Ts/Tv (similar to kappa/2) =%9.4f",mr/(1-mr));
2044
 
      matout (fout, Q, 4, 4);
 
2116
      for(j=0; j<n; j++) fprintf(fout," %8.5f", pi[j]);
 
2117
      if(n==4)
 
2118
         fprintf (fout,"\nrate matrix Q, Average Ts/Tv (similar to kappa/2) =%9.4f", ts/(1-ts));
 
2119
      else 
 
2120
         fprintf (fout,"\nrate matrix Q");
 
2121
      matout (fout, Q, n, n);
2045
2122
   }
2046
2123
   return (0);
2047
2124
}
2424
2501
         if(!isdigit(line[j])) { isname=1; break; }  
2425
2502
   }
2426
2503
   if(isname==0 || isname==2) {
2427
 
      sscanf(line,"%d",&k);
 
2504
      sscanf(line, "%d", &k);
2428
2505
      if(k<1||k>ns) {
2429
2506
         printf("species number %d outside range.", k);
2430
2507
         exit(-1);
2458
2535
                 node number etc.
2459
2536
 
2460
2537
   isname = 0:   species number; 1: species name; 2:both number and name
 
2538
 
 
2539
   Ziheng note (18/12/2011): I have changed the code so that sequence number is not used 
 
2540
   anymore.  isname = 1 always.
2461
2541
*/
2462
2542
   int cnode, cfather=-1;  /* current node and father */
2463
2543
   int inodeb=0;  /* node number that will have the next branch length */
2464
 
   int i,j,k, level=0, isname, ch=' ', icurspecies=0;
 
2544
   int cladeLabels=0, i,j,k, level=0, isname, ch=' ', icurspecies=0;
2465
2545
   char check[NS], delimiters[]="(),:#$=@><;", quote[]="\"\'";
2466
2546
   int lline=32000;
2467
2547
   char line[32000], *pch;
2560
2640
#endif
2561
2641
         }
2562
2642
      }
2563
 
      else if (ch=='#'||ch=='<') { *haslabel=1; fscanf(ftree,"%lf",&nodes[inodeb].label); }
2564
 
      else if (ch=='$')          { *haslabel=1; fscanf(ftree,"%d",&CladeLabel[inodeb]); }
2565
 
      else if (ch=='@'||ch=='=') { 
2566
 
         fscanf(ftree,"%lf",&nodes[inodeb].age);
 
2643
      else if (ch=='#' || ch=='<') { *haslabel=1; fscanf(ftree, "%lf", &nodes[inodeb].label); }
 
2644
      else if (ch=='$')            { *haslabel=1; fscanf(ftree, "%d",  &CladeLabel[inodeb]); }
 
2645
      else if (ch=='@' || ch=='=') { 
 
2646
         fscanf(ftree,"%lf", &nodes[inodeb].age);
2567
2647
#if (defined(BASEML) || defined(CODEML))
2568
2648
         if(com.clock) nodes[inodeb].fossil = 1;
2569
2649
#endif
2583
2663
         }
2584
2664
         for(j=i-1;j>0;j--) /* trim spaces*/
2585
2665
            if(isgraph(line[j])) break; else line[j]=0;
2586
 
         isname = IsNameNumber(line);
 
2666
 
 
2667
         if(FullSeqNames)
 
2668
            isname = 1;   /* numbers are part of names. */
 
2669
         else
 
2670
            isname = IsNameNumber(line);
2587
2671
 
2588
2672
         if (isname==2) {       /* both number and name */
2589
2673
            sscanf(line, "%d", &cnode);   cnode--;
2626
2710
         printf("\nSeq #%d occurs more than once in the tree\n",i+1); exit(-1); 
2627
2711
      }
2628
2712
      else if(check[i]<1) {
2629
 
         printf("\nSeq #%d (%s) is missing in the tree\n",i+1,com.spname[i]); exit(-1); 
 
2713
         printf("\nSeq #%d (%s) is missing in the tree\n", i+1, com.spname[i]);
 
2714
         exit(-1); 
2630
2715
      }
2631
2716
   }
2632
2717
   if(tree.nbranch>2*com.ns-2) { 
2642
2727
                         : tree.nbranch;
2643
2728
*/
2644
2729
 
2645
 
#if(defined(BASEML) || defined(CODEML))
 
2730
 
2646
2731
   /* check and convert clade labels $ */
2647
 
   if(com.clock>1 || (com.seqtype==1 && com.model>=2)) {
 
2732
#if(defined(BASEML) || defined(CODEML))
 
2733
#if(defined(BASEML))
 
2734
   if(com.seqtype==0 && com.nhomo==5) cladeLabels = 1;
 
2735
#endif
 
2736
   if(com.clock>1 || (com.seqtype==1 && com.model>=2)) cladeLabels = 1;
 
2737
   if(cladeLabels) {
2648
2738
      for(i=0,j=0; i<tree.nnode; i++) {
2649
2739
         if(CladeLabel[i] != -1) j++;
2650
2740
      }
2651
 
      if(j) {/* j is number of clade labels */
 
2741
      if(j) {  /* j is number of clade labels */
2652
2742
         DownTreeCladeLabel(tree.root, 0);
2653
2743
      }
2654
 
      else
2655
 
         for(i=0; i<tree.nnode; i++) 
2656
 
            if(i!=tree.root && nodes[i].label==-1) nodes[i].label = 0;
 
2744
 
 
2745
      /*** Somehow some labels are still -1 after this, so I changed this.  Needs checking.  ***/
 
2746
      for(i=0; i<tree.nnode; i++) 
 
2747
         if(i!=tree.root && nodes[i].label==-1) nodes[i].label = 0;
2657
2748
 
2658
2749
      /* OutTreeN(F0,1,PrBranch|PrNodeNum);  FPN(F0); */
2659
 
      /* FPN(F0);  OutTreeN(F0,1,PrLabel);  FPN(F0);  */
 
2750
      /* FPN(F0);  OutTreeN(F0,1,PrLabel);   FPN(F0); */
2660
2751
 
2661
2752
      for(i=0,com.nbtype=0; i<tree.nnode; i++) { 
2662
2753
         if(i == tree.root) continue;
2663
2754
         j = (int)nodes[i].label;
2664
 
         if(j+1 > com.nbtype)  com.nbtype=j+1;
 
2755
         if(j+1 > com.nbtype)  com.nbtype = j+1;
2665
2756
         if(j<0 || j>tree.nbranch-1)  
2666
2757
            error2("branch label in the tree (note labels start from 0 and are consecutive)");
2667
2758
      }
2679
2770
#endif
2680
2771
 
2681
2772
   }
2682
 
#endif
 
2773
#endif
2683
2774
 
2684
2775
   free(CladeLabel);
2685
2776
   return (0);
2712
2803
   if((printopt & PrLabel) && nodes[inode].label>0)
2713
2804
      fprintf(fout, labelfmt, nodes[inode].label);
2714
2805
   if((printopt & PrAge) && nodes[inode].age) 
2715
 
      fprintf(fout, " @%.3f", nodes[inode].age);
 
2806
      fprintf(fout, " @%.6f", nodes[inode].age);
2716
2807
 
2717
2808
/*  Add branch labels to be read by Rod Page's TreeView. */
2718
2809
#if (defined CODEML)
2724
2815
#endif
2725
2816
 
2726
2817
   if((printopt & PrBranch) && (inode!=tree.root || nodes[inode].branch>0))
2727
 
      fprintf(fout, ": %.6f", nodes[inode].branch);
 
2818
      fprintf(fout, ": %.5f", nodes[inode].branch);
2728
2819
   if(nsib == 0)            /* root */
2729
2820
      fputc(';', fout);
2730
2821
   else if (inode == nodes[dad].sons[nsib-1])  /* last sib */
2789
2880
*/
2790
2881
   int i,j, ison, father=nodes[inode].father, nson0=nodes[inode].nson;
2791
2882
 
2792
 
if(inode==10)
2793
 
printf("");
2794
 
 
2795
2883
   nodes[inode].label = 0;
2796
2884
   for(i=0; i<nson0; i++)
2797
2885
      PruneSubTreeN(nodes[inode].sons[i], keep);
2884
2972
   double *branch0;
2885
2973
   int debug=0;
2886
2974
 
2887
 
   if(debug) { FOR(i,com.ns) printf("%-15s %2d\n", com.spname[i], keep[i]); }
2888
 
   for(i=0,nsnew=0;i<com.ns;i++)
 
2975
   if(debug) { FOR(i,com.ns) printf("%-30s %2d\n", com.spname[i], keep[i]); }
 
2976
   for(i=0,nsnew=0; i<com.ns; i++)
2889
2977
      if(keep[i]) { nsnew++; sumnumber+=keep[i]; }
2890
2978
   if(nsnew<2)  return(-1);
2891
2979
 
2915
3003
 
2916
3004
   /* to renumber the nodes */
2917
3005
   if(sumnumber>nsnew) {
2918
 
      if(sumnumber!=nsnew*(nsnew+1)/2) error2("keep[] not right in GetSubTreeN");
 
3006
      if(sumnumber != nsnew*(nsnew+1)/2)
 
3007
         error2("keep[] not right in GetSubTreeN");
 
3008
 
2919
3009
      if((branch0=(double*)malloc(nnode0*sizeof(double)))==NULL) error2("oom#");
2920
3010
      FOR(i,nnode0) branch0[i] = nodes[i].branch;
2921
3011
      FOR(i,nnode0) newnodeNO[i] = -1;
2991
3081
   int i;
2992
3082
 
2993
3083
   for (i=com.ntime; i<np; i++) {
2994
 
      if (x[i]<xb[i][0]*1.05) x[i]=xb[i][0]*1.05;
2995
 
      if (x[i]>xb[i][1]/1.05) x[i]=xb[i][1]/1.05;
 
3084
      if (x[i]<xb[i][0]*1.005) x[i]=xb[i][0]*1.05;
 
3085
      if (x[i]>xb[i][1]/1.005) x[i]=xb[i][1]/1.05;
2996
3086
   }
2997
3087
   for (i=0; i<com.np; i++) {
2998
3088
      if (x[i]<xb[i][0]) x[i]=xb[i][0]*1.2;
4062
4152
            FOR(i,com.np) treebest.x[i]=x[i];
4063
4153
            best=itree+1;   improve=1;
4064
4154
         }
4065
 
         if (noisy) printf("%6d%2c %+8.6f", 
4066
 
                       NFunCall,(status?'?':'X'),treestar.tree.lnL-tree.lnL);
 
4155
         if (noisy) 
 
4156
            printf("%6d%2c %+8.6f", NFunCall,(status?'?':'X'),treestar.tree.lnL-tree.lnL);
4067
4157
         if (fsum) {
4068
4158
            fprintf(fsum, "%6d%2c", NFunCall, (status?'?':'X'));
4069
4159
            for (i=com.ntime; i<com.np; i++)  fprintf(fsum, "%7.3f", x[i]);
4137
4227
   if(com.ndata>1) error2("multiple data sets & multiple genes?");
4138
4228
 
4139
4229
   ngene0=com.ngene;  npatt0=com.npatt;
4140
 
   FOR (ig, ngene0)   lgene0[ig]=com.lgene[ig];
4141
 
   FOR (ig, ngene0+1) posG0[ig]=com.posG[ig];
 
4230
   for(ig=0; ig<ngene0; ig++)   lgene0[ig]=com.lgene[ig];
 
4231
   for(ig=0; ig<ngene0+1; ig++) posG0[ig]=com.posG[ig];
4142
4232
 
4143
4233
   ig=0;
4144
4234
/*
4643
4733
   jackknife if(lsb<com.ls && com.ngene==1).
4644
4734
   gmark[start+19] marks the position of the 19th site in that gene.
4645
4735
*/
4646
 
   int iboot,nboot=com.bootstrap, h,is,ig,lg[NGENE]={0},j, start;
 
4736
   int iboot,nboot=com.bootstrap, h, is, ig, lg[NGENE]={0}, j, start;
4647
4737
   int lsb=com.ls, n31=1,gap=10, gpos[NGENE];
4648
4738
   int *sites=(int*)malloc(com.ls*sizeof(int)), *gmark=NULL;
4649
 
   FILE *fseq=(FILE*)gfopen(seqf,"w");
 
4739
   FILE *fseq=(FILE*)gfopen(seqf, "w");
4650
4740
   enum {PAML=0, PAUP};
4651
 
   char *dt=(com.seqtype==AAseq?"protein":"dna");
4652
 
   char *paupstart="paupstart",*paupblock="paupblock",*paupend="paupend";
4653
 
   int format=0;  /* 0: paml-phylip; 1:paup-nexus */
 
4741
   char *datatype = (com.seqtype==AAseq?"protein":"dna");
 
4742
   char *paupstart="paupstart", *paupblock="paupblock", *paupend="paupend";
 
4743
   int  format=0;  /* 0: paml-phylip; 1:paup-nexus */
4654
4744
 
4655
4745
   if(com.readpattern) error2("work on bootstrapping pattern data.");
4656
4746
 
4658
4748
   if(format==PAUP) {
4659
4749
      printf("%s, %s, & %s will be appended if existent.\n",
4660
4750
         paupstart,paupblock,paupend);
4661
 
      appendfile(fseq,paupstart);
 
4751
      appendfile(fseq, paupstart);
4662
4752
   }
4663
4753
 
4664
4754
   if(com.seqtype==CODONseq||com.seqtype==CODON2AAseq) { n31=3; gap=1; }
4674
4764
      for(j=1; j<com.ngene; j++) com.lgene[j] += com.lgene[j-1];
4675
4765
 
4676
4766
      if(noisy && com.ngene>1) {
4677
 
         printf("Bootstrap uses stratefied sampling for %d partitions.",com.ngene);
 
4767
         printf("Bootstrap uses stratefied sampling for %d partitions.", com.ngene);
4678
4768
         printf("\nnumber of sites in each partition: ");
4679
 
         FOR(ig,com.ngene) printf(" %4d", lg[ig]);
 
4769
         for(ig=0; ig<com.ngene; ig++) printf(" %4d", lg[ig]);
4680
4770
         FPN(F0);
4681
4771
      }
4682
4772
 
4703
4793
         fprintf(fseq,"\n\n[Replicate # %d]\n", iboot+1);
4704
4794
         fprintf(fseq,"\nbegin data;\n");
4705
4795
         fprintf(fseq,"   dimensions ntax=%d nchar=%d;\n", com.ns, lsb*n31);
4706
 
         fprintf(fseq,"   format datatype=%s missing=? gap=-;\n   matrix\n",dt);
 
4796
         fprintf(fseq,"   format datatype=%s missing=? gap=-;\n   matrix\n", datatype);
4707
4797
 
4708
4798
         for(is=0;is<com.ns;is++,FPN(fseq)) {
4709
4799
            fprintf(fseq,"%-20s  ", com.spname[is]);
4719
4809
            fprintf(fseq, "\n\nbegin paup;\n");
4720
4810
            for(ig=0; ig<com.ngene; ig++)
4721
4811
               fprintf(fseq, "   charset partition%-2d = %-4d - %-4d;\n", 
4722
 
                  ig+1, (ig==0?1:com.lgene[ig-1]+1),com.lgene[ig]);
 
4812
                  ig+1, (ig==0 ? 1 : com.lgene[ig-1]+1), com.lgene[ig]);
4723
4813
            fprintf(fseq, "end;\n");
4724
4814
         }
4725
4815
         appendfile(fseq, paupblock);
4733
4823
               fprintf(fseq," %4d", lg[ig]);
4734
4824
            fprintf(fseq,"\n\n");
4735
4825
         }
4736
 
         for(is=0;is<com.ns;is++,FPN(fseq)) {
 
4826
         for(is=0; is<com.ns; is++,FPN(fseq)) {
4737
4827
            fprintf(fseq,"%-20s  ", com.spname[is]);
4738
4828
            for(h=0; h<lsb; h++) {
4739
 
               for(j=0; h<n31; h++)
 
4829
               for(j=0; j<n31; j++)
4740
4830
                  fprintf(fseq,"%c", com.z[is][sites[h]*n31+j]);
4741
4831
               if((h+1)%gap==0) fprintf(fseq," ");
4742
4832
            }
4746
4836
      if(noisy && (iboot+1)%10==0) printf("\rdid sample #%d", iboot+1);
4747
4837
   }  /* for(iboot) */
4748
4838
   free(sites);  if(com.ngene>1) free(gmark);
 
4839
   fclose(fseq);
4749
4840
   return(0);
4750
4841
}
4751
4842
 
5199
5290
   Deals with node scaling to avoid underflows.  See above 
5200
5291
   (Z. Yang, 2 Sept 2001)
5201
5292
*/
5202
 
   char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs)), *zanc;
5203
 
   char str[4]="",codon[2][4]={"   ","   "}, aa[4]="";
 
5293
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
 
5294
   char *zanc, str[4]="",codon[2][4]={"   ","   "}, aa[4]="";
5204
5295
   char *sitepatt=(com.readpattern?"pattern":"site");
5205
5296
   int n=com.ncode, inode, ic=0,b[3],i,j,k1=-1,k2=-1,c1,c2,k3, lsc=com.ls;
5206
5297
   int lst=(com.readpattern?com.npatt:com.ls);
5517
5608
   nonsynonymous changes are counted separately.
5518
5609
   Added in Nov 2000.
5519
5610
*/
5520
 
   char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs));
 
5611
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
5521
5612
   char codon[2][4]={"   ","   "};
5522
5613
   int  h,hp,inode,k1,k2,d, ls1=(com.readpattern?com.npatt:com.ls);
5523
5614
   double S,N,Sd,Nd, S1,N1,Sd1,Nd1, b,btotal=0, p,C;
5580
5671
         fprintf(frst,"%4d ",h+1);
5581
5672
         for(inode=0,d=0;inode<tree.nnode;inode++) {
5582
5673
            if(inode==tree.root) continue;
5583
 
            k1=pch[(int) zanc[(nodes[inode].father-com.ns)*com.npatt+hp] ];
 
5674
            k1 = pch[(int) zanc[(nodes[inode].father-com.ns)*com.npatt+hp] ];
5584
5675
            if(inode<com.ns)
5585
5676
               k2 = pch[com.z[inode][hp]];
5586
5677
            else  
5660
5751
   }
5661
5752
 
5662
5753
   if(inode!=tree.root) {    /* calculate log{P(t)} from father to inode */
5663
 
      t=nodes[inode].branch*_rateSite;
 
5754
      t = nodes[inode].branch*_rateSite;
5664
5755
      if(com.clock<5) {
5665
5756
         if(com.clock)  t *= GetBranchRate(igene,(int)nodes[inode].label,x,NULL);
5666
5757
         else           t *= com.rgene[igene];
5720
5811
         indexing(combScore, n*ncomb, combIndex, 0, combIndex+n*ncomb);
5721
5812
         if(debug) { printf("index "); for(i=0;i<n*ncomb; i++) printf(" %4d",combIndex[i]); FPN(F0); }
5722
5813
 
5723
 
 
5724
5814
         /* print out reconstructions at the site if inode is root. */
5725
5815
         if(inode==tree.root) {
5726
5816
            fprintf(fanc,"%4d ", h+1);
5732
5822
         psum1site=0;  /* used if inode is root */
5733
5823
 
5734
5824
         for(j=0; j<(inode!=tree.root ? NBESTANC : n*ncomb); j++) {
5735
 
            jchar=(it=combIndex[j])/ncomb; it%=ncomb;
 
5825
            jchar = (it=combIndex[j])/ncomb; it%=ncomb;
5736
5826
            if(j<NBESTANC) {
5737
 
               lnPanc[j][(inode-com.ns)*com.npatt*n+h*n+ichar]=combScore[combIndex[j]];
5738
 
               charNode[j][(inode-com.ns)*com.npatt*n+h*n+ichar]=jchar;
 
5827
               lnPanc[j][(inode-com.ns)*com.npatt*n+h*n+ichar] = combScore[combIndex[j]];
 
5828
               charNode[j][(inode-com.ns)*com.npatt*n+h*n+ichar] = jchar;
5739
5829
            }
5740
5830
            if(debug) printf("\t#%d: %6.1f %c ", j+1, combScore[combIndex[j]], BASEs[jchar]);
5741
5831
 
5776
5866
      printf("\r\tUp pass for gene %d node %d.", igene+1,inode+1);
5777
5867
}
5778
5868
 
5779
 
 
5780
5869
void DownPassPPSG2000OneSite (int h, int inode, int inodestate, int ipath)
5781
5870
{
5782
5871
/* this puts the state in ancState1site[nintern], using 
5802
5891
void PrintAncState1site (char ancState1site[], double prob)
5803
5892
{
5804
5893
   int i;
5805
 
   char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs)),codon[4]="";
 
5894
   char codon[4]="";
 
5895
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
5806
5896
 
5807
5897
   for(i=0; i<tree.nnode-com.ns; i++) {
5808
5898
      if(com.seqtype==1) {
5811
5901
#endif   
5812
5902
      }
5813
5903
      else
5814
 
         fprintf(fanc,"%c",pch[(int)ancState1site[i]]);
 
5904
         fprintf(fanc, "%c", pch[(int)ancState1site[i]]);
5815
5905
   }
5816
5906
   fprintf(fanc," (%5.3f) ", prob);
5817
5907
}
5850
5940
   This outputs results by pattern.  I tried to print results by sites, but gave up as 
5851
5941
   some variables use the same memory (e.g., combIndex) for different site patterns.
5852
5942
*/
5853
 
   char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs)),codon[4]="";
 
5943
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
 
5944
   char codon[4]="";
5854
5945
   int n=com.ncode,nintern=tree.nnode-com.ns, i,j,k,h,hp,igene;
5855
5946
   int maxnson=0, maxncomb, lst=(com.readpattern?com.npatt:com.ls);
5856
5947
   char *sitepatt=(com.readpattern?"pattern":"site");
5907
5998
   for(igene=0; igene<com.ngene; igene++) {
5908
5999
      if(com.Mgene>1) SetPGene(igene,1,1,0,x);
5909
6000
      for(i=0; i<com.ns; i++) {
5910
 
         t=nodes[i].branch*_rateSite;
 
6001
         t = nodes[i].branch*_rateSite;
5911
6002
         if(com.clock<5) {
5912
6003
            if(com.clock)  t *= GetBranchRate(igene,(int)nodes[i].label,x,NULL);
5913
6004
            else           t *= com.rgene[igene];
5975
6066
 
5976
6067
   if(com.alpha)
5977
6068
      puts("Rates are variable among sites, marginal reconstructions only.");
5978
 
   if(!com.verbose) fputs("Constant sites not listed for verbose=0\n", fout);
5979
6069
   if(!com.cleandata) fputs("Unreliable at sites with alignment gaps\n", fout);
5980
6070
 
5981
6071
   if(com.ncatG<=1 || com.method!=1)
6084
6174
 
6085
6175
   for(h=pos0; h<pos1; h++) {
6086
6176
      for(j=0,t=0;j<n;j++)
6087
 
         if(nodes[inode].conP[h*n+j]>t) t=nodes[inode].conP[h*n+j];
 
6177
         if(nodes[inode].conP[h*n+j]>t)
 
6178
            t = nodes[inode].conP[h*n+j];
6088
6179
 
6089
6180
      if(t<1e-300) {
6090
 
         for(j=0;j<n;j++)  nodes[inode].conP[h*n+j]=1;  /* both 0 and 1 fine */
6091
 
         com.nodeScaleF[k*com.npatt+h]=-800;  /* this is problematic? */
 
6181
         for(j=0;j<n;j++)
 
6182
            nodes[inode].conP[h*n+j]=1;  /* both 0 and 1 fine */
 
6183
         com.nodeScaleF[k*com.npatt+h] = -800;  /* this is problematic? */
6092
6184
      }
6093
6185
      else {  
6094
 
         for(j=0;j<n;j++)  nodes[inode].conP[h*n+j]/=t;
6095
 
         com.nodeScaleF[k*com.npatt+h]=log(t);
 
6186
         for(j=0;j<n;j++)
 
6187
            nodes[inode].conP[h*n+j]/=t;
 
6188
         com.nodeScaleF[k*com.npatt+h] = log(t);
6096
6189
      }
6097
6190
   }
6098
6191
   return(0);
6369
6462
/* P(t) for branch leading to inode, called by routines ConditionalPNode()
6370
6463
   and AncestralSeq() in baseml and codeml.  x[] is not used by baseml.
6371
6464
*/
6372
 
   double space[16] = {0};
 
6465
   int n=com.ncode, i;
 
6466
   double space[NCODE*NCODE*3] = {0};
6373
6467
 
6374
6468
   if (com.model<=K80)
6375
 
      PMatK80(Pt, t, (com.nhomo==2?nodes[inode].kappa:com.kappa));
 
6469
      PMatK80(Pt, t, (com.nhomo==2 ? *nodes[inode].pkappa : com.kappa));
6376
6470
   else {
6377
6471
      if (com.nhomo==2)
6378
 
         EigenTN93(com.model,nodes[inode].kappa, -1, com.pi,&nR,Root,Cijk);
6379
 
      else if (com.nhomo>2) /* need kappa on each node if fix_kappa ==0 */
6380
 
         EigenTN93(com.model,nodes[inode].kappa, -1, nodes[inode].pi,&nR,Root,Cijk);
 
6472
         eigenTN93(com.model, *nodes[inode].pkappa, -1, com.pi, &nR, Root, Cijk);
 
6473
      else if (com.nhomo>2 && com.model<=TN93)
 
6474
         eigenTN93(com.model, *nodes[inode].pkappa, *(nodes[inode].pkappa+1), nodes[inode].pi, &nR, Root, Cijk);
 
6475
      else if (com.nhomo>2 && com.model==REV)
 
6476
         eigenQREVbase(NULL, nodes[inode].pkappa, nodes[inode].pi, &nR, Root, Cijk);
 
6477
 
6381
6478
      if(com.model<=REV||com.model==REVu)  
6382
 
         PMatCijk(Pt,t);
 
6479
         PMatCijk(Pt, t);
6383
6480
      else {
6384
6481
         QUNREST(NULL, Pt, x+com.ntime+com.nrgene, com.pi);
6385
 
         matexp (Pt, t, 4, 10, space);
 
6482
         for(i=0; i<n*n; i++) Pt[i] *= t;
 
6483
         matexp (Pt, n, 7, 5, space);
6386
6484
      }
6387
6485
   }
6388
6486
   return(0);
6396
6494
   and AncestralSeq() in baseml and codeml.
6397
6495
 
6398
6496
   Qfactor in branch & site models (model = 2 or 3 and NSsites = 2 or 3):
6399
 
   Qfactor scaling is applied here and not inside EigenQcodon().
 
6497
   Qfactor scaling is applied here and not inside eigenQcodon().
6400
6498
*/
6401
6499
   int iUVR=0, nUVR=NBTYPE+2, ib = (int)nodes[inode].label, updateUVR=0;
6402
6500
   double *pkappa, w, mr=1, Qfactor=1;
6422
6520
         number of synonymous substitutions per codon) at their MLEs."
6423
6521
      */
6424
6522
      w = com.pomega[ib];
6425
 
      EigenQcodon(0,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, w, Pt);
 
6523
      eigenQcodon(1,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, w, Pt);
6426
6524
      Qfactor = Qfactor_NS_branch[ib];
6427
6525
   }
6428
6526
   else if (com.seqtype==CODONseq && (com.model==1 ||com.model==2) && com.nbtype<=nUVR) { 
6434
6532
      mr = 0;
6435
6533
      if(com.aaDist==AAClasses) { /* AAClass model */
6436
6534
         com.pomega = PointOmega(x+com.ntime, -1, inode, -1);
6437
 
         EigenQcodon(0,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, -1, Pt);
 
6535
         eigenQcodon(1,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, -1, Pt);
6438
6536
      }
6439
6537
      else if(com.nbtype>nUVR) {  /* branch models, with more than 8 omega */
6440
 
         EigenQcodon(0,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, nodes[inode].omega, Pt);
 
6538
         eigenQcodon(1,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, nodes[inode].omega, Pt);
6441
6539
      }
6442
6540
   }
6443
6541
 
6542
6640
            else if(com.aaDist==11) com.pomega = xcom + k + com.ncatG - 1 + 4*iclass;
6543
6641
            else if(com.aaDist==12) com.pomega = xcom + k + com.ncatG - 1 + 5*iclass;
6544
6642
         }
6545
 
         EigenQcodon(0,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, com.rK[iclass], PMat);
 
6643
         eigenQcodon(1,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, com.rK[iclass], PMat);
6546
6644
      }
6547
6645
   }
6548
6646
#endif
6671
6769
   site in the original data file or the h-th pattern.  The data are coded.
6672
6770
   naa > 1 if the codon codes for more than one amino acid.
6673
6771
*/
6674
 
   char *pch = (com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs)), compatibleAAs[20]="";
 
6772
   char *pch=(com.seqtype==0 ? BASEs : (com.seqtype==2 ? AAs: (com.seqtype==5?BASEs5:BINs)));
 
6773
   char compatibleAAs[20]="";
6675
6774
   int n=com.ncode, i, b, aa=0;
6676
6775
 
6677
6776
   for(i=0; i<com.ns; i++) {
6941
7040
   Note: At the end of the routine, nodes[].conP are not updated.
6942
7041
*/
6943
7042
   int ib,oldroot=tree.root, a,b;
6944
 
   int icycle, maxcycle=1000, icycleb, ncycleb=10, i;
 
7043
   int icycle, maxcycle=500, icycleb, ncycleb=10, i;
6945
7044
   double lnL, lnL0=0, l0,l,dl,ddl=-1, t,t0,t00, p,step=1, small=1e-20,y;
6946
7045
   double tb[2]={1e-8,50}, e=e_minbranches, *space=space_minbranches;
6947
7046
   double *xcom=x+com.ntime;  /* this is incorrect as com.ntime=0 */
7018
7117
 
7019
7118
#if (CODEML)
7020
7119
   nUVR = NBTYPE+2;
7021
 
   pkappa=(com.hkyREV||com.codonf==FMutSel ? xcom+com.nrgene : &com.kappa);
 
7120
   pkappa = (com.hkyREV||com.codonf==FMutSel ? xcom+com.nrgene : &com.kappa);
7022
7121
   if (com.seqtype==CODONseq && com.model) {
7023
7122
      if((com.model==NSbranchB || com.model==NSbranch2) && com.NSsites==0 && com.nbtype<=nUVR) {
7024
7123
         U = _UU[(int)nodes[a].label]; 
7026
7125
         Root = _Root[(int)nodes[a].label]; 
7027
7126
      }
7028
7127
      else {
7029
 
         EigenQcodon(0,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, nodes[a].omega, PMat);
 
7128
         eigenQcodon(1, -1, NULL, NULL, NULL, Root, U, V, &mr, pkappa, nodes[a].omega, PMat);
7030
7129
      }
7031
7130
   }
7032
7131
#endif
7033
7132
 
7034
7133
#if (BASEML)
7035
7134
   if (com.nhomo==2)
7036
 
      EigenTN93(com.model,nodes[a].kappa,1,com.pi,&nR,Root,Cijk);
 
7135
      eigenTN93(com.model, *nodes[a].pkappa, 1, com.pi, &nR, Root, Cijk);
7037
7136
   nroot = nR;
7038
7137
#endif
7039
7138
 
7105
7204
         Root = _Root[(int)nodes[a].label]; 
7106
7205
      }
7107
7206
      else {
7108
 
         EigenQcodon(0,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, nodes[a].omega, PMat);
 
7207
         eigenQcodon(1,-1,NULL,NULL,NULL,Root,U,V, &mr, pkappa, nodes[a].omega, PMat);
7109
7208
      }
7110
7209
   }
7111
7210
#endif
7112
7211
 
7113
7212
#if(BASEML)
7114
7213
   if (com.nhomo==2)
7115
 
      EigenTN93(com.model,nodes[a].kappa,1,com.pi,&nR,Root,Cijk);
 
7214
      eigenTN93(com.model, *nodes[a].pkappa, 1, com.pi, &nR, Root, Cijk);
7116
7215
   nroot=nR;
7117
7216
#endif
7118
7217
   *l = *dl = *ddl = 0;
7189
7288
 
7190
7289
#if (BASEML)
7191
7290
   if (com.nhomo==2)
7192
 
      EigenTN93(com.model,nodes[a].kappa,1,com.pi,&nR,Root,Cijk);
 
7291
      eigenTN93(com.model, *nodes[a].pkappa,1,com.pi,&nR,Root,Cijk);
7193
7292
   nroot=nR;
7194
7293
#endif
7195
7294
 
7312
7411
 
7313
7412
#if (BASEML)
7314
7413
   if (com.nhomo==2)
7315
 
      EigenTN93(com.model, nodes[a].kappa, 1, com.pi, &nR, Root, Cijk);
 
7414
      eigenTN93(com.model, *nodes[a].pkappa, 1, com.pi, &nR, Root, Cijk);
7316
7415
   nroot=nR;
7317
7416
#endif
7318
7417
   if(com.NnodeScale==0)
7549
7648
   (ii) for ReadTreeN to be able to read other node labels such as #, $, either with
7550
7649
   or without ' '.
7551
7650
*/
7552
 
   int i,j,k, nfossiltype=8;
 
7651
   int i,j,k, nfossiltype=7;
7553
7652
   char *pch;
7554
7653
   double tailL=0.025, tailR=0.025, p_LOWERBOUND=0.1, c_LOWERBOUND=1.0;
7555
7654
 
7796
7895
int GenerateGtree (int locus)
7797
7896
{
7798
7897
/* construct the gene tree at locus by pruning tips in the master species 
7799
 
   tree.  com.spname[] have names of species at the current locus and 
7800
 
   the routine use them to compare with sptree.nodes[].name to decide which 
7801
 
   species to keep for the locus.  See GetSubTreeN() for more info.
 
7898
   tree.  com.spname[] have names of species at the current locus (probably read 
 
7899
   from the sequence alignment at the locus).  They are used by the routine to compare 
 
7900
   with sptree.nodes[].name to decide which species to keep for the locus.  
 
7901
   See GetSubTreeN() for more details.
7802
7902
*/
7803
7903
   int ns=data.ns[locus], i,j, ipop[NS], keep[NS], newnodeNO[2*NS-1];
7804
7904
 
7805
 
   for(j=0;j<sptree.nspecies;j++) keep[j]=0;
 
7905
   for(j=0; j<sptree.nspecies; j++) keep[j]=0;
7806
7906
   for(i=0;i<ns;i++) {
7807
7907
      for(j=0;j<sptree.nspecies;j++)
7808
7908
         if(!strcmp(com.spname[i], sptree.nodes[j].name)) break;
7810
7910
         printf("species %s not found in master tree\n", com.spname[i]);
7811
7911
         exit(-1);
7812
7912
      }
7813
 
      keep[j]=i+1; ipop[i]=j;
 
7913
      if(keep[j]) {
 
7914
         printf("\nspecies %s occurs twice in locus %d", com.spname[i], locus+1);
 
7915
         error2("\ngiving up...");
 
7916
      }
 
7917
      keep[j] = i+1;  ipop[i] = j;  /* seq j in alignment is species i in master tree. */
7814
7918
      free(com.spname[i]);
7815
7919
   }
 
7920
 
7816
7921
   /* copy master species tree and then prune it. */
7817
7922
   copySptree();
7818
7923
   GetSubTreeN(keep, newnodeNO);
7912
8017
int ResetFinetuneSteps(FILE *fout, double Pjump[], double finetune[], int nsteps)
7913
8018
{
7914
8019
   int j, verybadstep=0;
7915
 
   double PjumpOpt = 0.30; /* this is the optimum for the 2NormalSym move. */
7916
 
   double P0Change[2]={0.2, 0.5}; /* no change if in interval? */
7917
 
 
7918
 
   printf("\n\nCurrent Pjump:    ");
7919
 
   for(j=0; j<nsteps; j++)
7920
 
      printf(" %8.5f", Pjump[j]);
7921
 
   printf("\nCurrent finetune: ");
7922
 
   for(j=0; j<nsteps; j++)
7923
 
      printf(" %8.5f", finetune[j]);
7924
 
 
 
8020
   double PjumpOpt = 0.30; /* this is the optimum for the Bactrian move. */
 
8021
 
 
8022
   if(noisy>=3) {
 
8023
      printf("\n\nCurrent Pjump:    ");
 
8024
      for(j=0; j<nsteps; j++)
 
8025
         printf(" %8.5f", Pjump[j]);
 
8026
      printf("\nCurrent finetune: ");
 
8027
      for(j=0; j<nsteps; j++)
 
8028
         printf(" %8.5f", finetune[j]);
 
8029
   }
7925
8030
   if(fout) {
7926
8031
      fprintf(fout, "\nCurrent Pjump:    ");
7927
8032
      for(j=0; j<nsteps; j++)
7941
8046
         verybadstep = 1;
7942
8047
      }
7943
8048
      else
7944
 
         finetune[j] *= tan(Pi/2*Pjump[j])/tan(Pi/2*PjumpOpt);
7945
 
      /*
7946
 
      else if(Pjump[j] < P0Change[0] || Pjump[j] > P0Change[1])
7947
 
         finetune[j] *= tan(Pi/2*Pjump[j])/tan(Pi/2*PjumpOpt);
7948
 
      */
7949
 
   }
7950
 
 
7951
 
   printf("\nNew     finetune: ");
7952
 
   for(j=0; j<nsteps; j++)
7953
 
      printf(" %8.5f", finetune[j]);
7954
 
   printf("\n\n");
7955
 
 
 
8049
         finetune[j] *= tan(Pi/2*Pjump[j]) / tan(Pi/2*PjumpOpt);
 
8050
   }
 
8051
 
 
8052
   if(noisy>=3) {
 
8053
      printf("\nNew     finetune: ");
 
8054
      for(j=0; j<nsteps; j++)
 
8055
         printf(" %8.5f", finetune[j]);
 
8056
      printf("\n\n");
 
8057
   }
7956
8058
   if(fout) {
7957
8059
      fprintf(fout, "\nNew     finetune: ");
7958
8060
      for(j=0; j<nsteps; j++)
8011
8113
         GetDaa(NULL, com.daa);
8012
8114
         if(com.model==Empirical_F) 
8013
8115
            xtoy(data.pi[locus], com.pi, nc);
8014
 
         EigenQaa(NULL, _Root[locus], _UU[locus], _VV[locus], NULL);
8015
 
 
8016
 
printf("Protein # %2d uses %-20s\n", locus+1,data.daafile[locus]);
8017
 
matout(F0, com.pi, 1, nc);
8018
 
matout(F0, _Root[locus], 1, nc);
 
8116
         eigenQaa(NULL, _Root[locus], _UU[locus], _VV[locus], NULL);
8019
8117
      }
8020
8118
   }
8021
8119
   else if(com.seqtype==1 && com.fix_kappa & com.fix_omega) {
8027
8125
         com.kappa=data.kappa[locus];
8028
8126
         com.omega=data.omega[locus];
8029
8127
         xtoy(data.pi[locus], com.pi, com.ncode);
8030
 
         EigenQcodon(0,-1,NULL,NULL,NULL, _Root[locus], _UU[locus], _VV[locus], &mr,
 
8128
         eigenQcodon(1,-1,NULL,NULL,NULL, _Root[locus], _UU[locus], _VV[locus], &mr,
8031
8129
            &com.kappa, com.omega, PMat);
8032
8130
      }
8033
8131
   }
8076
8174
      if(com.seqtype==0 && com.model!=0 && com.model!=1)
8077
8175
         xtoy(data.pi[locus], com.pi, com.ncode);
8078
8176
      if(com.model<=TN93)
8079
 
         EigenTN93(com.model, com.kappa, com.kappa, com.pi, &nR, Root, Cijk);
 
8177
         eigenTN93(com.model, com.kappa, com.kappa, com.pi, &nR, Root, Cijk);
8080
8178
      else if (com.model==REV)
8081
 
         EigenQREVbase (NULL, &com.kappa, com.pi, &nR, Root, Cijk);
 
8179
         eigenQREVbase (NULL, &com.kappa, com.pi, &nR, Root, Cijk);
8082
8180
#else
8083
8181
      if((com.seqtype==1 && com.codonf) || (com.seqtype==2 && com.model==3))
8084
8182
         xtoy(data.pi[locus], com.pi, com.ncode);
8088
8186
         Root=_Root[locus]; U=_UU[locus];  V=_VV[locus];
8089
8187
      }
8090
8188
      else {
8091
 
         EigenQcodon(0,-1,NULL,NULL,NULL,Root,U,V, &mr, &com.kappa, com.omega,PMat);
 
8189
         eigenQcodon(1,-1,NULL,NULL,NULL,Root,U,V, &mr, &com.kappa, com.omega,PMat);
8092
8190
      }
8093
8191
 
8094
8192
#endif