From 7252076cc8e83c93d78eac8446c3af09595bbc86 Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Tue, 11 Nov 2025 21:52:01 +0100 Subject: [PATCH 1/5] ggm compiles but runtime has a weird error --- R/RcppExports.R | 8 ++ src/RcppExports.cpp | 36 +++++++ src/adaptiveMetropolis.h | 55 ++++++++++ src/base_model.cpp | 0 src/base_model.h | 56 ++++++++++ src/chainResultNew.h | 32 ++++++ src/cholupdate.cpp | 129 +++++++++++++++++++++++ src/cholupdate.h | 6 ++ src/ggm_model.cpp | 215 +++++++++++++++++++++++++++++++++++++++ src/ggm_model.h | 134 ++++++++++++++++++++++++ src/sample_ggm.cpp | 185 +++++++++++++++++++++++++++++++++ test_ggm.R | 72 +++++++++++++ 12 files changed, 928 insertions(+) create mode 100644 src/adaptiveMetropolis.h create mode 100644 src/base_model.cpp create mode 100644 src/base_model.h create mode 100644 src/chainResultNew.h create mode 100644 src/cholupdate.cpp create mode 100644 src/cholupdate.h create mode 100644 src/ggm_model.cpp create mode 100644 src/ggm_model.h create mode 100644 src/sample_ggm.cpp create mode 100644 test_ggm.R diff --git a/R/RcppExports.R b/R/RcppExports.R index 6e04e6f1..6c7fab37 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,6 +9,10 @@ run_bgm_parallel <- function(observations, num_categories, pairwise_scale, edge_ .Call(`_bgms_run_bgm_parallel`, observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type) } +chol_update_arma <- function(R, u, downdate = FALSE, eps = 1e-12) { + .Call(`_bgms_chol_update_arma`, R, u, downdate, eps) +} + get_explog_switch <- function() { .Call(`_bgms_get_explog_switch`) } @@ -29,6 +33,10 @@ 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) +} + compute_Vn_mfm_sbm <- function(no_variables, dirichlet_alpha, t_max, lambda) { .Call(`_bgms_compute_Vn_mfm_sbm`, no_variables, dirichlet_alpha, t_max, lambda) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f3fc9664..65e29044 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -101,6 +101,20 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// chol_update_arma +arma::mat chol_update_arma(arma::mat& R, arma::vec& u, bool downdate, double eps); +RcppExport SEXP _bgms_chol_update_arma(SEXP RSEXP, SEXP uSEXP, SEXP downdateSEXP, SEXP epsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat& >::type R(RSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type u(uSEXP); + Rcpp::traits::input_parameter< bool >::type downdate(downdateSEXP); + Rcpp::traits::input_parameter< double >::type eps(epsSEXP); + rcpp_result_gen = Rcpp::wrap(chol_update_arma(R, u, downdate, eps)); + return rcpp_result_gen; +END_RCPP +} // get_explog_switch Rcpp::String get_explog_switch(); RcppExport SEXP _bgms_get_explog_switch() { @@ -167,6 +181,26 @@ BEGIN_RCPP return rcpp_result_gen; 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, 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) { +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 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); + Rcpp::traits::input_parameter< const int >::type no_warmup(no_warmupSEXP); + Rcpp::traits::input_parameter< const int >::type no_chains(no_chainsSEXP); + Rcpp::traits::input_parameter< const bool >::type edge_selection(edge_selectionSEXP); + Rcpp::traits::input_parameter< const int >::type seed(seedSEXP); + Rcpp::traits::input_parameter< const int >::type no_threads(no_threadsSEXP); + Rcpp::traits::input_parameter< 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)); + return rcpp_result_gen; +END_RCPP +} // compute_Vn_mfm_sbm arma::vec compute_Vn_mfm_sbm(arma::uword no_variables, double dirichlet_alpha, arma::uword t_max, double lambda); RcppExport SEXP _bgms_compute_Vn_mfm_sbm(SEXP no_variablesSEXP, SEXP dirichlet_alphaSEXP, SEXP t_maxSEXP, SEXP lambdaSEXP) { @@ -185,11 +219,13 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_bgms_run_bgmCompare_parallel", (DL_FUNC) &_bgms_run_bgmCompare_parallel, 36}, {"_bgms_run_bgm_parallel", (DL_FUNC) &_bgms_run_bgm_parallel, 34}, + {"_bgms_chol_update_arma", (DL_FUNC) &_bgms_chol_update_arma, 4}, {"_bgms_get_explog_switch", (DL_FUNC) &_bgms_get_explog_switch, 0}, {"_bgms_rcpp_ieee754_exp", (DL_FUNC) &_bgms_rcpp_ieee754_exp, 1}, {"_bgms_rcpp_ieee754_log", (DL_FUNC) &_bgms_rcpp_ieee754_log, 1}, {"_bgms_sample_omrf_gibbs", (DL_FUNC) &_bgms_sample_omrf_gibbs, 6}, {"_bgms_sample_bcomrf_gibbs", (DL_FUNC) &_bgms_sample_bcomrf_gibbs, 8}, + {"_bgms_sample_ggm", (DL_FUNC) &_bgms_sample_ggm, 10}, {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4}, {NULL, NULL, 0} }; diff --git a/src/adaptiveMetropolis.h b/src/adaptiveMetropolis.h new file mode 100644 index 00000000..f8df5c63 --- /dev/null +++ b/src/adaptiveMetropolis.h @@ -0,0 +1,55 @@ +#pragma once + +#include +#include + +class AdaptiveProposal { + +public: + + AdaptiveProposal(size_t num_params, size_t adaption_window = 50, double target_accept = 0.44) { + proposal_sds_ = arma::vec(num_params, arma::fill::ones) * 0.25; // Initial SD + // acceptance_counts_ = arma::ivec(num_params, arma::fill::zeros); + // total_proposals_ = arma::ivec(num_params, arma::fill::zeros); + adaptation_window_ = adaption_window; + target_accept_ = target_accept; + } + + double get_proposal_sd(size_t param_index) const { + validate_index(param_index); + return proposal_sds_[param_index]; + } + + void update_proposal_sd(size_t param_index, double alpha) { + + if (!adapting_) { + return; + } + + double current_sd = get_proposal_sd(param_index); + double updated_sd = current_sd + std::pow(1.0 / iterations_, 0.6) * (alpha - target_accept_); + // proposal_sds_[param_index] = std::min(20.0, std::max(1.0 / std::sqrt(n), updated_sd)); + proposal_sds_(param_index) = std::min(20.0, updated_sd); + } + + void increment_iteration() { + iterations_++; + if (iterations_ >= adaptation_window_) { + adapting_ = false; + } + } + +private: + arma::vec proposal_sds_; + int iterations_ = 0, + adaptation_window_; + double target_accept_ = 0.44; + bool adapting_ = true; + + void validate_index(size_t index) const { + if (index >= proposal_sds_.n_elem) { + throw std::out_of_range("Parameter index out of range"); + } + } + +}; diff --git a/src/base_model.cpp b/src/base_model.cpp new file mode 100644 index 00000000..e69de29b diff --git a/src/base_model.h b/src/base_model.h new file mode 100644 index 00000000..20cd4930 --- /dev/null +++ b/src/base_model.h @@ -0,0 +1,56 @@ +#pragma once + +#include +#include +#include + +class BaseModel { +public: + virtual ~BaseModel() = default; + + // Capability queries + virtual bool has_gradient() const { return false; } + virtual bool has_adaptive_mh() const { return false; } + + // Core methods (to be overridden by derived classes) + virtual double logp(const std::vector& parameters) = 0; + + virtual arma::vec gradient(const std::vector& parameters) { + if (!has_gradient()) { + throw std::runtime_error("Gradient not implemented for this model"); + } + throw std::runtime_error("Gradient method must be implemented in derived class"); + } + + virtual std::pair logp_and_gradient( + const std::vector& parameters) { + if (!has_gradient()) { + throw std::runtime_error("Gradient not implemented for this model"); + } + return {logp(parameters), gradient(parameters)}; + } + + // For Metropolis-Hastings (model handles parameter groups internally) + virtual void do_one_mh_step() { + throw std::runtime_error("do_one_mh_step method must be implemented in derived class"); + } + + virtual arma::vec get_vectorized_parameters() { + throw std::runtime_error("get_vectorized_parameters method must be implemented in derived class"); + } + + // Return dimensionality of the parameter space + virtual size_t parameter_dimension() const = 0; + + virtual void set_seed(int seed) { + throw std::runtime_error("set_seed method must be implemented in derived class"); + } + + virtual std::unique_ptr clone() const { + throw std::runtime_error("clone method must be implemented in derived class"); + } + + +protected: + BaseModel() = default; +}; diff --git a/src/chainResultNew.h b/src/chainResultNew.h new file mode 100644 index 00000000..5a6dd855 --- /dev/null +++ b/src/chainResultNew.h @@ -0,0 +1,32 @@ +#pragma once + +#include +#include + +class ChainResultNew { + +public: + ChainResultNew() {} + + bool error; + std::string error_msg; + int chain_id; + bool userInterrupt; + + arma::mat samples; + + void reserve(size_t param_dim, size_t n_iter) { + samples.set_size(param_dim, n_iter); + } + void store_sample(size_t iter, const arma::vec& sample) { + samples.col(iter) = sample; + } + + // arma::imat indicator_samples; + + // other samples + // arma::ivec treedepth_samples; + // arma::ivec divergent_samples; + // arma::vec energy_samples; + // arma::imat allocation_samples; +}; diff --git a/src/cholupdate.cpp b/src/cholupdate.cpp new file mode 100644 index 00000000..5ab8f6eb --- /dev/null +++ b/src/cholupdate.cpp @@ -0,0 +1,129 @@ +#include "cholupdate.h" + +extern "C" { + +// from mgcv: https://github.com/cran/mgcv/blob/1b6a4c8374612da27e36420b4459e93acb183f2d/src/mat.c#L1876-L1883 +static inline double hypote(double x, double y) { +/* stable computation of sqrt(x^2 + y^2) */ + double t; + x = fabs(x);y=fabs(y); + if (y>x) { t = x;x = y; y = t;} + if (x==0) return(y); else t = y/x; + return(x*sqrt(1+t*t)); +} /* hypote */ + +// from mgcv: https://github.com/cran/mgcv/blob/1b6a4c8374612da27e36420b4459e93acb183f2d/src/mat.c#L1956 +void chol_up(double *R,double *u, int *n,int *up,double *eps) { +/* Rank 1 update of a cholesky factor. Works as follows: + + [up=1] R'R + uu' = [u,R'][u,R']' = [u,R']Q'Q[u,R']', and then uses Givens rotations to + construct Q such that Q[u,R']' = [0,R1']'. Hence R1'R1 = R'R + uu'. The construction + operates from first column to last. + + [up=0] uses an almost identical sequence, but employs hyperbolic rotations + in place of Givens. See Golub and van Loan (2013, 4e 6.5.4) + + Givens rotations are of form [c,-s] where c = cos(theta), s = sin(theta). + [s,c] + + Assumes R upper triangular, and that it is OK to use first two columns + below diagonal as temporary strorage for Givens rotations (the storage is + needed to ensure algorithm is column oriented). + + For downdate returns a negative value in R[1] (R[1,0]) if not +ve definite. +*/ + double c0,s0,*c,*s,z,*x,z0,*c1; + int j,j1,n1; + n1 = *n - 1; + if (*up) for (j1=-1,j=0;j<*n;j++,u++,j1++) { /* loop over columns of R */ + z = *u; /* initial element of u */ + x = R + *n * j; /* current column */ + c = R + 2;s = R + *n + 2; /* Storage for first n-2 Givens rotations */ + for (c1=c+j1;c R[j,j] */ + z0 = hypote(z,*x); /* sqrt(z^2+R[j,j]^2) */ + c0 = *x/z0; s0 = z/z0; /* need to zero z */ + /* now apply this rotation and this column is finished (so no need to update z) */ + *x = s0 * z + c0 * *x; + } else for (j1=-1,j=0;j<*n;j++,u++,j1++) { /* loop over columns of R for down-dating */ + z = *u; /* initial element of u */ + x = R + *n * j; /* current column */ + c = R + 2;s = R + *n + 2; /* Storage for first n-2 hyperbolic rotations */ + for (c1=c+j1;c R[j,j] */ + z0 = z / *x; /* sqrt(z^2+R[j,j]^2) */ + if (fabs(z0)>=1) { /* downdate not +ve def */ + //Rprintf("j = %d d = %g ",j,z0); + if (*n>1) R[1] = -2.0; + return; /* signals error */ + } + if (z0 > 1 - *eps) z0 = 1 - *eps; + c0 = 1/sqrt(1-z0*z0);s0 = c0 * z0; + /* now apply this rotation and this column is finished (so no need to update z) */ + *x = -s0 * z + c0 * *x; + } + + /* now zero c and s storage */ + c = R + 2;s = R + *n + 2; + for (x = c + *n - 2;c + +void cholesky_update( arma::mat& R, arma::vec& u, double eps = 1e-12); +void cholesky_downdate(arma::mat& R, arma::vec& u, double eps = 1e-12); diff --git a/src/ggm_model.cpp b/src/ggm_model.cpp new file mode 100644 index 00000000..eee45bf7 --- /dev/null +++ b/src/ggm_model.cpp @@ -0,0 +1,215 @@ +#include "ggm_model.h" +#include "adaptiveMetropolis.h" +#include "rng_utils.h" +#include "cholupdate.h" + +double GGMModel::compute_inv_submatrix_i(const arma::mat& A, const size_t i, const size_t ii, const size_t jj) const { + return(A(ii, jj) - A(ii, i) * A(jj, i) / A(i, i)); +} + +void GGMModel::get_constants(size_t i, size_t j) { + + // TODO: helper function? + double logdet_omega = 0.0; + for (size_t i = 0; i < p_; i++) { + logdet_omega += std::log(phi_(i, i)); + } + + double log_adj_omega_ii = logdet_omega + log(abs(inv_omega_(i, i))); + double log_adj_omega_ij = logdet_omega + log(abs(inv_omega_(i, j))); + double log_adj_omega_jj = logdet_omega + log(abs(inv_omega_(j, j))); + + double inv_omega_sub_j1j1 = compute_inv_submatrix_i(inv_omega_, i, j, j); + double log_abs_inv_omega_sub_jj = log_adj_omega_ii + log(abs(inv_omega_sub_j1j1)); + + double Phi_q1q = (-1 * (2 * std::signbit(inv_omega_(i, j)) - 1)) * std::exp( + (log_adj_omega_ij - (log_adj_omega_jj + log_abs_inv_omega_sub_jj) / 2) + ); + double Phi_q1q1 = std::exp((log_adj_omega_jj - log_abs_inv_omega_sub_jj) / 2); + + constants_[1] = Phi_q1q; + constants_[2] = Phi_q1q1; + constants_[3] = omega_(i, j) - Phi_q1q * Phi_q1q1; + constants_[4] = Phi_q1q1; + constants_[5] = omega_(j, j) - Phi_q1q * Phi_q1q; + constants_[6] = constants_[5] + constants_[3] * constants_[3] / (constants_[4] * constants_[4]); + +} + +double GGMModel::R(const double x) const { + if (x == 0) { + return constants_[6]; + } else { + return constants_[3] + std::pow((x - constants_[3]) / constants_[4], 2); + } +} + +void GGMModel::update_edge_parameter(size_t i, size_t j) { + + if (edge_indicators_(i, j) == 0) { + return; // Edge is not included; skip update + } + + get_constants(i, j); + double Phi_q1q = constants_[1]; + double Phi_q1q1 = constants_[2]; + + size_t e = i * (i + 1) / 2 + j; // parameter index in vectorized form + double proposal_sd = proposal_.get_proposal_sd(e); + + double phi_prop = rnorm(rng_, Phi_q1q, proposal_sd); + double omega_prop_q1q = constants_[3] + constants_[4] * phi_prop; + double omega_prop_qq = R(omega_prop_q1q); + + // form full proposal matrix for Omega + omega_prop_ = omega_; // TODO: needs to be a copy! + omega_prop_(i, j) = omega_prop_q1q; + omega_prop_(j, i) = omega_prop_q1q; + omega_prop_(j, j) = omega_prop_qq; + + double ln_alpha = log_density(omega_prop_) - log_density(); + ln_alpha += R::dcauchy(omega_prop_(i, j), 0.0, 2.5, true); + ln_alpha -= R::dcauchy(omega_(i, j), 0.0, 2.5, true); + + double u = runif(rng_); + if (ln_alpha > log(u)) { + // accept proposal + + double omega_ij = omega_(i, j); + double omega_jj = omega_(j, j); + + omega_(i, j) = omega_prop_q1q; + omega_(j, i) = omega_prop_q1q; + omega_(j, j) = omega_prop_qq; + + // TODO: preallocate? + // find v for low rank update + arma::vec v1 = {0, -1}; + arma::vec v2 = {omega_ij - omega_prop_(i, j), (omega_jj - omega_prop_(j, j)) / 2}; + + arma::vec vf1 = arma::zeros(p_); + arma::vec vf2 = arma::zeros(p_); + vf1[i] = v1[1]; + vf1[j] = v1[2]; + vf2[i] = v2[1]; + vf2[j] = v2[2]; + + // we now have + // aOmega_prop - (aOmega + vf1 %*% t(vf2) + vf2 %*% t(vf1)) + + arma::vec u1 = (vf1 + vf2) / sqrt(2); + arma::vec u2 = (vf1 - vf2) / sqrt(2); + + // we now have + // omega_prop_ - (aOmega + u1 %*% t(u1) - u2 %*% t(u2)) + // and also + // aOmega_prop - (aOmega + cbind(vf1, vf2) %*% matrix(c(0, 1, 1, 0), 2, 2) %*% t(cbind(vf1, vf2))) + + // update phi + cholesky_update(phi_, u1); + cholesky_downdate(phi_, u2); + + // update inverse + inv_omega_ = phi_.t() * phi_; + + } + + double alpha = std::min(1.0, std::exp(ln_alpha)); + proposal_.update_proposal_sd(e, alpha); +} + +void GGMModel::do_one_mh_step() { + + // Update off-diagonals (upper triangle) + for (size_t i = 0; i < p_ - 1; ++i) { + for (size_t j = i + 1; j < p_; ++j) { + Rcpp::Rcout << "Updating edge parameter (" << i << ", " << j << ")" << std::endl; + update_edge_parameter(i, j); + } + } + + // Update diagonals + for (size_t i = 0; i < p_; ++i) { + Rcpp::Rcout << "Updating diagonal parameter " << i << std::endl; + update_diagonal_parameter(i); + } + + // if (edge_selection_) { + // for (size_t i = 0; i < p_ - 1; ++i) { + // for (size_t j = i + 1; j < p_; ++j) { + // update_edge_indicator_parameter_pair(i, j); + // } + // } + // } + proposal_.increment_iteration(); +} + +double GGMModel::log_density_impl(const arma::mat& omega, const arma::mat& phi) const { + double logdet_omega = 0.0; + for (size_t i = 0; i < p_; i++) { + logdet_omega += std::log(phi(i, i)); + } + + // TODO: does this allocate? + double trace_prod = arma::accu(omega % suf_stat_); + + double log_likelihood = n_ * (p_ * log(2 * arma::datum::pi) / 2 + logdet_omega / 2) - trace_prod / 2; + + return log_likelihood; +} + + +void GGMModel::update_diagonal_parameter(size_t i) { + // Implementation of diagonal parameter update + // 1-3) from before + double logdet_omega = 0.0; + for (size_t i = 0; i < p_; i++) { + logdet_omega += std::log(phi_(i, i)); + } + + double logdet_omega_sub_ii = logdet_omega + std::log(inv_omega_(i, i)); + + size_t e = i * (i + 1) / 2 + i; // parameter index in vectorized form + double proposal_sd = proposal_.get_proposal_sd(e); + + double theta_curr = (logdet_omega - logdet_omega_sub_ii) / 2; + double theta_prop = rnorm(rng_, theta_curr, proposal_sd); + + //4) Replace and rebuild omega + omega_prop_ = omega_; + omega_prop_(i, i) = omega_(i, i) - std::exp(theta_curr) * std::exp(theta_curr) + std::exp(theta_prop) * std::exp(theta_prop); + + // 5) Acceptance ratio + double ln_alpha = log_density(omega_prop_) - log_density(); + ln_alpha += R::dgamma(exp(theta_prop), 1.0, 1.0, true); + ln_alpha -= R::dgamma(exp(theta_curr), 1.0, 1.0, true); + ln_alpha += theta_prop - theta_curr; // Jacobian adjustment ? + + + if (log(runif(rng_)) < ln_alpha) { + + double omega_ii = omega_(i, i); + + arma::vec u(p_, arma::fill::zeros); + double delta = omega_ii - omega_prop_(i, i); + bool s = delta > 0; + u(i) = sqrt(abs(delta)); + + omega_(i, i) = omega_prop_(i, i); + + if (!s) + cholesky_update(phi_, u); + else + cholesky_downdate(phi_, u); + + inv_omega_ = phi_.t() * phi_; + + } + + double alpha = std::min(1.0, std::exp(ln_alpha)); + proposal_.update_proposal_sd(e, alpha); +} + +void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { + // Implementation of edge indicator parameter pair update +} diff --git a/src/ggm_model.h b/src/ggm_model.h new file mode 100644 index 00000000..05081db0 --- /dev/null +++ b/src/ggm_model.h @@ -0,0 +1,134 @@ +#pragma once + +#include +#include "base_model.h" +#include "adaptiveMetropolis.h" +#include "rng_utils.h" + +class GGMModel : public BaseModel { +public: + + GGMModel( + const arma::mat& X, + const arma::mat& prior_inclusion_prob, + const arma::imat& initial_edge_indicators, + const bool edge_selection = true + ) : n_(X.n_rows), + p_(X.n_cols), + dim_((p_ * (p_ + 1)) / 2), + suf_stat_(X.t() * X), + prior_inclusion_prob_(prior_inclusion_prob), + edge_selection_(edge_selection), + proposal_(AdaptiveProposal(dim_, 500)), + omega_(arma::eye(p_, p_)), + phi_(arma::eye(p_, p_)), + inv_omega_(arma::eye(p_, p_)), + edge_indicators_(initial_edge_indicators), + vectorized_parameters_(dim_), + constants_(6) + {} + + GGMModel(const GGMModel& other) + : BaseModel(other), + dim_(other.dim_), + suf_stat_(other.suf_stat_), + n_(other.n_), + p_(other.p_), + prior_inclusion_prob_(other.prior_inclusion_prob_), + edge_selection_(other.edge_selection_), + omega_(other.omega_), + phi_(other.phi_), + inv_omega_(other.inv_omega_), + edge_indicators_(other.edge_indicators_), + vectorized_parameters_(other.vectorized_parameters_), + proposal_(other.proposal_), + rng_(other.rng_), + omega_prop_(other.omega_prop_), + constants_(other.constants_) + {} + + // // rng_ = SafeRNG(123); + + // } + + void set_adaptive_proposal(AdaptiveProposal proposal) { + proposal_ = proposal; + } + + virtual bool has_gradient() const { return false; } + virtual bool has_adaptive_mh() const override { return true; } + + double logp(const std::vector& parameters) override { + // Implement log probability computation + return 0.0; + } + + // TODO: this can be done more efficiently, no need for the Cholesky! + double log_density(const arma::mat& omega) const { return log_density_impl(omega, arma::chol(omega)); }; + double log_density() const { return log_density_impl(omega_, phi_); } + + void do_one_mh_step() override; + + size_t parameter_dimension() const override { + Rcpp::Rcout << "GGMModel::parameter_dimension() returning: " << dim_ << std::endl; + return dim_; + } + + void set_seed(int seed) override { + rng_ = SafeRNG(seed); + } + + arma::vec get_vectorized_parameters() override { + size_t e = 0; + for (size_t j = 0; j < p_; ++j) { + for (size_t i = 0; i <= j; ++i) { + vectorized_parameters_(e) = omega_(i, j); + e++; + } + } + return vectorized_parameters_; + } + + std::unique_ptr clone() const override { + return std::make_unique(*this); // uses copy constructor + } + +private: + // data + size_t dim_; + arma::mat suf_stat_; + size_t n_; + size_t p_; + arma::mat prior_inclusion_prob_; + bool edge_selection_ = true; + + // parameters + arma::mat omega_, phi_, inv_omega_; + arma::imat edge_indicators_; + arma::vec vectorized_parameters_; + + + AdaptiveProposal proposal_; + SafeRNG rng_; + + // internal helper variables + arma::mat omega_prop_; + arma::vec constants_; // Phi_q1q, Phi_q1q1, c[1], c[2], c[3], c[4] + + // Parameter group updates with optimized likelihood evaluations + void update_edge_parameter(size_t i, size_t j); + void update_diagonal_parameter(size_t i); + void update_edge_indicator_parameter_pair(size_t i, size_t j); + + // Helper methods + void get_constants(size_t i, size_t j); + double compute_inv_submatrix_i(const arma::mat& A, const size_t i, const size_t ii, const size_t jj) const; + double R(const double x) const; + + double log_density_impl(const arma::mat& omega, const arma::mat& phi) const; + + // double find_reasonable_step_size_edge(const arma::mat& omega, size_t i, size_t j); + // double find_reasonable_step_size_diag(const arma::mat& omega, size_t i); + // 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); +}; diff --git a/src/sample_ggm.cpp b/src/sample_ggm.cpp new file mode 100644 index 00000000..0d9de749 --- /dev/null +++ b/src/sample_ggm.cpp @@ -0,0 +1,185 @@ +#include +#include +#include +#include +#include + +#include "ggm_model.h" +#include "progress_manager.h" +#include "chainResultNew.h" + +void run_mcmc_sampler_single_thread( + ChainResultNew& chain_result, + BaseModel& model, + const int no_iter, + const int no_warmup, + const int chain_id, + ProgressManager& pm +) { + + size_t i = 0; + for (size_t iter = 0; iter < no_iter + no_warmup; ++iter) { + + model.do_one_mh_step(); + + if (iter >= no_warmup) { + + chain_result.store_sample(i, model.get_vectorized_parameters()); + i++; + } + + pm.update(chain_id); + if (pm.shouldExit()) { + chain_result.userInterrupt = true; + break; + } + } +} + +struct GGMChainRunner : public RcppParallel::Worker { + std::vector& results_; + std::vector>& models_; + size_t no_iter_; + size_t no_warmup_; + int seed_; + ProgressManager& pm_; + + GGMChainRunner( + std::vector& results, + std::vector>& models, + const size_t no_iter, + const size_t no_warmup, + const int seed, + ProgressManager& pm + ) : + results_(results), + models_(models), + no_iter_(no_iter), + no_warmup_(no_warmup), + seed_(seed), + pm_(pm) + {} + + void operator()(std::size_t begin, std::size_t end) { + for (std::size_t i = begin; i < end; ++i) { + + ChainResultNew& chain_result = results_[i]; + chain_result.chain_id = static_cast(i + 1); + chain_result.error = false; + BaseModel& model = *models_[i]; + model.set_seed(seed_ + i); + try { + + run_mcmc_sampler_single_thread(chain_result, model, no_iter_, no_warmup_, i, pm_); + + } catch (std::exception& e) { + chain_result.error = true; + chain_result.error_msg = e.what(); + } catch (...) { + chain_result.error = true; + chain_result.error_msg = "Unknown error"; + } + } + } +}; + +void run_mcmc_sampler_threaded( + std::vector& results, + std::vector>& models, + const int no_iter, + const int no_warmup, + const int seed, + const int no_threads, + ProgressManager& pm +) { + + GGMChainRunner runner(results, models, no_iter, no_warmup, seed, pm); + tbb::global_control control(tbb::global_control::max_allowed_parallelism, no_threads); + RcppParallel::parallelFor(0, results.size(), runner); +} + + +std::vector run_mcmc_sampler( + BaseModel& model, + const int no_iter, + const int no_warmup, + const int no_chains, + const int seed, + const int no_threads, + ProgressManager& pm +) { + + Rcpp::Rcout << "Allocating results objects..." << std::endl; + std::vector results(no_chains); + for (size_t c = 0; c < no_chains; ++c) { + results[c].reserve(model.parameter_dimension(), no_iter); + } + + if (no_threads > 1) { + + Rcpp::Rcout << "Running multi-threaded MCMC sampling..." << std::endl; + std::vector> models; + models.reserve(no_chains); + for (size_t c = 0; c < no_chains; ++c) { + models.push_back(model.clone()); // deep copy via virtual clone + } + run_mcmc_sampler_threaded(results, models, no_iter, no_warmup, seed, no_threads, pm); + } else { + model.set_seed(seed); + Rcpp::Rcout << "Running single-threaded MCMC sampling..." << std::endl; + for (size_t c = 0; c < no_chains; ++c) { + run_mcmc_sampler_single_thread(results[c], model, no_iter, no_warmup, c, pm); + } + } + return results; +} + +Rcpp::List convert_sampler_output_to_ggm_result(const std::vector& results) { + + Rcpp::List output(results.size()); + for (size_t i = 0; i < results.size(); ++i) { + + Rcpp::List chain_i; + chain_i["chain_id"] = results[i].chain_id; + if (results[i].error) { + chain_i["error"] = results[i].error_msg; + } else { + chain_i["samples"] = results[i].samples; + chain_i["userInterrupt"] = results[i].userInterrupt; + + } + output[i] = chain_i; + } + return output; +} + +// [[Rcpp::export]] +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, + int progress_type +) { + + // should be done dynamically + // also adaptation method should be specified differently + GGMModel model(X, prior_inclusion_prob, initial_edge_indicators, edge_selection); + + Rcpp::Rcout << "GGMModel::parameter_dimension() returning: " << model.parameter_dimension() << std::endl; + + ProgressManager pm(no_chains, no_iter, no_warmup, 50, progress_type); + + std::vector output = run_mcmc_sampler(model, no_iter, no_warmup, no_chains, seed, no_threads, pm); + + Rcpp::List ggm_result = convert_sampler_output_to_ggm_result(output); + + pm.finish(); + + return ggm_result; +} \ No newline at end of file diff --git a/test_ggm.R b/test_ggm.R new file mode 100644 index 00000000..60c5ef0f --- /dev/null +++ b/test_ggm.R @@ -0,0 +1,72 @@ +library(bgms) + +# Dimension and true precision +p <- 10 + +adj <- matrix(0, nrow = p, ncol = p) +adj[lower.tri(adj)] <- rbinom(p * (p - 1) / 2, size = 1, prob = 0.3) +adj <- adj + t(adj) +# qgraph::qgraph(adj) +Omega <- BDgraph::rgwish(1, adj = adj, b = p + sample(0:p, 1), D = diag(p)) +Sigma <- solve(Omega) +zapsmall(Omega) + +# Data +n <- 1e2 +x <- mvtnorm::rmvnorm(n = n, mean = rep(0, p), sigma = Sigma) + + +# ---- Run MCMC with warmup and sampling ------------------------------------ + +# debugonce(mbgms:::bgm_gaussian) +sampling_results <- bgms:::sample_ggm( + X = x, + prior_inclusion_prob = matrix(.5, p, p), + initial_edge_indicators = adj, + no_iter = 4000, + no_warmup = 4000, + no_chains = 3, + edge_selection = FALSE, + no_threads = 1, + seed = 123, + progress_type = 2 +) + +profvis::profvis({ + sampling_results <- bgm_gaussian( + x = x, + n = n, + n_iter = 400, + n_warmup = 400, + n_phases = 10 + ) +}) + +# Extract results +aveOmega <- sampling_results$aveOmega +aveGamma <- sampling_results$aveGamma +aOmega <- sampling_results$aOmega +aGamma <- sampling_results$aGamma +prob <- sampling_results$prob +proposal_sd <- sampling_results$proposal_sd + +library(patchwork) +library(ggplot2) +df <- data.frame( + true = aveOmega[lower.tri(aveOmega)], + Omega[lower.tri(Omega)], + estimated = aveOmega[lower.tri(aveOmega)], + p_inclusion = aveGamma[lower.tri(aveGamma)] +) +p1 <- ggplot(df, aes(x = true, y = estimated)) + + geom_point(size = 5, alpha = 0.8, shape = 21, fill = "grey") + + geom_abline(slope = 1, intercept = 0, color = "grey") + + labs(x = "True Values Omega", y = "Estimated Values Omega (Posterior Mean)") +p2 <- ggplot(df, aes(x = estimated, y = p_inclusion)) + + geom_point(size = 5, alpha = 0.8, shape = 21, fill = "grey") + + labs( + x = "Estimated Values Omega (Posterior Mean)", + y = "Estimated Inclusion Probabilities" + ) +(p1 + p2) + plot_layout(ncol = 1) & theme_bw(base_size = 20) + From f719bdd475f954e82a51081b254c2568c5d7ebb0 Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Wed, 12 Nov 2025 16:49:42 +0100 Subject: [PATCH 2/5] almost functional and O(P^2) per update --- src/RcppExports.cpp | 4 +- src/adaptiveMetropolis.h | 37 +++- src/base_model.h | 10 +- src/chainResultNew.h | 12 +- src/ggm_model.cpp | 460 ++++++++++++++++++++++++++++++++------- src/ggm_model.h | 40 +++- src/sample_ggm.cpp | 14 +- test_ggm.R | 50 ++++- 8 files changed, 510 insertions(+), 117 deletions(-) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 65e29044..0b7e900f 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -182,7 +182,7 @@ 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, int progress_type); +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) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -196,7 +196,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const bool >::type edge_selection(edge_selectionSEXP); Rcpp::traits::input_parameter< const int >::type seed(seedSEXP); Rcpp::traits::input_parameter< const int >::type no_threads(no_threadsSEXP); - Rcpp::traits::input_parameter< int >::type progress_type(progress_typeSEXP); + 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)); return rcpp_result_gen; END_RCPP diff --git a/src/adaptiveMetropolis.h b/src/adaptiveMetropolis.h index f8df5c63..d8d9cb6c 100644 --- a/src/adaptiveMetropolis.h +++ b/src/adaptiveMetropolis.h @@ -8,9 +8,8 @@ class AdaptiveProposal { public: AdaptiveProposal(size_t num_params, size_t adaption_window = 50, double target_accept = 0.44) { - proposal_sds_ = arma::vec(num_params, arma::fill::ones) * 0.25; // Initial SD - // acceptance_counts_ = arma::ivec(num_params, arma::fill::zeros); - // total_proposals_ = arma::ivec(num_params, arma::fill::zeros); + proposal_sds_ = arma::vec(num_params, arma::fill::ones) * 0.25; // Initial SD, need to tweak this somehow? + acceptance_counts_ = arma::ivec(num_params, arma::fill::zeros); adaptation_window_ = adaption_window; target_accept_ = target_accept; } @@ -20,16 +19,26 @@ class AdaptiveProposal { return proposal_sds_[param_index]; } - void update_proposal_sd(size_t param_index, double alpha) { + void update_proposal_sd(size_t param_index) { if (!adapting_) { return; } double current_sd = get_proposal_sd(param_index); - double updated_sd = current_sd + std::pow(1.0 / iterations_, 0.6) * (alpha - target_accept_); - // proposal_sds_[param_index] = std::min(20.0, std::max(1.0 / std::sqrt(n), updated_sd)); - proposal_sds_(param_index) = std::min(20.0, updated_sd); + double observed_acceptance_probability = acceptance_counts_[param_index] / static_cast(iterations_ + 1); + double rm_weight = std::pow(iterations_, -decay_rate_); + + // Robbins-Monro update step + double updated_sd = current_sd + (observed_acceptance_probability - target_accept_) * rm_weight; + updated_sd = std::clamp(updated_sd, rm_lower_bound, rm_upper_bound); + + proposal_sds_(param_index) = updated_sd; + } + + void increment_accepts(size_t param_index) { + validate_index(param_index); + acceptance_counts_[param_index]++; } void increment_iteration() { @@ -40,11 +49,15 @@ class AdaptiveProposal { } private: - arma::vec proposal_sds_; - int iterations_ = 0, - adaptation_window_; - double target_accept_ = 0.44; - bool adapting_ = true; + arma::vec proposal_sds_; + arma::ivec acceptance_counts_; + int iterations_ = 0, + adaptation_window_; + double target_accept_ = 0.44, + decay_rate_ = 0.75, + rm_lower_bound = 0.001, + rm_upper_bound = 2.0; + bool adapting_ = true; void validate_index(size_t index) const { if (index >= proposal_sds_.n_elem) { diff --git a/src/base_model.h b/src/base_model.h index 20cd4930..10f443dc 100644 --- a/src/base_model.h +++ b/src/base_model.h @@ -13,9 +13,9 @@ class BaseModel { virtual bool has_adaptive_mh() const { return false; } // Core methods (to be overridden by derived classes) - virtual double logp(const std::vector& parameters) = 0; + virtual double logp(const arma::vec& parameters) = 0; - virtual arma::vec gradient(const std::vector& parameters) { + virtual arma::vec gradient(const arma::vec& parameters) { if (!has_gradient()) { throw std::runtime_error("Gradient not implemented for this model"); } @@ -23,7 +23,7 @@ class BaseModel { } virtual std::pair logp_and_gradient( - const std::vector& parameters) { + const arma::vec& parameters) { if (!has_gradient()) { throw std::runtime_error("Gradient not implemented for this model"); } @@ -39,6 +39,10 @@ class BaseModel { throw std::runtime_error("get_vectorized_parameters method must be implemented in derived class"); } + virtual arma::ivec get_vectorized_indicator_parameters() { + throw std::runtime_error("get_vectorized_indicator_parameters method must be implemented in derived class"); + } + // Return dimensionality of the parameter space virtual size_t parameter_dimension() const = 0; diff --git a/src/chainResultNew.h b/src/chainResultNew.h index 5a6dd855..fa269a54 100644 --- a/src/chainResultNew.h +++ b/src/chainResultNew.h @@ -8,17 +8,17 @@ class ChainResultNew { public: ChainResultNew() {} - bool error; + bool error = false, + userInterrupt = false; std::string error_msg; - int chain_id; - bool userInterrupt; + int chain_id; - arma::mat samples; + arma::mat samples; - void reserve(size_t param_dim, size_t n_iter) { + void reserve(const size_t param_dim, const size_t n_iter) { samples.set_size(param_dim, n_iter); } - void store_sample(size_t iter, const arma::vec& sample) { + void store_sample(const size_t iter, const arma::vec& sample) { samples.col(iter) = sample; } diff --git a/src/ggm_model.cpp b/src/ggm_model.cpp index eee45bf7..0ffc19b2 100644 --- a/src/ggm_model.cpp +++ b/src/ggm_model.cpp @@ -10,10 +10,7 @@ double GGMModel::compute_inv_submatrix_i(const arma::mat& A, const size_t i, con void GGMModel::get_constants(size_t i, size_t j) { // TODO: helper function? - double logdet_omega = 0.0; - for (size_t i = 0; i < p_; i++) { - logdet_omega += std::log(phi_(i, i)); - } + double logdet_omega = get_log_det(phi_); double log_adj_omega_ii = logdet_omega + log(abs(inv_omega_(i, i))); double log_adj_omega_ij = logdet_omega + log(abs(inv_omega_(i, j))); @@ -22,7 +19,7 @@ void GGMModel::get_constants(size_t i, size_t j) { double inv_omega_sub_j1j1 = compute_inv_submatrix_i(inv_omega_, i, j, j); double log_abs_inv_omega_sub_jj = log_adj_omega_ii + log(abs(inv_omega_sub_j1j1)); - double Phi_q1q = (-1 * (2 * std::signbit(inv_omega_(i, j)) - 1)) * std::exp( + double Phi_q1q = (2 * std::signbit(inv_omega_(i, j)) - 1) * std::exp( (log_adj_omega_ij - (log_adj_omega_jj + log_abs_inv_omega_sub_jj) / 2) ); double Phi_q1q1 = std::exp((log_adj_omega_jj - log_abs_inv_omega_sub_jj) / 2); @@ -40,10 +37,100 @@ double GGMModel::R(const double x) const { if (x == 0) { return constants_[6]; } else { - return constants_[3] + std::pow((x - constants_[3]) / constants_[4], 2); + return constants_[5] + std::pow((x - constants_[3]) / constants_[4], 2); } } +double GGMModel::get_log_det(arma::mat triangular_A) const { + // assume A is an (upper) triangular cholesky factor + // returns the log determinant of A'A + + // TODO: should we just do + // log_det(val, sign, trimatu(A))? + return 2 * arma::accu(arma::log(triangular_A.diag())); +} + +double GGMModel::log_density_impl(const arma::mat& omega, const arma::mat& phi) const { + + double logdet_omega = get_log_det(phi); + // TODO: why not just dot(omega, suf_stat_)? + double trace_prod = arma::accu(omega % suf_stat_); + + double log_likelihood = n_ * (p_ * log(2 * arma::datum::pi) / 2 + logdet_omega / 2) - trace_prod / 2; + + return log_likelihood; +} + +double GGMModel::log_density_impl_edge(size_t i, size_t j) const { + + // this is the log likelihood ratio, not the full log likelihood like GGMModel::log_density_impl + + double Ui2 = omega_(i, j) - omega_prop_(i, j); + // only reached from R + // if (omega_(j, j) == omega_prop_(j, j)) { + // k = i; + // i = j; + // j = k; + // } + double Uj2 = (omega_(j, j) - omega_prop_(j, j)) / 2; + + + // W <- matrix(c(0, 1, 1, 0), 2, 2) + // U0 <- matrix(c(0, -1, Ui2, Uj2)) + // U <- matrix(0, nrow(aOmega), 2) + // U[c(i, j), 1] <- c(0, -1) + // U[c(i, j), 2] <- c(Ui2, Uj2) + // aOmega_prop - (aOmega + U %*% W %*% t(U)) + // det(aOmega_prop) - det(aOmega + U %*% W %*% t(U)) + // det(aOmega_prop) - det(W + t(U) %*% inv_aOmega %*% U) * det(W) * det(aOmega) + // below computes logdet(W + t(U) %*% inv_aOmega %*% U) directly (this is a 2x2 matrix) + + double cc11 = 0 + inv_omega_(j, j); + double cc12 = 1 - (inv_omega_(i, j) * Ui2 + inv_omega_(j, j) * Uj2); + double cc22 = 0 + Ui2 * Ui2 * inv_omega_(i, i) + 2 * Ui2 * Uj2 * inv_omega_(i, j) + Uj2 * Uj2 * inv_omega_(j, j); + + double logdet = std::log(std::abs(cc11 * cc22 - cc12 * cc12)); + // logdet - (logdet(aOmega_prop) - logdet(aOmega)) + + double trace_prod = -2 * (suf_stat_(j, j) * Uj2 + suf_stat_(i, j) * Ui2); + + // This function uses the fact that the determinant doesn't change during edge updates. + // double trace_prod = 0.0; + // // TODO: we only need one of the two lines below, but it's not entirely clear which one + // trace_prod += suf_stat_(j, j) * (omega_prop(j, j) - omega(j, j)); + // trace_prod += suf_stat_(i, i) * (omega_prop(i, i) - omega(i, i)); + // trace_prod += 2 * suf_stat_(i, j) * (omega_prop(i, j) - omega(i, j)); + // trace_prod - sum((aOmega_prop - aOmega) * SufStat) + + double log_likelihood_ratio = (n_ * logdet - trace_prod) / 2; + return log_likelihood_ratio; + +} + +double GGMModel::log_density_impl_diag(size_t j) const { + // same as above but for i == j, so Ui2 = 0 + double Uj2 = (omega_(j, j) - omega_prop_(j, j)) / 2; + + double cc11 = 0 + inv_omega_(j, j); + double cc12 = 1 - inv_omega_(j, j) * Uj2; + double cc22 = 0 + Uj2 * Uj2 * inv_omega_(j, j); + + double logdet = log(abs(cc11 * cc22 - cc12 * cc12)); + double trace_prod = -2 * suf_stat_(j, j) * Uj2; + + // This function uses the fact that the determinant doesn't change during edge updates. + // double trace_prod = 0.0; + // // TODO: we only need one of the two lines below, but it's not entirely clear which one + // trace_prod += suf_stat_(j, j) * (omega_prop(j, j) - omega(j, j)); + // trace_prod += suf_stat_(i, i) * (omega_prop(i, i) - omega(i, i)); + // trace_prod += 2 * suf_stat_(i, j) * (omega_prop(i, j) - omega(i, j)); + // trace_prod - sum((aOmega_prop - aOmega) * SufStat) + + double log_likelihood_ratio = (n_ * logdet - trace_prod) / 2; + return log_likelihood_ratio; + +} + void GGMModel::update_edge_parameter(size_t i, size_t j) { if (edge_indicators_(i, j) == 0) { @@ -67,13 +154,44 @@ void GGMModel::update_edge_parameter(size_t i, size_t j) { omega_prop_(j, i) = omega_prop_q1q; omega_prop_(j, j) = omega_prop_qq; - double ln_alpha = log_density(omega_prop_) - log_density(); + // Rcpp::Rcout << "i: " << i << ", j: " << j << + // ", proposed phi: " << phi_prop << + // ", proposal_sd omega_ij: " << proposal_sd << + // ", proposed omega_ij: " << omega_prop_q1q << + // ", proposed omega_jj: " << omega_prop_qq << std::endl; + // constants_.print(Rcpp::Rcout, "Constants:"); + // omega_prop_.print(Rcpp::Rcout, "Proposed omega:"); + + // arma::vec eigval = eig_sym(omega_prop_); + // if (arma::any(eigval <= 0)) { + // Rcpp::Rcout << "Warning: omega_prop_ is not positive definite for edge (" << i << ", " << j << ")" << std::endl; + + // Rcpp::Rcout << + // ", proposed phi: " << phi_prop << + // ", proposal_sd omega_ij: " << proposal_sd << + // ", proposed omega_ij: " << omega_prop_q1q << + // ", proposed omega_jj: " << omega_prop_qq << std::endl; + // constants_.print(Rcpp::Rcout, "Constants:"); + // omega_prop_.print(Rcpp::Rcout, "Proposed omega:"); + // omega_.print(Rcpp::Rcout, "Current omega:"); + // phi_.print(Rcpp::Rcout, "Phi:"); + // inv_omega_.print(Rcpp::Rcout, "Inv(Omega):"); + + // } + + // double ln_alpha = log_density(omega_prop_) - log_density(); + double ln_alpha = log_density_impl_edge(i, j); + + if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { + Rcpp::Rcout << "Warning: log density implementations do not match for edge (" << i << ", " << j << ")" << std::endl; + } + ln_alpha += R::dcauchy(omega_prop_(i, j), 0.0, 2.5, true); ln_alpha -= R::dcauchy(omega_(i, j), 0.0, 2.5, true); - double u = runif(rng_); - if (ln_alpha > log(u)) { + if (std::log(runif(rng_)) < ln_alpha) { // accept proposal + proposal_.increment_accepts(e); double omega_ij = omega_(i, j); double omega_jj = omega_(j, j); @@ -89,10 +207,10 @@ void GGMModel::update_edge_parameter(size_t i, size_t j) { arma::vec vf1 = arma::zeros(p_); arma::vec vf2 = arma::zeros(p_); - vf1[i] = v1[1]; - vf1[j] = v1[2]; - vf2[i] = v2[1]; - vf2[j] = v2[2]; + vf1[i] = v1[0]; + vf1[j] = v1[1]; + vf2[i] = v2[0]; + vf2[j] = v2[1]; // we now have // aOmega_prop - (aOmega + vf1 %*% t(vf2) + vf2 %*% t(vf1)) @@ -105,68 +223,24 @@ void GGMModel::update_edge_parameter(size_t i, size_t j) { // and also // aOmega_prop - (aOmega + cbind(vf1, vf2) %*% matrix(c(0, 1, 1, 0), 2, 2) %*% t(cbind(vf1, vf2))) - // update phi + // update phi (2x O(p^2)) cholesky_update(phi_, u1); cholesky_downdate(phi_, u2); - // update inverse - inv_omega_ = phi_.t() * phi_; + // update inverse (2x O(p^2)) + arma::inv(inv_phi_, arma::trimatu(phi_)); + inv_omega_ = inv_phi_ * inv_phi_.t(); } - double alpha = std::min(1.0, std::exp(ln_alpha)); - proposal_.update_proposal_sd(e, alpha); -} - -void GGMModel::do_one_mh_step() { - - // Update off-diagonals (upper triangle) - for (size_t i = 0; i < p_ - 1; ++i) { - for (size_t j = i + 1; j < p_; ++j) { - Rcpp::Rcout << "Updating edge parameter (" << i << ", " << j << ")" << std::endl; - update_edge_parameter(i, j); - } - } - - // Update diagonals - for (size_t i = 0; i < p_; ++i) { - Rcpp::Rcout << "Updating diagonal parameter " << i << std::endl; - update_diagonal_parameter(i); - } - - // if (edge_selection_) { - // for (size_t i = 0; i < p_ - 1; ++i) { - // for (size_t j = i + 1; j < p_; ++j) { - // update_edge_indicator_parameter_pair(i, j); - // } - // } - // } - proposal_.increment_iteration(); -} - -double GGMModel::log_density_impl(const arma::mat& omega, const arma::mat& phi) const { - double logdet_omega = 0.0; - for (size_t i = 0; i < p_; i++) { - logdet_omega += std::log(phi(i, i)); - } - - // TODO: does this allocate? - double trace_prod = arma::accu(omega % suf_stat_); - - double log_likelihood = n_ * (p_ * log(2 * arma::datum::pi) / 2 + logdet_omega / 2) - trace_prod / 2; - - return log_likelihood; + proposal_.update_proposal_sd(e); } void GGMModel::update_diagonal_parameter(size_t i) { // Implementation of diagonal parameter update // 1-3) from before - double logdet_omega = 0.0; - for (size_t i = 0; i < p_; i++) { - logdet_omega += std::log(phi_(i, i)); - } - + double logdet_omega = get_log_det(phi_); double logdet_omega_sub_ii = logdet_omega + std::log(inv_omega_(i, i)); size_t e = i * (i + 1) / 2 + i; // parameter index in vectorized form @@ -179,37 +253,277 @@ void GGMModel::update_diagonal_parameter(size_t i) { omega_prop_ = omega_; omega_prop_(i, i) = omega_(i, i) - std::exp(theta_curr) * std::exp(theta_curr) + std::exp(theta_prop) * std::exp(theta_prop); + // Rcpp::Rcout << "i: " << i << + // ", current theta: " << theta_curr << + // ", proposed theta: " << theta_prop << + // ", proposal_sd: " << proposal_sd << std::endl; + // omega_prop_.print(Rcpp::Rcout, "Proposed omega:"); + // 5) Acceptance ratio - double ln_alpha = log_density(omega_prop_) - log_density(); + // double ln_alpha = log_density(omega_prop_) - log_density(); + double ln_alpha = log_density_impl_diag(i); + if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { + Rcpp::Rcout << "Warning: log density implementations do not match for diag (" << i << ", " << i << ")" << std::endl; + } + ln_alpha += R::dgamma(exp(theta_prop), 1.0, 1.0, true); ln_alpha -= R::dgamma(exp(theta_curr), 1.0, 1.0, true); ln_alpha += theta_prop - theta_curr; // Jacobian adjustment ? + if (std::log(runif(rng_)) < ln_alpha) { - if (log(runif(rng_)) < ln_alpha) { + proposal_.increment_accepts(e); double omega_ii = omega_(i, i); arma::vec u(p_, arma::fill::zeros); double delta = omega_ii - omega_prop_(i, i); bool s = delta > 0; - u(i) = sqrt(abs(delta)); + u(i) = std::sqrt(std::abs(delta)); omega_(i, i) = omega_prop_(i, i); - if (!s) - cholesky_update(phi_, u); - else + if (s) cholesky_downdate(phi_, u); + else + cholesky_update(phi_, u); + + // update inverse (2x O(p^2)) + arma::inv(inv_phi_, arma::trimatu(phi_)); + inv_omega_ = inv_phi_ * inv_phi_.t(); - inv_omega_ = phi_.t() * phi_; } - double alpha = std::min(1.0, std::exp(ln_alpha)); - proposal_.update_proposal_sd(e, alpha); + proposal_.update_proposal_sd(e); } void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { - // Implementation of edge indicator parameter pair update + + size_t e = i * (i + 1) / 2 + j; // parameter index in vectorized form + double proposal_sd = proposal_.get_proposal_sd(e); + + if (edge_indicators_(i, j) == 1) { + // Propose to turn OFF the edge + omega_prop_ = omega_; + omega_prop_(i, j) = 0.0; + omega_prop_(j, i) = 0.0; + + // Update diagonal using R function with omega_ij = 0 + get_constants(i, j); + omega_prop_(j, j) = R(0.0); + + // double ln_alpha = log_density(omega_prop_) - log_density(); + double ln_alpha = log_density_impl_edge(i, j); + if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { + Rcpp::Rcout << "Warning: log density indicator implementations do not match for edge (" << i << ", " << j << ")" << std::endl; + } + + ln_alpha += std::log(1.0 - prior_inclusion_prob_(i, j)) - std::log(prior_inclusion_prob_(i, j)); + + ln_alpha += R::dnorm(omega_(i, j) / constants_[4], 0.0, proposal_sd, true) - std::log(constants_[4]); + ln_alpha -= R::dcauchy(omega_(i, j), 0.0, 2.5, true); + + if (std::log(runif(rng_)) < ln_alpha) { + + // Store old values for Cholesky update + double omega_ij_old = omega_(i, j); + double omega_jj_old = omega_(j, j); + + // Update omega + omega_(i, j) = 0.0; + omega_(j, i) = 0.0; + omega_(j, j) = omega_prop_(j, j); + + // Update edge indicator + edge_indicators_(i, j) = 0; + edge_indicators_(j, i) = 0; + + // Cholesky update vectors + arma::vec v1 = {0, -1}; + arma::vec v2 = {omega_ij_old - 0.0, (omega_jj_old - omega_(j, j)) / 2}; + + arma::vec vf1 = arma::zeros(p_); + arma::vec vf2 = arma::zeros(p_); + vf1[i] = v1[0]; + vf1[j] = v1[1]; + vf2[i] = v2[0]; + vf2[j] = v2[1]; + + arma::vec u1 = (vf1 + vf2) / sqrt(2); + arma::vec u2 = (vf1 - vf2) / sqrt(2); + + // Update Cholesky factor + cholesky_update(phi_, u1); + cholesky_downdate(phi_, u2); + + // Update inverse + arma::inv(inv_phi_, arma::trimatu(phi_)); + inv_omega_ = inv_phi_ * inv_phi_.t(); + } + + } else { + // Propose to turn ON the edge + double epsilon = rnorm(rng_, 0.0, proposal_sd); + + // Get constants for current state (with edge OFF) + get_constants(i, j); + double omega_prop_ij = constants_[4] * epsilon; + double omega_prop_jj = R(omega_prop_ij); + + omega_prop_ = omega_; + omega_prop_(i, j) = omega_prop_ij; + omega_prop_(j, i) = omega_prop_ij; + omega_prop_(j, j) = omega_prop_jj; + + // double ln_alpha = log_density(omega_prop_) - log_density(); + double ln_alpha = log_density_impl_edge(i, j); + if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { + Rcpp::Rcout << "Warning: log density indicator implementations do not match for edge (" << i << ", " << j << ")" << std::endl; + } + ln_alpha += std::log(prior_inclusion_prob_(i, j)) - std::log(1.0 - prior_inclusion_prob_(i, j)); + + // Prior change: add slab (Cauchy prior) + ln_alpha += R::dcauchy(omega_prop_ij, 0.0, 2.5, true); + + // Proposal term: proposed edge value given it was generated from truncated normal + ln_alpha -= R::dnorm(omega_prop_ij / constants_[4], 0.0, proposal_sd, true) - std::log(constants_[4]); + + // TODO: this can be factored out? + if (std::log(runif(rng_)) < ln_alpha) { + // Accept: turn ON the edge + proposal_.increment_accepts(e); + + // Store old values for Cholesky update + double omega_ij_old = omega_(i, j); + double omega_jj_old = omega_(j, j); + + // Update omega + omega_(i, j) = omega_prop_ij; + omega_(j, i) = omega_prop_ij; + omega_(j, j) = omega_prop_jj; + + // Update edge indicator + edge_indicators_(i, j) = 1; + edge_indicators_(j, i) = 1; + + // Cholesky update vectors + arma::vec v1 = {0, -1}; + arma::vec v2 = {omega_ij_old - omega_(i, j), (omega_jj_old - omega_(j, j)) / 2}; + + arma::vec vf1 = arma::zeros(p_); + arma::vec vf2 = arma::zeros(p_); + vf1[i] = v1[0]; + vf1[j] = v1[1]; + vf2[i] = v2[0]; + vf2[j] = v2[1]; + + arma::vec u1 = (vf1 + vf2) / sqrt(2); + arma::vec u2 = (vf1 - vf2) / sqrt(2); + + // Update Cholesky factor + cholesky_update(phi_, u1); + cholesky_downdate(phi_, u2); + + // Update inverse + arma::inv(inv_phi_, arma::trimatu(phi_)); + inv_omega_ = inv_phi_ * inv_phi_.t(); + } + } +} + +void GGMModel::do_one_mh_step() { + + // Update off-diagonals (upper triangle) + for (size_t i = 0; i < p_ - 1; ++i) { + for (size_t j = i + 1; j < p_; ++j) { + // Rcpp::Rcout << "Updating edge parameter (" << i << ", " << j << ")" << std::endl; + update_edge_parameter(i, j); + // if (!arma:: approx_equal(omega_ * inv_omega_, arma::eye(p_, p_), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Omega * Inv(Omega) not equal to identity after updating edge (" << i << ", " << j << ")" << std::endl; + // (omega_ * inv_omega_).print(Rcpp::Rcout, "Omega * Inv(Omega):"); + // (phi_ * inv_phi_).print(Rcpp::Rcout, "Phi * Inv(Phi):"); + // } + // if (!arma:: approx_equal(omega_, phi_.t() * phi_, "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Omega not equal to Phi.t() * Phi after updating edge (" << i << ", " << j << ")" << std::endl; + // omega_.print(Rcpp::Rcout, "Omega:"); + // phi_.print(Rcpp::Rcout, "Phi:"); + // } + // if (!arma:: approx_equal(phi_ * inv_phi_, arma::eye(p_, p_), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Phi * Inv(Phi) not equal to identity after updating edge (" << i << ", " << j << ")" << std::endl; + // (omega_ * inv_omega_).print(Rcpp::Rcout, "Omega * Inv(Omega):"); + // (phi_ * inv_phi_).print(Rcpp::Rcout, "Phi * Inv(Phi):"); + // } + // if (!arma:: approx_equal(inv_omega_, inv_phi_ * inv_phi_.t(), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Inv(Omega) not equal to Inv(Phi) * Inv(Phi).t() after updating edge (" << i << ", " << j << ")" << std::endl; + // omega_.print(Rcpp::Rcout, "Omega:"); + // phi_.print(Rcpp::Rcout, "Phi:"); + // inv_omega_.print(Rcpp::Rcout, "Inv(Omega):"); + // inv_phi_.print(Rcpp::Rcout, "Inv(Phi):"); + // } + } + } + + // Update diagonals + for (size_t i = 0; i < p_; ++i) { + // Rcpp::Rcout << "Updating diagonal parameter " << i << std::endl; + update_diagonal_parameter(i); + + // if (!arma:: approx_equal(omega_ * inv_omega_, arma::eye(p_, p_), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Omega * Inv(Omega) not equal to identity after updating diagonal " << i << std::endl; + // (omega_ * inv_omega_).print(Rcpp::Rcout, "Omega * Inv(Omega):"); + // (phi_ * inv_phi_).print(Rcpp::Rcout, "Phi * Inv(Phi):"); + // } + // if (!arma:: approx_equal(omega_, phi_.t() * phi_, "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Omega not equal to Phi.t() * Phi after updating diagonal " << i << std::endl; + // omega_.print(Rcpp::Rcout, "Omega:"); + // phi_.print(Rcpp::Rcout, "Phi:"); + // } + // if (!arma:: approx_equal(phi_ * inv_phi_, arma::eye(p_, p_), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Phi * Inv(Phi) not equal to identity after updating diagonal " << i << std::endl; + // (omega_ * inv_omega_).print(Rcpp::Rcout, "Omega * Inv(Omega):"); + // (phi_ * inv_phi_).print(Rcpp::Rcout, "Phi * Inv(Phi):"); + // } + // if (!arma:: approx_equal(inv_omega_, inv_phi_ * inv_phi_.t(), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Inv(Omega) not equal to Inv(Phi) * Inv(Phi).t() after updating diagonal " << i << std::endl; + // omega_.print(Rcpp::Rcout, "Omega:"); + // phi_.print(Rcpp::Rcout, "Phi:"); + // inv_omega_.print(Rcpp::Rcout, "Inv(Omega):"); + // inv_phi_.print(Rcpp::Rcout, "Inv(Phi):"); + // } + } + + if (edge_selection_) { + for (size_t i = 0; i < p_ - 1; ++i) { + for (size_t j = i + 1; j < p_; ++j) { + // Rcpp::Rcout << "Between model move for edge (" << i << ", " << j << ")" << std::endl; + update_edge_indicator_parameter_pair(i, j); + // if (!arma:: approx_equal(omega_ * inv_omega_, arma::eye(p_, p_), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Omega * Inv(Omega) not equal to identity after updating edge (" << i << ", " << j << ")" << std::endl; + // (omega_ * inv_omega_).print(Rcpp::Rcout, "Omega * Inv(Omega):"); + // (phi_ * inv_phi_).print(Rcpp::Rcout, "Phi * Inv(Phi):"); + // } + // if (!arma:: approx_equal(omega_, phi_.t() * phi_, "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Omega not equal to Phi.t() * Phi after updating edge (" << i << ", " << j << ")" << std::endl; + // omega_.print(Rcpp::Rcout, "Omega:"); + // phi_.print(Rcpp::Rcout, "Phi:"); + // } + // if (!arma:: approx_equal(phi_ * inv_phi_, arma::eye(p_, p_), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Phi * Inv(Phi) not equal to identity after updating edge (" << i << ", " << j << ")" << std::endl; + // (omega_ * inv_omega_).print(Rcpp::Rcout, "Omega * Inv(Omega):"); + // (phi_ * inv_phi_).print(Rcpp::Rcout, "Phi * Inv(Phi):"); + // } + // if (!arma:: approx_equal(inv_omega_, inv_phi_ * inv_phi_.t(), "absdiff", 1e-6)) { + // Rcpp::Rcout << "Warning: Inv(Omega) not equal to Inv(Phi) * Inv(Phi).t() after updating edge (" << i << ", " << j << ")" << std::endl; + // omega_.print(Rcpp::Rcout, "Omega:"); + // phi_.print(Rcpp::Rcout, "Phi:"); + // inv_omega_.print(Rcpp::Rcout, "Inv(Omega):"); + // inv_phi_.print(Rcpp::Rcout, "Inv(Phi):"); + // } + } + } + } + + // could also be called in the main MCMC loop + proposal_.increment_iteration(); } diff --git a/src/ggm_model.h b/src/ggm_model.h index 05081db0..48158588 100644 --- a/src/ggm_model.h +++ b/src/ggm_model.h @@ -22,9 +22,12 @@ class GGMModel : public BaseModel { proposal_(AdaptiveProposal(dim_, 500)), omega_(arma::eye(p_, p_)), phi_(arma::eye(p_, p_)), + inv_phi_(arma::eye(p_, p_)), inv_omega_(arma::eye(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) {} @@ -38,9 +41,11 @@ class GGMModel : public BaseModel { edge_selection_(other.edge_selection_), omega_(other.omega_), phi_(other.phi_), + inv_phi_(other.inv_phi_), inv_omega_(other.inv_omega_), edge_indicators_(other.edge_indicators_), vectorized_parameters_(other.vectorized_parameters_), + vectorized_indicator_parameters_(other.vectorized_indicator_parameters_), proposal_(other.proposal_), rng_(other.rng_), omega_prop_(other.omega_prop_), @@ -55,10 +60,10 @@ class GGMModel : public BaseModel { proposal_ = proposal; } - virtual bool has_gradient() const { return false; } - virtual bool has_adaptive_mh() const override { return true; } + bool has_gradient() const { return false; } + bool has_adaptive_mh() const override { return true; } - double logp(const std::vector& parameters) override { + double logp(const arma::vec& parameters) override { // Implement log probability computation return 0.0; } @@ -70,7 +75,6 @@ class GGMModel : public BaseModel { void do_one_mh_step() override; size_t parameter_dimension() const override { - Rcpp::Rcout << "GGMModel::parameter_dimension() returning: " << dim_ << std::endl; return dim_; } @@ -79,33 +83,47 @@ class GGMModel : public BaseModel { } arma::vec get_vectorized_parameters() override { + // upper triangle of omega_ size_t e = 0; for (size_t j = 0; j < p_; ++j) { for (size_t i = 0; i <= j; ++i) { vectorized_parameters_(e) = omega_(i, j); - e++; + ++e; } } return vectorized_parameters_; } + arma::ivec get_vectorized_indicator_parameters() override { + // upper triangle of omega_ + size_t e = 0; + for (size_t j = 0; j < p_; ++j) { + for (size_t i = 0; i <= j; ++i) { + vectorized_indicator_parameters_(e) = edge_indicators_(i, j); + ++e; + } + } + return vectorized_indicator_parameters_; + } + std::unique_ptr clone() const override { return std::make_unique(*this); // uses copy constructor } private: // data - size_t dim_; - arma::mat suf_stat_; size_t n_; size_t p_; + size_t dim_; + arma::mat suf_stat_; arma::mat prior_inclusion_prob_; - bool edge_selection_ = true; + bool edge_selection_; // parameters - arma::mat omega_, phi_, inv_omega_; + arma::mat omega_, phi_, inv_phi_, inv_omega_; arma::imat edge_indicators_; arma::vec vectorized_parameters_; + arma::ivec vectorized_indicator_parameters_; AdaptiveProposal proposal_; @@ -126,7 +144,9 @@ class GGMModel : public BaseModel { double R(const double x) const; double log_density_impl(const arma::mat& omega, const arma::mat& phi) const; - + double log_density_impl_edge(size_t i, size_t j) const; + double log_density_impl_diag(size_t j) const; + double get_log_det(arma::mat triangular_A) const; // double find_reasonable_step_size_edge(const arma::mat& omega, size_t i, size_t j); // double find_reasonable_step_size_diag(const arma::mat& omega, size_t i); // double edge_log_ratio(const arma::mat& omega, size_t i, size_t j, double proposal); diff --git a/src/sample_ggm.cpp b/src/sample_ggm.cpp index 0d9de749..8ca13568 100644 --- a/src/sample_ggm.cpp +++ b/src/sample_ggm.cpp @@ -17,6 +17,7 @@ void run_mcmc_sampler_single_thread( ProgressManager& pm ) { + chain_result.chain_id = chain_id + 1; size_t i = 0; for (size_t iter = 0; iter < no_iter + no_warmup; ++iter) { @@ -25,7 +26,7 @@ void run_mcmc_sampler_single_thread( if (iter >= no_warmup) { chain_result.store_sample(i, model.get_vectorized_parameters()); - i++; + ++i; } pm.update(chain_id); @@ -64,8 +65,6 @@ struct GGMChainRunner : public RcppParallel::Worker { for (std::size_t i = begin; i < end; ++i) { ChainResultNew& chain_result = results_[i]; - chain_result.chain_id = static_cast(i + 1); - chain_result.error = false; BaseModel& model = *models_[i]; model.set_seed(seed_ + i); try { @@ -124,12 +123,17 @@ std::vector run_mcmc_sampler( models.push_back(model.clone()); // deep copy via virtual clone } run_mcmc_sampler_threaded(results, models, no_iter, no_warmup, seed, no_threads, pm); + } else { + model.set_seed(seed); Rcpp::Rcout << "Running single-threaded MCMC sampling..." << std::endl; + // TODO: this is actually not correct, each chain should have its own model object + // now chain 2 continues from chain 1 state for (size_t c = 0; c < no_chains; ++c) { run_mcmc_sampler_single_thread(results[c], model, no_iter, no_warmup, c, pm); } + } return results; } @@ -164,15 +168,13 @@ Rcpp::List sample_ggm( const bool edge_selection, const int seed, const int no_threads, - int progress_type + const int progress_type ) { // should be done dynamically // also adaptation method should be specified differently GGMModel model(X, prior_inclusion_prob, initial_edge_indicators, edge_selection); - Rcpp::Rcout << "GGMModel::parameter_dimension() returning: " << model.parameter_dimension() << std::endl; - ProgressManager pm(no_chains, no_iter, no_warmup, 50, progress_type); std::vector output = run_mcmc_sampler(model, no_iter, no_warmup, no_chains, seed, no_threads, pm); diff --git a/test_ggm.R b/test_ggm.R index 60c5ef0f..5d0ca2c4 100644 --- a/test_ggm.R +++ b/test_ggm.R @@ -1,7 +1,7 @@ library(bgms) # Dimension and true precision -p <- 10 +p <- 50 adj <- matrix(0, nrow = p, ncol = p) adj[lower.tri(adj)] <- rbinom(p * (p - 1) / 2, size = 1, prob = 0.3) @@ -12,7 +12,7 @@ Sigma <- solve(Omega) zapsmall(Omega) # Data -n <- 1e2 +n <- 1e3 x <- mvtnorm::rmvnorm(n = n, mean = rep(0, p), sigma = Sigma) @@ -23,15 +23,55 @@ sampling_results <- bgms:::sample_ggm( X = x, prior_inclusion_prob = matrix(.5, p, p), initial_edge_indicators = adj, - no_iter = 4000, - no_warmup = 4000, + no_iter = 500, + no_warmup = 500, no_chains = 3, edge_selection = FALSE, no_threads = 1, seed = 123, - progress_type = 2 + progress_type = 1 ) +true_values <- zapsmall(Omega[upper.tri(Omega, TRUE)]) +posterior_means <- rowMeans(sampling_results[[2]]$samples) + +plot(true_values, posterior_means) +abline(0, 1) + +sampling_results2 <- bgms:::sample_ggm( + X = x, + prior_inclusion_prob = matrix(.5, p, p), + initial_edge_indicators = adj, + no_iter = 500, + no_warmup = 500, + no_chains = 3, + edge_selection = TRUE, + no_threads = 1, + seed = 123, + progress_type = 1 +) + +true_values <- zapsmall(Omega[upper.tri(Omega, TRUE)]) +posterior_means <- rowMeans(sampling_results2[[2]]$samples) + +plot(true_values, posterior_means) +abline(0, 1) + +plot(posterior_means, rowMeans(sampling_results2[[2]]$samples != 0)) + + +mmm <- matrix(c( + 1.6735, 0, 0, 0, 0, + 0, 1.0000, 0, 0, -3.4524, + 0, 0, 1.0000, 0, 0, + 0, 0, 0, 1.0000, 0, + 0, -3.4524, 0, 0, 9.6674 +), p, p) +mmm +chol(mmm) +base::isSymmetric(mmm) +eigen(mmm) + profvis::profvis({ sampling_results <- bgm_gaussian( x = x, From 81c651195db5fb8e129c0ca96819bfb1c4e3c6eb Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Thu, 13 Nov 2025 09:16:24 +0100 Subject: [PATCH 3/5] can be tested --- src/ggm_model.cpp | 60 ++++++++++++++++++++++++++++++----------------- test_ggm.R | 3 ++- 2 files changed, 41 insertions(+), 22 deletions(-) diff --git a/src/ggm_model.cpp b/src/ggm_model.cpp index 0ffc19b2..db148fba 100644 --- a/src/ggm_model.cpp +++ b/src/ggm_model.cpp @@ -94,14 +94,6 @@ double GGMModel::log_density_impl_edge(size_t i, size_t j) const { double trace_prod = -2 * (suf_stat_(j, j) * Uj2 + suf_stat_(i, j) * Ui2); - // This function uses the fact that the determinant doesn't change during edge updates. - // double trace_prod = 0.0; - // // TODO: we only need one of the two lines below, but it's not entirely clear which one - // trace_prod += suf_stat_(j, j) * (omega_prop(j, j) - omega(j, j)); - // trace_prod += suf_stat_(i, i) * (omega_prop(i, i) - omega(i, i)); - // trace_prod += 2 * suf_stat_(i, j) * (omega_prop(i, j) - omega(i, j)); - // trace_prod - sum((aOmega_prop - aOmega) * SufStat) - double log_likelihood_ratio = (n_ * logdet - trace_prod) / 2; return log_likelihood_ratio; @@ -115,7 +107,7 @@ double GGMModel::log_density_impl_diag(size_t j) const { double cc12 = 1 - inv_omega_(j, j) * Uj2; double cc22 = 0 + Uj2 * Uj2 * inv_omega_(j, j); - double logdet = log(abs(cc11 * cc22 - cc12 * cc12)); + double logdet = std::log(std::abs(cc11 * cc22 - cc12 * cc12)); double trace_prod = -2 * suf_stat_(j, j) * Uj2; // This function uses the fact that the determinant doesn't change during edge updates. @@ -182,9 +174,15 @@ void GGMModel::update_edge_parameter(size_t i, size_t j) { // double ln_alpha = log_density(omega_prop_) - log_density(); double ln_alpha = log_density_impl_edge(i, j); - if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { - Rcpp::Rcout << "Warning: log density implementations do not match for edge (" << i << ", " << j << ")" << std::endl; - } + // { + // double ln_alpha_ref = log_density(omega_prop_) - log_density(); + // if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) { + // Rcpp::Rcout << "Warning: log density implementations do not match for edge (" << i << ", " << j << ")" << std::endl; + // omega_.print(Rcpp::Rcout, "Current omega:"); + // omega_prop_.print(Rcpp::Rcout, "Proposed omega:"); + // Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl; + // } + // } ln_alpha += R::dcauchy(omega_prop_(i, j), 0.0, 2.5, true); ln_alpha -= R::dcauchy(omega_(i, j), 0.0, 2.5, true); @@ -262,9 +260,16 @@ void GGMModel::update_diagonal_parameter(size_t i) { // 5) Acceptance ratio // double ln_alpha = log_density(omega_prop_) - log_density(); double ln_alpha = log_density_impl_diag(i); - if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { - Rcpp::Rcout << "Warning: log density implementations do not match for diag (" << i << ", " << i << ")" << std::endl; - } + // { + // double ln_alpha_ref = log_density(omega_prop_) - log_density(); + // if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) { + // Rcpp::Rcout << "Warning: log density implementations do not match for diag (" << i << ", " << i << ")" << std::endl; + // // omega_.print(Rcpp::Rcout, "Current omega:"); + // // omega_prop_.print(Rcpp::Rcout, "Proposed omega:"); + // Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl; + // Rcpp::Rcout << "1e4 * diff: " << 10000 * (ln_alpha - ln_alpha_ref) << std::endl; + // } + // } ln_alpha += R::dgamma(exp(theta_prop), 1.0, 1.0, true); ln_alpha -= R::dgamma(exp(theta_curr), 1.0, 1.0, true); @@ -315,9 +320,16 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { // double ln_alpha = log_density(omega_prop_) - log_density(); double ln_alpha = log_density_impl_edge(i, j); - if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { - Rcpp::Rcout << "Warning: log density indicator implementations do not match for edge (" << i << ", " << j << ")" << std::endl; - } + // { + // double ln_alpha_ref = log_density(omega_prop_) - log_density(); + // if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) { + // Rcpp::Rcout << "Warning: log density implementations do not match for edge indicator (" << i << ", " << j << ")" << std::endl; + // omega_.print(Rcpp::Rcout, "Current omega:"); + // omega_prop_.print(Rcpp::Rcout, "Proposed omega:"); + // Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl; + // } + // } + ln_alpha += std::log(1.0 - prior_inclusion_prob_(i, j)) - std::log(prior_inclusion_prob_(i, j)); @@ -378,9 +390,15 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { // double ln_alpha = log_density(omega_prop_) - log_density(); double ln_alpha = log_density_impl_edge(i, j); - if (std::abs(ln_alpha - (log_density(omega_prop_) - log_density())) > 1e-6) { - Rcpp::Rcout << "Warning: log density indicator implementations do not match for edge (" << i << ", " << j << ")" << std::endl; - } + // { + // double ln_alpha_ref = log_density(omega_prop_) - log_density(); + // if (std::abs(ln_alpha - ln_alpha_ref) > 1e-6) { + // Rcpp::Rcout << "Warning: log density implementations do not match for edge indicator (" << i << ", " << j << ")" << std::endl; + // omega_.print(Rcpp::Rcout, "Current omega:"); + // omega_prop_.print(Rcpp::Rcout, "Proposed omega:"); + // Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl; + // } + // } ln_alpha += std::log(prior_inclusion_prob_(i, j)) - std::log(1.0 - prior_inclusion_prob_(i, j)); // Prior change: add slab (Cauchy prior) diff --git a/test_ggm.R b/test_ggm.R index 5d0ca2c4..5f1995ca 100644 --- a/test_ggm.R +++ b/test_ggm.R @@ -1,7 +1,7 @@ library(bgms) # Dimension and true precision -p <- 50 +p <- 10 adj <- matrix(0, nrow = p, ncol = p) adj[lower.tri(adj)] <- rbinom(p * (p - 1) / 2, size = 1, prob = 0.3) @@ -34,6 +34,7 @@ sampling_results <- bgms:::sample_ggm( true_values <- zapsmall(Omega[upper.tri(Omega, TRUE)]) posterior_means <- rowMeans(sampling_results[[2]]$samples) +cbind(true_values, posterior_means) plot(true_values, posterior_means) abline(0, 1) From 447c64bea3dd183ee704b79cc43dcce40016f3db Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Wed, 19 Nov 2025 12:55:17 +0100 Subject: [PATCH 4/5] reduce allocations and factor out common cholesky stuff --- src/ggm_model.cpp | 219 ++++++++++++++++++++++++++++++---------------- src/ggm_model.h | 9 ++ 2 files changed, 151 insertions(+), 77 deletions(-) diff --git a/src/ggm_model.cpp b/src/ggm_model.cpp index db148fba..9d475176 100644 --- a/src/ggm_model.cpp +++ b/src/ggm_model.cpp @@ -12,13 +12,12 @@ void GGMModel::get_constants(size_t i, size_t j) { // TODO: helper function? double logdet_omega = get_log_det(phi_); - double log_adj_omega_ii = logdet_omega + log(abs(inv_omega_(i, i))); - double log_adj_omega_ij = logdet_omega + log(abs(inv_omega_(i, j))); - double log_adj_omega_jj = logdet_omega + log(abs(inv_omega_(j, j))); + double log_adj_omega_ii = logdet_omega + std::log(std::abs(inv_omega_(i, i))); + double log_adj_omega_ij = logdet_omega + std::log(std::abs(inv_omega_(i, j))); + double log_adj_omega_jj = logdet_omega + std::log(std::abs(inv_omega_(j, j))); double inv_omega_sub_j1j1 = compute_inv_submatrix_i(inv_omega_, i, j, j); - double log_abs_inv_omega_sub_jj = log_adj_omega_ii + log(abs(inv_omega_sub_j1j1)); - + double log_abs_inv_omega_sub_jj = log_adj_omega_ii + std::log(std::abs(inv_omega_sub_j1j1)); double Phi_q1q = (2 * std::signbit(inv_omega_(i, j)) - 1) * std::exp( (log_adj_omega_ij - (log_adj_omega_jj + log_abs_inv_omega_sub_jj) / 2) ); @@ -191,49 +190,89 @@ void GGMModel::update_edge_parameter(size_t i, size_t j) { // accept proposal proposal_.increment_accepts(e); - double omega_ij = omega_(i, j); - double omega_jj = omega_(j, j); + double omega_ij_old = omega_(i, j); + double omega_jj_old = omega_(j, j); + omega_(i, j) = omega_prop_q1q; omega_(j, i) = omega_prop_q1q; omega_(j, j) = omega_prop_qq; - // TODO: preallocate? - // find v for low rank update - arma::vec v1 = {0, -1}; - arma::vec v2 = {omega_ij - omega_prop_(i, j), (omega_jj - omega_prop_(j, j)) / 2}; + cholesky_update_after_edge(omega_ij_old, omega_jj_old, i, j); + + // // TODO: preallocate? + // // find v for low rank update + // arma::vec v1 = {0, -1}; + // arma::vec v2 = {omega_ij - omega_prop_(i, j), (omega_jj - omega_prop_(j, j)) / 2}; - arma::vec vf1 = arma::zeros(p_); - arma::vec vf2 = arma::zeros(p_); - vf1[i] = v1[0]; - vf1[j] = v1[1]; - vf2[i] = v2[0]; - vf2[j] = v2[1]; + // arma::vec vf1 = arma::zeros(p_); + // arma::vec vf2 = arma::zeros(p_); + // vf1[i] = v1[0]; + // vf1[j] = v1[1]; + // vf2[i] = v2[0]; + // vf2[j] = v2[1]; - // we now have - // aOmega_prop - (aOmega + vf1 %*% t(vf2) + vf2 %*% t(vf1)) + // // we now have + // // aOmega_prop - (aOmega + vf1 %*% t(vf2) + vf2 %*% t(vf1)) - arma::vec u1 = (vf1 + vf2) / sqrt(2); - arma::vec u2 = (vf1 - vf2) / sqrt(2); + // arma::vec u1 = (vf1 + vf2) / sqrt(2); + // arma::vec u2 = (vf1 - vf2) / sqrt(2); - // we now have - // omega_prop_ - (aOmega + u1 %*% t(u1) - u2 %*% t(u2)) - // and also - // aOmega_prop - (aOmega + cbind(vf1, vf2) %*% matrix(c(0, 1, 1, 0), 2, 2) %*% t(cbind(vf1, vf2))) + // // we now have + // // omega_prop_ - (aOmega + u1 %*% t(u1) - u2 %*% t(u2)) + // // and also + // // aOmega_prop - (aOmega + cbind(vf1, vf2) %*% matrix(c(0, 1, 1, 0), 2, 2) %*% t(cbind(vf1, vf2))) - // update phi (2x O(p^2)) - cholesky_update(phi_, u1); - cholesky_downdate(phi_, u2); + // // update phi (2x O(p^2)) + // cholesky_update(phi_, u1); + // cholesky_downdate(phi_, u2); - // update inverse (2x O(p^2)) - arma::inv(inv_phi_, arma::trimatu(phi_)); - inv_omega_ = inv_phi_ * inv_phi_.t(); + // // update inverse (2x O(p^2)) + // arma::inv(inv_phi_, arma::trimatu(phi_)); + // inv_omega_ = inv_phi_ * inv_phi_.t(); } proposal_.update_proposal_sd(e); } +void GGMModel::cholesky_update_after_edge(double omega_ij_old, double omega_jj_old, size_t i, size_t j) +{ + + v2_[0] = omega_ij_old - omega_prop_(i, j); + v2_[1] = (omega_jj_old - omega_prop_(j, j)) / 2; + + vf1_[i] = v1_[0]; + vf1_[j] = v1_[1]; + vf2_[i] = v2_[0]; + vf2_[j] = v2_[1]; + + // we now have + // aOmega_prop - (aOmega + vf1 %*% t(vf2) + vf2 %*% t(vf1)) + + u1_ = (vf1_ + vf2_) / sqrt(2); + u2_ = (vf1_ - vf2_) / sqrt(2); + + // we now have + // omega_prop_ - (aOmega + u1 %*% t(u1) - u2 %*% t(u2)) + // and also + // aOmega_prop - (aOmega + cbind(vf1, vf2) %*% matrix(c(0, 1, 1, 0), 2, 2) %*% t(cbind(vf1, vf2))) + + // update phi (2x O(p^2)) + cholesky_update(phi_, u1_); + cholesky_downdate(phi_, u2_); + + // update inverse (2x O(p^2)) + arma::inv(inv_phi_, arma::trimatu(phi_)); + inv_omega_ = inv_phi_ * inv_phi_.t(); + + // reset for next iteration + vf1_[i] = 0.0; + vf1_[j] = 0.0; + vf2_[i] = 0.0; + vf2_[j] = 0.0; + +} void GGMModel::update_diagonal_parameter(size_t i) { // Implementation of diagonal parameter update @@ -280,22 +319,24 @@ void GGMModel::update_diagonal_parameter(size_t i) { proposal_.increment_accepts(e); double omega_ii = omega_(i, i); + omega_(i, i) = omega_prop_(i, i); - arma::vec u(p_, arma::fill::zeros); - double delta = omega_ii - omega_prop_(i, i); - bool s = delta > 0; - u(i) = std::sqrt(std::abs(delta)); + cholesky_update_after_diag(omega_ii, i); - omega_(i, i) = omega_prop_(i, i); + // arma::vec u(p_, arma::fill::zeros); + // double delta = omega_ii - omega_prop_(i, i); + // bool s = delta > 0; + // u(i) = std::sqrt(std::abs(delta)); - if (s) - cholesky_downdate(phi_, u); - else - cholesky_update(phi_, u); - // update inverse (2x O(p^2)) - arma::inv(inv_phi_, arma::trimatu(phi_)); - inv_omega_ = inv_phi_ * inv_phi_.t(); + // if (s) + // cholesky_downdate(phi_, u); + // else + // cholesky_update(phi_, u); + + // // update inverse (2x O(p^2)) + // arma::inv(inv_phi_, arma::trimatu(phi_)); + // inv_omega_ = inv_phi_ * inv_phi_.t(); } @@ -303,6 +344,28 @@ void GGMModel::update_diagonal_parameter(size_t i) { proposal_.update_proposal_sd(e); } +void GGMModel::cholesky_update_after_diag(double omega_ii_old, size_t i) +{ + + double delta = omega_ii_old - omega_prop_(i, i); + + bool s = delta > 0; + vf1_(i) = std::sqrt(std::abs(delta)); + + if (s) + cholesky_downdate(phi_, vf1_); + else + cholesky_update(phi_, vf1_); + + // update inverse (2x O(p^2)) + arma::inv(inv_phi_, arma::trimatu(phi_)); + inv_omega_ = inv_phi_ * inv_phi_.t(); + + // reset for next iteration + vf1_(i) = 0.0; +} + + void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { size_t e = i * (i + 1) / 2 + j; // parameter index in vectorized form @@ -351,27 +414,28 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { edge_indicators_(i, j) = 0; edge_indicators_(j, i) = 0; - // Cholesky update vectors - arma::vec v1 = {0, -1}; - arma::vec v2 = {omega_ij_old - 0.0, (omega_jj_old - omega_(j, j)) / 2}; + cholesky_update_after_edge(omega_ij_old, omega_jj_old, i, j); + // // Cholesky update vectors + // arma::vec v1 = {0, -1}; + // arma::vec v2 = {omega_ij_old - 0.0, (omega_jj_old - omega_(j, j)) / 2}; - arma::vec vf1 = arma::zeros(p_); - arma::vec vf2 = arma::zeros(p_); - vf1[i] = v1[0]; - vf1[j] = v1[1]; - vf2[i] = v2[0]; - vf2[j] = v2[1]; + // arma::vec vf1 = arma::zeros(p_); + // arma::vec vf2 = arma::zeros(p_); + // vf1[i] = v1[0]; + // vf1[j] = v1[1]; + // vf2[i] = v2[0]; + // vf2[j] = v2[1]; - arma::vec u1 = (vf1 + vf2) / sqrt(2); - arma::vec u2 = (vf1 - vf2) / sqrt(2); + // arma::vec u1 = (vf1 + vf2) / sqrt(2); + // arma::vec u2 = (vf1 - vf2) / sqrt(2); - // Update Cholesky factor - cholesky_update(phi_, u1); - cholesky_downdate(phi_, u2); + // // Update Cholesky factor + // cholesky_update(phi_, u1); + // cholesky_downdate(phi_, u2); - // Update inverse - arma::inv(inv_phi_, arma::trimatu(phi_)); - inv_omega_ = inv_phi_ * inv_phi_.t(); + // // Update inverse + // arma::inv(inv_phi_, arma::trimatu(phi_)); + // inv_omega_ = inv_phi_ * inv_phi_.t(); } } else { @@ -425,27 +489,28 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { edge_indicators_(i, j) = 1; edge_indicators_(j, i) = 1; - // Cholesky update vectors - arma::vec v1 = {0, -1}; - arma::vec v2 = {omega_ij_old - omega_(i, j), (omega_jj_old - omega_(j, j)) / 2}; + cholesky_update_after_edge(omega_ij_old, omega_jj_old, i, j); + // // Cholesky update vectors + // arma::vec v1 = {0, -1}; + // arma::vec v2 = {omega_ij_old - omega_(i, j), (omega_jj_old - omega_(j, j)) / 2}; - arma::vec vf1 = arma::zeros(p_); - arma::vec vf2 = arma::zeros(p_); - vf1[i] = v1[0]; - vf1[j] = v1[1]; - vf2[i] = v2[0]; - vf2[j] = v2[1]; + // arma::vec vf1 = arma::zeros(p_); + // arma::vec vf2 = arma::zeros(p_); + // vf1[i] = v1[0]; + // vf1[j] = v1[1]; + // vf2[i] = v2[0]; + // vf2[j] = v2[1]; - arma::vec u1 = (vf1 + vf2) / sqrt(2); - arma::vec u2 = (vf1 - vf2) / sqrt(2); + // arma::vec u1 = (vf1 + vf2) / sqrt(2); + // arma::vec u2 = (vf1 - vf2) / sqrt(2); - // Update Cholesky factor - cholesky_update(phi_, u1); - cholesky_downdate(phi_, u2); + // // Update Cholesky factor + // cholesky_update(phi_, u1); + // cholesky_downdate(phi_, u2); - // Update inverse - arma::inv(inv_phi_, arma::trimatu(phi_)); - inv_omega_ = inv_phi_ * inv_phi_.t(); + // // Update inverse + // arma::inv(inv_phi_, arma::trimatu(phi_)); + // inv_omega_ = inv_phi_ * inv_phi_.t(); } } } diff --git a/src/ggm_model.h b/src/ggm_model.h index 48158588..48af799e 100644 --- a/src/ggm_model.h +++ b/src/ggm_model.h @@ -133,6 +133,13 @@ class GGMModel : public BaseModel { arma::mat omega_prop_; arma::vec constants_; // Phi_q1q, Phi_q1q1, c[1], c[2], c[3], c[4] + arma::vec v1_ = {0, -1}; + arma::vec v2_ = {0, 0}; + arma::vec vf1_ = arma::zeros(p_); + arma::vec vf2_ = arma::zeros(p_); + arma::vec u1_ = arma::zeros(p_); + arma::vec u2_ = arma::zeros(p_); + // Parameter group updates with optimized likelihood evaluations void update_edge_parameter(size_t i, size_t j); void update_diagonal_parameter(size_t i); @@ -147,6 +154,8 @@ class GGMModel : public BaseModel { double log_density_impl_edge(size_t i, size_t j) const; double log_density_impl_diag(size_t j) const; double get_log_det(arma::mat triangular_A) const; + void cholesky_update_after_edge(double omega_ij_old, double omega_jj_old, size_t i, size_t j); + void cholesky_update_after_diag(double omega_ii_old, size_t i); // double find_reasonable_step_size_edge(const arma::mat& omega, size_t i, size_t j); // double find_reasonable_step_size_diag(const arma::mat& omega, size_t i); // double edge_log_ratio(const arma::mat& omega, size_t i, size_t j, double proposal); From 7f74ab0fc8214fff4afef3cba618b300866f5915 Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Tue, 2 Dec 2025 09:36:01 +0100 Subject: [PATCH 5/5] R interface for raw data and sufficient statistics --- R/RcppExports.R | 4 ++-- src/RcppExports.cpp | 8 ++++---- src/ggm_model.cpp | 32 ++++++++++++++++++++++++++++++++ src/ggm_model.h | 33 +++++++++++++++++++++++++++++++++ src/sample_ggm.cpp | 5 +++-- 5 files changed, 74 insertions(+), 8 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 6c7fab37..e67adf97 100644 --- a/R/RcppExports.R +++ b/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) { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 0b7e900f..d0cf3f0f 100644 --- a/src/RcppExports.cpp +++ b/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 } diff --git a/src/ggm_model.cpp b/src/ggm_model.cpp index 9d475176..6b8e0202 100644 --- a/src/ggm_model.cpp +++ b/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(inputFromR["n"]); + arma::mat suf_stat = Rcpp::as(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(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'."); + } + +} diff --git a/src/ggm_model.h b/src/ggm_model.h index 48af799e..dcf4a35e 100644 --- a/src/ggm_model.h +++ b/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(p_, p_)), + phi_(arma::eye(p_, p_)), + inv_phi_(arma::eye(p_, p_)), + inv_omega_(arma::eye(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 +); diff --git a/src/sample_ggm.cpp b/src/sample_ggm.cpp index 8ca13568..4462e464 100644 --- a/src/sample_ggm.cpp +++ b/src/sample_ggm.cpp @@ -159,7 +159,7 @@ Rcpp::List convert_sampler_output_to_ggm_result(const std::vector