RcppArmadillo: arma::cube 的向量

Posted 2023-02-16

【中文标题】RcppArmadillo: arma::cube 的向量

【英文标题】：RcppArmadillo: vector of arma::cube's

【发布时间】：2018-04-19 16:05:18

【问题描述】：

我正在使用 RcppArmadillo 在 R / Rcpp 中开发应用程序，并且我需要使用 arma::cube 对象的向量。以下示例运行良好。

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace std;
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
bool posterior(int n, int L, NumericVector N, int TIME) 
    vector<cube> A(TIME);
    for (int t = 0; t < TIME; t++) A[t] = cube(n, L, max(N), fill::zeros);

    Rprintf("*** %.2f ***\n", A[3].at(5, 1, 2));

    return true;

这里有一些 R 代码来测试它。

library(Rcpp)

sourceCpp("VectorOfCubes.cpp")

posterior(200, 3, c(10, 5, 2), 10^4)

我的问题是：可以去掉上面C++函数中的for，直接初始化向量A吗？

【参考方案1】：

为什么不使用默认的vector 构造函数？

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
bool posterior_default(int n, int L, Rcpp::NumericVector N, int TIME) 
  std::vector<arma::cube> A(TIME, arma::cube(n, L, max(N), arma::fill::zeros));

  return true;

您可以使用std::fill 算法，c.f.

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
bool posterior(int n, int L, Rcpp::NumericVector N, int TIME) 
  std::vector<arma::cube> A(TIME);

  std::fill(A.begin(), A.end(), arma::cube(n, L, max(N), arma::fill::zeros));

  return true;

不过，一个更简单的变体是使用arma::field 来存储值而不是std::vector，然后使用.fill() 成员函数。

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::field<arma::cube> posterior_field(int n, int L, Rcpp::NumericVector N, int TIME) 
  arma::field<arma::cube> A(TIME);

  A.fill(arma::cube(n, L, max(N), arma::fill::zeros));

  return A;

输出：

posterior_field(3, 4, c(1,2,3), 10)
#      [,1]      
# [1,] Numeric,36
# [2,] Numeric,36
# [3,] Numeric,36
# [4,] Numeric,36
# [5,] Numeric,36
# [6,] Numeric,36
# [7,] Numeric,36
# [8,] Numeric,36
# [9,] Numeric,36
# [10,] Numeric,36

【讨论】：

很好的答案！非常感谢！

还有一件事。使用您的第二个建议arma::field<arma::cube> A(TIME); A.fill(arma::cube(n, L, max(N), arma::fill::zeros));，我以A[1].at(5, 1, 2) 访问这些元素。但是at() 不会检查边界，对吗？我有办法访问检查边界的元素吗？

正确的.at() 不做边界检查。要使用 边界检查的 2D 或 3D 子集，您需要使用 (x,y) 或 (x,y,z)。例如double b = A[1](5, 1, 2);

完美运行！再次感谢您。