diff --git a/src/models/ggm/ggm_model.cpp b/src/models/ggm/ggm_model.cpp index 03fa5029..501d7564 100644 --- a/src/models/ggm/ggm_model.cpp +++ b/src/models/ggm/ggm_model.cpp @@ -482,8 +482,7 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { cholesky_update_after_edge(omega_ij_old, omega_jj_old, i, j); - constraint_dirty_ = true; - theta_valid_ = false; + invalidate_gradient_cache(); } } else { @@ -539,8 +538,7 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { cholesky_update_after_edge(omega_ij_old, omega_jj_old, i, j); - constraint_dirty_ = true; - theta_valid_ = false; + invalidate_gradient_cache(); } } } @@ -641,8 +639,7 @@ void GGMModel::tune_proposal_sd(int iteration, const WarmupSchedule& schedule) { } // Invalidate gradient cache after MH updates - constraint_dirty_ = true; - theta_valid_ = false; + invalidate_gradient_cache(); } void GGMModel::refresh_cholesky() { diff --git a/src/models/ggm/ggm_model.h b/src/models/ggm/ggm_model.h index 2cbd3417..2e96f034 100644 --- a/src/models/ggm/ggm_model.h +++ b/src/models/ggm/ggm_model.h @@ -47,30 +47,22 @@ class GGMModel : public BaseModel { std::unique_ptr interaction_prior, std::unique_ptr diagonal_prior, const bool na_impute = false - ) : n_(observations.n_rows - 1), // centered data has n-1 effective df - p_(observations.n_cols), - dim_((p_ * (p_ + 1)) / 2), - suf_stat_(observations.t() * observations), - inclusion_probability_(inclusion_probability), - edge_selection_(edge_selection), - interaction_prior_(std::move(interaction_prior)), - diagonal_prior_(std::move(diagonal_prior)), - precision_matrix_(arma::eye(p_, p_)), - cholesky_of_precision_(arma::eye(p_, p_)), - inv_cholesky_of_precision_(arma::eye(p_, p_)), - covariance_matrix_(arma::eye(p_, p_)), - edge_indicators_(initial_edge_indicators), - vectorized_parameters_(dim_), - vectorized_indicator_parameters_(edge_selection_ ? dim_ : 0), - proposal_sds_(arma::mat(dim_, 1, arma::fill::ones) * 0.25), - num_pairwise_(p_ * (p_ - 1) / 2), - observations_(na_impute ? observations : arma::mat()), - precision_proposal_(arma::mat(p_, p_, arma::fill::none)) + ) // Delegate to the sufficient-statistics constructor: the raw-data path + // differs only in deriving n (centered data has n-1 effective df) and + // S = X'X from the observations, and in retaining the observations for + // missing-data imputation. All member init and the MLE warm-start live + // in the delegated constructor. + : GGMModel(observations.n_rows - 1, + observations.t() * observations, + inclusion_probability, + initial_edge_indicators, + edge_selection, + std::move(interaction_prior), + std::move(diagonal_prior)) { - int num_edges = arma::accu(edge_indicators_) / 2; - int max_edges = static_cast(p_ * (p_ - 1) / 2); - has_sparse_graph_ = !edge_selection_ && (num_edges < max_edges); - initialize_precision_from_mle(); + if (na_impute) { + observations_ = observations; + } } /** @@ -680,6 +672,21 @@ class GGMModel : public BaseModel { mutable arma::vec theta_; public: + /** + * Mark the NUTS gradient caches stale after a change to the precision + * matrix or the edge set. Both the constraint structure (active dimension + * depends on the edge indicators) and the cached theta parameterization + * become invalid, so the MH within-K, edge-selection, and proposal-SD + * tuning paths all funnel their invalidation through here rather than + * repeating the two flag writes. (set_determinant_tilt is intentionally + * NOT routed through this: it only needs a gradient-engine rebuild, not a + * theta invalidation.) + */ + void invalidate_gradient_cache() { + constraint_dirty_ = true; + theta_valid_ = false; + } + /** * Rebuild the constraint structure and gradient engine from current * edge indicators. Called lazily before gradient evaluation.