[med-svn] [Git][med-team/iqtree][upstream] New upstream version 1.6.11+dfsg
Andreas Tille
gitlab at salsa.debian.org
Sat Jul 27 22:33:33 BST 2019
Andreas Tille pushed to branch upstream at Debian Med / iqtree
Commits:
93437d2c by Andreas Tille at 2019-07-27T20:09:24Z
New upstream version 1.6.11+dfsg
- - - - -
18 changed files:
- CMakeLists.txt
- alignment/alignment.cpp
- alignment/superalignment.cpp
- main/main.cpp
- model/modeldna.cpp
- model/modelfactory.cpp
- model/modelprotein.cpp
- tree/iqtree.cpp
- tree/iqtree.h
- tree/phylokernelnew.h
- tree/phylokernelnonrev.h
- tree/phylosupertree.cpp
- tree/phylosupertree.h
- tree/phylotree.h
- tree/phylotreepars.cpp
- utils/MPIHelper.cpp
- utils/MPIHelper.h
- 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 "9")
+set (iqtree_VERSION_PATCH "11")
set(BUILD_SHARED_LIBS OFF)
@@ -277,9 +277,14 @@ if (NOT IQTREE_FLAGS MATCHES "single")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -pthread")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp -pthread")
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")
+ if (APPLE)
+ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -pthread")
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xpreprocessor -fopenmp -pthread")
+ set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Xpreprocessor -fopenmp -lomp")
+ else()
+ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -pthread")
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp -pthread")
+ endif()
endif()
else()
message("OpenMP : NONE")
@@ -509,12 +514,10 @@ add_library(kernelfma tree/phylokernelfma.cpp)
endif()
endif()
-if (IQTREE_FLAGS MATCHES "mpi")
- add_library(mympi utils/TreeCollection.cpp utils/ObjectStream.cpp)
-endif()
-
if (NOT IQTREE_FLAGS MATCHES "single")
- if (BINARY32)
+ if (APPLE)
+ link_directories(${PROJECT_SOURCE_DIR}/libmac)
+ elseif (BINARY32)
link_directories(${PROJECT_SOURCE_DIR}/lib32)
else()
link_directories(${PROJECT_SOURCE_DIR}/lib)
@@ -602,7 +605,6 @@ endif()
# MPI libraries
if (IQTREE_FLAGS MATCHES "mpi")
- target_link_libraries(iqtree mympi)
if (NOT CMAKE_CXX_COMPILER MATCHES "mpi")
target_link_libraries(iqtree ${MPI_CXX_LIBRARIES})
endif()
=====================================
alignment/alignment.cpp
=====================================
@@ -2943,13 +2943,22 @@ void convert_range(const char *str, int &lower, int &upper, int &step_size, char
if (*endptr != '-') return;
// parse the upper bound of the range
- str = endptr+1;
+ endptr++;
+ // skip blank chars
+ for (; *endptr == ' '; endptr++) {}
+ str = endptr;
d = strtol(str, &endptr, 10);
if ((d == 0 && endptr == str) || abs(d) == HUGE_VALL) {
- string err = "Expecting integer, but found \"";
- err += str;
- err += "\" instead";
- throw err;
+ if (str[0] == '.') {
+ // 2019-06-03: special character '.' for whatever ending position
+ d = lower-1;
+ endptr++;
+ } else {
+ string err = "Expecting integer, but found \"";
+ err += str;
+ err += "\" instead";
+ throw err;
+ }
}
//lower = d_save;
@@ -2982,6 +2991,9 @@ void extractSiteID(Alignment *aln, const char* spec, IntVector &site_id) {
for (; *str != 0; ) {
int lower, upper, step;
convert_range(str, lower, upper, step, str);
+ // 2019-06-03: special '.' character
+ if (upper == lower-1)
+ upper = aln->getNSite();
lower--;
upper--;
nchars += (upper-lower+1)/step;
=====================================
alignment/superalignment.cpp
=====================================
@@ -935,6 +935,7 @@ void SuperAlignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern
auto site_pattern = boot_aln->site_pattern;
for (int j = 1; j < gene_freq[i]; j++)
boot_aln->site_pattern.insert(boot_aln->site_pattern.end(), site_pattern.begin(), site_pattern.end());
+ boot_aln->countConstSite();
}
partitions.push_back(boot_aln);
}
=====================================
main/main.cpp
=====================================
@@ -68,9 +68,6 @@
#include "vectorclass/instrset.h"
#include "utils/MPIHelper.h"
-#ifdef _IQTREE_MPI
-#include <mpi.h>
-#endif
#ifdef _OPENMP
#include <omp.h>
@@ -1721,6 +1718,7 @@ void processNCBITree(Params ¶ms) {
class outstreambuf : public streambuf {
public:
outstreambuf* open( const char* name, ios::openmode mode = ios::out);
+ bool is_open();
outstreambuf* close();
~outstreambuf() { close(); }
streambuf *get_fout_buf() {
@@ -1753,6 +1751,10 @@ outstreambuf* outstreambuf::open( const char* name, ios::openmode mode) {
return this;
}
+bool outstreambuf::is_open() {
+ return fout.is_open();
+}
+
outstreambuf* outstreambuf::close() {
cout.rdbuf(cout_buf);
if ( fout.is_open()) {
@@ -1850,8 +1852,8 @@ extern "C" void startLogFile(bool append_log) {
}
extern "C" void endLogFile() {
- _out_buf.close();
-
+ if (_out_buf.is_open())
+ _out_buf.close();
}
void funcExit(void) {
@@ -1862,6 +1864,7 @@ void funcExit(void) {
}
endLogFile();
+ MPIHelper::getInstance().finalize();
}
extern "C" void funcAbort(int signal_number)
@@ -2218,20 +2221,9 @@ int main(int argc, char *argv[]) {
*/
int instruction_set;
-#ifdef _IQTREE_MPI
- double time_initial, time_current;
- int n_tasks, task_id;
- if (MPI_Init(&argc, &argv) != MPI_SUCCESS) {
- outError("MPI initialization failed!");
- }
- MPI_Comm_size(MPI_COMM_WORLD, &n_tasks);
- MPI_Comm_rank(MPI_COMM_WORLD, &task_id);
- MPIHelper::getInstance().setNumProcesses(n_tasks);
- MPIHelper::getInstance().setProcessID(task_id);
- MPIHelper::getInstance().setNumTreeReceived(0);
- MPIHelper::getInstance().setNumTreeSent(0);
- MPIHelper::getInstance().setNumNNISearch(0);
-#endif
+ MPIHelper::getInstance().init(argc, argv);
+
+ atexit(funcExit);
/*************************/
{ /* local scope */
@@ -2289,18 +2281,23 @@ int main(int argc, char *argv[]) {
if (checkpoint->hasKey("finished")) {
if (checkpoint->getBool("finished")) {
if (Params::getInstance().force_unfinished) {
- cout << "NOTE: Continue analysis although a previous run already finished" << endl;
+ if (MPIHelper::getInstance().isMaster())
+ cout << "NOTE: Continue analysis although a previous run already finished" << endl;
} else {
- outError("Checkpoint (" + filename + ") indicates that a previous run successfully finished\n" +
- "Use `-redo` option if you really want to redo the analysis and overwrite all output files.");
delete checkpoint;
- return EXIT_FAILURE;
+ if (MPIHelper::getInstance().isMaster())
+ outError("Checkpoint (" + filename + ") indicates that a previous run successfully finished\n" +
+ "Use `-redo` option if you really want to redo the analysis and overwrite all output files.");
+ else
+ exit(EXIT_SUCCESS);
+ exit(EXIT_FAILURE);
}
} else {
append_log = true;
}
} else {
- outWarning("Ignore invalid checkpoint file " + filename);
+ if (MPIHelper::getInstance().isMaster())
+ outWarning("Ignore invalid checkpoint file " + filename);
checkpoint->clear();
}
}
@@ -2318,26 +2315,10 @@ int main(int argc, char *argv[]) {
cout << endl << "******************************************************"
<< endl << "CHECKPOINT: Resuming analysis from " << filename << endl << endl;
}
-#ifdef _IQTREE_MPI
- cout << "************************************************" << endl;
- cout << "* START TREE SEARCH USING MPI WITH " << MPIHelper::getInstance().getNumProcesses() << " PROCESSES *" << endl;
- cout << "************************************************" << endl;
- unsigned int rndSeed;
- if (MPIHelper::getInstance().isMaster()) {
- rndSeed = Params::getInstance().ran_seed;
- cout << "Random seed of master = " << rndSeed << endl;
- }
- // Broadcast random seed
- MPI_Bcast(&rndSeed, 1, MPI_INT, PROC_MASTER, MPI_COMM_WORLD);
- if (MPIHelper::getInstance().isWorker()) {
-// Params::getInstance().ran_seed = rndSeed + task_id * 100000;
- Params::getInstance().ran_seed = rndSeed;
-// printf("Process %d: random_seed = %d\n", task_id, Params::getInstance().ran_seed);
- }
-#endif
- atexit(funcExit);
- signal(SIGABRT, &funcAbort);
+ MPIHelper::getInstance().syncRandomSeed();
+
+ signal(SIGABRT, &funcAbort);
signal(SIGFPE, &funcAbort);
signal(SIGILL, &funcAbort);
signal(SIGSEGV, &funcAbort);
@@ -2469,8 +2450,10 @@ int main(int argc, char *argv[]) {
}
#endif
-
-#ifndef _IQTREE_MPI
+#ifdef _IQTREE_MPI
+ cout << endl << "MPI: " << MPIHelper::getInstance().getNumProcesses() << " processes";
+#endif
+
int num_procs = countPhysicalCPUCores();
#ifdef _OPENMP
if (num_procs > 1 && Params::getInstance().num_threads == 1) {
@@ -2480,7 +2463,6 @@ int main(int argc, char *argv[]) {
#else
if (num_procs > 1)
cout << endl << endl << "NOTE: Consider using the multicore version because your CPU has " << num_procs << " cores!";
-#endif
#endif
//cout << "sizeof(int)=" << sizeof(int) << endl;
@@ -2638,8 +2620,5 @@ int main(int argc, char *argv[]) {
finish_random();
-#ifdef _IQTREE_MPI
- MPI_Finalize();
-#endif
return EXIT_SUCCESS;
}
=====================================
model/modeldna.cpp
=====================================
@@ -377,8 +377,14 @@ bool ModelDNA::setRateType(string rate_str) {
cout << rates[i] << " ";
cout << endl;
}
- param_fixed.resize(num_params+1, false);
- param_fixed[0] = true; // fix the last entry
+ if (param_fixed.size() == num_params + 1) {
+ num_params = 0;
+ for (auto p : param_fixed)
+ if (!p) num_params++;
+ } else {
+ param_fixed.resize(num_params+1, false);
+ param_fixed[0] = true; // fix the last entry
+ }
delete [] num_rates;
delete [] avg_rates;
return true;
=====================================
model/modelfactory.cpp
=====================================
@@ -579,7 +579,8 @@ ModelFactory::ModelFactory(Params ¶ms, string &model_name, PhyloTree *tree,
outError("Invalid use of +ASC because of " + convertIntToString(tree->aln->frac_invariant_sites*tree->aln->getNSite()) +
" invariant sites in the alignment");
}
- cout << "Ascertainment bias correction: " << unobserved_ptns.size() << " unobservable constant patterns"<< endl;
+ if (verbose_mode >= VB_MED)
+ cout << "Ascertainment bias correction: " << unobserved_ptns.size() << " unobservable constant patterns"<< endl;
rate_str = rate_str.substr(0, posasc) + rate_str.substr(posasc+4);
} else {
tree->aln->buildSeqStates(false);
=====================================
model/modelprotein.cpp
=====================================
@@ -75,7 +75,9 @@ model DAYHOFF=
0 201 23 0 0 0 0 0 27 0 46 0 0 76 0 75 0
24 8 95 0 96 0 22 0 127 37 28 13 0 698 0 34 42 61
208 24 15 18 49 35 37 54 44 889 175 10 258 12 48 30 157 0 28
-0.087127 0.040904 0.040432 0.046872 0.033474 0.038255 0.049530 0.088612 0.033618 0.036886 0.085357 0.080482 0.014753 0.039772 0.050680 0.069577 0.058542 0.010494 0.029916 0.064717;
+0.08712691 0.04090396 0.04043196 0.04687195 0.03347397 0.03825496 0.04952995 0.08861191 0.03361897 0.03688596 0.08535691 0.08048092 0.01475299 0.03977196 0.05067995 0.06957693 0.05854194 0.01049399 0.02991597 0.06471794;
+[ NOTE 2019-06-04: normalised from original Dayhoff freqs, which do not sum to 1.0
+ https://www.ebi.ac.uk/goldman-srv/dayhoff/dayhoff-dcmut.dat ]
model DCMUT=
26.7828
@@ -97,7 +99,9 @@ model DCMUT=
0.0000 200.1375 22.4968 0.0000 0.0000 0.0000 0.0000 0.0000 27.0564 0.0000 46.1776 0.0000 0.0000 76.2354 0.0000 74.0819 0.0000
24.4139 7.8012 94.6940 0.0000 95.3164 0.0000 21.4717 0.0000 126.5400 37.4834 28.6572 13.2142 0.0000 695.2629 0.0000 33.6289 41.7839 60.8070
205.9564 24.0368 15.8067 17.8316 48.4678 34.6983 36.7250 53.8165 43.8715 881.0038 174.5156 10.3850 256.5955 12.3606 48.5026 30.3836 156.1997 0.0000 27.9379
-0.087127 0.040904 0.040432 0.046872 0.033474 0.038255 0.049530 0.088612 0.033619 0.036886 0.085357 0.080481 0.014753 0.039772 0.050680 0.069577 0.058542 0.010494 0.029916 0.064717;
+ 0.08712691 0.04090396 0.04043196 0.04687195 0.03347397 0.03825496 0.04952995 0.08861191 0.03361897 0.03688596 0.08535691 0.08048092 0.01475299 0.03977196 0.05067995 0.06957693 0.05854194 0.01049399 0.02991597 0.06471794;
+[ NOTE 2019-06-04: normalised from original Dayhoff-DCMUT freqs, which do not sum to 1.0
+ https://www.ebi.ac.uk/goldman-srv/dayhoff/dayhoff-dcmut.dat ]
model JTT=
58
@@ -119,7 +123,8 @@ model JTT=
9 126 8 4 115 18 10 55 8 9 52 10 24 53 6 35 12
11 20 70 46 209 24 7 8 573 32 24 8 18 536 10 63 21 71
298 17 16 31 62 20 45 47 11 961 180 14 323 62 23 38 112 25 16
-0.076748 0.051691 0.042645 0.051544 0.019803 0.040752 0.061830 0.073152 0.022944 0.053761 0.091904 0.058676 0.023826 0.040126 0.050901 0.068765 0.058565 0.014261 0.032102 0.066004;
+0.07674792 0.05169095 0.04264496 0.05154395 0.01980298 0.04075196 0.06182994 0.07315193 0.02294398 0.05376095 0.09190391 0.05867594 0.02382598 0.04012596 0.05090095 0.06876493 0.05856494 0.01426099 0.03210197 0.06600493;
+[ NOTE 2019-06-04: original JTT freqs do not sum to 1.0, taken from PAML package ]
model MTREV=
23.18
@@ -142,7 +147,9 @@ model MTREV=
6.48 1.90 191.36 21.21 254.77 38.82 13.12 3.21 670.14 25.01 44.15 51.17 39.96 465.58 16.21 64.92 38.73 26.25
195.06 7.64 1.90 1.90 1.90 19.00 21.14 2.53 1.90 1222.94 91.67 1.90 387.54 6.35 8.23 1.90 204.54 5.37 1.90
0.072 0.019 0.039 0.019 0.006 0.025 0.024 0.056 0.028 0.088 0.169 0.023 0.054 0.061 0.054 0.072 0.086 0.029 0.033 0.043;
-
+[ NOTE 2019-06-04: The PI's used to sum to 0.999 and I (Z. Yang) changed one of the freq from 0.168
+ into 0.169 so that the sum is 1. Prepared by Z. Yang ]
+
model WAG=
55.15710
50.98480 63.53460
@@ -163,8 +170,8 @@ model WAG=
11.31330 116.39200 7.19167 12.97670 71.70700 21.57370 15.65570 33.69830 26.25690 21.24830 66.53090 13.75050 51.57060 152.96400 13.94050 52.37420 11.08640
24.07350 38.15330 108.60000 32.57110 54.38330 22.77100 19.63030 10.36040 387.34400 42.01700 39.86180 13.32640 42.84370 645.42800 21.60460 78.69930 29.11480 248.53900
200.60100 25.18490 19.62460 15.23350 100.21400 30.12810 58.87310 18.72470 11.83580 782.13000 180.03400 30.54340 205.84500 64.98920 31.48870 23.27390 138.82300 36.53690 31.47300
-
-0.0866279 0.043972 0.0390894 0.0570451 0.0193078 0.0367281 0.0580589 0.0832518 0.0244313 0.048466 0.086209 0.0620286 0.0195027 0.0384319 0.0457631 0.0695179 0.0610127 0.0143859 0.0352742 0.0708957;
+0.08662791 0.043972 0.0390894 0.05704511 0.0193078 0.0367281 0.05805891 0.08325181 0.0244313 0.048466 0.08620901 0.06202861 0.0195027 0.0384319 0.0457631 0.06951791 0.06101271 0.0143859 0.0352742 0.07089561;
+[ NOTE 2019-06-04: normalised from original WAG freqs, which do not sum to 1.0 ]
model RTREV=
34
@@ -208,7 +215,8 @@ model CPREV=
14 230 40 18 435 53 63 82 69 42 159 10 86 468 49 73 29
56 323 754 281 1466 391 142 10 1971 89 189 247 215 2370 97 522 71 346
968 92 83 75 592 54 200 91 25 4797 865 249 475 317 122 167 760 10 119
-0.0760 0.0620 0.0410 0.0370 0.0090 0.0380 0.0490 0.0840 0.0250 0.0810 0.1010 0.0500 0.0220 0.0510 0.0430 0.0620 0.0540 0.0180 0.0310 0.0660;
+0.0755 0.0621 0.0410 0.0371 0.0091 0.0382 0.0495 0.0838 0.0246 0.0806 0.1011 0.0504 0.0220 0.0506 0.0431 0.0622 0.0543 0.0181 0.0307 0.0660;
+[ NOTE 2019-06-04: CPREV freqs taken from PAML package with higher precision ]
model VT=
1.2412691067876198
@@ -296,7 +304,9 @@ model LG=
0.180717 0.593607 0.045376 0.029890 0.670128 0.236199 0.077852 0.268491 0.597054 0.111660 0.619632 0.049906 0.696175 2.457121 0.095131 0.248862 0.140825
0.218959 0.314440 0.612025 0.135107 1.165532 0.257336 0.120037 0.054679 5.306834 0.232523 0.299648 0.131932 0.481306 7.803902 0.089613 0.400547 0.245841 3.151815
2.547870 0.170887 0.083688 0.037967 1.959291 0.210332 0.245034 0.076701 0.119013 10.649107 1.702745 0.185202 1.898718 0.654683 0.296501 0.098369 2.188158 0.189510 0.249313
-0.079066 0.055941 0.041977 0.053052 0.012937 0.040767 0.071586 0.057337 0.022355 0.062157 0.099081 0.064600 0.022951 0.042302 0.044040 0.061197 0.053287 0.012066 0.034155 0.069146;
+0.07906592 0.05594094 0.04197696 0.05305195 0.01293699 0.04076696 0.07158593 0.05733694 0.02235498 0.06215694 0.0990809 0.06459994 0.02295098 0.04230196 0.04403996 0.06119694 0.05328695 0.01206599 0.03415497 0.06914693;
+[ NOTE 2019-06-04: normalised from original LG freqs, which do not sum to 1.0
+ http://www.atgc-montpellier.fr/download/datasets/models/lg_LG.PAML.txt ]
model MTART=
0.2
@@ -340,7 +350,8 @@ model MTZOA=
3.1 16.9 6.4 0.2 36.1 6.1 3.5 12.3 4.5 9.7 27.2 6.6 48.7 58.2 1.3 10.3 3.6
2.1 13.8 141.6 13.9 76.7 52.3 10.0 4.3 266.5 13.1 5.7 45.0 41.4 590.5 4.2 29.7 29.0 79.8
321.9 5.1 7.1 3.7 243.8 9.0 16.3 23.7 0.3 1710.6 126.1 11.1 279.6 59.6 17.9 49.5 396.4 13.7 15.6
-0.069 0.021 0.030 0.020 0.010 0.019 0.025 0.072 0.027 0.085 0.157 0.019 0.051 0.082 0.045 0.081 0.056 0.028 0.037 0.066;
+0.06887993 0.02103698 0.03038997 0.02069598 0.00996599 0.01862298 0.02498898 0.07196793 0.02681397 0.08507191 0.15671684 0.01927598 0.05065195 0.08171192 0.04480296 0.08053492 0.05638594 0.02799797 0.03740396 0.06608293;
+[ NOTE 2019-06-04: original mtzoa freqs do not sum to 1.0, modified from PAML package ]
model PMB=
0.674995699
@@ -362,7 +373,10 @@ model PMB=
0.350110510 0.618778365 0.422407388 0.362495245 0.445669347 0.720384740 0.261258229 0.378748270 0.724367510 0.516260502 0.794797115 0.433409620 0.768395107 3.295193440 0.499869138 0.496334956 0.383723610
0.573154753 0.628599063 0.720013799 0.436220437 0.556261630 0.728970584 0.507200030 0.284727562 2.210952064 0.570562395 0.811019594 0.664884513 0.932536060 5.894735673 0.433748126 0.593795813 0.523549536 2.996248013
2.063050067 0.388680158 0.474418852 0.275658381 0.998911631 0.634408285 0.527640634 0.314700907 0.305792277 8.002789424 2.113077156 0.526184203 1.737356217 0.983844803 0.551333603 0.507506011 1.899650790 0.429570747 0.716795463
-0.076 0.054 0.038 0.045 0.028 0.034 0.053 0.078 0.030 0.060 0.096 0.052 0.022 0.045 0.042 0.068 0.056 0.016 0.036 0.071;
+
+ 0.07559244 0.05379462 0.03769623 0.04469553 0.02849715 0.03389661 0.05349465 0.0779922 0.029997 0.05989401 0.09579042 0.0519948 0.02189781 0.0449955 0.0419958 0.06819318 0.05639436 0.01569843 0.0359964 0.07149285;
+ [ NOTE 2019-06-04: normalised from original PMB freqs, which do not sum to 1.0:
+ 0.0756 0.0538 0.0377 0.0447 0.0285 0.0339 0.0535 0.0780 0.0300 0.0599 0.0958 0.0520 0.0219 0.0450 0.0420 0.0682 0.0564 0.0157 0.0360 0.0715 ]
model HIVB=
0.30750700
@@ -428,7 +442,9 @@ model JTTDCMUT=
0.084329 1.257961 0.027700 0.057466 1.104181 0.172206 0.114381 0.544180 0.128193 0.134510 0.530324 0.089134 0.201334 0.537922 0.069965 0.310927 0.080556
0.139492 0.235601 0.700693 0.453952 2.114852 0.254745 0.063452 0.052500 5.848400 0.303445 0.241094 0.087904 0.189870 5.484236 0.113850 0.628608 0.201094 0.747889
2.924161 0.171995 0.164525 0.315261 0.621323 0.179771 0.465271 0.470140 0.121827 9.533943 1.761439 0.124066 3.038533 0.593478 0.211561 0.408532 1.143980 0.239697 0.165473
-0.077 0.051 0.043 0.051 0.020 0.041 0.062 0.075 0.023 0.053 0.091 0.059 0.024 0.040 0.051 0.068 0.059 0.014 0.032 0.066;
+0.07686192 0.05105695 0.04254596 0.05126895 0.02027898 0.04106096 0.06181994 0.07471393 0.02298298 0.05256895 0.09111091 0.05949794 0.02341398 0.04052996 0.05053195 0.06822493 0.05851794 0.01433599 0.03230297 0.06637393;
+[ NOTE 2019-06-04: normalised from original JTTDCMUTT freqs, which do not sum to 1.0
+ https://www.ebi.ac.uk/goldman-srv/dayhoff/jtt-dcmut.dat ]
model FLU=
0.138658765
@@ -450,8 +466,9 @@ model FLU=
0.196000000 1.369429408 0.000536000 0.000014900 0.094106680 0.044000000 0.155245492 0.196486447 0.022400000 0.032132150 0.431277663 0.000049800 0.070460039 0.814753094 0.000431021 0.099835753 0.207066206
0.018289288 0.099855497 0.373101927 0.525398543 0.601692431 0.072205935 0.104092870 0.074814997 6.448954446 0.273934263 0.340058468 0.012416222 0.874272175 5.393924245 0.000182000 0.392552240 0.124898020 0.427755430
3.532005270 0.103964386 0.010257517 0.297123975 0.054904564 0.406697814 0.285047948 0.337229619 0.098631355 14.394052190 0.890598579 0.073127930 4.904842235 0.592587985 0.058971975 0.088256423 0.654109108 0.256900461 0.167581647
-0.0471 0.0509 0.0742 0.0479 0.0250 0.0333 0.0546 0.0764 0.0200 0.0671 0.0715 0.0568 0.0181 0.0305 0.0507 0.0884 0.0743 0.0185 0.0315 0.0632;
-
+0.04707195 0.05090995 0.07421393 0.04785995 0.02502197 0.03330397 0.05458695 0.07637292 0.01996398 0.06713393 0.07149793 0.05678494 0.01815098 0.03049597 0.05065595 0.08840891 0.07433893 0.01852398 0.03147397 0.06322894;
+[ NOTE 2019-06-04: normalised from FLU freqs in PhyML, which do not sum 1.0 ]
+
model MTMET=
0.058078177576542
0.032893910131824 0.141364232590718
@@ -567,7 +584,7 @@ void ModelProtein::init(const char *model_name, string model_params, StateFreqTy
double sum = 0.0;
for (i = 0; i < num_states; i++)
sum += (double) state_freq[i];
- if (round(sum*1e8) != 1e8) {
+ if (fabs(sum-1.0) > 1e-7) {
cout.precision(7);
cout << "WARNING: " << name_upper << " state frequencies do not sum up to 1: " << sum << endl;
}
=====================================
tree/iqtree.cpp
=====================================
@@ -35,10 +35,6 @@
#include "utils/MPIHelper.h"
#include "utils/pllnni.h"
-#ifdef _IQTREE_MPI
-#include <mpi.h>
-#endif
-
Params *globalParams;
Alignment *globalAlignment;
extern StringIntMap pllTreeCounter;
@@ -830,25 +826,6 @@ void IQTree::initCandidateTreeSet(int nParTrees, int nNNITrees) {
cout << "Current best score: " << candidateTrees.getBestScore() << endl;
}
-/*
- //---- NON-BLOCKING COMMUNICATION
-#ifdef _IQTREE_MPI
- vector<string> trees;
- vector<double> scores;
- // FIX BUG: send candidateTrees instead of intermediateTrees
- candidateTrees.getAllTrees(trees, scores, WT_TAXON_ID + WT_BR_LEN + WT_BR_LEN_SHORT);
- // Send all trees to other processes
- MPIHelper::getInstance().distributeTrees(trees, scores, TREE_TAG);
-
- // Get trees from other nodes
- cout << "Getting initial trees from other processes ... " << endl;
- int maxNumTrees = (nParTrees + 2) * (MPIHelper::getInstance().getNumProcesses() - 1);
- MPI_CollectTrees(true, maxNumTrees, false);
-
- MPI_Barrier(MPI_COMM_WORLD);
-#endif
-*/
-
//---- BLOCKING COMMUNICATION
syncCandidateTrees(nNNITrees, false);
@@ -875,10 +852,6 @@ void IQTree::initCandidateTreeSet(int nParTrees, int nNNITrees) {
addTreeToCandidateSet(treeString, curScore, true, MPIHelper::getInstance().getProcessID());
if (Params::getInstance().writeDistImdTrees)
intermediateTrees.update(treeString, curScore);
-//#ifdef _IQTREE_MPI
-// MPIHelper::getInstance().distributeTree(getTreeString(), curScore, TREE_TAG);
-// MPI_CollectTrees(false, maxNumTrees, true);
-//#endif
}
// TODO turning this
@@ -902,32 +875,6 @@ void IQTree::initCandidateTreeSet(int nParTrees, int nNNITrees) {
//---- BLOCKING COMMUNICATION
syncCandidateTrees(Params::getInstance().numSupportTrees, true);
-/*
-#ifdef _IQTREE_MPI
- //------ NON-BLOCKING COMMUNICATION
- // FIX BUG: send candidateTrees instead of intermediateTrees
- candidateTrees.getAllTrees(trees, scores, WT_TAXON_ID + WT_BR_LEN + WT_BR_LEN_SHORT);
- // Send all trees to other processes
- MPIHelper::getInstance().distributeTrees(trees, scores, TREE_TAG);
-
- // Get trees from other nodes
- cout << "Getting top candidate trees from other processes ... " << endl;
- MPI_CollectTrees(true, maxNumTrees, true);
- MPI_Barrier(MPI_COMM_WORLD);
-#endif
-*/
-
-// #ifdef _IQTREE_MPI
-// // Send trees
-// MPIHelper::getInstance().distributeTrees(nniTrees, nniScores, TREE_TAG);
-// MPI_Barrier(MPI_COMM_WORLD);
-// // Receive trees
-// maxNumTrees = treesPerProc * (MPIHelper::getInstance().getNumProcesses() - 1);
-// MPI_CollectTrees(true,maxNumTrees,true);
-// #endif
-// if (params->fixStableSplits && candidateTrees.size() > 1) {
-// candidateTrees.computeSplitOccurences(Params::getInstance().stableSplitThreshold, Params::getInstance().numSupportTrees);
-// }
}
string IQTree::generateParsimonyTree(int randomSeed) {
@@ -2212,43 +2159,6 @@ string IQTree::optimizeBranches(int maxTraversal) {
return tree;
}
-void IQTree::collectBootTrees() {
-#ifdef _IQTREE_MPI
- if (boot_trees.size() == 0)
- return;
- // send UFBoot trees between processes
- if (MPIHelper::getInstance().isMaster()) {
- MPIHelper::getInstance().sendMsg(BOOT_TAG, "BOOT TREES PLEASE!");
- TreeCollection trees;
- int count = 0;
- do {
- int source = MPIHelper::getInstance().receiveTrees(trees, BOOT_TREE_TAG);
- if (source > 0) {
- count++;
- ASSERT(trees.getNumTrees() == boot_trees.size());
- int better_trees = 0;
- for (int id = 0; id < trees.getNumTrees(); id++)
- if (trees.getScores()[id] > boot_logl[id]) {
- boot_trees[id] = trees.getTreeStrings()[id];
- boot_logl[id] = trees.getScores()[id];
- better_trees++;
- }
- trees.clear();
- cout << better_trees << " better bootstrap trees from process " << source << endl;
- }
- } while (count < MPIHelper::getInstance().getNumProcesses()-1);
- } else {
- // worker
- if (MPIHelper::getInstance().checkMsg(BOOT_TAG))
- MPIHelper::getInstance().sendTrees(PROC_MASTER, boot_trees, boot_logl, BOOT_TREE_TAG);
- string msg;
- if (MPIHelper::getInstance().checkMsg(LOGL_CUTOFF_TAG, msg)) {
- logl_cutoff = convert_double(msg.c_str());
- cout << "Log-likelihood cutoff on original alignment: " << logl_cutoff << endl;
- }
- }
-#endif
-}
double IQTree::doTreeSearch() {
cout << "--------------------------------------------------------------------" << endl;
@@ -2344,24 +2254,6 @@ double IQTree::doTreeSearch() {
while (!stop_rule.meetStopCondition(stop_rule.getCurIt(), cur_correlation)) {
-/*
-#ifdef _IQTREE_MPI
- // check stopping rule
- if (MPIHelper::getInstance().isMaster()) {
- if (stop_rule.meetStopCondition(stop_rule.getCurIt(), cur_correlation)) {
- MPIHelper::getInstance().sendMsg(STOP_TAG, "STOP!");
- break;
- }
- } else {
- if(MPIHelper::getInstance().checkMsg(STOP_TAG)) {
- break;
- }
- }
-#else
- if (stop_rule.meetStopCondition(stop_rule.getCurIt(), cur_correlation))
- break;
-#endif
-*/
searchinfo.curIter = stop_rule.getCurIt();
// estimate logl_cutoff for bootstrap
if (!boot_orig_logl.empty())
@@ -2393,26 +2285,6 @@ double IQTree::doTreeSearch() {
if (MPIHelper::getInstance().isWorker() || MPIHelper::getInstance().gotMessage())
syncCurrentTree();
-/*
-#ifdef _IQTREE_MPI
- //----------- NON-BLOCKING COMMUNICATION ---------//
- int maxNumTrees = (MPIHelper::getInstance().getNumProcesses() - 1) * 2;
- if (MPIHelper::getInstance().isMaster()) {
- // master: receive tree from WORKERS
- bool candidateset_changed = MPI_CollectTrees(false, maxNumTrees, true);
- if (candidateset_changed) {
- vector<string> bestTrees = candidateTrees.getBestTreeStrings(Params::getInstance().popSize);
- vector<double> bestScores = candidateTrees.getBestScores(Params::getInstance().popSize);
- MPIHelper::getInstance().distributeTrees(bestTrees, bestScores, TREE_TAG);
- }
- } else {
- // worker: always send tree to MASTER
- MPIHelper::getInstance().sendTree(PROC_MASTER, getTreeString(), curScore, TREE_TAG);
- MPI_CollectTrees(false, maxNumTrees, true);
- }
-#endif
-*/
-
// TODO: cannot check yet, need to somehow return treechanged
// if (nni_count == 0 && params->snni && numPerturb > 0 && treechanged) {
@@ -2453,14 +2325,10 @@ double IQTree::doTreeSearch() {
* convergence criterion for ultrafast bootstrap
*---------------------------------------*/
- // workers send bootstrap trees, TODO: blocking communication
-// if (params->stop_condition == SC_BOOTSTRAP_CORRELATION && MPIHelper::getInstance().isWorker())
-// collectBootTrees();
// MASTER receives bootstrap trees and perform stop convergence test
if ((stop_rule.getCurIt()) >= ufboot_count &&
params->stop_condition == SC_BOOTSTRAP_CORRELATION && MPIHelper::getInstance().isMaster()) {
-// collectBootTrees();
ufboot_count += params->step_iterations/2;
// compute split support every half step
SplitGraph *sg = new SplitGraph;
@@ -2501,6 +2369,14 @@ double IQTree::doTreeSearch() {
}
+ // 2019-06-03: check convergence here to avoid effect of refineBootTrees
+ if (boot_splits.size() >= 2 && MPIHelper::getInstance().isMaster()) {
+ // check the stopping criterion for ultra-fast bootstrap
+ if (computeBootstrapCorrelation() < params->min_correlation)
+ cout << "WARNING: bootstrap analysis did not converge. You should rerun with higher number of iterations (-nm option)" << endl;
+
+ }
+
if(params->ufboot2corr) refineBootTrees();
if (optimization_looped)
@@ -2527,10 +2403,6 @@ double IQTree::doTreeSearch() {
cout << "Total number of trees sent: " << MPIHelper::getInstance().getNumTreeSent() << endl;
cout << "Total number of NNI searches done by myself: " << MPIHelper::getInstance().getNumNNISearch() << endl;
MPIHelper::getInstance().resetNumbers();
-// MPI_Finalize();
-// if (MPIHelper::getInstance().getProcessID() != MASTER) {
-// exit(0);
-// }
#endif
@@ -2774,6 +2646,8 @@ void IQTree::refineBootTrees() {
// 2018-08-17: delete duplicated memory
deleteAllPartialLh();
+ ModelsBlock *models_block = readModelsDefinition(*params);
+
// do bootstrap analysis
for (int sample = refined_samples; sample < boot_trees.size(); sample++) {
// create bootstrap alignment
@@ -2814,9 +2688,18 @@ void IQTree::refineBootTrees() {
boot_tree->setParams(params);
- // copy model
- boot_tree->setModelFactory(getModelFactory());
+ // 2019-06-03: bug fix setting part_info properly
+ if (boot_tree->isSuperTree())
+ ((PhyloSuperTree*)boot_tree)->setPartInfo((PhyloSuperTree*)this);
+ // copy model
+ // BQM 2019-05-31: bug fix with -bsam option
+ boot_tree->initializeModel(*params, aln->model_name, models_block);
+ boot_tree->getModelFactory()->setCheckpoint(getCheckpoint());
+ if (isSuperTree())
+ ((PartitionModel*)boot_tree->getModelFactory())->PartitionModel::restoreCheckpoint();
+ else
+ boot_tree->getModelFactory()->restoreCheckpoint();
// set likelihood kernel
boot_tree->setParams(params);
@@ -2824,7 +2707,15 @@ void IQTree::refineBootTrees() {
boot_tree->setNumThreads(num_threads);
// load the current ufboot tree
- boot_tree->readTreeString(boot_trees[sample]);
+ // 2019-02-06: fix crash with -sp and -bnni
+ boot_tree->PhyloTree::readTreeString(boot_trees[sample]);
+
+ if (boot_tree->isSuperTree() && params->partition_type == BRLEN_OPTIMIZE) {
+ if (((PhyloSuperTree*)boot_tree)->size() > 1) {
+ // re-initialize branch lengths for unlinked model
+ boot_tree->wrapperFixNegativeBranch(true);
+ }
+ }
// TODO: check if this resolves the crash in reorientPartialLh()
boot_tree->initializeAllPartialLh();
@@ -2857,7 +2748,7 @@ void IQTree::refineBootTrees() {
// delete memory
- boot_tree->setModelFactory(NULL);
+ //boot_tree->setModelFactory(NULL);
boot_tree->save_all_trees = 2;
bootstrap_alignment = boot_tree->aln;
@@ -2878,6 +2769,8 @@ void IQTree::refineBootTrees() {
checkpoint->dump();
}
+
+ delete models_block;
cout << "Total " << refined_trees << " ufboot trees refined" << endl;
@@ -2953,48 +2846,6 @@ void IQTree::printIterationInfo(int sourceProcID) {
// }
//}
-#ifdef _IQTREE_MPI
-bool IQTree::MPI_CollectTrees(bool allTrees, int maxNumTrees, bool updateStopRule) {
- if (MPIHelper::getInstance().getNumProcesses() == 1)
- return false;
- TreeCollection outTrees;
- double start = getRealTime();
- MPIHelper::getInstance().receiveTrees(allTrees, maxNumTrees, outTrees, TREE_TAG);
- double commTime = getRealTime() - start;
- if (verbose_mode >= VB_MED && outTrees.getNumTrees()> 0) {
- cout << outTrees.getNumTrees() << " trees received from other processes in ";
- cout << commTime << " seconds" << endl;
- }
- if (commTime > 1.0) {
- cout << "WARNING: Communication time (" << commTime << " sec) is too slow. Please increase MP_BUFFER_MEM and MP_EAGER_LIMIT" << endl;
- }
-
-// PhyloTree phyloTree;
-// phyloTree.aln = this->aln;
-// phyloTree.setParams(&(Params::getInstance()));
-
- bool candidateset_changed = false;
-
- for (int i = 0; i < outTrees.getNumTrees(); i++) {
- pair<string, double> tree = outTrees.getTree(i);
- if (tree.first == "notree") {
- if (updateStopRule) {
- stop_rule.setCurIt(stop_rule.getCurIt() + 1);
- curScore = tree.second;
- cout << "Bad tree with score: " << tree.second << " skipped" << endl;
- printIterationInfo(outTrees.getSourceProcID()[i]);
- }
- } else {
-// phyloTree.readTreeString(tree.first, true);
-// string treeString = phyloTree.getTreeString();
- int pos = addTreeToCandidateSet(tree.first, tree.second, updateStopRule, outTrees.getSourceProcID()[i]);
- if (pos >= 0 && pos < params->popSize)
- candidateset_changed = true;
- }
- }
- return candidateset_changed;
-}
-#endif
double IQTree::doTreePerturbation() {
if (iqp_assess_quartet == IQP_BOOTSTRAP) {
@@ -3810,13 +3661,6 @@ void IQTree::summarizeBootstrap(Params ¶ms, MTreeSet &trees) {
if (verbose_mode >= VB_MED)
cout << sg.size() << " splits found" << endl;
- if (!boot_splits.empty()) {
- // check the stopping criterion for ultra-fast bootstrap
- if (computeBootstrapCorrelation() < params.min_correlation)
- cout << "WARNING: bootstrap analysis did not converge. You should rerun with higher number of iterations (-nm option)" << endl;
-
- }
-
sg.scaleWeight(1.0 / trees.sumTreeWeights(), false, 4);
string out_file;
out_file = params.out_prefix;
@@ -4333,7 +4177,7 @@ int PhyloTree::testNumThreads() {
#ifndef _OPENMP
return 1;
#else
- int max_procs = min(countPhysicalCPUCores(), params->num_threads_max);
+ int max_procs = min(countPhysicalCPUCores()/MPIHelper::getInstance().countSameHost(), params->num_threads_max);
cout << "Measuring multi-threading efficiency up to " << max_procs << " CPU cores" << endl;
DoubleVector runTimes;
int bestProc = 0;
=====================================
tree/iqtree.h
=====================================
@@ -310,8 +310,6 @@ public:
*/
double swapTaxa(PhyloNode *node1, PhyloNode *node2);
- /** collect boostrap trees from workers to master */
- void collectBootTrees();
/**
perform tree search
@@ -961,22 +959,6 @@ protected:
string generateParsimonyTree(int randomSeed);
-#ifdef _IQTREE_MPI
- /**
- * Receive trees from other processes and add them to the candidate set
- *
- * @param allTrees
- * If true, wait for tree from every node
- * If false, only collect trees that have been sent
- * @param maxNumTrees
- * Only received up to maxNumTrees to prevent the function to block because it can constantly receive
- * new trees
- * @param updateStopRule
- * To update the stop rule or not
- */
- bool MPI_CollectTrees(bool allTrees, int maxNumTrees, bool updateStopRule);
-#endif
-
double doTreePerturbation();
void estimateLoglCutoffBS();
=====================================
tree/phylokernelnew.h
=====================================
@@ -826,7 +826,7 @@ inline void scaleLikelihood(VectorClass &lh_max, double *invar, double *dad_part
// now do the likelihood scaling
double *partial_lh = &dad_partial_lh[x];
for (i = 0; i < nstates; i++)
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
dad_scale_num[x*ncat_mix] += 1;
}
}
@@ -840,7 +840,7 @@ inline void scaleLikelihood(VectorClass &lh_max, double *invar, double *dad_part
double *partial_lh = &dad_partial_lh[x];
// now do the likelihood scaling
for (i = 0; i < block; i++) {
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
}
dad_scale_num[x] += 1;
}
@@ -1512,7 +1512,7 @@ void PhyloTree::computePartialLikelihoodGenericSIMD(TraversalInfo &info, size_t
// now do the likelihood scaling
double *partial_lh = (double*)partial_lh_tmp + (x);
for (i = 0; i < nstates; i++)
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
dad_branch->scale_num[(ptn+x)*ncat_mix+c] += 1;
}
}
@@ -1530,7 +1530,7 @@ void PhyloTree::computePartialLikelihoodGenericSIMD(TraversalInfo &info, size_t
double *partial_lh = (double*)partial_lh_all + (x);
// now do the likelihood scaling
for (i = 0; i < block; i++) {
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
}
// sum_scale += LOG_SCALING_THRESHOLD * ptn_freq[ptn+x];
dad_branch->scale_num[ptn+x] += 1;
@@ -1729,7 +1729,7 @@ void PhyloTree::computePartialLikelihoodGenericSIMD(TraversalInfo &info, size_t
// now do the likelihood scaling
double *partial_lh = dad_branch->partial_lh + (ptn*block + c*nstates*VectorClass::size() + x);
for (i = 0; i < nstates; i++)
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
dad_branch->scale_num[(ptn+x)*ncat_mix+c] += 1;
}
}
@@ -1792,7 +1792,7 @@ void PhyloTree::computePartialLikelihoodGenericSIMD(TraversalInfo &info, size_t
// now do the likelihood scaling
double *partial_lh = dad_branch->partial_lh + (ptn*block + c*nstates*VectorClass::size() + x);
for (i = 0; i < nstates; i++)
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
dad_branch->scale_num[(ptn+x)*ncat_mix+c] += 1;
}
}
@@ -1811,7 +1811,7 @@ void PhyloTree::computePartialLikelihoodGenericSIMD(TraversalInfo &info, size_t
double *partial_lh = dad_branch->partial_lh + (ptn*block + x);
// now do the likelihood scaling
for (i = 0; i < block; i++) {
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
}
// sum_scale += LOG_SCALING_THRESHOLD * ptn_freq[ptn+x];
dad_branch->scale_num[ptn+x] += 1;
@@ -1915,7 +1915,7 @@ void PhyloTree::computePartialLikelihoodGenericSIMD(TraversalInfo &info, size_t
// now do the likelihood scaling
double *partial_lh = dad_branch->partial_lh + (ptn*block + c*nstates*VectorClass::size() + x);
for (i = 0; i < nstates; i++)
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
scale_dad[x*ncat_mix] += 1;
}
scale_dad++;
@@ -1936,7 +1936,7 @@ void PhyloTree::computePartialLikelihoodGenericSIMD(TraversalInfo &info, size_t
double *partial_lh = dad_branch->partial_lh + (ptn*block + x);
// now do the likelihood scaling
for (i = 0; i < block; i++) {
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
}
// sum_scale += LOG_SCALING_THRESHOLD * ptn_freq[ptn+x];
dad_branch->scale_num[ptn+x] += 1;
@@ -2464,7 +2464,7 @@ void PhyloTree::computeLikelihoodDervGenericSIMD(PhyloNeighbor *dad_branch, Phyl
// mixed branch length model
for (i = 0; i < nmixlen; i++) {
df[i] = horizontal_add(all_dfvec[i]);
- ASSERT(!std::isnan(df[i]) && !std::isinf(df[i]) && "Numerical underflow for lh-derivative");
+ ASSERT(std::isfinite(df[i]) && "Numerical underflow for lh-derivative");
}
for (i = 0; i < nmixlen2; i++)
ddf[i] = horizontal_add(all_ddfvec[i]);
@@ -2476,15 +2476,15 @@ void PhyloTree::computeLikelihoodDervGenericSIMD(PhyloNeighbor *dad_branch, Phyl
// normal joint branch length model
*df = horizontal_add(all_df);
*ddf = horizontal_add(all_ddf);
- if (std::isnan(*df) || std::isinf(*df)) {
+ if (!std::isfinite(*df)) {
getModel()->writeInfo(cout);
getRate()->writeInfo(cout);
}
- if (!SAFE_NUMERIC && (std::isnan(*df) || std::isinf(*df)))
+ if (!SAFE_NUMERIC && !std::isfinite(*df))
outError("Numerical underflow (lh-derivative). Run again with the safe likelihood kernel via `-safe` option");
-// ASSERT(!std::isnan(*df) && !std::isinf(*df) && "Numerical underflow for lh-derivative");
+// ASSERT(std::isfinite(*df) && "Numerical underflow for lh-derivative");
if (isASC) {
double prob_const = 0.0, df_const = 0.0, ddf_const = 0.0;
@@ -2500,7 +2500,7 @@ void PhyloTree::computeLikelihoodDervGenericSIMD(PhyloNeighbor *dad_branch, Phyl
*ddf += nsites *(ddf_frac + df_frac*df_frac);
}
- if (std::isnan(*df) || std::isinf(*df)) {
+ if (!std::isfinite(*df)) {
cout << "WARNING: Numerical underflow for lh-derivative" << endl;
*df = *ddf = 0.0;
}
@@ -2949,17 +2949,17 @@ double PhyloTree::computeLikelihoodBranchGenericSIMD(PhyloNeighbor *dad_branch,
tree_lh += horizontal_add(all_tree_lh);
- if (!SAFE_NUMERIC && (std::isnan(tree_lh)))
+ if (!SAFE_NUMERIC && !std::isfinite(tree_lh))
outError("Numerical underflow (lh-branch). Run again with the safe likelihood kernel via `-safe` option");
- if (std::isnan(tree_lh))
+ if (!std::isfinite(tree_lh))
outWarning("Numerical underflow for lh-branch");
// arbitrarily fix tree_lh if underflown for some sites
- if (std::isnan(tree_lh)) {
+ if (!std::isfinite(tree_lh)) {
tree_lh = 0.0;
for (ptn = 0; ptn < orig_nptn; ptn++) {
- if (std::isnan(_pattern_lh[ptn])) {
+ if (!std::isfinite(_pattern_lh[ptn])) {
_pattern_lh[ptn] = LOG_SCALING_THRESHOLD*4; // log(2^(-1024))
}
tree_lh += _pattern_lh[ptn] * ptn_freq[ptn];
@@ -2986,7 +2986,7 @@ double PhyloTree::computeLikelihoodBranchGenericSIMD(PhyloNeighbor *dad_branch,
(VectorClass().load_a(&_pattern_lh[ptn])-prob_const).store_a(&_pattern_lh[ptn]);
// _pattern_lh[ptn] -= prob_const;
tree_lh -= aln->getNSite()*prob_const;
- ASSERT(!std::isnan(tree_lh) && !std::isinf(tree_lh));
+ ASSERT(std::isfinite(tree_lh));
}
return tree_lh;
@@ -3153,16 +3153,16 @@ double PhyloTree::computeLikelihoodFromBufferGenericSIMD()
double tree_lh = horizontal_add(all_tree_lh);
- if (!safe_numeric && (std::isnan(tree_lh) || std::isinf(tree_lh)))
+ if (!safe_numeric && !std::isfinite(tree_lh))
outError("Numerical underflow (lh-from-buffer). Run again with the safe likelihood kernel via `-safe` option");
- ASSERT(!std::isnan(tree_lh) && "Numerical underflow for lh-from-buffer");
+ ASSERT(std::isfinite(tree_lh) && "Numerical underflow for lh-from-buffer");
// arbitrarily fix tree_lh if underflown for some sites
- if (std::isinf(tree_lh)) {
+ if (!std::isfinite(tree_lh)) {
tree_lh = 0.0;
for (ptn = 0; ptn < orig_nptn; ptn++) {
- if (std::isinf(_pattern_lh[ptn])) {
+ if (!std::isfinite(_pattern_lh[ptn])) {
_pattern_lh[ptn] = LOG_SCALING_THRESHOLD*4; // log(2^(-1024))
}
tree_lh += _pattern_lh[ptn] * ptn_freq[ptn];
@@ -3189,7 +3189,7 @@ double PhyloTree::computeLikelihoodFromBufferGenericSIMD()
(VectorClass().load_a(&_pattern_lh[ptn])-prob_const).store_a(&_pattern_lh[ptn]);
// _pattern_lh[ptn] -= prob_const;
tree_lh -= aln->getNSite()*prob_const;
- ASSERT(!std::isnan(tree_lh) && !std::isinf(tree_lh));
+ ASSERT(std::isfinite(tree_lh));
}
return tree_lh;
@@ -3387,7 +3387,7 @@ void PhyloTree::computeLikelihoodDervMixlenGenericSIMD(PhyloNeighbor *dad_branch
df = horizontal_add(all_df);
ddf = horizontal_add(all_ddf);
- if (!SAFE_NUMERIC && (std::isnan(df) || std::isinf(df)))
+ if (!SAFE_NUMERIC && !std::isfinite(df))
outError("Numerical underflow (lh-derivative-mixlen). Run again with the safe likelihood kernel via `-safe` option");
if (isASC) {
@@ -3403,7 +3403,7 @@ void PhyloTree::computeLikelihoodDervMixlenGenericSIMD(PhyloNeighbor *dad_branch
ddf += nsites * (ddf_const + df_const*df_const);
}
- if (std::isnan(df) || std::isinf(df)) {
+ if (!std::isfinite(df)) {
cout << "WARNING: Numerical underflow for lh-derivative-mixlen" << endl;
df = ddf = 0.0;
}
=====================================
tree/phylokernelnonrev.h
=====================================
@@ -172,7 +172,7 @@ void PhyloTree::computeNonrevPartialLikelihoodGenericSIMD(TraversalInfo &info, s
double *partial_lh = dad_branch->partial_lh + (ptn*block + x);
// now do the likelihood scaling
for (i = 0; i < block; i++) {
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
}
dad_branch->scale_num[ptn+x] += 1;
}
@@ -360,7 +360,7 @@ void PhyloTree::computeNonrevPartialLikelihoodGenericSIMD(TraversalInfo &info, s
double *partial_lh = dad_branch->partial_lh + (ptn*block + x);
// now do the likelihood scaling
for (i = 0; i < block; i++) {
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
}
dad_branch->scale_num[ptn+x] += 1;
}
@@ -408,7 +408,7 @@ void PhyloTree::computeNonrevPartialLikelihoodGenericSIMD(TraversalInfo &info, s
double *partial_lh = dad_branch->partial_lh + (ptn*block + x);
// now do the likelihood scaling
for (i = 0; i < block; i++) {
- partial_lh[i*VectorClass::size()] *= SCALING_THRESHOLD_INVER;
+ partial_lh[i*VectorClass::size()] = ldexp(partial_lh[i*VectorClass::size()], SCALING_THRESHOLD_EXP);
}
dad_branch->scale_num[ptn+x] += 1;
}
@@ -738,7 +738,7 @@ void PhyloTree::computeNonrevLikelihoodDervGenericSIMD(PhyloNeighbor *dad_branch
*df = horizontal_add(all_df);
*ddf = horizontal_add(all_ddf);
- ASSERT(!std::isnan(*df) && !std::isinf(*df) && "Numerical underflow for non-rev lh-derivative");
+ ASSERT(std::isfinite(*df) && "Numerical underflow for non-rev lh-derivative");
if (isASC) {
double prob_const = 0.0, df_const = 0.0, ddf_const = 0.0;
@@ -1075,7 +1075,7 @@ double PhyloTree::computeNonrevLikelihoodBranchGenericSIMD(PhyloNeighbor *dad_br
tree_lh = horizontal_add(all_tree_lh);
- if (std::isnan(tree_lh) || std::isinf(tree_lh)) {
+ if (!std::isfinite(tree_lh)) {
model->writeInfo(cout);
site_rate->writeInfo(cout);
ASSERT(0 && "Numerical underflow for non-rev lh-branch");
@@ -1094,7 +1094,7 @@ double PhyloTree::computeNonrevLikelihoodBranchGenericSIMD(PhyloNeighbor *dad_br
for (ptn = 0; ptn < orig_nptn; ptn+=VectorClass::size())
(VectorClass().load_a(&_pattern_lh[ptn])-prob_const).store_a(&_pattern_lh[ptn]);
tree_lh -= aln->getNSite()*prob_const;
- ASSERT(!std::isnan(tree_lh) && !std::isinf(tree_lh));
+ ASSERT(std::isfinite(tree_lh));
}
return tree_lh;
=====================================
tree/phylosupertree.cpp
=====================================
@@ -95,6 +95,33 @@ void PhyloSuperTree::setModelFactory(ModelFactory *model_fac) {
}
}
+void PhyloSuperTree::setPartInfo(PhyloSuperTree *tree) {
+ part_info = tree->part_info;
+ int part = 0;
+ for (iterator it = begin(); it != end(); it++, part++) {
+ part_info[part].evalNNIs = 0.0;
+ }
+ if (!params->bootstrap_spec) {
+ return;
+ }
+ // 2018-06-03: for -bsam GENE, number of partitions might be reduced
+ if (part_info.size() <= size())
+ return;
+ part_info.erase(part_info.begin()+size(), part_info.end());
+ for (part = 0; part < part_info.size(); part++) {
+ bool found = false;
+ for (int p = 0; p < tree->size(); p++)
+ if (tree->at(p)->aln->name == at(part)->aln->name) {
+ part_info[part] = tree->part_info[p];
+ part_info[part].evalNNIs = 0.0;
+ found = true;
+ break;
+ }
+ ASSERT(found);
+ }
+}
+
+
void PhyloSuperTree::setSuperAlignment(Alignment *alignment) {
PhyloTree::setAlignment(alignment);
=====================================
tree/phylosupertree.h
=====================================
@@ -75,6 +75,12 @@ public:
*/
virtual void setModelFactory(ModelFactory *model_fac);
+ /**
+ 2019-06-03: copy part_info from tree, taking into account -bsam option
+ @param tree input super tree
+ */
+ void setPartInfo(PhyloSuperTree *tree);
+
/**
Set the alignment, important to compute parsimony or likelihood score
Assing taxa ids according to their position in the alignment
=====================================
tree/phylotree.h
=====================================
@@ -64,9 +64,14 @@ const double TOL_LIKELIHOOD_PARAMOPT = 0.001; // BQM: newly introduced for Model
//const static double LOG_SCALING_THRESHOLD = log(SCALING_THRESHOLD);
//#include "pll/pll.h"
// 2^256
-#define SCALING_THRESHOLD_INVER 115792089237316195423570985008687907853269984665640564039457584007913129639936.0
-#define SCALING_THRESHOLD (1.0/SCALING_THRESHOLD_INVER)
-#define LOG_SCALING_THRESHOLD log(SCALING_THRESHOLD)
+//#define SCALING_THRESHOLD_INVER 115792089237316195423570985008687907853269984665640564039457584007913129639936.0
+#define SCALING_THRESHOLD_EXP 256
+//#define SCALING_THRESHOLD (1.0/SCALING_THRESHOLD_INVER)
+// 2^{-256}
+#define SCALING_THRESHOLD 8.636168555094444625386e-78
+//#define SCALING_THRESHOLD ldexp(1.0, -256)
+//#define LOG_SCALING_THRESHOLD log(SCALING_THRESHOLD)
+#define LOG_SCALING_THRESHOLD -177.4456782233459932741
const int SPR_DEPTH = 2;
@@ -444,9 +449,9 @@ public:
*/
virtual void setAlignment(Alignment *alignment);
- /** set the root by name
- @param my_root root node name
- @param multi_taxa TRUE if my_root is a comma-separated list of nodes
+ /** set the root by name
+ @param my_root root node name
+ @param multi_taxa TRUE if my_root is a comma-separated list of nodes
*/
void setRootNode(const char *my_root, bool multi_taxa = false);
@@ -1759,8 +1764,8 @@ public:
virtual void changeLikelihoodKernel(LikelihoodKernel lk);
virtual void setLikelihoodKernel(LikelihoodKernel lk);
-
- virtual void setNumThreads(int num_threads);
+
+ virtual void setNumThreads(int num_threads);
#if defined(BINARY32) || defined(__NOAVX__)
void setLikelihoodKernelAVX() {}
@@ -1790,9 +1795,9 @@ public:
* Variance matrix
*/
double *var_matrix;
-
- /** distance matrix file */
- string dist_file;
+
+ /** distance matrix file */
+ string dist_file;
/**
TRUE if you want to optimize branch lengths by Newton-Raphson method
@@ -1958,13 +1963,13 @@ public:
@param partid partition ID as first column of the line. -1 to omit it
*/
virtual void writeSiteLh(ostream &out, SiteLoglType wsl, int partid = -1);
-
- /**
- write branches into a csv file
- Feature requested by Rob Lanfear
- @param out output stream
- */
- virtual void writeBranches(ostream &out);
+
+ /**
+ write branches into a csv file
+ Feature requested by Rob Lanfear
+ @param out output stream
+ */
+ virtual void writeBranches(ostream &out);
protected:
=====================================
tree/phylotreepars.cpp
=====================================
@@ -397,6 +397,11 @@ int PhyloTree::computeParsimonyTree(const char *out_prefix, Alignment *alignment
UINT *tmp_partial_pars;
tmp_partial_pars = newBitsBlock();
+ best_pars_score = 0;
+ if (leafNum == size) {
+ outWarning("Constraint tree has all taxa and is bifurcating, which strictly enforces final tree!");
+ }
+
// stepwise adding the next taxon for the remaining taxa
for (; leafNum < size; leafNum++) {
if (verbose_mode >= VB_MAX)
=====================================
utils/MPIHelper.cpp
=====================================
@@ -18,254 +18,70 @@ MPIHelper& MPIHelper::getInstance() {
return instance;
}
-void MPIHelper::distributeTrees(vector<string> &treeStrings, vector<double> &scores, int tag) {
- if (getNumProcesses() == 1)
- return;
+void MPIHelper::init(int argc, char *argv[]) {
#ifdef _IQTREE_MPI
- vector<int> sourceProcID;
- sourceProcID.insert(sourceProcID.end(), scores.size(), getProcessID());
- TreeCollection outTrees(treeStrings, scores, sourceProcID);
- cleanUpMessages();
- for (int i = 0; i < getNumProcesses(); i++) {
- if (i != getProcessID()) {
- MPI_Request *request = new MPI_Request;
- ObjectStream *os = new ObjectStream(outTrees);
- MPI_Isend(os->getObjectData(), os->getDataLength(), MPI_CHAR, i, tag, MPI_COMM_WORLD, request);
- sentMessages.push_back(make_pair(request, os));
- int flag = 0;
- MPI_Status status;
- MPI_Test(request, &flag, &status);
- }
+ int n_tasks, task_id;
+ if (MPI_Init(&argc, &argv) != MPI_SUCCESS) {
+ outError("MPI initialization failed!");
}
- //numTreeSent += treeStrings.size();
+ MPI_Comm_size(MPI_COMM_WORLD, &n_tasks);
+ MPI_Comm_rank(MPI_COMM_WORLD, &task_id);
+ setNumProcesses(n_tasks);
+ setProcessID(task_id);
+ setNumTreeReceived(0);
+ setNumTreeSent(0);
+ setNumNNISearch(0);
#endif
}
-void MPIHelper::distributeTree(string treeString, double score, int tag) {
- if (getNumProcesses() == 1)
- return;
+void MPIHelper::finalize() {
#ifdef _IQTREE_MPI
- double start = getRealTime();
- vector<string> trees;
- vector<double> scores;
- trees.push_back(treeString);
- scores.push_back(score);
- distributeTrees(trees, scores, tag);
- if (verbose_mode >= VB_MED)
- cout << "Sent tree to other processes in " << getRealTime() - start << " seconds" << endl;
- numTreeSent++;
+ MPI_Finalize();
#endif
}
-void MPIHelper::sendTrees(int dest, vector<string> &treeStrings, vector<double> &scores, int tag) {
- if (getNumProcesses() == 1 || dest == getProcessID())
- return;
+void MPIHelper::syncRandomSeed() {
#ifdef _IQTREE_MPI
- vector<int> sourceProcID;
- sourceProcID.insert(sourceProcID.end(), scores.size(), getProcessID());
- TreeCollection outTrees(treeStrings, scores, sourceProcID);
- cleanUpMessages();
- MPI_Request *request = new MPI_Request;
- ObjectStream *os = new ObjectStream(outTrees);
- MPI_Isend(os->getObjectData(), os->getDataLength(), MPI_CHAR, dest, tag, MPI_COMM_WORLD, request);
- sentMessages.push_back(make_pair(request, os));
- numTreeSent += treeStrings.size();
-
- int flag = 0;
- MPI_Status status;
- MPI_Test(request, &flag, &status);
-#endif
-}
-
-void MPIHelper::sendTree(int dest, string treeString, double score, int tag) {
- if (getNumProcesses() == 1 || dest == getProcessID())
- return;
-#ifdef _IQTREE_MPI
- StrVector treeStrings;
- treeStrings.push_back(treeString);
- DoubleVector scores;
- scores.push_back(score);
- sendTrees(dest, treeStrings, scores, tag);
-#endif
-}
-
-int MPIHelper::sendRecvTrees(int dest, vector<string> &treeStrings, vector<double> &scores, int tag) {
- if (getNumProcesses() == 1 || dest == getProcessID())
- return tag;
-#ifdef _IQTREE_MPI
- double beginTime = getRealTime();
- // prepare message
- vector<int> sourceProcID;
- sourceProcID.insert(sourceProcID.end(), scores.size(), getProcessID());
- TreeCollection outTrees(treeStrings, scores, sourceProcID);
- ObjectStream *os = new ObjectStream(outTrees);
-
- // blocking send
- MPI_Send(os->getObjectData(), os->getDataLength(), MPI_CHAR, dest, tag, MPI_COMM_WORLD);
- numTreeSent += treeStrings.size();
- delete os;
-
- // blocking probe
- MPI_Status status;
- MPI_Probe(dest, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
- int msgCount;
- MPI_Get_count(&status, MPI_CHAR, &msgCount);
-
- // receive the message
- char *recvBuffer = new char[msgCount];
- MPI_Recv(recvBuffer, msgCount, MPI_CHAR, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, &status);
- treeStrings.clear();
- scores.clear();
-
- if (status.MPI_TAG != STOP_TAG) {
- os = new ObjectStream(recvBuffer, msgCount);
- TreeCollection curTrees = os->getTreeCollection();
- treeStrings = curTrees.getTreeStrings();
- scores = curTrees.getScores();
- numTreeReceived += treeStrings.size();
+ unsigned int rndSeed;
+ if (MPIHelper::getInstance().isMaster()) {
+ rndSeed = Params::getInstance().ran_seed;
+ }
+ // Broadcast random seed
+ MPI_Bcast(&rndSeed, 1, MPI_INT, PROC_MASTER, MPI_COMM_WORLD);
+ if (MPIHelper::getInstance().isWorker()) {
+ // Params::getInstance().ran_seed = rndSeed + task_id * 100000;
+ Params::getInstance().ran_seed = rndSeed;
+ // printf("Process %d: random_seed = %d\n", task_id, Params::getInstance().ran_seed);
}
- delete [] recvBuffer;
-
- double endTime = getRealTime();
- cout << "INFO: " << endTime - beginTime << " seconds for " << __func__ << endl;
-
- return status.MPI_TAG;
-#else
- return tag;
#endif
}
-int MPIHelper::recvSendTrees(vector<string> &treeStrings, vector<double> &scores, vector<bool> &should_send, int tag) {
- if (getNumProcesses() == 1)
- return 0;
+int MPIHelper::countSameHost() {
#ifdef _IQTREE_MPI
- double beginTime = getRealTime();
- // blocking probe
- MPI_Status status;
- MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
- int msgCount;
- MPI_Get_count(&status, MPI_CHAR, &msgCount);
- int dest = status.MPI_SOURCE;
-
- // receive the message
- char *recvBuffer = new char[msgCount];
- MPI_Recv(recvBuffer, msgCount, MPI_CHAR, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, &status);
-
- // now send message
- if (!should_send[dest]) {
- treeStrings.resize(1, "notree");
- scores.resize(1, -DBL_MAX);
- }
- IntVector sourceProcID;
- sourceProcID.insert(sourceProcID.end(), scores.size(), getProcessID());
- TreeCollection outTrees(treeStrings, scores, sourceProcID);
- ObjectStream *os = new ObjectStream(outTrees);
-
- // blocking send
- MPI_Send(os->getObjectData(), os->getDataLength(), MPI_CHAR, dest, tag, MPI_COMM_WORLD);
- numTreeSent += treeStrings.size();
- delete os;
-
- // now extract trees from received buffer
- treeStrings.clear();
- scores.clear();
- os = new ObjectStream(recvBuffer, msgCount);
- TreeCollection curTrees = os->getTreeCollection();
- treeStrings = curTrees.getTreeStrings();
- scores = curTrees.getScores();
- delete [] recvBuffer;
- numTreeReceived += treeStrings.size();
+ // detect if processes are in the same host
+ char host_name[MPI_MAX_PROCESSOR_NAME];
+ int resultlen;
+ int pID = MPIHelper::getInstance().getProcessID();
+ MPI_Get_processor_name(host_name, &resultlen);
+ char *host_names;
+ host_names = new char[MPI_MAX_PROCESSOR_NAME * MPIHelper::getInstance().getNumProcesses()];
- should_send[dest] = false;
-
- double endTime = getRealTime();
- if (endTime - beginTime > 1)
- cout << "WARNING: " << endTime - beginTime << " seconds for " << __func__ << endl;
-
- return dest;
+ MPI_Allgather(host_name, resultlen+1, MPI_CHAR, host_names, MPI_MAX_PROCESSOR_NAME, MPI_CHAR,
+ MPI_COMM_WORLD);
+
+ int count = 0;
+ for (int i = 0; i < MPIHelper::getInstance().getNumProcesses(); i++)
+ if (strcmp(&host_names[i*MPI_MAX_PROCESSOR_NAME], host_name) == 0)
+ count++;
+ delete [] host_names;
+ if (count>1)
+ cout << "NOTE: " << count << " processes are running on the same host " << host_name << endl;
+ return count;
#else
- return 0;
-#endif
-}
-
-void MPIHelper::gatherTrees(TreeCollection &trees) {
- if (getNumProcesses() == 1)
- return;
-#ifdef _IQTREE_MPI
- double beginTime = getRealTime();
-
- if (isMaster()) {
- trees.clear();
- // Master: receive from all Workers
- for (int w = 1; w < getNumProcesses(); w++) {
- // blocking probe
- MPI_Status status;
- MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
- int msgCount;
- MPI_Get_count(&status, MPI_CHAR, &msgCount);
- // receive the message
- char *recvBuffer = new char[msgCount];
- MPI_Recv(recvBuffer, msgCount, MPI_CHAR, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, &status);
- ObjectStream *os = new ObjectStream(recvBuffer, msgCount);
- TreeCollection curTrees = os->getTreeCollection();
- trees.addTrees(curTrees);
- numTreeReceived += curTrees.getNumTrees();
- delete [] recvBuffer;
- }
- cout << trees.getNumTrees() << " trees gathered from workers in ";
- } else {
- // Worker: send trees to Master
- ObjectStream *os = new ObjectStream(trees);
- // blocking send
- MPI_Send(os->getObjectData(), os->getDataLength(), MPI_CHAR, PROC_MASTER, TREE_TAG, MPI_COMM_WORLD);
- numTreeSent += trees.getNumTrees();
- delete os;
- cout << trees.getNumTrees() << " trees sent to master in ";
- }
-
- double endTime = getRealTime();
- cout << endTime - beginTime << " seconds" << endl;
-#endif
-}
-
-void MPIHelper::broadcastTrees(TreeCollection &trees) {
- if (getNumProcesses() == 1)
- return;
-#ifdef _IQTREE_MPI
- double beginTime = getRealTime();
-
- // prepare data from Master
- ObjectStream *os;
- int msgCount = 0;
- if (isMaster()) {
- os = new ObjectStream(trees);
- msgCount = os->getDataLength();
- }
-
- // broadcast the count for workers
- MPI_Bcast(&msgCount, 1, MPI_INT, PROC_MASTER, MPI_COMM_WORLD);
-
- char *recvBuffer = new char[msgCount];
- if (isMaster())
- memcpy(recvBuffer, os->getObjectData(), msgCount);
-
- // broadcast trees to workers
- MPI_Bcast(recvBuffer, msgCount, MPI_CHAR, PROC_MASTER, MPI_COMM_WORLD);
-
- if (isWorker()) {
- os = new ObjectStream(recvBuffer, msgCount);
- trees = os->getTreeCollection();
- }
- delete os;
- delete [] recvBuffer;
-
- double endTime = getRealTime();
- cout << trees.getNumTrees() << " trees broadcasted to workers in " << endTime - beginTime << " seconds" << endl;
-
+ return 1;
#endif
}
-
bool MPIHelper::gotMessage() {
// Check for incoming messages
if (getNumProcesses() == 1)
@@ -283,170 +99,7 @@ bool MPIHelper::gotMessage() {
#endif
}
-void MPIHelper::sendMsg(int tag, string msg) {
- if (getNumProcesses() == 1)
- return;
-#ifdef _IQTREE_MPI
- if (tag == STOP_TAG)
- cleanUpMessages();
- for (int i = 0; i < getNumProcesses(); i++) {
- if (i != getProcessID()) {
- MPI_Request *request = new MPI_Request;
- ObjectStream *os = new ObjectStream(msg.c_str(), msg.size()+1);
- MPI_Isend(os->getObjectData(), os->getDataLength(), MPI_CHAR, i, tag, MPI_COMM_WORLD, request);
- sentMessages.push_back(make_pair(request, os));
- int flag = 0;
- MPI_Status status;
- MPI_Test(request, &flag, &status);
- }
- }
-#endif
-}
-
-bool MPIHelper::checkMsg(int tag, string &msg) {
- if (getNumProcesses() == 1)
- return true;
-#ifdef _IQTREE_MPI
- int flag=0;
- MPI_Status status;
- char *recvBuffer;
- int numBytes;
- // Check for incoming messages
- MPI_Iprobe(PROC_MASTER, tag, MPI_COMM_WORLD, &flag, &status);
- // flag == true if there is a message
- if (flag) {
- MPI_Get_count(&status, MPI_CHAR, &numBytes);
- recvBuffer = new char[numBytes];
- MPI_Recv(recvBuffer, numBytes, MPI_CHAR, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, &status);
- msg = recvBuffer;
- delete[] recvBuffer;
- return true;
- }
-#endif
- return false;
-}
-
-bool MPIHelper::checkMsg(int tag) {
- if (getNumProcesses() == 1) {
- return false;
- }
-#ifdef _IQTREE_MPI
- string msg;
- if (checkMsg(tag, msg)) {
- cout << "Worker " << getProcessID() << " gets message " << msg << endl;
- return true;
- }
-#endif
- return false;
-}
-
-
-void MPIHelper::receiveTrees(bool fromAll, int maxNumTrees, TreeCollection &trees, int tag) {
- if (getNumProcesses() == 1) {
- return;
- }
-#ifdef _IQTREE_MPI
- int flag = 0;
- int minNumTrees = 0;
- bool nodes[getNumProcesses()];
- if (fromAll)
- minNumTrees = getNumProcesses() - 1;
- for (int i = 0; i < getNumProcesses(); i++)
- nodes[i] = false;
- nodes[getProcessID()] = true;
- // Process all pending messages
- MPI_Status status;
- size_t totalMsgSize = 0;
- do {
- char* recvBuffer;
- int numBytes;
- flag = 0;
- // Check for incoming messages
- MPI_Iprobe(MPI_ANY_SOURCE, tag, MPI_COMM_WORLD, &flag, &status);
- // flag == true if there is a message
- if (flag) {
- //cout << "Getting messages from node " << status.MPI_SOURCE << endl;
- MPI_Get_count(&status, MPI_CHAR, &numBytes);
- totalMsgSize += numBytes;
- recvBuffer = new char[numBytes];
- MPI_Recv(recvBuffer, numBytes, MPI_CHAR, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, &status);
- ObjectStream os(recvBuffer, numBytes);
- if (status.MPI_TAG == STOP_TAG) {
- cout << os.getObjectData() << endl;
- MPI_Finalize();
- exit(0);
- }
- TreeCollection curTrees = os.getTreeCollection();
- trees.addTrees(curTrees);
- if (trees.getNumTrees() >= maxNumTrees) {
- break;
- }
- if (fromAll && !nodes[status.MPI_SOURCE]) {
- nodes[status.MPI_SOURCE] = true;
- minNumTrees--;
- }
- delete [] recvBuffer;
- }
- } while (minNumTrees > 0 || flag);
- numTreeReceived += trees.getNumTrees();
- if (trees.getNumTrees() > 0) {
- cout << "Proc " << getProcessID() << ": " << trees.getNumTrees() << " trees received from other processes (" << totalMsgSize << " bytes)" << endl;
- }
-#endif
-}
-int MPIHelper::receiveTrees(TreeCollection &trees, int tag) {
- if (getNumProcesses() == 1) {
- return -1;
- }
-#ifdef _IQTREE_MPI
- int flag = 0;
- // Process all pending messages
- MPI_Status status;
- char* recvBuffer;
- int numBytes;
- // Check for incoming messages
- MPI_Iprobe(MPI_ANY_SOURCE, tag, MPI_COMM_WORLD, &flag, &status);
- // flag == true if there is a message
- if (flag) {
- //cout << "Getting messages from node " << status.MPI_SOURCE << endl;
- MPI_Get_count(&status, MPI_CHAR, &numBytes);
- recvBuffer = new char[numBytes];
- MPI_Recv(recvBuffer, numBytes, MPI_CHAR, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, &status);
- ObjectStream os(recvBuffer, numBytes);
- TreeCollection curTrees = os.getTreeCollection();
- trees.addTrees(curTrees);
- delete [] recvBuffer;
- return status.MPI_SOURCE;
- }
-#endif
- return -1;
-}
-
-int MPIHelper::cleanUpMessages() {
-#ifdef _IQTREE_MPI
- int numMsgCleaned = 0;
- // change iterator to index to avoid iterator being invalidated after erase()
- for (int i = 0; i < sentMessages.size(); ) {
- int flag = 0;
- MPI_Status status;
- MPI_Test(sentMessages[i].first, &flag, &status);
- if (flag) {
- delete sentMessages[i].first;
- delete sentMessages[i].second;
- numMsgCleaned++;
- sentMessages.erase(sentMessages.begin()+i);
- } else {
- i++;
- }
- }
- if (verbose_mode >= VB_MED && numMsgCleaned)
- cout << numMsgCleaned << " messages sent and cleaned up" << endl;
- return numMsgCleaned;
-#else
- return 0;
-#endif
-}
#ifdef _IQTREE_MPI
void MPIHelper::sendString(string &str, int dest, int tag) {
=====================================
utils/MPIHelper.h
=====================================
@@ -24,14 +24,11 @@
#include <string>
#include <vector>
-#include "tools.h"
-#include "TreeCollection.h"
-#include "ObjectStream.h"
+#include "utils/tools.h"
+#include "utils/checkpoint.h"
#ifdef _IQTREE_MPI
-
#include <mpi.h>
-
#endif
#define PROC_MASTER 0
@@ -50,6 +47,12 @@ public:
*/
static MPIHelper &getInstance();
+ /** initialize MPI */
+ void init(int argc, char *argv[]);
+
+ /** finalize MPI */
+ void finalize();
+
/**
destructor
*/
@@ -79,112 +82,15 @@ public:
MPIHelper::processID = processID;
}
+ /** synchronize random seed from master to all workers */
+ void syncRandomSeed();
+
+ /** count the number of host with the same name as the current host */
+ int countSameHost();
+
/** @return true if got any message from another process */
bool gotMessage();
- /**
- * Receive trees that sent to the current process
- *
- * @param fromAll
- * wait until at least one tree from each remaining process has been received
- * @param maxNumTrees
- * Only received up to maxNumTrees to prevent the function to block because it can constantly receive
- * new trees
- * @param trees[OUT]
- * Trees received from other processes
- * @param tag MPI tag
- */
- void receiveTrees(bool fromAll, int maxNumTrees, TreeCollection &trees, int tag);
-
-
- /**
- * Receive trees that sent to the current process
- *
- * @param trees[OUT]
- * Trees received from other processes
- * @param tag MPI tag
- * @return source process ID
- */
- int receiveTrees(TreeCollection &trees, int tag);
-
- /**
- * Send trees to all other processes
- * @param treeStrings vector of trees
- * @param scores vector containing scores of the trees with same order as in treeStrings
- * @param tag used to classified the message
- */
- void distributeTrees(vector<string> &treeStrings, vector<double> &scores, int tag = TREE_TAG);
-
- /**
- * Similar to distributeTrees but only 1 tree is sent
- * @param treeString
- * @param score
- * @param tag
- */
- void distributeTree(string treeString, double score, int tag);
-
- /**
- * Send trees to a dest process
- * @param dest MPI rank of destination process
- * @param treeStrings vector of trees
- * @param scores vector containing scores of the trees with same order as in treeStrings
- * @param tag used to classified the message
- */
- void sendTrees(int dest, vector<string> &treeStrings, vector<double> &scores, int tag);
-
- /**
- * Send one tree to a dest process
- * @param dest MPI rank of destination process
- * @param treeString NEWICK tree string
- * @param score its score
- * @param tag used to classified the message
- */
- void sendTree(int dest, string treeString, double score, int tag);
-
- /**
- * Blocking Send and then receive trees with a dest process
- * @param dest MPI rank of destination process
- * @param[in,out] treeString NEWICK tree string
- * @param[in,out] score its score
- * @param tag used to classified the message
- * return the message tag
- */
- int sendRecvTrees(int dest, vector<string> &treeStrings, vector<double> &scores, int tag);
-
- /**
- * Blocking receive and then send trees with a dest process
- * @param dest MPI rank of destination process
- * @param[in,out] treeString NEWICK tree string
- * @param[in,out] score its score
- * @param tag used to classified the message
- * return the message tag
- */
- int recvSendTrees(vector<string> &treeStrings, vector<double> &scores, vector<bool> &should_send, int tag);
-
- /**
- gather trees from workers to master
- */
- void gatherTrees(TreeCollection &trees);
-
- /**
- broadcase trees from master to works
- */
- void broadcastTrees(TreeCollection &trees);
-
- /**
- * Send a message to other process, e.g. STOP_TAG
- */
- void sendMsg(int tag, string msg);
-
- /**
- * Check if a message is received, e.g. STOP_TAG
- */
- bool checkMsg(int tag);
-
- /**
- * Check if a message is received, e.g. STOP_TAG
- */
- bool checkMsg(int tag, string &msg);
/** wrapper for MPI_Send a string
@param str string to send
@@ -293,12 +199,6 @@ public:
private:
int numNNISearch;
-#ifdef _IQTREE_MPI
- // A list storing messages and the corresponding requests that have been sent from the current process.
- // When a message has been successfully received, it will be deleted from the list
- vector< pair<MPI_Request *, ObjectStream *> > sentMessages;
-#endif
-
};
=====================================
utils/tools.cpp
=====================================
@@ -1891,7 +1891,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
if (cnt >= argc)
throw "Use -n <#iterations>";
if (params.gbo_replicates != 0) {
- outError("Ultrafast bootstrap does not work with -n option");
+ throw("Ultrafast bootstrap does not work with -n option");
}
params.min_iterations = convert_int(argv[cnt]);
params.stop_condition = SC_FIXED_ITERATION;
@@ -2244,11 +2244,11 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
throw "Use -blmin <min_branch_length>";
params.min_branch_length = convert_double(argv[cnt]);
if (params.min_branch_length < 0.0)
- outError("Negative -blmin not allowed!");
+ throw("Negative -blmin not allowed!");
if (params.min_branch_length == 0.0)
- outError("Zero -blmin is not allowed due to numerical problems");
+ throw("Zero -blmin is not allowed due to numerical problems");
if (params.min_branch_length > 0.1)
- outError("-blmin must be < 0.1");
+ throw("-blmin must be < 0.1");
continue;
}
@@ -2258,7 +2258,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
throw "Use -blmax <max_branch_length>";
params.max_branch_length = convert_double(argv[cnt]);
if (params.max_branch_length < 0.5)
- outError("-blmax smaller than 0.5 is not allowed");
+ throw("-blmax smaller than 0.5 is not allowed");
continue;
}
if (strcmp(argv[cnt], "--show-lh") == 0) {
@@ -2746,7 +2746,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
if (cnt >= argc)
throw "Use -bb <#replicates>";
if (params.stop_condition == SC_FIXED_ITERATION) {
- outError("Ultrafast bootstrap does not work with -fast, -te or -n option");
+ throw("Ultrafast bootstrap does not work with -fast, -te or -n option");
}
params.gbo_replicates = convert_int(argv[cnt]);
// params.avoid_duplicated_trees = true;
@@ -2929,10 +2929,11 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
if (strcmp(argv[cnt], "-fast") == 0) {
// fast search option to resemble FastTree
if (params.gbo_replicates != 0) {
- outError("Ultrafast bootstrap (-bb) does not work with -fast option");
+ throw("Ultrafast bootstrap (-bb) does not work with -fast option");
}
params.numInitTrees = 2;
- params.min_iterations = 2;
+ if (params.min_iterations == -1)
+ params.min_iterations = 2;
params.stop_condition = SC_FIXED_ITERATION;
params.modelEps = 0.05;
continue;
@@ -2963,7 +2964,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
if (strcmp(argv[cnt], "--adt-pert") == 0) {
if (params.tabu == true) {
- outError("option -tabu and --adt-pert cannot be combined");
+ throw("option -tabu and --adt-pert cannot be combined");
}
params.adaptPertubation = true;
params.stableSplitThreshold = 1.0;
@@ -3057,7 +3058,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
continue;
}
if (strcmp(argv[cnt], "-pll") == 0) {
- outError("-pll option is discontinued.");
+ throw("-pll option is discontinued.");
params.pll = true;
continue;
}
@@ -3146,7 +3147,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
throw "Use -nni-eval <num_evaluation>";
params.nni5_num_eval = convert_int(argv[cnt]);
if (params.nni5_num_eval < 1)
- outError("Positive -nni-eval expected");
+ throw("Positive -nni-eval expected");
continue;
}
@@ -3156,7 +3157,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
throw "Use -bl-eval <num_evaluation>";
params.brlen_num_traversal = convert_int(argv[cnt]);
if (params.brlen_num_traversal < 1)
- outError("Positive -bl-eval expected");
+ throw("Positive -bl-eval expected");
continue;
}
@@ -3332,7 +3333,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
if (strcmp(argv[cnt], "-t") == 0 || strcmp(argv[cnt], "-te") == 0) {
if (strcmp(argv[cnt], "-te") == 0) {
if (params.gbo_replicates != 0) {
- outError("Ultrafast bootstrap does not work with -te option");
+ throw("Ultrafast bootstrap does not work with -te option");
}
params.min_iterations = 0;
params.stop_condition = SC_FIXED_ITERATION;
@@ -3399,7 +3400,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
throw "Use -mixlen <number of mixture branch lengths for heterotachy model>";
params.num_mixlen = convert_int(argv[cnt]);
if (params.num_mixlen < 1)
- outError("-mixlen must be >= 1");
+ throw("-mixlen must be >= 1");
continue;
}
@@ -3482,19 +3483,49 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
if (params.root != NULL && params.is_rooted)
throw "Not allowed to specify both -o <taxon> and -root";
+ if (params.model_test_and_tree && params.partition_type != BRLEN_OPTIMIZE)
+ throw("-mtree not allowed with edge-linked partition model (-spp or -q)");
+
+ if (params.do_au_test && params.topotest_replicates == 0)
+ throw("For AU test please specify number of bootstrap replicates via -zb option");
+
+ if (params.lh_mem_save == LM_MEM_SAVE && params.partition_file)
+ throw("-mem option does not work with partition models yet");
+
+ if (params.gbo_replicates && params.num_bootstrap_samples)
+ throw("UFBoot (-bb) and standard bootstrap (-b) must not be specified together");
+
+ if ((params.model_name.find("ONLY") != string::npos || (params.model_name.substr(0,2) == "MF" && params.model_name.substr(0,3) != "MFP")) && (params.gbo_replicates || params.num_bootstrap_samples))
+ throw("ModelFinder only cannot be combined with bootstrap analysis");
+
+ if (params.num_runs > 1 && params.treeset_file)
+ throw("Can't combine --runs and -z options");
+
+ if (params.num_runs > 1 && params.lmap_num_quartets >= 0)
+ throw("Can't combine --runs and -lmap options");
+
}
// try
catch (const char *str) {
- outError(str);
+ if (MPIHelper::getInstance().isMaster())
+ outError(str);
+ else
+ exit(EXIT_SUCCESS);
//} catch (char *str) {
//outError(str);
} catch (string str) {
- outError(str);
+ if (MPIHelper::getInstance().isMaster())
+ outError(str);
+ else
+ exit(EXIT_SUCCESS);
} catch (...) {
string err = "Unknown argument \"";
err += argv[cnt];
err += "\"";
- outError(err);
+ if (MPIHelper::getInstance().isMaster())
+ outError(err);
+ else
+ exit(EXIT_SUCCESS);
}
} // for
@@ -3510,26 +3541,6 @@ void parseArg(int argc, char *argv[], Params ¶ms) {
// 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)");
-
- if (params.do_au_test && params.topotest_replicates == 0)
- outError("For AU test please specify number of bootstrap replicates via -zb option");
-
- if (params.lh_mem_save == LM_MEM_SAVE && params.partition_file)
- outError("-mem option does not work with partition models yet");
-
- if (params.gbo_replicates && params.num_bootstrap_samples)
- outError("UFBoot (-bb) and standard bootstrap (-b) must not be specified together");
-
- if ((params.model_name.find("ONLY") != string::npos || (params.model_name.substr(0,2) == "MF" && params.model_name.substr(0,3) != "MFP")) && (params.gbo_replicates || params.num_bootstrap_samples))
- outError("ModelFinder only cannot be combined with bootstrap analysis");
-
- if (params.num_runs > 1 && params.treeset_file)
- outError("Can't combine --runs and -z options");
-
- if (params.num_runs > 1 && params.lmap_num_quartets >= 0)
- outError("Can't combine --runs and -lmap options");
// Diep:
if(params.ufboot2corr == true){
View it on GitLab: https://salsa.debian.org/med-team/iqtree/commit/93437d2cffc19f0180ade8a2003742306d8cd09e
--
View it on GitLab: https://salsa.debian.org/med-team/iqtree/commit/93437d2cffc19f0180ade8a2003742306d8cd09e
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/20190727/d6756ddc/attachment-0001.html>
More information about the debian-med-commit
mailing list