bgmCompare: persistent sweep state and current-state normalizer cache

src/mcmc/algorithms/metropolis.cpp

@@ -23,3 +23,23 @@ StepResult metropolis_step(
 
   return {State, accept_prob};
 }
+
+
+StepResult metropolis_step_cached(
+    double current_state,
+    double step_size,
+    double log_post_current,
+    const std::function<double(double)>& log_post_proposed,
+    SafeRNG& rng
+) {
+  double proposed_state = rnorm(rng, current_state, step_size);
+  double log_accept = log_post_proposed(proposed_state) - log_post_current;
+  double accept_prob = std::min(1.0, MY_EXP(log_accept));
+
+  double state = (runif(rng) < accept_prob) ? proposed_state : current_state;
+
+  arma::vec State(1);
+  State[0] = state;
+
+  return {State, accept_prob};
+}

src/mcmc/algorithms/metropolis.h

@@ -24,3 +24,26 @@ StepResult metropolis_step(
     const std::function<double(double)>& log_post,
     SafeRNG& rng
 );
+
+/**
+ * Performs one Random Walk Metropolis step with a cached current-state value
+ *
+ * Variant of metropolis_step() for callers that maintain the current-state
+ * log-posterior across proposals: only the proposed state is evaluated. The
+ * random-number stream matches metropolis_step() (one normal draw, one
+ * uniform draw).
+ *
+ * @param current_state     Current scalar parameter value
+ * @param step_size         Standard deviation of the Gaussian proposal
+ * @param log_post_current  Log-posterior value at current_state
+ * @param log_post_proposed Log-posterior function evaluated at the proposal
+ * @param rng               Thread-safe random number generator
+ * @return StepResult with accepted state (1-element vector) and acceptance probability
+ */
+StepResult metropolis_step_cached(
+    double current_state,
+    double step_size,
+    double log_post_current,
+    const std::function<double(double)>& log_post_proposed,
+    SafeRNG& rng
+);

src/models/bgmCompare/bgmCompare_helper.cpp

@@ -1,10 +1,74 @@
 #include <RcppArmadillo.h>
 #include <cmath>
 #include "models/bgmCompare/bgmCompare_helper.h"
+#include "models/bgmCompare/bgmCompare_state.h"
 #include "utils/common_helpers.h"
 
 
 
+// Converts the integer observations to double, whole-matrix and per group.
+// Weight and residual members are left untouched; call
+// rebuild_sweep_state_weights() afterwards to make the state consistent.
+void initialize_sweep_state_observations(
+    CompareSweepState& state,
+    const arma::imat& observations,
+    const arma::imat& group_indices,
+    const int num_groups
+) {
+  state.obs_double_all = arma::conv_to<arma::mat>::from(observations);
+  state.obs_double.resize(num_groups);
+  for (int g = 0; g < num_groups; g++) {
+    const int r0 = group_indices(g, 0);
+    const int r1 = group_indices(g, 1);
+    state.obs_double[g] = state.obs_double_all.rows(r0, r1);
+  }
+
+  const int num_variables = observations.n_cols;
+  state.log_normalizer.zeros(num_variables, num_groups);
+  state.normalizer_valid.zeros(num_variables);
+}
+
+
+
+// Recomputes the per-group effective pairwise weights from the current
+// pairwise effects and inclusion indicators, and the residual matrices as
+// one matrix product per group. Invalidates the normalizer cache.
+void rebuild_sweep_state_weights(
+    CompareSweepState& state,
+    const arma::mat& pairwise_effects,
+    const arma::imat& pairwise_effect_indices,
+    const arma::imat& inclusion_indicator,
+    const arma::mat& projection,
+    const int num_groups
+) {
+  const int num_variables = inclusion_indicator.n_rows;
+  state.pairwise_group.resize(num_groups);
+  state.residual.resize(num_groups);
+
+  for (int g = 0; g < num_groups; g++) {
+    const arma::vec proj_g = projection.row(g).t();
+
+    arma::mat& pairwise_g = state.pairwise_group[g];
+    pairwise_g.zeros(num_variables, num_variables);
+    for (int v = 0; v < num_variables - 1; v++) {
+      for (int u = v + 1; u < num_variables; u++) {
+        double w = compute_group_pairwise_effects(
+          v, u, num_groups, pairwise_effects, pairwise_effect_indices,
+          inclusion_indicator, proj_g
+        );
+        pairwise_g(v, u) = w;
+        pairwise_g(u, v) = w;
+      }
+    }
+
+    state.residual[g] = state.obs_double[g] * pairwise_g;
+  }
+
+  state.normalizer_valid.zeros();
+}
+
+
+
 // Computes group-specific main effects for a given variable (bgmCompare model).
 //
 // For a variable, the main-effect parameters are stored with: