Skip to content

Remove uses of arma::schur, that have a different semantics in Armadillo 6.100 #15

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Remove uses of arma::schur, that have a different semantics in Armadillo 6.100 #15

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 0 additions & 51 deletions cppbugs/distributions/mcmc.exponential.censored.hpp
View file
Open in desktop

This file was deleted.

88 changes: 0 additions & 88 deletions cppbugs/mcmc.arma.extensions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,94 +121,6 @@ namespace arma {
return eOpCube<T1, eop_factln>(A.get_ref());
}

// cube
//! element-wise multiplication of BaseCube objects with same element type
template<typename T1, typename T2>
arma_inline
const eGlueCube<T1, T2, eglue_schur>
schur
(
const BaseCube<typename T1::elem_type,T1>& X,
const BaseCube<typename T1::elem_type,T2>& Y
)
{
arma_extra_debug_sigprint();
return eGlueCube<T1, T2, eglue_schur>(X.get_ref(), Y.get_ref());
}

//! element-wise multiplication of BaseCube objects with different element types
template<typename T1, typename T2>
inline
const mtGlueCube<typename promote_type<typename T1::elem_type, typename T2::elem_type>::result, T1, T2, glue_mixed_schur>
schur
(
const BaseCube< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T1_result, T1>& X,
const BaseCube< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T2_result, T2>& Y
)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT1;
typedef typename T2::elem_type eT2;
typedef typename promote_type<eT1,eT2>::result out_eT;
promote_type<eT1,eT2>::check();
return mtGlueCube<out_eT, T1, T2, glue_mixed_schur>( X.get_ref(), Y.get_ref() );
}

// matrix
template<typename T1, typename T2>
arma_inline
const eGlue<T1, T2, eglue_schur>
schur(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y) {
arma_extra_debug_sigprint();
return eGlue<T1, T2, eglue_schur>(X.get_ref(), Y.get_ref());
}

//! element-wise multiplication of Base objects with different element types
template<typename T1, typename T2>
inline
const mtGlue<typename promote_type<typename T1::elem_type, typename T2::elem_type>::result, T1, T2, glue_mixed_schur>
schur
(
const Base< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T1_result, T1>& X,
const Base< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T2_result, T2>& Y
)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT1;
typedef typename T2::elem_type eT2;
typedef typename promote_type<eT1,eT2>::result out_eT;
promote_type<eT1,eT2>::check();
return mtGlue<out_eT, T1, T2, glue_mixed_schur>( X.get_ref(), Y.get_ref() );
}


//! Base * scalar
template<typename T1>
arma_inline
const eOp<T1, eop_scalar_times>
schur
(const Base<typename T1::elem_type,T1>& X, const typename T1::elem_type k)
{
arma_extra_debug_sigprint();
return eOp<T1, eop_scalar_times>(X.get_ref(),k);
}

//! scalar * Base
template<typename T1>
arma_inline
const eOp<T1, eop_scalar_times>
schur
(const typename T1::elem_type k, const Base<typename T1::elem_type,T1>& X)
{
arma_extra_debug_sigprint();
return eOp<T1, eop_scalar_times>(X.get_ref(),k); // NOTE: order is swapped
}

double schur(const int x, const double y) { return x * y; }
double schur(const double x, const int y) { return x * y; }
double schur(const double& x, const double& y) { return x * y; }
double schur(const int& x, const int& y) { return x * y; }

