New function to update norm const (all)

include/barry/model-bones.hpp

@@ -60,9 +60,9 @@ class Model {
      */
     ///@{
     std::vector< std::vector< double > > stats_support;          ///< Sufficient statistics of the model (support)
-    std::vector< size_t >                  stats_support_n_arrays; ///< Number of arrays included per support.
+    std::vector< size_t >                stats_support_n_arrays; ///< Number of arrays included per support.
     std::vector< std::vector< double > > stats_target;           ///< Target statistics of the model
-    std::vector< size_t >                  arrays2support;
+    std::vector< size_t >                arrays2support;
     ///@}
 
     /**
@@ -125,6 +125,14 @@ class Model {
     std::vector< std::string > transform_model_term_names;
 
 public:
+
+    /**
+     * @brief Computes the normalizing constant for a given set of parameters
+     * @details This function will compute the normalizing constant for a given
+     * set of parameters. It will also update the `normalizing_constants` member
+     * variable.
+    */
+    void update_normalizing_constants(const std::vector< double > & params);
 
     void set_rengine(std::mt19937 * rengine_, bool delete_ = false) {
 
@@ -240,34 +248,39 @@ class Model {
     double likelihood(
         const std::vector<double> & params,
         const size_t & i,
-        bool as_log = false
+        bool as_log = false,
+        bool no_update_normconst = false
     );
 
     double likelihood(
         const std::vector<double> & params,
         const Array_Type & Array_,
         int i = -1,
-        bool as_log = false
+        bool as_log = false,
+        bool no_update_normconst = false
     );
 
     double likelihood(
         const std::vector<double> & params,
         const std::vector<double> & target_,
         const size_t & i,
-        bool as_log = false
+        bool as_log = false,
+        bool no_update_normconst = false
     );
 
     double likelihood(
         const std::vector<double> & params,
         const double * target_,
         const size_t & i,
-        bool as_log = false
+        bool as_log = false,
+        bool no_update_normconst = false
     );
 
     double likelihood_total(
         const std::vector<double> & params,
         bool as_log = false,
-        BARRY_NCORES_ARG(=2)
+        BARRY_NCORES_ARG(=2),
+        bool no_update_normconst = false
     );
     ///@}
 

include/barry/model-meat.hpp

@@ -150,6 +150,37 @@ inline double likelihood_(
 
 }
 
+template <
+    typename Array_Type,
+    typename Data_Counter_Type,
+    typename Data_Rule_Type,
+    typename Data_Rule_Dyn_Type
+    >
+inline void Model<Array_Type, Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_Type>::update_normalizing_constants(
+    const std::vector< double > & params
+) {
+
+    #if defined(__OPENMP) || defined(_OPENMP)
+    #pragma omp parallel for firstprivate(params) \
+        shared(stats_support, normalizing_constants, first_calc_done)
+    #endif
+    for (size_t i = 0u; i < stats_support.size(); ++i)
+    {
+
+        size_t k = params.size() + 1u;
+        size_t n = stats_support[i].size() / k;
+
+        first_calc_done[i] = true;
+        normalizing_constants[i] = update_normalizing_constant(
+            params, &stats_support[i][0u], k, n
+        );
+
+    }
+
+    return;
+
+}
+
 template <
     typename Array_Type,
     typename Data_Counter_Type,
@@ -600,7 +631,8 @@ template <typename Array_Type, typename Data_Counter_Type, typename Data_Rule_Ty
 inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_Type>::likelihood(
     const std::vector<double> & params,
     const size_t & i,
-    bool as_log
+    bool as_log,
+    bool no_update_normconst
 ) {
 
     // Checking if the index exists
@@ -614,7 +646,7 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
         return as_log ? -std::numeric_limits<double>::infinity() : 0.0;
 
     // Checking if we have updated the normalizing constant or not
-    if (!first_calc_done[idx] || !vec_equal_approx(params, params_last[idx]) )
+    if (!no_update_normconst && (!first_calc_done[idx] || !vec_equal_approx(params, params_last[idx])))
     {
 
         first_calc_done[idx] = true;
@@ -645,7 +677,8 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     const std::vector<double> & params,
     const Array_Type & Array_,
     int i,
-    bool as_log
+    bool as_log,
+    bool no_update_normconst
 ) {
 
     // Key of the support set to use
@@ -691,7 +724,7 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
         target_ = transform_model_fun(&target_[0u], target_.size());
 
     // Checking if we have updated the normalizing constant or not
-    if (!first_calc_done[loc] || !vec_equal_approx(params, params_last[loc]) )
+    if (!no_update_normconst && (!first_calc_done[loc] || !vec_equal_approx(params, params_last[loc])) )
     {
 
         first_calc_done[loc] = true;
@@ -726,7 +759,8 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     const std::vector<double> & params,
     const std::vector<double> & target_,
     const size_t & i,
-    bool as_log
+    bool as_log,
+    bool no_update_normconst
 ) {
 
     // Checking if the index exists
@@ -759,7 +793,7 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     }
 
     // Checking if we have updated the normalizing constant or not
-    if (!first_calc_done[loc] || !vec_equal_approx(params, params_last[loc]) ) {
+    if (!no_update_normconst && (!first_calc_done[loc] || !vec_equal_approx(params, params_last[loc])) ) {
 
         first_calc_done[loc] = true;
 
@@ -789,7 +823,8 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     const std::vector<double> & params,
     const double * target_,
     const size_t & i,
-    bool as_log
+    bool as_log,
+    bool no_update_normconst
 ) {
 
     // Checking if the index exists
@@ -828,7 +863,7 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     }
 
     // Checking if we have updated the normalizing constant or not
-    if (!first_calc_done[loc] || !vec_equal_approx(params, params_last[loc]) ) {
+    if (!no_update_normconst && (!first_calc_done[loc] || !vec_equal_approx(params, params_last[loc]) )) {
 
         first_calc_done[loc] = true;
 
@@ -857,32 +892,42 @@ template <typename Array_Type, typename Data_Counter_Type, typename Data_Rule_Ty
 inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_Type>::likelihood_total(
     const std::vector<double> & params,
     bool as_log,
-    BARRY_NCORES_ARG()
+    BARRY_NCORES_ARG(),
+    bool no_update_normconst
 ) {
 
     size_t params_last_size = params_last.size();
 
     // #if defined(__OPENMP) || defined(_OPENMP)
     // #pragma omp parallel for num_threads(ncores)
     // #endif
-    for (size_t i = 0u; i < params_last_size; ++i)
-    {
 
-        if (!first_calc_done[i] || !vec_equal_approx(params, params_last[i]) )
+    if (!no_update_normconst)
+    {
+        #if defined(__OPENMP) || defined(_OPENMP)
+        #pragma omp parallel for num_threads(ncores) \
+            shared(normalizing_constants, params_last, first_calc_done, stats_support) \
+            firstprivate(params)
+        #endif
+        for (size_t i = 0u; i < params_last_size; ++i)
         {
 
-            size_t k = params.size() + 1u;
-            size_t n = stats_support[i].size() / k;
-
-            first_calc_done[i] = true;
-            normalizing_constants[i] = update_normalizing_constant(
-                params, &stats_support[i][0u], k, n
-            );
-
-            params_last[i] = params;
-
-        }
+            if (!first_calc_done[i] || !vec_equal_approx(params, params_last[i]) )
+            {
+
+                size_t k = params.size() + 1u;
+                size_t n = stats_support[i].size() / k;
+
+                first_calc_done[i] = true;
+                normalizing_constants[i] = update_normalizing_constant(
+                    params, &stats_support[i][0u], k, n
+                );
+
+                params_last[i] = params;
+
+            }
 
+        }
     }
 
     double res = 0.0;

include/barry/models/geese/flock-bones.hpp

@@ -19,7 +19,7 @@ class Flock {
     bool initialized = false;
 
     // Common components
-    std::mt19937              rengine;
+    std::mt19937 rengine;
     PhyloModel model = PhyloModel();
 
     Flock() {};
@@ -37,8 +37,8 @@ class Flock {
     size_t add_data(
         std::vector< std::vector<size_t> > & annotations,
         std::vector< size_t > &              geneid,
-        std::vector< int > &                       parent,
-        std::vector< bool > &                      duplication
+        std::vector< int > &                 parent,
+        std::vector< bool > &                duplication
     );
 
     /**

include/barry/models/geese/flock-meat.hpp

@@ -147,17 +147,20 @@ inline double Flock::likelihood_joint(
 
     double ans = as_log ? 0.0: 1.0;
 
+    std::vector< double > par0(par.begin(), par.end() - nfunctions);
+    model.update_normalizing_constants(par0);
+
     if (as_log) {
 
         for (auto& d : this->dat) 
-            ans += d.likelihood(par, as_log, use_reduced_sequence, ncores);
+            ans += d.likelihood(par, as_log, use_reduced_sequence, ncores, true);
 
     }
     else
     {
 
         for (auto& d : this->dat) 
-            ans *= d.likelihood(par, as_log, use_reduced_sequence, ncores);
+            ans *= d.likelihood(par, as_log, use_reduced_sequence, ncores, true);
 
     }
 

include/barry/models/geese/geese-bones.hpp

@@ -216,7 +216,8 @@ class Geese {
         const std::vector< double > & par,
         bool as_log = false,
         bool use_reduced_sequence = true,
-        size_t ncores = 1u
+        size_t ncores = 1u,
+        bool no_update_normalizing_constant = false
         );
 
     double likelihood_exhaust(const std::vector< double > & par);

include/barry/models/geese/geese-meat-likelihood.hpp

@@ -7,9 +7,11 @@ inline double Geese::likelihood(
     const std::vector< double > & par,
     bool as_log,
     bool use_reduced_sequence,
-    size_t ncores
+    size_t ncores,
+    bool no_update_normalizing_constant
 ) {
 
+    // Checking whether the model is initialized
     INITIALIZED()
 
     // Splitting the probabilities
@@ -22,6 +24,10 @@ inline double Geese::likelihood(
 
     double ll = 0.0;
 
+    // Updating normalizing constants
+    if (!no_update_normalizing_constant)
+        model->update_normalizing_constants(par0);
+
     // Following the prunning sequence
     std::vector< size_t > * preseq;
 
@@ -157,7 +163,8 @@ inline double Geese::likelihood(
                         par0,
                         temp_stats,
                         array_id,
-                        false
+                        false,
+                        true
                     );
                 } catch (std::exception & e) {
 

tests/10-geese-predict.cpp

@@ -1,3 +1,4 @@
+#undef _OPENMP
 #include "tests.hpp"
 
 #include "../include/barry/models/geese.hpp"