diff --git a/include/barry/model-bones.hpp b/include/barry/model-bones.hpp index 80ecfbac6..b5afddd12 100644 --- a/include/barry/model-bones.hpp +++ b/include/barry/model-bones.hpp @@ -240,32 +240,28 @@ class Model { double likelihood( const std::vector & params, const size_t & i, - bool as_log = false, - BARRY_NCORES_ARG(=2) + bool as_log = false ); double likelihood( const std::vector & params, const Array_Type & Array_, int i = -1, - bool as_log = false, - BARRY_NCORES_ARG(=2) + bool as_log = false ); double likelihood( const std::vector & params, const std::vector & target_, const size_t & i, - bool as_log = false, - BARRY_NCORES_ARG(=2) + bool as_log = false ); double likelihood( const std::vector & params, const double * target_, const size_t & i, - bool as_log = false, - BARRY_NCORES_ARG(=2) + bool as_log = false ); double likelihood_total( diff --git a/include/barry/model-meat.hpp b/include/barry/model-meat.hpp index f1617272c..7cfa7caf9 100644 --- a/include/barry/model-meat.hpp +++ b/include/barry/model-meat.hpp @@ -20,11 +20,11 @@ inline double update_normalizing_constant( std::vector< double > resv(n, 0.0); - #if defined(__OPENMP) || defined(_OPENMP) - #pragma omp parallel for shared(resv) firstprivate(params, n, k) - #elif defined(__GNUC__) && !defined(__clang__) - #pragma GCC ivdep - #endif + // #if defined(__OPENMP) || defined(_OPENMP) + // #pragma omp parallel for shared(resv) firstprivate(params, n, k) + // #elif defined(__GNUC__) && !defined(__clang__) + // #pragma GCC ivdep + // #endif for (size_t j = 0u; j < (k - 1u); ++j) { @@ -600,13 +600,8 @@ template ::likelihood( const std::vector & params, const size_t & i, - bool as_log, - BARRY_NCORES_ARG() + bool as_log ) { - - #if defined(__OPENMP) || defined(_OPENMP) - omp_set_num_threads(ncores); - #endif // Checking if the index exists if (i >= arrays2support.size()) @@ -650,13 +645,8 @@ inline double Model & params, const Array_Type & Array_, int i, - bool as_log, - BARRY_NCORES_ARG() + bool as_log ) { - - #if defined(__OPENMP) || defined(_OPENMP) - omp_set_num_threads(ncores); - #endif // Key of the support set to use int loc; @@ -736,13 +726,8 @@ inline double Model & params, const std::vector & target_, const size_t & i, - bool as_log, - BARRY_NCORES_ARG() + bool as_log ) { - - #if defined(__OPENMP) || defined(_OPENMP) - omp_set_num_threads(ncores); - #endif // Checking if the index exists if (i >= arrays2support.size()) @@ -804,13 +789,8 @@ inline double Model & params, const double * target_, const size_t & i, - bool as_log, - BARRY_NCORES_ARG() + bool as_log ) { - - #if defined(__OPENMP) || defined(_OPENMP) - omp_set_num_threads(ncores); - #endif // Checking if the index exists if (i >= arrays2support.size()) @@ -879,13 +859,12 @@ inline double Model * preseq; @@ -60,38 +58,47 @@ inline double Geese::likelihood( { // Starting the prob - double totprob = 0.0; + size_t array_id = node.narray[s]; // Retrieving the sets of arrays - const std::vector< PhyloArray > * psets = - model->get_pset(node.narray[s]); + const std::vector< PhyloArray > & psets = + *(model->get_pset(array_id)); - const std::vector * psets_stats = - model->get_pset_stats(node.narray[s]); + const std::vector & psets_stats = + *(model->get_pset_stats(array_id)); std::vector< std::vector< size_t > > & locations = pset_loc[ - arrays2support->operator[](node.narray[s]) + arrays2support->operator[](array_id) ]; // Summation over all possible values of X + const auto & node_offspring = node.offspring; + std::vector< double > totprob_n(psets.size(), 0.0); size_t nstate = 0u; - size_t narray = 0u; - for (auto x = psets->begin(); x != psets->end(); ++x) + #if defined(_OPENMP) || defined(__OPENMP) + #pragma omp parallel for num_threads(ncores) \ + shared(locations, psets, psets_stats, totprob_n) \ + firstprivate(nfunctions, par0, node_offspring, array_id) + #endif + for (size_t n = 0u; n < psets.size(); ++n) // x = psets->begin(); x != psets->end(); ++x) { - if (!x->is_dense()) + // Retrieving the pset + const auto & x = psets[n]; + + if (!x.is_dense()) throw std::logic_error("This is only supported for dense arrays."); - std::vector< size_t > & location_x = locations[narray++]; + const std::vector< size_t > & location_x = locations[n]; // Extracting the possible values of each offspring double off_mult = 1.0; - for (auto o = 0u; o < x->ncol(); ++o) + for (auto o = 0u; o < x.ncol(); ++o) { // Setting the node - n_off = node.offspring[o]; + const Node * n_off = node_offspring[o]; // In the case that the offspring is a leaf, then we need to // check whether the state makes sense. @@ -102,7 +109,7 @@ inline double Geese::likelihood( if (n_off->annotations[f] != 9u) { - if (x->operator()(f, o) != n_off->annotations[f]) + if (x(f, o) != n_off->annotations[f]) { off_mult = -1.0; @@ -123,7 +130,7 @@ inline double Geese::likelihood( } // Retrieving the location to the respective set of probabilities - off_mult *= node.offspring[o]->subtree_prob[location_x[o]]; + off_mult *= n_off->subtree_prob[location_x[o]]; } @@ -140,7 +147,7 @@ inline double Geese::likelihood( std::vector< double > temp_stats; temp_stats.reserve(par0.size()); for (auto p = 0u; p < par0.size(); ++p) - temp_stats.push_back(psets_stats->operator[](par0.size() * nstate + p)); + temp_stats.push_back(psets_stats[par0.size() * nstate + p]); nstate++; @@ -149,9 +156,8 @@ inline double Geese::likelihood( off_mult *= model->likelihood( par0, temp_stats, - node.narray[s], - false, - ncores + array_id, + false ); } catch (std::exception & e) { @@ -171,12 +177,13 @@ inline double Geese::likelihood( } // Adding to the total probabilities - totprob += off_mult; + totprob_n[n] = off_mult; } // Setting the probability at the node - node.subtree_prob[s] = totprob; + for (size_t n = 0u; n < psets.size(); ++n) + node.subtree_prob[s] += totprob_n[n]; }