犰狳矩阵列的非连续访问

Posted 2023-02-16

【中文标题】犰狳矩阵列的非连续访问

【英文标题】：noncontiguous access of Armadillo matrix columns

【发布时间】：2020-11-25 03:13:36

【问题描述】：

我有一个矩阵的索引存储在一个列表中，我想访问 Rcpp 中的这些列。

library(Rcpp)
library(RcppArmadillo)
cppFunction('arma::mat tmp(arma::mat mat1, List list1) 
            arma::mat mat2=mat1.row(1);
return mat2.cols(as<arma::uvec>(list1[1]));

            ',depends = "RcppArmadillo")

但是，上述函数会导致致命错误（“R sessions Aborted”）。例如，

tmp(matrix(1:4,2),list(list(c(0))))

纯 rcpp 中的类似函数会返回此错误。

note: in instantiation of function template specialization 'Rcpp::Vector<19, PreserveStorage>::create<Rcpp::traits::named_object<arma::subview_row<double> >, Rcpp::traits::named_object<arma::subview_elem2<double, arma::Mat<unsigned int>, arma::Mat<unsigned int> > >, Rcpp::traits::named_object<Rcpp::internal::generic_proxy<19, PreserveStorage> > >' requested here

通过documentation，我应该可以通过arma:uvec 访问这些列。

X.cols( vector_of_column_indices )

我不知道为什么它仍然要求我使用这种格式X.cols( first_col, last_col )。 非常感谢任何帮助！

【参考方案1】：

这个问题不是特别清楚（可能会受益于一些编辑/澄清），但这里有一些关键点：

R 以及无数其他应用程序的内存模型很简单，因为向量需要连续内存 矩阵实际上只是一个带有二维索引的向量，因此同样成立：需要连续的内存 R 中的列表不同于矩阵，因为它们实际上是列列表 data.frames 也是如此，它们在某种程度上是子类列表

所以有些事情你可以做，也有些事情不能做。 Armadillo documentation 很好，可以帮助将部分向量堆叠和组合成更大的向量。

编辑：这是一个非常简单的例子，展示了如何解决这个问题。因为 Armadillo 显然想要一个 matrix 来作为子集，而 uvec 带有索引，我们……只需给它一个矩阵和一个（无符号）整数向量（RcppArmadillo 为我们映射）。所以你的主要问题真的是你如何从List 到 uvec --- 它不能以你所拥有的形式工作。我建议更明确地说明如何创建向量并正确实例化一个。

代码

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat subsetOfMatrix(arma::mat X, arma::uvec v) 
  return X.cols(v);


/*** R
M <- matrix(1:16,4,4)
subsetOfMatrix(M, c(0L, 2L))
*/

输出

R> Rcpp::sourceCpp("~/git/***/64998018/answer.cpp")

R> M <- matrix(1:16,4,4)

R> subsetOfMatrix(M, c(0L, 2L))
     [,1] [,2]
[1,]    1    9
[2,]    2   10
[3,]    3   11
[4,]    4   12
R>