From fc700918dc7a1e043bf236e9ff73600045b2d34c Mon Sep 17 00:00:00 2001 From: dimim Date: Sat, 16 Nov 2019 15:36:00 +0100 Subject: [PATCH 1/5] Add .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..75b1567 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.swp +*.o From 0c74ea9eaa5b455e239503ebd7d1d8d5b9e9757b Mon Sep 17 00:00:00 2001 From: dimim Date: Sat, 16 Nov 2019 15:55:57 +0100 Subject: [PATCH 2/5] Add args checker Remove obsolete code/comments --- src/mmp_c.cpp | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/mmp_c.cpp b/src/mmp_c.cpp index 17d7005..58319e6 100644 --- a/src/mmp_c.cpp +++ b/src/mmp_c.cpp @@ -8,9 +8,6 @@ #include "mmp_c.h" -// [[Rcpp::plugins("cpp11")]] -// [[Rcpp::depends("RcppArmadillo")]] - #define DEBUG 0 #define db_print(...) \ do { if (DEBUG) Rprintf(__VA_ARGS__); } while (0) @@ -31,6 +28,18 @@ static double STAT_PVALUE[2]; Rcpp::List res; unsigned int kv_length; +static void check_args(const double thres, const int max_k, const std::string method) { + if (max_k < 1) { + Rcpp::stop("Invalid max_k argument provided.\nExiting...\n"); + } + else if (thres < 0 || thres >= 1) { + Rcpp::stop("Invalid thres argument provided.\nExiting...\n"); + } + else if (method.compare("pearson") && method.compare("spearman")) { + Rcpp::stop("Invalid method name provided.\nExiting...\n"); + } +} + static arma::mat calc_resid(arma::mat& cor_ds, arma::mat& sol_ds, arma::uvec& idxs_a, arma::uvec& idxs_b) { arma::mat lh_sub = cor_ds.submat(idxs_a, idxs_a); @@ -560,15 +569,8 @@ Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int max_k, hash_on ? kv_length = stats_kv[".length"] : kv_length = 0; db_print("Call: adj_med_NAs\n"); adj_med_NAs(ds); - if (max_k < 1 || thres < 0 || thres >= 1) { - db_print("True: max_k < 1 || thres < 0 || thres >= 1\n"); - if (max_k < 1) { - Rcpp::stop("Invalid max_k argument provided.\nExiting...\n"); - } - else { - Rcpp::stop("Invalud thres argument provided.\nExiting...\n"); - } - } + check_args(thres, max_k, method); + const unsigned int var_size = ds.n_cols; if (max_k > (int) var_size) { db_print("True: max_k > var_size\n"); From d91c1d6f193d0688bfe308c4637ec188e9e5fe6b Mon Sep 17 00:00:00 2001 From: dimim Date: Sun, 17 Nov 2019 01:48:52 +0100 Subject: [PATCH 3/5] Add spearman option --- src/mmp_c.cpp | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/mmp_c.cpp b/src/mmp_c.cpp index 58319e6..6f30540 100644 --- a/src/mmp_c.cpp +++ b/src/mmp_c.cpp @@ -40,6 +40,13 @@ static void check_args(const double thres, const int max_k, const std::string me } } +static void upd_args(arma::vec& target_vars, arma::mat& ds, const std::string method) { + if (!method.compare("spearman")) { + target_vars = rank_mean(target_vars, false); + ds = calc_rank(ds); + } +} + static arma::mat calc_resid(arma::mat& cor_ds, arma::mat& sol_ds, arma::uvec& idxs_a, arma::uvec& idxs_b) { arma::mat lh_sub = cor_ds.submat(idxs_a, idxs_a); @@ -380,15 +387,20 @@ static arma::vec calc_univ_pvalues(arma::vec& stats, const double dof) { } static void calc_univs(arma::vec& target_vars, arma::mat& ds, const std::string method, Rcpp::List& univs) { + arma::mat cor_mat = arma::cor(target_vars, ds); + arma::vec cor_vec = to_vec(cor_mat); + const unsigned int dof = ds.n_rows - 3; + arma::vec stats; if (!method.compare("pearson")) { db_print("True: !method.compare(\"pearson\")\n"); - arma::mat cor_mat = arma::cor(target_vars, ds); - arma::vec cor_vec = to_vec(cor_mat); - const unsigned int dof = ds.n_rows - 3; - arma::vec stats = 0.5 * arma::log((1 + cor_vec) / (1 - cor_vec)) * std::sqrt(dof); - univs["stat"] = stats; - univs["pvalue"] = calc_univ_pvalues(stats, dof); + stats = 0.5 * arma::log((1 + cor_vec) / (1 - cor_vec)) * std::sqrt(dof); } + else if (!method.compare("spearman")) { + db_print("True: !method.compare(\"spearman\")\n"); + stats = 0.5 * arma::log((1 + cor_vec) / (1 - cor_vec)) * std::sqrt(dof) / 1.029563; + } + univs["stat"] = stats; + univs["pvalue"] = calc_univ_pvalues(stats, dof); } static double get_time_spent(const clock_t begin) { @@ -569,8 +581,11 @@ Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int max_k, hash_on ? kv_length = stats_kv[".length"] : kv_length = 0; db_print("Call: adj_med_NAs\n"); adj_med_NAs(ds); + db_print("Call: check_args\n"); check_args(thres, max_k, method); + db_print("Call: upd_args\n"); + upd_args(target_vars, ds, method); const unsigned int var_size = ds.n_cols; if (max_k > (int) var_size) { db_print("True: max_k > var_size\n"); From 87d4199031efb3d8bb299efc1ea2903d88c779a0 Mon Sep 17 00:00:00 2001 From: dimim Date: Sun, 17 Nov 2019 02:09:19 +0100 Subject: [PATCH 4/5] Change *int to *int*_t Update .gitignore --- .gitignore | 1 + src/cts.cpp | 306 +++++++++++++++++++++++++------------------------- src/cts.h | 62 +++++----- src/mmp_c.cpp | 126 ++++++++++----------- src/mmp_c.h | 2 +- src/sf.cpp | 4 +- 6 files changed, 251 insertions(+), 250 deletions(-) diff --git a/.gitignore b/.gitignore index 75b1567..cb93b88 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ *.swp *.o +*.o.tmp diff --git a/src/cts.cpp b/src/cts.cpp index f35a9ec..9e8ecf0 100644 --- a/src/cts.cpp +++ b/src/cts.cpp @@ -7,8 +7,8 @@ #define db_print(...) \ do { if (DEBUG) Rprintf(__VA_ARGS__); } while (0) -static unsigned int skip_ahead(arma::uvec& vals, const unsigned int curr) { - unsigned int i; +static uint32_t skip_ahead(arma::uvec& vals, const uint32_t curr) { + uint32_t i; for (i = curr + 1; i < vals.size() && vals(i) == vals(curr); ++i) { } return i; @@ -16,14 +16,14 @@ static unsigned int skip_ahead(arma::uvec& vals, const unsigned int curr) { double calc_med(arma::vec& vals); -double calc_med(arma::vec& vals, const unsigned int size); +double calc_med(arma::vec& vals, const uint32_t size); bool adj_med_NAs(arma::mat& ds) { bool found_NA = false; - for (unsigned int i = 0; i < ds.n_rows; i++) { + for (uint32_t i = 0; i < ds.n_rows; i++) { bool found_NA_col = false; double med = 0; - for (unsigned int j = 0; j < ds.n_cols; j++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { if (!arma::is_finite(ds.at(i, j))) { if (!found_NA) { found_NA = true; } if (!found_NA_col) { @@ -39,7 +39,7 @@ bool adj_med_NAs(arma::mat& ds) { } double calc_med(arma::vec& vals) { - const unsigned int mid = vals.size() / 2 - 1; + const uint32_t mid = vals.size() / 2 - 1; double med = 0; if (vals.size() % 2 == 0) { std::nth_element(vals.begin(), vals.begin() + mid, vals.end()); @@ -56,10 +56,10 @@ double find_freq(arma::vec& vals); bool adj_freq_NAs(arma::mat& ds) { bool found_NA = false; - for (unsigned int i = 0; i < ds.n_rows; i++) { + for (uint32_t i = 0; i < ds.n_rows; i++) { bool found_NA_col = false; double freq = 0; - for (unsigned int j = 0; j < ds.n_cols; j++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { if (!arma::is_finite(ds.at(i, j))) { if (!found_NA) { found_NA = true; } if (!found_NA_col) { @@ -78,9 +78,9 @@ double find_freq(arma::vec& vals) { std::sort(vals.begin(), vals.end()); double prev = vals[0]; double freq = vals[0]; - int curr_count = 1; - int max_count = 1; - for (unsigned int i = 1; i < vals.size(); i++) { + int32_t curr_count = 1; + int32_t max_count = 1; + for (uint32_t i = 1; i < vals.size(); i++) { if (vals[i] == prev) { curr_count++; } @@ -101,7 +101,7 @@ double find_freq(arma::vec& vals) { arma::mat calc_rank(arma::mat& ds) { arma::mat rds(ds.n_rows, ds.n_cols); - for (unsigned int i = 0; i < ds.n_cols; ++i) { + for (uint32_t i = 0; i < ds.n_cols; ++i) { arma::vec curr_col = ds.col(i); rds.col(i) = rank_mean(curr_col, false); } @@ -109,16 +109,16 @@ arma::mat calc_rank(arma::mat& ds) { } arma::uvec sub_col_max_min(arma::mat& ds, const bool cont) { - int extra_val = 0; + int32_t extra_val = 0; if (!cont) { extra_val = 1; } arma::uvec max_min(ds.n_cols); - for (unsigned int j = 0; j < ds.n_cols; j++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { double max = ds(0, j); double min = ds(0, j); - for (unsigned int i = 1; i < ds.n_rows; i++) { - const unsigned int curr = ds(i, j); + for (uint32_t i = 1; i < ds.n_rows; i++) { + const uint32_t curr = ds(i, j); if (curr > max) { max = curr; } @@ -132,21 +132,21 @@ arma::uvec sub_col_max_min(arma::mat& ds, const bool cont) { } arma::mat calc_pt(arma::mat& ds, - const int df, const bool lower_tail, const bool log_p, + const int32_t df, const bool lower_tail, const bool log_p, const double add) { arma::mat pt_ds(ds.n_rows, ds.n_cols); - for (unsigned int i = 0; i < ds.n_rows; i++) { - for (unsigned int j = 0; j < ds.n_cols; j++) { + for (uint32_t i = 0; i < ds.n_rows; i++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { pt_ds(i, j) = add + R::pt(ds(i, j), df, lower_tail, log_p); } } return pt_ds; } -std::vector det_cols(arma::umat& ds, const unsigned int val) { - std::vector cols; - for (unsigned int i = 0; i < ds.n_rows; ++i) { - for (unsigned int j = 0; j < ds.n_cols; ++j) { +std::vector det_cols(arma::umat& ds, const uint32_t val) { + std::vector cols; + for (uint32_t i = 0; i < ds.n_rows; ++i) { + for (uint32_t j = 0; j < ds.n_cols; ++j) { if (ds(i, j) == val) { cols.push_back(j); } @@ -156,39 +156,39 @@ std::vector det_cols(arma::umat& ds, const unsigned int val) { return cols; } -arma::mat ext_cols(arma::mat& ds, const unsigned int col_a, - const unsigned int col_b) { - const unsigned int ncols = 2; +arma::mat ext_cols(arma::mat& ds, const uint32_t col_a, + const uint32_t col_b) { + const uint32_t ncols = 2; arma::mat ext_ds(ds.n_rows, ncols); - const unsigned int ext_col_a = 0; - const unsigned int ext_col_b = 1; - for (unsigned int i = 0; i < ds.n_rows; i++) { + const uint32_t ext_col_a = 0; + const uint32_t ext_col_b = 1; + for (uint32_t i = 0; i < ds.n_rows; i++) { ext_ds(i, ext_col_a) = ds(i, col_a); ext_ds(i, ext_col_b) = ds(i, col_b); } return ext_ds; } -void cp_lt(arma::mat& src, arma::mat& dst, const int val) { - for (unsigned int i = 0; i < src.n_rows; i++) { - for (unsigned int j = 0; j < src.n_cols; j++) { +void cp_lt(arma::mat& src, arma::mat& dst, const int32_t val) { + for (uint32_t i = 0; i < src.n_rows; i++) { + for (uint32_t j = 0; j < src.n_cols; j++) { i > j ? dst(i, j) = val : dst(i, j) = src(i, j); } } } void adj_diag(arma::mat& ds, const double val) { - for (unsigned int i = 0; i < ds.n_rows && i < ds.n_cols; i++) { + for (uint32_t i = 0; i < ds.n_rows && i < ds.n_cols; i++) { ds(i, i) = val; } } arma::mat cbind_tran_mat(arma::mat& ds1, arma::mat& ds2) { - const unsigned int nrows = ds1.n_cols; - const unsigned int ncols = ds1.n_rows + ds2.n_rows; + const uint32_t nrows = ds1.n_cols; + const uint32_t ncols = ds1.n_rows + ds2.n_rows; arma::mat ds(nrows, ncols); - for (unsigned int i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { - for (unsigned int j = 0; j < ds1.n_cols && j < ds2.n_cols; j++) { + for (uint32_t i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { + for (uint32_t j = 0; j < ds1.n_cols && j < ds2.n_cols; j++) { ds(j, i) = ds1(i, j); ds(j, i + ds1.n_rows) = ds2(i, j); } @@ -196,23 +196,23 @@ arma::mat cbind_tran_mat(arma::mat& ds1, arma::mat& ds2) { return ds; } -bool is_dupl_row(arma::mat& ds, const unsigned int lrow); +bool is_dupl_row(arma::mat& ds, const uint32_t lrow); -std::vector get_dupl_rows_pos(arma::mat& ds); +std::vector get_dupl_rows_pos(arma::mat& ds); arma::mat rm_dupl_rows(arma::mat& src); arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, const bool ass1, const bool ass2) { - const unsigned int nrows = ds1.n_rows + ds2.n_rows; - unsigned int ncols; + const uint32_t nrows = ds1.n_rows + ds2.n_rows; + uint32_t ncols; ds1.n_cols > ds2.n_cols ? ncols = ds1.n_cols : ncols = ds2.n_cols; arma::mat ds(nrows, ncols, arma::fill::zeros); - unsigned int row = 0; - unsigned int col = 0; + uint32_t row = 0; + uint32_t col = 0; if (ass1) { - for (unsigned int i = 0; i < ds1.n_rows; i++) { - for (unsigned int j = 0; j < ds1.n_cols; j++) { + for (uint32_t i = 0; i < ds1.n_rows; i++) { + for (uint32_t j = 0; j < ds1.n_cols; j++) { ds(row, col) = ds1(i, j); col++; } @@ -224,8 +224,8 @@ arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, row = ds1.n_rows; } if (ass2) { - for (unsigned int i = 0; i < ds2.n_rows; i++) { - for (unsigned int j = 0; j < ds2.n_cols; j++) { + for (uint32_t i = 0; i < ds2.n_rows; i++) { + for (uint32_t j = 0; j < ds2.n_cols; j++) { ds(row, col) = ds2(i, j); col++; } @@ -237,29 +237,29 @@ arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, } arma::mat rm_dupl_rows(arma::mat& src) { - std::vector dupls_pos = get_dupl_rows_pos(src); + std::vector dupls_pos = get_dupl_rows_pos(src); if (dupls_pos.empty()) { return src; } - const unsigned int nrows = src.n_rows - dupls_pos.size(); + const uint32_t nrows = src.n_rows - dupls_pos.size(); arma::mat dst(nrows, src.n_cols); - unsigned int dupls_pos_cntr = 0; - for (unsigned int i = 0, src_row = 0; i < dst.n_rows; i++, src_row++) { + uint32_t dupls_pos_cntr = 0; + for (uint32_t i = 0, src_row = 0; i < dst.n_rows; i++, src_row++) { while (dupls_pos_cntr < dupls_pos.size() && dupls_pos.at(dupls_pos_cntr) == src_row) { src_row++; dupls_pos_cntr++; } - for (unsigned int j = 0; j < dst.n_cols; j++) { + for (uint32_t j = 0; j < dst.n_cols; j++) { dst(i, j) = src(src_row, j); } } return dst; } -std::vector get_dupl_rows_pos(arma::mat& ds) { - std::vector dupls_pos; - for (unsigned int i = 1; i < ds.n_rows; i++) { +std::vector get_dupl_rows_pos(arma::mat& ds) { + std::vector dupls_pos; + for (uint32_t i = 1; i < ds.n_rows; i++) { if (is_dupl_row(ds, i)) { dupls_pos.push_back(i); } @@ -267,9 +267,9 @@ std::vector get_dupl_rows_pos(arma::mat& ds) { return dupls_pos; } -bool is_dupl_row(arma::mat& ds, const unsigned int lrow) { - for (unsigned int i = 0; i < lrow; i++) { - for (unsigned int j = 0; +bool is_dupl_row(arma::mat& ds, const uint32_t lrow) { + for (uint32_t i = 0; i < lrow; i++) { + for (uint32_t j = 0; ; j++) { if (ds(i, j) != ds(lrow, j)) { @@ -285,17 +285,17 @@ bool is_dupl_row(arma::mat& ds, const unsigned int lrow) { arma::vec to_vec(arma::mat& ds) { arma::vec vals(ds.n_rows * ds.n_cols); - unsigned int k = 0; - for (unsigned int i = 0; i < ds.n_rows; i++) { - for (unsigned int j = 0; j < ds.n_cols; j++) { + uint32_t k = 0; + for (uint32_t i = 0; i < ds.n_rows; i++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { vals(k++) = ds(i, j); } } return vals; } -bool found_match(const unsigned int x, arma::uvec& vals) { - for (unsigned int i = 0; i < vals.size(); ++i) { +bool found_match(const uint32_t x, arma::uvec& vals) { + for (uint32_t i = 0; i < vals.size(); ++i) { if (x == vals[i]) { return true; } @@ -309,11 +309,11 @@ bool are_equal(arma::mat& ds, arma::vec& vals, (!by_col_idx && (ds.n_rows * ds.n_cols) != vals.size())) { return false; } - unsigned int j = 0; - unsigned int k = 0; + uint32_t j = 0; + uint32_t k = 0; by_col_idx ? j = col : j = 0; for (; (!by_col_idx && j < ds.n_cols) || (by_col_idx && col == j); ++j) { - for (unsigned int i = 0; i < ds.n_rows; ++i) { + for (uint32_t i = 0; i < ds.n_rows; ++i) { if (ds(i, j) != vals[k++]) { return false; } @@ -322,9 +322,9 @@ bool are_equal(arma::mat& ds, arma::vec& vals, return true; } -std::vector rm_lt_nan(arma::uvec& idxs, const unsigned int limit) { - std::vector idxs_adj; - for (unsigned int i = 0; i < idxs.size(); ++i) { +std::vector rm_lt_nan(arma::uvec& idxs, const uint32_t limit) { + std::vector idxs_adj; + for (uint32_t i = 0; i < idxs.size(); ++i) { if (arma::is_finite(idxs[i]) || idxs[i] >= limit) { idxs_adj.push_back(idxs[i]); } @@ -332,9 +332,9 @@ std::vector rm_lt_nan(arma::uvec& idxs, const unsigned int limit) return idxs_adj; } -std::unordered_set get_diff(std::vector& lh, const unsigned int val) { - std::unordered_set diff; - for (unsigned int i = 0; i < lh.size(); ++i) { +std::unordered_set get_diff(std::vector& lh, const uint32_t val) { + std::unordered_set diff; + for (uint32_t i = 0; i < lh.size(); ++i) { if (lh[i] != val) { diff.insert(lh[i]); } @@ -342,24 +342,24 @@ std::unordered_set get_diff(std::vector& lh, const u return diff; } -arma::umat nchoosek(std::vector& idxs, const unsigned int k) { +arma::umat nchoosek(std::vector& idxs, const uint32_t k) { if (idxs.size() == 1) { arma::umat ret(1, 1); ret(0, 0) = R::choose(idxs[0], k); return ret; } - return find_combn>(idxs, k); + return find_combn>(idxs, k); } std::vector inter_helper(arma::vec& vals1, arma::vec& vals2) { std::sort(vals1.begin(), vals1.end()); std::sort(vals2.begin(), vals2.end()); std::vector inter_vals; - for (unsigned int i = 0, j = 0; i < vals1.size() && j < vals2.size();) { + for (uint32_t i = 0, j = 0; i < vals1.size() && j < vals2.size();) { const double curr1 = vals1(i); const double curr2 = vals2(j); if (curr1 == curr2) { - const unsigned int size = inter_vals.size(); + const uint32_t size = inter_vals.size(); if (!size || (inter_vals.back() != curr1)) { inter_vals.push_back(curr1); } @@ -383,11 +383,11 @@ std::vector inter_helper(arma::vec& vals1, arma::vec& vals2) { } // Alters row_idxs -void append_rows(arma::mat& ds, const double val, std::vector& row_idxs); +void append_rows(arma::mat& ds, const double val, std::vector& row_idxs); -std::vector index_row_eq(arma::mat& ds, std::vector& vals) { - std::vector row_idxs; - for (unsigned int i = 0; i < vals.size(); i++) { +std::vector index_row_eq(arma::mat& ds, std::vector& vals) { + std::vector row_idxs; + for (uint32_t i = 0; i < vals.size(); i++) { append_rows(ds, vals.at(i), row_idxs); } std::sort(row_idxs.begin(), row_idxs.end()); @@ -395,9 +395,9 @@ std::vector index_row_eq(arma::mat& ds, std::vector& vals) return row_idxs; } -void append_rows(arma::mat& ds, const double val, std::vector& row_idxs) { - for (unsigned int i = 0; i < ds.n_rows; i++) { - for (unsigned int j = 0; j < ds.n_cols; j++) { +void append_rows(arma::mat& ds, const double val, std::vector& row_idxs) { + for (uint32_t i = 0; i < ds.n_rows; i++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { if (ds(i, j) == val) { row_idxs.push_back(i); } @@ -406,17 +406,17 @@ void append_rows(arma::mat& ds, const double val, std::vector& row } arma::mat rm_rows(arma::mat& src, arma::uvec& rows) { - unsigned int dst_nrows = src.n_rows - rows.size(); - unsigned int dst_ncols = src.n_cols; + uint32_t dst_nrows = src.n_rows - rows.size(); + uint32_t dst_ncols = src.n_cols; arma::mat dst(dst_nrows, dst_ncols); - unsigned int src_row = 0; - unsigned int rows_idx = 0; - for (unsigned int dst_row = 0; dst_row < dst_nrows; dst_row++) { + uint32_t src_row = 0; + uint32_t rows_idx = 0; + for (uint32_t dst_row = 0; dst_row < dst_nrows; dst_row++) { while (src_row < src.n_rows && rows_idx < rows.size() && src_row == rows(rows_idx)) { src_row++; rows_idx = skip_ahead(rows, rows_idx); } - for (unsigned int col = 0; col < dst_ncols; col++) { + for (uint32_t col = 0; col < dst_ncols; col++) { dst(dst_row, col) = src(src_row, col); } src_row++; @@ -424,20 +424,20 @@ arma::mat rm_rows(arma::mat& src, arma::uvec& rows) { return dst; } -unsigned int skip_ahead_std(std::vector rows, const unsigned int curr); +uint32_t skip_ahead_std(std::vector rows, const uint32_t curr); -arma::mat rm_rows_std(arma::mat& src, std::vector& rows) { - unsigned int dst_nrows = src.n_rows - rows.size(); - unsigned int dst_ncols = src.n_cols; +arma::mat rm_rows_std(arma::mat& src, std::vector& rows) { + uint32_t dst_nrows = src.n_rows - rows.size(); + uint32_t dst_ncols = src.n_cols; arma::mat dst(dst_nrows, dst_ncols); - unsigned int src_row = 0; - unsigned int rows_idx = 0; - for (unsigned int dst_row = 0; dst_row < dst_nrows; dst_row++) { + uint32_t src_row = 0; + uint32_t rows_idx = 0; + for (uint32_t dst_row = 0; dst_row < dst_nrows; dst_row++) { while (src_row < src.n_rows && rows_idx < rows.size() && src_row == rows.at(rows_idx)) { src_row++; rows_idx = skip_ahead_std(rows, rows_idx); } - for (unsigned int col = 0; col < dst_ncols; col++) { + for (uint32_t col = 0; col < dst_ncols; col++) { dst(dst_row, col) = src(src_row, col); } src_row++; @@ -445,27 +445,27 @@ arma::mat rm_rows_std(arma::mat& src, std::vector& rows) { return dst; } -unsigned int skip_ahead_std(std::vector rows, const unsigned int curr) { - unsigned int i; +uint32_t skip_ahead_std(std::vector rows, const uint32_t curr) { + uint32_t i; for (i = curr + 1; i < rows.size() && rows.at(i) == rows.at(curr); i++) { } return i; } // No sorting -arma::umat rm_cols(arma::umat src, std::vector cols) { - unsigned int dst_nrows = src.n_rows; - unsigned int dst_ncols = src.n_cols - cols.size(); +arma::umat rm_cols(arma::umat src, std::vector cols) { + uint32_t dst_nrows = src.n_rows; + uint32_t dst_ncols = src.n_cols - cols.size(); arma::umat dst(dst_nrows, dst_ncols); - unsigned int src_col = 0; - unsigned int cols_idx = 0; - for (unsigned int dst_col = 0; dst_col < dst_ncols; ++dst_col) { + uint32_t src_col = 0; + uint32_t cols_idx = 0; + for (uint32_t dst_col = 0; dst_col < dst_ncols; ++dst_col) { while (src_col < src.n_cols && cols_idx < cols.size() && src_col == cols.at(cols_idx)) { ++src_col; ++cols_idx; } - for (unsigned int row = 0; row < dst_nrows; ++row) { + for (uint32_t row = 0; row < dst_nrows; ++row) { dst(row, dst_col) = src(row, src_col); } ++src_col; @@ -473,12 +473,12 @@ arma::umat rm_cols(arma::umat src, std::vector cols) { return dst; } -arma::mat order_col(arma::mat& ds, const unsigned int col) { +arma::mat order_col(arma::mat& ds, const uint32_t col) { arma::mat ordered_ds(ds.n_rows, ds.n_cols); arma::uvec order_idxs = arma::sort_index(ds.col(col)); - for (unsigned int i = 0; i < ds.n_rows; i++) { - const unsigned int idx = order_idxs(i); - for (unsigned int j = 0; j < ds.n_cols; j++) { + for (uint32_t i = 0; i < ds.n_rows; i++) { + const uint32_t idx = order_idxs(i); + for (uint32_t j = 0; j < ds.n_cols; j++) { ordered_ds(i, j) = ds(idx, j); } } @@ -487,8 +487,8 @@ arma::mat order_col(arma::mat& ds, const unsigned int col) { arma::mat form_cmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols) { arma::mat formed_ds(cols.size(), rows.size()); - for (unsigned int i = 0; i < rows.size(); i++) { - for (unsigned int j = 0; j < cols.size(); j++) { + for (uint32_t i = 0; i < rows.size(); i++) { + for (uint32_t j = 0; j < cols.size(); j++) { formed_ds(j, i) = ds(rows(i), cols(j)); } } @@ -497,19 +497,19 @@ arma::mat form_cmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols) { arma::mat form_rmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols) { arma::mat formed_ds(rows.size(), cols.size()); - for (unsigned int i = 0; i < rows.size(); i++) { - for (unsigned int j = 0; j < cols.size(); j++) { + for (uint32_t i = 0; i < rows.size(); i++) { + for (uint32_t j = 0; j < cols.size(); j++) { formed_ds(i, j) = ds(rows(i), cols(j)); } } return formed_ds; } -arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, - std::vector& cols) { +arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, + std::vector& cols) { arma::mat formed_ds(rows.size(), cols.size()); - for (unsigned int i = 0; i < rows.size(); i++) { - for (unsigned int j = 0; j < cols.size(); j++) { + for (uint32_t i = 0; i < rows.size(); i++) { + for (uint32_t j = 0; j < cols.size(); j++) { formed_ds(i, j) = ds(rows.at(i), cols.at(j)); } } @@ -518,8 +518,8 @@ arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, arma::mat merge_cols(arma::mat& ds, arma::uvec& idxs) { arma::mat merged_ds(ds.n_rows, idxs.size()); - for (unsigned int i = 0; i < idxs.size(); i++) { - for (unsigned int j = 0; j < ds.n_rows; j++) { + for (uint32_t i = 0; i < idxs.size(); i++) { + for (uint32_t j = 0; j < ds.n_rows; j++) { merged_ds(j, i) = ds(j, idxs(i)); } } @@ -527,12 +527,12 @@ arma::mat merge_cols(arma::mat& ds, arma::uvec& idxs) { } arma::mat form_ncolcmat(arma::vec& vals, arma::mat& ds) { - const unsigned int nrows = vals.size(); - const unsigned int ncols = 1 + ds.n_cols; + const uint32_t nrows = vals.size(); + const uint32_t ncols = 1 + ds.n_cols; arma::mat formed_ds(nrows, ncols); - for (unsigned int i = 0; i < nrows; i++) { + for (uint32_t i = 0; i < nrows; i++) { formed_ds(i, 0) = vals(i); - for (unsigned int j = 1; j < ncols; j++) { + for (uint32_t j = 1; j < ncols; j++) { formed_ds(i, j) = ds(i, (j - 1)); } } @@ -540,49 +540,49 @@ arma::mat form_ncolcmat(arma::vec& vals, arma::mat& ds) { } arma::mat form_c2mat(arma::vec& vals1, arma::vec& vals2) { - const unsigned int nrows = vals1.size(); - const unsigned int ncols = 2; + const uint32_t nrows = vals1.size(); + const uint32_t ncols = 2; arma::mat formed_ds(nrows, ncols); - for (unsigned int i = 0; i < nrows; i++) { + for (uint32_t i = 0; i < nrows; i++) { formed_ds(i, 0) = vals1(i); formed_ds(i, 1) = vals2(i); } return formed_ds; } -arma::uvec form_vec(arma::mat& ds, const unsigned int row, arma::uvec& cols) { +arma::uvec form_vec(arma::mat& ds, const uint32_t row, arma::uvec& cols) { arma::uvec vals(cols.size()); - for (unsigned int i = 0; i < cols.size(); i++) { + for (uint32_t i = 0; i < cols.size(); i++) { vals(i) = ds(row, cols(i)); } return vals; } -arma::vec form_vec_wvals(arma::mat& ds, const unsigned int row, arma::uvec& cols, +arma::vec form_vec_wvals(arma::mat& ds, const uint32_t row, arma::uvec& cols, arma::vec& vals) { arma::vec formed_vec(cols.size() + vals.size()); - unsigned int i; + uint32_t i; for (i = 0; i < cols.size(); i++) { formed_vec(i) = ds(row, cols(i)); } - for (unsigned int j = 0; i < formed_vec.size(); i++, j++) { + for (uint32_t j = 0; i < formed_vec.size(); i++, j++) { formed_vec(i) = vals(j); } return formed_vec; } -arma::mat append_row(arma::mat& ds, const unsigned int row, arma::vec& vals) { - for (unsigned int j = 0; j < ds.n_cols; j++) { +arma::mat append_row(arma::mat& ds, const uint32_t row, arma::vec& vals) { + for (uint32_t j = 0; j < ds.n_cols; j++) { ds(row, j) = vals(j); } return ds; } -std::vector rsum_gt_zero_idxs(arma::mat& ds) { - std::vector idxs; - for (unsigned int i = 0; i < ds.n_rows; i++) { +std::vector rsum_gt_zero_idxs(arma::mat& ds) { + std::vector idxs; + for (uint32_t i = 0; i < ds.n_rows; i++) { double sum = 0; - for (unsigned int j = 0; j < ds.n_cols; j++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { sum += ds(i, j); } if (sum > 0) { @@ -594,40 +594,40 @@ std::vector rsum_gt_zero_idxs(arma::mat& ds) { arma::vec form_cmat_vec(arma::mat& ds, arma::rowvec& vals) { arma::vec concat_vals(ds.n_rows * ds.n_cols + vals.size()); - unsigned int cd_cntr = 0; - for (unsigned int j = 0; j < ds.n_cols; j++) { - for (unsigned int i = 0; i < ds.n_rows; i++) { + uint32_t cd_cntr = 0; + for (uint32_t j = 0; j < ds.n_cols; j++) { + for (uint32_t i = 0; i < ds.n_rows; i++) { concat_vals(cd_cntr++) = ds(i, j); } } - for (unsigned int i = 0; i < vals.size(); i++) { + for (uint32_t i = 0; i < vals.size(); i++) { concat_vals[cd_cntr++] = vals(i); } return concat_vals; } arma::mat cbind_mat(arma::mat& ds1, arma::mat& ds2) { - const unsigned int nrows = ds1.n_rows; - const unsigned int ncols = ds1.n_cols + ds2.n_cols; + const uint32_t nrows = ds1.n_rows; + const uint32_t ncols = ds1.n_cols + ds2.n_cols; arma::mat ds(nrows, ncols); - for (unsigned int i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { - for (unsigned int j = 0; j < ds1.n_cols; j++) { + for (uint32_t i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { + for (uint32_t j = 0; j < ds1.n_cols; j++) { ds(i, j) = ds1(i, j); } - for (unsigned int j = 0; j < ds2.n_cols; j++) { + for (uint32_t j = 0; j < ds2.n_cols; j++) { ds(i, j + ds1.n_cols) = ds2(i, j); } } return ds; } -arma::mat adj_cols(arma::mat& src, const unsigned int dst_ncols) { - const unsigned int dst_nrows = src.n_rows * src.n_cols / dst_ncols; +arma::mat adj_cols(arma::mat& src, const uint32_t dst_ncols) { + const uint32_t dst_nrows = src.n_rows * src.n_cols / dst_ncols; arma::mat dst(dst_nrows, dst_ncols); - unsigned int src_row = 0; - unsigned int src_col = 0; - unsigned int dst_row = 0; - unsigned int dst_col = 0; + uint32_t src_row = 0; + uint32_t src_col = 0; + uint32_t dst_row = 0; + uint32_t dst_col = 0; while (src_col < src.n_cols && dst_col < dst_ncols) { while (src_row < src.n_rows && dst_row < dst_nrows) { dst(dst_row++, dst_col) = src(src_row++, src_col); diff --git a/src/cts.h b/src/cts.h index 564fc75..5a7ef11 100644 --- a/src/cts.h +++ b/src/cts.h @@ -13,32 +13,32 @@ // [[Rcpp::depends("RcppArmadillo")]] template -static void combn(T2& vals, const int n, const unsigned int start_idx, - std::vector& combn_data, T1& combn_ds, unsigned int& combn_col); +static void combn(T2& vals, const int32_t n, const uint32_t start_idx, + std::vector& combn_data, T1& combn_ds, uint32_t& combn_col); template -T1 find_combn(T2& vals, const int n) { - static unsigned int combn_col = 0; - const unsigned int nrows = n; - const unsigned int ncols = std::round(R::choose(vals.size(), n)); +T1 find_combn(T2& vals, const int32_t n) { + static uint32_t combn_col = 0; + const uint32_t nrows = n; + const uint32_t ncols = std::round(R::choose(vals.size(), n)); T1 combn_ds(nrows, ncols); std::vector combn_data(nrows); - const unsigned int start_idx = 0; + const uint32_t start_idx = 0; combn_col = 0; combn(vals, n, start_idx, combn_data, combn_ds, combn_col); return combn_ds; } template -static void combn(T2& vals, const int n, const unsigned int start_idx, - std::vector& combn_data, T1& combn_ds, unsigned int& combn_col) { +static void combn(T2& vals, const int32_t n, const uint32_t start_idx, + std::vector& combn_data, T1& combn_ds, uint32_t& combn_col) { if (!n) { - for (unsigned int i = 0; i < combn_ds.n_rows && combn_col < combn_ds.n_cols; i++) { + for (uint32_t i = 0; i < combn_ds.n_rows && combn_col < combn_ds.n_cols; i++) { combn_ds(i, combn_col) = combn_data.at(i); } combn_col++; return; } - for (unsigned int i = start_idx; i <= (vals.size() - n); i++) { + for (uint32_t i = start_idx; i <= (vals.size() - n); i++) { combn_data.at(combn_ds.n_rows - n) = vals[i]; combn(vals, n - 1, i + 1, combn_data, combn_ds, combn_col); } @@ -55,16 +55,16 @@ arma::mat calc_rank(arma::mat& ds); arma::uvec sub_col_max_min(arma::mat& ds, const bool cont); arma::mat calc_pt(arma::mat& ds, - const int df, const bool lower_tail, const bool log_p, + const int32_t df, const bool lower_tail, const bool log_p, const double add); -std::vector det_cols(arma::umat& ds, const unsigned int val); +std::vector det_cols(arma::umat& ds, const uint32_t val); -arma::mat ext_cols(arma::mat& ds, const unsigned int col_a, - const unsigned int col_b); +arma::mat ext_cols(arma::mat& ds, const uint32_t col_a, + const uint32_t col_b); // Alters dst -void cp_lt(arma::mat& src, arma::mat& dst, const int val); +void cp_lt(arma::mat& src, arma::mat& dst, const int32_t val); // Alters ds void adj_diag(arma::mat& ds, const double val); @@ -76,35 +76,35 @@ arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, arma::vec to_vec(arma::mat& ds); -bool found_match(const unsigned int x, arma::uvec& vals); +bool found_match(const uint32_t x, arma::uvec& vals); bool are_equal(arma::mat& ds, arma::vec& vals, const bool by_col_idx = false, const unsigned col = 0); -std::vector rm_lt_nan(arma::uvec& idxs, const unsigned int limit); +std::vector rm_lt_nan(arma::uvec& idxs, const uint32_t limit); -std::unordered_set get_diff(std::vector& lh, const unsigned int val); +std::unordered_set get_diff(std::vector& lh, const uint32_t val); -arma::umat nchoosek(std::vector& idxs, const unsigned int k); +arma::umat nchoosek(std::vector& idxs, const uint32_t k); std::vector inter_helper(arma::vec& vals1, arma::vec& vals2); -std::vector index_row_eq(arma::mat& ds, std::vector& vals); +std::vector index_row_eq(arma::mat& ds, std::vector& vals); arma::mat rm_rows(arma::mat& src, arma::uvec& rows); -arma::mat rm_rows_std(arma::mat& src, std::vector& rows); +arma::mat rm_rows_std(arma::mat& src, std::vector& rows); -arma::umat rm_cols(arma::umat src, std::vector cols); +arma::umat rm_cols(arma::umat src, std::vector cols); -arma::mat order_col(arma::mat& ds, const unsigned int col); +arma::mat order_col(arma::mat& ds, const uint32_t col); arma::mat form_cmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols); arma::mat form_rmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols); -arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, - std::vector& cols); +arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, + std::vector& cols); arma::mat merge_cols(arma::mat& ds, arma::uvec& idxs); @@ -112,20 +112,20 @@ arma::mat form_ncolcmat(arma::vec& vals, arma::mat& ds); arma::mat form_c2mat(arma::vec& vals1, arma::vec& vals2); -arma::uvec form_vec(arma::mat& ds, const unsigned int row, arma::uvec& cols); +arma::uvec form_vec(arma::mat& ds, const uint32_t row, arma::uvec& cols); -arma::vec form_vec_wvals(arma::mat& ds, const unsigned int row, arma::uvec& cols, +arma::vec form_vec_wvals(arma::mat& ds, const uint32_t row, arma::uvec& cols, arma::vec& vals); // Alters ds -arma::mat append_row(arma::mat& ds, const unsigned int row, arma::vec& vals); +arma::mat append_row(arma::mat& ds, const uint32_t row, arma::vec& vals); -std::vector rsum_gt_zero_idxs(arma::mat& ds); +std::vector rsum_gt_zero_idxs(arma::mat& ds); arma::vec form_cmat_vec(arma::mat& ds, arma::rowvec& vals); arma::mat cbind_mat(arma::mat& ds1, arma::mat& ds2); -arma::mat adj_cols(arma::mat& src, const unsigned int dst_ncols); +arma::mat adj_cols(arma::mat& src, const uint32_t dst_ncols); #endif diff --git a/src/mmp_c.cpp b/src/mmp_c.cpp index 6f30540..1477c62 100644 --- a/src/mmp_c.cpp +++ b/src/mmp_c.cpp @@ -21,14 +21,14 @@ static arma::uvec DEF_UVEC; static Rcpp::List DEF_LIST; static Rcpp::Environment DEF_ENV; -static const unsigned int STAT_POS = 0; -static const unsigned int PVALUE_POS = 1; +static const uint32_t STAT_POS = 0; +static const uint32_t PVALUE_POS = 1; static double STAT_PVALUE[2]; Rcpp::List res; -unsigned int kv_length; +uint32_t kv_length; -static void check_args(const double thres, const int max_k, const std::string method) { +static void check_args(const double thres, const int32_t max_k, const std::string method) { if (max_k < 1) { Rcpp::stop("Invalid max_k argument provided.\nExiting...\n"); } @@ -63,22 +63,22 @@ static arma::mat calc_sol(arma::mat& ds, arma::uvec& idxs_a, arma::uvec& idxs_b) return res; } -static void upd_col(arma::vec& src, arma::mat& dst, const unsigned int col) { - for (unsigned int i = 0; i < src.size(); ++i) { +static void upd_col(arma::vec& src, arma::mat& dst, const uint32_t col) { + for (uint32_t i = 0; i < src.size(); ++i) { dst(i, col) = src[i]; } } static arma::mat col_bind(arma::vec& src_a, arma::vec& src_b, arma::mat& src_c) { - const unsigned int nrows = src_c.n_rows; - const unsigned int ncols = src_c.n_cols + 2; + const uint32_t nrows = src_c.n_rows; + const uint32_t ncols = src_c.n_cols + 2; arma::mat dst(nrows, ncols); upd_col(src_a, dst, 0); upd_col(src_b, dst, 1); - for (unsigned int src_col = 0, dst_col = 2; + for (uint32_t src_col = 0, dst_col = 2; src_col < src_c.n_cols && dst_col < dst.n_cols; ++src_col, ++dst_col) { - for (unsigned int i = 0; i < nrows; ++i) { + for (uint32_t i = 0; i < nrows; ++i) { dst(i, dst_col) = src_c(i, src_col); } } @@ -92,26 +92,26 @@ static void form_ret(const double stat, const double pvalue) { } -static std::string form_key(const unsigned int val, std::vector& vals) { +static std::string form_key(const uint32_t val, std::vector& vals) { std::string key = std::to_string(val + 1); - for (unsigned int i = 0; i < vals.size(); ++i) { + for (uint32_t i = 0; i < vals.size(); ++i) { key += " " + std::to_string(vals[i] + 1); } return key; } void calc_pearson(arma::vec& target_vars, arma::mat& ds, - const unsigned int x, arma::uvec& idxs, Rcpp::List& univs, + const uint32_t x, arma::uvec& idxs, Rcpp::List& univs, const bool hash_on, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { const double def_stat = 0; const double def_pvalue = std::log(1); db_print("Call: rm_lt_nan\n"); - std::vector idxs_adj = rm_lt_nan(idxs, 0); + std::vector idxs_adj = rm_lt_nan(idxs, 0); std::string key; if (hash_on) { db_print("True: hash_on\n"); - std::vector idxs_nz(idxs_adj); + std::vector idxs_nz(idxs_adj); db_print("Call: std::sort\n"); std::sort(idxs_nz.begin(), idxs_nz.end()); db_print("Call: form_key\n"); @@ -161,7 +161,7 @@ void calc_pearson(arma::vec& target_vars, arma::mat& ds, else { db_print("True: tmp_ds.n_cols != 1\n"); db_print("Checking tmp_ds, tmp_vec by column equality.\n"); - for (unsigned int j = 0; j < tmp_ds.n_cols; ++j) { + for (uint32_t j = 0; j < tmp_ds.n_cols; ++j) { if (are_equal(tmp_ds, tmp_vec, true, j)) { if (hash_on) { if (!stats_kv.exists(key)) { @@ -201,7 +201,7 @@ void calc_pearson(arma::vec& target_vars, arma::mat& ds, arma::uvec idxs_a(2); std::iota(idxs_a.begin(), idxs_a.end(), 0); arma::uvec idxs_b(idxs.size()); - for (unsigned int i = 0; i < idxs_b.size(); ++i) { + for (uint32_t i = 0; i < idxs_b.size(); ++i) { idxs_b[i] = i + 2; } db_print("Call: calc_sol\n"); @@ -241,16 +241,16 @@ bool cmp_pvalues(const double stat_lh, const double stat_rh, } // Alters ds -static void rbind_mat(arma::umat& ds, const unsigned int val) { +static void rbind_mat(arma::umat& ds, const uint32_t val) { ds.resize(ds.n_rows + 1, ds.n_cols); - for (unsigned int j = 0; j < ds.n_cols; ++j) { + for (uint32_t j = 0; j < ds.n_cols; ++j) { ds(ds.n_rows - 1, j) = val; } } -static std::vector keep_val(arma::uvec& vals, const unsigned int val_to_keep) { - std::vector vals_adj; - for (unsigned int i = 0; i < vals.size(); ++i) { +static std::vector keep_val(arma::uvec& vals, const uint32_t val_to_keep) { + std::vector vals_adj; + for (uint32_t i = 0; i < vals.size(); ++i) { if (vals[i] == val_to_keep) { vals_adj.push_back(i); } @@ -259,18 +259,18 @@ static std::vector keep_val(arma::uvec& vals, const unsigned int v } void assoc_min(arma::vec& target_vars, arma::mat& ds, const std::string method, - const int max_k, const unsigned int sel_idx, arma::uvec& sel_idxs, arma::vec& stats, + const int32_t max_k, const uint32_t sel_idx, arma::uvec& sel_idxs, arma::vec& stats, arma::vec& pvalues, arma::uvec& sel_ord_idxs, const bool hash_on, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { double sel_stat = stats[sel_idx]; double sel_pvalue = pvalues[sel_idx]; - std::vector sel_idxs_adj = keep_val(sel_idxs, 1); - const unsigned int k = std::min(max_k, (int) sel_idxs_adj.size()); - unsigned int curr_k = 1; + std::vector sel_idxs_adj = keep_val(sel_idxs, 1); + const uint32_t k = std::min(max_k, (int) sel_idxs_adj.size()); + uint32_t curr_k = 1; while (curr_k <= k) { - const unsigned int max_idx = arma::index_max(sel_ord_idxs); - std::unordered_set tmp_idxs = get_diff(sel_idxs_adj, max_idx); - std::vector tmp_idxs_adj(tmp_idxs.begin(), tmp_idxs.end()); + const uint32_t max_idx = arma::index_max(sel_ord_idxs); + std::unordered_set tmp_idxs = get_diff(sel_idxs_adj, max_idx); + std::vector tmp_idxs_adj(tmp_idxs.begin(), tmp_idxs.end()); arma::umat subsetcsk; if (curr_k == 1) { subsetcsk = arma::umat(1, 1); @@ -280,7 +280,7 @@ void assoc_min(arma::vec& target_vars, arma::mat& ds, const std::string method, subsetcsk = nchoosek(tmp_idxs_adj, curr_k - 1); rbind_mat(subsetcsk, max_idx); } - for (unsigned int i = 0; i < subsetcsk.n_cols; ++i) { + for (uint32_t i = 0; i < subsetcsk.n_cols; ++i) { arma::uvec subsetcsk_col = subsetcsk.col(i); Rcpp::List univs; calc_pearson(target_vars, ds, sel_idx, subsetcsk_col, univs, hash_on, stats_kv, pvalues_kv); @@ -296,15 +296,15 @@ void assoc_min(arma::vec& target_vars, arma::mat& ds, const std::string method, form_ret(sel_stat, sel_pvalue); } -unsigned int assoc_max_min(arma::vec& target_vars, arma::mat& ds, const std::string method, - const double thres, const int max_k, arma::uvec& sel_idxs, +uint32_t assoc_max_min(arma::vec& target_vars, arma::mat& ds, const std::string method, + const double thres, const int32_t max_k, arma::uvec& sel_idxs, arma::vec& stats, arma::vec& pvalues, arma::uvec& rem_idxs, arma::uvec& sel_ord_idxs, const bool hash_on, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { - int sel_idx = -1; + int32_t sel_idx = -1; double sel_pvalue = 2; double sel_stat = 0; - std::vector rem_idxs_adj = keep_val(rem_idxs, 1); - for (unsigned int i = 0; i < rem_idxs_adj.size(); ++i) { + std::vector rem_idxs_adj = keep_val(rem_idxs, 1); + for (uint32_t i = 0; i < rem_idxs_adj.size(); ++i) { assoc_min(target_vars, ds, method, max_k, rem_idxs_adj[i], sel_idxs, stats, pvalues, sel_ord_idxs, hash_on, stats_kv, pvalues_kv); const double tmp_pvalue = STAT_PVALUE[PVALUE_POS]; @@ -323,16 +323,16 @@ unsigned int assoc_max_min(arma::vec& target_vars, arma::mat& ds, const std::str static arma::uvec get_sort_idxs(arma::vec& vals, arma::uvec& idxs) { std::vector vals_tmp; - for (unsigned int i = 0; i < idxs.size(); ++i) { + for (uint32_t i = 0; i < idxs.size(); ++i) { vals_tmp.push_back(vals[idxs[i]]); } arma::vec vals_adj(vals_tmp); return arma::sort_index(vals_adj); } -static arma::uvec get_idxs_eq(arma::uvec vals, const unsigned int val) { - std::vector vals_adj; - for (unsigned int i = 0; i < vals.size(); ++i) { +static arma::uvec get_idxs_eq(arma::uvec vals, const uint32_t val) { + std::vector vals_adj; + for (uint32_t i = 0; i < vals.size(); ++i) { if (vals[i] == val) { vals_adj.push_back(i); } @@ -340,8 +340,8 @@ static arma::uvec get_idxs_eq(arma::uvec vals, const unsigned int val) { return arma::uvec(vals_adj); } -static void neg_zero_idxs(arma::uvec& src, arma::uvec& base, const unsigned int val) { - for (unsigned int i = 0; i < src.size() && i < base.size(); ++i) { +static void neg_zero_idxs(arma::uvec& src, arma::uvec& base, const uint32_t val) { + for (uint32_t i = 0; i < src.size() && i < base.size(); ++i) { if (!base[i]) { src[i] = val; } @@ -349,7 +349,7 @@ static void neg_zero_idxs(arma::uvec& src, arma::uvec& base, const unsigned int } static bool is_true(arma::uvec& vals) { - for (unsigned int i = 0; i < vals.size(); ++i) { + for (uint32_t i = 0; i < vals.size(); ++i) { if (vals[i]) { return true; } @@ -358,7 +358,7 @@ static bool is_true(arma::uvec& vals) { } static void adj_gt_vals(arma::uvec& rem_idxs, arma::vec& pvalues, const double thres, const double val) { - for (unsigned int i = 0; i < rem_idxs.size(); ++i) { + for (uint32_t i = 0; i < rem_idxs.size(); ++i) { if (pvalues[i] > thres) { rem_idxs[i] = 0; } @@ -371,7 +371,7 @@ static void copy_vecs(Rcpp::List& src, arma::vec& dst1, arma::vec& dst2) { dst1.set_size(src1.size()); dst2.set_size(src2.size()); - for (unsigned int i = 0; i < src1.size(); ++i) { + for (uint32_t i = 0; i < src1.size(); ++i) { dst1[i] = src1[i]; dst2[i] = src2[i]; } @@ -380,7 +380,7 @@ static void copy_vecs(Rcpp::List& src, arma::vec& dst1, arma::vec& dst2) { static arma::vec calc_univ_pvalues(arma::vec& stats, const double dof) { arma::vec pvalues(stats.size()); static const double tmp_log = std::log(2); - for (unsigned int i = 0; i < stats.size(); ++i) { + for (uint32_t i = 0; i < stats.size(); ++i) { pvalues[i] = tmp_log + R::pt(std::abs(stats[i]), dof, false, true); } return pvalues; @@ -389,7 +389,7 @@ static arma::vec calc_univ_pvalues(arma::vec& stats, const double dof) { static void calc_univs(arma::vec& target_vars, arma::mat& ds, const std::string method, Rcpp::List& univs) { arma::mat cor_mat = arma::cor(target_vars, ds); arma::vec cor_vec = to_vec(cor_mat); - const unsigned int dof = ds.n_rows - 3; + const uint32_t dof = ds.n_rows - 3; arma::vec stats; if (!method.compare("pearson")) { db_print("True: !method.compare(\"pearson\")\n"); @@ -408,9 +408,9 @@ static double get_time_spent(const clock_t begin) { return (double) (end - begin) / CLOCKS_PER_SEC; } -void inter_mmp_c(arma::vec& target_vars, arma::mat& ds, const int max_k, +void inter_mmp_c(arma::vec& target_vars, arma::mat& ds, const int32_t max_k, const double thres, const std::string method, Rcpp::List inits, const bool hash_on, - const unsigned int var_size, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { + const uint32_t var_size, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { db_print("Initializing clock.\n"); clock_t begin = clock(); @@ -450,7 +450,7 @@ void inter_mmp_c(arma::vec& target_vars, arma::mat& ds, const int max_k, db_print("Initializing idxs.\n"); arma::uvec sel_idxs(var_size, arma::fill::zeros); arma::uvec sel_ord_idxs(var_size, arma::fill::zeros); - unsigned int sel_idx = arma::index_min(pvalues); + uint32_t sel_idx = arma::index_min(pvalues); sel_idxs[sel_idx] = 1; sel_ord_idxs[sel_idx] = 1; @@ -490,14 +490,14 @@ void inter_mmp_c(arma::vec& target_vars, arma::mat& ds, const int max_k, inits.size() ? res["n.tests"] = 0 : res["n.tests"] = stats.n_rows * stats.n_cols + kv_length; } -Rcpp::List calc_mmp_c_bp(arma::vec& target_vars, arma::mat& ds, const int limit_k, const double thres, +Rcpp::List calc_mmp_c_bp(arma::vec& target_vars, arma::mat& ds, const int32_t limit_k, const double thres, const std::string method) { db_print("Initializing values\n"); const double thres_log = std::log(thres); arma::uvec idxs(ds.n_cols); std::iota(idxs.begin(), idxs.end(), 0); - unsigned int cntr = 0; + uint32_t cntr = 0; std::vector pvalues; if (ds.n_cols == 1) { db_print("True: ds.n_cols == 1\n"); @@ -535,18 +535,18 @@ Rcpp::List calc_mmp_c_bp(arma::vec& target_vars, arma::mat& ds, const int limit_ } else { db_print("True: !ds.n_cols || ds.n_cols > 2\n"); - std::vector combn_vals(ds.n_cols); + std::vector combn_vals(ds.n_cols); std::iota(combn_vals.begin(), combn_vals.end(), 1); - for (unsigned int i = 0; i < ds.n_cols; ++i) { - int k = 0; + for (uint32_t i = 0; i < ds.n_cols; ++i) { + int32_t k = 0; double tmp_pvalue = -5.0; std::vector tmp_pvalues; - while (k < (int) ds.n_cols - 1 && k < limit_k && tmp_pvalue < thres_log) { - arma::umat combns = find_combn>(combn_vals, ++k); + while (k < (int32_t) ds.n_cols - 1 && k < limit_k && tmp_pvalue < thres_log) { + arma::umat combns = find_combn>(combn_vals, ++k); combns -= 1; - std::vector cols_rem = det_cols(combns, i); + std::vector cols_rem = det_cols(combns, i); arma::umat combns_adj = rm_cols(combns, cols_rem); - unsigned int j = 0; + uint32_t j = 0; while (j < combns_adj.n_cols && tmp_pvalue < thres_log) { arma::uvec tmp_idxs = combns_adj.col(j); calc_pearson(target_vars, ds, i, tmp_idxs, @@ -568,13 +568,13 @@ Rcpp::List calc_mmp_c_bp(arma::vec& target_vars, arma::mat& ds, const int limit_ static std::vector upd_vals(std::vector& src_pvalues, arma::uvec& idxs, std::vector dst_pvalues) { - for (unsigned int i = 0; i < idxs.size(); ++i) { + for (uint32_t i = 0; i < idxs.size(); ++i) { dst_pvalues[idxs[i]] = src_pvalues[i]; } return dst_pvalues; } -Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int max_k, +Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int32_t max_k, const double thres, const std::string method, Rcpp::List& inits, const bool hash_on, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv, const bool bws_on) { @@ -586,8 +586,8 @@ Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int max_k, db_print("Call: upd_args\n"); upd_args(target_vars, ds, method); - const unsigned int var_size = ds.n_cols; - if (max_k > (int) var_size) { + const uint32_t var_size = ds.n_cols; + if (max_k > (int32_t) var_size) { db_print("True: max_k > var_size\n"); max_k = var_size; } @@ -607,7 +607,7 @@ Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int max_k, std::vector res_im_pvalues = res["pvalues"]; std::vector res_mb_pvalues = res_mb["pvalue"]; res["pvalues"] = upd_vals(res_mb_pvalues, res_im_idxs, res_im_pvalues); - inits.size() ? res["n.tests"] = 0 : res["n.tests"] = (int) res["n.tests"] + (int) res_mb["counter"]; + inits.size() ? res["n.tests"] = 0 : res["n.tests"] = (int32_t) res["n.tests"] + (int32_t) res_mb["counter"]; } if (hash_on) { stats_kv[".length"] = kv_length; diff --git a/src/mmp_c.h b/src/mmp_c.h index 04fe75d..971a21c 100644 --- a/src/mmp_c.h +++ b/src/mmp_c.h @@ -14,7 +14,7 @@ // [[Rcpp::plugins("cpp11")]] // [[Rcpp::depends("RcppArmadillo")]] -Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int max_k, +Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int32_t max_k, const double thres, const std::string method, Rcpp::List& inits, const bool hash_on, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv, const bool bws_on); diff --git a/src/sf.cpp b/src/sf.cpp index 1a52b4c..5eafae4 100644 --- a/src/sf.cpp +++ b/src/sf.cpp @@ -21,7 +21,7 @@ END_RCPP } // [[Rcpp::export]] -Rcpp::List mmp_c(arma::vec target_vars, arma::mat ds, int max_k, +Rcpp::List mmp_c(arma::vec target_vars, arma::mat ds, int32_t max_k, const double thres, const std::string method, Rcpp::List inits, const bool hash_on, Rcpp::Environment stats_kv, Rcpp::Environment pvalues_kv, const bool bws_on) { @@ -35,7 +35,7 @@ BEGIN_RCPP RNGScope __rngScope; traits::input_parameter< arma::vec >::type target_vars(target_varsSEXP); traits::input_parameter< arma::mat >::type ds(dsSEXP); - traits::input_parameter< int >::type max_k(max_kSEXP); + traits::input_parameter< int32_t >::type max_k(max_kSEXP); traits::input_parameter< const double >::type thres(thresSEXP); traits::input_parameter< const std::string >::type method(methodSEXP); traits::input_parameter< List >::type inits(initsSEXP); From 475c5abbc26dfb363363a2dd387b7ce1a6fcf8f8 Mon Sep 17 00:00:00 2001 From: dimim Date: Sun, 17 Nov 2019 13:56:58 +0100 Subject: [PATCH 5/5] Apply formatting --- src/cts.cpp | 1009 +++++++++++++++++++++++------------------------ src/cts.h | 81 ++-- src/mmp_c.cpp | 1038 +++++++++++++++++++++++++------------------------ src/mmp_c.h | 21 +- src/sf.cpp | 60 +-- 5 files changed, 1124 insertions(+), 1085 deletions(-) diff --git a/src/cts.cpp b/src/cts.cpp index 9e8ecf0..3f542a8 100644 --- a/src/cts.cpp +++ b/src/cts.cpp @@ -4,14 +4,16 @@ #include "cts.h" #define DEBUG 0 -#define db_print(...) \ - do { if (DEBUG) Rprintf(__VA_ARGS__); } while (0) +#define db_print(...) \ + do { \ + if (DEBUG) Rprintf(__VA_ARGS__); \ + } while (0) static uint32_t skip_ahead(arma::uvec& vals, const uint32_t curr) { - uint32_t i; - for (i = curr + 1; i < vals.size() && vals(i) == vals(curr); ++i) { - } - return i; + uint32_t i; + for (i = curr + 1; i < vals.size() && vals(i) == vals(curr); ++i) { + } + return i; } double calc_med(arma::vec& vals); @@ -19,23 +21,25 @@ double calc_med(arma::vec& vals); double calc_med(arma::vec& vals, const uint32_t size); bool adj_med_NAs(arma::mat& ds) { - bool found_NA = false; - for (uint32_t i = 0; i < ds.n_rows; i++) { - bool found_NA_col = false; - double med = 0; - for (uint32_t j = 0; j < ds.n_cols; j++) { - if (!arma::is_finite(ds.at(i, j))) { - if (!found_NA) { found_NA = true; } - if (!found_NA_col) { - found_NA_col = true; - arma::vec tmp = ds.col(j); - med = calc_med(tmp); - } - ds.at(i, j) = med; - } - } - } - return found_NA; + bool found_NA = false; + for (uint32_t i = 0; i < ds.n_rows; i++) { + bool found_NA_col = false; + double med = 0; + for (uint32_t j = 0; j < ds.n_cols; j++) { + if (!arma::is_finite(ds.at(i, j))) { + if (!found_NA) { + found_NA = true; + } + if (!found_NA_col) { + found_NA_col = true; + arma::vec tmp = ds.col(j); + med = calc_med(tmp); + } + ds.at(i, j) = med; + } + } + } + return found_NA; } double calc_med(arma::vec& vals) { @@ -43,9 +47,10 @@ double calc_med(arma::vec& vals) { double med = 0; if (vals.size() % 2 == 0) { std::nth_element(vals.begin(), vals.begin() + mid, vals.end()); - med = (vals[mid] + *std::min_element(vals.begin() + (mid + 1), vals.end())) / 2.0; - } - else { + med = (vals[mid] + + *std::min_element(vals.begin() + (mid + 1), vals.end())) / + 2.0; + } else { std::nth_element(vals.begin(), vals.begin() + (mid + 1), vals.end()); med = vals[mid + 1]; } @@ -55,145 +60,145 @@ double calc_med(arma::vec& vals) { double find_freq(arma::vec& vals); bool adj_freq_NAs(arma::mat& ds) { - bool found_NA = false; - for (uint32_t i = 0; i < ds.n_rows; i++) { - bool found_NA_col = false; - double freq = 0; - for (uint32_t j = 0; j < ds.n_cols; j++) { - if (!arma::is_finite(ds.at(i, j))) { - if (!found_NA) { found_NA = true; } - if (!found_NA_col) { - found_NA_col = true; - arma::vec tmp = ds.col(j); - freq = find_freq(tmp); - } - ds.at(i, j) = freq; - } - } - } - return found_NA; + bool found_NA = false; + for (uint32_t i = 0; i < ds.n_rows; i++) { + bool found_NA_col = false; + double freq = 0; + for (uint32_t j = 0; j < ds.n_cols; j++) { + if (!arma::is_finite(ds.at(i, j))) { + if (!found_NA) { + found_NA = true; + } + if (!found_NA_col) { + found_NA_col = true; + arma::vec tmp = ds.col(j); + freq = find_freq(tmp); + } + ds.at(i, j) = freq; + } + } + } + return found_NA; } double find_freq(arma::vec& vals) { - std::sort(vals.begin(), vals.end()); - double prev = vals[0]; - double freq = vals[0]; - int32_t curr_count = 1; - int32_t max_count = 1; - for (uint32_t i = 1; i < vals.size(); i++) { - if (vals[i] == prev) { - curr_count++; - } - else { - if (curr_count > max_count) { - freq = vals[i - 1]; - max_count = curr_count; - } - prev = vals[i]; - curr_count = 1; - } - } - if (curr_count == 1 && max_count == 1) { - return *std::min_element(vals.begin(), vals.end()); - } - return curr_count > max_count ? vals[vals.size() - 1] : freq; + std::sort(vals.begin(), vals.end()); + double prev = vals[0]; + double freq = vals[0]; + int32_t curr_count = 1; + int32_t max_count = 1; + for (uint32_t i = 1; i < vals.size(); i++) { + if (vals[i] == prev) { + curr_count++; + } else { + if (curr_count > max_count) { + freq = vals[i - 1]; + max_count = curr_count; + } + prev = vals[i]; + curr_count = 1; + } + } + if (curr_count == 1 && max_count == 1) { + return *std::min_element(vals.begin(), vals.end()); + } + return curr_count > max_count ? vals[vals.size() - 1] : freq; } arma::mat calc_rank(arma::mat& ds) { - arma::mat rds(ds.n_rows, ds.n_cols); - for (uint32_t i = 0; i < ds.n_cols; ++i) { - arma::vec curr_col = ds.col(i); - rds.col(i) = rank_mean(curr_col, false); - } - return rds; -} - + arma::mat rds(ds.n_rows, ds.n_cols); + for (uint32_t i = 0; i < ds.n_cols; ++i) { + arma::vec curr_col = ds.col(i); + rds.col(i) = + rank_mean(curr_col, false); + } + return rds; +} + arma::uvec sub_col_max_min(arma::mat& ds, const bool cont) { - int32_t extra_val = 0; - if (!cont) { - extra_val = 1; - } - arma::uvec max_min(ds.n_cols); - for (uint32_t j = 0; j < ds.n_cols; j++) { - double max = ds(0, j); - double min = ds(0, j); - for (uint32_t i = 1; i < ds.n_rows; i++) { - const uint32_t curr = ds(i, j); - if (curr > max) { - max = curr; - } - if (curr < min) { - min = curr; - } - } - max_min(j) = (max - min) + extra_val; - } - return max_min; -} - -arma::mat calc_pt(arma::mat& ds, - const int32_t df, const bool lower_tail, const bool log_p, - const double add) { - arma::mat pt_ds(ds.n_rows, ds.n_cols); - for (uint32_t i = 0; i < ds.n_rows; i++) { - for (uint32_t j = 0; j < ds.n_cols; j++) { - pt_ds(i, j) = add + R::pt(ds(i, j), df, lower_tail, log_p); - } - } - return pt_ds; + int32_t extra_val = 0; + if (!cont) { + extra_val = 1; + } + arma::uvec max_min(ds.n_cols); + for (uint32_t j = 0; j < ds.n_cols; j++) { + double max = ds(0, j); + double min = ds(0, j); + for (uint32_t i = 1; i < ds.n_rows; i++) { + const uint32_t curr = ds(i, j); + if (curr > max) { + max = curr; + } + if (curr < min) { + min = curr; + } + } + max_min(j) = (max - min) + extra_val; + } + return max_min; +} + +arma::mat calc_pt(arma::mat& ds, const int32_t df, const bool lower_tail, + const bool log_p, const double add) { + arma::mat pt_ds(ds.n_rows, ds.n_cols); + for (uint32_t i = 0; i < ds.n_rows; i++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { + pt_ds(i, j) = add + R::pt(ds(i, j), df, lower_tail, log_p); + } + } + return pt_ds; } std::vector det_cols(arma::umat& ds, const uint32_t val) { - std::vector cols; - for (uint32_t i = 0; i < ds.n_rows; ++i) { - for (uint32_t j = 0; j < ds.n_cols; ++j) { - if (ds(i, j) == val) { - cols.push_back(j); - } - } - } - std::sort(cols.begin(), cols.end()); - return cols; -} - -arma::mat ext_cols(arma::mat& ds, const uint32_t col_a, - const uint32_t col_b) { - const uint32_t ncols = 2; - arma::mat ext_ds(ds.n_rows, ncols); - const uint32_t ext_col_a = 0; - const uint32_t ext_col_b = 1; - for (uint32_t i = 0; i < ds.n_rows; i++) { - ext_ds(i, ext_col_a) = ds(i, col_a); - ext_ds(i, ext_col_b) = ds(i, col_b); - } - return ext_ds; + std::vector cols; + for (uint32_t i = 0; i < ds.n_rows; ++i) { + for (uint32_t j = 0; j < ds.n_cols; ++j) { + if (ds(i, j) == val) { + cols.push_back(j); + } + } + } + std::sort(cols.begin(), cols.end()); + return cols; +} + +arma::mat ext_cols(arma::mat& ds, const uint32_t col_a, const uint32_t col_b) { + const uint32_t ncols = 2; + arma::mat ext_ds(ds.n_rows, ncols); + const uint32_t ext_col_a = 0; + const uint32_t ext_col_b = 1; + for (uint32_t i = 0; i < ds.n_rows; i++) { + ext_ds(i, ext_col_a) = ds(i, col_a); + ext_ds(i, ext_col_b) = ds(i, col_b); + } + return ext_ds; } void cp_lt(arma::mat& src, arma::mat& dst, const int32_t val) { - for (uint32_t i = 0; i < src.n_rows; i++) { - for (uint32_t j = 0; j < src.n_cols; j++) { - i > j ? dst(i, j) = val : dst(i, j) = src(i, j); - } - } + for (uint32_t i = 0; i < src.n_rows; i++) { + for (uint32_t j = 0; j < src.n_cols; j++) { + i > j ? dst(i, j) = val : dst(i, j) = src(i, j); + } + } } void adj_diag(arma::mat& ds, const double val) { - for (uint32_t i = 0; i < ds.n_rows && i < ds.n_cols; i++) { - ds(i, i) = val; - } + for (uint32_t i = 0; i < ds.n_rows && i < ds.n_cols; i++) { + ds(i, i) = val; + } } arma::mat cbind_tran_mat(arma::mat& ds1, arma::mat& ds2) { - const uint32_t nrows = ds1.n_cols; - const uint32_t ncols = ds1.n_rows + ds2.n_rows; - arma::mat ds(nrows, ncols); - for (uint32_t i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { - for (uint32_t j = 0; j < ds1.n_cols && j < ds2.n_cols; j++) { - ds(j, i) = ds1(i, j); - ds(j, i + ds1.n_rows) = ds2(i, j); - } - } - return ds; + const uint32_t nrows = ds1.n_cols; + const uint32_t ncols = ds1.n_rows + ds2.n_rows; + arma::mat ds(nrows, ncols); + for (uint32_t i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { + for (uint32_t j = 0; j < ds1.n_cols && j < ds2.n_cols; j++) { + ds(j, i) = ds1(i, j); + ds(j, i + ds1.n_rows) = ds2(i, j); + } + } + return ds; } bool is_dupl_row(arma::mat& ds, const uint32_t lrow); @@ -202,445 +207,447 @@ std::vector get_dupl_rows_pos(arma::mat& ds); arma::mat rm_dupl_rows(arma::mat& src); -arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, - const bool ass1, const bool ass2) { - const uint32_t nrows = ds1.n_rows + ds2.n_rows; - uint32_t ncols; - ds1.n_cols > ds2.n_cols ? ncols = ds1.n_cols : ncols = ds2.n_cols; - arma::mat ds(nrows, ncols, arma::fill::zeros); - uint32_t row = 0; - uint32_t col = 0; - if (ass1) { - for (uint32_t i = 0; i < ds1.n_rows; i++) { - for (uint32_t j = 0; j < ds1.n_cols; j++) { - ds(row, col) = ds1(i, j); - col++; - } - row++; - col = 0; - } - } - else { - row = ds1.n_rows; - } - if (ass2) { - for (uint32_t i = 0; i < ds2.n_rows; i++) { - for (uint32_t j = 0; j < ds2.n_cols; j++) { - ds(row, col) = ds2(i, j); - col++; - } - row++; - col = 0; - } - } - return rm_dupl_rows(ds); +arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, const bool ass1, + const bool ass2) { + const uint32_t nrows = ds1.n_rows + ds2.n_rows; + uint32_t ncols; + ds1.n_cols > ds2.n_cols ? ncols = ds1.n_cols : ncols = ds2.n_cols; + arma::mat ds(nrows, ncols, arma::fill::zeros); + uint32_t row = 0; + uint32_t col = 0; + if (ass1) { + for (uint32_t i = 0; i < ds1.n_rows; i++) { + for (uint32_t j = 0; j < ds1.n_cols; j++) { + ds(row, col) = ds1(i, j); + col++; + } + row++; + col = 0; + } + } else { + row = ds1.n_rows; + } + if (ass2) { + for (uint32_t i = 0; i < ds2.n_rows; i++) { + for (uint32_t j = 0; j < ds2.n_cols; j++) { + ds(row, col) = ds2(i, j); + col++; + } + row++; + col = 0; + } + } + return rm_dupl_rows(ds); } arma::mat rm_dupl_rows(arma::mat& src) { - std::vector dupls_pos = get_dupl_rows_pos(src); - if (dupls_pos.empty()) { - return src; - } - const uint32_t nrows = src.n_rows - dupls_pos.size(); - arma::mat dst(nrows, src.n_cols); - uint32_t dupls_pos_cntr = 0; - for (uint32_t i = 0, src_row = 0; i < dst.n_rows; i++, src_row++) { - while (dupls_pos_cntr < dupls_pos.size() && - dupls_pos.at(dupls_pos_cntr) == src_row) { - src_row++; - dupls_pos_cntr++; - } - for (uint32_t j = 0; j < dst.n_cols; j++) { - dst(i, j) = src(src_row, j); - } - } - return dst; + std::vector dupls_pos = get_dupl_rows_pos(src); + if (dupls_pos.empty()) { + return src; + } + const uint32_t nrows = src.n_rows - dupls_pos.size(); + arma::mat dst(nrows, src.n_cols); + uint32_t dupls_pos_cntr = 0; + for (uint32_t i = 0, src_row = 0; i < dst.n_rows; i++, src_row++) { + while (dupls_pos_cntr < dupls_pos.size() && + dupls_pos.at(dupls_pos_cntr) == src_row) { + src_row++; + dupls_pos_cntr++; + } + for (uint32_t j = 0; j < dst.n_cols; j++) { + dst(i, j) = src(src_row, j); + } + } + return dst; } std::vector get_dupl_rows_pos(arma::mat& ds) { - std::vector dupls_pos; - for (uint32_t i = 1; i < ds.n_rows; i++) { - if (is_dupl_row(ds, i)) { - dupls_pos.push_back(i); - } - } - return dupls_pos; + std::vector dupls_pos; + for (uint32_t i = 1; i < ds.n_rows; i++) { + if (is_dupl_row(ds, i)) { + dupls_pos.push_back(i); + } + } + return dupls_pos; } bool is_dupl_row(arma::mat& ds, const uint32_t lrow) { - for (uint32_t i = 0; i < lrow; i++) { - for (uint32_t j = 0; - ; - j++) { - if (ds(i, j) != ds(lrow, j)) { - break; - } - if (j == ds.n_cols - 1) { - return true; - } - } - } - return false; + for (uint32_t i = 0; i < lrow; i++) { + for (uint32_t j = 0;; j++) { + if (ds(i, j) != ds(lrow, j)) { + break; + } + if (j == ds.n_cols - 1) { + return true; + } + } + } + return false; } arma::vec to_vec(arma::mat& ds) { - arma::vec vals(ds.n_rows * ds.n_cols); - uint32_t k = 0; - for (uint32_t i = 0; i < ds.n_rows; i++) { - for (uint32_t j = 0; j < ds.n_cols; j++) { - vals(k++) = ds(i, j); - } - } - return vals; + arma::vec vals(ds.n_rows * ds.n_cols); + uint32_t k = 0; + for (uint32_t i = 0; i < ds.n_rows; i++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { + vals(k++) = ds(i, j); + } + } + return vals; } bool found_match(const uint32_t x, arma::uvec& vals) { - for (uint32_t i = 0; i < vals.size(); ++i) { - if (x == vals[i]) { - return true; - } - } - return false; -} - -bool are_equal(arma::mat& ds, arma::vec& vals, - const bool by_col_idx, const unsigned col) { - if ((by_col_idx && ds.n_rows != vals.size() && ds.n_cols != vals.size()) || - (!by_col_idx && (ds.n_rows * ds.n_cols) != vals.size())) { - return false; - } - uint32_t j = 0; - uint32_t k = 0; - by_col_idx ? j = col : j = 0; - for (; (!by_col_idx && j < ds.n_cols) || (by_col_idx && col == j); ++j) { - for (uint32_t i = 0; i < ds.n_rows; ++i) { - if (ds(i, j) != vals[k++]) { - return false; - } - } - } - return true; + for (uint32_t i = 0; i < vals.size(); ++i) { + if (x == vals[i]) { + return true; + } + } + return false; +} + +bool are_equal(arma::mat& ds, arma::vec& vals, const bool by_col_idx, + const unsigned col) { + if ((by_col_idx && ds.n_rows != vals.size() && ds.n_cols != vals.size()) || + (!by_col_idx && (ds.n_rows * ds.n_cols) != vals.size())) { + return false; + } + uint32_t j = 0; + uint32_t k = 0; + by_col_idx ? j = col : j = 0; + for (; (!by_col_idx && j < ds.n_cols) || (by_col_idx && col == j); ++j) { + for (uint32_t i = 0; i < ds.n_rows; ++i) { + if (ds(i, j) != vals[k++]) { + return false; + } + } + } + return true; } std::vector rm_lt_nan(arma::uvec& idxs, const uint32_t limit) { - std::vector idxs_adj; - for (uint32_t i = 0; i < idxs.size(); ++i) { - if (arma::is_finite(idxs[i]) || idxs[i] >= limit) { - idxs_adj.push_back(idxs[i]); - } - } - return idxs_adj; + std::vector idxs_adj; + for (uint32_t i = 0; i < idxs.size(); ++i) { + if (arma::is_finite(idxs[i]) || idxs[i] >= limit) { + idxs_adj.push_back(idxs[i]); + } + } + return idxs_adj; } -std::unordered_set get_diff(std::vector& lh, const uint32_t val) { - std::unordered_set diff; - for (uint32_t i = 0; i < lh.size(); ++i) { - if (lh[i] != val) { - diff.insert(lh[i]); - } - } - return diff; +std::unordered_set get_diff(std::vector& lh, + const uint32_t val) { + std::unordered_set diff; + for (uint32_t i = 0; i < lh.size(); ++i) { + if (lh[i] != val) { + diff.insert(lh[i]); + } + } + return diff; } arma::umat nchoosek(std::vector& idxs, const uint32_t k) { - if (idxs.size() == 1) { - arma::umat ret(1, 1); - ret(0, 0) = R::choose(idxs[0], k); - return ret; - } - return find_combn>(idxs, k); + if (idxs.size() == 1) { + arma::umat ret(1, 1); + ret(0, 0) = R::choose(idxs[0], k); + return ret; + } + return find_combn>(idxs, k); } std::vector inter_helper(arma::vec& vals1, arma::vec& vals2) { - std::sort(vals1.begin(), vals1.end()); - std::sort(vals2.begin(), vals2.end()); - std::vector inter_vals; - for (uint32_t i = 0, j = 0; i < vals1.size() && j < vals2.size();) { - const double curr1 = vals1(i); - const double curr2 = vals2(j); - if (curr1 == curr2) { - const uint32_t size = inter_vals.size(); - if (!size || (inter_vals.back() != curr1)) { - inter_vals.push_back(curr1); - } - i++; - j++; - } - else if (curr1 < curr2) { - if (vals1(vals1.size() - 1) < curr2) { - break; - } - i++; - } - else { - if (vals2(vals2.size() - 1) < curr1) { - break; - } - j++; - } - } - return inter_vals; + std::sort(vals1.begin(), vals1.end()); + std::sort(vals2.begin(), vals2.end()); + std::vector inter_vals; + for (uint32_t i = 0, j = 0; i < vals1.size() && j < vals2.size();) { + const double curr1 = vals1(i); + const double curr2 = vals2(j); + if (curr1 == curr2) { + const uint32_t size = inter_vals.size(); + if (!size || (inter_vals.back() != curr1)) { + inter_vals.push_back(curr1); + } + i++; + j++; + } else if (curr1 < curr2) { + if (vals1(vals1.size() - 1) < curr2) { + break; + } + i++; + } else { + if (vals2(vals2.size() - 1) < curr1) { + break; + } + j++; + } + } + return inter_vals; } // Alters row_idxs -void append_rows(arma::mat& ds, const double val, std::vector& row_idxs); +void append_rows(arma::mat& ds, const double val, + std::vector& row_idxs); std::vector index_row_eq(arma::mat& ds, std::vector& vals) { - std::vector row_idxs; - for (uint32_t i = 0; i < vals.size(); i++) { - append_rows(ds, vals.at(i), row_idxs); - } - std::sort(row_idxs.begin(), row_idxs.end()); - row_idxs.erase(std::unique(row_idxs.begin(), row_idxs.end()), row_idxs.end()); - return row_idxs; -} - -void append_rows(arma::mat& ds, const double val, std::vector& row_idxs) { - for (uint32_t i = 0; i < ds.n_rows; i++) { - for (uint32_t j = 0; j < ds.n_cols; j++) { - if (ds(i, j) == val) { - row_idxs.push_back(i); - } - } - } + std::vector row_idxs; + for (uint32_t i = 0; i < vals.size(); i++) { + append_rows(ds, vals.at(i), row_idxs); + } + std::sort(row_idxs.begin(), row_idxs.end()); + row_idxs.erase(std::unique(row_idxs.begin(), row_idxs.end()), + row_idxs.end()); + return row_idxs; +} + +void append_rows(arma::mat& ds, const double val, + std::vector& row_idxs) { + for (uint32_t i = 0; i < ds.n_rows; i++) { + for (uint32_t j = 0; j < ds.n_cols; j++) { + if (ds(i, j) == val) { + row_idxs.push_back(i); + } + } + } } arma::mat rm_rows(arma::mat& src, arma::uvec& rows) { - uint32_t dst_nrows = src.n_rows - rows.size(); - uint32_t dst_ncols = src.n_cols; - arma::mat dst(dst_nrows, dst_ncols); - uint32_t src_row = 0; - uint32_t rows_idx = 0; - for (uint32_t dst_row = 0; dst_row < dst_nrows; dst_row++) { - while (src_row < src.n_rows && rows_idx < rows.size() && src_row == rows(rows_idx)) { - src_row++; - rows_idx = skip_ahead(rows, rows_idx); - } - for (uint32_t col = 0; col < dst_ncols; col++) { - dst(dst_row, col) = src(src_row, col); - } - src_row++; - } - return dst; + uint32_t dst_nrows = src.n_rows - rows.size(); + uint32_t dst_ncols = src.n_cols; + arma::mat dst(dst_nrows, dst_ncols); + uint32_t src_row = 0; + uint32_t rows_idx = 0; + for (uint32_t dst_row = 0; dst_row < dst_nrows; dst_row++) { + while (src_row < src.n_rows && rows_idx < rows.size() && + src_row == rows(rows_idx)) { + src_row++; + rows_idx = skip_ahead(rows, rows_idx); + } + for (uint32_t col = 0; col < dst_ncols; col++) { + dst(dst_row, col) = src(src_row, col); + } + src_row++; + } + return dst; } uint32_t skip_ahead_std(std::vector rows, const uint32_t curr); arma::mat rm_rows_std(arma::mat& src, std::vector& rows) { - uint32_t dst_nrows = src.n_rows - rows.size(); - uint32_t dst_ncols = src.n_cols; - arma::mat dst(dst_nrows, dst_ncols); - uint32_t src_row = 0; - uint32_t rows_idx = 0; - for (uint32_t dst_row = 0; dst_row < dst_nrows; dst_row++) { - while (src_row < src.n_rows && rows_idx < rows.size() && src_row == rows.at(rows_idx)) { - src_row++; - rows_idx = skip_ahead_std(rows, rows_idx); - } - for (uint32_t col = 0; col < dst_ncols; col++) { - dst(dst_row, col) = src(src_row, col); - } - src_row++; - } - return dst; + uint32_t dst_nrows = src.n_rows - rows.size(); + uint32_t dst_ncols = src.n_cols; + arma::mat dst(dst_nrows, dst_ncols); + uint32_t src_row = 0; + uint32_t rows_idx = 0; + for (uint32_t dst_row = 0; dst_row < dst_nrows; dst_row++) { + while (src_row < src.n_rows && rows_idx < rows.size() && + src_row == rows.at(rows_idx)) { + src_row++; + rows_idx = skip_ahead_std(rows, rows_idx); + } + for (uint32_t col = 0; col < dst_ncols; col++) { + dst(dst_row, col) = src(src_row, col); + } + src_row++; + } + return dst; } uint32_t skip_ahead_std(std::vector rows, const uint32_t curr) { - uint32_t i; - for (i = curr + 1; i < rows.size() && rows.at(i) == rows.at(curr); i++) { - } - return i; + uint32_t i; + for (i = curr + 1; i < rows.size() && rows.at(i) == rows.at(curr); i++) { + } + return i; } // No sorting arma::umat rm_cols(arma::umat src, std::vector cols) { - uint32_t dst_nrows = src.n_rows; - uint32_t dst_ncols = src.n_cols - cols.size(); - arma::umat dst(dst_nrows, dst_ncols); - - uint32_t src_col = 0; - uint32_t cols_idx = 0; - for (uint32_t dst_col = 0; dst_col < dst_ncols; ++dst_col) { - while (src_col < src.n_cols && cols_idx < cols.size() && src_col == cols.at(cols_idx)) { - ++src_col; - ++cols_idx; - } - for (uint32_t row = 0; row < dst_nrows; ++row) { - dst(row, dst_col) = src(row, src_col); - } - ++src_col; - } - return dst; + uint32_t dst_nrows = src.n_rows; + uint32_t dst_ncols = src.n_cols - cols.size(); + arma::umat dst(dst_nrows, dst_ncols); + + uint32_t src_col = 0; + uint32_t cols_idx = 0; + for (uint32_t dst_col = 0; dst_col < dst_ncols; ++dst_col) { + while (src_col < src.n_cols && cols_idx < cols.size() && + src_col == cols.at(cols_idx)) { + ++src_col; + ++cols_idx; + } + for (uint32_t row = 0; row < dst_nrows; ++row) { + dst(row, dst_col) = src(row, src_col); + } + ++src_col; + } + return dst; } arma::mat order_col(arma::mat& ds, const uint32_t col) { - arma::mat ordered_ds(ds.n_rows, ds.n_cols); - arma::uvec order_idxs = arma::sort_index(ds.col(col)); - for (uint32_t i = 0; i < ds.n_rows; i++) { - const uint32_t idx = order_idxs(i); - for (uint32_t j = 0; j < ds.n_cols; j++) { - ordered_ds(i, j) = ds(idx, j); - } - } - return ordered_ds; + arma::mat ordered_ds(ds.n_rows, ds.n_cols); + arma::uvec order_idxs = arma::sort_index(ds.col(col)); + for (uint32_t i = 0; i < ds.n_rows; i++) { + const uint32_t idx = order_idxs(i); + for (uint32_t j = 0; j < ds.n_cols; j++) { + ordered_ds(i, j) = ds(idx, j); + } + } + return ordered_ds; } arma::mat form_cmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols) { - arma::mat formed_ds(cols.size(), rows.size()); - for (uint32_t i = 0; i < rows.size(); i++) { - for (uint32_t j = 0; j < cols.size(); j++) { - formed_ds(j, i) = ds(rows(i), cols(j)); - } - } - return formed_ds; + arma::mat formed_ds(cols.size(), rows.size()); + for (uint32_t i = 0; i < rows.size(); i++) { + for (uint32_t j = 0; j < cols.size(); j++) { + formed_ds(j, i) = ds(rows(i), cols(j)); + } + } + return formed_ds; } arma::mat form_rmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols) { - arma::mat formed_ds(rows.size(), cols.size()); - for (uint32_t i = 0; i < rows.size(); i++) { - for (uint32_t j = 0; j < cols.size(); j++) { - formed_ds(i, j) = ds(rows(i), cols(j)); - } - } - return formed_ds; -} - -arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, - std::vector& cols) { - arma::mat formed_ds(rows.size(), cols.size()); - for (uint32_t i = 0; i < rows.size(); i++) { - for (uint32_t j = 0; j < cols.size(); j++) { - formed_ds(i, j) = ds(rows.at(i), cols.at(j)); - } - } - return formed_ds; + arma::mat formed_ds(rows.size(), cols.size()); + for (uint32_t i = 0; i < rows.size(); i++) { + for (uint32_t j = 0; j < cols.size(); j++) { + formed_ds(i, j) = ds(rows(i), cols(j)); + } + } + return formed_ds; +} + +arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, + std::vector& cols) { + arma::mat formed_ds(rows.size(), cols.size()); + for (uint32_t i = 0; i < rows.size(); i++) { + for (uint32_t j = 0; j < cols.size(); j++) { + formed_ds(i, j) = ds(rows.at(i), cols.at(j)); + } + } + return formed_ds; } arma::mat merge_cols(arma::mat& ds, arma::uvec& idxs) { - arma::mat merged_ds(ds.n_rows, idxs.size()); - for (uint32_t i = 0; i < idxs.size(); i++) { - for (uint32_t j = 0; j < ds.n_rows; j++) { - merged_ds(j, i) = ds(j, idxs(i)); - } - } - return merged_ds; + arma::mat merged_ds(ds.n_rows, idxs.size()); + for (uint32_t i = 0; i < idxs.size(); i++) { + for (uint32_t j = 0; j < ds.n_rows; j++) { + merged_ds(j, i) = ds(j, idxs(i)); + } + } + return merged_ds; } arma::mat form_ncolcmat(arma::vec& vals, arma::mat& ds) { - const uint32_t nrows = vals.size(); - const uint32_t ncols = 1 + ds.n_cols; - arma::mat formed_ds(nrows, ncols); - for (uint32_t i = 0; i < nrows; i++) { - formed_ds(i, 0) = vals(i); - for (uint32_t j = 1; j < ncols; j++) { - formed_ds(i, j) = ds(i, (j - 1)); - } - } - return formed_ds; + const uint32_t nrows = vals.size(); + const uint32_t ncols = 1 + ds.n_cols; + arma::mat formed_ds(nrows, ncols); + for (uint32_t i = 0; i < nrows; i++) { + formed_ds(i, 0) = vals(i); + for (uint32_t j = 1; j < ncols; j++) { + formed_ds(i, j) = ds(i, (j - 1)); + } + } + return formed_ds; } arma::mat form_c2mat(arma::vec& vals1, arma::vec& vals2) { - const uint32_t nrows = vals1.size(); - const uint32_t ncols = 2; - arma::mat formed_ds(nrows, ncols); - for (uint32_t i = 0; i < nrows; i++) { - formed_ds(i, 0) = vals1(i); - formed_ds(i, 1) = vals2(i); - } - return formed_ds; + const uint32_t nrows = vals1.size(); + const uint32_t ncols = 2; + arma::mat formed_ds(nrows, ncols); + for (uint32_t i = 0; i < nrows; i++) { + formed_ds(i, 0) = vals1(i); + formed_ds(i, 1) = vals2(i); + } + return formed_ds; } arma::uvec form_vec(arma::mat& ds, const uint32_t row, arma::uvec& cols) { - arma::uvec vals(cols.size()); - for (uint32_t i = 0; i < cols.size(); i++) { - vals(i) = ds(row, cols(i)); - } - return vals; + arma::uvec vals(cols.size()); + for (uint32_t i = 0; i < cols.size(); i++) { + vals(i) = ds(row, cols(i)); + } + return vals; } arma::vec form_vec_wvals(arma::mat& ds, const uint32_t row, arma::uvec& cols, - arma::vec& vals) { - arma::vec formed_vec(cols.size() + vals.size()); - uint32_t i; - for (i = 0; i < cols.size(); i++) { - formed_vec(i) = ds(row, cols(i)); - } - for (uint32_t j = 0; i < formed_vec.size(); i++, j++) { - formed_vec(i) = vals(j); - } - return formed_vec; + arma::vec& vals) { + arma::vec formed_vec(cols.size() + vals.size()); + uint32_t i; + for (i = 0; i < cols.size(); i++) { + formed_vec(i) = ds(row, cols(i)); + } + for (uint32_t j = 0; i < formed_vec.size(); i++, j++) { + formed_vec(i) = vals(j); + } + return formed_vec; } arma::mat append_row(arma::mat& ds, const uint32_t row, arma::vec& vals) { - for (uint32_t j = 0; j < ds.n_cols; j++) { - ds(row, j) = vals(j); - } - return ds; + for (uint32_t j = 0; j < ds.n_cols; j++) { + ds(row, j) = vals(j); + } + return ds; } std::vector rsum_gt_zero_idxs(arma::mat& ds) { - std::vector idxs; - for (uint32_t i = 0; i < ds.n_rows; i++) { - double sum = 0; - for (uint32_t j = 0; j < ds.n_cols; j++) { - sum += ds(i, j); - } - if (sum > 0) { - idxs.push_back(i); - } - } - return idxs; + std::vector idxs; + for (uint32_t i = 0; i < ds.n_rows; i++) { + double sum = 0; + for (uint32_t j = 0; j < ds.n_cols; j++) { + sum += ds(i, j); + } + if (sum > 0) { + idxs.push_back(i); + } + } + return idxs; } arma::vec form_cmat_vec(arma::mat& ds, arma::rowvec& vals) { - arma::vec concat_vals(ds.n_rows * ds.n_cols + vals.size()); - uint32_t cd_cntr = 0; - for (uint32_t j = 0; j < ds.n_cols; j++) { - for (uint32_t i = 0; i < ds.n_rows; i++) { - concat_vals(cd_cntr++) = ds(i, j); - } - } - for (uint32_t i = 0; i < vals.size(); i++) { - concat_vals[cd_cntr++] = vals(i); - } - return concat_vals; + arma::vec concat_vals(ds.n_rows * ds.n_cols + vals.size()); + uint32_t cd_cntr = 0; + for (uint32_t j = 0; j < ds.n_cols; j++) { + for (uint32_t i = 0; i < ds.n_rows; i++) { + concat_vals(cd_cntr++) = ds(i, j); + } + } + for (uint32_t i = 0; i < vals.size(); i++) { + concat_vals[cd_cntr++] = vals(i); + } + return concat_vals; } arma::mat cbind_mat(arma::mat& ds1, arma::mat& ds2) { - const uint32_t nrows = ds1.n_rows; - const uint32_t ncols = ds1.n_cols + ds2.n_cols; - arma::mat ds(nrows, ncols); - for (uint32_t i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { - for (uint32_t j = 0; j < ds1.n_cols; j++) { - ds(i, j) = ds1(i, j); - } - for (uint32_t j = 0; j < ds2.n_cols; j++) { - ds(i, j + ds1.n_cols) = ds2(i, j); - } - } - return ds; + const uint32_t nrows = ds1.n_rows; + const uint32_t ncols = ds1.n_cols + ds2.n_cols; + arma::mat ds(nrows, ncols); + for (uint32_t i = 0; i < ds1.n_rows && i < ds2.n_rows; i++) { + for (uint32_t j = 0; j < ds1.n_cols; j++) { + ds(i, j) = ds1(i, j); + } + for (uint32_t j = 0; j < ds2.n_cols; j++) { + ds(i, j + ds1.n_cols) = ds2(i, j); + } + } + return ds; } arma::mat adj_cols(arma::mat& src, const uint32_t dst_ncols) { - const uint32_t dst_nrows = src.n_rows * src.n_cols / dst_ncols; - arma::mat dst(dst_nrows, dst_ncols); - uint32_t src_row = 0; - uint32_t src_col = 0; - uint32_t dst_row = 0; - uint32_t dst_col = 0; - while (src_col < src.n_cols && dst_col < dst_ncols) { - while (src_row < src.n_rows && dst_row < dst_nrows) { - dst(dst_row++, dst_col) = src(src_row++, src_col); - } - if (src_row >= src.n_rows) { - src_row = 0; - src_col++; - } - if (dst_row >= dst_nrows) { - dst_row = 0; - dst_col++; - } - } - return dst; + const uint32_t dst_nrows = src.n_rows * src.n_cols / dst_ncols; + arma::mat dst(dst_nrows, dst_ncols); + uint32_t src_row = 0; + uint32_t src_col = 0; + uint32_t dst_row = 0; + uint32_t dst_col = 0; + while (src_col < src.n_cols && dst_col < dst_ncols) { + while (src_row < src.n_rows && dst_row < dst_nrows) { + dst(dst_row++, dst_col) = src(src_row++, src_col); + } + if (src_row >= src.n_rows) { + src_row = 0; + src_col++; + } + if (dst_row >= dst_nrows) { + dst_row = 0; + dst_col++; + } + } + return dst; } diff --git a/src/cts.h b/src/cts.h index 5a7ef11..6027532 100644 --- a/src/cts.h +++ b/src/cts.h @@ -1,47 +1,51 @@ #ifndef _cts_h_ #define _cts_h_ -#include +#include #include -#include +#include #include #include -#include +#include #include "templates.h" // [[Rcpp::plugins("cpp11")]] // [[Rcpp::depends("RcppArmadillo")]] template -static void combn(T2& vals, const int32_t n, const uint32_t start_idx, - std::vector& combn_data, T1& combn_ds, uint32_t& combn_col); +static void combn(T2& vals, const int32_t n, const uint32_t start_idx, + std::vector& combn_data, T1& combn_ds, + uint32_t& combn_col); template T1 find_combn(T2& vals, const int32_t n) { - static uint32_t combn_col = 0; - const uint32_t nrows = n; - const uint32_t ncols = std::round(R::choose(vals.size(), n)); - T1 combn_ds(nrows, ncols); - std::vector combn_data(nrows); - const uint32_t start_idx = 0; - combn_col = 0; combn(vals, n, start_idx, combn_data, combn_ds, combn_col); - return combn_ds; + static uint32_t combn_col = 0; + const uint32_t nrows = n; + const uint32_t ncols = std::round(R::choose(vals.size(), n)); + T1 combn_ds(nrows, ncols); + std::vector combn_data(nrows); + const uint32_t start_idx = 0; + combn_col = 0; + combn(vals, n, start_idx, combn_data, combn_ds, combn_col); + return combn_ds; } template -static void combn(T2& vals, const int32_t n, const uint32_t start_idx, - std::vector& combn_data, T1& combn_ds, uint32_t& combn_col) { - if (!n) { - for (uint32_t i = 0; i < combn_ds.n_rows && combn_col < combn_ds.n_cols; i++) { - combn_ds(i, combn_col) = combn_data.at(i); - } - combn_col++; - return; - } - for (uint32_t i = start_idx; i <= (vals.size() - n); i++) { - combn_data.at(combn_ds.n_rows - n) = vals[i]; - combn(vals, n - 1, i + 1, combn_data, combn_ds, combn_col); - } +static void combn(T2& vals, const int32_t n, const uint32_t start_idx, + std::vector& combn_data, T1& combn_ds, + uint32_t& combn_col) { + if (!n) { + for (uint32_t i = 0; i < combn_ds.n_rows && combn_col < combn_ds.n_cols; + i++) { + combn_ds(i, combn_col) = combn_data.at(i); + } + combn_col++; + return; + } + for (uint32_t i = start_idx; i <= (vals.size() - n); i++) { + combn_data.at(combn_ds.n_rows - n) = vals[i]; + combn(vals, n - 1, i + 1, combn_data, combn_ds, combn_col); + } } // Alters ds @@ -54,14 +58,12 @@ arma::mat calc_rank(arma::mat& ds); arma::uvec sub_col_max_min(arma::mat& ds, const bool cont); -arma::mat calc_pt(arma::mat& ds, - const int32_t df, const bool lower_tail, const bool log_p, - const double add); +arma::mat calc_pt(arma::mat& ds, const int32_t df, const bool lower_tail, + const bool log_p, const double add); std::vector det_cols(arma::umat& ds, const uint32_t val); -arma::mat ext_cols(arma::mat& ds, const uint32_t col_a, - const uint32_t col_b); +arma::mat ext_cols(arma::mat& ds, const uint32_t col_a, const uint32_t col_b); // Alters dst void cp_lt(arma::mat& src, arma::mat& dst, const int32_t val); @@ -71,19 +73,20 @@ void adj_diag(arma::mat& ds, const double val); arma::mat cbind_tran_mat(arma::mat& ds1, arma::mat& ds2); -arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, - const bool ass1, const bool ass2); +arma::mat rbind_uniq(arma::mat& ds1, arma::mat& ds2, const bool ass1, + const bool ass2); arma::vec to_vec(arma::mat& ds); bool found_match(const uint32_t x, arma::uvec& vals); -bool are_equal(arma::mat& ds, arma::vec& vals, - const bool by_col_idx = false, const unsigned col = 0); +bool are_equal(arma::mat& ds, arma::vec& vals, const bool by_col_idx = false, + const unsigned col = 0); std::vector rm_lt_nan(arma::uvec& idxs, const uint32_t limit); -std::unordered_set get_diff(std::vector& lh, const uint32_t val); +std::unordered_set get_diff(std::vector& lh, + const uint32_t val); arma::umat nchoosek(std::vector& idxs, const uint32_t k); @@ -103,8 +106,8 @@ arma::mat form_cmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols); arma::mat form_rmat(arma::mat& ds, arma::uvec& rows, arma::uvec& cols); -arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, - std::vector& cols); +arma::mat form_rmat_std(arma::mat& ds, std::vector& rows, + std::vector& cols); arma::mat merge_cols(arma::mat& ds, arma::uvec& idxs); @@ -115,7 +118,7 @@ arma::mat form_c2mat(arma::vec& vals1, arma::vec& vals2); arma::uvec form_vec(arma::mat& ds, const uint32_t row, arma::uvec& cols); arma::vec form_vec_wvals(arma::mat& ds, const uint32_t row, arma::uvec& cols, - arma::vec& vals); + arma::vec& vals); // Alters ds arma::mat append_row(arma::mat& ds, const uint32_t row, arma::vec& vals); diff --git a/src/mmp_c.cpp b/src/mmp_c.cpp index 1477c62..4083fc9 100644 --- a/src/mmp_c.cpp +++ b/src/mmp_c.cpp @@ -9,12 +9,16 @@ #include "mmp_c.h" #define DEBUG 0 -#define db_print(...) \ - do { if (DEBUG) Rprintf(__VA_ARGS__); } while (0) +#define db_print(...) \ + do { \ + if (DEBUG) Rprintf(__VA_ARGS__); \ + } while (0) #define DEBUG_L 0 -#define dbl_print(...) \ - do { if (DEBUG_L) Rprintf(__VA_ARGS__); } while (0) +#define dbl_print(...) \ + do { \ + if (DEBUG_L) Rprintf(__VA_ARGS__); \ + } while (0) static arma::vec DEF_VEC; static arma::uvec DEF_UVEC; @@ -28,591 +32,611 @@ static double STAT_PVALUE[2]; Rcpp::List res; uint32_t kv_length; -static void check_args(const double thres, const int32_t max_k, const std::string method) { +static void check_args(const double thres, const int32_t max_k, + const std::string method) { if (max_k < 1) { Rcpp::stop("Invalid max_k argument provided.\nExiting...\n"); - } - else if (thres < 0 || thres >= 1) { + } else if (thres < 0 || thres >= 1) { Rcpp::stop("Invalid thres argument provided.\nExiting...\n"); - } - else if (method.compare("pearson") && method.compare("spearman")) { + } else if (method.compare("pearson") && method.compare("spearman")) { Rcpp::stop("Invalid method name provided.\nExiting...\n"); } } -static void upd_args(arma::vec& target_vars, arma::mat& ds, const std::string method) { +static void upd_args(arma::vec& target_vars, arma::mat& ds, + const std::string method) { if (!method.compare("spearman")) { - target_vars = rank_mean(target_vars, false); + target_vars = + rank_mean(target_vars, false); ds = calc_rank(ds); } } -static arma::mat calc_resid(arma::mat& cor_ds, arma::mat& sol_ds, - arma::uvec& idxs_a, arma::uvec& idxs_b) { - arma::mat lh_sub = cor_ds.submat(idxs_a, idxs_a); - arma::mat rh_sub = cor_ds.submat(idxs_a, idxs_b); - arma::mat mult_res = rh_sub * sol_ds; - arma::mat sub_res = lh_sub - mult_res; - return sub_res; +static arma::mat calc_resid(arma::mat& cor_ds, arma::mat& sol_ds, + arma::uvec& idxs_a, arma::uvec& idxs_b) { + arma::mat lh_sub = cor_ds.submat(idxs_a, idxs_a); + arma::mat rh_sub = cor_ds.submat(idxs_a, idxs_b); + arma::mat mult_res = rh_sub * sol_ds; + arma::mat sub_res = lh_sub - mult_res; + return sub_res; } -static arma::mat calc_sol(arma::mat& ds, arma::uvec& idxs_a, arma::uvec& idxs_b) { - arma::mat lh = ds.submat(idxs_b, idxs_b); - arma::mat rh = ds.submat(idxs_b, idxs_a); - arma::mat res = arma::solve(lh, rh); - return res; +static arma::mat calc_sol(arma::mat& ds, arma::uvec& idxs_a, + arma::uvec& idxs_b) { + arma::mat lh = ds.submat(idxs_b, idxs_b); + arma::mat rh = ds.submat(idxs_b, idxs_a); + arma::mat res = arma::solve(lh, rh); + return res; } static void upd_col(arma::vec& src, arma::mat& dst, const uint32_t col) { - for (uint32_t i = 0; i < src.size(); ++i) { - dst(i, col) = src[i]; - } + for (uint32_t i = 0; i < src.size(); ++i) { + dst(i, col) = src[i]; + } } -static arma::mat col_bind(arma::vec& src_a, arma::vec& src_b, arma::mat& src_c) { - const uint32_t nrows = src_c.n_rows; - const uint32_t ncols = src_c.n_cols + 2; - arma::mat dst(nrows, ncols); - upd_col(src_a, dst, 0); - upd_col(src_b, dst, 1); - for (uint32_t src_col = 0, dst_col = 2; - src_col < src_c.n_cols && dst_col < dst.n_cols; - ++src_col, ++dst_col) { - for (uint32_t i = 0; i < nrows; ++i) { - dst(i, dst_col) = src_c(i, src_col); - } - } - return dst; +static arma::mat col_bind(arma::vec& src_a, arma::vec& src_b, + arma::mat& src_c) { + const uint32_t nrows = src_c.n_rows; + const uint32_t ncols = src_c.n_cols + 2; + arma::mat dst(nrows, ncols); + upd_col(src_a, dst, 0); + upd_col(src_b, dst, 1); + for (uint32_t src_col = 0, dst_col = 2; + src_col < src_c.n_cols && dst_col < dst.n_cols; ++src_col, ++dst_col) { + for (uint32_t i = 0; i < nrows; ++i) { + dst(i, dst_col) = src_c(i, src_col); + } + } + return dst; } - static void form_ret(const double stat, const double pvalue) { - STAT_PVALUE[STAT_POS] = stat; - STAT_PVALUE[PVALUE_POS] = pvalue; + STAT_PVALUE[STAT_POS] = stat; + STAT_PVALUE[PVALUE_POS] = pvalue; } - static std::string form_key(const uint32_t val, std::vector& vals) { - std::string key = std::to_string(val + 1); - for (uint32_t i = 0; i < vals.size(); ++i) { - key += " " + std::to_string(vals[i] + 1); - } - return key; + std::string key = std::to_string(val + 1); + for (uint32_t i = 0; i < vals.size(); ++i) { + key += " " + std::to_string(vals[i] + 1); + } + return key; } -void calc_pearson(arma::vec& target_vars, arma::mat& ds, - const uint32_t x, arma::uvec& idxs, Rcpp::List& univs, - const bool hash_on, Rcpp::Environment& stats_kv, - Rcpp::Environment& pvalues_kv) { - const double def_stat = 0; - const double def_pvalue = std::log(1); - db_print("Call: rm_lt_nan\n"); - std::vector idxs_adj = rm_lt_nan(idxs, 0); - std::string key; - if (hash_on) { - db_print("True: hash_on\n"); - std::vector idxs_nz(idxs_adj); - db_print("Call: std::sort\n"); - std::sort(idxs_nz.begin(), idxs_nz.end()); - db_print("Call: form_key\n"); - key = form_key(x, idxs_nz); - if (stats_kv.exists(key)) { - db_print("True: stats_kv.exists(key)\n"); - form_ret(stats_kv[key], pvalues_kv[key]); - return; - } - } - db_print("Call: found_match\n"); - if (found_match(x, idxs)) { - db_print("True: found_match(x, idxs)\n"); - if (hash_on) { - db_print("True: hash_on\n"); - if (!stats_kv.exists(key)) { - ++kv_length; - } - stats_kv[key] = def_stat; - pvalues_kv[key] = def_pvalue; - } - form_ret(def_stat, def_pvalue); - return; - } - db_print("Call: arma::unique\n"); - idxs = arma::unique(idxs); - arma::mat tmp_ds = ds.cols(idxs); - arma::vec tmp_vec = ds.col(x); - if (!tmp_ds.is_empty()) { - db_print("True: !tmp_ds.is_empty()\n"); - if (tmp_ds.n_cols == 1) { - db_print("True: tmp_ds.n_cols == 1\n"); - if (are_equal(tmp_ds, tmp_vec)) { - db_print("True: are_equal(tmp_ds, tmp_vec)\n"); - if (hash_on) { - db_print("True: hash_on\n"); - if (!stats_kv.exists(key)) { - ++kv_length; - } - stats_kv[key] = def_stat; - pvalues_kv[key] = def_pvalue; - } - form_ret(def_stat, def_pvalue); - return; - } - } - else { - db_print("True: tmp_ds.n_cols != 1\n"); - db_print("Checking tmp_ds, tmp_vec by column equality.\n"); - for (uint32_t j = 0; j < tmp_ds.n_cols; ++j) { - if (are_equal(tmp_ds, tmp_vec, true, j)) { - if (hash_on) { - if (!stats_kv.exists(key)) { - ++kv_length; - } - stats_kv[key] = def_stat; - pvalues_kv[key] = def_pvalue; - } - form_ret(def_stat, def_pvalue); - return; - } - } - } - } - double stat = 0; - double pvalue = 0; - if (tmp_ds.is_empty()) { - db_print("True: tmp_ds.is_empty()\n"); - if (univs.size()) { - db_print("True: univs\n"); - arma::mat tmp_stats = univs["stat"]; - arma::mat tmp_pvalues = univs["pvalue"]; - form_ret(tmp_stats[x], tmp_pvalues[x]); - return; - } - db_print("Call: arma::cor\n"); - arma::mat tmp = arma::cor(tmp_vec, target_vars); - stat = tmp(0, 0); - } - else { - db_print("True: !tmp_ds.is_empty()\n"); - db_print("Call: col_bind\n"); - arma::mat tmp_bind = col_bind(tmp_vec, target_vars, tmp_ds); - db_print("Call: arma::cor\n"); - arma::mat cor_ds = arma::cor(tmp_bind); - db_print("Forming indices.\n"); - arma::uvec idxs_a(2); - std::iota(idxs_a.begin(), idxs_a.end(), 0); - arma::uvec idxs_b(idxs.size()); - for (uint32_t i = 0; i < idxs_b.size(); ++i) { - idxs_b[i] = i + 2; - } - db_print("Call: calc_sol\n"); - arma::mat sol_ds = calc_sol(cor_ds, idxs_a, idxs_b); - db_print("Call: calc_resid\n"); - arma::mat resid_ds = calc_resid(cor_ds, sol_ds, idxs_a, idxs_b); - stat = resid_ds(0, 1) / std::sqrt(resid_ds(0, 0) * resid_ds(1, 1)); - } - db_print("Forming final stat, pvalue.\n"); - const double lh = 0.5 * std::log((1 + stat) / (1 - stat)); - const double rh = target_vars.size() - idxs.size() - 3; - stat = std::sqrt(rh) * std::abs(lh); - pvalue = std::log(2) + R::pt(stat, rh, false, true); - if (hash_on) { - db_print("True: hash_on\n"); - if (!stats_kv.exists(key)) { - ++kv_length; - } - stats_kv[key] = stat; - pvalues_kv[key] = pvalue; - } - form_ret(stat, pvalue); +void calc_pearson(arma::vec& target_vars, arma::mat& ds, const uint32_t x, + arma::uvec& idxs, Rcpp::List& univs, const bool hash_on, + Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { + const double def_stat = 0; + const double def_pvalue = std::log(1); + db_print("Call: rm_lt_nan\n"); + std::vector idxs_adj = rm_lt_nan(idxs, 0); + std::string key; + if (hash_on) { + db_print("True: hash_on\n"); + std::vector idxs_nz(idxs_adj); + db_print("Call: std::sort\n"); + std::sort(idxs_nz.begin(), idxs_nz.end()); + db_print("Call: form_key\n"); + key = form_key(x, idxs_nz); + if (stats_kv.exists(key)) { + db_print("True: stats_kv.exists(key)\n"); + form_ret(stats_kv[key], pvalues_kv[key]); + return; + } + } + db_print("Call: found_match\n"); + if (found_match(x, idxs)) { + db_print("True: found_match(x, idxs)\n"); + if (hash_on) { + db_print("True: hash_on\n"); + if (!stats_kv.exists(key)) { + ++kv_length; + } + stats_kv[key] = def_stat; + pvalues_kv[key] = def_pvalue; + } + form_ret(def_stat, def_pvalue); + return; + } + db_print("Call: arma::unique\n"); + idxs = arma::unique(idxs); + arma::mat tmp_ds = ds.cols(idxs); + arma::vec tmp_vec = ds.col(x); + if (!tmp_ds.is_empty()) { + db_print("True: !tmp_ds.is_empty()\n"); + if (tmp_ds.n_cols == 1) { + db_print("True: tmp_ds.n_cols == 1\n"); + if (are_equal(tmp_ds, tmp_vec)) { + db_print("True: are_equal(tmp_ds, tmp_vec)\n"); + if (hash_on) { + db_print("True: hash_on\n"); + if (!stats_kv.exists(key)) { + ++kv_length; + } + stats_kv[key] = def_stat; + pvalues_kv[key] = def_pvalue; + } + form_ret(def_stat, def_pvalue); + return; + } + } else { + db_print("True: tmp_ds.n_cols != 1\n"); + db_print("Checking tmp_ds, tmp_vec by column equality.\n"); + for (uint32_t j = 0; j < tmp_ds.n_cols; ++j) { + if (are_equal(tmp_ds, tmp_vec, true, j)) { + if (hash_on) { + if (!stats_kv.exists(key)) { + ++kv_length; + } + stats_kv[key] = def_stat; + pvalues_kv[key] = def_pvalue; + } + form_ret(def_stat, def_pvalue); + return; + } + } + } + } + double stat = 0; + double pvalue = 0; + if (tmp_ds.is_empty()) { + db_print("True: tmp_ds.is_empty()\n"); + if (univs.size()) { + db_print("True: univs\n"); + arma::mat tmp_stats = univs["stat"]; + arma::mat tmp_pvalues = univs["pvalue"]; + form_ret(tmp_stats[x], tmp_pvalues[x]); + return; + } + db_print("Call: arma::cor\n"); + arma::mat tmp = arma::cor(tmp_vec, target_vars); + stat = tmp(0, 0); + } else { + db_print("True: !tmp_ds.is_empty()\n"); + db_print("Call: col_bind\n"); + arma::mat tmp_bind = col_bind(tmp_vec, target_vars, tmp_ds); + db_print("Call: arma::cor\n"); + arma::mat cor_ds = arma::cor(tmp_bind); + db_print("Forming indices.\n"); + arma::uvec idxs_a(2); + std::iota(idxs_a.begin(), idxs_a.end(), 0); + arma::uvec idxs_b(idxs.size()); + for (uint32_t i = 0; i < idxs_b.size(); ++i) { + idxs_b[i] = i + 2; + } + db_print("Call: calc_sol\n"); + arma::mat sol_ds = calc_sol(cor_ds, idxs_a, idxs_b); + db_print("Call: calc_resid\n"); + arma::mat resid_ds = calc_resid(cor_ds, sol_ds, idxs_a, idxs_b); + stat = resid_ds(0, 1) / std::sqrt(resid_ds(0, 0) * resid_ds(1, 1)); + } + db_print("Forming final stat, pvalue.\n"); + const double lh = 0.5 * std::log((1 + stat) / (1 - stat)); + const double rh = target_vars.size() - idxs.size() - 3; + stat = std::sqrt(rh) * std::abs(lh); + pvalue = std::log(2) + R::pt(stat, rh, false, true); + if (hash_on) { + db_print("True: hash_on\n"); + if (!stats_kv.exists(key)) { + ++kv_length; + } + stats_kv[key] = stat; + pvalues_kv[key] = pvalue; + } + form_ret(stat, pvalue); } -bool cmp_pvalues(const double stat_lh, const double stat_rh, - const double pvalue_lh, const double pvalue_rh) { - if (!arma::is_finite(stat_lh) || !arma::is_finite(stat_rh) || - !arma::is_finite(pvalue_lh) || !arma::is_finite(pvalue_rh)) { - return false; - } - if (pvalue_lh == pvalue_rh) { - return stat_lh > stat_rh; - } - else { - return pvalue_lh < pvalue_rh; - } +bool cmp_pvalues(const double stat_lh, const double stat_rh, + const double pvalue_lh, const double pvalue_rh) { + if (!arma::is_finite(stat_lh) || !arma::is_finite(stat_rh) || + !arma::is_finite(pvalue_lh) || !arma::is_finite(pvalue_rh)) { + return false; + } + if (pvalue_lh == pvalue_rh) { + return stat_lh > stat_rh; + } else { + return pvalue_lh < pvalue_rh; + } } // Alters ds static void rbind_mat(arma::umat& ds, const uint32_t val) { - ds.resize(ds.n_rows + 1, ds.n_cols); - for (uint32_t j = 0; j < ds.n_cols; ++j) { - ds(ds.n_rows - 1, j) = val; - } + ds.resize(ds.n_rows + 1, ds.n_cols); + for (uint32_t j = 0; j < ds.n_cols; ++j) { + ds(ds.n_rows - 1, j) = val; + } } -static std::vector keep_val(arma::uvec& vals, const uint32_t val_to_keep) { - std::vector vals_adj; - for (uint32_t i = 0; i < vals.size(); ++i) { - if (vals[i] == val_to_keep) { - vals_adj.push_back(i); - } - } - return vals_adj; +static std::vector keep_val(arma::uvec& vals, + const uint32_t val_to_keep) { + std::vector vals_adj; + for (uint32_t i = 0; i < vals.size(); ++i) { + if (vals[i] == val_to_keep) { + vals_adj.push_back(i); + } + } + return vals_adj; } -void assoc_min(arma::vec& target_vars, arma::mat& ds, const std::string method, - const int32_t max_k, const uint32_t sel_idx, arma::uvec& sel_idxs, arma::vec& stats, - arma::vec& pvalues, arma::uvec& sel_ord_idxs, const bool hash_on, - Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { - double sel_stat = stats[sel_idx]; - double sel_pvalue = pvalues[sel_idx]; - std::vector sel_idxs_adj = keep_val(sel_idxs, 1); - const uint32_t k = std::min(max_k, (int) sel_idxs_adj.size()); - uint32_t curr_k = 1; - while (curr_k <= k) { - const uint32_t max_idx = arma::index_max(sel_ord_idxs); - std::unordered_set tmp_idxs = get_diff(sel_idxs_adj, max_idx); - std::vector tmp_idxs_adj(tmp_idxs.begin(), tmp_idxs.end()); - arma::umat subsetcsk; - if (curr_k == 1) { - subsetcsk = arma::umat(1, 1); - subsetcsk(0, 0) = max_idx; - } - else { - subsetcsk = nchoosek(tmp_idxs_adj, curr_k - 1); - rbind_mat(subsetcsk, max_idx); - } - for (uint32_t i = 0; i < subsetcsk.n_cols; ++i) { - arma::uvec subsetcsk_col = subsetcsk.col(i); - Rcpp::List univs; - calc_pearson(target_vars, ds, sel_idx, subsetcsk_col, univs, hash_on, stats_kv, pvalues_kv); - if (!cmp_pvalues(STAT_PVALUE[STAT_POS], sel_stat, STAT_PVALUE[PVALUE_POS], sel_pvalue)) { - sel_stat = STAT_PVALUE[STAT_POS]; - stats[sel_idx] = sel_stat; - sel_pvalue = STAT_PVALUE[PVALUE_POS]; - pvalues[sel_idx] = sel_pvalue; - } - } - ++curr_k; - } - form_ret(sel_stat, sel_pvalue); +void assoc_min(arma::vec& target_vars, arma::mat& ds, const std::string method, + const int32_t max_k, const uint32_t sel_idx, + arma::uvec& sel_idxs, arma::vec& stats, arma::vec& pvalues, + arma::uvec& sel_ord_idxs, const bool hash_on, + Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { + double sel_stat = stats[sel_idx]; + double sel_pvalue = pvalues[sel_idx]; + std::vector sel_idxs_adj = keep_val(sel_idxs, 1); + const uint32_t k = std::min(max_k, (int)sel_idxs_adj.size()); + uint32_t curr_k = 1; + while (curr_k <= k) { + const uint32_t max_idx = arma::index_max(sel_ord_idxs); + std::unordered_set tmp_idxs = get_diff(sel_idxs_adj, max_idx); + std::vector tmp_idxs_adj(tmp_idxs.begin(), tmp_idxs.end()); + arma::umat subsetcsk; + if (curr_k == 1) { + subsetcsk = arma::umat(1, 1); + subsetcsk(0, 0) = max_idx; + } else { + subsetcsk = nchoosek(tmp_idxs_adj, curr_k - 1); + rbind_mat(subsetcsk, max_idx); + } + for (uint32_t i = 0; i < subsetcsk.n_cols; ++i) { + arma::uvec subsetcsk_col = subsetcsk.col(i); + Rcpp::List univs; + calc_pearson(target_vars, ds, sel_idx, subsetcsk_col, univs, + hash_on, stats_kv, pvalues_kv); + if (!cmp_pvalues(STAT_PVALUE[STAT_POS], sel_stat, + STAT_PVALUE[PVALUE_POS], sel_pvalue)) { + sel_stat = STAT_PVALUE[STAT_POS]; + stats[sel_idx] = sel_stat; + sel_pvalue = STAT_PVALUE[PVALUE_POS]; + pvalues[sel_idx] = sel_pvalue; + } + } + ++curr_k; + } + form_ret(sel_stat, sel_pvalue); } -uint32_t assoc_max_min(arma::vec& target_vars, arma::mat& ds, const std::string method, - const double thres, const int32_t max_k, arma::uvec& sel_idxs, - arma::vec& stats, arma::vec& pvalues, arma::uvec& rem_idxs, arma::uvec& sel_ord_idxs, - const bool hash_on, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { - int32_t sel_idx = -1; - double sel_pvalue = 2; - double sel_stat = 0; - std::vector rem_idxs_adj = keep_val(rem_idxs, 1); - for (uint32_t i = 0; i < rem_idxs_adj.size(); ++i) { - assoc_min(target_vars, ds, method, max_k, rem_idxs_adj[i], sel_idxs, - stats, pvalues, sel_ord_idxs, hash_on, stats_kv, pvalues_kv); - const double tmp_pvalue = STAT_PVALUE[PVALUE_POS]; - if (tmp_pvalue > thres) { - rem_idxs[rem_idxs_adj[i]] = 0; - } - if (cmp_pvalues(STAT_PVALUE[STAT_POS], sel_stat, STAT_PVALUE[PVALUE_POS], sel_pvalue)) { - sel_idx = rem_idxs_adj[i]; - sel_stat = STAT_PVALUE[STAT_POS]; - sel_pvalue = tmp_pvalue; - } - } - form_ret(sel_stat, sel_pvalue); - return sel_idx; +uint32_t assoc_max_min(arma::vec& target_vars, arma::mat& ds, + const std::string method, const double thres, + const int32_t max_k, arma::uvec& sel_idxs, + arma::vec& stats, arma::vec& pvalues, + arma::uvec& rem_idxs, arma::uvec& sel_ord_idxs, + const bool hash_on, Rcpp::Environment& stats_kv, + Rcpp::Environment& pvalues_kv) { + int32_t sel_idx = -1; + double sel_pvalue = 2; + double sel_stat = 0; + std::vector rem_idxs_adj = keep_val(rem_idxs, 1); + for (uint32_t i = 0; i < rem_idxs_adj.size(); ++i) { + assoc_min(target_vars, ds, method, max_k, rem_idxs_adj[i], sel_idxs, + stats, pvalues, sel_ord_idxs, hash_on, stats_kv, pvalues_kv); + const double tmp_pvalue = STAT_PVALUE[PVALUE_POS]; + if (tmp_pvalue > thres) { + rem_idxs[rem_idxs_adj[i]] = 0; + } + if (cmp_pvalues(STAT_PVALUE[STAT_POS], sel_stat, + STAT_PVALUE[PVALUE_POS], sel_pvalue)) { + sel_idx = rem_idxs_adj[i]; + sel_stat = STAT_PVALUE[STAT_POS]; + sel_pvalue = tmp_pvalue; + } + } + form_ret(sel_stat, sel_pvalue); + return sel_idx; } static arma::uvec get_sort_idxs(arma::vec& vals, arma::uvec& idxs) { - std::vector vals_tmp; - for (uint32_t i = 0; i < idxs.size(); ++i) { - vals_tmp.push_back(vals[idxs[i]]); - } - arma::vec vals_adj(vals_tmp); - return arma::sort_index(vals_adj); + std::vector vals_tmp; + for (uint32_t i = 0; i < idxs.size(); ++i) { + vals_tmp.push_back(vals[idxs[i]]); + } + arma::vec vals_adj(vals_tmp); + return arma::sort_index(vals_adj); } static arma::uvec get_idxs_eq(arma::uvec vals, const uint32_t val) { - std::vector vals_adj; - for (uint32_t i = 0; i < vals.size(); ++i) { - if (vals[i] == val) { - vals_adj.push_back(i); - } - } - return arma::uvec(vals_adj); + std::vector vals_adj; + for (uint32_t i = 0; i < vals.size(); ++i) { + if (vals[i] == val) { + vals_adj.push_back(i); + } + } + return arma::uvec(vals_adj); } -static void neg_zero_idxs(arma::uvec& src, arma::uvec& base, const uint32_t val) { - for (uint32_t i = 0; i < src.size() && i < base.size(); ++i) { - if (!base[i]) { - src[i] = val; - } - } +static void neg_zero_idxs(arma::uvec& src, arma::uvec& base, + const uint32_t val) { + for (uint32_t i = 0; i < src.size() && i < base.size(); ++i) { + if (!base[i]) { + src[i] = val; + } + } } static bool is_true(arma::uvec& vals) { - for (uint32_t i = 0; i < vals.size(); ++i) { - if (vals[i]) { - return true; - } - } - return false; + for (uint32_t i = 0; i < vals.size(); ++i) { + if (vals[i]) { + return true; + } + } + return false; } -static void adj_gt_vals(arma::uvec& rem_idxs, arma::vec& pvalues, const double thres, const double val) { - for (uint32_t i = 0; i < rem_idxs.size(); ++i) { - if (pvalues[i] > thres) { - rem_idxs[i] = 0; - } - } +static void adj_gt_vals(arma::uvec& rem_idxs, arma::vec& pvalues, + const double thres, const double val) { + for (uint32_t i = 0; i < rem_idxs.size(); ++i) { + if (pvalues[i] > thres) { + rem_idxs[i] = 0; + } + } } static void copy_vecs(Rcpp::List& src, arma::vec& dst1, arma::vec& dst2) { - arma::vec src1 = src["stat"]; - arma::vec src2 = src["pvalue"]; - dst1.set_size(src1.size()); - dst2.set_size(src2.size()); - - for (uint32_t i = 0; i < src1.size(); ++i) { - dst1[i] = src1[i]; - dst2[i] = src2[i]; - } + arma::vec src1 = src["stat"]; + arma::vec src2 = src["pvalue"]; + dst1.set_size(src1.size()); + dst2.set_size(src2.size()); + + for (uint32_t i = 0; i < src1.size(); ++i) { + dst1[i] = src1[i]; + dst2[i] = src2[i]; + } } static arma::vec calc_univ_pvalues(arma::vec& stats, const double dof) { - arma::vec pvalues(stats.size()); - static const double tmp_log = std::log(2); - for (uint32_t i = 0; i < stats.size(); ++i) { - pvalues[i] = tmp_log + R::pt(std::abs(stats[i]), dof, false, true); - } - return pvalues; + arma::vec pvalues(stats.size()); + static const double tmp_log = std::log(2); + for (uint32_t i = 0; i < stats.size(); ++i) { + pvalues[i] = tmp_log + R::pt(std::abs(stats[i]), dof, false, true); + } + return pvalues; } -static void calc_univs(arma::vec& target_vars, arma::mat& ds, const std::string method, Rcpp::List& univs) { +static void calc_univs(arma::vec& target_vars, arma::mat& ds, + const std::string method, Rcpp::List& univs) { arma::mat cor_mat = arma::cor(target_vars, ds); arma::vec cor_vec = to_vec(cor_mat); const uint32_t dof = ds.n_rows - 3; arma::vec stats; - if (!method.compare("pearson")) { - db_print("True: !method.compare(\"pearson\")\n"); - stats = 0.5 * arma::log((1 + cor_vec) / (1 - cor_vec)) * std::sqrt(dof); - } - else if (!method.compare("spearman")) { - db_print("True: !method.compare(\"spearman\")\n"); - stats = 0.5 * arma::log((1 + cor_vec) / (1 - cor_vec)) * std::sqrt(dof) / 1.029563; + if (!method.compare("pearson")) { + db_print("True: !method.compare(\"pearson\")\n"); + stats = 0.5 * arma::log((1 + cor_vec) / (1 - cor_vec)) * std::sqrt(dof); + } else if (!method.compare("spearman")) { + db_print("True: !method.compare(\"spearman\")\n"); + stats = 0.5 * arma::log((1 + cor_vec) / (1 - cor_vec)) * + std::sqrt(dof) / 1.029563; } univs["stat"] = stats; - univs["pvalue"] = calc_univ_pvalues(stats, dof); + univs["pvalue"] = calc_univ_pvalues(stats, dof); } static double get_time_spent(const clock_t begin) { - const clock_t end = clock(); - return (double) (end - begin) / CLOCKS_PER_SEC; + const clock_t end = clock(); + return (double)(end - begin) / CLOCKS_PER_SEC; } -void inter_mmp_c(arma::vec& target_vars, arma::mat& ds, const int32_t max_k, - const double thres, const std::string method, Rcpp::List inits, const bool hash_on, - const uint32_t var_size, Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { - db_print("Initializing clock.\n"); - clock_t begin = clock(); - - db_print("Declaring vars.\n"); - Rcpp::List univs; - arma::vec stats; - arma::vec pvalues; - - if (!inits.size()) { - db_print("True: !inits\n"); - calc_univs(target_vars, ds, method, univs); - } - else { - db_print("True: inits.size()\n"); - univs = inits; - } - if (univs.size()) { - db_print("True: univs.size()\n"); - db_print("Call: copy_vec\n"); - copy_vecs(univs, stats, pvalues); - } - db_print("Call: arma::min\n"); - const double min = arma::min(pvalues); - if (min > thres) { - db_print("True: min > thres\n"); - res["selected"] = R_NilValue; - res["hashobject"] = R_NilValue; - res["stats"] = stats; - res["pvalues"] = pvalues; - res["univ"] = univs; - res["max_k"] = max_k; - res["alpha"] = thres; - res["runtime"] = get_time_spent(begin); - inits.size() ? res["n.tests"] = 0 : res["n.tests"] = stats.size(); +void inter_mmp_c(arma::vec& target_vars, arma::mat& ds, const int32_t max_k, + const double thres, const std::string method, Rcpp::List inits, + const bool hash_on, const uint32_t var_size, + Rcpp::Environment& stats_kv, Rcpp::Environment& pvalues_kv) { + db_print("Initializing clock.\n"); + clock_t begin = clock(); + + db_print("Declaring vars.\n"); + Rcpp::List univs; + arma::vec stats; + arma::vec pvalues; + + if (!inits.size()) { + db_print("True: !inits\n"); + calc_univs(target_vars, ds, method, univs); + } else { + db_print("True: inits.size()\n"); + univs = inits; + } + if (univs.size()) { + db_print("True: univs.size()\n"); + db_print("Call: copy_vec\n"); + copy_vecs(univs, stats, pvalues); + } + db_print("Call: arma::min\n"); + const double min = arma::min(pvalues); + if (min > thres) { + db_print("True: min > thres\n"); + res["selected"] = R_NilValue; + res["hashobject"] = R_NilValue; + res["stats"] = stats; + res["pvalues"] = pvalues; + res["univ"] = univs; + res["max_k"] = max_k; + res["alpha"] = thres; + res["runtime"] = get_time_spent(begin); + inits.size() ? res["n.tests"] = 0 : res["n.tests"] = stats.size(); return; - } - db_print("Initializing idxs.\n"); - arma::uvec sel_idxs(var_size, arma::fill::zeros); - arma::uvec sel_ord_idxs(var_size, arma::fill::zeros); - uint32_t sel_idx = arma::index_min(pvalues); - sel_idxs[sel_idx] = 1; - sel_ord_idxs[sel_idx] = 1; - - arma::uvec rem_idxs(var_size, arma::fill::ones); - rem_idxs[sel_idx] = 0; - db_print("Call: adj_gt_vals\n"); - adj_gt_vals(rem_idxs, pvalues, thres, 0); - db_print("Entering main loop.\n"); - while (is_true(rem_idxs)) { - sel_idx = assoc_max_min(target_vars, ds, method, thres, max_k, - sel_idxs, stats, pvalues, rem_idxs, sel_ord_idxs, hash_on, stats_kv, pvalues_kv); - const double sel_pvalue = STAT_PVALUE[PVALUE_POS]; - if (sel_pvalue <= thres) { - sel_idxs[sel_idx] = 1; - sel_ord_idxs[sel_idx] = arma::max(sel_ord_idxs) + 1; - rem_idxs[sel_idx] = 0; - } - } - db_print("Call: neg_zero_idxs\n"); - neg_zero_idxs(sel_ord_idxs, sel_idxs, var_size); - db_print("Call: arma::sort\n"); - sel_ord_idxs = arma::sort(sel_ord_idxs); - - db_print("Call: get_idxs_eq\n"); - arma::uvec tmp_sel_idxs = get_idxs_eq(sel_idxs, 1); - db_print("Call: get_sort_idxs\n"); - arma::uvec tmp_sort_idxs = get_sort_idxs(pvalues, tmp_sel_idxs); - arma::uvec tmp_sel_ord_idxs = tmp_sel_idxs(tmp_sort_idxs); - db_print("Forming return list.\n"); - res["selected"] = tmp_sel_ord_idxs; - res["stats"] = stats; - res["pvalues"] = pvalues; - res["univ"] = univs; - res["max_k"] = max_k; - res["alpha"] = thres; - res["runtime"] = get_time_spent(begin); - inits.size() ? res["n.tests"] = 0 : res["n.tests"] = stats.n_rows * stats.n_cols + kv_length; + } + db_print("Initializing idxs.\n"); + arma::uvec sel_idxs(var_size, arma::fill::zeros); + arma::uvec sel_ord_idxs(var_size, arma::fill::zeros); + uint32_t sel_idx = arma::index_min(pvalues); + sel_idxs[sel_idx] = 1; + sel_ord_idxs[sel_idx] = 1; + + arma::uvec rem_idxs(var_size, arma::fill::ones); + rem_idxs[sel_idx] = 0; + db_print("Call: adj_gt_vals\n"); + adj_gt_vals(rem_idxs, pvalues, thres, 0); + db_print("Entering main loop.\n"); + while (is_true(rem_idxs)) { + sel_idx = assoc_max_min(target_vars, ds, method, thres, max_k, sel_idxs, + stats, pvalues, rem_idxs, sel_ord_idxs, hash_on, + stats_kv, pvalues_kv); + const double sel_pvalue = STAT_PVALUE[PVALUE_POS]; + if (sel_pvalue <= thres) { + sel_idxs[sel_idx] = 1; + sel_ord_idxs[sel_idx] = arma::max(sel_ord_idxs) + 1; + rem_idxs[sel_idx] = 0; + } + } + db_print("Call: neg_zero_idxs\n"); + neg_zero_idxs(sel_ord_idxs, sel_idxs, var_size); + db_print("Call: arma::sort\n"); + sel_ord_idxs = arma::sort(sel_ord_idxs); + + db_print("Call: get_idxs_eq\n"); + arma::uvec tmp_sel_idxs = get_idxs_eq(sel_idxs, 1); + db_print("Call: get_sort_idxs\n"); + arma::uvec tmp_sort_idxs = get_sort_idxs(pvalues, tmp_sel_idxs); + arma::uvec tmp_sel_ord_idxs = tmp_sel_idxs(tmp_sort_idxs); + db_print("Forming return list.\n"); + res["selected"] = tmp_sel_ord_idxs; + res["stats"] = stats; + res["pvalues"] = pvalues; + res["univ"] = univs; + res["max_k"] = max_k; + res["alpha"] = thres; + res["runtime"] = get_time_spent(begin); + inits.size() ? res["n.tests"] = 0 + : res["n.tests"] = stats.n_rows * stats.n_cols + kv_length; } -Rcpp::List calc_mmp_c_bp(arma::vec& target_vars, arma::mat& ds, const int32_t limit_k, const double thres, - const std::string method) { - db_print("Initializing values\n"); - const double thres_log = std::log(thres); - arma::uvec idxs(ds.n_cols); - std::iota(idxs.begin(), idxs.end(), 0); - - uint32_t cntr = 0; - std::vector pvalues; - if (ds.n_cols == 1) { - db_print("True: ds.n_cols == 1\n"); - calc_pearson(target_vars, ds, 0, DEF_UVEC, - DEF_LIST, false, DEF_ENV, DEF_ENV); - cntr = 1; - const double tmp_pvalue = STAT_PVALUE[PVALUE_POS]; - if (tmp_pvalue > thres_log) { - idxs.fill(0); - } - pvalues.push_back(tmp_pvalue); - } - else if (ds.n_cols == 2) { - db_print("True: ds.ncols == 2\n"); - arma::uvec tmp_idxs(1); - tmp_idxs[0] = 1; - calc_pearson(target_vars, ds, 0, tmp_idxs, - DEF_LIST, false, DEF_ENV, DEF_ENV); - const double tmp_pvalue1 = STAT_PVALUE[PVALUE_POS]; - tmp_idxs[0] = 0; - calc_pearson(target_vars, ds, 1, tmp_idxs, - DEF_LIST, false, DEF_ENV, DEF_ENV); - const double tmp_pvalue2 = STAT_PVALUE[PVALUE_POS]; - cntr = 2; - if (tmp_pvalue1 > thres_log) { - db_print("True: tmp_pvalue1 > thres_log\n"); - idxs[0] = 0; - } - if (tmp_pvalue2 > thres_log) { - db_print("True: tmp_pvalue2 > thres_log\n"); - idxs[1] = 0; - } - pvalues.push_back(tmp_pvalue1); - pvalues.push_back(tmp_pvalue2); - } - else { - db_print("True: !ds.n_cols || ds.n_cols > 2\n"); - std::vector combn_vals(ds.n_cols); - std::iota(combn_vals.begin(), combn_vals.end(), 1); - for (uint32_t i = 0; i < ds.n_cols; ++i) { - int32_t k = 0; - double tmp_pvalue = -5.0; - std::vector tmp_pvalues; - while (k < (int32_t) ds.n_cols - 1 && k < limit_k && tmp_pvalue < thres_log) { - arma::umat combns = find_combn>(combn_vals, ++k); - combns -= 1; - std::vector cols_rem = det_cols(combns, i); - arma::umat combns_adj = rm_cols(combns, cols_rem); - uint32_t j = 0; - while (j < combns_adj.n_cols && tmp_pvalue < thres_log) { - arma::uvec tmp_idxs = combns_adj.col(j); - calc_pearson(target_vars, ds, i, tmp_idxs, - DEF_LIST, false, DEF_ENV, DEF_ENV); - tmp_pvalue = STAT_PVALUE[PVALUE_POS]; - tmp_pvalues.push_back(tmp_pvalue); - ++j; ++cntr; - } - } - if (tmp_pvalue > thres_log) { - idxs[i] = 0; - } - pvalues.push_back(*std::max_element(tmp_pvalues.begin(), tmp_pvalues.end())); - } - } - Rcpp::List ret; ret["met"] = idxs; ret["counter"] = cntr; ret["pvalue"] = pvalues; - return ret; +Rcpp::List calc_mmp_c_bp(arma::vec& target_vars, arma::mat& ds, + const int32_t limit_k, const double thres, + const std::string method) { + db_print("Initializing values\n"); + const double thres_log = std::log(thres); + arma::uvec idxs(ds.n_cols); + std::iota(idxs.begin(), idxs.end(), 0); + + uint32_t cntr = 0; + std::vector pvalues; + if (ds.n_cols == 1) { + db_print("True: ds.n_cols == 1\n"); + calc_pearson(target_vars, ds, 0, DEF_UVEC, DEF_LIST, false, DEF_ENV, + DEF_ENV); + cntr = 1; + const double tmp_pvalue = STAT_PVALUE[PVALUE_POS]; + if (tmp_pvalue > thres_log) { + idxs.fill(0); + } + pvalues.push_back(tmp_pvalue); + } else if (ds.n_cols == 2) { + db_print("True: ds.ncols == 2\n"); + arma::uvec tmp_idxs(1); + tmp_idxs[0] = 1; + calc_pearson(target_vars, ds, 0, tmp_idxs, DEF_LIST, false, DEF_ENV, + DEF_ENV); + const double tmp_pvalue1 = STAT_PVALUE[PVALUE_POS]; + tmp_idxs[0] = 0; + calc_pearson(target_vars, ds, 1, tmp_idxs, DEF_LIST, false, DEF_ENV, + DEF_ENV); + const double tmp_pvalue2 = STAT_PVALUE[PVALUE_POS]; + cntr = 2; + if (tmp_pvalue1 > thres_log) { + db_print("True: tmp_pvalue1 > thres_log\n"); + idxs[0] = 0; + } + if (tmp_pvalue2 > thres_log) { + db_print("True: tmp_pvalue2 > thres_log\n"); + idxs[1] = 0; + } + pvalues.push_back(tmp_pvalue1); + pvalues.push_back(tmp_pvalue2); + } else { + db_print("True: !ds.n_cols || ds.n_cols > 2\n"); + std::vector combn_vals(ds.n_cols); + std::iota(combn_vals.begin(), combn_vals.end(), 1); + for (uint32_t i = 0; i < ds.n_cols; ++i) { + int32_t k = 0; + double tmp_pvalue = -5.0; + std::vector tmp_pvalues; + while (k < (int32_t)ds.n_cols - 1 && k < limit_k && + tmp_pvalue < thres_log) { + arma::umat combns = + find_combn>(combn_vals, + ++k); + combns -= 1; + std::vector cols_rem = det_cols(combns, i); + arma::umat combns_adj = rm_cols(combns, cols_rem); + uint32_t j = 0; + while (j < combns_adj.n_cols && tmp_pvalue < thres_log) { + arma::uvec tmp_idxs = combns_adj.col(j); + calc_pearson(target_vars, ds, i, tmp_idxs, DEF_LIST, false, + DEF_ENV, DEF_ENV); + tmp_pvalue = STAT_PVALUE[PVALUE_POS]; + tmp_pvalues.push_back(tmp_pvalue); + ++j; + ++cntr; + } + } + if (tmp_pvalue > thres_log) { + idxs[i] = 0; + } + pvalues.push_back( + *std::max_element(tmp_pvalues.begin(), tmp_pvalues.end())); + } + } + Rcpp::List ret; + ret["met"] = idxs; + ret["counter"] = cntr; + ret["pvalue"] = pvalues; + return ret; } -static std::vector upd_vals(std::vector& src_pvalues, arma::uvec& idxs, - std::vector dst_pvalues) { - for (uint32_t i = 0; i < idxs.size(); ++i) { - dst_pvalues[idxs[i]] = src_pvalues[i]; - } - return dst_pvalues; +static std::vector upd_vals(std::vector& src_pvalues, + arma::uvec& idxs, + std::vector dst_pvalues) { + for (uint32_t i = 0; i < idxs.size(); ++i) { + dst_pvalues[idxs[i]] = src_pvalues[i]; + } + return dst_pvalues; } -Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int32_t max_k, - const double thres, const std::string method, Rcpp::List& inits, - const bool hash_on, Rcpp::Environment& stats_kv, - Rcpp::Environment& pvalues_kv, const bool bws_on) { - hash_on ? kv_length = stats_kv[".length"] : kv_length = 0; - db_print("Call: adj_med_NAs\n"); - adj_med_NAs(ds); - db_print("Call: check_args\n"); +Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int32_t max_k, + const double thres, const std::string method, + Rcpp::List& inits, const bool hash_on, + Rcpp::Environment& stats_kv, + Rcpp::Environment& pvalues_kv, const bool bws_on) { + hash_on ? kv_length = stats_kv[".length"] : kv_length = 0; + db_print("Call: adj_med_NAs\n"); + adj_med_NAs(ds); + db_print("Call: check_args\n"); check_args(thres, max_k, method); - db_print("Call: upd_args\n"); + db_print("Call: upd_args\n"); upd_args(target_vars, ds, method); - const uint32_t var_size = ds.n_cols; - if (max_k > (int32_t) var_size) { - db_print("True: max_k > var_size\n"); - max_k = var_size; - } - db_print("Call: inter_mmp_c\n"); - inter_mmp_c(target_vars, ds, max_k, std::log(thres), method, inits, - hash_on, var_size, stats_kv, pvalues_kv); - if (bws_on && res["selected"] != R_NilValue) { + const uint32_t var_size = ds.n_cols; + if (max_k > (int32_t)var_size) { + db_print("True: max_k > var_size\n"); + max_k = var_size; + } + db_print("Call: inter_mmp_c\n"); + inter_mmp_c(target_vars, ds, max_k, std::log(thres), method, inits, hash_on, + var_size, stats_kv, pvalues_kv); + if (bws_on && res["selected"] != R_NilValue) { db_print("True: bws_on && res[\"selected\"] != R_NilValue\n"); arma::uvec res_im_idxs = res["selected"]; - arma::uvec res_im_idxs_order = res["selected"]; - arma::mat ds_of_idxs = ds.cols(res_im_idxs); - db_print("Call: calc_mmp_c_bp\n"); - Rcpp::List res_mb = calc_mmp_c_bp(target_vars, ds_of_idxs, max_k, thres, method); - arma::uvec res_mb_idxs = res_mb["met"]; - arma::uvec sel_idxs = res_im_idxs_order(res_mb_idxs); - res["selected"] = sel_idxs; - std::vector res_im_pvalues = res["pvalues"]; - std::vector res_mb_pvalues = res_mb["pvalue"]; - res["pvalues"] = upd_vals(res_mb_pvalues, res_im_idxs, res_im_pvalues); - inits.size() ? res["n.tests"] = 0 : res["n.tests"] = (int32_t) res["n.tests"] + (int32_t) res_mb["counter"]; - } - if (hash_on) { - stats_kv[".length"] = kv_length; - pvalues_kv[".length"] = kv_length; - } - return res; + arma::uvec res_im_idxs_order = res["selected"]; + arma::mat ds_of_idxs = ds.cols(res_im_idxs); + db_print("Call: calc_mmp_c_bp\n"); + Rcpp::List res_mb = + calc_mmp_c_bp(target_vars, ds_of_idxs, max_k, thres, method); + arma::uvec res_mb_idxs = res_mb["met"]; + arma::uvec sel_idxs = res_im_idxs_order(res_mb_idxs); + res["selected"] = sel_idxs; + std::vector res_im_pvalues = res["pvalues"]; + std::vector res_mb_pvalues = res_mb["pvalue"]; + res["pvalues"] = upd_vals(res_mb_pvalues, res_im_idxs, res_im_pvalues); + inits.size() ? res["n.tests"] = 0 + : res["n.tests"] = + (int32_t)res["n.tests"] + (int32_t)res_mb["counter"]; + } + if (hash_on) { + stats_kv[".length"] = kv_length; + pvalues_kv[".length"] = kv_length; + } + return res; } diff --git a/src/mmp_c.h b/src/mmp_c.h index 971a21c..14b2e12 100644 --- a/src/mmp_c.h +++ b/src/mmp_c.h @@ -1,22 +1,23 @@ #ifndef _mmp_c_h_ #define _mmp_c_h_ -#include -#include -#include +#include #include -#include -#include +#include #include -#include +#include +#include +#include +#include #include "cts.h" // [[Rcpp::plugins("cpp11")]] // [[Rcpp::depends("RcppArmadillo")]] -Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int32_t max_k, - const double thres, const std::string method, Rcpp::List& inits, - const bool hash_on, Rcpp::Environment& stats_kv, - Rcpp::Environment& pvalues_kv, const bool bws_on); +Rcpp::List calc_mmp_c(arma::vec& target_vars, arma::mat& ds, int32_t max_k, + const double thres, const std::string method, + Rcpp::List& inits, const bool hash_on, + Rcpp::Environment& stats_kv, + Rcpp::Environment& pvalues_kv, const bool bws_on); #endif diff --git a/src/sf.cpp b/src/sf.cpp index 5eafae4..204db28 100644 --- a/src/sf.cpp +++ b/src/sf.cpp @@ -2,48 +2,52 @@ // Contact: kmdimitriadis@gmail.com #include -#include "mmp_c.h" #include "cts.h" +#include "mmp_c.h" std::vector inter(arma::vec vals1, arma::vec vals2) { - return inter_helper(vals1, vals2); + return inter_helper(vals1, vals2); } -RcppExport SEXP Rfast2_inter(SEXP xSEXP,SEXP ySEXP) { -BEGIN_RCPP +RcppExport SEXP Rfast2_inter(SEXP xSEXP, SEXP ySEXP) { + BEGIN_RCPP RObject __result; RNGScope __rngScope; - traits::input_parameter< arma::vec >::type x(xSEXP); - traits::input_parameter< arma::vec >::type y(ySEXP); - __result = inter(x,y); + traits::input_parameter::type x(xSEXP); + traits::input_parameter::type y(ySEXP); + __result = inter(x, y); return __result; -END_RCPP + END_RCPP } // [[Rcpp::export]] -Rcpp::List mmp_c(arma::vec target_vars, arma::mat ds, int32_t max_k, - const double thres, const std::string method, Rcpp::List inits, - const bool hash_on, Rcpp::Environment stats_kv, - Rcpp::Environment pvalues_kv, const bool bws_on) { - return calc_mmp_c(target_vars, ds, max_k, thres, method, - inits, hash_on, stats_kv, pvalues_kv, bws_on); +Rcpp::List mmp_c(arma::vec target_vars, arma::mat ds, int32_t max_k, + const double thres, const std::string method, Rcpp::List inits, + const bool hash_on, Rcpp::Environment stats_kv, + Rcpp::Environment pvalues_kv, const bool bws_on) { + return calc_mmp_c(target_vars, ds, max_k, thres, method, inits, hash_on, + stats_kv, pvalues_kv, bws_on); } -RcppExport SEXP Rfast2_mmp_c(SEXP target_varsSEXP,SEXP dsSEXP,SEXP max_kSEXP,SEXP thresSEXP,SEXP methodSEXP,SEXP initsSEXP,SEXP hash_onSEXP,SEXP stats_kvSEXP,SEXP pvalues_kvSEXP,SEXP bws_onSEXP){ -BEGIN_RCPP +RcppExport SEXP Rfast2_mmp_c(SEXP target_varsSEXP, SEXP dsSEXP, SEXP max_kSEXP, + SEXP thresSEXP, SEXP methodSEXP, SEXP initsSEXP, + SEXP hash_onSEXP, SEXP stats_kvSEXP, + SEXP pvalues_kvSEXP, SEXP bws_onSEXP) { + BEGIN_RCPP RObject __result; RNGScope __rngScope; - traits::input_parameter< arma::vec >::type target_vars(target_varsSEXP); - traits::input_parameter< arma::mat >::type ds(dsSEXP); - traits::input_parameter< int32_t >::type max_k(max_kSEXP); - traits::input_parameter< const double >::type thres(thresSEXP); - traits::input_parameter< const std::string >::type method(methodSEXP); - traits::input_parameter< List >::type inits(initsSEXP); - traits::input_parameter< const bool >::type hash_on(hash_onSEXP); - traits::input_parameter< Environment >::type stats_kv(stats_kvSEXP); - traits::input_parameter< Environment >::type pvalues_kv(pvalues_kvSEXP); - traits::input_parameter< const bool >::type bws_on(bws_onSEXP); - __result = mmp_c(target_vars,ds,max_k,thres,method,inits,hash_on,stats_kv,pvalues_kv,bws_on); + traits::input_parameter::type target_vars(target_varsSEXP); + traits::input_parameter::type ds(dsSEXP); + traits::input_parameter::type max_k(max_kSEXP); + traits::input_parameter::type thres(thresSEXP); + traits::input_parameter::type method(methodSEXP); + traits::input_parameter::type inits(initsSEXP); + traits::input_parameter::type hash_on(hash_onSEXP); + traits::input_parameter::type stats_kv(stats_kvSEXP); + traits::input_parameter::type pvalues_kv(pvalues_kvSEXP); + traits::input_parameter::type bws_on(bws_onSEXP); + __result = mmp_c(target_vars, ds, max_k, thres, method, inits, hash_on, + stats_kv, pvalues_kv, bws_on); return __result; -END_RCPP + END_RCPP }