Fixing error in likelihood calc

USCbiostats · Sep 27, 2023 · e7bdf61 · e7bdf61
1 parent f0efb22
commit e7bdf61
Show file tree

Hide file tree

Showing 2 changed files with 44 additions and 5 deletions.
diff --git a/include/barry/model-meat.hpp b/include/barry/model-meat.hpp
@@ -87,9 +87,9 @@ inline double likelihood_(
         numerator += *(stats_target + j) * params[j];
 
     if (!log_)
-        numerator = exp(numerator BARRY_SAFE_EXP);
+        numerator = std::exp(numerator BARRY_SAFE_EXP);
     else
-        return numerator BARRY_SAFE_EXP - log(normalizing_constant);
+        return numerator BARRY_SAFE_EXP - std::log(normalizing_constant);
 
     double ans = numerator/normalizing_constant;
 
@@ -1309,8 +1309,7 @@ MODEL_TEMPLATE(Array_Type, sample)(
     } else { 
 
         probs.resize(pset_arrays[a].size());
-        std::vector< double > temp_stats;
-        temp_stats.reserve(params.size());
+        std::vector< double > temp_stats(params.size());
         const std::vector< double > & stats = pset_stats[a];
 
         int i_matches = -1;
@@ -1319,7 +1318,7 @@ MODEL_TEMPLATE(Array_Type, sample)(
 
             // Filling out the parameters
             for (auto p = 0u; p < params.size(); ++p)
-                temp_stats.push_back(stats[array * k + p]);
+                temp_stats[p] = stats[array * k + p];
 
             probs[array] = this->likelihood(params, temp_stats, i, false);
             cumprob += probs[array];

diff --git a/include/barry/models/geese/geese-types.hpp b/include/barry/models/geese/geese-types.hpp
@@ -1,5 +1,8 @@
 #ifndef GEESE_TYPES_HPP
 #define GEESE_TYPES_HPP
+
+#define POS(a,b) (b)*N + (a)
+
 /**
  * @name Convenient typedefs for Node objects.
  * */
@@ -114,4 +117,41 @@ typedef barry::Model<PhyloArray, PhyloCounterData, PhyloRuleData, PhyloRuleDynDa
 typedef barry::PowerSet<PhyloArray, PhyloRuleData> PhyloPowerSet;
 ///@}
 
+// template<> 
+// inline void PhyloArray::insert_cell(
+//     size_t i,
+//     size_t j,
+//     const Cell< size_t > & v,
+//     bool check_bounds,
+//     bool
+// ) {
+
+//     if (check_bounds)
+//         out_of_range(i,j); 
+
+//     auto & elptr = el[POS(i,j)];
+
+//     if (elptr == 0u)
+//     {
+
+//         el_rowsums[i] += v.value;
+//         el_colsums[j] += v.value;
+
+//     } 
+//     else
+//     {
+
+//         el_rowsums[i] += (v.value - elptr);
+//         el_colsums[j] += (v.value - elptr);
+
+//     }
+
+//     elptr = v.value;
+
+//     return;
+
+// }
+
+#undef POS
+
 #endif