refactor(s5): share GGM/Mixed precision kernels as cholesky_helpers free functions

src/math/cholesky_helpers.h

@@ -10,7 +10,9 @@
  */
 
 #include <RcppArmadillo.h>
+#include <array>
 #include <cmath>
+#include "math/explog_macros.h"
 
 namespace cholesky_helpers {
 
@@ -42,4 +44,90 @@ inline double compute_inv_submatrix_i(const arma::mat& A, size_t i,
     return A(ii, jj) - A(ii, i) * A(jj, i) / A(i, i);
 }
 
+/**
+ * Rank-2 matrix-determinant-lemma log-ratio log|K'| − log|K| for a precision
+ * change confined to entries (i,j), (j,i), (j,j). cc11/cc12/cc22 are the 2×2
+ * update-Gram-matrix entries I + Vᵀ Σ U.
+ *
+ * Shared by GGMModel::log_det_ratio_edge and MixedMRFModel::log_det_ratio_yy_edge,
+ * which differ only in where the current precision entries come from (GGM stores
+ * K directly; Mixed stores −½K), so they pass the resolved precision scalars in.
+ *
+ * @param covariance      Σ = K⁻¹.
+ * @param i, j            Off-diagonal/diagonal indices (i < j).
+ * @param precision_curr_ij, precision_curr_jj  Current K entries.
+ * @param precision_prop_ij, precision_prop_jj  Proposed K entries.
+ */
+inline double log_det_ratio_edge_kernel(
+        const arma::mat& covariance, size_t i, size_t j,
+        double precision_curr_ij, double precision_prop_ij,
+        double precision_curr_jj, double precision_prop_jj) {
+    double Ui2 = precision_curr_ij - precision_prop_ij;
+    double Uj2 = (precision_curr_jj - precision_prop_jj) / 2;
+
+    double cc11 = covariance(j, j);
+    double cc12 = 1.0 - (covariance(i, j) * Ui2 + covariance(j, j) * Uj2);
+    double cc22 = Ui2 * Ui2 * covariance(i, i)
+                + 2.0 * Ui2 * Uj2 * covariance(i, j)
+                + Uj2 * Uj2 * covariance(j, j);
+
+    return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12));
+}
+
+/**
+ * Rank-1 specialisation of log_det_ratio_edge_kernel (Ui2 = 0): a diagonal-only
+ * precision change at (j,j). Shared by GGMModel::log_det_ratio_diag and
+ * MixedMRFModel::log_det_ratio_yy_diag.
+ */
+inline double log_det_ratio_diag_kernel(
+        const arma::mat& covariance, size_t j,
+        double precision_curr_jj, double precision_prop_jj) {
+    double Uj2 = (precision_curr_jj - precision_prop_jj) / 2;
+
+    double cc11 = covariance(j, j);
+    double cc12 = 1.0 - covariance(j, j) * Uj2;
+    double cc22 = Uj2 * Uj2 * covariance(j, j);
+
+    return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12));
+}
+
+/**
+ * Extract the six reparameterization constants for a rank-2 off-diagonal MH
+ * proposal in precision space (Roverato move). Shared by GGMModel::get_constants
+ * and MixedMRFModel::get_precision_constants, which differ only in where the
+ * current precision entries come from (K vs −½K), passed in as scalars.
+ *
+ * @param cholesky      Upper-triangular Cholesky factor of K.
+ * @param covariance    Σ = K⁻¹.
+ * @param i, j          Off-diagonal/diagonal indices (i < j).
+ * @param precision_ij, precision_jj  Current K entries K(i,j), K(j,j).
+ * @return  {Phi_q1q, Phi_q1q1, c2, Phi_q1q1, c4, c5}.
+ */
+inline std::array<double, 6> precision_proposal_constants(
+        const arma::mat& cholesky, const arma::mat& covariance,
+        size_t i, size_t j, double precision_ij, double precision_jj) {
+    double logdet = get_log_det(cholesky);
+
+    double log_adj_ii = logdet + MY_LOG(std::abs(covariance(i, i)));
+    double log_adj_ij = logdet + MY_LOG(std::abs(covariance(i, j)));
+    double log_adj_jj = logdet + MY_LOG(std::abs(covariance(j, j)));
+
+    double inv_sub_jj = compute_inv_submatrix_i(covariance, i, j, j);
+    double log_abs_inv_sub_jj = log_adj_ii + MY_LOG(std::abs(inv_sub_jj));
+
+    double Phi_q1q  = (2 * std::signbit(covariance(i, j)) - 1) * MY_EXP(
+        (log_adj_ij - (log_adj_jj + log_abs_inv_sub_jj) / 2)
+    );
+    double Phi_q1q1 = MY_EXP((log_adj_jj - log_abs_inv_sub_jj) / 2);
+
+    std::array<double, 6> c{};
+    c[0] = Phi_q1q;
+    c[1] = Phi_q1q1;
+    c[2] = precision_ij - Phi_q1q * Phi_q1q1;
+    c[3] = Phi_q1q1;
+    c[4] = precision_jj - Phi_q1q * Phi_q1q;
+    c[5] = c[4] + c[2] * c[2] / (c[3] * c[3]);
+    return c;
+}
+
 } // namespace cholesky_helpers

