CMakeLists.txt

@@ -51,7 +51,7 @@ add_definitions(-DIQ_TREE)
 # The version number.
 set (iqtree_VERSION_MAJOR 1)
 set (iqtree_VERSION_MINOR 6)
-set (iqtree_VERSION_PATCH "7.1")
+set (iqtree_VERSION_PATCH "8")
 
 set(BUILD_SHARED_LIBS OFF)
 

alignment/alignment.cpp

@@ -3177,7 +3177,17 @@ void Alignment::createBootstrapAlignment(int *pattern_freq, const char *spec, in
     if (Params::getInstance().jackknife_prop > 0.0 && spec)
         outError((string)"Unsupported jackknife with " + spec);
 
-    if (!spec) {
+    if (spec && strncmp(spec, "SCALE=", 6) == 0) {
+        // multi-scale bootstrapping called by AU test
+        int orig_nsite = nsite;
+        double scale = convert_double(spec+6);
+        nsite = (int)round(scale * nsite);
+        for (site = 0; site < nsite; site++) {
+            int site_id = random_int(orig_nsite, rstream);
+            int ptn_id = getPatternID(site_id);
+            pattern_freq[ptn_id]++;
+        }
+    } else if (!spec) {
 
         int nptn = getNPattern();
 
@@ -4337,6 +4347,12 @@ void Alignment::getPatternFreq(IntVector &freq) {
 		freq[cnt] = (*it).frequency;
 }
 
+void Alignment::getPatternFreq(int *freq) {
+    int cnt = 0;
+    for (iterator it = begin(); it < end(); it++, cnt++)
+        freq[cnt] = (*it).frequency;
+}
+
 //added by MA
 void Alignment::multinomialProb(Alignment refAlign, double &prob)
 {

alignment/alignment.h

@@ -332,6 +332,11 @@ public:
      */
     virtual void getPatternFreq(IntVector &freq);
 
+    /**
+     * @param[out] freq vector of site-pattern frequencies
+     */
+    virtual void getPatternFreq(int *freq);
+
     /**
             @param i sequence index
             @return sequence name

alignment/superalignment.cpp

@@ -776,7 +776,7 @@ void SuperAlignment::getSitePatternIndex(IntVector &pattern_index) {
 
 void SuperAlignment::getPatternFreq(IntVector &pattern_freq) {
 	ASSERT(isSuperAlignment());
-	int offset = 0;
+	size_t offset = 0;
 	if (!pattern_freq.empty()) pattern_freq.resize(0);
 	for (vector<Alignment*>::iterator it = partitions.begin(); it != partitions.end(); it++) {
 		IntVector freq;
@@ -786,6 +786,15 @@ void SuperAlignment::getPatternFreq(IntVector &pattern_freq) {
 	}
 }
 
+void SuperAlignment::getPatternFreq(int *pattern_freq) {
+    ASSERT(isSuperAlignment());
+    size_t offset = 0;
+    for (vector<Alignment*>::iterator it = partitions.begin(); it != partitions.end(); it++) {
+        (*it)->getPatternFreq(pattern_freq + offset);
+        offset += (*it)->getNPattern();
+    }
+}
+
 void SuperAlignment::printSiteInfo(const char* filename) {
     try {
         ofstream out(filename);

alignment/superalignment.h

@@ -114,6 +114,11 @@ public:
 	*/
 	virtual void getPatternFreq(IntVector &pattern_freq);
 
+    /**
+     * @param[out] freq vector of site-pattern frequencies
+     */
+    virtual void getPatternFreq(int *freq);
+
     /**
         Print all site information to a file
         @param filename output file name

debian/changelog

+iqtree (1.6.8+dfsg-1) unstable; urgency=medium
+
+  * New upstream version
+  * Autopkgtest: Upstream now requires -redo option to allow overriding
+    existing output data
+
+ -- Andreas Tille <tille@debian.org>  Thu, 13 Dec 2018 20:07:24 +0100
+
 iqtree (1.6.7.1+dfsg-1) unstable; urgency=medium
 
   * New upstream version

debian/tests/run-unit-test

@@ -8,10 +8,12 @@ else
 fi
 
 pkg=iqtree
-if [ "$ADTTMP" = "" ] ; then
-  ADTTMP=`mktemp -d /tmp/${pkg}-test.XXXXXX`
+if [ "$AUTOPKGTEST_TMP" = "" ] ; then
+  AUTOPKGTEST_TMP=`mktemp -d /tmp/${pkg}-test.XXXXXX`
+  trap "rm -rf $AUTOPKGTEST_TMP" 0 INT QUIT ABRT PIPE TERM
 fi
-cd $ADTTMP
+cd $AUTOPKGTEST_TMP
+
 cp -a $EXAMPLEDIR/example.phy* $EXAMPLEDIR/example.nex* .
 find . -name "*.gz" -exec gunzip \{\} \;
 
@@ -27,32 +29,33 @@ fi
 #  PATHTOEXE=`find $CURDIR -name iqtree -type l`
 #fi
 echo "Executing tests using $PATHTOEXE ..."
+set -x
 time $PATHTOEXE -s example.phy
-time $PATHTOEXE -s example.phy -pre myprefix
-time $PATHTOEXE -s example.phy -nni1
-time $PATHTOEXE -s example.phy -m TEST
-time $PATHTOEXE -s example.phy -m TIM+I+G
-time $PATHTOEXE -s example.phy -m TESTONLY
-time $PATHTOEXE -s example.phy -m TIM+I+G -bb 1000
-time $PATHTOEXE -s example.phy -m TIM+I+G -b 100
-time $PATHTOEXE -s example.phy -m TIM+I+G -alrt 1000
-time $PATHTOEXE -s example.phy -m TIM+I+G -lbp 1000
-time $PATHTOEXE -s example.phy -m TIM+I+G -alrt 1000 -lbp 1000
-time $PATHTOEXE -s example.phy -m TIM+I+G -bb 1000 -alrt 1000 -lbp 1000
-time $PATHTOEXE -s example.phy -sp example.nex
-time $PATHTOEXEOMP -s example.phy -omp 2
+time $PATHTOEXE -redo -s example.phy -pre myprefix
+time $PATHTOEXE -redo -s example.phy -nni1
+time $PATHTOEXE -redo -s example.phy -m TEST
+time $PATHTOEXE -redo -s example.phy -m TIM+I+G
+time $PATHTOEXE -redo -s example.phy -m TESTONLY
+time $PATHTOEXE -redo -s example.phy -m TIM+I+G -bb 1000
+time $PATHTOEXE -redo -s example.phy -m TIM+I+G -b 100
+time $PATHTOEXE -redo -s example.phy -m TIM+I+G -alrt 1000
+time $PATHTOEXE -redo -s example.phy -m TIM+I+G -lbp 1000
+time $PATHTOEXE -redo -s example.phy -m TIM+I+G -alrt 1000 -lbp 1000
+time $PATHTOEXE -redo -s example.phy -m TIM+I+G -bb 1000 -alrt 1000 -lbp 1000
+time $PATHTOEXE -redo -s example.phy -sp example.nex
+time $PATHTOEXEOMP -redo -s example.phy -omp 2
 if [ $(nproc) -ge 3 ] ; then
-  time $PATHTOEXEOMP -s example.phy -omp 3
+  time $PATHTOEXEOMP -redo -s example.phy -omp 3
 fi
 if [ -e example.treels ] ; then
-  time $PATHTOEXE -s example.phy -z example.treels
-  time $PATHTOEXE -s example.phy -z example.treels -n 1
-  time $PATHTOEXE -s example.phy -z example.treels -n 1 -zb 1000
-  time $PATHTOEXE -s example.phy -z example.treels -n 1 -zb 1000 -zw
+  time $PATHTOEXE -redo -s example.phy -z example.treels
+  time $PATHTOEXE -redo -s example.phy -z example.treels -n 1
+  time $PATHTOEXE -redo -s example.phy -z example.treels -n 1 -zb 1000
+  time $PATHTOEXE -redo -s example.phy -z example.treels -n 1 -zb 1000 -zw
 fi
-time $PATHTOEXE -s example.phy -m 010010+G
+time $PATHTOEXE -redo -s example.phy -m 010010+G
 if [ -e mymodel ] ; then
-  time $PATHTOEXE -s example.phy -m mymodel+G
+  time $PATHTOEXE -redo -s example.phy -m mymodel+G
 fi
-time $PATHTOEXE -s example.phy -m 'TN{2.0,3.0}+G8{0.5}+I{0.15}'
-time $PATHTOEXE -s example.phy -m GTR+G+Fo
+time $PATHTOEXE -redo -s example.phy -m 'TN{2.0,3.0}+G8{0.5}+I{0.15}'
+time $PATHTOEXE -redo -s example.phy -m GTR+G+Fo

gsl/mygsl.h

@@ -46,6 +46,13 @@ double gsl_ran_ugaussian_pdf (const double x);
 */
 double gsl_cdf_ugaussian_P (const double x);
 
+/*
+ 1.0 - cumulative distribution function for standard normal distribution
+ @param x x-value
+ @return 1.0-CDF at x
+ */
+double gsl_cdf_ugaussian_Q (const double x);
+
 /*
     quantile function for standard normal distribution (or CDF-inverse function)
     @param P probability value

main/phyloanalysis.cpp

@@ -1190,6 +1190,8 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
 			out << endl << "USER TREES" << endl << "----------" << endl << endl;
 			out << "See " << params.out_prefix << ".trees for trees with branch lengths." << endl << endl;
 			if (params.topotest_replicates && info.size() > 1) {
+                if (params.do_au_test && params.topotest_replicates < 10000)
+                    out << "WARNING: Too few replicates for AU test. At least -zb 10000 for reliable results!" << endl << endl;
                 out << "Tree      logL    deltaL  bp-RELL    p-KH     p-SH    ";
 				if (params.do_weighted_test)
 					out << "p-WKH    p-WSH    ";
@@ -1219,17 +1221,19 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
 					out << " = tree " << distinct_trees[orig_id]+1 << endl;
 					continue;
 				}
-				out.precision(3);
+                out.unsetf(ios::fixed);
+				out.precision(10);
 				out.width(12);
 				out << info[tid].logl << " ";
 				out.width(7);
+                out.precision(5);
 				out << maxL - info[tid].logl;
 				if (!params.topotest_replicates || info.size() <= 1) {
 					out << endl;
 					tid++;
 					continue;
 				}
-				out.precision(4);
+				out.precision(3);
 				out << "  ";
 				out.width(6);
 				out << info[tid].rell_bp;
@@ -1264,21 +1268,22 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
 					else
 						out << " + ";
 				}
-				out.width(6);
-				out << info[tid].elw_value;
+				out.width(9);
+				out << right << info[tid].elw_value;
 				if (info[tid].elw_confident)
 					out << " + ";
 				else
 					out << " - ";
 
                 if (params.do_au_test) {
-                    out.width(6);
+                    out.width(8);
                     out << right << info[tid].au_pvalue;
                     if (info[tid].au_pvalue < 0.05)
                         out << " - ";
                     else
                         out << " + ";
                 }
+                out.setf(ios::fixed);
 
 				out << endl;
 				tid++;
@@ -1930,7 +1935,7 @@ void startTreeReconstruction(Params &params, IQTree* &iqtree, ModelCheckpoint &m
 			PhyloSuperTree *stree = (PhyloSuperTree*)iqtree;
 			for (PhyloSuperTree::iterator it = stree->begin(); it != stree->end(); it++)
 				if ((*it)->aln->seq_type != SEQ_DNA && (*it)->aln->seq_type != SEQ_PROTEIN)
-					params.start_tree = STT_BIONJ;
+					params.start_tree = STT_PARSIMONY;
 		} else if (iqtree->aln->seq_type != SEQ_DNA && iqtree->aln->seq_type != SEQ_PROTEIN)
 			params.start_tree = STT_PARSIMONY;
     }

main/phylotesting.cpp

@@ -1418,9 +1418,9 @@ string testOneModel(string &model_name, Params &params, Alignment *in_aln,
     else
         iqtree = new IQTree(in_aln);
     iqtree->setParams(&params);
-    iqtree->sse = params.SSE;
+    iqtree->setLikelihoodKernel(params.SSE);
     iqtree->optimize_by_newton = params.optimize_by_newton;
-    iqtree->num_threads = num_threads;
+    iqtree->setNumThreads(num_threads);
 
     iqtree->setCheckpoint(&model_info);
     iqtree->restoreCheckpoint();
@@ -1729,8 +1729,9 @@ void testPartitionModel(Params &params, PhyloSuperTree* in_tree, ModelCheckpoint
 	dfvec.resize(in_tree->size());
 	lenvec.resize(in_tree->size());
 
-    double *dist = new double[in_tree->size()*(in_tree->size()-1)/2];
-    pair<int,int> *distID = new pair<int,int>[in_tree->size()*(in_tree->size()-1)/2];
+    // 2018-10-23: +1 to avoid crash when size <= 2
+    double *dist = new double[in_tree->size()*(in_tree->size()-1)/2+1];
+    pair<int,int> *distID = new pair<int,int>[in_tree->size()*(in_tree->size()-1)/2+1];
     
     // sort partition by computational cost for OpenMP effciency
 	for (i = 0; i < in_tree->size(); i++) {
@@ -2671,7 +2672,7 @@ int mleloopmax=30;
 double mleeps=1e-10;
 int mlecoef(double *cnts, double *rr, double bb, int kk,
 	    double *coef0, /* set initinal value (size=2) */
-	    double *lrt, double *df, /* LRT statistic */
+	    double *lrt, int *df, /* LRT statistic */
         double *se
 	    )
 {
@@ -2727,22 +2728,22 @@ int mlecoef(double *cnts, double *rr, double bb, int kk,
   }
 
   /* calc log-likelihood */
-  *lrt=0.0; *df=0.0;
+  *lrt=0.0; *df=0;
   for(i=0;i<m;i++) {
     if(p[i]>0.0 && p[i]<1.0) {
-      *df+=1.0;
+      *df+=1;
       if(c[i]>0.0) a=c[i]*log(c[i]/b[i]/p[i]); else a=0.0;
       if(c[i]<b[i]) a+=(b[i]-c[i])*(log3(p[i])-log3(c[i]/b[i]));
       *lrt += a;
     }
   }
-  *lrt *= 2.0; *df -= 2.0;
+  *lrt *= 2.0; *df -= 2;
 
   /* write back the results */
   coef0[0]=coef[0]; coef0[1]=coef[1];
   *se = v11 + v22 - 2*v12;
 //  vmat[0][0]=v11;vmat[0][1]=vmat[1][0]=v12;vmat[1][1]=v22; 
-  if(loop==mleloopmax || *df< -0.01) i=1; else i=0;
+  if(loop==mleloopmax || *df< 0) i=1; else i=0;
   return i;
 }
 
@@ -2807,6 +2808,10 @@ void performAUTest(Params &params, PhyloTree *tree, double *pattern_lhs, vector<
     for (k = 0; k < nscales; k++) {
         string str = "SCALE=" + convertDoubleToString(r[k]);    
 		for (boot = 0; boot < nboot; boot++) {
+            if (r[k] == 1.0 && boot == 0)
+                // 2018-10-23: get one of the bootstrap sample as the original alignment
+                tree->aln->getPatternFreq(boot_sample);
+            else
                 tree->aln->createBootstrapAlignment(boot_sample, str.c_str(), rstream);
             for (ptn = 0; ptn < maxnptn; ptn++)
                 boot_sample_dbl[ptn] = boot_sample[ptn];
@@ -2843,14 +2848,14 @@ void performAUTest(Params &params, PhyloTree *tree, double *pattern_lhs, vector<
                 else
                     treelhs[(tid*nscales+k)*nboot + boot] = second_max_lh - max_lh;
 //            bp[k*ntrees+max_tid] += nboot_inv;
-        }
+        } // for boot
         
         // sort the replicates
         for (tid = 0; tid < ntrees; tid++) {
             quicksort<double,int>(treelhs + (tid*nscales+k)*nboot, 0, nboot-1);
         }
         
-    }
+    } // for scale
 
     aligned_free(boot_sample_dbl);
     aligned_free(boot_sample);
@@ -2914,6 +2919,15 @@ void performAUTest(Params &params, PhyloTree *tree, double *pattern_lhs, vector<
                 // first obtain d and c by weighted least square
                 doWeightedLeastSquare(nscales, w, rr, rr_inv, cc, d, c, se);
 
+                // maximum likelhood fit
+                double coef0[2] = {d, c};
+                int mlefail = mlecoef(this_bp, r, nboot, nscales, coef0, &rss, &df, &se);
+                
+                if (!mlefail) {
+                    d = coef0[0];
+                    c = coef0[1];
+                }
+
                 se = gsl_ran_ugaussian_pdf(d-c)*sqrt(se);
                 
                 // second, perform MLE estimate of d and c
@@ -2923,7 +2937,7 @@ void performAUTest(Params &params, PhyloTree *tree, double *pattern_lhs, vector<
     //            c = mle.c;
 
                 /* STEP 4: compute p-value according to Eq. 11 */
-                pval = 1.0 - gsl_cdf_ugaussian_P(d-c);
+                pval = gsl_cdf_ugaussian_Q(d-c);
                 z = -pval;
                 ze = se;
                 // compute sum of squared difference
@@ -2935,10 +2949,10 @@ void performAUTest(Params &params, PhyloTree *tree, double *pattern_lhs, vector<
                 
             } else {
                 // not enough data for WLS
-                double sum = 0.0;
+                int num0 = 0;
                 for (k = 0; k < nscales; k++)
-                    sum += cc[k];
-                if (sum >= 0.0) 
+                    if (this_bp[k] <= 0.0) num0++;
+                if (num0 > nscales/2)
                     pval = 0.0;
                 else
                     pval = 1.0;
@@ -2947,20 +2961,10 @@ void performAUTest(Params &params, PhyloTree *tree, double *pattern_lhs, vector<
                 rss = 0.0;
                 if (verbose_mode >= VB_MED)
                     cout << "   error in wls" << endl;
+                //info[tid].au_pvalue = pval;
+                //break;
             }
 
-            // maximum likelhood fit
-//            double coef0[2] = {d, c};
-//            double df;
-//            int mlefail = mlecoef(this_bp, r, nboot, nscales, coef0, &rss, &df, &se);
-//            
-//            if (!mlefail) {
-//                d = coef0[0];
-//                c = coef0[1];
-//                pval = 1.0 - gsl_cdf_ugaussian_P(d-c);
-//                z = -pval;
-//                ze = se;
-//            }
             
             if (verbose_mode >= VB_MED) {
                 cout.unsetf(ios::fixed);
@@ -2986,9 +2990,9 @@ void performAUTest(Params &params, PhyloTree *tree, double *pattern_lhs, vector<
             thp=x; 
             z0=z;
             ze0=ze;
-            idf0 = nscales-2;
+            idf0 = df;
             if(fabs(x-th)<1e-10) break;
-        }
+        } // for step
         
         if (failed && verbose_mode >= VB_MED)
             cout << "   degenerated" << endl;
@@ -3100,6 +3104,9 @@ void evaluateTrees(Params &params, IQTree *tree, vector<TreeInfo> &info, IntVect
         int *rstream = randstream;
 #endif
 		for (boot = 0; boot < params.topotest_replicates; boot++)
+            if (boot == 0)
+                tree->aln->getPatternFreq(boot_samples + (boot*nptn));
+            else
                 tree->aln->createBootstrapAlignment(boot_samples + (boot*nptn), params.bootstrap_spec, rstream);
 #ifdef _OPENMP
         finish_random(rstream);

model/modelmarkov.cpp

@@ -1001,7 +1001,8 @@ void ModelMarkov::readRates(istream &in) throw(const char*, string) {
 	if (str == "equalrate") {
 		for (int i = 0; i < nrates; i++)
 			rates[i] = 1.0;
-	} else {
+	} else if (is_reversible ){
+        // reversible model
 		try {
 			rates[0] = convert_double(str.c_str());
 		} catch (string &str) {
@@ -1015,6 +1016,36 @@ void ModelMarkov::readRates(istream &in) throw(const char*, string) {
 			if (rates[i] < 0.0)
 				throw "Negative rates not allowed";
 		}
+    } else {
+        // non-reversible model, read the whole rate matrix
+        int i = 0, row, col;
+        for (row = 0; row < num_states; row++) {
+            double row_sum = 0.0;
+            for (col = 0; col < num_states; col++)
+                if (row == 0 && col == 0) {
+                    // top-left element was already red
+                    try {
+                        row_sum = convert_double(str.c_str());
+                    } catch (string &str) {
+                        outError(str);
+                    }
+                } else if (row != col) {
+                    // non-diagonal element
+                    if (!(in >> rates[i]))
+                        throw name+string(": Rate entries could not be read");
+                    if (rates[i] < 0.0)
+                        throw "Negative rates found";
+                    row_sum += rates[i];
+                    i++;
+                } else {
+                    // diagonal element
+                    double d;
+                    in >> d;
+                    row_sum += d;
+                }
+            if (fabs(row_sum) > 1e-3)
+                throw "Row " + convertIntToString(row) + " does not sum to 0";
+        }
     }
 }
 

pda/hashsplitset.cpp

@@ -41,7 +41,8 @@ Split *SplitIntMap::findSplit(Split *sp, int &value) {
 
 int SplitIntMap::getValue(Split *sp) {
     int value;
-    ASSERT(findSplit(sp, value));
+    Split* findsp = findSplit(sp, value);
+    ASSERT(findsp);
     return value;
 }
 

pda/splitgraph.cpp

@@ -192,7 +192,8 @@ void SplitGraph::restoreCheckpoint() {
     for (int split = 0; split < nsplits; split++) {
         checkpoint->addListElement();
         string str;
-        ASSERT(checkpoint->getString("", str));
+        bool found = checkpoint->getString("", str);
+        ASSERT(found);
         stringstream ss(str);
         double weight;
         ss >> weight;

tree/mtree.cpp

@@ -842,7 +842,7 @@ void MTree::parseFile(istream &infile, char &ch, Node* &root, DoubleVector &bran
         renameString(seqname);
 //    seqname[seqlen] = 0;
     if (seqlen == 0 && root->isLeaf())
-        throw "A taxon has no name.";
+        throw "Redundant double-bracket ‘((…))’ with closing bracket ending at";
     if (seqlen > 0)
         root->name.append(seqname);
     if (root->isLeaf()) {

tree/phylokernel.h

@@ -1665,8 +1665,10 @@ void PhyloTree::computePartialParsimonyFastSIMD(PhyloNeighbor *dad_branch, Phylo
         // add dummy states
         if (site > 0 && site < NUM_BITS) {
             x += site/UINT_BITS;
+            if (site % UINT_BITS != 0) {
                 *x |= ~((1<<(site%UINT_BITS)) - 1);
                 x++;
+            }
             int max_sites = ((site+UINT_BITS-1)/UINT_BITS);
             memset(x, 255, (VCSIZE - max_sites)*sizeof(UINT));
         }

utils/tools.cpp

@@ -3505,8 +3505,8 @@ void parseArg(int argc, char *argv[], Params &params) {
 #endif
     }
 
-    if (params.do_au_test)
-        outError("The AU test is temporarily disabled due to numerical issue when bp-RELL=0");
+//    if (params.do_au_test)
+//        outError("The AU test is temporarily disabled due to numerical issue when bp-RELL=0");
     
     if (params.model_test_and_tree && params.partition_type != BRLEN_OPTIMIZE)
         outError("-mtree not allowed with edge-linked partition model (-spp or -q)");