From 800f1189e34568f09fb2818f2fe785ec2a432bc6 Mon Sep 17 00:00:00 2001 From: vandenman Date: Tue, 27 May 2025 13:54:38 +0200 Subject: [PATCH 01/10] custom exp and log --- .Rbuildignore | 2 + src/Makevars.win | 4 + src/gibbs_functions.cpp | 240 ++++++++++++++++++++++++++++++---------- 3 files changed, 189 insertions(+), 57 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index d5373344..7c2ceb5f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,3 +1,5 @@ +^renv$ +^renv\.lock$ ^.*\.Rproj$ ^\.Rproj\.user$ diff --git a/src/Makevars.win b/src/Makevars.win index 22ebc632..f45d7333 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1 +1,5 @@ +DEBUG=1 PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +PKG_CXXFLAGS+= -g +PKG_DLLFLAGS= +PKG_SHLIB_LDFLAGS= \ No newline at end of file diff --git a/src/gibbs_functions.cpp b/src/gibbs_functions.cpp index a5d8860b..823c71e9 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -5,6 +5,21 @@ #include using namespace Rcpp; +inline double fast_exp(double x) { + // 5th-order Pade approximant + return 1.0 + x * (1.0 + x * (0.5 + x * (1.0/6.0 + x * (1.0/24.0 + x * (1.0/120.0))))); +} + +inline double fast_log_bitwise(double x) { + int* xp = reinterpret_cast(&x); + int log2 = ((xp[1] >> 20) & 0x7FF) - 1023; // Extract exponent bits + double y = x / std::ldexp(1.0, log2); // Normalize mantissa + return log2 * 0.69314718 + (y - 1) - (y - 1)*(y - 1)/2; // linear/quadratic correction +} + +#define MY_EXP fast_exp +#define MY_LOG fast_log_bitwise + // ----------------------------------------------------------------------------| // Impute missing data from full-conditional // ----------------------------------------------------------------------------| @@ -50,7 +65,7 @@ List impute_missing_data( for(int category = 0; category < no_categories[variable]; category++) { exponent = thresholds(variable, category); exponent += (category + 1) * rest_score; - cumsum += std::exp(exponent); + cumsum += MY_EXP(exponent); probabilities[category + 1] = cumsum; } @@ -60,7 +75,7 @@ List impute_missing_data( exponent = thresholds(variable, 1) * reference_category[variable] * reference_category[variable]; - cumsum = std::exp(exponent); + cumsum = MY_EXP(exponent); probabilities[0] = cumsum; for(int category = 0; category < no_categories[variable]; category++) { exponent = thresholds(variable, 0) * (category + 1); @@ -68,7 +83,7 @@ List impute_missing_data( (category + 1 - reference_category[variable]) * (category + 1 - reference_category[variable]); exponent += (category + 1) * rest_score; - cumsum += std::exp(exponent); + cumsum += MY_EXP(exponent); probabilities[category + 1] = cumsum; } } @@ -145,7 +160,7 @@ void metropolis_thresholds_regular( for(int category = 0; category < no_categories[variable]; category++) { current_state = thresholds(variable, category); - exp_current = std::exp(current_state); + exp_current = MY_EXP(current_state); c = (threshold_alpha + threshold_beta) / (1 + exp_current); for(int person = 0; person < no_persons; person++) { g[person] = 1.0; @@ -153,11 +168,11 @@ void metropolis_thresholds_regular( rest_score = rest_matrix(person, variable); for(int cat = 0; cat < no_categories[variable]; cat++) { if(cat != category) { - g[person] += std::exp(thresholds(variable, cat) + + g[person] += MY_EXP(thresholds(variable, cat) + (cat + 1) * rest_score); } } - q[person] = std::exp((category + 1) * rest_score); + q[person] = MY_EXP((category + 1) * rest_score); c += q[person] / (g[person] + q[person] * exp_current); } c = c / ((no_persons + threshold_alpha + threshold_beta) - @@ -167,26 +182,26 @@ void metropolis_thresholds_regular( a = n_cat_obs(category + 1, variable) + threshold_alpha; b = no_persons + threshold_beta - n_cat_obs(category + 1, variable); tmp = R::rbeta(a, b); - proposed_state = std::log(tmp / (1 - tmp) / c); + proposed_state = MY_LOG(tmp / (1 - tmp) / c); exp_proposed = exp(proposed_state); //Compute log_acceptance probability for Metropolis. //First, we use g and q above to compute the ratio of pseudolikelihoods log_prob = 0; for(int person = 0; person < no_persons; person++) { - log_prob += std::log(g[person] + q[person] * exp_current); - log_prob -= std::log(g[person] + q[person] * exp_proposed); + log_prob += MY_LOG(g[person] + q[person] * exp_current); + log_prob -= MY_LOG(g[person] + q[person] * exp_proposed); } //Second, we add the ratio of prior probabilities log_prob -= (threshold_alpha + threshold_beta) * - std::log(1 + exp_proposed); + MY_LOG(1 + exp_proposed); log_prob += (threshold_alpha + threshold_beta) * - std::log(1 + exp_current); + MY_LOG(1 + exp_current); //Third, we add the ratio of proposals - log_prob -= (a + b) * std::log(1 + c * exp_current); - log_prob += (a + b) * std::log(1 + c * exp_proposed); + log_prob -= (a + b) * MY_LOG(1 + c * exp_current); + log_prob += (a + b) * MY_LOG(1 + c * exp_proposed); - U = std::log(R::unif_rand()); + U = MY_LOG(R::unif_rand()); if(U < log_prob) { thresholds(variable, category) = proposed_state; } @@ -262,24 +277,24 @@ void metropolis_thresholds_blumecapel( } else { bound = lbound; } - numerator = std::exp(constant_numerator[0] - bound); - denominator = std::exp(constant_denominator[0] - bound); + numerator = MY_EXP(constant_numerator[0] - bound); + denominator = MY_EXP(constant_denominator[0] - bound); for(int category = 0; category < no_categories[variable]; category ++) { exponent = (category + 1) * rest_score - bound; - numerator += std::exp(constant_numerator[category + 1] + exponent); - denominator += std::exp(constant_denominator[category + 1] + exponent); + numerator += MY_EXP(constant_numerator[category + 1] + exponent); + denominator += MY_EXP(constant_denominator[category + 1] + exponent); } - log_prob += std::log(numerator); - log_prob -= std::log(denominator); + log_prob += MY_LOG(numerator); + log_prob -= MY_LOG(denominator); } log_prob += (threshold_alpha + threshold_beta) * - std::log(1 + std::exp(current_state)); + MY_LOG(1 + MY_EXP(current_state)); log_prob -= (threshold_alpha + threshold_beta) * - std::log(1 + std::exp(proposed_state)); + MY_LOG(1 + MY_EXP(proposed_state)); U = R::unif_rand(); - if(std::log(U) < log_prob) { + if(MY_LOG(U) < log_prob) { thresholds(variable, 0) = proposed_state; } @@ -287,11 +302,11 @@ void metropolis_thresholds_blumecapel( if(log_prob > 0) { log_prob = 1; } else { - log_prob = std::exp(log_prob); + log_prob = MY_EXP(log_prob); } double update_proposal_sd = proposal_sd_blumecapel(variable, 0) + - (log_prob - target_ar) * std::exp(-log(t) * phi); + (log_prob - target_ar) * MY_EXP(-log(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -342,25 +357,25 @@ void metropolis_thresholds_blumecapel( bound = lbound; } - numerator = std::exp(constant_numerator[0] - bound); - denominator = std::exp(constant_denominator[0] - bound); + numerator = MY_EXP(constant_numerator[0] - bound); + denominator = MY_EXP(constant_denominator[0] - bound); for(int category = 0; category < no_categories[variable]; category ++) { exponent = (category + 1) * rest_score - bound; - numerator += std::exp(constant_numerator[category + 1] + exponent); - denominator += std::exp(constant_denominator[category + 1] + exponent); + numerator += MY_EXP(constant_numerator[category + 1] + exponent); + denominator += MY_EXP(constant_denominator[category + 1] + exponent); } - log_prob += std::log(numerator); - log_prob -= std::log(denominator); + log_prob += MY_LOG(numerator); + log_prob -= MY_LOG(denominator); } log_prob += (threshold_alpha + threshold_beta) * - std::log(1 + std::exp(current_state)); + MY_LOG(1 + MY_EXP(current_state)); log_prob -= (threshold_alpha + threshold_beta) * - std::log(1 + std::exp(proposed_state)); + MY_LOG(1 + MY_EXP(proposed_state)); U = R::unif_rand(); - if(std::log(U) < log_prob) { + if(MY_LOG(U) < log_prob) { thresholds(variable, 1) = proposed_state; } @@ -368,11 +383,11 @@ void metropolis_thresholds_blumecapel( if(log_prob > 0) { log_prob = 1; } else { - log_prob = std::exp(log_prob); + log_prob = MY_EXP(log_prob); } update_proposal_sd = proposal_sd_blumecapel(variable, 1) + - (log_prob - target_ar) * std::exp(-log(t) * phi); + (log_prob - target_ar) * MY_EXP(-log(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -383,6 +398,114 @@ void metropolis_thresholds_blumecapel( } +double log_pseudolikelihood_ratio( + const NumericMatrix& interactions, + const NumericMatrix& thresholds, + const IntegerMatrix& observations, + const IntegerVector& no_categories, + const int no_persons, + const int variable1, + const int variable2, + const double proposed_state, + const double current_state, + const NumericMatrix& rest_matrix, + const LogicalVector& variable_bool, + const IntegerVector& reference_category +) { + double rest_score, bound; + double pseudolikelihood_ratio = 0.0; + double denominator_prop, denominator_curr, exponent; + int score, obs_score1, obs_score2; + + double delta_state = proposed_state - current_state; + + int nc1 = no_categories[variable1]; + int nc2 = no_categories[variable2]; + bool is_ord1 = variable_bool[variable1]; + bool is_ord2 = variable_bool[variable2]; + int ref_cat1 = reference_category[variable1]; + int ref_cat2 = reference_category[variable2]; + + for (int person = 0; person < no_persons; person++) { + obs_score1 = observations(person, variable1); + obs_score2 = observations(person, variable2); + + pseudolikelihood_ratio += 2 * obs_score1 * obs_score2 * delta_state; + + // Variable 1 + rest_score = rest_matrix(person, variable1) - + obs_score2 * interactions(variable2, variable1); + bound = rest_score > 0.0 ? nc1 * rest_score : 0.0; + + double obs2_prop_state = obs_score2 * proposed_state; + double obs2_curr_state = obs_score2 * current_state; + + if (is_ord1) { + denominator_prop = MY_EXP(-bound); + denominator_curr = MY_EXP(-bound); + int index = variable1; + for (int category = 0; category < nc1; category++) { + score = category + 1; + exponent = thresholds[index] + score * rest_score - bound; + index += thresholds.nrow(); + double exp_term = exponent; + denominator_prop += MY_EXP(exp_term + score * obs2_prop_state); + denominator_curr += MY_EXP(exp_term + score * obs2_curr_state); + } + } else { + denominator_prop = 0.0; + denominator_curr = 0.0; + for (int category = 0; category <= nc1; category++) { + int diff = category - ref_cat1; + exponent = thresholds(variable1, 0) * category + + thresholds(variable1, 1) * diff * diff + + category * rest_score - bound; + denominator_prop += MY_EXP(exponent + category * obs2_prop_state); + denominator_curr += MY_EXP(exponent + category * obs2_curr_state); + } + } + + pseudolikelihood_ratio += MY_LOG(denominator_curr / denominator_prop); + + // Variable 2 + rest_score = rest_matrix(person, variable2) - + obs_score1 * interactions(variable1, variable2); + bound = rest_score > 0.0 ? nc2 * rest_score : 0.0; + + double obs1_prop_state = obs_score1 * proposed_state; + double obs1_curr_state = obs_score1 * current_state; + + if (is_ord2) { + denominator_prop = MY_EXP(-bound); + denominator_curr = MY_EXP(-bound); + int index = variable2; + for (int category = 0; category < nc2; category++) { + score = category + 1; + exponent = thresholds[index] + score * rest_score - bound; + index += thresholds.nrow(); + denominator_prop += MY_EXP(exponent + score * obs1_prop_state); + denominator_curr += MY_EXP(exponent + score * obs1_curr_state); + } + } else { + denominator_prop = 0.0; + denominator_curr = 0.0; + for (int category = 0; category <= nc2; category++) { + int diff = category - ref_cat2; + exponent = thresholds(variable2, 0) * category + + thresholds(variable2, 1) * diff * diff + + category * rest_score - bound; + denominator_prop += MY_EXP(exponent + category * obs1_prop_state); + denominator_curr += MY_EXP(exponent + category * obs1_curr_state); + } + } + + pseudolikelihood_ratio += MY_LOG(denominator_curr / denominator_prop); + } + + return pseudolikelihood_ratio; +} + +/* // ----------------------------------------------------------------------------| // The log pseudolikelihood ratio [proposed against current] for an interaction // ----------------------------------------------------------------------------| @@ -425,17 +548,17 @@ double log_pseudolikelihood_ratio( if(variable_bool[variable1] == true) { //Regular binary or ordinal MRF variable --------------------------------- - denominator_prop = std::exp(-bound); - denominator_curr = std::exp(-bound); + denominator_prop = MY_EXP(-bound); + denominator_curr = MY_EXP(-bound); for(int category = 0; category < no_categories[variable1]; category++) { score = category + 1; exponent = thresholds(variable1, category) + score * rest_score - bound; denominator_prop += - std::exp(exponent + score * obs_score2 * proposed_state); + MY_EXP(exponent + score * obs_score2 * proposed_state); denominator_curr += - std::exp(exponent + score * obs_score2 * current_state); + MY_EXP(exponent + score * obs_score2 * current_state); } } else { //Blume-Capel ordinal MRF variable --------------------------------------- @@ -448,13 +571,13 @@ double log_pseudolikelihood_ratio( (category - reference_category[variable1]); exponent+= category * rest_score - bound; denominator_prop += - std::exp(exponent + category * obs_score2 * proposed_state); + MY_EXP(exponent + category * obs_score2 * proposed_state); denominator_curr += - std::exp(exponent + category * obs_score2 * current_state); + MY_EXP(exponent + category * obs_score2 * current_state); } } - pseudolikelihood_ratio -= std::log(denominator_prop); - pseudolikelihood_ratio += std::log(denominator_curr); + pseudolikelihood_ratio -= MY_LOG(denominator_prop); + pseudolikelihood_ratio += MY_LOG(denominator_curr); //variable 2 log pseudolikelihood ratio rest_score = rest_matrix(person, variable2) - @@ -468,17 +591,18 @@ double log_pseudolikelihood_ratio( if(variable_bool[variable2] == true) { //Regular binary or ordinal MRF variable --------------------------------- - denominator_prop = std::exp(-bound); - denominator_curr = std::exp(-bound); + denominator_prop = MY_EXP(-bound); + denominator_curr = MY_EXP(-bound); + for(int category = 0; category < no_categories[variable2]; category++) { score = category + 1; exponent = thresholds(variable2, category) + score * rest_score - bound; denominator_prop += - std::exp(exponent + score * obs_score1 * proposed_state); + MY_EXP(exponent + score * obs_score1 * proposed_state); denominator_curr += - std::exp(exponent + score * obs_score1 * current_state); + MY_EXP(exponent + score * obs_score1 * current_state); } } else { //Blume-Capel ordinal MRF variable --------------------------------------- @@ -491,17 +615,17 @@ double log_pseudolikelihood_ratio( (category - reference_category[variable2]); exponent+= category * rest_score - bound; denominator_prop += - std::exp(exponent + category * obs_score1 * proposed_state); + MY_EXP(exponent + category * obs_score1 * proposed_state); denominator_curr += - std::exp(exponent + category * obs_score1 * current_state); + MY_EXP(exponent + category * obs_score1 * current_state); } } - pseudolikelihood_ratio -= std::log(denominator_prop); - pseudolikelihood_ratio += std::log(denominator_curr); + pseudolikelihood_ratio -= MY_LOG(denominator_prop); + pseudolikelihood_ratio += MY_LOG(denominator_curr); } return pseudolikelihood_ratio; } - +*/ // ----------------------------------------------------------------------------| // MH algorithm to sample from the full-conditional of the active interaction // parameters for Bayesian edge selection @@ -552,7 +676,7 @@ void metropolis_interactions( log_prob -= R::dcauchy(current_state, 0.0, interaction_scale, true); U = R::unif_rand(); - if(std::log(U) < log_prob) { + if(MY_LOG(U) < log_prob) { double state_diff = proposed_state - current_state; interactions(variable1, variable2) = proposed_state; interactions(variable2, variable1) = proposed_state; @@ -569,11 +693,11 @@ void metropolis_interactions( if(log_prob > 0) { log_prob = 1; } else { - log_prob = std::exp(log_prob); + log_prob = MY_EXP(log_prob); } double update_proposal_sd = proposal_sd(variable1, variable2) + - (log_prob - target_ar) * std::exp(-log(t) * phi); + (log_prob - target_ar) * MY_EXP(-log(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -587,6 +711,7 @@ void metropolis_interactions( } } + // ----------------------------------------------------------------------------| // MH algorithm to sample from the full-conditional of an edge + interaction // pair for Bayesian edge selection @@ -659,7 +784,7 @@ void metropolis_edge_interaction_pair( } U = R::unif_rand(); - if(std::log(U) < log_prob) { + if(MY_LOG(U) < log_prob) { indicator(variable1, variable2) = 1 - indicator(variable1, variable2); indicator(variable2, variable1) = 1 - indicator(variable2, variable1); @@ -679,6 +804,7 @@ void metropolis_edge_interaction_pair( } } + // ----------------------------------------------------------------------------| // A Gibbs step for graphical model parameters for Bayesian edge selection // ----------------------------------------------------------------------------| From 4970d4d931ed9e314db9273da90974ca864123be Mon Sep 17 00:00:00 2001 From: vandenman Date: Tue, 27 May 2025 14:01:48 +0200 Subject: [PATCH 02/10] uncomment debug stuff in Makevars --- src/Makevars.win | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Makevars.win b/src/Makevars.win index f45d7333..3c8dc8e2 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,5 +1,5 @@ -DEBUG=1 +#DEBUG=1 PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -PKG_CXXFLAGS+= -g -PKG_DLLFLAGS= -PKG_SHLIB_LDFLAGS= \ No newline at end of file +#PKG_CXXFLAGS+= -g +#PKG_DLLFLAGS= +#PKG_SHLIB_LDFLAGS= \ No newline at end of file From ea9eb6c9237c9db5cd9a793a92f499dff466d059 Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Tue, 27 May 2025 14:10:35 +0200 Subject: [PATCH 03/10] replace a few more calls --- src/gibbs_functions.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/gibbs_functions.cpp b/src/gibbs_functions.cpp index 823c71e9..aa0e2a8c 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -183,7 +183,7 @@ void metropolis_thresholds_regular( b = no_persons + threshold_beta - n_cat_obs(category + 1, variable); tmp = R::rbeta(a, b); proposed_state = MY_LOG(tmp / (1 - tmp) / c); - exp_proposed = exp(proposed_state); + exp_proposed = MY_EXP(proposed_state); //Compute log_acceptance probability for Metropolis. //First, we use g and q above to compute the ratio of pseudolikelihoods @@ -306,7 +306,7 @@ void metropolis_thresholds_blumecapel( } double update_proposal_sd = proposal_sd_blumecapel(variable, 0) + - (log_prob - target_ar) * MY_EXP(-log(t) * phi); + (log_prob - target_ar) * MY_EXP(-MY_LOG(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -387,7 +387,7 @@ void metropolis_thresholds_blumecapel( } update_proposal_sd = proposal_sd_blumecapel(variable, 1) + - (log_prob - target_ar) * MY_EXP(-log(t) * phi); + (log_prob - target_ar) * MY_EXP(-MY_LOG(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -697,7 +697,7 @@ void metropolis_interactions( } double update_proposal_sd = proposal_sd(variable1, variable2) + - (log_prob - target_ar) * MY_EXP(-log(t) * phi); + (log_prob - target_ar) * MY_EXP(-MY_LOG(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -772,7 +772,7 @@ void metropolis_edge_interaction_pair( proposal_sd(variable1, variable2), true); - log_prob += log(theta(variable1, variable2) / (1 - theta(variable1, variable2))); + log_prob += MY_LOG(theta(variable1, variable2) / (1 - theta(variable1, variable2))); } else { log_prob -= R::dcauchy(current_state, 0.0, interaction_scale, true); log_prob += R::dnorm(current_state, @@ -780,7 +780,7 @@ void metropolis_edge_interaction_pair( proposal_sd(variable1, variable2), true); - log_prob -= log(theta(variable1, variable2) / (1 - theta(variable1, variable2))); + log_prob -= MY_LOG(theta(variable1, variable2) / (1 - theta(variable1, variable2))); } U = R::unif_rand(); From acecd4320559f73d435a2b9b15b1f3753d977563 Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Thu, 5 Jun 2025 11:12:55 +0200 Subject: [PATCH 04/10] use exp from openlibm --- R/RcppExports.R | 4 + bgms.Rproj | 1 + src/RcppExports.cpp | 12 ++ src/e_exp.cpp | 303 ++++++++++++++++++++++++++++++++++++++++ src/e_exp.h | 6 + src/gibbs_functions.cpp | 3 +- 6 files changed, 328 insertions(+), 1 deletion(-) create mode 100644 src/e_exp.cpp create mode 100644 src/e_exp.h diff --git a/R/RcppExports.R b/R/RcppExports.R index 0d64db97..13e7cba1 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,6 +9,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) } +rcpp_ieee754_exp <- function(x) { + .Call(`_bgms_rcpp_ieee754_exp`, x) +} + gibbs_sampler <- function(observations, indicator, interactions, thresholds, no_categories, interaction_scale, proposal_sd, proposal_sd_blumecapel, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, n_cat_obs, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, variable_bool, reference_category, save = FALSE, display_progress = FALSE, edge_selection = TRUE) { .Call(`_bgms_gibbs_sampler`, observations, indicator, interactions, thresholds, no_categories, interaction_scale, proposal_sd, proposal_sd_blumecapel, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, n_cat_obs, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, variable_bool, reference_category, save, display_progress, edge_selection) } diff --git a/bgms.Rproj b/bgms.Rproj index 6291372e..73e9df0c 100644 --- a/bgms.Rproj +++ b/bgms.Rproj @@ -16,5 +16,6 @@ StripTrailingWhitespace: Yes BuildType: Package PackageUseDevtools: Yes +PackageCleanBeforeInstall: No PackageInstallArgs: --no-multiarch --with-keep.source PackageCheckArgs: --as-cran diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3b71f51e..947e3a98 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -45,6 +45,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// rcpp_ieee754_exp +Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x); +RcppExport SEXP _bgms_rcpp_ieee754_exp(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_exp(x)); + return rcpp_result_gen; +END_RCPP +} // gibbs_sampler List gibbs_sampler(IntegerMatrix& observations, IntegerMatrix& indicator, NumericMatrix& interactions, NumericMatrix& thresholds, const IntegerVector& no_categories, const double interaction_scale, NumericMatrix& proposal_sd, NumericMatrix& proposal_sd_blumecapel, const String& edge_prior, NumericMatrix& theta, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double dirichlet_alpha, const double lambda, const IntegerMatrix& Index, const int iter, const int burnin, IntegerMatrix& n_cat_obs, IntegerMatrix& sufficient_blume_capel, const double threshold_alpha, const double threshold_beta, const bool na_impute, const IntegerMatrix& missing_index, const LogicalVector& variable_bool, const IntegerVector& reference_category, const bool save, const bool display_progress, bool edge_selection); RcppExport SEXP _bgms_gibbs_sampler(SEXP observationsSEXP, SEXP indicatorSEXP, SEXP interactionsSEXP, SEXP thresholdsSEXP, SEXP no_categoriesSEXP, SEXP interaction_scaleSEXP, SEXP proposal_sdSEXP, SEXP proposal_sd_blumecapelSEXP, SEXP edge_priorSEXP, SEXP thetaSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP IndexSEXP, SEXP iterSEXP, SEXP burninSEXP, SEXP n_cat_obsSEXP, SEXP sufficient_blume_capelSEXP, SEXP threshold_alphaSEXP, SEXP threshold_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP variable_boolSEXP, SEXP reference_categorySEXP, SEXP saveSEXP, SEXP display_progressSEXP, SEXP edge_selectionSEXP) { @@ -145,6 +156,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_bgms_sample_omrf_gibbs", (DL_FUNC) &_bgms_sample_omrf_gibbs, 6}, {"_bgms_sample_bcomrf_gibbs", (DL_FUNC) &_bgms_sample_bcomrf_gibbs, 8}, + {"_bgms_rcpp_ieee754_exp", (DL_FUNC) &_bgms_rcpp_ieee754_exp, 1}, {"_bgms_gibbs_sampler", (DL_FUNC) &_bgms_gibbs_sampler, 28}, {"_bgms_compare_anova_gibbs_sampler", (DL_FUNC) &_bgms_compare_anova_gibbs_sampler, 34}, {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4}, diff --git a/src/e_exp.cpp b/src/e_exp.cpp new file mode 100644 index 00000000..8a8ae9a6 --- /dev/null +++ b/src/e_exp.cpp @@ -0,0 +1,303 @@ + +/* @(#)e_exp.c 1.6 04/04/22 */ +/* + * ==================================================== + * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. + * + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + + +// this is a stripped version of +// https://github.com/JuliaMath/openlibm/blob/c9c6fd6eade60ba38427c675a7323b0d2bae6521/src/e_exp.c +// specifically, the header files have been reduced to include only the minimially necessary things +// throughout the file EDIT indicates a modification relative to the original source code + +// EDIT: omitted header +// #include "cdefs-compat.h" + +//__FBSDID("$FreeBSD: src/lib/msun/src/e_exp.c,v 1.14 2011/10/21 06:26:38 das Exp $"); + +/* __ieee754_exp(x) + * Returns the exponential of x. + * + * Method + * 1. Argument reduction: + * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. + * Given x, find r and integer k such that + * + * x = k*ln2 + r, |r| <= 0.5*ln2. + * + * Here r will be represented as r = hi-lo for better + * accuracy. + * + * 2. Approximation of exp(r) by a special rational function on + * the interval [0,0.34658]: + * Write + * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... + * We use a special Remes algorithm on [0,0.34658] to generate + * a polynomial of degree 5 to approximate R. The maximum error + * of this polynomial approximation is bounded by 2**-59. In + * other words, + * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 + * (where z=r*r, and the values of P1 to P5 are listed below) + * and + * | 5 | -59 + * | 2.0+P1*z+...+P5*z - R(z) | <= 2 + * | | + * The computation of exp(r) thus becomes + * 2*r + * exp(r) = 1 + ------- + * R - r + * r*R1(r) + * = 1 + r + ----------- (for better accuracy) + * 2 - R1(r) + * where + * 2 4 10 + * R1(r) = r - (P1*r + P2*r + ... + P5*r ). + * + * 3. Scale back to obtain exp(x): + * From step 1, we have + * exp(x) = 2^k * exp(r) + * + * Special cases: + * exp(INF) is INF, exp(NaN) is NaN; + * exp(-INF) is 0, and + * for finite argument, only exp(0)=1 is exact. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Misc. info. + * For IEEE double + * if x > 7.09782712893383973096e+02 then exp(x) overflow + * if x < -7.45133219101941108420e+02 then exp(x) underflow + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +// EDIT: omitted headers, the necessary definitions are pasted below +// #include +// #include "math_private.h" + +/* Definitions provided directly by GCC and Clang. */ +#if !(defined(__BYTE_ORDER__) && defined(__ORDER_LITTLE_ENDIAN__) && defined(__ORDER_BIG_ENDIAN__)) + +#if defined(__GLIBC__) + +#include +#include +#define __ORDER_LITTLE_ENDIAN__ __LITTLE_ENDIAN +#define __ORDER_BIG_ENDIAN__ __BIG_ENDIAN +#define __BYTE_ORDER__ __BYTE_ORDER + +#elif defined(__APPLE__) + +#include +#define __ORDER_LITTLE_ENDIAN__ LITTLE_ENDIAN +#define __ORDER_BIG_ENDIAN__ BIG_ENDIAN +#define __BYTE_ORDER__ BYTE_ORDER + +#elif defined(_WIN32) + +#define __ORDER_LITTLE_ENDIAN__ 1234 +#define __ORDER_BIG_ENDIAN__ 4321 +#define __BYTE_ORDER__ __ORDER_LITTLE_ENDIAN__ + +#endif + +#endif /* __BYTE_ORDER__, __ORDER_LITTLE_ENDIAN__ and __ORDER_BIG_ENDIAN__ */ + +#ifndef __FLOAT_WORD_ORDER__ +#define __FLOAT_WORD_ORDER__ __BYTE_ORDER__ +#endif + + +#include +// #include // EDIT: for int32_t, u_int32_t +#include // EDIT: for u_int64_t + +//EDIT: types-compat.h +typedef uint8_t u_int8_t; +typedef uint16_t u_int16_t; +typedef uint32_t u_int32_t; +typedef uint64_t u_int64_t; + + +/* + * A union which permits us to convert between a double and two 32 bit + * ints. + */ + + +#if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__ + +typedef union +{ + double value; + struct + { + u_int32_t msw; + u_int32_t lsw; + } parts; + struct + { + u_int64_t w; + } xparts; +} ieee_double_shape_type; + +#endif + +#if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__ + +typedef union +{ + double value; + struct + { + u_int32_t lsw; + u_int32_t msw; + } parts; + struct + { + u_int64_t w; + } xparts; +} ieee_double_shape_type; + +#endif + +/* Get the more significant 32 bit int from a double. */ + +#define GET_HIGH_WORD(i,d) \ +do { \ + ieee_double_shape_type gh_u; \ + gh_u.value = (d); \ + (i) = gh_u.parts.msw; \ +} while (0) + +/* Get the less significant 32 bit int from a double. */ + +#define GET_LOW_WORD(i,d) \ +do { \ + ieee_double_shape_type gl_u; \ + gl_u.value = (d); \ + (i) = gl_u.parts.lsw; \ +} while (0) + +/* Set a double from two 32 bit ints. */ + +#define INSERT_WORDS(d,ix0,ix1) \ +do { \ + ieee_double_shape_type iw_u; \ + iw_u.parts.msw = (ix0); \ + iw_u.parts.lsw = (ix1); \ + (d) = iw_u.value; \ +} while (0) + + +static const double +one = 1.0, +halF[2] = {0.5,-0.5,}, +huge = 1.0e+300, +o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ +u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ +ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ + -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ +ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ + -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ +invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ +P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ +P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ +P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ +P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ +P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ + +static volatile double +twom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/ + +double __ieee754_exp(double x) /* default IEEE double exp */ +{ + double y,hi=0.0,lo=0.0,c,t,twopk; + int32_t k=0,xsb; + u_int32_t hx; + + GET_HIGH_WORD(hx,x); + xsb = (hx>>31)&1; /* sign bit of x */ + hx &= 0x7fffffff; /* high word of |x| */ + + /* filter out non-finite argument */ + if(hx >= 0x40862E42) { /* if |x|>=709.78... */ + if(hx>=0x7ff00000) { + u_int32_t lx; + GET_LOW_WORD(lx,x); + if(((hx&0xfffff)|lx)!=0) + return x+x; /* NaN */ + else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ + } + if(x > o_threshold) return huge*huge; /* overflow */ + if(x < u_threshold) return twom1000*twom1000; /* underflow */ + } + + /* this implementation gives 2.7182818284590455 for exp(1.0), + which is well within the allowable error. however, + 2.718281828459045 is closer to the true value so we prefer that + answer, given that 1.0 is such an important argument value. */ + if (x == 1.0) + return 2.718281828459045235360; + + /* argument reduction */ + if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ + if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ + hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; + } else { + k = (int)(invln2*x+halF[xsb]); + t = k; + hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ + lo = t*ln2LO[0]; + } + //STRICT_ASSIGN(double, x, hi - lo); + x = hi - lo; // EDIT: this is equivalent to the line above on platforms we care about + } + else if(hx < 0x3e300000) { /* when |x|<2**-28 */ + if(huge+x>one) return one+x;/* trigger inexact */ + } + else k = 0; + + /* x is now in primary range */ + t = x*x; + if(k >= -1021) + INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0); + else + INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0); + c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); + if(k==0) return one-((x*c)/(c-2.0)-x); + else y = one-((lo-(x*c)/(2.0-c))-hi); + if(k >= -1021) { + if (k==1024) return y*2.0*0x1p1023; + return y*twopk; + } else { + return y*twopk*twom1000; + } +} + +// #if (LDBL_MANT_DIG == 53) +// openlibm_weak_reference(exp, expl); +// #endif + +#include "Rcpp.h" +// [[Rcpp::export]] +Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x) { + Rcpp::NumericVector y(x.size()); + for (int i = 0; i < x.size(); i++) { + y[i] = __ieee754_exp(x[i]); + } + return y; +} \ No newline at end of file diff --git a/src/e_exp.h b/src/e_exp.h new file mode 100644 index 00000000..1bfd0c05 --- /dev/null +++ b/src/e_exp.h @@ -0,0 +1,6 @@ +#ifndef _E_EXP_H_ +#define _E_EXP_H_ + +double __ieee754_exp(double x); + +#endif diff --git a/src/gibbs_functions.cpp b/src/gibbs_functions.cpp index aa0e2a8c..8df9c931 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -1,6 +1,7 @@ // [[Rcpp::depends(RcppProgress)]] #include #include "gibbs_functions_edge_prior.h" +#include "e_exp.h" #include #include using namespace Rcpp; @@ -17,7 +18,7 @@ inline double fast_log_bitwise(double x) { return log2 * 0.69314718 + (y - 1) - (y - 1)*(y - 1)/2; // linear/quadratic correction } -#define MY_EXP fast_exp +#define MY_EXP __ieee754_exp//fast_exp #define MY_LOG fast_log_bitwise // ----------------------------------------------------------------------------| From c5315ce9f124e1fc6d830b6ecaf8f30671d32aeb Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Thu, 5 Jun 2025 11:32:00 +0200 Subject: [PATCH 05/10] use log from openlibm --- R/RcppExports.R | 4 ++ src/RcppExports.cpp | 12 +++++ src/e_exp.cpp | 106 ++++++++++++++++++++++++++++++++++++++++- src/e_exp.h | 1 + validate_exp_and_log.R | 57 ++++++++++++++++++++++ 5 files changed, 179 insertions(+), 1 deletion(-) create mode 100644 validate_exp_and_log.R diff --git a/R/RcppExports.R b/R/RcppExports.R index 13e7cba1..e5feda93 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,6 +13,10 @@ rcpp_ieee754_exp <- function(x) { .Call(`_bgms_rcpp_ieee754_exp`, x) } +rcpp_ieee754_log <- function(x) { + .Call(`_bgms_rcpp_ieee754_log`, x) +} + gibbs_sampler <- function(observations, indicator, interactions, thresholds, no_categories, interaction_scale, proposal_sd, proposal_sd_blumecapel, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, n_cat_obs, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, variable_bool, reference_category, save = FALSE, display_progress = FALSE, edge_selection = TRUE) { .Call(`_bgms_gibbs_sampler`, observations, indicator, interactions, thresholds, no_categories, interaction_scale, proposal_sd, proposal_sd_blumecapel, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, n_cat_obs, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, variable_bool, reference_category, save, display_progress, edge_selection) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 947e3a98..7cf15ca4 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -56,6 +56,17 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// rcpp_ieee754_log +Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x); +RcppExport SEXP _bgms_rcpp_ieee754_log(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_log(x)); + return rcpp_result_gen; +END_RCPP +} // gibbs_sampler List gibbs_sampler(IntegerMatrix& observations, IntegerMatrix& indicator, NumericMatrix& interactions, NumericMatrix& thresholds, const IntegerVector& no_categories, const double interaction_scale, NumericMatrix& proposal_sd, NumericMatrix& proposal_sd_blumecapel, const String& edge_prior, NumericMatrix& theta, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double dirichlet_alpha, const double lambda, const IntegerMatrix& Index, const int iter, const int burnin, IntegerMatrix& n_cat_obs, IntegerMatrix& sufficient_blume_capel, const double threshold_alpha, const double threshold_beta, const bool na_impute, const IntegerMatrix& missing_index, const LogicalVector& variable_bool, const IntegerVector& reference_category, const bool save, const bool display_progress, bool edge_selection); RcppExport SEXP _bgms_gibbs_sampler(SEXP observationsSEXP, SEXP indicatorSEXP, SEXP interactionsSEXP, SEXP thresholdsSEXP, SEXP no_categoriesSEXP, SEXP interaction_scaleSEXP, SEXP proposal_sdSEXP, SEXP proposal_sd_blumecapelSEXP, SEXP edge_priorSEXP, SEXP thetaSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP IndexSEXP, SEXP iterSEXP, SEXP burninSEXP, SEXP n_cat_obsSEXP, SEXP sufficient_blume_capelSEXP, SEXP threshold_alphaSEXP, SEXP threshold_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP variable_boolSEXP, SEXP reference_categorySEXP, SEXP saveSEXP, SEXP display_progressSEXP, SEXP edge_selectionSEXP) { @@ -157,6 +168,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bgms_sample_omrf_gibbs", (DL_FUNC) &_bgms_sample_omrf_gibbs, 6}, {"_bgms_sample_bcomrf_gibbs", (DL_FUNC) &_bgms_sample_bcomrf_gibbs, 8}, {"_bgms_rcpp_ieee754_exp", (DL_FUNC) &_bgms_rcpp_ieee754_exp, 1}, + {"_bgms_rcpp_ieee754_log", (DL_FUNC) &_bgms_rcpp_ieee754_log, 1}, {"_bgms_gibbs_sampler", (DL_FUNC) &_bgms_gibbs_sampler, 28}, {"_bgms_compare_anova_gibbs_sampler", (DL_FUNC) &_bgms_compare_anova_gibbs_sampler, 34}, {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4}, diff --git a/src/e_exp.cpp b/src/e_exp.cpp index 8a8ae9a6..702d2cc1 100644 --- a/src/e_exp.cpp +++ b/src/e_exp.cpp @@ -88,6 +88,7 @@ // #include // #include "math_private.h" +// EDIT: fpmath.h /* Definitions provided directly by GCC and Clang. */ #if !(defined(__BYTE_ORDER__) && defined(__ORDER_LITTLE_ENDIAN__) && defined(__ORDER_BIG_ENDIAN__)) @@ -137,7 +138,7 @@ typedef uint64_t u_int64_t; * ints. */ - +// the defines below are from "math_private.h" #if __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__ typedef union @@ -202,7 +203,28 @@ do { \ (d) = iw_u.value; \ } while (0) +/* Set the more significant 32 bits of a double from an int. */ + +#define SET_HIGH_WORD(d,v) \ +do { \ + ieee_double_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts.msw = (v); \ + (d) = sh_u.value; \ +} while (0) \ + +/* Get two 32 bit ints from a double. */ + +#define EXTRACT_WORDS(ix0,ix1,d) \ +do { \ + ieee_double_shape_type ew_u; \ + ew_u.value = (d); \ + (ix0) = ew_u.parts.msw; \ + (ix1) = ew_u.parts.lsw; \ +} while (0) + +// EDIT: e_exp.c static const double one = 1.0, halF[2] = {0.5,-0.5,}, @@ -292,6 +314,78 @@ double __ieee754_exp(double x) /* default IEEE double exp */ // openlibm_weak_reference(exp, expl); // #endif +// EDIT: e_log.c +static const double + ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ + ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ + two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ + Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ + Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ + Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ + Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ + Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ + Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ + Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + + static const double zero = 0.0; + +double __ieee754_log(double x) + { + double hfsq,f,s,z,R,w,t1,t2,dk; + int32_t k,hx,i,j; + u_int32_t lx; + + EXTRACT_WORDS(hx,lx,x); + + k=0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx&0x7fffffff)|lx)==0) + return -two54/zero; /* log(+-0)=-inf */ + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 54; x *= two54; /* subnormal number, scale up x */ + GET_HIGH_WORD(hx,x); + } + if (hx >= 0x7ff00000) return x+x; + k += (hx>>20)-1023; + hx &= 0x000fffff; + i = (hx+0x95f64)&0x100000; + SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ + k += (i>>20); + f = x-1.0; + if((0x000fffff&(2+hx))<3) { /* -2**-20 <= f < 2**-20 */ + if(f==zero) { + if(k==0) { + return zero; + } else { + dk=(double)k; + return dk*ln2_hi+dk*ln2_lo; + } + } + R = f*f*(0.5-0.33333333333333333*f); + if(k==0) return f-R; else {dk=(double)k; + return dk*ln2_hi-((R-dk*ln2_lo)-f);} + } + s = f/(2.0+f); + dk = (double)k; + z = s*s; + i = hx-0x6147a; + w = z*z; + j = 0x6b851-hx; + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + i |= j; + R = t2+t1; + if(i>0) { + hfsq=0.5*f*f; + if(k==0) return f-(hfsq-s*(hfsq+R)); else + return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); + } else { + if(k==0) return f-s*(f-R); else + return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); + } +} + + #include "Rcpp.h" // [[Rcpp::export]] Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x) { @@ -300,4 +394,14 @@ Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x) { y[i] = __ieee754_exp(x[i]); } return y; +} + +#include "Rcpp.h" +// [[Rcpp::export]] +Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x) { + Rcpp::NumericVector y(x.size()); + for (int i = 0; i < x.size(); i++) { + y[i] = __ieee754_log(x[i]); + } + return y; } \ No newline at end of file diff --git a/src/e_exp.h b/src/e_exp.h index 1bfd0c05..b26f5ae9 100644 --- a/src/e_exp.h +++ b/src/e_exp.h @@ -2,5 +2,6 @@ #define _E_EXP_H_ double __ieee754_exp(double x); +double __ieee754_log(double x); #endif diff --git a/validate_exp_and_log.R b/validate_exp_and_log.R new file mode 100644 index 00000000..2cc7043e --- /dev/null +++ b/validate_exp_and_log.R @@ -0,0 +1,57 @@ + +values <- matrix(c( + -1e6, 1e6, + -1e4, 1e4, + -1e3, 1e3, + -1e2, 1e2, + -1e1, 1e1, + -1e-3, 1e-3 +), ncol = 2, byrow = TRUE) +size <- 1E6L + +for (i in seq_len(nrow(values))) { + from <- values[i, 1] + to <- values[i, 2] + cat("i = ", i, ", Testing range: ", from, " to ", to, "\n") + test_values <- seq(from, to, length.out = size) + y_R <- exp(test_values) + y_R2 <- bgms:::rcpp_ieee754_exp(test_values) + result <- all.equal(y_R, y_R2) + + abs_error = max(abs(y_R - y_R2), na.rm = TRUE) + rel_error = max(abs((y_R - y_R2) / y_R), na.rm = TRUE) + + cat("i = ", i, ", All.equal = ", isTRUE(result), + " max abs error:", sprintf("%.8g", abs_error), + " max rel error:", sprintf("%.8g", rel_error), + "\n") +} + +values <- matrix(c( + 1e10, 1e30, + 1e04, 1e09, + 1e01, 1e03, + 1e-2, 1e01, + 1e-5, 1e-3, + 1e-20, 1e-6 +), ncol = 2, byrow = TRUE) +size <- 1E6L + +for (i in seq_len(nrow(values))) { + from <- values[i, 1] + to <- values[i, 2] + cat("i = ", i, ", Testing range: ", from, " to ", to, "\n") + test_values <- seq(from, to, length.out = size) + y_R <- log(test_values) + y_R2 <- bgms:::rcpp_ieee754_log(test_values) + result <- all.equal(y_R, y_R2) + + abs_error = max(abs(y_R - y_R2), na.rm = TRUE) + rel_error = max(abs((y_R - y_R2) / y_R), na.rm = TRUE) + + cat("i = ", i, ", All.equal = ", isTRUE(result), + " max abs error:", sprintf("%.8g", abs_error), + " max rel error:", sprintf("%.8g", rel_error), + "\n") +} + From effeacf85cd64b15341bba5acd4416e5739ae6aa Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Thu, 5 Jun 2025 09:36:27 +0200 Subject: [PATCH 06/10] use log --- src/gibbs_functions.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gibbs_functions.cpp b/src/gibbs_functions.cpp index 8df9c931..ae3a3afe 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -19,7 +19,7 @@ inline double fast_log_bitwise(double x) { } #define MY_EXP __ieee754_exp//fast_exp -#define MY_LOG fast_log_bitwise +#define MY_LOG __ieee754_log // ----------------------------------------------------------------------------| // Impute missing data from full-conditional From 94073135df819098f396b0bacf034604a118c22d Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Tue, 17 Jun 2025 13:37:57 +0200 Subject: [PATCH 07/10] use it everywhere --- R/RcppExports.R | 16 +- src/RcppExports.cpp | 59 ++++--- src/custom_exp.cpp | 29 ++++ src/data_simulation.cpp | 7 +- src/e_exp.cpp | 21 --- src/explog_switch.h | 41 +++++ src/gibbs_functions.cpp | 16 +- src/gibbs_functions_compare.cpp | 239 +++++++++++++++-------------- src/gibbs_functions_edge_prior.cpp | 31 ++-- 9 files changed, 256 insertions(+), 203 deletions(-) create mode 100644 src/custom_exp.cpp create mode 100644 src/explog_switch.h diff --git a/R/RcppExports.R b/R/RcppExports.R index e5feda93..35aa5389 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,12 +1,8 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -sample_omrf_gibbs <- function(no_states, no_variables, no_categories, interactions, thresholds, iter) { - .Call(`_bgms_sample_omrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, iter) -} - -sample_bcomrf_gibbs <- function(no_states, no_variables, no_categories, interactions, thresholds, variable_type, reference_category, iter) { - .Call(`_bgms_sample_bcomrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, variable_type, reference_category, iter) +get_explog_switch <- function() { + .Call(`_bgms_get_explog_switch`) } rcpp_ieee754_exp <- function(x) { @@ -17,6 +13,14 @@ rcpp_ieee754_log <- function(x) { .Call(`_bgms_rcpp_ieee754_log`, x) } +sample_omrf_gibbs <- function(no_states, no_variables, no_categories, interactions, thresholds, iter) { + .Call(`_bgms_sample_omrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, iter) +} + +sample_bcomrf_gibbs <- function(no_states, no_variables, no_categories, interactions, thresholds, variable_type, reference_category, iter) { + .Call(`_bgms_sample_bcomrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, variable_type, reference_category, iter) +} + gibbs_sampler <- function(observations, indicator, interactions, thresholds, no_categories, interaction_scale, proposal_sd, proposal_sd_blumecapel, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, n_cat_obs, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, variable_bool, reference_category, save = FALSE, display_progress = FALSE, edge_selection = TRUE) { .Call(`_bgms_gibbs_sampler`, observations, indicator, interactions, thresholds, no_categories, interaction_scale, proposal_sd, proposal_sd_blumecapel, edge_prior, theta, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, Index, iter, burnin, n_cat_obs, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, variable_bool, reference_category, save, display_progress, edge_selection) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 7cf15ca4..95dac6a1 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,6 +11,38 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// get_explog_switch +Rcpp::String get_explog_switch(); +RcppExport SEXP _bgms_get_explog_switch() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(get_explog_switch()); + return rcpp_result_gen; +END_RCPP +} +// rcpp_ieee754_exp +Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x); +RcppExport SEXP _bgms_rcpp_ieee754_exp(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_exp(x)); + return rcpp_result_gen; +END_RCPP +} +// rcpp_ieee754_log +Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x); +RcppExport SEXP _bgms_rcpp_ieee754_log(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_log(x)); + return rcpp_result_gen; +END_RCPP +} // sample_omrf_gibbs IntegerMatrix sample_omrf_gibbs(int no_states, int no_variables, IntegerVector no_categories, NumericMatrix interactions, NumericMatrix thresholds, int iter); RcppExport SEXP _bgms_sample_omrf_gibbs(SEXP no_statesSEXP, SEXP no_variablesSEXP, SEXP no_categoriesSEXP, SEXP interactionsSEXP, SEXP thresholdsSEXP, SEXP iterSEXP) { @@ -45,28 +77,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// rcpp_ieee754_exp -Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x); -RcppExport SEXP _bgms_rcpp_ieee754_exp(SEXP xSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP); - rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_exp(x)); - return rcpp_result_gen; -END_RCPP -} -// rcpp_ieee754_log -Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x); -RcppExport SEXP _bgms_rcpp_ieee754_log(SEXP xSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::NumericVector >::type x(xSEXP); - rcpp_result_gen = Rcpp::wrap(rcpp_ieee754_log(x)); - return rcpp_result_gen; -END_RCPP -} // gibbs_sampler List gibbs_sampler(IntegerMatrix& observations, IntegerMatrix& indicator, NumericMatrix& interactions, NumericMatrix& thresholds, const IntegerVector& no_categories, const double interaction_scale, NumericMatrix& proposal_sd, NumericMatrix& proposal_sd_blumecapel, const String& edge_prior, NumericMatrix& theta, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double dirichlet_alpha, const double lambda, const IntegerMatrix& Index, const int iter, const int burnin, IntegerMatrix& n_cat_obs, IntegerMatrix& sufficient_blume_capel, const double threshold_alpha, const double threshold_beta, const bool na_impute, const IntegerMatrix& missing_index, const LogicalVector& variable_bool, const IntegerVector& reference_category, const bool save, const bool display_progress, bool edge_selection); RcppExport SEXP _bgms_gibbs_sampler(SEXP observationsSEXP, SEXP indicatorSEXP, SEXP interactionsSEXP, SEXP thresholdsSEXP, SEXP no_categoriesSEXP, SEXP interaction_scaleSEXP, SEXP proposal_sdSEXP, SEXP proposal_sd_blumecapelSEXP, SEXP edge_priorSEXP, SEXP thetaSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP IndexSEXP, SEXP iterSEXP, SEXP burninSEXP, SEXP n_cat_obsSEXP, SEXP sufficient_blume_capelSEXP, SEXP threshold_alphaSEXP, SEXP threshold_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP variable_boolSEXP, SEXP reference_categorySEXP, SEXP saveSEXP, SEXP display_progressSEXP, SEXP edge_selectionSEXP) { @@ -165,10 +175,11 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_bgms_sample_omrf_gibbs", (DL_FUNC) &_bgms_sample_omrf_gibbs, 6}, - {"_bgms_sample_bcomrf_gibbs", (DL_FUNC) &_bgms_sample_bcomrf_gibbs, 8}, + {"_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_gibbs_sampler", (DL_FUNC) &_bgms_gibbs_sampler, 28}, {"_bgms_compare_anova_gibbs_sampler", (DL_FUNC) &_bgms_compare_anova_gibbs_sampler, 34}, {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4}, diff --git a/src/custom_exp.cpp b/src/custom_exp.cpp new file mode 100644 index 00000000..dfe2c1d7 --- /dev/null +++ b/src/custom_exp.cpp @@ -0,0 +1,29 @@ +#include "Rcpp.h" +#include "explog_switch.h" + +// [[Rcpp::export]] +Rcpp::String get_explog_switch() { +#if USE_CUSTOM_LOG + return "custom"; + #else + return "standard"; + #endif +} + +// [[Rcpp::export]] +Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x) { + Rcpp::NumericVector y(x.size()); + for (int i = 0; i < x.size(); i++) { + y[i] = MY_EXP(x[i]); + } + return y; +} + +// [[Rcpp::export]] +Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x) { + Rcpp::NumericVector y(x.size()); + for (int i = 0; i < x.size(); i++) { + y[i] = MY_LOG(x[i]); + } + return y; +} \ No newline at end of file diff --git a/src/data_simulation.cpp b/src/data_simulation.cpp index 8686a821..b4aa539a 100644 --- a/src/data_simulation.cpp +++ b/src/data_simulation.cpp @@ -1,4 +1,5 @@ #include +#include "explog_switch.h" using namespace Rcpp; // [[Rcpp::export]] @@ -53,7 +54,7 @@ IntegerMatrix sample_omrf_gibbs(int no_states, for(int category = 0; category < no_categories[variable]; category++) { exponent = thresholds(variable, category); exponent += (category + 1) * rest_score; - cumsum += std::exp(exponent); + cumsum += MY_EXP(exponent); probabilities[category + 1] = cumsum; } @@ -132,7 +133,7 @@ IntegerMatrix sample_bcomrf_gibbs(int no_states, (category - reference_category[variable]); //The pairwise interactions exponent += category * rest_score; - cumsum += std::exp(exponent); + cumsum += MY_EXP(exponent); probabilities[category] = cumsum; } } else { @@ -141,7 +142,7 @@ IntegerMatrix sample_bcomrf_gibbs(int no_states, for(int category = 0; category < no_categories[variable]; category++) { exponent = thresholds(variable, category); exponent += (category + 1) * rest_score; - cumsum += std::exp(exponent); + cumsum += MY_EXP(exponent); probabilities[category + 1] = cumsum; } } diff --git a/src/e_exp.cpp b/src/e_exp.cpp index 702d2cc1..c39172ac 100644 --- a/src/e_exp.cpp +++ b/src/e_exp.cpp @@ -384,24 +384,3 @@ double __ieee754_log(double x) return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); } } - - -#include "Rcpp.h" -// [[Rcpp::export]] -Rcpp::NumericVector rcpp_ieee754_exp(Rcpp::NumericVector x) { - Rcpp::NumericVector y(x.size()); - for (int i = 0; i < x.size(); i++) { - y[i] = __ieee754_exp(x[i]); - } - return y; -} - -#include "Rcpp.h" -// [[Rcpp::export]] -Rcpp::NumericVector rcpp_ieee754_log(Rcpp::NumericVector x) { - Rcpp::NumericVector y(x.size()); - for (int i = 0; i < x.size(); i++) { - y[i] = __ieee754_log(x[i]); - } - return y; -} \ No newline at end of file diff --git a/src/explog_switch.h b/src/explog_switch.h new file mode 100644 index 00000000..e4a7e2fd --- /dev/null +++ b/src/explog_switch.h @@ -0,0 +1,41 @@ +#ifndef _EXPLOG_SWITCH_H_ +#define _EXPLOG_SWITCH_H_ + +/* + * To switch between these, try in R: + * Sys.setenv("PKG_CPPFLAGS" = "-DCUSTOM_EXP_LOG=0") # 0 for not using it, 1 for using it + * renv::install(".") # or whatever way you use to install bgms + */ + +#if (defined(_WIN32) && (!defined(CUSTOM_EXP_LOG) || CUSTOM_EXP_LOG == 0)) \ + || (!defined(_WIN32) && defined(CUSTOM_EXP_LOG) && CUSTOM_EXP_LOG != 0) + +#define USE_CUSTOM_LOG 1 + +#else + +#define USE_CUSTOM_LOG 0 + +#endif + +#if USE_CUSTOM_LOG + +#include "e_exp.h" + +#define MY_EXP __ieee754_exp +#define MY_LOG __ieee754_log + +// TODO: add and use these +// #define MY_EXPM1 std::expm1 +// #define MY_LOG1P std::log1p +// #define MY_LOG1PEXP + + +#else + +#define MY_EXP std::exp +#define MY_LOG std::log + +#endif + +#endif \ No newline at end of file diff --git a/src/gibbs_functions.cpp b/src/gibbs_functions.cpp index ae3a3afe..4a4e238b 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -1,25 +1,11 @@ // [[Rcpp::depends(RcppProgress)]] #include #include "gibbs_functions_edge_prior.h" -#include "e_exp.h" #include #include using namespace Rcpp; -inline double fast_exp(double x) { - // 5th-order Pade approximant - return 1.0 + x * (1.0 + x * (0.5 + x * (1.0/6.0 + x * (1.0/24.0 + x * (1.0/120.0))))); -} - -inline double fast_log_bitwise(double x) { - int* xp = reinterpret_cast(&x); - int log2 = ((xp[1] >> 20) & 0x7FF) - 1023; // Extract exponent bits - double y = x / std::ldexp(1.0, log2); // Normalize mantissa - return log2 * 0.69314718 + (y - 1) - (y - 1)*(y - 1)/2; // linear/quadratic correction -} - -#define MY_EXP __ieee754_exp//fast_exp -#define MY_LOG __ieee754_log +#include "explog_switch.h" // ----------------------------------------------------------------------------| // Impute missing data from full-conditional diff --git a/src/gibbs_functions_compare.cpp b/src/gibbs_functions_compare.cpp index 7b43de90..93aab08c 100644 --- a/src/gibbs_functions_compare.cpp +++ b/src/gibbs_functions_compare.cpp @@ -6,6 +6,7 @@ #include #include +#include "explog_switch.h" using namespace Rcpp; @@ -34,7 +35,7 @@ double update_with_robbins_monro ( // Normalize the acceptance probability double observed_acceptance_probability = 1.0; if(observed_log_acceptance_probability < 0.0) { - observed_acceptance_probability = std::exp(observed_log_acceptance_probability); + observed_acceptance_probability = MY_EXP(observed_log_acceptance_probability); } // Update the proposal standard deviation @@ -168,7 +169,7 @@ List impute_missing_data_for_anova_model( for(int category = 1; category <= num_categories(variable, gr); category++) { exponent = GroupThresholds(category - 1); exponent += category * rest_score; - cumsum += std::exp(exponent); + cumsum += MY_EXP(exponent); category_response_probabilities[category] = cumsum; } } else { @@ -180,7 +181,7 @@ List impute_missing_data_for_anova_model( (category - baseline_category[variable]) * (category - baseline_category[variable]); exponent += category * rest_score; - cumsum += std::exp(exponent); + cumsum += MY_EXP(exponent); category_response_probabilities[category] = cumsum; } } @@ -330,13 +331,13 @@ double log_pseudolikelihood_ratio_interaction( // Compute pseudo-likelihood terms if(is_ordinal_variable[var]) { // Regular binary or ordinal MRF variable - denominator_prop += std::exp(-bound); - denominator_curr += std::exp(-bound); + denominator_prop += MY_EXP(-bound); + denominator_curr += MY_EXP(-bound); for(int cat = 0; cat < n_cats; cat++) { score = cat + 1; exponent = GroupThresholds[cat] + score * rest_score - bound; - denominator_prop += std::exp(exponent + score * obs_proposed); - denominator_curr += std::exp(exponent + score * obs_current); + denominator_prop += MY_EXP(exponent + score * obs_proposed); + denominator_curr += MY_EXP(exponent + score * obs_current); } } else { // Blume-Capel ordinal MRF variable @@ -346,14 +347,14 @@ double log_pseudolikelihood_ratio_interaction( (cat - baseline_category[var]) * (cat - baseline_category[var]); exponent += cat * rest_score - bound; - denominator_prop += std::exp(exponent + cat * obs_proposed); - denominator_curr += std::exp(exponent + cat * obs_current); + denominator_prop += MY_EXP(exponent + cat * obs_proposed); + denominator_curr += MY_EXP(exponent + cat * obs_current); } } // Update the pseudo-likelihood ratio - pseudolikelihood_ratio -= std::log(denominator_prop); - pseudolikelihood_ratio += std::log(denominator_curr); + pseudolikelihood_ratio -= MY_LOG(denominator_prop); + pseudolikelihood_ratio += MY_LOG(denominator_curr); } } } @@ -443,7 +444,7 @@ void metropolis_interaction( // Metropolis-Hastings acceptance step U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { // Update the interaction parameter state_difference = proposed_state - current_state; pairwise_effects(int_index, 0) = proposed_state; @@ -570,13 +571,13 @@ double log_pseudolikelihood_ratio_pairwise_difference( // Compute denominators based on whether the variable is ordinal if (is_ordinal_variable[var]) { // Binary or ordinal MRF variable - denominator_prop = std::exp(-bound); - denominator_curr = std::exp(-bound); + denominator_prop = MY_EXP(-bound); + denominator_curr = MY_EXP(-bound); for (int cat = 0; cat < n_cats; cat++) { score = cat + 1; exponent = GroupThresholds[cat] + score * rest_score - bound; - denominator_prop += std::exp(exponent + score * obs_proposed_p); - denominator_curr += std::exp(exponent + score * obs_current_p); + denominator_prop += MY_EXP(exponent + score * obs_proposed_p); + denominator_curr += MY_EXP(exponent + score * obs_current_p); } } else { // Blume-Capel ordinal MRF variable @@ -587,14 +588,14 @@ double log_pseudolikelihood_ratio_pairwise_difference( GroupThresholds[1] * (cat - baseline_category[var]) * (cat - baseline_category[var]) + cat * rest_score - bound; - denominator_prop += std::exp(exponent + cat * obs_proposed_p); - denominator_curr += std::exp(exponent + cat * obs_current_p); + denominator_prop += MY_EXP(exponent + cat * obs_proposed_p); + denominator_curr += MY_EXP(exponent + cat * obs_current_p); } } // Update the pseudo-likelihood ratio - pseudolikelihood_ratio -= std::log(denominator_prop); - pseudolikelihood_ratio += std::log(denominator_curr); + pseudolikelihood_ratio -= MY_LOG(denominator_prop); + pseudolikelihood_ratio += MY_LOG(denominator_curr); } } } @@ -688,7 +689,7 @@ void metropolis_pairwise_difference( // Metropolis-Hastings acceptance step U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { pairwise_effects(int_index, h) = proposed_state; // Update the rest matrix to reflect the new pairwise difference @@ -820,14 +821,14 @@ double log_pseudolikelihood_ratio_pairwise_differences( // Compute denominators if (is_ordinal_variable[var]) { - denominator_prop += std::exp(-bound); - denominator_curr += std::exp(-bound); + denominator_prop += MY_EXP(-bound); + denominator_curr += MY_EXP(-bound); for (int cat = 0; cat < n_cats; cat++) { int score = cat + 1; double exponent = GroupThresholds[cat] + rest_score * score - bound; - denominator_curr += std::exp(exponent + score * obs_current_p); - denominator_prop += std::exp(exponent + score * obs_proposed_p); + denominator_curr += MY_EXP(exponent + score * obs_current_p); + denominator_prop += MY_EXP(exponent + score * obs_proposed_p); } } else { for (int cat = 0; cat <= n_cats; cat++) { @@ -835,13 +836,13 @@ double log_pseudolikelihood_ratio_pairwise_differences( GroupThresholds[1] * (cat - baseline_category[var]) * (cat - baseline_category[var]) + rest_score * cat - bound; - denominator_curr += std::exp(exponent + cat * obs_current_p); - denominator_prop += std::exp(exponent + cat * obs_proposed_p); + denominator_curr += MY_EXP(exponent + cat * obs_current_p); + denominator_prop += MY_EXP(exponent + cat * obs_proposed_p); } } // Update log pseudo-likelihood ratio - pseudolikelihood_ratio += std::log(denominator_curr) - std::log(denominator_prop); + pseudolikelihood_ratio += MY_LOG(denominator_curr) - MY_LOG(denominator_prop); } } } @@ -940,11 +941,11 @@ void metropolis_pairwise_difference_between_model( // Update log probability for the inclusion inclusion_indicator if (inclusion_indicator(variable1, variable2) == 1) { - log_acceptance_probability -= std::log(inclusion_probability_difference(variable1, variable2)); - log_acceptance_probability += std::log(1 - inclusion_probability_difference(variable1, variable2)); + log_acceptance_probability -= MY_LOG(inclusion_probability_difference(variable1, variable2)); + log_acceptance_probability += MY_LOG(1 - inclusion_probability_difference(variable1, variable2)); } else { - log_acceptance_probability += std::log(inclusion_probability_difference(variable1, variable2)); - log_acceptance_probability -= std::log(1 - inclusion_probability_difference(variable1, variable2)); + log_acceptance_probability += MY_LOG(inclusion_probability_difference(variable1, variable2)); + log_acceptance_probability -= MY_LOG(1 - inclusion_probability_difference(variable1, variable2)); } // Compute log pseudo-likelihood ratio @@ -956,7 +957,7 @@ void metropolis_pairwise_difference_between_model( // Metropolis-Hastings acceptance step double U = R::unif_rand(); - if (std::log(U) < log_acceptance_probability) { + if (MY_LOG(U) < log_acceptance_probability) { // Update inclusion inclusion_indicator inclusion_indicator(variable1, variable2) = 1 - inclusion_indicator(variable1, variable2); inclusion_indicator(variable2, variable1) = inclusion_indicator(variable1, variable2); @@ -1036,7 +1037,7 @@ void metropolis_threshold_regular( // Iterate over each category threshold for the variable for(int category = 0; category < n_cats; category++) { double current_state = main_effects(cat_index + category, 0); // Current state of the threshold - double exp_current = std::exp(current_state); // Exponentiated current threshold + double exp_current = MY_EXP(current_state); // Exponentiated current threshold double c = (prior_threshold_alpha + prior_threshold_beta) / (1 + exp_current); // Initial value for c // Compute group-specific thresholds and contributions to `q` and `r` @@ -1060,11 +1061,11 @@ void metropolis_threshold_regular( for (int cat = 0; cat < n_cats; ++cat) { if (cat != category) { double exponent = GroupThresholds[cat] + (cat + 1) * rest_score; - q_person += std::exp(exponent); + q_person += MY_EXP(exponent); } } double exponent_r = GroupThresholds[category] + (category + 1) * rest_score; - double r_person = std::exp(exponent_r); + double r_person = MY_EXP(exponent_r); q[person] = q_person; r[person] = r_person; c += r_person / (q_person + r_person * exp_current); @@ -1086,8 +1087,8 @@ void metropolis_threshold_regular( // Sample from the beta distribution double tmp_beta = R::rbeta(a, b); - double proposed_state = std::log(tmp_beta / (1 - tmp_beta) / c); - double exp_proposed = std::exp(proposed_state); + double proposed_state = MY_LOG(tmp_beta / (1 - tmp_beta) / c); + double exp_proposed = MY_EXP(proposed_state); // Compute the log acceptance probability double log_acceptance_probability = 0.0; @@ -1095,21 +1096,21 @@ void metropolis_threshold_regular( // Compute pseudo-likelihood ratio for (int gr = 0; gr < num_groups; ++gr) { for (int person = group_indices(gr, 0); person <= group_indices(gr, 1); ++person) { - log_acceptance_probability += std::log(q[person] + r[person] * exp_current); - log_acceptance_probability -= std::log(q[person] + r[person] * exp_proposed); + log_acceptance_probability += MY_LOG(q[person] + r[person] * exp_current); + log_acceptance_probability -= MY_LOG(q[person] + r[person] * exp_proposed); } } // Add prior density ratio - log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + exp_proposed); - log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + exp_current); + log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + exp_proposed); + log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + exp_current); // Add proposal density ratio - log_acceptance_probability -= (a + b) * std::log(1 + c * exp_current); - log_acceptance_probability += (a + b) * std::log(1 + c * exp_proposed); + log_acceptance_probability -= (a + b) * MY_LOG(1 + c * exp_current); + log_acceptance_probability += (a + b) * MY_LOG(1 + c * exp_proposed); // Perform Metropolis-Hastings acceptance step - double U = std::log(R::unif_rand()); + double U = MY_LOG(R::unif_rand()); if(U < log_acceptance_probability) { main_effects(cat_index + category, 0) = proposed_state; // Accept the proposal } @@ -1195,17 +1196,17 @@ double log_pseudolikelihood_ratio_main_difference( double bound = (rest_score > 0) ? n_cats * rest_score : 0.0; // Compute the denominators for proposed and current thresholds - double denominator_proposed = std::exp(-bound); - double denominator_current = std::exp(-bound); + double denominator_proposed = MY_EXP(-bound); + double denominator_current = MY_EXP(-bound); for(int cat = 0; cat < n_cats; cat++) { double exponent = (cat + 1) * rest_score - bound; - denominator_proposed += std::exp(exponent + proposed_thresholds[cat]); - denominator_current += std::exp(exponent + current_thresholds[cat]); + denominator_proposed += MY_EXP(exponent + proposed_thresholds[cat]); + denominator_current += MY_EXP(exponent + current_thresholds[cat]); } // Update the pseudo-likelihood ratio with log-likelihood differences - pseudolikelihood_ratio -= std::log(denominator_proposed); - pseudolikelihood_ratio += std::log(denominator_current); + pseudolikelihood_ratio -= MY_LOG(denominator_proposed); + pseudolikelihood_ratio += MY_LOG(denominator_current); } } @@ -1293,7 +1294,7 @@ void metropolis_main_difference_regular( // Metropolis-Hastings acceptance step double U = R::unif_rand(); - if (std::log(U) < log_acceptance_probability) { + if (MY_LOG(U) < log_acceptance_probability) { main_effects(cat_index, h) = proposed_state; // Accept the proposed state } @@ -1418,11 +1419,11 @@ double log_pseudolikelihood_ratio_thresholds_blumecapel( denominator = 0.0; for(int category = 0; category <= num_categories[variable]; category ++) { exponent = category * rest_score - bound; - numerator += std::exp(constant_numerator[category] + exponent); - denominator += std::exp(constant_denominator[category] + exponent); + numerator += MY_EXP(constant_numerator[category] + exponent); + denominator += MY_EXP(constant_denominator[category] + exponent); } - pseudolikelihood_ratio += std::log(numerator); - pseudolikelihood_ratio -= std::log(denominator); + pseudolikelihood_ratio += MY_LOG(numerator); + pseudolikelihood_ratio -= MY_LOG(denominator); } } @@ -1498,12 +1499,12 @@ void metropolis_threshold_blumecapel( // Add prior ratio to log acceptance probability log_acceptance_probability += prior_threshold_alpha * (proposed_state - current_state); - log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(current_state)); - log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(proposed_state)); + log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(current_state)); + log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(proposed_state)); // Perform Metropolis-Hastings acceptance step double U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { main_effects(cat_index, 0) = proposed_state; } @@ -1530,12 +1531,12 @@ void metropolis_threshold_blumecapel( // Add prior ratio to log acceptance probability log_acceptance_probability += prior_threshold_alpha * (proposed_state - current_state); - log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(current_state)); - log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(proposed_state)); + log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(current_state)); + log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(proposed_state)); // Perform Metropolis-Hastings acceptance step U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { main_effects(cat_index + 1, 0) = proposed_state; } @@ -1672,13 +1673,13 @@ double log_pseudolikelihood_ratio_main_difference_blumecapel( // Compute category-specific contributions for(int category = 0; category <= num_categories[variable]; category ++) { exponent = category * rest_score - bound; - numerator += std::exp(constant_numerator[category] + exponent); - denominator += std::exp(constant_denominator[category] + exponent); + numerator += MY_EXP(constant_numerator[category] + exponent); + denominator += MY_EXP(constant_denominator[category] + exponent); } // Update the pseudo-likelihood ratio with log probabilities - pseudolikelihood_ratio += std::log(numerator); - pseudolikelihood_ratio -= std::log(denominator); + pseudolikelihood_ratio += MY_LOG(numerator); + pseudolikelihood_ratio -= MY_LOG(denominator); } } @@ -1776,7 +1777,7 @@ void metropolis_main_difference_blumecapel( // Metropolis-Hastings acceptance step U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { main_effects(cat_index, h) = proposed_state; // Accept proposed value } @@ -1808,7 +1809,7 @@ void metropolis_main_difference_blumecapel( // Metropolis-Hastings acceptance step U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { main_effects(cat_index + 1, h) = proposed_state; } @@ -1885,7 +1886,7 @@ void metropolis_thresholds_regular_free( // Loop over categories for(int category = 0; category < n_cats; category++) { double current_state = main_effects(base_cat_index + category, group); - double exp_current = std::exp(current_state); + double exp_current = MY_EXP(current_state); // Initialize scaling factor `c` for the generalized beta-prime proposal double c = (prior_threshold_alpha + prior_threshold_beta) / (1 + exp_current); @@ -1898,12 +1899,12 @@ void metropolis_thresholds_regular_free( double q_person = 1.0; for(int cat = 0; cat < n_cats; cat++) { if(cat != category) { - q_person += std::exp(main_effects(base_cat_index + cat, group) + (cat + 1) * rest_score); + q_person += MY_EXP(main_effects(base_cat_index + cat, group) + (cat + 1) * rest_score); } } // Calculate pseudo-likelihood component `r` - double r_person = std::exp((category + 1) * rest_score); + double r_person = MY_EXP((category + 1) * rest_score); // Store results for each person q[person] = q_person; @@ -1920,8 +1921,8 @@ void metropolis_thresholds_regular_free( double a = num_obs_categories_gr(category, variable) + prior_threshold_alpha; double b = num_persons + prior_threshold_beta - num_obs_categories_gr(category, variable); double tmp = R::rbeta(a, b); - double proposed_state = std::log(tmp / ((1 - tmp) * c)); - double exp_proposed = std::exp(proposed_state); + double proposed_state = MY_LOG(tmp / ((1 - tmp) * c)); + double exp_proposed = MY_EXP(proposed_state); // Compute the log acceptance probability @@ -1929,20 +1930,20 @@ void metropolis_thresholds_regular_free( // Add pseudo-likelihood ratio contributions for (int person = 0; person < num_persons; ++person) { - log_acceptance_probability += std::log(q[person] + r[person] * exp_current); - log_acceptance_probability -= std::log(q[person] + r[person] * exp_proposed); + log_acceptance_probability += MY_LOG(q[person] + r[person] * exp_current); + log_acceptance_probability -= MY_LOG(q[person] + r[person] * exp_proposed); } // Add prior ratio contributions - log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + exp_proposed); - log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + exp_current); + log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + exp_proposed); + log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + exp_current); // Add proposal ratio contributions - log_acceptance_probability -= (a + b) * std::log(1 + c * exp_current); - log_acceptance_probability += (a + b) * std::log(1 + c * exp_proposed); + log_acceptance_probability -= (a + b) * MY_LOG(1 + c * exp_current); + log_acceptance_probability += (a + b) * MY_LOG(1 + c * exp_proposed); // Perform the Metropolis-Hastings acceptance step - double U = std::log(R::unif_rand()); + double U = MY_LOG(R::unif_rand()); if(U < log_acceptance_probability) { // Update the main effects matrix if the proposal is accepted main_effects(base_cat_index + category, group) = proposed_state; @@ -2061,24 +2062,24 @@ void metropolis_thresholds_blumecapel_free( } // Compute likelihood numerator and denominator - double numerator = std::exp(constant_numerator[0] - bound); - double denominator = std::exp(constant_denominator[0] - bound); + double numerator = MY_EXP(constant_numerator[0] - bound); + double denominator = MY_EXP(constant_denominator[0] - bound); for(int score = 1; score <= num_categories(variable, group); score++) { double exponent = score * rest_score - bound; - numerator += std::exp(constant_numerator[score] + exponent); - denominator += std::exp(constant_denominator[score] + exponent); + numerator += MY_EXP(constant_numerator[score] + exponent); + denominator += MY_EXP(constant_denominator[score] + exponent); } - log_acceptance_probability += std::log(numerator); - log_acceptance_probability -= std::log(denominator); + log_acceptance_probability += MY_LOG(numerator); + log_acceptance_probability -= MY_LOG(denominator); } // Add prior ratio to log acceptance probability - log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(current_state)); - log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(proposed_state)); + log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(current_state)); + log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(proposed_state)); // Metropolis acceptance step double U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { main_effects(cat_index, group) = proposed_state; } @@ -2126,23 +2127,23 @@ void metropolis_thresholds_blumecapel_free( bound += num_categories[variable] * rest_score; } - double numerator = std::exp(constant_numerator[0] - bound); - double denominator = std::exp(constant_denominator[0] - bound); + double numerator = MY_EXP(constant_numerator[0] - bound); + double denominator = MY_EXP(constant_denominator[0] - bound); for(int score = 1; score <= num_categories(variable, group); score ++) { double exponent = score * rest_score - bound; - numerator += std::exp(constant_numerator[score] + exponent); - denominator += std::exp(constant_denominator[score] + exponent); + numerator += MY_EXP(constant_numerator[score] + exponent); + denominator += MY_EXP(constant_denominator[score] + exponent); } - log_acceptance_probability += std::log(numerator); - log_acceptance_probability -= std::log(denominator); + log_acceptance_probability += MY_LOG(numerator); + log_acceptance_probability -= MY_LOG(denominator); } - log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(current_state)); - log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * std::log(1 + std::exp(proposed_state)); + log_acceptance_probability += (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(current_state)); + log_acceptance_probability -= (prior_threshold_alpha + prior_threshold_beta) * MY_LOG(1 + MY_EXP(proposed_state)); U = R::unif_rand(); - if(std::log(U) < log_acceptance_probability) { + if(MY_LOG(U) < log_acceptance_probability) { main_effects(cat_index + 1, group) = proposed_state; } @@ -2221,19 +2222,19 @@ double log_pseudolikelihood_ratio_main_difference_regular_between_model( double bound = (rest_score > 0) ? num_cats * rest_score : 0.0; // Initialize denominator terms for the proposed and current pseudolikelihoods - double denominator_proposed = std::exp(-bound); - double denominator_current = std::exp(-bound); + double denominator_proposed = MY_EXP(-bound); + double denominator_current = MY_EXP(-bound); // Add contributions from each category to the denominators for(int cat = 0; cat < num_cats; cat++) { double exponent = (cat + 1) * rest_score - bound; - denominator_proposed += std::exp(exponent + proposed_thresholds[cat]); - denominator_current += std::exp(exponent + current_thresholds[cat]); + denominator_proposed += MY_EXP(exponent + proposed_thresholds[cat]); + denominator_current += MY_EXP(exponent + current_thresholds[cat]); } // Update the pseudo-likelihood ratio with log-likelihood differences - pseudolikelihood_ratio -= std::log(denominator_proposed); - pseudolikelihood_ratio += std::log(denominator_current); + pseudolikelihood_ratio -= MY_LOG(denominator_proposed); + pseudolikelihood_ratio += MY_LOG(denominator_current); } } @@ -2337,16 +2338,16 @@ void metropolis_main_difference_regular_between_model( // Add prior inclusion odds contributions if(current_indicator == 1) { - log_prob -= std::log(inclusion_probability_difference(variable, variable)); - log_prob += std::log(1 - inclusion_probability_difference(variable, variable)); + log_prob -= MY_LOG(inclusion_probability_difference(variable, variable)); + log_prob += MY_LOG(1 - inclusion_probability_difference(variable, variable)); } else { - log_prob += std::log(inclusion_probability_difference(variable, variable)); - log_prob -= std::log(1 - inclusion_probability_difference(variable, variable)); + log_prob += MY_LOG(inclusion_probability_difference(variable, variable)); + log_prob -= MY_LOG(1 - inclusion_probability_difference(variable, variable)); } // Perform Metropolis-Hastings step double U = R::unif_rand(); - if(std::log(U) < log_prob) { + if(MY_LOG(U) < log_prob) { inclusion_indicator(variable, variable) = proposed_indicator; for(int cat = 0; cat < num_cats; cat++) { for(int h = 1; h < num_groups; h++) { @@ -2456,13 +2457,13 @@ double log_pseudolikelihood_ratio_main_difference_blume_capel_between_model( // Compute category-specific contributions for(int cat = 0; cat <= num_cats; cat++) { double exponent = cat * rest_score - bound; - denominator_proposed += std::exp(proposed_denominator_constant[cat] + exponent); - denominator_current += std::exp(current_denominator_constant[cat] + exponent); + denominator_proposed += MY_EXP(proposed_denominator_constant[cat] + exponent); + denominator_current += MY_EXP(current_denominator_constant[cat] + exponent); } // Update the pseudo-likelihood ratio - pseudolikelihood_ratio -= std::log(denominator_proposed); - pseudolikelihood_ratio += std::log(denominator_current); + pseudolikelihood_ratio -= MY_LOG(denominator_proposed); + pseudolikelihood_ratio += MY_LOG(denominator_current); } } @@ -2590,16 +2591,16 @@ void metropolis_main_difference_blume_capel_between_model( // Add prior odds contribution if(current_inclusion_indicator == 1) { - log_prob -= std::log(inclusion_probability_difference(variable, variable)); - log_prob += std::log(1-inclusion_probability_difference(variable, variable)); + log_prob -= MY_LOG(inclusion_probability_difference(variable, variable)); + log_prob += MY_LOG(1-inclusion_probability_difference(variable, variable)); } else { - log_prob += std::log(inclusion_probability_difference(variable, variable)); - log_prob -= std::log(1-inclusion_probability_difference(variable, variable)); + log_prob += MY_LOG(inclusion_probability_difference(variable, variable)); + log_prob -= MY_LOG(1-inclusion_probability_difference(variable, variable)); } // Perform Metropolis-Hastings acceptance step double U = R::unif_rand(); - if(std::log(U) < log_prob) { + if(MY_LOG(U) < log_prob) { inclusion_indicator(variable, variable) = proposed_inclusion_indicator; main_effects.rows(start_index, stop_index).cols(1, num_groups - 1) = proposed_main_effects.cols(1, num_groups - 1); } @@ -2701,7 +2702,7 @@ List gibbs_step_gm( const arma::imat& index ) { // Precompute the Robbins-Monro decay term - double exp_neg_log_t_rm_adaptation_rate = std::exp(-std::log(t) * rm_adaptation_rate); + double exp_neg_log_t_rm_adaptation_rate = MY_EXP(-MY_LOG(t) * rm_adaptation_rate); // Step 1: Update pairwise interaction parameters metropolis_interaction( diff --git a/src/gibbs_functions_edge_prior.cpp b/src/gibbs_functions_edge_prior.cpp index fef05628..236f8994 100644 --- a/src/gibbs_functions_edge_prior.cpp +++ b/src/gibbs_functions_edge_prior.cpp @@ -1,4 +1,5 @@ #include +#include "explog_switch.h" using namespace Rcpp; // ----------------------------------------------------------------------------| @@ -96,10 +97,10 @@ NumericVector compute_Vn_mfm_sbm(int no_variables, for(int k = t; k < 500; k++) { tmp = 0.0; for(int tt = 1 - t; tt < 1; tt ++) { - tmp += std::log(dirichlet_alpha * k + tt); + tmp += MY_LOG(dirichlet_alpha * k + tt); } for(int n = 0; n < no_variables; n++) { - tmp -= std::log(dirichlet_alpha * k + n); + tmp -= MY_LOG(dirichlet_alpha * k + n); } tmp -= std::lgamma(k + 1); @@ -109,12 +110,12 @@ NumericVector compute_Vn_mfm_sbm(int no_variables, // Compute the maximum between r and tmp if (tmp > r) { - r = std::log(std::exp(r - tmp) + 1) + tmp; + r = MY_LOG(MY_EXP(r - tmp) + 1) + tmp; } else { - r = std::log(1 + std::exp(tmp - r)) + r; + r = MY_LOG(1 + MY_EXP(tmp - r)) + r; } } - log_Vn(t-1) = r - std::log(std::exp(1) - 1); + log_Vn(t-1) = r - MY_LOG(MY_EXP(1) - 1); } return log_Vn; } @@ -133,14 +134,14 @@ double log_likelihood_mfm_sbm(IntegerVector cluster_assign, if(j != node) { if(j < node) { output += indicator(j, node) * - std::log(cluster_probs(cluster_assign[j], cluster_assign[node])); + MY_LOG(cluster_probs(cluster_assign[j], cluster_assign[node])); output += (1 - indicator(j, node)) * - std::log(1 - cluster_probs(cluster_assign[j], cluster_assign[node])); + MY_LOG(1 - cluster_probs(cluster_assign[j], cluster_assign[node])); } else { output += indicator(node, j) * - std::log(cluster_probs(cluster_assign[node], cluster_assign[j])); + MY_LOG(cluster_probs(cluster_assign[node], cluster_assign[j])); output += (1 - indicator(node, j)) * - std::log(1 - cluster_probs(cluster_assign[node], cluster_assign[j])); + MY_LOG(1 - cluster_probs(cluster_assign[node], cluster_assign[j])); } } } @@ -264,7 +265,7 @@ IntegerVector block_allocations_mfm_sbm(IntegerVector cluster_assign, no_variables); prob = (dirichlet_alpha + cluster_size_node[c]) * - std::exp(loglike); + MY_EXP(loglike); } else { logmarg = log_marginal_mfm_sbm(cluster_assign_tmp, indicator, @@ -274,8 +275,8 @@ IntegerVector block_allocations_mfm_sbm(IntegerVector cluster_assign, beta_bernoulli_beta); prob = dirichlet_alpha * - std::exp(logmarg) * - std::exp(log_Vn[no_clusters - 1] - log_Vn[no_clusters - 2]); + MY_EXP(logmarg) * + MY_EXP(log_Vn[no_clusters - 1] - log_Vn[no_clusters - 2]); } cluster_prob[c] = prob; @@ -316,7 +317,7 @@ IntegerVector block_allocations_mfm_sbm(IntegerVector cluster_assign, no_variables); prob = (dirichlet_alpha + cluster_size_node[c]) * - std::exp(loglike); + MY_EXP(loglike); } else { logmarg = log_marginal_mfm_sbm(cluster_assign_tmp, indicator, @@ -326,8 +327,8 @@ IntegerVector block_allocations_mfm_sbm(IntegerVector cluster_assign, beta_bernoulli_beta); prob = dirichlet_alpha * - std::exp(logmarg) * - std::exp(log_Vn[no_clusters] - log_Vn[no_clusters - 1]); + MY_EXP(logmarg) * + MY_EXP(log_Vn[no_clusters] - log_Vn[no_clusters - 1]); } cluster_prob[c] = prob; From adb3c76672d638b0ff035d90687b09b650cc0c3a Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Tue, 17 Jun 2025 13:39:02 +0200 Subject: [PATCH 08/10] cleanup --- src/Makevars.win | 6 +---- validate_exp_and_log.R | 57 ------------------------------------------ 2 files changed, 1 insertion(+), 62 deletions(-) delete mode 100644 validate_exp_and_log.R diff --git a/src/Makevars.win b/src/Makevars.win index 3c8dc8e2..a00f5b8f 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,5 +1 @@ -#DEBUG=1 -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -#PKG_CXXFLAGS+= -g -#PKG_DLLFLAGS= -#PKG_SHLIB_LDFLAGS= \ No newline at end of file +PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) \ No newline at end of file diff --git a/validate_exp_and_log.R b/validate_exp_and_log.R deleted file mode 100644 index 2cc7043e..00000000 --- a/validate_exp_and_log.R +++ /dev/null @@ -1,57 +0,0 @@ - -values <- matrix(c( - -1e6, 1e6, - -1e4, 1e4, - -1e3, 1e3, - -1e2, 1e2, - -1e1, 1e1, - -1e-3, 1e-3 -), ncol = 2, byrow = TRUE) -size <- 1E6L - -for (i in seq_len(nrow(values))) { - from <- values[i, 1] - to <- values[i, 2] - cat("i = ", i, ", Testing range: ", from, " to ", to, "\n") - test_values <- seq(from, to, length.out = size) - y_R <- exp(test_values) - y_R2 <- bgms:::rcpp_ieee754_exp(test_values) - result <- all.equal(y_R, y_R2) - - abs_error = max(abs(y_R - y_R2), na.rm = TRUE) - rel_error = max(abs((y_R - y_R2) / y_R), na.rm = TRUE) - - cat("i = ", i, ", All.equal = ", isTRUE(result), - " max abs error:", sprintf("%.8g", abs_error), - " max rel error:", sprintf("%.8g", rel_error), - "\n") -} - -values <- matrix(c( - 1e10, 1e30, - 1e04, 1e09, - 1e01, 1e03, - 1e-2, 1e01, - 1e-5, 1e-3, - 1e-20, 1e-6 -), ncol = 2, byrow = TRUE) -size <- 1E6L - -for (i in seq_len(nrow(values))) { - from <- values[i, 1] - to <- values[i, 2] - cat("i = ", i, ", Testing range: ", from, " to ", to, "\n") - test_values <- seq(from, to, length.out = size) - y_R <- log(test_values) - y_R2 <- bgms:::rcpp_ieee754_log(test_values) - result <- all.equal(y_R, y_R2) - - abs_error = max(abs(y_R - y_R2), na.rm = TRUE) - rel_error = max(abs((y_R - y_R2) / y_R), na.rm = TRUE) - - cat("i = ", i, ", All.equal = ", isTRUE(result), - " max abs error:", sprintf("%.8g", abs_error), - " max rel error:", sprintf("%.8g", rel_error), - "\n") -} - From d840ac0ac8efaa89012897656f41f32dcbcfd126 Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Wed, 9 Jul 2025 10:39:05 +0200 Subject: [PATCH 09/10] remove duplicated function --- src/gibbs_functions.cpp | 109 +--------------------------------------- 1 file changed, 1 insertion(+), 108 deletions(-) diff --git a/src/gibbs_functions.cpp b/src/gibbs_functions.cpp index 4a4e238b..d85173fb 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -385,114 +385,7 @@ void metropolis_thresholds_blumecapel( } -double log_pseudolikelihood_ratio( - const NumericMatrix& interactions, - const NumericMatrix& thresholds, - const IntegerMatrix& observations, - const IntegerVector& no_categories, - const int no_persons, - const int variable1, - const int variable2, - const double proposed_state, - const double current_state, - const NumericMatrix& rest_matrix, - const LogicalVector& variable_bool, - const IntegerVector& reference_category -) { - double rest_score, bound; - double pseudolikelihood_ratio = 0.0; - double denominator_prop, denominator_curr, exponent; - int score, obs_score1, obs_score2; - - double delta_state = proposed_state - current_state; - - int nc1 = no_categories[variable1]; - int nc2 = no_categories[variable2]; - bool is_ord1 = variable_bool[variable1]; - bool is_ord2 = variable_bool[variable2]; - int ref_cat1 = reference_category[variable1]; - int ref_cat2 = reference_category[variable2]; - - for (int person = 0; person < no_persons; person++) { - obs_score1 = observations(person, variable1); - obs_score2 = observations(person, variable2); - pseudolikelihood_ratio += 2 * obs_score1 * obs_score2 * delta_state; - - // Variable 1 - rest_score = rest_matrix(person, variable1) - - obs_score2 * interactions(variable2, variable1); - bound = rest_score > 0.0 ? nc1 * rest_score : 0.0; - - double obs2_prop_state = obs_score2 * proposed_state; - double obs2_curr_state = obs_score2 * current_state; - - if (is_ord1) { - denominator_prop = MY_EXP(-bound); - denominator_curr = MY_EXP(-bound); - int index = variable1; - for (int category = 0; category < nc1; category++) { - score = category + 1; - exponent = thresholds[index] + score * rest_score - bound; - index += thresholds.nrow(); - double exp_term = exponent; - denominator_prop += MY_EXP(exp_term + score * obs2_prop_state); - denominator_curr += MY_EXP(exp_term + score * obs2_curr_state); - } - } else { - denominator_prop = 0.0; - denominator_curr = 0.0; - for (int category = 0; category <= nc1; category++) { - int diff = category - ref_cat1; - exponent = thresholds(variable1, 0) * category + - thresholds(variable1, 1) * diff * diff + - category * rest_score - bound; - denominator_prop += MY_EXP(exponent + category * obs2_prop_state); - denominator_curr += MY_EXP(exponent + category * obs2_curr_state); - } - } - - pseudolikelihood_ratio += MY_LOG(denominator_curr / denominator_prop); - - // Variable 2 - rest_score = rest_matrix(person, variable2) - - obs_score1 * interactions(variable1, variable2); - bound = rest_score > 0.0 ? nc2 * rest_score : 0.0; - - double obs1_prop_state = obs_score1 * proposed_state; - double obs1_curr_state = obs_score1 * current_state; - - if (is_ord2) { - denominator_prop = MY_EXP(-bound); - denominator_curr = MY_EXP(-bound); - int index = variable2; - for (int category = 0; category < nc2; category++) { - score = category + 1; - exponent = thresholds[index] + score * rest_score - bound; - index += thresholds.nrow(); - denominator_prop += MY_EXP(exponent + score * obs1_prop_state); - denominator_curr += MY_EXP(exponent + score * obs1_curr_state); - } - } else { - denominator_prop = 0.0; - denominator_curr = 0.0; - for (int category = 0; category <= nc2; category++) { - int diff = category - ref_cat2; - exponent = thresholds(variable2, 0) * category + - thresholds(variable2, 1) * diff * diff + - category * rest_score - bound; - denominator_prop += MY_EXP(exponent + category * obs1_prop_state); - denominator_curr += MY_EXP(exponent + category * obs1_curr_state); - } - } - - pseudolikelihood_ratio += MY_LOG(denominator_curr / denominator_prop); - } - - return pseudolikelihood_ratio; -} - -/* // ----------------------------------------------------------------------------| // The log pseudolikelihood ratio [proposed against current] for an interaction // ----------------------------------------------------------------------------| @@ -612,7 +505,7 @@ double log_pseudolikelihood_ratio( } return pseudolikelihood_ratio; } -*/ + // ----------------------------------------------------------------------------| // MH algorithm to sample from the full-conditional of the active interaction // parameters for Bayesian edge selection From c249e46a49cd8c390a5b117d874939f1dcd736ba Mon Sep 17 00:00:00 2001 From: Don van den Bergh Date: Wed, 9 Jul 2025 10:41:49 +0200 Subject: [PATCH 10/10] remove newlines --- src/gibbs_functions.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/gibbs_functions.cpp b/src/gibbs_functions.cpp index d85173fb..3761c863 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -591,7 +591,6 @@ void metropolis_interactions( } } - // ----------------------------------------------------------------------------| // MH algorithm to sample from the full-conditional of an edge + interaction // pair for Bayesian edge selection @@ -684,7 +683,6 @@ void metropolis_edge_interaction_pair( } } - // ----------------------------------------------------------------------------| // A Gibbs step for graphical model parameters for Bayesian edge selection // ----------------------------------------------------------------------------|