在rcpp犰狳中将向量与double进行比较

Posted 2023-02-16

【中文标题】在rcpp犰狳中将向量与double进行比较

【英文标题】：Compare vector to double in rcpp armadillo

【发布时间】：2021-01-07 16:07:36

【问题描述】：

考虑这个 Rcpp Armadillo 函数：

// [[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;

现在运行以下 R 脚本我得到不一致的结果：

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

所以不知何故，比较失败了。对于如何解决这个问题的任何建议，我都会很高兴。来自 R，我认为循环对于这项任务来说太复杂了。也许我错了。

【问题讨论】：

在 C++ 中不会自动回收。将标量 double 展开为向量，然后进行比较——Rcpp Sugar 涵盖了这一点。

对不起：“计算”不是“比较”。我们有矢量化操作，但你必须准备输入。

【参考方案1】：

几个简短的陈述：

Armadillo C++ 矩阵线性代数库在向量到标量运算之外没有自动回收功能。 C++ 确实有更严格的类型控制，从导致 unsigned int 的比较到使用 double 的减法都有问题。

第二部分确实是麻烦所在。为了说明，我们将通过编写如下调试函数系统地逐步执行每个操作：

// [[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);

运行函数给出：

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

输出与预期相反：

(1.0*(y < x) - z)

通过进行显式类型转换，从uvec 到vec，计算再次可行：

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

注意，显式转换请求是通过arma::conv_to<>::from()对计算的向量部分完成的。