debian/changelog

+paml (4.9i+dfsg-1) UNRELEASED; urgency=medium
+
+  * New upstream version
+  * debhelper-compat 12
+  * Standards-Version: 4.4.0
+
+ -- Andreas Tille <tille@debian.org>  Thu, 01 Aug 2019 19:19:13 +0200
+
 paml (4.9h+dfsg-1) unstable; urgency=medium
 
   * New upstream version

debian/compat

-11

debian/control

@@ -5,8 +5,8 @@ Uploaders: Pjotr Prins <pjotr.debian@thebird.nl>,
            Andreas Tille <tille@debian.org>
 Section: science
 Priority: optional
-Build-Depends: debhelper (>= 11~)
-Standards-Version: 4.1.4
+Build-Depends: debhelper-compat (= 12)
+Standards-Version: 4.4.0
 Vcs-Browser: https://salsa.debian.org/med-team/paml
 Vcs-Git: https://salsa.debian.org/med-team/paml.git
 Homepage: http://abacus.gene.ucl.ac.uk/software/paml.html

doc/pamlHistory.txt

@@ -8,6 +8,41 @@ discussion site
 https://groups.google.com/forum/#!forum/pamlsoftware.
 
 
+Version 4.9i, March 2019
+
+(*) mcmctree: I have added an option (duplication = 1) for dating a
+tree with both speciations and gene duplications, so that some nodes
+on the tree share divergence times.  Nodes sharing ages are identified
+using labels in the tree file: #1, #2, ....  I have yet to update the
+document about specification of the model.
+
+(*) mcmctree: The TipDate option was written for one locus or
+partition and never worked for more than two loci/partitions.  I have
+edited the code so that it works for multiple partitions, some of
+which may be molecular and the others morphological.
+
+(*) codeml: The option estFreq = 0 when codonFreq = 6 (FMutSel0) and 7
+(FMutSel) is not working in versions 4.9g and 4.9h.  This is fixed
+now.  This option uses the observed codon or amino acid frequencies
+for the mutation-selection models of codon usage.  Instead the program
+estimates the frequencies using maximum likeihood, which is what the
+option estFreq = 1 does.  Look at the README file in the
+examples/mtCDNAape/ folder.
+
+(*) codeml clade model D: The bounds for the w (dN/dS) ratios in the
+first site classes are set tp (0.0001, 0.5) for w0 and (0.5, 1.5) for
+w1, in versions 4.9b,c,d,e,f,g, since I added the BEB calculation for
+clade model D in 4.9b.  The motivation for the bounds is that site
+class 0 represents strong purifying selection with a small w0, while
+site class 1 should include sites under weak purifying selection with
+a larger w1.  However the bounds are arbitrary.  In some datasets, the
+MLEs are found to be at the bounds, making the interpretation awkward.
+I have changed the bounds to the following: 
+w0b[]={0.0001, 1.0}, w1b[]={0.01, 1.5}.
+This means that the user should swap the estimates of w0 and w2 if w0 > w1.
+
+
+
 Version 4.9h, March 2018
 
 (*) mcmctree: gamma-Dirichlet versus conditional i.i.d. priors for

examples/DatingSoftBound/mtCDNApri.trees

  7  1
 
