Fix for stack imbalance bgmCompare parallel

Author: MaartenMarsman

NAMESPACE

@@ -16,6 +16,7 @@ S3method(print,summary.bgms)
 S3method(summary,bgms)
 export(bgm)
 export(bgmCompare)
+export(bgmCompare2)
 export(extract_arguments)
 export(extract_category_thresholds.bgms)
 export(extract_edge_indicators)

R/bgmCompare2.R

@@ -196,7 +196,7 @@ bgmCompare2 = function(
       stop("Argument 'seed' must be a non-negative integer or vector of non-negative integers.")
     }
     # Force to integer type
-    seed <- as.integer(seed)
+    seed = as.integer(seed)
     dqrng::dqset.seed(seed)
   }
 
@@ -225,44 +225,13 @@ bgmCompare2 = function(
     nuts_max_depth = nuts_max_depth,
     learn_mass_matrix = learn_mass_matrix,
     projection = projection,
-    group_membership = sorted_group - 1,  ######################################
+    group_membership = sorted_group - 1,
     group_indices = group_indices,
     interaction_index_matrix = Index,
     inclusion_probability = model$inclusion_probability_difference,
     num_chains = chains, nThreads = cores
   )
 
-  # out = run_gibbs_sampler_for_bgmCompare(
-  #   chain_id = 1,
-  #   observations = observations,
-  #   num_groups = num_groups,
-  #   num_obs_categories = num_obs_categories,
-  #   sufficient_blume_capel = sufficient_blume_capel,
-  #   sufficient_pairwise = sufficient_pairwise,
-  #   num_categories = num_categories[, 1],
-  #   main_alpha = threshold_alpha,
-  #   main_beta = threshold_beta,
-  #   pairwise_scale = interaction_scale,
-  #   difference_scale = difference_scale,
-  #   difference_selection_alpha = beta_bernoulli_alpha,
-  #   difference_selection_beta = beta_bernoulli_beta,
-  #   difference_prior = model$difference_prior, iter = iter, burnin = burnin,
-  #   na_impute = na_impute, missing_data_indices = missing_index,
-  #   is_ordinal_variable = ordinal_variable,
-  #   baseline_category = reference_category,
-  #   difference_selection = difference_selection,
-  #   main_effect_indices = main_effect_indices,
-  #   pairwise_effect_indices = pairwise_effect_indices,
-  #   target_accept = target_accept,
-  #   nuts_max_depth = nuts_max_depth,
-  #   learn_mass_matrix = learn_mass_matrix,
-  #   projection = projection,
-  #   group_membership = sorted_group - 1,  ######################################
-  #   group_indices = group_indices,
-  #   interaction_index_matrix = Index,
-  #   inclusion_probability = model$inclusion_probability_difference)
-
-
   # Main output handler in the wrapper function
   # output = prepare_output_bgmCompare2(
   #   out = out, ...

R/data_utils.R

