From 4591a44a7085d18b3486f5b698d795ba5e188e51 Mon Sep 17 00:00:00 2001 From: brunoguindani Date: Wed, 18 May 2022 10:18:17 +0200 Subject: [PATCH] Comments cleanup and README updates --- CMakeLists.txt | 7 -- README.md | 9 +- docs/conf.py | 2 - src/hierarchies/fa_hierarchy.h | 3 +- .../likelihoods/multi_norm_likelihood.cc | 4 +- .../likelihoods/states/CMakeLists.txt | 2 +- src/hierarchies/nnw_hierarchy.h | 1 + src/hierarchies/priors/fa_prior_model.cc | 2 - .../updaters/semi_conjugate_updater.h | 2 +- src/utils/covariates_getter.h | 2 - test/hierarchies.cc | 38 --------- test/likelihoods.cc | 2 +- test/lpdf.cc | 84 ------------------- test/write_proto.cc | 24 ------ 14 files changed, 10 insertions(+), 172 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a4d7cea61..e79f07f43 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -171,7 +171,6 @@ if (NOT DISABLE_DOCS) add_subdirectory(docs) endif() -# add_subdirectory(examples) if (NOT DISABLE_PLOTS) include(FetchContent) @@ -210,9 +209,3 @@ endif() if (NOT DISABLE_EXAMPLES) add_subdirectory(examples) endif() - -# Test MH updater -# add_executable(test_mh $ test_mh_updater.cpp) -# target_include_directories(test_mh PUBLIC ${INCLUDE_PATHS}) -# target_link_libraries(test_mh PUBLIC ${LINK_LIBRARIES}) -# target_compile_options(test_mh PUBLIC ${COMPILE_OPTIONS}) diff --git a/README.md b/README.md index 019625c81..7cb6972a4 100644 --- a/README.md +++ b/README.md @@ -16,13 +16,7 @@ Current state of the software: --> -where P is either the Dirichlet process or the Pitman--Yor process - -- We currently support univariate and multivariate location-scale mixture of Gaussian densities - -- Inference is carried out using algorithms such as Algorithm 2 in [Neal (2000)](http://www.stat.columbia.edu/npbayes/papers/neal_sampling.pdf) - -- Serialization of the MCMC chains is possible using Google's [Protocol Buffers](https://developers.google.com/protocol-buffers) aka `protobuf` +For descriptions of the models supported in our library, discussion of software design, and examples, please refer to the following paper: https://arxiv.org/abs/2205.08144 # Installation @@ -101,3 +95,4 @@ Documentation is available at https://bayesmix.readthedocs.io. # Contributions are welcome! Please check out [CONTRIBUTING.md](CONTRIBUTING.md) for details on how to collaborate with us. +You can also head to our [issues page](https://github.com/bayesmix-dev/bayesmix/issues) to check for useful enhancements needed. diff --git a/docs/conf.py b/docs/conf.py index af16e0bdc..9acbf232a 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -57,8 +57,6 @@ def configureDoxyfile(input_dir, output_dir): html_theme = 'haiku' -# html_static_path = ['_static'] - highlight_language = 'cpp' imgmath_latex = 'latex' diff --git a/src/hierarchies/fa_hierarchy.h b/src/hierarchies/fa_hierarchy.h index 17da43564..fc7ce3079 100644 --- a/src/hierarchies/fa_hierarchy.h +++ b/src/hierarchies/fa_hierarchy.h @@ -62,10 +62,11 @@ class FAHierarchy } }; -//! Empirical-Bayes hyperparameters initialization for the FA HIerarchy. +//! Empirical-Bayes hyperparameters initialization for the FAHierarchy. //! Sets the hyperparameters in `hier` starting from the data on which the user //! wants to fit the model. inline void set_fa_hyperparams_from_data(FAHierarchy* hier) { + // TODO test this function auto dataset_ptr = std::static_pointer_cast(hier->get_likelihood()) ->get_dataset(); diff --git a/src/hierarchies/likelihoods/multi_norm_likelihood.cc b/src/hierarchies/likelihoods/multi_norm_likelihood.cc index ae51f3d47..f0cfae90d 100644 --- a/src/hierarchies/likelihoods/multi_norm_likelihood.cc +++ b/src/hierarchies/likelihoods/multi_norm_likelihood.cc @@ -12,8 +12,8 @@ double MultiNormLikelihood::compute_lpdf( void MultiNormLikelihood::update_sum_stats(const Eigen::RowVectorXd &datum, bool add) { - // Check if dim is not defined yet (usually not happens if hierarchy is - // initialized) + // Check if dim is not defined yet (this usually doesn't happen if the + // hierarchy is initialized) if (!dim) set_dim(datum.size()); // Updates if (add) { diff --git a/src/hierarchies/likelihoods/states/CMakeLists.txt b/src/hierarchies/likelihoods/states/CMakeLists.txt index 4f850a603..933c337b4 100644 --- a/src/hierarchies/likelihoods/states/CMakeLists.txt +++ b/src/hierarchies/likelihoods/states/CMakeLists.txt @@ -4,5 +4,5 @@ target_sources(bayesmix PUBLIC uni_ls_state.h multi_ls_state.h uni_lin_reg_ls_state.h - # fa_state.h + fa_state.h ) diff --git a/src/hierarchies/nnw_hierarchy.h b/src/hierarchies/nnw_hierarchy.h index 7d8bd330d..2cf36464a 100644 --- a/src/hierarchies/nnw_hierarchy.h +++ b/src/hierarchies/nnw_hierarchy.h @@ -65,6 +65,7 @@ class NNWHierarchy //! @return The evaluation of the lpdf double marg_lpdf(ProtoHypersPtr hier_params, const Eigen::RowVectorXd &datum) const override { + // TODO check Bayes rule for this hierarchy HyperParams pred_params = get_predictive_t_parameters(hier_params); Eigen::VectorXd diag = pred_params.scale_chol.diagonal(); double logdet = 2 * log(diag.array()).sum(); diff --git a/src/hierarchies/priors/fa_prior_model.cc b/src/hierarchies/priors/fa_prior_model.cc index d75442795..fa402b1d1 100644 --- a/src/hierarchies/priors/fa_prior_model.cc +++ b/src/hierarchies/priors/fa_prior_model.cc @@ -7,8 +7,6 @@ double FAPriorModel::lpdf(const google::protobuf::Message &state_) { // Proto2Eigen conversion Eigen::VectorXd mu = bayesmix::to_eigen(state.mu()); Eigen::VectorXd psi = bayesmix::to_eigen(state.psi()); - - // Eigen::MatrixXd eta = bayesmix::to_eigen(state.eta()); Eigen::MatrixXd lambda = bayesmix::to_eigen(state.lambda()); // Initialize lpdf value diff --git a/src/hierarchies/updaters/semi_conjugate_updater.h b/src/hierarchies/updaters/semi_conjugate_updater.h index 41e3bc23a..5609bf1b8 100644 --- a/src/hierarchies/updaters/semi_conjugate_updater.h +++ b/src/hierarchies/updaters/semi_conjugate_updater.h @@ -66,7 +66,7 @@ void SemiConjugateUpdater::draw( auto& likecast = downcast_likelihood(like); auto& priorcast = downcast_prior(prior); // Sample from the full conditional of a semi-conjugate hierarchy - bool set_card = true; /*, use_post_hypers=true;*/ + bool set_card = true; if (likecast.get_card() == 0) { likecast.set_state(priorcast.sample(), !set_card); } else { diff --git a/src/utils/covariates_getter.h b/src/utils/covariates_getter.h index 46adb9802..530590182 100644 --- a/src/utils/covariates_getter.h +++ b/src/utils/covariates_getter.h @@ -2,8 +2,6 @@ #define BAYESMIX_SRC_UTILS_COVARIATES_GETTER_H #include -// #include "src/hierarchies/likelihoods/abstract_likelihood.h" -// #include "src/hierarchies/priors/abstract_prior_model.h" class covariates_getter { protected: diff --git a/test/hierarchies.cc b/test/hierarchies.cc index db3090f0c..0a4a94c25 100644 --- a/test/hierarchies.cc +++ b/test/hierarchies.cc @@ -335,44 +335,6 @@ TEST(fa_hierarchy, draw) { ASSERT_TRUE(clusval->DebugString() != clusval2->DebugString()); } -// TEST(fahierarchy, draw_auto) { -// auto hier = std::make_shared(); -// bayesmix::FAPrior prior; -// Eigen::VectorXd mutilde(0); -// bayesmix::Vector mutilde_proto; -// bayesmix::to_proto(mutilde, &mutilde_proto); -// int q = 2; -// double phi = 1.0; -// double alpha0 = 5.0; -// Eigen::VectorXd beta(0); -// bayesmix::Vector beta_proto; -// bayesmix::to_proto(beta, &beta_proto); -// Eigen::MatrixXd dataset(5, 5); -// dataset << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, -// 19, -// 20, 1, 5, 7, 8, 9; -// hier->set_dataset(&dataset); -// *prior.mutable_fixed_values()->mutable_mutilde() = mutilde_proto; -// prior.mutable_fixed_values()->set_phi(phi); -// prior.mutable_fixed_values()->set_alpha0(alpha0); -// prior.mutable_fixed_values()->set_q(q); -// *prior.mutable_fixed_values()->mutable_beta() = beta_proto; -// hier->get_mutable_prior()->CopyFrom(prior); -// hier->initialize(); - -// auto hier2 = hier->clone(); -// hier2->sample_prior(); - -// bayesmix::AlgorithmState out; -// bayesmix::AlgorithmState::ClusterState* clusval = -// out.add_cluster_states(); bayesmix::AlgorithmState::ClusterState* clusval2 -// = out.add_cluster_states(); hier->write_state_to_proto(clusval); -// hier2->write_state_to_proto(clusval2); - -// ASSERT_TRUE(clusval->DebugString() != clusval2->DebugString()) -// << clusval->DebugString() << clusval2->DebugString(); -// } - TEST(fa_hierarchy, sample_given_data) { auto hier = std::make_shared(); bayesmix::FAPrior prior; diff --git a/test/likelihoods.cc b/test/likelihoods.cc index d0fccc775..4cb2a7d5e 100644 --- a/test/likelihoods.cc +++ b/test/likelihoods.cc @@ -380,7 +380,7 @@ TEST(laplace_likelihood, eval_lpdf_unconstrained) { lpdf += like->lpdf(data.row(i)); } - like->set_dataset(&data); // Questa cosa รจ sempre garantita?? + like->set_dataset(&data); double clus_lpdf = like->cluster_lpdf_from_unconstrained(unconstrained_params); ASSERT_DOUBLE_EQ(lpdf, clus_lpdf); diff --git a/test/lpdf.cc b/test/lpdf.cc index 0f41f69ba..fe1cde610 100644 --- a/test/lpdf.cc +++ b/test/lpdf.cc @@ -7,7 +7,6 @@ #include "algorithm_state.pb.h" #include "src/hierarchies/lin_reg_uni_hierarchy.h" #include "src/hierarchies/nnig_hierarchy.h" -// #include "src/hierarchies/nnw_hierarchy.h" #include "src/utils/proto_utils.h" TEST(lpdf, nnig) { @@ -54,89 +53,6 @@ TEST(lpdf, nnig) { ASSERT_DOUBLE_EQ(sum, marg); } -// TEST(lpdf, nnw) { // TODO -// using namespace stan::math; -// NNWHierarchy hier; -// bayesmix::NNWPrior hier_prior; -// Eigen::Vector2d mu0; mu0 << 5.5, 5.5; -// bayesmix::Vector mu0_proto; -// bayesmix::to_proto(mu0, &mu0_proto); -// double lambda0 = 0.2; -// double nu0 = 5.0; -// Eigen::Matrix2d tau0 = Eigen::Matrix2d::Identity() / nu0; -// bayesmix::Matrix tau0_proto; -// bayesmix::to_proto(tau0, &tau0_proto); -// *hier_prior.mutable_fixed_values()->mutable_mean() = mu0_proto; -// hier_prior.mutable_fixed_values()->set_var_scaling(lambda0); -// hier_prior.mutable_fixed_values()->set_deg_free(nu0); -// *hier_prior.mutable_fixed_values()->mutable_scale() = tau0_proto; -// hier.set_prior(hier_prior); -// hier.initialize(); -// -// Eigen::VectorXd mu = mu0; -// Eigen::MatrixXd tau = lambda0 * Eigen::Matrix2d::Identity(); -// -// Eigen::RowVectorXd datum(2); -// datum << 4.5, 4.5; -// -// // Compute prior parameters -// Eigen::MatrixXd tau_pr = lambda0 * tau0; -// -// // Compute posterior parameters -// double mu_n = (lambda0 * mu0 + datum(0)) / (lambda0 + 1); -// double alpha_n = alpha0 + 0.5; -// double lambda_n = lambda0 + 1; -// double nu_n = nu0 + 0.5; -// Eigen::VectorXd mu_n = -// (lambda0 * mu0 + datum.transpose()) / (lambda0 + 1); -// Eigen::MatrixXd tau_temp = -// stan::math::inverse_spd(tau0) + (0.5 * lambda0 / (lambda0 + 1)) * -// (datum.transpose() - mu0) * -// (datum - mu0.transpose()); -// Eigen::MatrixXd tau_n = stan::math::inverse_spd(tau_temp); -// Eigen::MatrixXd tau_post = lambda_n * tau_n; -// -// // Compute pieces -// double prior1 = stan::math::wishart_lpdf(tau, nu0, tau0); -// double prior2 = stan::math::multi_normal_prec_lpdf(mu, mu0, tau_pr); -// double prior = prior1 + prior2; -// double like = hier.get_like_lpdf(datum); -// double post1 = stan::math::wishart_lpdf(tau, nu_n, tau_post); -// double post2 = stan::math::multi_normal_prec_lpdf(mu, mu0, tau_post); -// double post = post1 + post2; -// -// // Bayes: logmarg(x) = logprior(phi) + loglik(x|phi) - logpost(phi|x) -// double sum = prior + like - post; -// double marg = hier.prior_pred_lpdf(false, datum); -// -// // Compute logdet's -// Eigen::MatrixXd tauchol0 = -// Eigen::LLT(tau0).matrixL().transpose(); -// double logdet0 = 2 * log(tauchol0.diagonal().array()).sum(); -// Eigen::MatrixXd tauchol_n = -// Eigen::LLT(tau_n).matrixL().transpose(); -// double logdet_n = 2 * log(tauchol_n.diagonal().array()).sum(); -// -// // lmgamma(dim, x) -// int dim = 2; -// double marg_murphy = lmgamma(dim, 0.5 * nu_n) + 0.5 * nu_n * logdet_n + -// 0.5 * dim * log(lambda0) + dim * NEG_LOG_SQRT_TWO_PI -// - lmgamma(dim, 0.5 * nu0) - 0.5 * nu0 * logdet0 - 0.5 -// * dim * log(lambda_n); -// -// // std::cout << "prior1=" << prior1 << std::endl; -// // std::cout << "prior2=" << prior2 << std::endl; -// // std::cout << "prior =" << prior << std::endl; -// // std::cout << "like =" << like << std::endl; -// // std::cout << "post1 =" << post1 << std::endl; -// // std::cout << "post2 =" << post2 << std::endl; -// // std::cout << "post =" << post << std::endl; -// std::cout << "sum =" << sum << std::endl; -// std::cout << "marg =" << marg << std::endl; -// std::cout << "murphy=" << marg_murphy << std::endl; -// ASSERT_DOUBLE_EQ(marg, marg_murphy); -// } - TEST(lpdf, lin_reg_uni) { // Create hierarchy objects LinRegUniHierarchy hier; diff --git a/test/write_proto.cc b/test/write_proto.cc index 3be320ec9..da866e40a 100644 --- a/test/write_proto.cc +++ b/test/write_proto.cc @@ -3,7 +3,6 @@ #include "algorithm_state.pb.h" #include "ls_state.pb.h" #include "src/hierarchies/nnig_hierarchy.h" -// #include "src/hierarchies/nnw_hierarchy.h" #include "src/utils/proto_utils.h" TEST(set_state, uni_ls) { @@ -46,26 +45,3 @@ TEST(write_proto, uni_ls) { ASSERT_EQ(mean, out_mean); ASSERT_EQ(var, out_var); } - -// TEST(set_state, multi_ls) { -// Eigen::VectorXd mean = Eigen::VectorXd::Ones(5); -// Eigen::MatrixXd prec = Eigen::MatrixXd::Identity(5, 5); -// prec(1, 1) = 10.0; - -// bayesmix::MultiLSState curr; -// bayesmix::to_proto(mean, curr.mutable_mean()); -// bayesmix::to_proto(prec, curr.mutable_prec()); - -// ASSERT_EQ(curr.mean().data(0), 1.0); -// ASSERT_EQ(curr.prec().data(0), 1.0); -// ASSERT_EQ(curr.prec().data(6), 10.0); - -// bayesmix::AlgorithmState::ClusterState clusval_in; -// clusval_in.mutable_multi_ls_state()->CopyFrom(curr); -// NNWHierarchy cluster; -// cluster.set_state_from_proto(clusval_in); - -// ASSERT_EQ(curr.mean().data(0), cluster.get_state().mean(0)); -// ASSERT_EQ(curr.prec().data(0), cluster.get_state().prec(0, 0)); -// ASSERT_EQ(curr.prec().data(6), cluster.get_state().prec(1, 1)); -// }