rpy2 和标准 R 矩阵之间的差异

Discrepancy between rpy2 and standard R matrices

我的目标

所以我试图从 preprocessCore 脚本中使用 rpy2 包在一个 巨大的 矩阵上(10 gb+ 文件)。我的记忆力几乎是无限的。当我 运行 通过 R 本身使用以下内容时,我能够完成规范化并打印输出:

require(preprocessCore);
all <- data.matrix(read.table("data_table.txt",sep="\t",header=TRUE));
all[,6:57]=normalize.quantiles(all[,6:57]);
write.table(all,"QN_data_table.txt",sep="\t",row.names=FALSE);

我正在尝试将其构建到 python 脚本中,该脚本还使用 rpy2 python 包执行其他操作,但我在构建方式上遇到了问题矩阵。示例如下:

matrix = sample_list  # My 2-d python array containing the data.
v = robjects.FloatVector([ element for col in matrix for element in col ])
m = robjects.r['matrix'](v, ncol = len(matrix), byrow=False)
print("Performing quantile normalization.")
Rnormalized_matrix = preprocessCore.normalize_quantiles(m)
norm_matrix = np.array(Rnormalized_matrix)

return header, pos_list, norm_matrix

问题

这适用于较小的文件,但是当我 运行 它用于我的大文件时,它会因错误而死:rpy2.rinterface.RRuntimeError: Error: cannot allocate vector of size 9.7 Gb

我知道 R 的向量的最大大小是 8 Gb,这解释了为什么会抛出上述错误。 rpy2 docs 说:

"A Matrix is a special case of Array. As with arrays, one must remember that this is just a vector with dimension attributes (number of rows, number of columns)."

我有点想知道它是如何严格遵守这一点的,所以我更改了我的代码以初始化我想要的大小的矩阵,然后遍历并将元素分配给值:

matrix = sample_list  # My 2-d python array of data.
m_count = 1
m = robjects.r['matrix'](0.0, ncol=len(matrix), nrow=len(matrix[0]))
for samp in matrix:
    i_count = 1
    for entry in samp:
        m.rx[i_count, m_count] = entry  # Assign the data to the element.
        i_count += 1
    m_count += 1

print("Performing quantile normalization.")

Rnormalized_matrix = preprocessCore.normalize_quantiles(m)
norm_matrix = np.array(Rnormalized_matrix)

return header, pos_list, norm_matrix

同样,这适用于较小的文件,但会因与之前相同的错误而崩溃。

所以我的问题是 允许在 R 中分配巨大的矩阵但在 rpy 中导致问题的根本区别是什么?我需要用不同的方法来解决这个问题吗?我应该接受它并在 R 中做吗?或者有什么方法可以避免我遇到的问题?

从根本上说,R 是一种函数式语言。这意味着在 R

中做
m[i, j] <- 123

发生的事情是这样的

assign_to_cell <- `[<-`
m <- assign_to_cell(m, i, j, 123)

其中参数按值传递。

这意味着应该返回一个新矩阵 m,并且 (i,j) 处的单元格包含新值。为每个赋值制作一个 m 的副本将是相当低效的,尤其是当你正在经历它时对于较大的对象,因此 R 解释器有一个巧妙的技巧(有关详细信息,请参阅 R 的 C 源代码):左侧将表达式与表达式的右侧进行比较,如果对象相同,解释器将知道不需要复制,可以进行修改 "in place".

现在使用 rpy2 还有两点需要考虑:虽然 Python 是(主要)通过引用传递参数,但嵌入式 R 引擎无法知道左边发生了什么Python 表达式的一侧。

表达式

m.rx[i_count, m_count] = entry 

正在忠实地构建一个 R 调用,例如

m <- assign_to_cell(m, i, j, entry)

但是 R 可以在表达式的左侧向前看。结果是制作了 m 的副本。每次修改。

但是,rpy2 可以轻松地将 R 中定义的向量、矩阵和数组移动到 Python 的引用传递世界。例如,这些 R 对象可以别名为相应的 numpy 对象(使用 asarra - 请参阅 http://rpy.sourceforge.net/rpy2/doc-2.0/html/numpy.html#low-level-interface)。记住 R 数组是按列优先顺序排列的,还可以计算索引并跳过别名到 numpy 数组并就地修改单元格:

m[idx] = entry 

注:

我认为由于 R 索引是 32 位整数(如果我没记错的话)导致的 8Gb 限制在几年前就已经取消了。 您的无限内存可能比您想象的要少。系统上的物理内存并不一定意味着可以在一个连续的块中分配所有物理内存。