R interface for raw data and sufficient statistics

Author: vandenman

R/RcppExports.R

@@ -33,8 +33,8 @@ sample_bcomrf_gibbs <- function(no_states, no_variables, no_categories, interact
     .Call(`_bgms_sample_bcomrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, variable_type, reference_category, iter)
 }
 
-sample_ggm <- function(X, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type) {
-    .Call(`_bgms_sample_ggm`, X, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type)
+sample_ggm <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type) {
+    .Call(`_bgms_sample_ggm`, inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type)
 }
 
 compute_Vn_mfm_sbm <- function(no_variables, dirichlet_alpha, t_max, lambda) {

src/RcppExports.cpp

@@ -182,12 +182,12 @@ BEGIN_RCPP
 END_RCPP
 }
 // sample_ggm
-Rcpp::List sample_ggm(const arma::mat& X, const arma::mat& prior_inclusion_prob, const arma::imat& initial_edge_indicators, const int no_iter, const int no_warmup, const int no_chains, const bool edge_selection, const int seed, const int no_threads, const int progress_type);
-RcppExport SEXP _bgms_sample_ggm(SEXP XSEXP, SEXP prior_inclusion_probSEXP, SEXP initial_edge_indicatorsSEXP, SEXP no_iterSEXP, SEXP no_warmupSEXP, SEXP no_chainsSEXP, SEXP edge_selectionSEXP, SEXP seedSEXP, SEXP no_threadsSEXP, SEXP progress_typeSEXP) {
+Rcpp::List sample_ggm(const Rcpp::List& inputFromR, const arma::mat& prior_inclusion_prob, const arma::imat& initial_edge_indicators, const int no_iter, const int no_warmup, const int no_chains, const bool edge_selection, const int seed, const int no_threads, const int progress_type);
+RcppExport SEXP _bgms_sample_ggm(SEXP inputFromRSEXP, SEXP prior_inclusion_probSEXP, SEXP initial_edge_indicatorsSEXP, SEXP no_iterSEXP, SEXP no_warmupSEXP, SEXP no_chainsSEXP, SEXP edge_selectionSEXP, SEXP seedSEXP, SEXP no_threadsSEXP, SEXP progress_typeSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::List& >::type inputFromR(inputFromRSEXP);
     Rcpp::traits::input_parameter< const arma::mat& >::type prior_inclusion_prob(prior_inclusion_probSEXP);
     Rcpp::traits::input_parameter< const arma::imat& >::type initial_edge_indicators(initial_edge_indicatorsSEXP);
     Rcpp::traits::input_parameter< const int >::type no_iter(no_iterSEXP);
@@ -197,7 +197,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const int >::type seed(seedSEXP);
     Rcpp::traits::input_parameter< const int >::type no_threads(no_threadsSEXP);
     Rcpp::traits::input_parameter< const int >::type progress_type(progress_typeSEXP);
-    rcpp_result_gen = Rcpp::wrap(sample_ggm(X, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type));
+    rcpp_result_gen = Rcpp::wrap(sample_ggm(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type));
     return rcpp_result_gen;
 END_RCPP
 }

src/ggm_model.cpp

@@ -610,3 +610,35 @@ void GGMModel::do_one_mh_step() {
     // could also be called in the main MCMC loop
     proposal_.increment_iteration();
 }
+
+
+GGMModel createGGMFromR(
+    const Rcpp::List& inputFromR,
+    const arma::mat& prior_inclusion_prob,
+    const arma::imat& initial_edge_indicators,
+    const bool edge_selection
+) {
+
+    if (inputFromR.containsElementNamed("n") && inputFromR.containsElementNamed("suf_stat")) {
+        int n = Rcpp::as<int>(inputFromR["n"]);
+        arma::mat suf_stat = Rcpp::as<arma::mat>(inputFromR["suf_stat"]);
+        return GGMModel(
+            n,
+            suf_stat,
+            prior_inclusion_prob,
+            initial_edge_indicators,
+            edge_selection
+        );
+    } else if (inputFromR.containsElementNamed("X")) {
+        arma::mat X = Rcpp::as<arma::mat>(inputFromR["X"]);
+        return GGMModel(
+            X,
+            prior_inclusion_prob,
+            initial_edge_indicators,
+            edge_selection
+        );
+    } else {
+        throw std::invalid_argument("Input list must contain either 'X' or both 'n' and 'suf_stat'.");
+    }
+
+}

src/ggm_model.h

@@ -5,6 +5,7 @@
 #include "adaptiveMetropolis.h"
 #include "rng_utils.h"
 
+
 class GGMModel : public BaseModel {
 public:
 
@@ -31,6 +32,30 @@ class GGMModel : public BaseModel {
         constants_(6)
     {}
 
+    GGMModel(
+            const int n,
+            const arma::mat& suf_stat,
+            const arma::mat& prior_inclusion_prob,
+            const arma::imat& initial_edge_indicators,
+            const bool edge_selection = true
+    ) : n_(n),
+        p_(suf_stat.n_cols),
+        dim_((p_ * (p_ + 1)) / 2),
+        suf_stat_(suf_stat),
+        prior_inclusion_prob_(prior_inclusion_prob),
+        edge_selection_(edge_selection),
+        proposal_(AdaptiveProposal(dim_, 500)),
+        omega_(arma::eye<arma::mat>(p_, p_)),
+        phi_(arma::eye<arma::mat>(p_, p_)),
+        inv_phi_(arma::eye<arma::mat>(p_, p_)),
+        inv_omega_(arma::eye<arma::mat>(p_, p_)),
+        edge_indicators_(initial_edge_indicators),
+        vectorized_parameters_(dim_),
+        vectorized_indicator_parameters_(edge_selection_ ? dim_ : 0),
+        omega_prop_(arma::mat(p_, p_, arma::fill::none)),
+        constants_(6)
+    {}
+
     GGMModel(const GGMModel& other)
         : BaseModel(other),
           dim_(other.dim_),
@@ -161,3 +186,11 @@ class GGMModel : public BaseModel {
     // double edge_log_ratio(const arma::mat& omega, size_t i, size_t j, double proposal);
     // double diag_log_ratio(const arma::mat& omega, size_t i, double proposal);
 };
+
+
+GGMModel createGGMFromR(
+    const Rcpp::List& inputFromR,
+    const arma::mat& prior_inclusion_prob,
+    const arma::imat& initial_edge_indicators,
+    const bool edge_selection = true
+);

src/sample_ggm.cpp

@@ -159,7 +159,7 @@ Rcpp::List convert_sampler_output_to_ggm_result(const std::vector<ChainResultNew
 
 // [[Rcpp::export]]
 Rcpp::List sample_ggm(
-    const arma::mat& X,
+    const Rcpp::List& inputFromR,
     const arma::mat& prior_inclusion_prob,
     const arma::imat& initial_edge_indicators,
     const int no_iter,
@@ -173,7 +173,8 @@ Rcpp::List sample_ggm(
 
     // should be done dynamically
     // also adaptation method should be specified differently
-    GGMModel model(X, prior_inclusion_prob, initial_edge_indicators, edge_selection);
+    // GGMModel model(X, prior_inclusion_prob, initial_edge_indicators, edge_selection);
+    GGMModel model = createGGMFromR(inputFromR, prior_inclusion_prob, initial_edge_indicators, edge_selection);
 
     ProgressManager pm(no_chains, no_iter, no_warmup, 50, progress_type);