Some barry speed improvements

USCbiostats · Sep 27, 2023 · a6adb22 · a6adb22
1 parent 2c18a6d
commit a6adb22
Show file tree

Hide file tree

Showing 7 changed files with 86 additions and 43 deletions.
diff --git a/inst/include/barry/barraydense-bones.hpp b/inst/include/barry/barraydense-bones.hpp
@@ -194,9 +194,9 @@ class BArrayDense {
 
     void rm_cell(size_t i, size_t j, bool check_bounds = true, bool check_exists = true);
 
-    void insert_cell(size_t i, size_t j, const Cell< Cell_Type > & v, bool check_bounds, bool check_exists);
+    void insert_cell(size_t i, size_t j, const Cell< Cell_Type > & v, bool check_bounds, bool);
     // void insert_cell(size_t i, size_t j, Cell< Cell_Type > && v, bool check_bounds, bool check_exists);
-    void insert_cell(size_t i, size_t j, Cell_Type v, bool check_bounds, bool check_exists);
+    void insert_cell(size_t i, size_t j, Cell_Type v, bool check_bounds, bool);
 
     void swap_cells(
         size_t i0, size_t j0, size_t i1, size_t j1, bool check_bounds = true,

diff --git a/inst/include/barry/barraydense-meat.hpp b/inst/include/barry/barraydense-meat.hpp
@@ -153,13 +153,16 @@ inline BArrayDense<Cell_Type, Data_Type>:: BArrayDense(
 ) : N(Array_.N), M(Array_.M){
 
     // Dimensions
-    el.resize(0u);
-    el_rowsums.resize(0u);
-    el_colsums.resize(0u);
+    el = Array_.el;
+    el_rowsums = Array_.el_rowsums;
+    el_colsums = Array_.el_colsums;
+    // el.resize(0u);
+    // el_rowsums.resize(0u);
+    // el_colsums.resize(0u);
 
-    std::copy(Array_.el.begin(), Array_.el.end(), std::back_inserter(el));
-    std::copy(Array_.el_rowsums.begin(), Array_.el_rowsums.end(), std::back_inserter(el_rowsums));
-    std::copy(Array_.el_colsums.begin(), Array_.el_colsums.end(), std::back_inserter(el_colsums));
+    // std::copy(Array_.el.begin(), Array_.el.end(), std::back_inserter(el));
+    // std::copy(Array_.el_rowsums.begin(), Array_.el_rowsums.end(), std::back_inserter(el_rowsums));
+    // std::copy(Array_.el_colsums.begin(), Array_.el_colsums.end(), std::back_inserter(el_colsums));
 
     // this->NCells  = Array_.NCells;
     this->visited = Array_.visited;
@@ -196,14 +199,17 @@ inline BArrayDense<Cell_Type,Data_Type> & BArrayDense<Cell_Type, Data_Type>::ope
     if (this != &Array_)
     {
 
-        el.resize(0u);
-        el_rowsums.resize(0u);
-        el_colsums.resize(0u);
+        el = Array_.el;
+        el_rowsums = Array_.el_rowsums;
+        el_colsums = Array_.el_colsums;
+        // el.resize(0u);
+        // el_rowsums.resize(0u);
+        // el_colsums.resize(0u);
 
-        // Entries
-        std::copy(Array_.el.begin(), Array_.el.end(), std::back_inserter(el));
-        std::copy(Array_.el_rowsums.begin(), Array_.el_rowsums.end(), std::back_inserter(el_rowsums));
-        std::copy(Array_.el_colsums.begin(), Array_.el_colsums.end(), std::back_inserter(el_colsums));
+        // // Entries
+        // std::copy(Array_.el.begin(), Array_.el.end(), std::back_inserter(el));
+        // std::copy(Array_.el_rowsums.begin(), Array_.el_rowsums.end(), std::back_inserter(el_rowsums));
+        // std::copy(Array_.el_colsums.begin(), Array_.el_colsums.end(), std::back_inserter(el_colsums));
 
 
         // this->NCells = Array_.NCells;
@@ -406,9 +412,10 @@ inline std::vector< Cell_Type > BArrayDense<Cell_Type, Data_Type>::get_row_vec (
     if (check_bounds) 
         out_of_range(i, 0u);
 
-    std::vector< Cell_Type > ans(ncol(), static_cast< Cell_Type >(false));
+    std::vector< Cell_Type > ans;
+    ans.reserve(ncol());
     for (size_t j = 0u; j < M; ++j) 
-        ans[j] = el[POS(i, j)];
+        ans.push_back(el[POS(i, j)]);
 
     return ans;
 
@@ -425,7 +432,7 @@ template<typename Cell_Type, typename Data_Type> inline void BArrayDense<Cell_Ty
         out_of_range(i, 0u);
 
     for (size_t j = 0u; j < M; ++j) 
-        x->at(j) = el[POS(i, j)];
+        x->operator[](j) = el[POS(i, j)];
 
 }
 
@@ -438,9 +445,10 @@ template<typename Cell_Type, typename Data_Type> inline std::vector< Cell_Type >
     if (check_bounds) 
         out_of_range(0u, i);
 
-    std::vector< Cell_Type > ans(nrow(), static_cast< Cell_Type >(false));
+    std::vector< Cell_Type > ans;
+    ans.reserve(nrow()); 
     for (size_t j = 0u; j < N; ++j) 
-        ans[j] = el[POS(j, i)];
+        ans.push_back(el[POS(j, i)]);
 
     return ans;
 
@@ -457,7 +465,7 @@ template<typename Cell_Type, typename Data_Type> inline void BArrayDense<Cell_Ty
         out_of_range(0u, i);
 
     for (size_t j = 0u; j < N; ++j) 
-        x->at(j) = el[POS(j, i)];//this->get_cell(iter->first, i, false);
+        x->operator[](j) = el[POS(j, i)];//this->get_cell(iter->first, i, false);
 
 }
 template<typename Cell_Type, typename Data_Type>
@@ -655,7 +663,7 @@ inline void BArrayDense<Cell_Type, Data_Type>::rm_cell (
     if (check_bounds)
         out_of_range(i,j);
 
-    BARRY_UNUSED(check_exists)
+    // BARRY_UNUSED(check_exists)
 
     // Remove the pointer first (so it wont point to empty)
     el_rowsums[i] -= el[POS(i, j)];
@@ -672,13 +680,11 @@ inline void BArrayDense<Cell_Type, Data_Type>::insert_cell (
     size_t j,
     const Cell< Cell_Type> & v,
     bool check_bounds,
-    bool check_exists
+    bool
 ) { 
 
     if (check_bounds)
         out_of_range(i,j); 
-
-    BARRY_UNUSED(check_exists)
 
     if (el[POS(i,j)] == BARRY_ZERO_DENSE)
     {
@@ -708,13 +714,11 @@ template<typename Cell_Type, typename Data_Type> inline void BArrayDense<Cell_Ty
     size_t j,
     Cell_Type v,
     bool check_bounds,
-    bool check_exists
+    bool
 ) {
 
     if (check_bounds)
         out_of_range(i,j);
-
-    BARRY_UNUSED(check_exists)
 
     if (el[POS(i,j)] == BARRY_ZERO_DENSE)
     {

diff --git a/inst/include/barry/model-meat.hpp b/inst/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/inst/include/barry/models/geese/geese-bones.hpp b/inst/include/barry/models/geese/geese-bones.hpp
@@ -41,21 +41,20 @@ inline std::vector< double > keygen_full(
     // Baseline data: nrows and columns
     std::vector< double > dat = {
         static_cast<double>(array.nrow()) * 100000 +
-         static_cast<double>(array.ncol())
+         static_cast<double>(array.ncol()),
+         1000000.0, // state of the parent
+         array.D_ptr()->duplication ? 1.0 : 0.0 // type of the parent
     };
 
     // State of the parent
-    dat.push_back(1000000.0);
-    size_t count = 0u;
+    double pow10 = 1.0;
     for (bool i : array.D_ptr()->states) {
-        dat[dat.size() - 1u] += (i ? 1.0 : 0.0) * pow(10, static_cast<double>(count));
-        count++;
+        dat[1u] += (i ? 1.0 : 0.0) * pow10;
+        pow10 *= 10.0;
     }
 
-    // Type of the parent
-    dat.push_back(array.D_ptr()->duplication ? 1.0 : 0.0);
-
     return dat;
+
 }
 
 inline bool vec_diff(

diff --git a/inst/include/barry/models/geese/geese-meat-predict.hpp b/inst/include/barry/models/geese/geese-meat-predict.hpp
@@ -170,7 +170,7 @@ inline std::vector< std::vector<double> > Geese::predict_backend(
                         // Getting the offspring state, and how it maps, only
                         // if it is not an offspring
                         const auto & off_state = array_p.get_col_vec(off);
-                        size_t loc = this->map_to_state_id[off_state];
+                        size_t loc = this->map_to_state_id.find(off_state)->second;
 
                         everything_below_p *= parent.offspring[off]->subtree_prob[loc];
 

diff --git a/inst/include/barry/models/geese/geese-types.hpp b/inst/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
diff --git a/src/Makevars b/src/Makevars
@@ -4,6 +4,7 @@
 PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS)
 PKG_CXXFLAGS= -I../inst/include/ \
 	-Dprintf_barry=Rprintf $(SHLIB_OPENMP_CXXFLAGS) \
-	-DBARRY_USER_INTERRUPT="Rcpp::checkUserInterrupt();"
+	-DBARRY_USER_INTERRUPT="Rcpp::checkUserInterrupt();" \
+	-DBARRY_USE_UNORDERED_MAP
 
 #CXXFLAGS=-O3 -Wall -pedantic