Skip to content
Commits on Source (6)
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ e5190a45
......@@ -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)
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=libomp -pthread")
set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -fopenmp=libomp")
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
View file @ e5190a45
......@@ -2943,14 +2943,23 @@ 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) {
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;
upper = d;
......@@ -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
View file @ e5190a45
......@@ -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);
}
......
debian/changelog
View file @ e5190a45
iqtree (1.6.11+dfsg-1) unstable; urgency=medium
* New upstream version
* debhelper-compat 12
* Standards-Version: 4.4.0
-- Andreas Tille <tille@debian.org> Sat, 27 Jul 2019 22:09:33 +0200
iqtree (1.6.9+dfsg-1) unstable; urgency=medium
* New upstream version
......
debian/control
View file @ e5190a45
......@@ -4,7 +4,7 @@ Uploaders: Andreas Tille <tille@debian.org>,
Kevin Murray <kdmfoss@gmail.com>
Section: science
Priority: optional
Build-Depends: debhelper (>= 11~),
Build-Depends: debhelper-compat (= 12),
cmake,
libeigen3-dev,
libsprng2-dev,
......@@ -13,7 +13,7 @@ Build-Depends: debhelper (>= 11~),
help2man,
time,
chrpath
Standards-Version: 4.3.0
Standards-Version: 4.4.0
Vcs-Browser: https://salsa.debian.org/med-team/iqtree
Vcs-Git: https://salsa.debian.org/med-team/iqtree.git
Homepage: http://www.cibiv.at/software/iqtree/
......
main/main.cpp
View file @ e5190a45
......@@ -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 &params) {
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() {
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,17 +2281,22 @@ int main(int argc, char *argv[]) {
if (checkpoint->hasKey("finished")) {
if (checkpoint->getBool("finished")) {
if (Params::getInstance().force_unfinished) {
if (MPIHelper::getInstance().isMaster())
cout << "NOTE: Continue analysis although a previous run already finished" << endl;
} else {
delete checkpoint;
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.");
delete checkpoint;
return EXIT_FAILURE;
else
exit(EXIT_SUCCESS);
exit(EXIT_FAILURE);
}
} else {
append_log = true;
}
} else {
if (MPIHelper::getInstance().isMaster())
outWarning("Ignore invalid checkpoint file " + filename);
checkpoint->clear();
}
......@@ -2318,25 +2315,9 @@ 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);
MPIHelper::getInstance().syncRandomSeed();
signal(SIGABRT, &funcAbort);
signal(SIGFPE, &funcAbort);
signal(SIGILL, &funcAbort);
......@@ -2469,8 +2450,10 @@ int main(int argc, char *argv[]) {
}
#endif
#ifdef _IQTREE_MPI
cout << endl << "MPI: " << MPIHelper::getInstance().getNumProcesses() << " processes";
#endif
#ifndef _IQTREE_MPI
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
View file @ e5190a45
......@@ -377,8 +377,14 @@ bool ModelDNA::setRateType(string rate_str) {
cout << rates[i] << " ";
cout << endl;
}
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
View file @ e5190a45
......@@ -579,6 +579,7 @@ ModelFactory::ModelFactory(Params &params, 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");
}
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 {
......
model/modelprotein.cpp
View file @ e5190a45
......@@ -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,6 +147,8 @@ 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
......@@ -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,7 +466,8 @@ 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
......@@ -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
View file @ e5190a45
......@@ -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;
......@@ -2879,6 +2770,8 @@ void IQTree::refineBootTrees() {
}
delete models_block;
cout << "Total " << refined_trees << " ufboot trees refined" << endl;
// restore randstream
......@@ -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 &params, 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
View file @ e5190a45
......@@ -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
View file @ e5190a45
......@@ -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
View file @ e5190a45
......@@ -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
View file @ e5190a45
......@@ -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
View file @ e5190a45
......@@ -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
View file @ e5190a45
......@@ -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;
......
tree/phylotreepars.cpp
View file @ e5190a45
......@@ -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
View file @ e5190a45
......@@ -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
unsigned int rndSeed;
if (MPIHelper::getInstance().isMaster()) {
rndSeed = Params::getInstance().ran_seed;
}
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
// 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);
}
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();
}
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;
// 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()];
// receive the message
char *recvBuffer = new char[msgCount];
MPI_Recv(recvBuffer, msgCount, MPI_CHAR, status.MPI_SOURCE, status.MPI_TAG, MPI_COMM_WORLD, &status);
MPI_Allgather(host_name, resultlen+1, MPI_CHAR, host_names, MPI_MAX_PROCESSOR_NAME, MPI_CHAR,
MPI_COMM_WORLD);
// 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();
should_send[dest] = false;
double endTime = getRealTime();
if (endTime - beginTime > 1)
cout << "WARNING: " << endTime - beginTime << " seconds for " << __func__ << endl;
return dest;
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;
return 1;
#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;
#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
View file @ e5190a45
......@@ -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;
}
/** @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);
/** synchronize random seed from master to all workers */
void syncRandomSeed();
/**
* Send a message to other process, e.g. STOP_TAG
*/
void sendMsg(int tag, string msg);
/** count the number of host with the same name as the current host */
int countSameHost();
/**
* Check if a message is received, e.g. STOP_TAG
*/
bool checkMsg(int tag);
/** @return true if got any message from another process */
bool gotMessage();
/**
* 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
};
......