10

我写了这个R函数，给定任意数量的向量 ( ...)，通过根据它们的名称对各自的元素值求和来组合它们。

add_vectors <- function(...) {
  a <- list(...)
  nms <- sort(unique(unlist(lapply(a, names))))
  out <- numeric(length(nms))
  names(out) <- nms
  for (v in a) out[names(v)] <- out[names(v)] + v

  out
}

例子：

v1 <- c(a=2,b=3,e=4)
v2 <- c(b=1,c=6,d=0,a=4)
add_vectors(v1, v2)
#
a b c d e 
6 4 6 0 4

我正在尝试编写一个更快的等效函数。

不幸的是，目前我不知道如何实现这一点，R所以我想Rcpp。但是，为了在Rcpp这个函数中进行转换，我错过了一些概念：

  1. 如何管理...参数。List使用in 类型的参数Rcpp？
  2. 如何迭代...参数中的向量。
  3. 如何通过名称访问（然后求和）向量的元素（这在 中非常简单R，但我不知道如何在 中进行操作Rcpp）。

因此，我正在寻找可以帮助我提高此功能的性能的人（在Ror中Rcpp，或两者兼而有之）。

任何帮助表示赞赏，谢谢。

4

5 回答 5

6

我会使用这样的东西：

#include <Rcpp.h>
using namespace Rcpp; 

// [[Rcpp::export]]
NumericVector add_all(List vectors){
    RCPP_UNORDERED_MAP<std::string,double> out ; 
    int n = vectors.size() ;
    for( int i=0; i<n; i++){
        NumericVector x = vectors[i] ;
        CharacterVector names = x.attr("names") ;
        int m = x.size() ;

        for( int j=0; j<m; j++){
            String name = names[j] ;
            out[ name ] += x[j] ;   
        }
    }
    return wrap(out) ;
}

使用以下包装器：

add_vectors_cpp <- function(...){
    add_all( list(...) )
}

RCPP_UNORDERED_MAP只是一个 typedef to unordered_map，std::取决于std::tr1::你的编译器，等等......

这里的技巧是...使用经典创建一个常规列表list(...)。

如果你真的想直接...在 C++ 中传递并在内部处理它，你将不得不使用.External接口。这很少使用，因此 Rcpp 属性不支持该.External接口。

使用.External，它看起来像这样（未经测试）：

SEXP add_vectors(SEXP args){
    RCPP_UNORDERED_MAP<std::string,double> out ; 
    args = CDR(args) ;
    while( args != R_NilValue ){
        NumericVector x = CAR(args) ;

        CharacterVector names = x.attr("names") ;
        int m = x.size() ;

        for( int j=0; j<m; j++){
            String name = names[j] ;
            out[ name ] += x[j] ;   
        }        
        args = CDR(args) ;
    }   
    return wrap(out) ;
}
于 2013-04-02T18:00:09.730 回答
3

使用编译器包编译成字节码会给你一些改进。此软件包随 R 一起提供。

library(compiler)
library(microbenchmark)

add_vectors_cmp <- cmpfun(add_vectors)

set.seed(1)
v <- rpois(length(letters), 10)
names(v) <- letters
vs <- replicate(150, v, simplify=FALSE)

not_compiled <- function(l) do.call(add_vectors, l)
compiled <- function(l) do.call(add_vectors_cmp, l)

plot(microbenchmark(not_compiled(vs), compiled(vs)))

在此处输入图像描述

于 2013-04-02T14:30:57.727 回答
3

我刚刚在Rcpp.

我不知道如何使用...参数（以及如何对其进行迭代），Rcpp所以我将此函数封装在一个简单的R函数中。

解决方案

library(Rcpp)
cppFunction(
  code = '
  NumericVector add_vectors_cpp(NumericVector v1, NumericVector v2) {
    // merging names, sorting them and removing duplicates
    std::vector<std::string> nms1 = v1.names();
    std::vector<std::string> nms2 = v2.names();
    std::vector<std::string> nms;
    nms.resize(nms1.size() + nms2.size());
    std::merge(nms1.begin(), nms1.end(), nms2.begin(), nms2.end(), nms.begin());
    std::sort(nms.begin(), nms.end());
    nms.erase(std::unique(nms.begin(), nms.end()), nms.end());
    // summing vector elements by their names and storing them in an associative data structure
    int num_names = nms.size();
    std::tr1::unordered_map<std::string, double> map(num_names);
    for (std::vector<int>::size_type i1 = 0; i1 != nms1.size(); i1++) {
        map[nms1[i1]] += v1[i1];
    }
    for (std::vector<int>::size_type i2 = 0; i2 != nms2.size(); i2++) {
        map[nms2[i2]] += v2[i2];
    }
    // extracting map values (to use as result vector) and keys (to use as result vector names)
    NumericVector vals(map.size());
    for (unsigned r = 0; r < num_names; ++r) {
        vals[r] = map[nms[r]];
    }
    vals.names() = nms;
    return vals;
  }',
  includes = '
  #include <vector>
  #include <tr1/unordered_map>
  #include <algorithm>'
)