@@ -490,7 +490,7 @@ compute_sufficient_blume_capel = function(x, reference_category, ordinal_variabl
 
 # Helper function for computing sufficient statistics for pairwise interactions
 compute_sufficient_pairwise <- function(x, group) {
-  result <- vector("list", length(group))
+  result <- vector("list", length(unique(group)))
 
   for (g in unique(group)) {
     obs <- x[group == g, , drop = FALSE]

R/function_input_utils.R

@@ -737,8 +737,10 @@ check_compare2_model = function(
         if(difference_inclusion_probability >= 1)
           stop("The inclusion probability for differences cannot equal or exceed the value one.")
 
-       inclusion_probability_difference = inclusion_probability_difference +
-         difference_probability
+       inclusion_probability_difference = matrix(difference_probability,
+                                                 nrow = ncol(x),
+                                                 ncol = ncol(x))
+
       } else {
         if(!inherits(difference_probability, what = "matrix") &&
            !inherits(difference_probability, what = "data.frame"))
@@ -765,8 +767,9 @@ check_compare2_model = function(
           stop("One or more inclusion probabilities for differences are one or larger.")
       }
     } else {
-      inclusion_probability_difference = inclusion_probability_difference + 0.5
-
+      inclusion_probability_difference = matrix(0.5,
+                                                nrow = ncol(x),
+                                                ncol = ncol(x))
       if(beta_bernoulli_alpha <= 0 || beta_bernoulli_beta <= 0)
         stop("The scale parameters of the beta distribution for the differences need to be positive.")
       if(!is.finite(beta_bernoulli_alpha) || !is.finite(beta_bernoulli_beta))

man/bgm.Rd

@@ -1,4 +1,4 @@
-CXX_STD = CXX14
+CXX_STD = CXX20
 
 ## Pull in the include-paths for RcppParallel
 PKG_CPPFLAGS = \
@@ -11,4 +11,4 @@ PKG_CPPFLAGS = \
 PKG_LIBS = \
   $(shell "${R_HOME}/bin/Rscript" -e "cat(RcppParallel::LdFlags())") \
   -Wl,-rpath,$(shell "${R_HOME}/bin/Rscript" -e "cat(system.file('lib',package='RcppParallel'))") \
-  $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+  $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

man/summarySBM.Rd

@@ -12,43 +12,43 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
 // run_bgmCompare_parallel
-Rcpp::List run_bgmCompare_parallel(const arma::imat& observations, const int num_groups, const std::vector<arma::imat>& num_obs_categories, const std::vector<arma::imat>& sufficient_blume_capel, const std::vector<arma::mat>& sufficient_pairwise, const arma::ivec& num_categories, const double main_alpha, const double main_beta, const double pairwise_scale, const double difference_scale, const double difference_selection_alpha, const double difference_selection_beta, const std::string& difference_prior, const int iter, const int burnin, const bool na_impute, const arma::imat& missing_data_indices, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, const bool difference_selection, const arma::imat& main_effect_indices, const arma::imat& pairwise_effect_indices, const double target_accept, const int nuts_max_depth, const bool learn_mass_matrix, const arma::mat& projection, const arma::ivec& group_membership, const arma::imat& group_indices, const arma::imat& interaction_index_matrix, const arma::mat& inclusion_probability, const int num_chains, const int nThreads);
+Rcpp::List run_bgmCompare_parallel(const arma::imat& observations, int num_groups, const std::vector<arma::imat>& num_obs_categories, const std::vector<arma::imat>& sufficient_blume_capel, const std::vector<arma::mat>& sufficient_pairwise, const arma::ivec& num_categories, double main_alpha, double main_beta, double pairwise_scale, double difference_scale, double difference_selection_alpha, double difference_selection_beta, const std::string& difference_prior, int iter, int burnin, bool na_impute, const arma::imat& missing_data_indices, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, bool difference_selection, const arma::imat& main_effect_indices, const arma::imat& pairwise_effect_indices, double target_accept, int nuts_max_depth, bool learn_mass_matrix, const arma::mat& projection, const arma::ivec& group_membership, const arma::imat& group_indices, const arma::imat& interaction_index_matrix, const arma::mat& inclusion_probability, int num_chains, int nThreads);
 RcppExport SEXP _bgms_run_bgmCompare_parallel(SEXP observationsSEXP, SEXP num_groupsSEXP, SEXP num_obs_categoriesSEXP, SEXP sufficient_blume_capelSEXP, SEXP sufficient_pairwiseSEXP, SEXP num_categoriesSEXP, SEXP main_alphaSEXP, SEXP main_betaSEXP, SEXP pairwise_scaleSEXP, SEXP difference_scaleSEXP, SEXP difference_selection_alphaSEXP, SEXP difference_selection_betaSEXP, SEXP difference_priorSEXP, SEXP iterSEXP, SEXP burninSEXP, SEXP na_imputeSEXP, SEXP missing_data_indicesSEXP, SEXP is_ordinal_variableSEXP, SEXP baseline_categorySEXP, SEXP difference_selectionSEXP, SEXP main_effect_indicesSEXP, SEXP pairwise_effect_indicesSEXP, SEXP target_acceptSEXP, SEXP nuts_max_depthSEXP, SEXP learn_mass_matrixSEXP, SEXP projectionSEXP, SEXP group_membershipSEXP, SEXP group_indicesSEXP, SEXP interaction_index_matrixSEXP, SEXP inclusion_probabilitySEXP, SEXP num_chainsSEXP, SEXP nThreadsSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const arma::imat& >::type observations(observationsSEXP);
-    Rcpp::traits::input_parameter< const int >::type num_groups(num_groupsSEXP);
+    Rcpp::traits::input_parameter< int >::type num_groups(num_groupsSEXP);
     Rcpp::traits::input_parameter< const std::vector<arma::imat>& >::type num_obs_categories(num_obs_categoriesSEXP);
     Rcpp::traits::input_parameter< const std::vector<arma::imat>& >::type sufficient_blume_capel(sufficient_blume_capelSEXP);
     Rcpp::traits::input_parameter< const std::vector<arma::mat>& >::type sufficient_pairwise(sufficient_pairwiseSEXP);
     Rcpp::traits::input_parameter< const arma::ivec& >::type num_categories(num_categoriesSEXP);
-    Rcpp::traits::input_parameter< const double >::type main_alpha(main_alphaSEXP);
-    Rcpp::traits::input_parameter< const double >::type main_beta(main_betaSEXP);
-    Rcpp::traits::input_parameter< const double >::type pairwise_scale(pairwise_scaleSEXP);
-    Rcpp::traits::input_parameter< const double >::type difference_scale(difference_scaleSEXP);
-    Rcpp::traits::input_parameter< const double >::type difference_selection_alpha(difference_selection_alphaSEXP);
-    Rcpp::traits::input_parameter< const double >::type difference_selection_beta(difference_selection_betaSEXP);
+    Rcpp::traits::input_parameter< double >::type main_alpha(main_alphaSEXP);
+    Rcpp::traits::input_parameter< double >::type main_beta(main_betaSEXP);
+    Rcpp::traits::input_parameter< double >::type pairwise_scale(pairwise_scaleSEXP);
+    Rcpp::traits::input_parameter< double >::type difference_scale(difference_scaleSEXP);
+    Rcpp::traits::input_parameter< double >::type difference_selection_alpha(difference_selection_alphaSEXP);
+    Rcpp::traits::input_parameter< double >::type difference_selection_beta(difference_selection_betaSEXP);
     Rcpp::traits::input_parameter< const std::string& >::type difference_prior(difference_priorSEXP);
-    Rcpp::traits::input_parameter< const int >::type iter(iterSEXP);
-    Rcpp::traits::input_parameter< const int >::type burnin(burninSEXP);
-    Rcpp::traits::input_parameter< const bool >::type na_impute(na_imputeSEXP);
+    Rcpp::traits::input_parameter< int >::type iter(iterSEXP);
+    Rcpp::traits::input_parameter< int >::type burnin(burninSEXP);
+    Rcpp::traits::input_parameter< bool >::type na_impute(na_imputeSEXP);
     Rcpp::traits::input_parameter< const arma::imat& >::type missing_data_indices(missing_data_indicesSEXP);
     Rcpp::traits::input_parameter< const arma::uvec& >::type is_ordinal_variable(is_ordinal_variableSEXP);
     Rcpp::traits::input_parameter< const arma::ivec& >::type baseline_category(baseline_categorySEXP);
-    Rcpp::traits::input_parameter< const bool >::type difference_selection(difference_selectionSEXP);
+    Rcpp::traits::input_parameter< bool >::type difference_selection(difference_selectionSEXP);
     Rcpp::traits::input_parameter< const arma::imat& >::type main_effect_indices(main_effect_indicesSEXP);
     Rcpp::traits::input_parameter< const arma::imat& >::type pairwise_effect_indices(pairwise_effect_indicesSEXP);
-    Rcpp::traits::input_parameter< const double >::type target_accept(target_acceptSEXP);
-    Rcpp::traits::input_parameter< const int >::type nuts_max_depth(nuts_max_depthSEXP);
-    Rcpp::traits::input_parameter< const bool >::type learn_mass_matrix(learn_mass_matrixSEXP);
+    Rcpp::traits::input_parameter< double >::type target_accept(target_acceptSEXP);
+    Rcpp::traits::input_parameter< int >::type nuts_max_depth(nuts_max_depthSEXP);
+    Rcpp::traits::input_parameter< bool >::type learn_mass_matrix(learn_mass_matrixSEXP);
     Rcpp::traits::input_parameter< const arma::mat& >::type projection(projectionSEXP);
     Rcpp::traits::input_parameter< const arma::ivec& >::type group_membership(group_membershipSEXP);
     Rcpp::traits::input_parameter< const arma::imat& >::type group_indices(group_indicesSEXP);
     Rcpp::traits::input_parameter< const arma::imat& >::type interaction_index_matrix(interaction_index_matrixSEXP);
     Rcpp::traits::input_parameter< const arma::mat& >::type inclusion_probability(inclusion_probabilitySEXP);
-    Rcpp::traits::input_parameter< const int >::type num_chains(num_chainsSEXP);
-    Rcpp::traits::input_parameter< const int >::type nThreads(nThreadsSEXP);
+    Rcpp::traits::input_parameter< int >::type num_chains(num_chainsSEXP);
+    Rcpp::traits::input_parameter< int >::type nThreads(nThreadsSEXP);
     rcpp_result_gen = Rcpp::wrap(run_bgmCompare_parallel(observations, num_groups, num_obs_categories, sufficient_blume_capel, sufficient_pairwise, num_categories, main_alpha, main_beta, pairwise_scale, difference_scale, difference_selection_alpha, difference_selection_beta, difference_prior, iter, burnin, na_impute, missing_data_indices, is_ordinal_variable, baseline_category, difference_selection, main_effect_indices, pairwise_effect_indices, target_accept, nuts_max_depth, learn_mass_matrix, projection, group_membership, group_indices, interaction_index_matrix, inclusion_probability, num_chains, nThreads));
     return rcpp_result_gen;
 END_RCPP

src/Makevars

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <RcppArmadillo.h>
+#include "rng_utils.h"
 
 
 
@@ -65,4 +66,40 @@ arma::vec inv_mass_active(
     const arma::imat& main_effect_indices,
     const arma::imat& pairwise_effect_indices,
     const bool& selection
-);
+);
+
+inline void initialise_graph(
+    arma::imat& indicator,
+    arma::mat& main,
+    arma::mat& pairwise,
+    const arma::imat& main_indices,
+    const arma::imat& pairwise_indices,
+    const arma::mat& incl_prob,
+    dqrng::xoshiro256plus& rng
+) {
+  int V = indicator.n_rows;
+  int G = main.n_cols;
+  for (int i = 0; i < V-1; ++i) {
+    for (int j = i+1; j < V; ++j) {
+      double p = incl_prob(i,j);
+      int draw = (runif(rng) < p) ? 1 : 0;
+      indicator(i,j) = indicator(j,i) = draw;
+      if (!draw) {
+        int row = pairwise_indices(i, j);
+        pairwise.row(row).cols(1, G-1) = 0.0;
+      }
+    }
+  }
+  for(int i = 0; i < V; i++) {
+    double p = incl_prob(i,i);
+    int draw = (runif(rng) < p) ? 1 : 0;
+    indicator(i,i) = draw;
+    if(!draw) {
+      int start = main_indices(i,0);
+      int end = main_indices(i,1);
+      for(int row = start; row < end; row++) {
+        main.row(row).cols(1, G-1) = 0.0;
+      }
+    }
+  }
+};