src/models/ggm/ggm_model.cpp

@@ -162,27 +162,10 @@ arma::vec GGMModel::get_active_inv_mass() const {
 
 
 void GGMModel::get_constants(size_t i, size_t j) {
-
-    double logdet_omega = cholesky_helpers::get_log_det(cholesky_of_precision_);
-
-    double log_adj_omega_ii = logdet_omega + MY_LOG(std::abs(covariance_matrix_(i, i)));
-    double log_adj_omega_ij = logdet_omega + MY_LOG(std::abs(covariance_matrix_(i, j)));
-    double log_adj_omega_jj = logdet_omega + MY_LOG(std::abs(covariance_matrix_(j, j)));
-
-    double inv_omega_sub_j1j1 = cholesky_helpers::compute_inv_submatrix_i(covariance_matrix_, i, j, j);
-    double log_abs_inv_omega_sub_jj = log_adj_omega_ii + MY_LOG(std::abs(inv_omega_sub_j1j1));
-    double Phi_q1q  = (2 * std::signbit(covariance_matrix_(i, j)) - 1) * MY_EXP(
-        (log_adj_omega_ij - (log_adj_omega_jj + log_abs_inv_omega_sub_jj) / 2)
-    );
-    double Phi_q1q1 = MY_EXP((log_adj_omega_jj - log_abs_inv_omega_sub_jj) / 2);
-
-    constants_[0] = Phi_q1q;
-    constants_[1] = Phi_q1q1;
-    constants_[2] = precision_matrix_(i, j) - Phi_q1q * Phi_q1q1;
-    constants_[3] = Phi_q1q1;
-    constants_[4] = precision_matrix_(j, j) - Phi_q1q * Phi_q1q;
-    constants_[5] = constants_[4] + constants_[2] * constants_[2] / (constants_[3] * constants_[3]);
-
+    // GGM stores K directly, so the precision entries are precision_matrix_.
+    constants_ = cholesky_helpers::precision_proposal_constants(
+        cholesky_of_precision_, covariance_matrix_, i, j,
+        precision_matrix_(i, j), precision_matrix_(j, j));
 }
 
 double GGMModel::constrained_diagonal(const double x) const {
@@ -205,31 +188,18 @@ double GGMModel::log_density_impl(const arma::mat& omega, const arma::mat& phi)
 
 double GGMModel::log_det_ratio_edge(size_t i, size_t j) const {
     // Rank-2 matrix-determinant lemma: log|K_prop| - log|K_curr| where K_prop
-    // differs from K_curr at entries (i,j), (j,i), and (j,j). cc11, cc12, cc22
-    // are the entries of the 2x2 update Gram matrix I + V^T Sigma U used by
-    // the lemma; |det(...)| is its determinant.
-    double Ui2 = precision_matrix_(i, j) - precision_proposal_(i, j);
-    double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2;
-
-    double cc11 = covariance_matrix_(j, j);
-    double cc12 = 1 - (covariance_matrix_(i, j) * Ui2
-                       + covariance_matrix_(j, j) * Uj2);
-    double cc22 = Ui2 * Ui2 * covariance_matrix_(i, i)
-                + 2 * Ui2 * Uj2 * covariance_matrix_(i, j)
-                + Uj2 * Uj2 * covariance_matrix_(j, j);
-
-    return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12));
+    // differs from K_curr at entries (i,j), (j,i), and (j,j). GGM stores K
+    // directly, so the precision scalars are precision_matrix_ entries.
+    return cholesky_helpers::log_det_ratio_edge_kernel(
+        covariance_matrix_, i, j,
+        precision_matrix_(i, j), precision_proposal_(i, j),
+        precision_matrix_(j, j), precision_proposal_(j, j));
 }
 
 double GGMModel::log_det_ratio_diag(size_t j) const {
-    // Rank-1 specialisation of log_det_ratio_edge (Ui2 = 0).
-    double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2;
-
-    double cc11 = covariance_matrix_(j, j);
-    double cc12 = 1 - covariance_matrix_(j, j) * Uj2;
-    double cc22 = Uj2 * Uj2 * covariance_matrix_(j, j);
-
-    return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12));
+    return cholesky_helpers::log_det_ratio_diag_kernel(
+        covariance_matrix_, j,
+        precision_matrix_(j, j), precision_proposal_(j, j));
 }
 
 double GGMModel::log_density_impl_edge(size_t i, size_t j) const {

src/models/mixed/mixed_mrf_metropolis.cpp

@@ -182,31 +182,14 @@ double MixedMRFModel::update_pairwise_discrete(int i, int j, std::optional<doubl
 // =============================================================================
 
 void MixedMRFModel::get_precision_constants(int i, int j) {
-    double logdet = cholesky_helpers::get_log_det(cholesky_of_precision_);
-
-    double log_adj_ii = logdet + MY_LOG(std::abs(covariance_continuous_(i, i)));
-    double log_adj_ij = logdet + MY_LOG(std::abs(covariance_continuous_(i, j)));
-    double log_adj_jj = logdet + MY_LOG(std::abs(covariance_continuous_(j, j)));
-
-    double inv_sub_jj = cholesky_helpers::compute_inv_submatrix_i(covariance_continuous_, i, j, j);
-    double log_abs_inv_sub_jj = log_adj_ii + MY_LOG(std::abs(inv_sub_jj));
-
-    double Phi_q1q  = (2 * std::signbit(covariance_continuous_(i, j)) - 1) * MY_EXP(
-        (log_adj_ij - (log_adj_jj + log_abs_inv_sub_jj) / 2)
-    );
-    double Phi_q1q1 = MY_EXP((log_adj_jj - log_abs_inv_sub_jj) / 2);
-
-    // Read precision = -2 * pairwise_effects_continuous_ for constants extraction
-    double precision_ij = -2.0 * pairwise_effects_continuous_(i, j);
-    double precision_jj = -2.0 * pairwise_effects_continuous_(j, j);
-
-    cont_constants_[0] = Phi_q1q;
-    cont_constants_[1] = Phi_q1q1;
-    cont_constants_[2] = precision_ij - Phi_q1q * Phi_q1q1;
-    cont_constants_[3] = Phi_q1q1;
-    cont_constants_[4] = precision_jj - Phi_q1q * Phi_q1q;
-    cont_constants_[5] = cont_constants_[4] +
-        cont_constants_[2] * cont_constants_[2] / (cont_constants_[3] * cont_constants_[3]);
+    // Kyy = -2 * pairwise_effects_continuous_; delegate to the shared kernel,
+    // passing the resolved precision entries.
+    size_t ui = static_cast<size_t>(i);
+    size_t uj = static_cast<size_t>(j);
+    cont_constants_ = cholesky_helpers::precision_proposal_constants(
+        cholesky_of_precision_, covariance_continuous_, ui, uj,
+        -2.0 * pairwise_effects_continuous_(ui, uj),
+        -2.0 * pairwise_effects_continuous_(uj, uj));
 }
 
 double MixedMRFModel::precision_constrained_diagonal(double x) const {
@@ -231,37 +214,20 @@ double MixedMRFModel::precision_constrained_diagonal(double x) const {
 // =============================================================================
 
 double MixedMRFModel::log_det_ratio_yy_edge(int i, int j) const {
-    // Rank-2 matrix-determinant lemma applied to the Kyy block. Mirrors
-    // GGMModel::log_det_ratio_edge but uses pairwise_effects_continuous_
-    // (Kyy = -2 * pairwise_effects_continuous_) and covariance_continuous_.
+    // Kyy = -2 * pairwise_effects_continuous_; delegate to the shared kernel.
     size_t ui = static_cast<size_t>(i);
     size_t uj = static_cast<size_t>(j);
-    double precision_ij_curr = -2.0 * pairwise_effects_continuous_(ui, uj);
-    double precision_jj_curr = -2.0 * pairwise_effects_continuous_(uj, uj);
-    double Ui = precision_ij_curr - precision_proposal_(ui, uj);
-    double Uj = (precision_jj_curr - precision_proposal_(uj, uj)) / 2.0;
-
-    double cc11 = covariance_continuous_(uj, uj);
-    double cc12 = 1.0 - (covariance_continuous_(ui, uj) * Ui +
-                         covariance_continuous_(uj, uj) * Uj);
-    double cc22 = Ui * Ui * covariance_continuous_(ui, ui) +
-                  2.0 * Ui * Uj * covariance_continuous_(ui, uj) +
-                  Uj * Uj * covariance_continuous_(uj, uj);
-
-    return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12));
+    return cholesky_helpers::log_det_ratio_edge_kernel(
+        covariance_continuous_, ui, uj,
+        -2.0 * pairwise_effects_continuous_(ui, uj), precision_proposal_(ui, uj),
+        -2.0 * pairwise_effects_continuous_(uj, uj), precision_proposal_(uj, uj));
 }
 
 double MixedMRFModel::log_det_ratio_yy_diag(int i) const {
-    // Rank-1 specialisation of log_det_ratio_yy_edge (Ui = 0).
     size_t ui = static_cast<size_t>(i);
-    double precision_ii_curr = -2.0 * pairwise_effects_continuous_(ui, ui);
-    double Uj = (precision_ii_curr - precision_proposal_(ui, ui)) / 2.0;
-
-    double cc11 = covariance_continuous_(ui, ui);
-    double cc12 = 1.0 - covariance_continuous_(ui, ui) * Uj;
-    double cc22 = Uj * Uj * covariance_continuous_(ui, ui);
-
-    return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12));
+    return cholesky_helpers::log_det_ratio_diag_kernel(
+        covariance_continuous_, ui,
+        -2.0 * pairwise_effects_continuous_(ui, ui), precision_proposal_(ui, ui));
 }