[med-svn] [Git][med-team/iqtree][master] 7 commits: New upstream version 1.6.8+dfsg
Andreas Tille
gitlab at salsa.debian.org
Thu Dec 13 20:54:06 GMT 2018
Andreas Tille pushed to branch master at Debian Med / iqtree
Commits:
3ef1bdb2 by Andreas Tille at 2018-12-13T14:08:11Z
New upstream version 1.6.8+dfsg
- - - - -
a4ad90bf by Andreas Tille at 2018-12-13T14:08:16Z
Update upstream source from tag 'upstream/1.6.8+dfsg'
Update to upstream version '1.6.8+dfsg'
with Debian dir 59ea322a4df6c49bbb29b21873c01dacf2c91219
- - - - -
099d62f2 by Andreas Tille at 2018-12-13T14:08:16Z
New upstream version
- - - - -
c7f7944b by Andreas Tille at 2018-12-13T15:34:12Z
Autopkgtest: Upstream now requires -redo option to allow overriding existing output data
- - - - -
4142649b by Andreas Tille at 2018-12-13T17:17:41Z
Use AUTOPKGTEST_TMP instead of ADTTMP
- - - - -
7e2b7283 by Andreas Tille at 2018-12-13T18:48:03Z
Add one missing -redo
- - - - -
31fe0ddb by Andreas Tille at 2018-12-13T19:25:38Z
Upload to unstable
- - - - -
17 changed files:
- CMakeLists.txt
- alignment/alignment.cpp
- alignment/alignment.h
- alignment/superalignment.cpp
- alignment/superalignment.h
- debian/changelog
- debian/tests/run-unit-test
- gsl/mygsl.h
- main/phyloanalysis.cpp
- main/phylotesting.cpp
- model/modelmarkov.cpp
- pda/hashsplitset.cpp
- pda/splitgraph.cpp
- tree/mtree.cpp
- tree/phylokernel.h
- tree/phylotree.cpp
- utils/tools.cpp
Changes:
=====================================
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)
@@ -70,7 +70,7 @@ if (CMAKE_BUILD_TYPE STREQUAL "Release")
endif()
if (CMAKE_GENERATOR MATCHES "Xcode")
- set(CMAKE_XCODE_ATTRIBUTE_DEBUG_INFORMATION_FORMAT "dwarf-with-dsym")
+ set(CMAKE_XCODE_ATTRIBUTE_DEBUG_INFORMATION_FORMAT "dwarf-with-dsym")
set(CMAKE_XCODE_ATTRIBUTE_COMPILER_INDEX_STORE_ENABLE "No")
endif()
@@ -279,7 +279,7 @@ if (NOT IQTREE_FLAGS MATCHES "single")
elseif (CLANG)
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -pthread")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp=libomp -pthread")
- set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -fopenmp=libomp")
+ set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -fopenmp=libomp")
endif()
else()
message("OpenMP : NONE")
@@ -437,8 +437,8 @@ check_function_exists (gettimeofday HAVE_GETTIMEOFDAY)
check_function_exists (getrusage HAVE_GETRUSAGE)
check_function_exists (GlobalMemoryStatusEx HAVE_GLOBALMEMORYSTATUSEX)
check_function_exists (strndup HAVE_STRNDUP)
-check_function_exists (strtok_r HAVE_STRTOK_R)
-
+check_function_exists (strtok_r HAVE_STRTOK_R)
+
find_package(Backtrace)
# configure a header file to pass some of the CMake settings
@@ -512,21 +512,21 @@ endif()
if (IQTREE_FLAGS MATCHES "mpi")
add_library(mympi utils/TreeCollection.cpp utils/ObjectStream.cpp)
endif()
-
-if (NOT IQTREE_FLAGS MATCHES "single")
- if (BINARY32)
- link_directories(${PROJECT_SOURCE_DIR}/lib32)
- else()
- link_directories(${PROJECT_SOURCE_DIR}/lib)
- endif()
-endif()
+
+if (NOT IQTREE_FLAGS MATCHES "single")
+ if (BINARY32)
+ link_directories(${PROJECT_SOURCE_DIR}/lib32)
+ else()
+ link_directories(${PROJECT_SOURCE_DIR}/lib)
+ endif()
+endif()
add_executable(iqtree
main/main.cpp
main/phyloanalysis.cpp
main/phyloanalysis.h
main/phylotesting.cpp
-main/phylotesting.h
+main/phylotesting.h
obsolete/parsmultistate.cpp
)
=====================================
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
=====================================
@@ -1,3 +1,11 @@
+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 at 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,12 +1190,14 @@ void reportPhyloAnalysis(Params ¶ms, 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 ";
- out << "c-ELW";
+ out << " c-ELW";
if (params.do_au_test)
- out << " p-AU";
+ out << " p-AU";
out << endl << "------------------------------------------------------------------";
if (params.do_weighted_test)
@@ -1219,17 +1221,19 @@ void reportPhyloAnalysis(Params ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, Alignment *in_aln,
else
iqtree = new IQTree(in_aln);
iqtree->setParams(¶ms);
- 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 ¶ms, 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,7 +2808,11 @@ void performAUTest(Params ¶ms, PhyloTree *tree, double *pattern_lhs, vector<
for (k = 0; k < nscales; k++) {
string str = "SCALE=" + convertDoubleToString(r[k]);
for (boot = 0; boot < nboot; boot++) {
- tree->aln->createBootstrapAlignment(boot_sample, str.c_str(), rstream);
+ 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];
double max_lh = -DBL_MAX, second_max_lh = -DBL_MAX;
@@ -2843,14 +2848,14 @@ void performAUTest(Params ¶ms, 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);
@@ -2913,7 +2918,16 @@ void performAUTest(Params ¶ms, PhyloTree *tree, double *pattern_lhs, vector<
if (num_k >= 2) {
// 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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,7 +3104,10 @@ void evaluateTrees(Params ¶ms, IQTree *tree, vector<TreeInfo> &info, IntVect
int *rstream = randstream;
#endif
for (boot = 0; boot < params.topotest_replicates; boot++)
- tree->aln->createBootstrapAlignment(boot_samples + (boot*nptn), params.bootstrap_spec, rstream);
+ 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,7 +1016,37 @@ 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";
+ }
+ }
}
void ModelMarkov::readRates(string str) throw(const char*) {
=====================================
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;
- *x |= ~((1<<(site%UINT_BITS)) - 1);
- x++;
+ 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));
}
=====================================
tree/phylotree.cpp
=====================================
@@ -99,7 +99,7 @@ void PhyloTree::init() {
dist_matrix = NULL;
var_matrix = NULL;
params = NULL;
- setLikelihoodKernel(LK_SSE2); // FOR TUNG: you forgot to initialize this variable!
+ setLikelihoodKernel(LK_SSE2); // FOR TUNG: you forgot to initialize this variable!
setNumThreads(1);
num_threads = 0;
max_lh_slots = 0;
@@ -383,41 +383,41 @@ void PhyloTree::setAlignment(Alignment *alignment) {
}
void PhyloTree::setRootNode(const char *my_root, bool multi_taxa) {
- if (rooted) {
- computeBranchDirection();
- return;
- }
-
- if (!my_root) {
- root = findNodeName(aln->getSeqName(0));
- ASSERT(root);
- return;
- }
-
- if (strchr(my_root, ',') == NULL) {
- string root_name = my_root;
- root = findNodeName(root_name);
- ASSERT(root);
- return;
- }
-
- // my_root is a list of taxa
- StrVector taxa;
- convert_string_vec(my_root, taxa);
- root = findNodeName(taxa[0]);
- ASSERT(root);
- if (!multi_taxa) {
- return;
- }
- unordered_set<string> taxa_set;
- for (auto it = taxa.begin(); it != taxa.end(); it++)
- taxa_set.insert(*it);
- pair<Node*,Neighbor*> res = {NULL, NULL};
- findNodeNames(taxa_set, res, root->neighbors[0]->node, root);
- if (res.first)
- root = res.first;
- else
- outWarning("Branch separating outgroup is not found");
+ if (rooted) {
+ computeBranchDirection();
+ return;
+ }
+
+ if (!my_root) {
+ root = findNodeName(aln->getSeqName(0));
+ ASSERT(root);
+ return;
+ }
+
+ if (strchr(my_root, ',') == NULL) {
+ string root_name = my_root;
+ root = findNodeName(root_name);
+ ASSERT(root);
+ return;
+ }
+
+ // my_root is a list of taxa
+ StrVector taxa;
+ convert_string_vec(my_root, taxa);
+ root = findNodeName(taxa[0]);
+ ASSERT(root);
+ if (!multi_taxa) {
+ return;
+ }
+ unordered_set<string> taxa_set;
+ for (auto it = taxa.begin(); it != taxa.end(); it++)
+ taxa_set.insert(*it);
+ pair<Node*,Neighbor*> res = {NULL, NULL};
+ findNodeNames(taxa_set, res, root->neighbors[0]->node, root);
+ if (res.first)
+ root = res.first;
+ else
+ outWarning("Branch separating outgroup is not found");
}
//void PhyloTree::setParams(Params* params) {
@@ -548,10 +548,10 @@ void PhyloTree::setModelFactory(ModelFactory *model_fac) {
model_factory = model_fac;
if (model_fac) {
model = model_factory->model;
- site_rate = model_factory->site_rate;
- if (!isSuperTree()) {
- setLikelihoodKernel(sse);
- setNumThreads(num_threads);
+ site_rate = model_factory->site_rate;
+ if (!isSuperTree()) {
+ setLikelihoodKernel(sse);
+ setNumThreads(num_threads);
}
} else {
model = NULL;
@@ -2352,24 +2352,24 @@ const double TOL_TREE_LENGTH_SCALE = 0.001;
double PhyloTree::optimizeTreeLengthScaling(double min_scaling, double &scaling, double max_scaling, double gradient_epsilon) {
is_opt_scaling = true;
current_scaling = scaling;
- double negative_lh, ferror;
- // 2018-08-20: make sure that max and min branch lengths do not go over bounds
- vector<DoubleVector> brlens;
- brlens.resize(branchNum);
- getBranchLengths(brlens);
- double min_brlen = params->max_branch_length, max_brlen = params->min_branch_length;
- for (auto brlenvec = brlens.begin(); brlenvec != brlens.end(); brlenvec++) {
- for (auto brlen = brlenvec->begin(); brlen != brlenvec->end(); brlen++) {
- max_brlen = max(max_brlen, *brlen);
- min_brlen = min(min_brlen, *brlen);
- }
- }
- if (min_brlen <= 0.0) min_brlen = params->min_branch_length;
- if (max_scaling > params->max_branch_length / max_brlen)
- max_scaling = params->max_branch_length / max_brlen;
- if (min_scaling < params->min_branch_length / min_brlen)
- min_scaling = params->min_branch_length / min_brlen;
-
+ double negative_lh, ferror;
+ // 2018-08-20: make sure that max and min branch lengths do not go over bounds
+ vector<DoubleVector> brlens;
+ brlens.resize(branchNum);
+ getBranchLengths(brlens);
+ double min_brlen = params->max_branch_length, max_brlen = params->min_branch_length;
+ for (auto brlenvec = brlens.begin(); brlenvec != brlens.end(); brlenvec++) {
+ for (auto brlen = brlenvec->begin(); brlen != brlenvec->end(); brlen++) {
+ max_brlen = max(max_brlen, *brlen);
+ min_brlen = min(min_brlen, *brlen);
+ }
+ }
+ if (min_brlen <= 0.0) min_brlen = params->min_branch_length;
+ if (max_scaling > params->max_branch_length / max_brlen)
+ max_scaling = params->max_branch_length / max_brlen;
+ if (min_scaling < params->min_branch_length / min_brlen)
+ min_scaling = params->min_branch_length / min_brlen;
+
scaling = minimizeOneDimen(min(scaling, min_scaling), scaling, max(max_scaling, scaling), max(TOL_TREE_LENGTH_SCALE, gradient_epsilon), &negative_lh, &ferror);
if (scaling != current_scaling) {
scaleLength(scaling / current_scaling);
@@ -4340,7 +4340,7 @@ void PhyloTree::resampleLh(double **pat_lh, double *lh_new, int *rstream) {
lh_new[0] += boot_freq[i] * pat_lh[0][i];
lh_new[1] += boot_freq[i] * pat_lh[1][i];
lh_new[2] += boot_freq[i] * pat_lh[2][i];
- }
+ }
aligned_free(boot_freq);
}
@@ -4523,14 +4523,14 @@ double PhyloTree::testOneBranch(double best_score, double *pattern_lh, int reps,
double lh[NUM_NNI];
double *pat_lh[NUM_NNI];
lh[0] = best_score;
- pat_lh[0] = pattern_lh;
+ pat_lh[0] = pattern_lh;
int nptn = getAlnNPattern();
pat_lh[1] = new double[nptn];
pat_lh[2] = new double[nptn];
- int tmp = save_all_trees;
- save_all_trees = 0;
+ int tmp = save_all_trees;
+ save_all_trees = 0;
computeNNIPatternLh(best_score, lh[1], pat_lh[1], lh[2], pat_lh[2], node1, node2);
- save_all_trees = tmp;
+ save_all_trees = tmp;
double aLRT;
if (lh[1] > lh[2])
aLRT = (lh[0] - lh[1]);
@@ -4546,7 +4546,7 @@ double PhyloTree::testOneBranch(double best_score, double *pattern_lh, int reps,
aBayes_support = 1.0 / (1.0 + exp(lh[1]-lh[0]) + exp(lh[2]-lh[0]));
- int SH_aLRT_support = 0;
+ int SH_aLRT_support = 0;
int lbp_support_int = 0;
int times = max(reps, lbp_reps);
@@ -4554,22 +4554,22 @@ double PhyloTree::testOneBranch(double best_score, double *pattern_lh, int reps,
if (max(lh[1],lh[2]) == -DBL_MAX) {
SH_aLRT_support = times;
outWarning("Branch where both NNIs violate constraint tree will show 100% SH-aLRT support");
- } else
-#ifdef _OPENMP
-#pragma omp parallel
- {
- int *rstream;
- init_random(params->ran_seed + omp_get_thread_num(), false, &rstream);
-#pragma omp for reduction(+: lbp_support_int, SH_aLRT_support)
-#endif
+ } else
+#ifdef _OPENMP
+#pragma omp parallel
+ {
+ int *rstream;
+ init_random(params->ran_seed + omp_get_thread_num(), false, &rstream);
+#pragma omp for reduction(+: lbp_support_int, SH_aLRT_support)
+#endif
for (int i = 0; i < times; i++) {
double lh_new[NUM_NNI];
// resampling estimated log-likelihood (RELL)
-#ifdef _OPENMP
+#ifdef _OPENMP
resampleLh(pat_lh, lh_new, rstream);
-#else
- resampleLh(pat_lh, lh_new, randstream);
-#endif
+#else
+ resampleLh(pat_lh, lh_new, randstream);
+#endif
if (lh_new[0] > lh_new[1] && lh_new[0] > lh_new[2])
lbp_support_int++;
double cs[NUM_NNI], cs_best, cs_2nd_best;
@@ -4597,14 +4597,14 @@ double PhyloTree::testOneBranch(double best_score, double *pattern_lh, int reps,
}
if (aLRT > (cs_best - cs_2nd_best) + 0.05)
SH_aLRT_support++;
- }
-#ifdef _OPENMP
- finish_random(rstream);
- }
-#endif
+ }
+#ifdef _OPENMP
+ finish_random(rstream);
+ }
+#endif
delete[] pat_lh[2];
delete[] pat_lh[1];
-
+
lbp_support = 0.0;
if (times > 0)
lbp_support = ((double)lbp_support_int) / times;
@@ -5078,7 +5078,7 @@ void PhyloTree::convertToRooted() {
}
void PhyloTree::convertToUnrooted() {
- ASSERT(rooted);
+ ASSERT(rooted);
if (aln)
ASSERT(leafNum == aln->getNSeq()+1);
ASSERT(root);
@@ -5111,9 +5111,9 @@ void PhyloTree::convertToUnrooted() {
delete root;
// set a temporary taxon so that tree traversal works
- root = taxon;
-
- if (params)
+ root = taxon;
+
+ if (params)
setRootNode(params->root);
initializeTree();
@@ -5324,8 +5324,8 @@ void PhyloTree::writeSiteRates(ostream &out, int partid) {
cout << " " << ((double)count[i])/nsite;
cout << endl;
}
-
-void PhyloTree::writeBranches(ostream &out) {
- outError("Please only use this feature with partition model");
-}
-
+
+void PhyloTree::writeBranches(ostream &out) {
+ outError("Please only use this feature with partition model");
+}
+
=====================================
utils/tools.cpp
=====================================
@@ -3505,8 +3505,8 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
#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)");
View it on GitLab: https://salsa.debian.org/med-team/iqtree/compare/def15b03e5fca0a7baeb9235e4ba0cd3ee82cc25...31fe0ddb8a4814ee240ff733975e6c8010ee8bbc
--
View it on GitLab: https://salsa.debian.org/med-team/iqtree/compare/def15b03e5fca0a7baeb9235e4ba0cd3ee82cc25...31fe0ddb8a4814ee240ff733975e6c8010ee8bbc
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20181213/a208fece/attachment-0001.html>
More information about the debian-med-commit
mailing list