// insert an 'any' function for bools into the arma namespace
bool any(const bool x) {
return x;
Expand Down
43 changes: 37 additions & 6 deletions cppbugs/mcmc.math.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,37 @@
// Stochastic/Math related functions
namespace cppbugs {

template<typename T, typename U>
inline
auto schur_product(T&& x, U&& y) ->
decltype (arma::operator%(std::forward<T>(x), std::forward<U>(y))){
return arma::operator%(std::forward<T>(x), std::forward<U>(y));
}

template<typename T>
inline
auto schur_product(const typename T::elem_type& x, T&& y) ->
decltype (arma::operator*(x, std::forward<T>(y))){
return arma::operator*(x, std::forward<T>(y));
}

template<typename T>
inline
auto schur_product(T&& x, const typename T::elem_type& y) ->
decltype (arma::operator*(std::forward<T>(x), y)){
return arma::operator*(std::forward<T>(x), y);
}

double schur_product(const int x, const double y) { return x * y; }
double schur_product(const double x, const int y) { return x * y; }
double schur_product(const double x, const double y) { return x * y; }
double schur_product(const float x, const double y) { return x * y; }
double schur_product(const double x, const float y) { return x * y; }
float schur_product(const int x, const float y) { return x * y; }
float schur_product(const float x, const int y) { return x * y; }
float schur_product(const float x, const float y) { return x * y; }
int schur_product(const int x, const int y) { return x * y; }

static inline double square(double x) {
return x*x;
}
Expand Down Expand Up @@ -55,7 +86,7 @@ namespace cppbugs {

template<typename T, typename U, typename V>
double normal_logp(const T& x, const U& mu, const V& tau) {
return arma::accu(0.5*log_approx(0.5*tau/arma::math::pi()) - 0.5 * arma::schur(tau, square(x - mu)));
return arma::accu(0.5*log_approx(0.5*tau/arma::math::pi()) - 0.5 * schur_product(tau, square(x - mu)));
}

template<typename T, typename U, typename V>
Expand All @@ -67,15 +98,15 @@ namespace cppbugs {
double gamma_logp(const T& x, const U& alpha, const V& beta) {
return arma::any(arma::vectorise(x < 0)) ?
-std::numeric_limits<double>::infinity() :
arma::accu(arma::schur((alpha - 1.0),log_approx(x)) - arma::schur(beta,x) - lgamma(alpha) + arma::schur(alpha,log_approx(beta)));
arma::accu(schur_product((alpha - 1.0),log_approx(x)) - schur_product(beta,x) - lgamma(alpha) + schur_product(alpha,log_approx(beta)));
}

template<typename T, typename U, typename V>
double beta_logp(const T& x, const U& alpha, const V& beta) {
const double one = 1.0;
return arma::any(arma::vectorise(x <= 0)) || arma::any(arma::vectorise(x >= 1)) || arma::any(arma::vectorise(alpha <= 0)) || arma::any(arma::vectorise(beta <= 0)) ?
-std::numeric_limits<double>::infinity() :
arma::accu(lgamma(alpha+beta) - lgamma(alpha) - lgamma(beta) + arma::schur((alpha-one),log_approx(x)) + arma::schur((beta-one),log_approx(one-x)));
arma::accu(lgamma(alpha+beta) - lgamma(alpha) - lgamma(beta) + schur_product((alpha-one),log_approx(x)) + schur_product((beta-one),log_approx(one-x)));
}

double categorical_logp(const arma::ivec& x, const arma::mat& p) {
Expand Down Expand Up @@ -111,15 +142,15 @@ namespace cppbugs {
if(arma::any(arma::vectorise(p <= 0)) || arma::any(arma::vectorise(p >= 1)) || arma::any(arma::vectorise(x < 0)) || arma::any(arma::vectorise(x > n))) {
return -std::numeric_limits<double>::infinity();
}
return arma::accu(arma::schur(x,log_approx(p)) + arma::schur((n-x),log_approx(1-p)) + arma::factln(n) - arma::factln(x) - arma::factln(n-x));
return arma::accu(schur_product(x,log_approx(p)) + schur_product((n-x),log_approx(1-p)) + arma::factln(n) - arma::factln(x) - arma::factln(n-x));
}

template<typename T, typename U>
double bernoulli_logp(const T& x, const U& p) {
if( arma::any(arma::vectorise(p <= 0)) || arma::any(arma::vectorise(p >= 1)) || arma::any(arma::vectorise(x < 0)) || arma::any(arma::vectorise(x > 1)) ) {
return -std::numeric_limits<double>::infinity();
} else {
return arma::accu(arma::schur(x,log_approx(p)) + arma::schur((1-x), log_approx(1-p)));
return arma::accu(schur_product(x,log_approx(p)) + schur_product((1-x), log_approx(1-p)));
}
}

Expand All @@ -134,7 +165,7 @@ namespace cppbugs {

template<typename T, typename U>
double exponential_logp(const T& x, const U& lambda) {
return arma::accu(log_approx(lambda) - arma::schur(lambda, x));
return arma::accu(log_approx(lambda) - schur_product(lambda, x));
}

template<typename T, typename U>
Expand Down