然后封装在一个R函数中：

add_vectors_2 <- function(...) {
  Reduce(function(x, y) add_vectors_cpp(x, y), list(...))
}

请注意，此解决方案使用STL库。我不知道这是否是一个写得很好的 C++解决方案，或者是否可以（可能）编写一个更有效的解决方案，但可以肯定它是一个好的（和工作的）起点。

使用示例

v1 <- c(b = 1, d = 2, c = 3, a = 4, e = 6, f = 5)
v2 <- c(d = 2, c = 3, a = 4, e = 6, f = 5)
add_vectors(v1, v2, v1, v2)
#  a  b  c  d  e  f 
# 16  2 12  8 24 20
add_vectors_2(v1, v2, v1, v2)
#  a  b  c  d  e  f 
# 16  2 12  8 24 20

注意：此函数也适用于名称不是唯一的向量。

v1 <- c(b = 1, d = 2, c = 3, a = 4, e = 6, f = 5)
v2 <- c(d = 2, c = 3, a = 4, e = 6, f = 5, f = 10, a = 12)
add_vectors(v1, v2)
#  a  b  c  d  e  f 
# 16  1  6  4 12 15 
add_vectors_2(v1, v2)
#  a  b  c  d  e  f 
# 20  1  6  4 12 20

如上一个示例所示，即使输入向量具有非唯一名称，该解决方案也有效，将具有相同名称的相同向量的元素相加。

基准

R在最简单的情况下（两个向量），我的解决方案比解决方案快大约 3 倍。这是一个很好的改进，但可能还有更好的C++解决方案进行进一步小的改进的空间。

Unit: microseconds
                 expr    min     lq median      uq     max neval
  add_vectors(v1, v2) 65.460 68.569 70.913 73.5205 614.274   100
add_vectors_2(v1, v2) 20.743 23.389 25.142 26.9920 337.544   100

在此处输入图像描述

当将此函数应用于更多向量时，性能会有所下降（仅快 2 倍）。

Unit: microseconds
                                 expr     min       lq  median       uq     max neval
  add_vectors(v1, v2, v1, v2, v1, v1) 105.994 195.7565 205.174 212.5745 993.756   100
add_vectors_2(v1, v2, v1, v2, v1, v1)  66.168 125.2110 135.060 139.7725 666.975   100

所以现在的最后一个目标是删除R 直接使用....ListRcpp

我认为这是可能的，因为Rcpp糖具有与其相似的功能（例如功能的移植sapply），但希望得到一些反馈。

于 2013-04-02T16:55:15.350 回答
3

data.table 包非常擅长执行聚合和其他操作。我不是真正的专家，但

library(data.table)
add_vectors5 <- function(...)
{
    vals <- do.call(c, list(...))
    dt <- data.table(nm=names(vals), v=vals, key="nm")
    dt <- dt[,sum(v), by=nm]
    setNames(dt[[2]], dt[[1]])
}

似乎比其他纯 R 实现快约 2 倍。一个更神秘的实现是

add_vectors6 <- function(..., method="radix")
{
    vals <- do.call(c, list(...))
    ## order by name, but use integers for faster order algo
    idx <- match(names(vals), unique(names(vals)))
    o <- sort.list(idx, method=method, na.last=NA)

    ## cummulative sum of ordered values
    csum <- cumsum(vals[o])

    ## subset where ordering factor changes, and then diff
    idxo <- idx[o]
    diff(c(0, csum[idxo != c(idxo[-1], TRUE)]))
}

容易出现数值溢出；如果名称少于 100,000 个，则使用 method="radix"，如 上所暗示的那样?sort.list，否则使用 method="quick"。

于 2013-04-02T20:47:09.507 回答
1

我不认为你会得到太多的加速。我在 R 代码中采用了另一种方法，将所有输入组合成一个向量，然后按名称重新拆分，并使用vapply. 那里的所有函数或多或少都调用内部 C 代码，并且速度与大型向量的函数相当（在长度为 1e5 和 1e6 的向量上测试）。对于 3 或 4 个元素的玩具示例，它会慢一些。

add_vectors2 <- function(...) {
  y <- do.call(c, unname(list(...)))
  vapply(split(y, names(y)), sum, numeric(1))
}

#Longer sample vectors
m <- 1e3
n <- 1e6
v1 <- sample(m, n, replace = TRUE)
names(v1) <- sample(n)
v2 <- sample(m, n, replace = TRUE)
names(v2) <- sample(seq_len(n) + n / 2)  

#Timings
k <- 20
system.time(for(i in 1:k) add_vectors(v1, v2))   #5 or 6 seconds
system.time(for(i in 1:k) add_vectors2(v1, v2))  #same

编辑：矢量名称固定为唯一，反映了 Roland 的评论。我的解决方案现在比 OP 慢一点。

于 2013-04-02T15:24:44.743 回答