From d6c45559eea18e4c454ddea2deebaa57c6337a75 Mon Sep 17 00:00:00 2001 From: "George G. Vega Yon" Date: Mon, 2 Oct 2023 11:05:51 -0600 Subject: [PATCH] OpenMP finally working --- include/barry/model-bones.hpp | 6 +-- include/barry/model-meat.hpp | 40 ++-------------- .../models/geese/geese-meat-likelihood.hpp | 47 ++++++------------- tests/07-geese-flock.cpp | 1 - 4 files changed, 19 insertions(+), 75 deletions(-) diff --git a/include/barry/model-bones.hpp b/include/barry/model-bones.hpp index fd0dc1852..e454b86b3 100644 --- a/include/barry/model-bones.hpp +++ b/include/barry/model-bones.hpp @@ -293,11 +293,7 @@ class Model { * constant. */ ///@{ - double get_norm_const( - const std::vector< double > & params, - const size_t & i, - bool as_log = false - ); + std::vector< double > & get_normalizing_constants(); const std::vector< Array_Type > * get_pset( const size_t & i diff --git a/include/barry/model-meat.hpp b/include/barry/model-meat.hpp index 7a977bd2b..7aa8893ac 100644 --- a/include/barry/model-meat.hpp +++ b/include/barry/model-meat.hpp @@ -106,11 +106,6 @@ inline double likelihood_( double numerator = 0.0; // Computing the numerator - #if defined(__OPENMP) || defined(_OPENMP) - #pragma omp simd reduction(+:numerator) - #elif defined(__GNUC__) && !defined(__clang__) - #pragma GCC ivdep - #endif for (size_t j = 0u; j < params.size(); ++j) numerator += *(stats_target + j) * params[j]; @@ -970,39 +965,10 @@ inline double Model -inline double Model:: get_norm_const( - const std::vector & params, - const size_t & i, - bool as_log -) { - - // Checking if the index exists - if (i >= arrays2support.size()) - throw std::range_error("The requested support is out of range"); - - const auto id = arrays2support[i]; - - // Checking if we have updated the normalizing constant or not - if (!first_calc_done[id] || !vec_equal_approx(params, params_last[id]) ) - { - - first_calc_done[id] = true; - - size_t k = params.size() + 1u; - size_t n = stats_support[id].size() / k; - - normalizing_constants[id] = update_normalizing_constant( - params, &stats_support[id][0u], k, n - ); - - params_last[id] = params; - - } +inline std::vector< double > & +Model:: get_normalizing_constants() { - return as_log ? - std::log(normalizing_constants[id]) : - normalizing_constants[id] - ; + return normalizing_constants; } diff --git a/include/barry/models/geese/geese-meat-likelihood.hpp b/include/barry/models/geese/geese-meat-likelihood.hpp index 4b78e289a..d7fa7877c 100644 --- a/include/barry/models/geese/geese-meat-likelihood.hpp +++ b/include/barry/models/geese/geese-meat-likelihood.hpp @@ -29,27 +29,16 @@ inline double Geese::likelihood( model->update_normalizing_constants(par0); // Following the prunning sequence - std::vector< size_t > * preseq; - - if (use_reduced_sequence) - { - - preseq = &this->reduced_sequence; - - } - else - { - - preseq = &this->sequence; - - } + const std::vector< size_t > & preseq = use_reduced_sequence ? + this->reduced_sequence : this->sequence; // The first time it is called, it need to generate the corresponding // hashes of the columns so it is fast to access then (saves time // hashing and looking in the map.) - auto arrays2support = model->get_arrays2support(); + const auto arrays2support = model->get_arrays2support(); + const auto & normconst = model->get_normalizing_constants(); - for (auto& i : *preseq) + for (auto& i : preseq) { // We cannot compute probability at the leaf, we need to continue @@ -80,11 +69,10 @@ inline double Geese::likelihood( // 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; #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) + firstprivate(nfunctions, par0, node_offspring, array_id, normconst) #endif for (size_t n = 0u; n < psets.size(); ++n) // x = psets->begin(); x != psets->end(); ++x) { @@ -142,30 +130,25 @@ inline double Geese::likelihood( // Is this state valid? if (off_mult < 0.0) - { - - ++nstate; continue; - - } // Multiplying by P(x|x_n), the transition probability std::vector< double > temp_stats; temp_stats.reserve(par0.size()); for (auto p = 0u; p < par0.size(); ++p) - temp_stats.push_back(psets_stats[par0.size() * nstate + p]); - - nstate++; + temp_stats.push_back(psets_stats[par0.size() * n + p]); // Use try catch in the following line try { - off_mult *= model->likelihood( + + off_mult *= barry::likelihood_( + &temp_stats[0u], par0, - temp_stats, - array_id, - false, - true + normconst[arrays2support->operator[](array_id)], + par0.size(), + false ); + } catch (std::exception & e) { auto err = std::string(e.what()); @@ -221,7 +204,7 @@ inline double Geese::likelihood( // In the case that the sequence is empty, then it means // that we are looking at a completely unnanotated tree, // thus the likelihood should be one - if (preseq->size() == 0u) + if (preseq.size() == 0u) return as_log ? -std::numeric_limits::infinity() : 1.0; diff --git a/tests/07-geese-flock.cpp b/tests/07-geese-flock.cpp index dd984ce36..06ae7bf8d 100644 --- a/tests/07-geese-flock.cpp +++ b/tests/07-geese-flock.cpp @@ -1,4 +1,3 @@ -#undef _OPENMP #ifndef CATCH_CONFIG_MAIN // #include "/opt/intel/oneapi/advisor/2022.0.0/include/advisor-annotate.h" #include