beachmat 2.15.0
的beachmat包提供了一个c++ API来处理由DelayedArray的块处理机制。这使得Bioconductor包能够使用c++对存储在几乎任何矩阵表示中的数据进行高性能处理，从密集的普通矩阵，稀疏的矩阵到矩阵包、文件支持的矩阵或甚至基于云的矩阵(例如，HDF5Array，TileDBArray)．包使用DelayedArray将数据块以稀疏或密集矩阵的形式存储到内存中，然后可以很容易地被c++代码使用beachmat的API。其目的是在R中抽象出矩阵对象的特定类(甚至类型!)，以方便c++扩展的开发。
早在2017年，关于如何使各种c++函数与替代矩阵表示一起工作，我们有两种广泛的策略。的版本1和2中实现了第一种方法beachmat，这需要Rcpp: RObject在c++中，计算出底层的表示，然后选择适当的方法从给定的行或列中读取值的向量。这工作得相当好，并且允许人们以一种与矩阵表示无关的方式编写c++代码。最重要的是，它是一个即插即用的替代品Rcpp.column ()而且.row ()方法用于矩阵访问，这允许客户端包快速启动并与现有的c++代码一起运行。
第二种策略是在R中使用DelayedArray: blockApply．这将读入一个值块作为普通矩阵，然后将其传递到c++代码中;这意味着开发人员只需要担心在他们bob电竞体育官网的c++代码中处理普通的(column-major等)矩阵。问题是，实现密集矩阵会牺牲稀疏矩阵格式的一些效率。此外，实现这将涉及更多的重构，因为循环将在c++和R之间共享——在R中对每个块进行一次循环，在c++中对每个块进行另一次循环，以完成所需的计算。
方程现在有了适当的稀疏支持DelayedArray．块处理现在可以生成密集和稀疏的子矩阵(取决于底层的表示)，从而避免了以前从强制转换到密集格式的效率损失。使用基于r的区块循环，我们可以很容易地访问BiocParallel基于并行化;我们避免了从c++回调到R的讨厌的调用，这对于处理不属于已知类型的矩阵是必要的;我们得到了一种更便于用户调整的处理大型矩阵输出的方法RealizationSinks.块处理现在是处理跨行或跨列独立计算的大型矩阵的首选机制。
那么，如何beachmat适应这种新模式?在块处理循环中，我们可能希望调用一个c++函数，该函数接受块并对其进行操作。这个块可以是密集的，也可以是稀疏的，但是c++函数可能并不关心它是否只是想获取，例如，单个的行向量以供进一步处理。在这种情况下，beachmat可以继续提供来自已实现块的不可知数据访问，这样开发人员就不必担心在c++代码中处理不同的格式。bob电竞体育官网或者，如果稀疏格式适合更有效的算法，beachmat还提供专门的稀疏类，只对非零值进行操作。
的beachmatPackage目前有几个依赖项:
编译器应该支持c++ 11标准。这通常需要GCC 4.8.1或更高版本。您需要告诉构建系统使用c++ 11，通过修改SystemRequirements字段描述文件:
SystemRequirements: c++ 11Rcpp应安装。你还需要确保Rcpp在加载包时加载。这需要添加Rcpp到进口字段描述文件:
进口:Rcpp和一个对应的importFrom《名称空间文件:
importFrom (Rcpp sourceCpp)(要导入的确切函数并不重要，只要加载了名称空间。请查看Rcpp有关详细信息的文档。)
beachmat是一个头文件库，所以链接就像写一样简单:
链接到:Rcpp，海滩垫…在你的描述文件。(这确保了包含Rcpp头)。
在C或c++代码文件中，使用标准技术来包含头文件定义，例如:# include“beachmat3 / numeric_matrix.h”(见在这里有关详情)。头文件可在以下位置阅读(在R会话中输入):
系统。文件(包= " beachmat”,“包括”)## [1] "/tmp/RtmphAMjmS/Rinst7c87b6092a494/beachmat/include"如果你看到在“使用”之前预期的非限定id错误，这通常意味着用于编译的GCC版本beachmat是过时的，不完全支持c++ 11。
API本身全面的文档我们就不详细讲了。相反，我们将关注最常见的使用模式。假设我们想要从一个逻辑的、整数的或数字矩阵类对象中计算列和:
#include "beachmat3/beachmat.h" #include #include //在包上下文中不是必需的:// [[Rcpp::depends(beachmat)]] // [[Rcpp::export(rng=false)]] Rcpp::NumericVector column_sum (Rcpp::RObject mat) {auto ptr = beachmat::read_lin_block(mat);std::向量<二>工作区(ptr - > get_nrow ());Rcpp: NumericVector输出(ptr - > get_ncol ());For (size_t I = 0;I < ptr->get_ncol();++i) {auto colptr = ptr->get_col(i, workspace.data());Output [i] = std::accumulate(colptr, colptr + ptr->get_nrow()， 0.0);}返回输出;｝ 的read_lin_block ()函数设置一个指针，该指针指向可使用get_col ()方法(或get_row ()，如果我们对行感兴趣的话)。这可能涉及到内存分配，因此我们为该方法提供一个指向“工作区”的指针，以便它在其中放入值。的get_col ()方法然后返回一个指针，该指针指向每个请求列的数据值向量我，我们积累并储存在输出返回到R会话。为了演示:
column_sums(矩阵(rnorm (1000), ncol = 10))[1] 3.331864 3.855267 -20.649126 -9.746339 -1.470335 16.693330 ## [7] -2.675940 -3.898415 14.784445 8.219382此代码可用于生成的密集或稀疏块DelayedArray块处理。此外，它自动支持可能存在的所有数字类型数据垫，即逻辑，整数或双精度值。所有值都转换为工作空间的类型—在本例中为双精度向量，但同样可以通过创建int工作区。这意味着我们可以避免编写多个版本的代码来处理不同类型的输入。
# Integer column_sum (matrix(rpois(1000, 10)， ncol=10))## [1] 1038 971 1040 999 1010 944 1000 1024 1015 973# column_sum (matrix(rbinom(1000, 1, 0.5)==1, ncol=10))## [1] 47 48 50 53 45 45 47 54 47 46#稀疏column_sum (Matrix::rsparsematrix(100,10,0.1))## [1] -1.1600 2.2647 -7.2530 -0.1040 1.2450 -2.7100 2.4640 -1.5110 5.6700 ## [10] -0.1500返回的指针get_col ()可以指向工作空间，并且只在工作空间下一次使用之前有效。用户应该只依赖于返回的指针之前的任何后续调用get_col ()．此外，不能保证工作空间实际上是用数据值填充的，因为一些矩阵表示可能不需要复制过程来访问底层数据。如果需要复制到工作区，可以使用以下模式:
#include "beachmat3/beachmat.h" #include #include //在包上下文中不是必需的:// [[Rcpp::depends(beachmat)]] // [[Rcpp::export(rng=false)]] Rcpp::NumericVector column_sums_copy (Rcpp::RObject mat) {auto ptr = beachmat::read_lin_block(mat);std::向量<二>工作区(ptr - > get_nrow ());Rcpp: NumericVector输出(ptr - > get_ncol ());For (size_t I = 0;I < ptr->get_ncol();++i) {auto colptr = ptr->get_col(i, workspace.data());//检查colptr是否已经指向工作区。if (colptr != workspace.data()){//如果不是，强制将内容复制到工作区。Std::copy(colptr, colptr + ptr->get_nrow()， workspace.data());} output[i] = std::accumulate(workspace.begin()， workspace.end()， 0.0); } return output; } 对于某些应用程序，可以为仅使用非零项的稀疏数据实现更有效的算法。我们可以用read_lin_sparse_block ()来创建一个指向稀疏c++矩阵的指针，该指针将重载get_col ()而且get_row ()方法返回非零项的列表。下面演示了仅通过迭代非零项来计算列和。注意，这次我们需要两个工作区向量，一个用于非零值，另一个用于它们的索引。(同样，不能保证填充了工作空间，因此如果有必要，应该执行手动复制。)
#include "beachmat3/beachmat.h" #include #include //在包上下文中不是必需的:// [[Rcpp::depends(beachmat)]] // [[Rcpp::export(rng=false)]] Rcpp::NumericVector column_sums_sparse (Rcpp::RObject mat) {auto ptr = beachmat::read_lin_sparse_block(mat);std::向量<二> workspace_x (ptr - > get_nrow ());std::向量< int > workspace_i (ptr - > get_nrow ());Rcpp: NumericVector输出(ptr - > get_ncol ());For (size_t I = 0;I < ptr->get_ncol();++i) {auto indexes = ptr->get_col(i, workspace_x.data()， workspace_i.data());Auto XPTR = indexes .x;Auto iptr = index .i;//行索引，这里不使用。 auto nnzero = indices.n; output[i] = std::accumulate(xptr, xptr + nnzero, 0.0); } return output; } 为了演示:
column_sums_sparse(Matrix::rsparsematrix(100,10,0.1))## [1] 0.733 -0.341 -1.730 0.535 0.338 3.085 -1.950 -1.096 1.589 4.850try(column_sums_sparse(matrix(rpois(1000, 10)， ncol=10)))## column_sums_sparse(matrix(rpois(1000, 10)， ncol = 10)): ##对象没有class属性在实践中，我们希望c++函数能够同时处理来自块处理的稀疏和密集输入。我们可以通过检查矩阵is_sparse ()，如果是，则将其提升为稀疏实例，如下所示。这使得开发人员可以根据观bob电竞体育官网察到的数据表示顺畅地在算法之间切换。
#include "beachmat3/beachmat.h" #include #include //在包上下文中不是必需的:// [[Rcpp::depends(beachmat)]] // [[Rcpp::export(rng=false)]] Rcpp::NumericVector column_sums_flexible (Rcpp::RObject mat) {auto ptr = beachmat::read_lin_block(mat);Rcpp: NumericVector输出(ptr - > get_ncol ());std::向量<二>工作区(ptr - > get_nrow ());If (ptr->is_sparse()) {auto SPTR = beachmat::promote_to_sparse(ptr);//注意，ptr变成NULL。std::向量< int > workspace_i (sptr - > get_nrow ());For (size_t I = 0;I < sptr->get_ncol();++i) {auto indexes = sptr->get_col(i, workspace.data()， workspace_i.data());Auto XPTR = indexes .x; output[i] = std::accumulate(xptr, xptr + indices.n, 0.0); } } else { for (size_t i = 0; i < ptr->get_ncol(); ++i) { auto vec = ptr->get_col(i, workspace.data()); output[i] = std::accumulate(vec, vec + ptr->get_nrow(), 0.0); } } return output; } 为了演示:
column_sums_flexible(矩阵::rsparsematrix(100, 10, 0.1))## [1] 1.445 -0.473 0.190 3.580 -0.852 0.684 2.608 4.000 0.870 -4.547Column_sums_flexible(矩阵(rpois(1000, 10)， ncol=10))## [1] 1013 982 950 1012 963 999 1040 996 984 938在大多数情况下，任何输出矩阵都应该使用Rcpp:矩阵对象。如果预计输出量很大，我们建议使用DelayedArrayRealizationSink存储结果的基础设施，例如，跨块处理迭代的文件支持矩阵。
有时，可能需要生成稀疏输出beachmat提供一些有用的实用程序。当预先不知道非零项的数量时，第一种更通用的方法是相关的。这假设您可以将所有非零项收集到a中std::地图存储:
#include "beachmat3/beachmat.h" #include 为了演示:
generate_sparse_general ()## 7 x 10稀疏矩阵类“dgCMatrix”## ##[1，]。0.1 . . . . . . . .## [2，] . . . . . . . . . .## [3，] . . . . . . . . . .## [4，] 0.3 . . . . . . . . .## [5，] . . . . .0.4 . . . .## [6，] . .0.2 . . . . . . .## [7，] . . . . . . . . . .另一种更具体的方法是替换现有对象的非零值* gCMatrix对象的另一个向量对应于与原始值相同的顺序和位置的新值。这对于不改变非零项位置的操作更有效。
#include "beachmat3/beachmat.h" #include 为了演示:
library(Matrix) mat <- rsparsematrix(10,10,0.1) generate_sparse_specific(mat)## 10 × 10稀疏矩阵类“lgCMatrix”## ## [1，]. . . .| | . . . .##[2，]…|……|。## [3，] . . . . . . . . . .## [4，] . . . . .| . . . .## [5，] . . . . . . . . . .## [6，] . . | . . . . . . | ## [7,] . . . . . | . . . . ## [8,] . . . . . . . . . | ## [9,] . . . . . . . . . . ## [10,] . . . . . | . . . .一旦我们写出了我们想要的函数，我们就可以执行了DelayedArray对任何类似矩阵的输入进行块处理。可以直接使用blockApply ()函数，或者我们可以用beachmat的包装器来方便地将输入矩阵分割成连续的行或列块:
library(beachmat) blockedcolsum <- function(x，…){out <- colBlockApply(x, FUN=column_sum，…)unlist(out) #跨块组合结果。｝这个函数现在能够对任何可以进行块处理的矩阵类对象进行操作。粗略地说，如果一个物体有两个维度和一个as.array ()方法，即可正确处理。
#普通矩阵:x1 <- matrix(rnorm(1000)， ncol=20) blockedcolsum (x1)[1] -9.7145327 -1.0007139 -17.0405649 0.9958055 ## [7] -10.6014081 -0.1348006 -0.1747428 -2.6273769 -8.8778876 0.9401495 ## [13] -3.6048843 -2.8983927 2.9849290 14.3433742 -4.6471962 [19] -1.2189234 -4.4202641#适用于矩阵类:x2 <-矩阵(x1) blockedcolsum (x2)[1] -9.7145327 -1.0007139 -17.0405649 0.9958055 ## [7] -10.6014081 -0.1348006 -0.1747428 -2.6273769 -8.8778876 0.9401495 ## [13] -3.6048843 -2.8983927 2.9849290 14.3433742 -4.6471962 [19] -1.2189234 -4.4202641# DelayedArray对象:library(DelayedArray) x3 <- DelayedArray(x1) blockedcolsum (x3)[1] -9.7145327 -1.0007139 -17.0405649 0.9958055 ## [7] -10.6014081 -0.1348006 -0.1747428 -2.6273769 -8.8778876 0.9401495 ## [13] -3.6048843 -2.8983927 2.9849290 14.3433742 -4.6471962 [19] -1.2189234 -4.4202641如上所述，这个过程可以通过简单的传递轻松地并行化BPPARAM =到colBlockApply ()函数。请注意，这只对大型矩阵有好处，其中并行化的好处抵消了相关的开销。
library(BiocParallel) blockedcolsum (x3, BPPARAM=MulticoreParam(2))[1] -9.7145327 -1.0007139 -17.0405649 0.9958055 ## [7] -10.6014081 -0.1348006 -0.1747428 -2.6273769 -8.8778876 0.9401495 ## [13] -3.6048843 -2.8983927 2.9849290 14.3433742 -4.6471962 [19] -1.2189234 -4.4202641sessionInfo ()## R正在开发中(不稳定)(2022-10-25 r83175) ##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 22.04.1 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas。so ## LAPACK: /usr/lib/x86_64-linux-gnu/ LAPACK /liblapack.so.3.10.0 ## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基础包:## [1]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:## [1]BiocParallel_1.33.0 DelayedArray_0.25.0 IRanges_2.33.0 ## [4] S4Vectors_0.37.0 MatrixGenerics_1.11.0 matrixStats_0.62.0 ## [7] BiocGenerics_0.45.0 beachmat_2.15.0 Matrix_1.5-1 ## [10] knitr_1.40 BiocStyle_2.27.0 ## ##通过命名空间加载(且未附加):## [4] stringi_1.7.8 jsonlite_1.8.3 htmltools_0.5.3 ## [7] sass_0.4.2 rmarkdown_2.17 grid_4.3.0 ## [10] evaluate_0.17 jquerylib_0.1.4 fastmap_1.1.0 ## [13] yaml_2.3.6 bookdown_0.29 BiocManager_1.30.19 ## [19] rcpp_1 .4.1 compiler_4.3.0 codetools_0.2-18 ## [19] Rcpp_1.0.9 lattice_0.20-45 digest_0.6.30 ## [25] bslib_0.4.0 tools_4.3.0 cachem_1.0.6 ## [25] bslib_0.4.0