Compare vector to double in rcpp armadillo

Question

Consider this Rcpp Armadillo function:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;

// [[Rcpp::export]]
vec testfun(const vec &x,
            const double &y,
            const double &z)
{
    vec out = ((y < x) - z) % (x - y);
    return out;
}

Now running the following R Script I get inconsistent results:

Rcpp::sourceCpp("functions/test.cpp")

x <- 1:3
y <- 2
z <- 0.5

out_r <- ((y < x) - z) * (x - y)
out_cpp <- testfun(x, y, z)

print(out_r)
print(out_cpp)
[1] 0.5 0.0 0.5
     [,1]
[1,]    0
[2,]    0
[3,]    1

So somehow the comparison fails. I would be glad for any advice on how to solve this. Coming from R, I think a loop is too complicated for this task. Maybe I'm wrong.

Answer 1

A couple of quick statements:

1. The Armadillo C++ Matrix Linear Algebra Library does not have automatic recycling outside of vector-to-scalar operations.
2. C++ does have stricter type controls that make going from a comparison resulting in unsigned int to subtraction with a double problematic.

The second part is really where trouble resides. To illustrate, we'll systematically step through each operation by writing a debug function as follows:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;

// [[Rcpp::export]]
arma::vec debug_inline_statements(const vec &x, const double &y, const double &z)
{
    // Isolate the problem statement:
    // Ok
    Rcpp::Rcout << "(x - y)" << std::endl << (x - y) << std::endl;
    // Ok
    Rcpp::Rcout << "1.0*(y < x):" << std::endl << 1.0*(y < x) << std::endl;
    // Bad
    Rcpp::Rcout << "(1.0*(y < x) - z):" << std::endl << ((1.0*(y < x)) - z) << std::endl;
    // What went wrong? Conversion from Unsigned integer to double. 
    
    // Solution: Help the template expansion:
    vec bool_to_double = arma::conv_to<arma::vec>::from(y < x);
    Rcpp::Rcout << "(double(y < x) - z):" << std::endl << (bool_to_double - z) << std::endl;
    // Success!
    
    // All together now:
    Rcpp::Rcout << "(double(y < x) - z) % (x - y):" << std::endl << 
      (arma::conv_to<arma::vec>::from(y < x) - z) % (x - y) << std::endl;
    
    return (arma::conv_to<arma::vec>::from(y < x) - z) % (x - y);
}

Running the function gives:

x <- 1:3
y <- 2
z <- 0.5

out_cpp <- debug_inline_statements(x, y, z)
# (x - y)
#   -1.0000
#         0
#    1.0000
# 
# 1.0*(y < x):
#         0
#         0
#         1
# 
# (1.0*(y < x) - z):
#         0
#         0
#         1
# 
# (double(y < x) - z):
#   -0.5000
#   -0.5000
#    0.5000
# 
# (double(y < x) - z) % (x - y):
#    0.5000
#         0
#    0.5000

The output is against expectations at:

(1.0*(y < x) - z)

By making an explicit type conversion, from uvec to vec, the computation is viable again:

(arma::conv_to<arma::vec>::from(y < x) - z)

Note, the explicit conversion request was done on the vector portion of the computation through arma::conv_to<>::from().