-
-((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla),
-(orangutan, sumatran)) '>.12<.16', gibbon);
-
+((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);
 
 //end of file
 
 
-((((human, (chimpanzee, bonobo)) 'B(.06, .08)', gorilla),
-(orangutan, sumatran)) 'B(.12, .16)', gibbon);
+((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);
+
+((((human, (chimpanzee, bonobo)) 'B(.06, .08)', gorilla), (orangutan, sumatran)) 'B(.12, .16)', gibbon);

examples/TipDate.FluH1/baseml.ctl

@@ -18,7 +18,7 @@
     fix_kappa = 0   * 0: estimate kappa; 1: fix kappa at value below
         kappa = 2   * initial or fixed kappas
  
-    fix_alpha = 0   * 0: estimate alpha; 1: fix alpha at value below
+    fix_alpha = 1   * 0: estimate alpha; 1: fix alpha at value below
         alpha = 0.25   * initial or fixed alpha, 0:infinity (constant rate)
         ncatG = 5   * # of categories in the dG, AdG, or nparK models of rates
       fix_rho = 1   * 0: estimate rho; 1: fix rho at value below

examples/TipDate.FluH1/mcmctreeClock1.ctl

-          seed = -1357
+          seed = 1357
       seqfile = H1.txt
 *      seqfile = ../H1.CP123.txt
      treefile = H1.tre
@@ -7,7 +7,7 @@
        seqtype = 0  * 0: nucleotides; 1:codons; 2:AAs
          noisy = 3
 
-       usedata = 0 in.BV.HKYG5  * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
+       usedata = 2 in.BV.HKYG5  * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
 
          ndata = 1   * 
          clock = 1   * 1: global clock; 2: independent rates; 3: correlated rates

examples/TipDate.HIV2/FigTree.tre

+#NEXUS
+BEGIN TREES;
+
+	UTREE 1 = (((SL92eSM1992: 0.381307, (SL92fSM1992: 0.154397, SL92aSM1992: 0.154397) [&95%={0.0347455, 0.522379}]: 0.226910) [&95%={0.0904505, 0.907043}]: 0.504653, (((Lib1sm1988: 0.142534, FO784h1989: 0.152534) [&95%={0.0748715, 0.520538}]: 0.202566, SIVMNE1982: 0.285100) [&95%={0.152587, 0.820146}]: 0.357506, (SIVSTM1989: 0.615921, ((((((JA1990: 0.060047, ON1990: 0.060047) [&95%={0.0520286, 0.251389}]: 0.066570, EHO1988: 0.106617) [&95%={0.0769271, 0.358953}]: 0.073194, UC1h1988: 0.179811) [&95%={0.105804, 0.465565}]: 0.087350, D205h1986: 0.247161) [&95%={0.152152, 0.584476}]: 0.111152, 2238h1989: 0.388313) [&95%={0.223053, 0.723296}]: 0.156536, (((D194h1987: 0.111214, UC2h1988: 0.121214) [&95%={0.0835921, 0.423991}]: 0.141317, BEN1987: 0.252531) [&95%={0.124754, 0.605509}]: 0.202257, (CAM2CG1987: 0.390471, (ROD1985: 0.316609, ((P07h1995: 0.259167, (P09h1995: 0.168184, NIHZ1986: 0.078184) [&95%={0.0925188, 0.329415}]: 0.090983) [&95%={0.120916, 0.445566}]: 0.110770, (P08h1995: 0.318209, (P010h1995: 0.272484, ((P01h1995: 0.188251, ((P02h1995: 0.071261, P06h1995: 0.071261) [&95%={0.0031347, 0.190992}]: 0.077251, (P05h1995: 0.102844, ALI1989: 0.042844) [&95%={0.0612933, 0.200128}]: 0.045668) [&95%={0.0744984, 0.27015}]: 0.039739) [&95%={0.0946369, 0.324612}]: 0.043111, (P04h1995: 0.107655, P03h1995: 0.107655) [&95%={0.0046477, 0.272953}]: 0.123707) [&95%={0.121765, 0.381337}]: 0.041122) [&95%={0.149241, 0.434521}]: 0.045725) [&95%={0.182706, 0.493541}]: 0.051729) [&95%={0.223594, 0.558168}]: 0.046672) [&95%={0.255413, 0.624248}]: 0.053862) [&95%={0.292723, 0.698665}]: 0.064317) [&95%={0.337601, 0.788576}]: 0.070062) [&95%={0.379475, 0.895678}]: 0.071072) [&95%={0.41195, 1.01957}]: 0.096685) [&95%={0.44965, 1.207}]: 0.143354) [&95%={0.48519, 1.51936}]: 0.312910, (SL2h1993: 0.423337, (SL92bSM1992: 0.160871, SL92cSM1992: 0.160871) [&95%={0.0348426, 0.554373}]: 0.252466) [&95%={0.0928551, 1.06356}]: 0.785533) [&95%={0.51691, 1.97612}];
+
+END;
+
+[Note for FigTree: Under Time Scale, set Offset = 1995.0, Scale factor = -100.0
+Untick Scale Bar, & tick Tip Labels, Node Bars, Scale Axis, Reverse Axis, Show Grid.]

src/Makefile.MSVC → src/Makefile.MSC

 PRGS =  baseml codeml basemlg mcmctree pamp evolver yn00 chi2 
 CC = cl # cc, gcc, cl
-CFLAGS =  -O2 -Ot 
+CFLAGS =  -Ox 
 
 DEFINE = #
 LIBS = #-lm -lM

src/baseml.c

@@ -30,7 +30,7 @@ int Forestry(FILE *fout);
 void DetailOutput(FILE *fout, double x[], double var[]);
 int GetOptions(char *ctlf);
 int GetInitials(double x[], int *fromfile);
-int GetInitialsTime(double x[]);
+int GetInitialsTimes(double x[]);
 int SetxInitials(int np, double x[], double xb[][2]);
 int SetParameters(double x[]);
 int SetPGene(int igene, int _pi, int _UVRoot, int _alpha, double x[]);
@@ -52,7 +52,7 @@ int GroupDistances();
 
 struct CommonInfo {
    unsigned char *z[NS];
-   char *spname[NS], seqf[512], outf[512], treef[512], cleandata;
+   char *spname[NS], seqf[2048], outf[2048], treef[2048], cleandata;
    char oldconP[NNODE];  /* update conP for nodes to save computation (0 yes; 1 no) */
    int seqtype, ns, ls, ngene, posG[NGENE + 1], lgene[NGENE], *pose, npatt, readpattern;
    int np, ntime, nrgene, nrate, nalpha, npi, nhomo, ncatG, ncode, Mgene;
@@ -72,14 +72,16 @@ struct CommonInfo {
    int fix_omega;
    double omega;
 }  com;
+
 struct TREEB {
    int  nbranch, nnode, root, branches[NBRANCH][2];
    double lnL;
 }  tree;
+
 struct TREEN {
    int father, nson, sons[MAXNSONS], ibranch, ipop;
-   double branch, age, *pkappa, pi[4], *conP, label;
-   char *nodeStr, fossil;
+   double branch, age, *pkappa, pi[4], *conP, label, label2;
+   char fossil, *name, *annotation;
 }  *nodes, **gnodes, nodes_t[2 * NS - 1];
 
 
@@ -92,7 +94,7 @@ struct SPECIESTREE {
    struct TREESPN {
       char name[LSPNAME + 1], fossil, usefossil;  /* fossil: 0, 1, 2, 3 */
       int father, nson, sons[2];
-      double age, pfossil[7];   /* lower and upper bounds or alpha & beta */
+      double label, age, pfossil[7];   /* lower and upper bounds or alpha & beta */
       double *lnrates;          /* log rates for loci */
    } nodes[2 * NS - 1];
 }  stree;
@@ -133,7 +135,7 @@ int main(int argc, char *argv[])
 {
    FILE *fout, *fseq = NULL, *fpair[6];
    char pairfs[1][32] = { "2base.t" };
-   char rstf[512] = "rst", ctlf[512] = "baseml.ctl", timestr[64];
+   char rstf[2048] = "rst", ctlf[2048] = "baseml.ctl", timestr[64];
    char *Mgenestr[] = { "diff. rate", "separate data", "diff. rate & pi",
                      "diff. rate & kappa", "diff. rate & pi & kappa" };
    int getdistance = 1, i, idata;
@@ -266,7 +268,13 @@ int main(int argc, char *argv[])
             fprintf(fout, "\n\nPrinting out site pattern counts\n");
             printPatterns(fout);
          }
-         if (com.verbose == 2) {
+         if (com.verbose >= 2) {
+            fprintf(fout, "\nSite-to-pattern map: ");
+            for (i = 0; i < com.ls; i++)
+               fprintf(fout, " %2d", com.pose[i] + 1);
+            fprintf(fout, "\n");
+
+            fprintf(fout, "\n**** Alignment mutated for JC69-like models ****\n");
             fprintf(fout, "\n%6d %6d\n", com.ns, com.ls);
             for (i = 0; i < com.ns; i++) {
                fprintf(fout, "\n%-30s  ", com.spname[i]);
@@ -410,9 +418,11 @@ int Forestry(FILE *fout)
    fprintf(flnf, "%6d %6d %6d\n", ntree, com.ls, com.npatt);
 
    for (itree = 0; ntree == -1 || itree < ntree; itree++, iteration = 1) {
-      if (ReadTreeN(ftree, &haslength, &i, 0, 1)) {
+      if (ReadTreeN(ftree, &haslength, 0, 1)) {
          puts("\nend of tree file.");   break;
       }
+      ProcessNodeAnnotation(&i);
+
       if (noisy) printf("\nTREE # %2d: ", itree + 1);
       fprintf(fout, "\nTREE # %2d:  ", itree + 1);
       fprintf(frub, "\n\nTREE # %2d\n", itree + 1);
@@ -490,7 +500,9 @@ int Forestry(FILE *fout)
          else                                  status = 0;
       }
 
-      if (itree == 0) { lnL0 = lnLbest = lnL; btree = 0; }
+      if (itree == 0) {
+         lnL0 = lnLbest = lnL; btree = 0;
+      }
       else if (lnL < lnLbest) {
          lnLbest = lnL;  btree = itree;
       }
@@ -623,12 +635,12 @@ int Forestry(FILE *fout)
 
 void PrintBestTree(FILE *fout, FILE *ftree, int btree)
 {
-   int itree, ntree, i, j;
+   int itree, ntree, i;
 
    rewind(ftree);
    GetTreeFileType(ftree, &ntree, &i, 0);
    for (itree = 0; ntree == -1 || itree < ntree; itree++) {
-      if (ReadTreeN(ftree, &i, &j, 0, 1)) {
+      if (ReadTreeN(ftree, &i, 0, 1)) {
          puts("\nend of tree file."); break;
       }
       if (itree == btree)
@@ -1473,28 +1485,27 @@ int ConditionalPNode(int inode, int igene, double x[])
    return (0);
 }
 
-
 int PMatCijk(double P[], double t)
 {
 /* P(t)ij = SUM Cijk * exp{Root*t}
 */
    int i, j, k, n = com.ncode, nr = nR;
-   double expt[5], pij;
-
-   if (t < -.001 && noisy>3)
-      printf("\nt = %.5f in PMatCijk", t);
-   if (t < 1e-200) { identity(P, n); return(0); }
+   double exptm1[5] = {0};
 
-   for (k = 1; k < nr; k++) expt[k] = exp(t*Root[k]);
-   for (i = 0; i < n; i++) for (j = 0; j < n; j++) {
-      for (k = 1, pij = Cijk[i*n*nr + j*nr + 0]; k < nr; k++)
-         pij += Cijk[i*n*nr + j*nr + k] * expt[k];
-      P[i*n + j] = (pij > 0 ? pij : 0);
+   memset(P, 0, n*n * sizeof(double));
+   for (k = 1; k < nr; k++) exptm1[k] = expm1(t*Root[k]);
+   for (i = 0; i < n; i++) {
+      for (j = 0; j < n; j++) {
+         for (k = 0; k < nr; k++)
+            P[i*n + j] += Cijk[i*n*nr + j*nr + k] * exptm1[k];
+      }
+      P[i*n + i] ++;
    }
    return (0);
 }
 
 
+
 #ifdef UNDEFINED
 
 int CollapsSite(FILE *fout, int nsep, int ns, int *ncat, int SiteCat[])

src/basemlg.c

@@ -40,7 +40,7 @@ int lfunG_d (double x[], double *lnL, double dl[], int np);
 int lfunG_dd (double x[], double *lnL, double dl[], double ddl[], int np);
 
 struct CommonInfo {
-   unsigned char *z[NS], *spname[NS], seqf[96],outf[96],treef[96];
+   unsigned char *z[NS], *spname[NS], seqf[2048],outf[2048],treef[2048];
    int  seqtype, ns, ls, ngene, posG[NGENE+1],lgene[NGENE],*pose,npatt, readpattern;
    int  clock,fix_alpha,fix_kappa,fix_rgene,Malpha,print,verbose;
    int  model, runmode, cleandata, ndata;
@@ -58,8 +58,8 @@ struct TREEB {
 
 struct TREEN {
    int father, nson, sons[MAXNSONS], ibranch;
-   double branch, age, label, *conP;
-   char *nodeStr, fossil;
+   double branch, age, label, label2, *conP;
+   char *annotation, fossil;
 }  nodes[2*NS-1];
 
 static int nR=4, CijkIs0[64];
@@ -83,7 +83,7 @@ int main (int argc, char *argv[])
 {
    FILE *fout, *fseq=NULL, *fpair[6];
    char pairfs[1][32]={"2base.t"};
-   char ctlf[96]="baseml.ctl";
+   char ctlf[2048]="baseml.ctl";
    double  space[50000];
 
    noisy=2;  com.cleandata=1;  /* works with clean data only */
@@ -158,7 +158,7 @@ int Forestry (FILE *fout, double space[])
       fprintf(frub,"\n\nTREE # %2d", itree+1);
       fprintf(frate,"\nTREE # %2d\n", itree+1);
 
-      if(ReadTreeN(ftree,&haslength,&i,0,1))
+      if (ReadTreeN(ftree, &haslength, 0, 1))
            { puts("err or end of tree file."); break; }
 
       com.ntime = com.clock ? tree.nnode-com.ns: tree.nbranch;
@@ -248,7 +248,7 @@ int SetBranch (double x[])
 int testx (double x[], int np)
 {
    int i,k;
-   double tb[]={1e-5,4}, rgeneb[]={.01,20}, kappab[]={0,80}, alphab[]={.01,99};
+   double tb[]={1e-5,99}, rgeneb[]={.01,20}, kappab[]={0,80}, alphab[]={.01,999};
 
    if (SetBranch(x)) return (-1);
    FOR (i,com.ntime)   if (x[i]<tb[0] || x[i]>tb[1])   return (-1);
@@ -266,8 +266,8 @@ int SetxBound (int np, double xb[][2])
 /* sets lower and upper bounds for variables during iteration
 */
    int i,j, k=com.ntime+com.nrgene+com.nrate;
-   double tb[]={.4e-5, 20}, rgeneb[]={1e-4,99}, rateb[]={1e-5,999};
-   double alphab[]={0.005, 99}, pb[2]={1e-5,1-1e-5};
+   double tb[]={.4e-5, 99}, rgeneb[]={1e-4,99}, rateb[]={1e-5,999};
+   double alphab[]={0.005, 999}, pb[2]={1e-5,1-1e-5};
 
    if(com.clock) {
       xb[0][0]=tb[0];  xb[0][1]=tb[1];

src/codeml.c

@@ -106,7 +106,7 @@ double CDFLogis(double x, double m, double s);
 
 struct common_info {
    unsigned char *z[NS];
-   char *spname[NS], seqf[512], outf[512], treef[512], daafile[512], cleandata;
+   char *spname[NS], seqf[2048], outf[2048], treef[2048], daafile[2048], cleandata;
    char oldconP[NNODE];       /* update conP for nodes? to save computation */
    int seqtype, ns, ls, ngene, posG[NGENE + 1], lgene[NGENE], npatt, *pose, readpattern;
    int runmode, clock, verbose, print, codonf, aaDist, model, NSsites;
@@ -132,14 +132,16 @@ struct common_info {
      & eigenQcodon.  Try to remove them? */
    double *pomega, pkappa[5], *ppi;
 }  com;
+
 struct TREEB {
    int  nbranch, nnode, root, branches[NBRANCH][2];
    double lnL;
 }  tree;
+
 struct TREEN {
    int father, nson, sons[MAXNSONS], ibranch, ipop;
-   double branch, age, omega, *conP, label;
-   char *nodeStr, fossil, usefossil;
+   double branch, age, omega, *conP, label, label2;
+   char fossil, usefossil, *name, *annotation;
 }  *nodes, **gnodes, nodes_t[2 * NS - 1];
 
 
@@ -152,7 +154,7 @@ struct SPECIESTREE {
    struct TREESPN {
       char name[LSPNAME + 1], fossil, usefossil;  /* fossil: 0, 1, 2, 3 */
       int father, nson, sons[2];
-      double age, pfossil[7];   /* lower and upper bounds or alpha & beta */
+      double label, age, pfossil[7];   /* lower and upper bounds or alpha & beta */
       double *lnrates;          /* log rates for loci */
    } nodes[2 * NS - 1];
 }  stree;
@@ -643,10 +645,11 @@ int Forestry(FILE *fout)
    }
 
    for (itree = 0; ntree == -1 || itree < ntree; itree++, iteration = 1) {
-      if (ReadTreeN(ftree, &haslength, &i, 0, 1)) {
+      if (ReadTreeN(ftree, &haslength, 0, 1)) {
          puts("end of tree file."); break;
       }
-
+      if (com.seqtype == 1 && com.model)
+         ProcessNodeAnnotation(&k);
       printf("\nTREE # %2d\n", itree + 1);
       fprintf(fout, "\n\nTREE # %2d:  ", itree + 1);
       fprintf(flnf, "\n\n%2d\n", itree + 1);
@@ -903,15 +906,13 @@ double *PointOmega(double xcom[], int igene, int inode, int isiteclass)
       com.NSsites and igene.  This sometimes points to com.omega or com.rK[].
       This is called by SetParameters(), DetailOutput(), etc.
 
+      Ziheng Notes: 8 August 2003.
       Difficulties in using this with lfunt() etc.
-
       Trying to remove global variables com.pomega and com.pkappa through
       PointOmega and PointKappa, but was unsuccessful when too many changes were
       made at the same time.  Perhaps look at this later.  Note that some
       variables are passed over the head from lfunt() etc. to eigenQcodon().
 
-      Ziheng Notes: 8 August 2003.
-
    */
    int k = com.nrgene + com.nkappa, backfore;
    int nka = (com.hkyREV ? 5 : 1), nw = (com.aaDist == AAClasses ? com.nOmegaType : 1);
@@ -1472,13 +1473,13 @@ int getnrate(int firsttime)
          if (com.aaDist) com.nrate++;
          if (!com.fix_omega)  com.nrate++;
          if (com.codonf) {
-            com.codonf = 0;  puts("CodonFreq=0 reset for model=6.");
+            com.codonf = 0;  puts("CodonFreq=0 reset for model=5.");
          }
          break;
       case (FromCodon):
          com.nrate = com.nkappa = (com.hkyREV ? 5 : !com.fix_kappa);
          if (com.aaDist) com.nrate++;
-         if(com.fix_omega) error2("fix_omega = 1?  omega is inestimable!");
+         if (com.fix_omega) error2("fix_omega = 1?  omega is not estimable!");
          com.omega = -1;
          if (com.codonf) {
             com.codonf = 0;  puts("CodonFreq=0 reset for model=6.");
@@ -1531,6 +1532,7 @@ int getnrate(int firsttime)
       /* pi_T, pi_C, pi_A are counted as frequency parameters pi. */
       if (com.codonf == 0)
          com.npi = 0;
+      if (firsttime) {
          if (com.codonf == FMutSel0)     /* FMutSel0: pi_TCA plus 20 AA freqs. */
             com.npi = 3 + (com.npi ? 20 - 1 : 0);
          else if (com.codonf == FMutSel) /* FMutSel:  pi_TCA plus 60 codon freqs. */
@@ -1540,6 +1542,7 @@ int getnrate(int firsttime)
             else if (com.codonf == F3x4 || com.codonf == F3x4MG)  com.npi = 9;
             else if (com.codonf == Fcodon)                        com.npi = com.ncode - 1;
          }
+      }
       com.nrate += com.npi;
 
       if (com.aaDist != AAClasses) {
@@ -1846,7 +1849,7 @@ int SetxBound(int np, double xb[][2])
    double alphab[] = { 0.02,49 }, betab[] = { 0.005,99 };
    double rhob[] = { 0.01,0.99 }, pb[] = { .00001,.99999 };
    double wb[] = { 0.0001,999 };
-   double w0b[] = { 0.0001, 0.5 }, w1b[] = { 0.5, 1.5 };  /* clade model D */
+   double w0b[] = { 0.0001, 1.0 }, w1b[] = { 0.01, 1.5 };  /* clade model D */
 
    SetxBoundTimes(xb);
    for (i = com.ntime; i < np; i++) for (j = 0; j < 2; j++) xb[i][j] = rateb[j];
@@ -5157,7 +5160,6 @@ int lfunNSsites_rate(FILE* frst, double x[], int np)
    char sitelabel[96], *colors[5] = { "darkblue", "lightblue", "purple", "pinkred", "red" };
    char *colorvalues[5] = { "[2,2,120]", "[133,57,240]", "[186,60,200]", "[200,60,160]", "[250,5,5]" };
 
-
    if (com.nparK) error2("lfunNSsites_rate to be done for HMM.");
 
    if ((meanw = (double*)malloc(com.npatt * sizeof(double))) == NULL)
@@ -6722,8 +6724,8 @@ int lfunNSsites_ACD(FILE* frst, double x[], int np)
       8 March 2016, added BEB for clade model D.  deleted ncatG=2 for D so that ncatG=3 for A C D.
    */
    int n1d = 10, debug = 0, site = 0;
-   double w0b[] = { 0.0, 0.5 };  /* bounds for w0 */
-   double w1b[] = { 0.5, 1.5 };  /* bounds for w1 (for model D only) */
+   double w0b[] = { 0.0, 1.0 };  /* bounds for w0 */
+   double w1b[] = { 0.01, 1.5 }; /* bounds for w1 (for model D only) */
    double wsb[] = { 1, 11 };     /* bounds for w2 in model A */
    double w2b[] = { 0, 3 };      /* bounds for w2 w3 w4... in clade models C & D */
    enum ModelsACD { mA, mC, mD };

src/evolver.c

 /* evolver.c
    Copyright, Ziheng Yang, April 1995.
 
-     cl -Ot -O2 evolver.c tools.c
-     cl -Ot -O2 -DCodonNSbranches    -FeevolverNSbranches.exe    evolver.c tools.c
-     cl -Ot -O2 -DCodonNSsites       -FeevolverNSsites.exe       evolver.c tools.c
-     cl -Ot -O2 -DCodonNSbranchsites -FeevolverNSbranchsites.exe evolver.c tools.c
+     cl -Ox evolver.c tools.c
+     cl -Ox -DCodonNSbranches    -FeevolverNSbranches.exe    evolver.c tools.c
+     cl -Ox -DCodonNSsites       -FeevolverNSsites.exe       evolver.c tools.c
+     cl -Ox -DCodonNSbranchsites -FeevolverNSbranchsites.exe evolver.c tools.c
 
-     cc -fast -o evolver evolver.c tools.c -lm
+     cc -O3 -o evolver evolver.c tools.c -lm
      cc -O4 -DCodonNSbranches -o evolverNSbranches evolver.c tools.c -lm
      cc -O4 -DCodonNSsites -o evolverNSsites evolver.c tools.c -lm
      cc -O4 -DCodonNSbranchsites -o evolverNSbranchsites evolver.c tools.c -lm
@@ -36,7 +36,7 @@
 
 struct CommonInfo {
    unsigned char *z[2 * NS - 1];
-   char spname[NS][LSPNAME + 1], daafile[512], cleandata, readpattern;
+   unsigned char *spname[NS], daafile[512], cleandata, readpattern;
    int ns, ls, npatt, np, ntime, ncode, clock, rooted, model, icode;
    int seqtype, *pose, ncatG, NSsites;
    int ngene, lgene[1], posG[1 + 1];  /* not used */
@@ -52,8 +52,8 @@ struct TREEB {
 }  tree;
 struct TREEN {
    int father, nson, sons[MAXNSONS], ibranch;
-   double branch, age, omega, label, *conP;
-   char *nodeStr, fossil;
+   double branch, age, omega, label, label2, *conP;
+   char *annotation, fossil;
 }  *nodes;
 
 extern char BASEs[];
@@ -107,7 +107,7 @@ int main(int argc, char*argv[])
    if (argc == 1)
       printf("Results for options 1-4 & 8 go into %s\n", outf);
    else if (option != 5 && option != 6 && option != 7 && option != 9) {
-      puts("Usage: \n\tevolver \n\tevolver option# MyDataFile"); exit(-1);
+      puts("Usage: \n\tevolver \n\tevolver option datafile"); exit(-1);
    }
    if (option >= 4 && option <= 6)
       MCctlf = argv[2];
@@ -117,6 +117,11 @@ int main(int argc, char*argv[])
       if (argc > 4) sscanf(argv[4], "%d", &pick1tree);
    }
 
+   for (i = 0; i < NS; i++) {
+      com.spname[i] = (char*)malloc(LSPNAME * sizeof(char*));
+      if (com.spname[i] == NULL) error2("oom");
+   }
+
 #if defined (CodonNSbranches)
    option = 6;  com.model = 1;
    MCctlf = (argc == 3 ? argv[2] : "MCcodonNSbranches.dat");
@@ -219,8 +224,9 @@ int main(int argc, char*argv[])
       Simulate(MCctlf ? MCctlf : MCctlf0[option - 5]);
    }
    else if (option == 9) {
+
       CladeSupport(fout, treefile, 1, mastertreefile, pick1tree);
-      /* CladeMrBayesProbabilities("/papers/BPPJC3sB/Karol.trees"); */
+      /* CladeMrBayesProbabilities("\data\Lakner2008SB/M1510YY03.nex.t"); */
    }
    return(0);
 }
@@ -271,12 +277,12 @@ void LabelClades(FILE *fout)
 
    i = (com.ns * 2 - 1) * sizeof(struct TREEN);
    if ((nodes = (struct TREEN*)malloc(i)) == NULL) error2("oom");
-   for (i = 0; i < com.ns * 2 - 1; i++)  nodes[i].nodeStr = NULL;
+   for (i = 0; i < com.ns * 2 - 1; i++)  nodes[i].annotation = NULL;
    for (i = 0; i < com.ns - 1; i++) {
       anc[i] = (int*)malloc((com.ns / SI + 1) * sizeof(int));
       if (anc[i] == NULL)  error2("oom");
    }
-   ReadTreeN(ftree, &haslength, &j, 1, 0);
+   ReadTreeN(ftree, &haslength, 1, 0);
    fclose(ftree);
    if (debug) { OutTreeN(F0, 1, PrNodeNum);  FPN(F0); }
 
@@ -288,7 +294,6 @@ void LabelClades(FILE *fout)
       for (i = 0; i < com.ns; i++)
          chosen[i] = '\0';
       
-
       k = strlen(key);
       for (i = 0; i < com.ns; i++) {
          if ((p = strstr(com.spname[i], key))
@@ -400,7 +405,7 @@ void TreeDistanceDistribution(FILE* fout)
    puts("\ntree #: mean prop of tree pairs with 0 1 2 ... shared bipartitions\n");
    fputs("\ntree #: prop of tree pairs with 0 1 2 ... shared bipartitions\n", fout);
    for (i = 0; i < ntree; i++) {
-      ReadTreeN(ftree, &j, &k, 0, 1);
+      ReadTreeN(ftree, &j, 0, 1);
       nib[i] = tree.nbranch - com.ns;
       Tree2Partition(partition + i*lpart);
    }
@@ -462,7 +467,7 @@ void TreeDistances(FILE* fout)
    if (k == 2) {    /* pairwise comparisons */
       fputs("Number of identical bi-partitions in pairwise comparisons\n", fout);
       for (i = 0; i < ntree; i++) {
-         ReadTreeN(ftree, &j, &k, 0, 1);
+         ReadTreeN(ftree, &j, 0, 1);
          nib[i] = tree.nbranch - com.ns;
          Tree2Partition(partition + i*lpart);
       }
@@ -477,7 +482,7 @@ void TreeDistances(FILE* fout)
       }
    }
    else {  /* first vs. others */
-      ReadTreeN(ftree, &j, &k, 0, 1);
+      ReadTreeN(ftree, &j, 0, 1);
       nib[0] = tree.nbranch - com.ns;
       if (nib[0] == 0) error2("1st tree is a star tree..");
       Tree2Partition(partition);
@@ -496,7 +501,7 @@ void TreeDistances(FILE* fout)
       fputs("\nCorrect internal branches compared with the 1st tree:\n", fout);
       FOR(k, nib[0]) nIBsame[k] = 0;
       for (i = 1, mp = vp = 0; i < ntree; i++, FPN(fout)) {
-         ReadTreeN(ftree, &j, &k, 0, 1);
+         ReadTreeN(ftree, &j, 0, 1);
          nib[1] = tree.nbranch - com.ns;
          Tree2Partition(partition + lpart);
          nsame = NSameBranch(partition, partition + lpart, nib[0], nib[1], IBsame);
@@ -823,7 +828,7 @@ void Simulate(char *ctlf)
 
    if (fixtree) {
       fscanf(fin, "%lf", &tlength);   fgets(line, lline, fin);
-      if (ReadTreeN(fin, &i, &j, 1, 1))  /* might overwrite spname */
+      if (ReadTreeN(fin, &i, 1, 1))  /* might overwrite spname */
          error2("err tree..");
 
       if (i == 0) error2("use : to specify branch lengths in tree");
@@ -897,9 +902,9 @@ void Simulate(char *ctlf)
          for (i = 0; i < tree.nnode; i++)
             blengthBS[i] = nodes[i].branch;
          for (k = 0; k < com.ncatG; k++) {
-            ReadTreeN(fin, &i, &j, 0, 1);
+            ReadTreeN(fin, &i, 0, 1);
             if (i) error2("do not include branch lengths except in the first tree.");
-            if (!j) error2("Use # to specify omega's for branches");
+            /* if (!j) error2("Use # to specify omega's for branches"); */
             for (i = 0; i < tree.nnode; i++)  com.omegaBS[i*com.ncatG + k] = nodes[i].label;
          }
          for (i = 0; i < tree.nnode; i++)
@@ -907,7 +912,7 @@ void Simulate(char *ctlf)
             nodes[i].branch = blengthBS[i];  nodes[i].label = nodes[i].omega = 0;
          }
          for (i = 0; i < tree.nnode; i++) {  /* print out omega as node labels. */
-            nodes[i].nodeStr = pc = (char*)malloc(20 * com.ncatG * sizeof(char));
+            nodes[i].annotation = pc = (char*)malloc(20 * com.ncatG * sizeof(char));
             sprintf(pc, "'[%.2f", com.omegaBS[i*com.ncatG + 0]);
             for (k = 1, pc += strlen(pc); k < com.ncatG; k++, pc += strlen(pc))
                sprintf(pc, ", %.2f", com.omegaBS[i*com.ncatG + k]);
@@ -1100,8 +1105,7 @@ void Simulate(char *ctlf)
             PatternWeightJC69like();
       }
       if (format == 2 || format == 3) fprintf(fseq, "\n\n[Replicate # %d]\n", ir + 1);
-      printSeqs(fseq, NULL, NULL, format); /* printsma not usable as it codes into 0,1,...,60. */
-
+      printSeqs(fseq, com.z, com.spname, com.ns, com.ls, com.npatt, com.fpatt, NULL, NULL, format); /* printsma not usable as it codes into 0,1,...,60. */
       if ((format == 2 || format == 3) && !fixtree) {
          fprintf(fseq, "\nbegin tree;\n   tree true_tree = [&U] ");
          OutTreeN(fseq, 1, 1); fputs(";\n", fseq);
@@ -1169,7 +1173,7 @@ void Simulate(char *ctlf)
    for (j = 0; j < com.ns * 2 - 1; j++) free(com.z[j]);
    free(space);
    if (com.model && com.NSsites) /* branch-site model */
-      for (i = 0; i < tree.nnode; i++)  free(nodes[i].nodeStr);
+      for (i = 0; i < tree.nnode; i++)  free(nodes[i].annotation);
    free(nodes);
    if (com.alpha || com.ncatG) {
       free(com.siterates);  com.siterates = NULL;
@@ -1243,8 +1247,7 @@ void CladeMrBayesProbabilities(char treefile[])
    6 -- ...........................*************   8001 1.000 0.005 (0.000)
    7 -- ....................********************   8001 1.000 0.006 (0.000)
 
-      Note 4 Jan 2014: This uses parti2B[], and is broken after i rewrote
-      Tree2Partition().
+   Note 4 Jan 2014: This uses parti2B[], and is broken after i rewrote Tree2Partition().
 */
    int lline = 100000, i, j, k, nib, inode, parti2B[NS];
    char line[100000], *partition, *p;
@@ -1266,13 +1269,13 @@ void CladeMrBayesProbabilities(char treefile[])
    if (i && i != com.ns) error2("do you mean to specify ns in the tree file?");
    i = (com.ns * 2 - 1) * sizeof(struct TREEN);
    if ((nodes = (struct TREEN*)malloc(i)) == NULL) error2("oom");
-   ReadTreeN(ftree, &i, &j, 0, 1);
+   ReadTreeN(ftree, &i, 0, 1);
 
    FPN(F0);  OutTreeN(F0, 0, 0);  FPN(F0);  FPN(F0);
    nib = tree.nbranch - com.ns;
    for (i = 0; i < tree.nnode; i++) {
-      nodes[i].nodeStr = NULL;
-      if (i > com.ns) nodes[i].nodeStr = (char*)malloc(100 * sizeof(char));
+      nodes[i].annotation = NULL;
+      if (i > com.ns) nodes[i].annotation = (char*)malloc(100 * sizeof(char));
    }
 
    partition = (char*)malloc(nib*com.ns * sizeof(char));
@@ -1294,19 +1297,19 @@ void CladeMrBayesProbabilities(char treefile[])
       for (k = 0; k < nmbfiles; k++) {
          if (GrepLine(fmb[k], cladestr, line, lline)) {
             p = strstr(line, cladestr);
-            sscanf(p + com.ns, "%lf%lf\0", &t, &Pclade[i*nmbfiles + k]);
+            sscanf(p + com.ns, "%lf%lf", &t, &Pclade[i*nmbfiles + k]);
          }
       }
       for (k = 0; k < nmbfiles; k++) printf("%6.2f", Pclade[i*nmbfiles + k]);
       FPN(F0);
-      for (k = 0, p = nodes[inode].nodeStr; k < nmbfiles; k++) {
+      for (k = 0, p = nodes[inode].annotation; k < nmbfiles; k++) {
          sprintf(p, "%3.0f%s", Pclade[i*nmbfiles + k] * 100, (k < nmbfiles - 1 ? "/" : ""));
          p += 4;
       }
    }
    FPN(F0);  OutTreeN(F0, 1, PrLabel);  FPN(F0);
 
-   for (i = 0; i < tree.nnode; i++) free(nodes[i].nodeStr);
+   for (i = 0; i < tree.nnode; i++) free(nodes[i].annotation);
    free(nodes); free(partition);  free(Pclade);
    fclose(ftree);
    for (k = 0; k < nmbfiles; k++) fclose(fmb[k]);

src/paml.h

@@ -13,6 +13,7 @@
 #include <limits.h>
 #include <float.h>
 #include <time.h>
+#include <assert.h>
 
 #define square(a) ((a)*(a))
 #define FOR(i,n) for(i=0; i<n; i++)
@@ -138,7 +139,7 @@ int f_mono_di (FILE *fout, char z[], int ls, int iring, double fb1[], double fb2
 int PickExtreme (FILE *fout, char z[], int ls, int iring, int lfrag, int ffrag[]);
 
 int print1seq (FILE*fout, unsigned char *z, int ls, int pose[]);
-void printSeqs(FILE *fout, int *pose, char keep[], int format);
+void printSeqs(FILE *fout, unsigned char *z[], unsigned char *spnames[], int ns, int ls, int npatt, double fpatt[], int *pose, char keep[], int format);
 int printPatterns(FILE *fout);
 void printSeqsMgenes (void);
 int printsma (FILE*fout, char*spname[], unsigned char*z[], int ns, int l, int lline, int gap, int seqtype, 
@@ -150,10 +151,8 @@ double SeqDivergence (double x[], int model, double alpha, double *kapa);
 int symtest (double freq[], int n, int nobs, double space[], double *chisym, 
     double* chihom);
 int dSdNNG1986(char *z1, char *z2, int lc, int icode, int transfed, double *dS, double *dN, double *Ssites, double *Nsites);
-int difcodonNG(char codon1[], char codon2[], double *SynSite,double *AsynSite, 
-    double *SynDif, double *AsynDif, int transfed, int icode);
-int difcodonLWL85 (char codon1[], char codon2[], double sites[3], double sdiff[3], 
-                   double vdiff[3], int transfed, int icode);
+int difcodonNG(char codon1[], char codon2[], double *SynSite,double *AsynSite, double *SynDif, double *AsynDif, int transfed, int icode);
+int difcodonLWL85 (char codon1[], char codon2[], double sites[3], double sdiff[3], double vdiff[3], int transfed, int icode);
 int testTransP (double P[], int n);
 double testDetailedBalance (double P[], double pi[], int n);
 void pijJC69 (double pij[2], double t);
@@ -163,12 +162,10 @@ int PMatTN93 (double P[], double a1t, double a2t, double bt, double pi[]);
 int PMatUVRoot (double P[],double t,int n,double U[],double V[],double Root[]);
 int PMatCijk (double PMat[], double t);
 int PMatQRev(double P[], double pi[], double t, int n, double space[]);
-int EvolveHKY85 (char source[], char target[], int ls, double t, 
-    double rates[], double pi[], double kapa, int isHKY85);
+int EvolveHKY85 (char source[], char target[], int ls, double t, double rates[], double pi[], double kapa, int isHKY85);
 double DistanceIJ (int is, int js, int model, double alpha, double *kappa);
 int DistanceMatNuc (FILE *fout, FILE*f2base, int model, double alpha);
-int EigenQGTRbase (FILE* fout, double kappa[], double pi[], 
-                   int *nR, double Root[], double Cijk[]);
+int EigenQGTRbase (FILE* fout, double kappa[], double pi[], int *nR, double Root[], double Cijk[]);
 int  DistanceMatNG86 (FILE *fout, FILE*fds, FILE*fdn, FILE*dt, double alpha);
 int  setmark_61_64 (void);
 
@@ -236,12 +233,11 @@ int eigenRealSym(double A[], int n, double Root[], double work[]);
 
 int MeanVar (double x[], int n, double *mean, double *var);
 int variance (double x[], int n, int nx, double mx[], double vx[]);
-int correl (double x[], double y[], int n, double *mx, double *my,
-            double *vxx, double *vxy, double *vyy, double *r);
+int correl (double x[], double y[], int n, double *mx, double *my, double *vxx, double *vxy, double *vyy, double *r);
 int comparefloat  (const void *a, const void *b);
 int comparedouble (const void *a, const void *b);
-double Eff_IntegratedCorrelationTime (double x[], int n, double *mx, double *varx, double *rho1);
-double Eff_IntegratedCorrelationTime2(double x[], int n, int nbatch, double *mx, double *varx);
+double Eff_IntegratedCorrelationTime (double x[], int n, double *mx, double *vx, double *rho1);
+double Eff_IntegratedCorrelationTime2(double x[], int n, int nbatch, double *mx, double *vx);
 
 int HPDinterval(double x[], int n, double HPD[2], double alpha);
 int DescriptiveStatistics(FILE *fout, char infile[], int nbin, int propternary, int SkipColumns);
@@ -249,12 +245,9 @@ int DescriptiveStatisticsSimple (FILE *fout, char infile[], int SkipColumns);
 int splitline(char line[], int nfields, int fields[]);
 int scanfile (FILE*fin, int *nrecords, int *nx, int *HasHeader, char line[], int ifields[]);
 
-double bound (int nx, double x0[], double p[], double x[],
-    int (*testx) (double x[], int nx));
-int gradient (int n, double x[], double f0, double g[], 
-    double (* fun)(double x[],int n), double space[], int Central);
-int Hessian (int nx, double x[], double f, double g[], double H[],
-    double (*fun)(double x[], int n), double space[]);
+double bound (int nx, double x0[], double p[], double x[], int (*testx) (double x[], int nx));
+int gradient (int n, double x[], double f0, double g[], double (* fun)(double x[],int n), double space[], int Central);
+int Hessian (int nx, double x[], double f, double g[], double H[], double (*fun)(double x[], int n), double space[]);
 int HessianSKT2004 (double xmle[], double lnLm, double g[], double H[]);
 
 int H_end (double x0[], double x1[], double f0, double f1, double e1, double e2, int n);
@@ -264,7 +257,6 @@ double LineSearch2 (double(*fun)(double x[],int n), double *f, double x0[],
 
 void xtoFreq(double x[], double freq[], int n);
 
-
 int SetxBound (int np, double xb[][2]);
 int ming1 (FILE *fout, double *f, double (* fun)(double x[], int n),
     int (*dfun) (double x[], double *f, double dx[], int n),
@@ -288,19 +280,19 @@ int nls2 (FILE *fout, double *sx, double * x0, int nx,
     int (*testx) (double x[], int nx),
     int ny, double e);
 
-int ResetFinetuneSteps(FILE *fout, double Pjump[], double finetune[], int nsteps);
+int ResetStepLengths(FILE *fout, double Pjump[], double finetune[], int nsteps);
 
 /* tree structure functions in treesub.c */
 void NodeToBranch (void);
 void BranchToNode (void);
 void ClearNode (int inode);
-int ReadTreeN (FILE *ftree, int *haslength, int *haslabel, int copyname, int popline);
+int ReadTreeN (FILE *ftree, int *haslength, int copyname, int popline);
 int ReadTreeB (FILE *ftree, int popline);
 int OutTreeN (FILE *fout, int spnames, int printopt);
 int OutTreeB (FILE *fout);
 int DeRoot (void);
 int GetSubTreeN (int keep[], int space[]);
-void printtree (int timebranches);
+void PrintTree (int timebranches);
 void PointconPnodes (void);
 int SetBranch (double x[]);
 int DistanceMat (FILE *fout, int ischeme, double alfa, double *kapa);
@@ -317,8 +309,7 @@ int PopEmptyLines (FILE* fseq, int lline, char line[]);
 int hasbase (char *str);
 int blankline (char *str);
 
-void BranchLengthBD(int rooted, double birth, double death, double sample, 
-     double mut);
+void BranchLengthBD(int rooted, double birth, double death, double sample, double mut);
 int RandomLHistory (int rooted, double space[]);
 
 void Tree2Partition (char partition[]);
@@ -359,8 +350,7 @@ int GetSubSeqs(int nsnew);
 
 /* functions for evolving sequences */
 int GenerateSeq (void);
-int Rates4Sites (double rates[],double alfa,int ncatG,int ls, int cdf,
-    double space[]);
+int Rates4Sites (double rates[],double alfa,int ncatG,int ls, int cdf, double space[]);
 void Evolve (int inode);
 void EvolveJC (int inode);
 
@@ -396,8 +386,6 @@ enum {PrBranch=1, PrNodeNum=2, PrLabel=4, PrNodeStr=8, PrAge=16, PrOmega=32} Out
 
 #define PAML_RELEASE      0
 
-#define FullSeqNames      0   /* 1: numbers at the beginning of sequence name are part of name */
-
-#define pamlVerStr "paml version 4.9h, March 2018"
+#define pamlVerStr "paml version 4.9i, September 2018"
 
 #endif

src/pamp.c

@@ -32,10 +32,10 @@ double lfunAlpha_YK96 (double x);
 
 struct CommonInfo {
    unsigned char *z[NS];
-   char *spname[NS], seqf[256],outf[256],treef[256];
+   char *spname[NS], seqf[2048], outf[2048], treef[2048];
    int seqtype, ns, ls, ngene, posG[NGENE + 1], lgene[NGENE], *pose, npatt, readpattern;
    int np, ntime, ncode, fix_kappa, fix_rgene, fix_alpha, clock, model, ncatG, cleandata;
-   int print, nhomo;
+   int print, nhomo, verbose;
    double *fpatt, *conP;
    /* not used */
    double lmax, pi[NCODE], kappa, alpha, rou, rgene[NGENE], piG[NGENE][NCODE];
@@ -46,8 +46,8 @@ struct TREEB {
 }  tree;
 struct TREEN {
    int father, nson, sons[MAXNSONS], ibranch;
-   double branch, age, label, *conP;
-   char *nodeStr, fossil;
+   double branch, age, label, label2, *conP;
+   char *annotation, fossil;
 }  *nodes;
 
 
@@ -68,7 +68,7 @@ int LASTROUND=0; /* no use for this */
 int main(int argc, char *argv[])
 {
    FILE *ftree, *fout, *fseq;
-   char ctlf[32]="pamp.ctl";
+   char ctlf[2048] = "pamp.ctl";
    char *Seqstr[] = { "nucleotide", "", "amino-acid", "Binary" };
    int itree, ntree, i, j, s3;
    double *space, *Ft;
@@ -120,7 +120,7 @@ int main (int argc, char *argv[])
       printf("\nTREE # %2d\n", itree + 1);
       fprintf(fout, "\nTREE # %2d\n", itree + 1);
 
-      if (ReadTreeN (ftree, &i,&j, 0, 1)) error2 ("err tree..");
+      if (ReadTreeN(ftree, &i, 0, 1)) error2("err tree..");
       OutTreeN(F0, 0, 0);    FPN(F0);
       OutTreeN(fout, 0, 0);  FPN(fout);
 
@@ -178,7 +178,9 @@ int GetOptions (char *ctlf)
                }
             }
             if (iopt == nopt)
-            { printf ("\nopt %s in %s\n", opt, ctlf);  exit (-1); }
+            {
+               printf("\nopt %s in %s\n", opt, ctlf);  exit(-1);
+            }
       }
       fclose(fctl);
    }
@@ -442,7 +444,9 @@ int PathwayMP1 (FILE *fout, int *maxchange, int NSiteChange[],
                sumpr += pr;
                FOR(i, nid) pnsite[i*n + nodeb[i + com.ns]] += pr;
                if (pr > bestpr)
-                  { bestpr=pr; FOR(i,nid) bestPath[i]=PATHWay[i];}
+               {
+                  bestpr = pr; FOR(i, nid) bestPath[i] = PATHWay[i];
+               }
             }
             else {
                for (i = 0, nchange0 = 0; i < tree.nbranch; i++) {