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/R/RcppExports.R b/R/RcppExports.R index 0d64db97..35aa5389 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,6 +1,18 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +get_explog_switch <- function() { + .Call(`_bgms_get_explog_switch`) +} + +rcpp_ieee754_exp <- function(x) { + .Call(`_bgms_rcpp_ieee754_exp`, x) +} + +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) } 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/Makevars.win b/src/Makevars.win index 22ebc632..a00f5b8f 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1 +1 @@ -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) \ No newline at end of file diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3b71f51e..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) { @@ -143,6 +175,9 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { + {"_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}, 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 new file mode 100644 index 00000000..c39172ac --- /dev/null +++ b/src/e_exp.cpp @@ -0,0 +1,386 @@ + +/* @(#)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" + +// EDIT: fpmath.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. + */ + +// the defines below are from "math_private.h" +#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) + +/* 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,}, +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 + +// 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); + } +} diff --git a/src/e_exp.h b/src/e_exp.h new file mode 100644 index 00000000..b26f5ae9 --- /dev/null +++ b/src/e_exp.h @@ -0,0 +1,7 @@ +#ifndef _E_EXP_H_ +#define _E_EXP_H_ + +double __ieee754_exp(double x); +double __ieee754_log(double x); + +#endif 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 a5d8860b..3761c863 100644 --- a/src/gibbs_functions.cpp +++ b/src/gibbs_functions.cpp @@ -5,6 +5,8 @@ #include using namespace Rcpp; +#include "explog_switch.h" + // ----------------------------------------------------------------------------| // Impute missing data from full-conditional // ----------------------------------------------------------------------------| @@ -50,7 +52,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 +62,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 +70,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 +147,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 +155,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 +169,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); - exp_proposed = exp(proposed_state); + proposed_state = MY_LOG(tmp / (1 - tmp) / c); + 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 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 +264,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 +289,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(-MY_LOG(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -342,25 +344,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 +370,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(-MY_LOG(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -383,6 +385,7 @@ void metropolis_thresholds_blumecapel( } + // ----------------------------------------------------------------------------| // The log pseudolikelihood ratio [proposed against current] for an interaction // ----------------------------------------------------------------------------| @@ -425,17 +428,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 +451,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 +471,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,13 +495,13 @@ 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; } @@ -552,7 +556,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 +573,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(-MY_LOG(t) * phi); if(std::isnan(update_proposal_sd) == true) { update_proposal_sd = 1.0; @@ -647,7 +651,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, @@ -655,11 +659,11 @@ 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(); - 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); 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;