找出R和Rcpp中最大的n个唯一值及其频率
Find the largest n unique values and their frequencies in R and Rcpp
我有一个数字向量v（已经省略了NA），希望获得第n个最大值及其各自的频率。
我发现http://gallery.rcpp.org/articles/top-elements-from-vectors-using-priority-queue/速度相当快。
// [[Rcpp::export]]
std::vector<int> top_i_pq(NumericVector v, unsigned int n)
{
typedef pair<double, int> Elt;
priority_queue< Elt, vector<Elt>, greater<Elt> > pq;
vector<int> result;
for (int i = 0; i != v.size(); ++i) {
if (pq.size() < n)
pq.push(Elt(v[i], i));
else {
Elt elt = Elt(v[i], i);
if (pq.top() < elt) {
pq.pop();
pq.push(elt);
}
}
}
result.reserve(pq.size());
while (!pq.empty()) {
result.push_back(pq.top().second + 1);
pq.pop();
}
return result ;
}
然而，这种关系不会得到尊重。事实上，我不需要索引，返回值也可以。
我想得到的是一个包含值和频率的列表，比如：
numv <- c(4.2, 4.2, 4.5, 0.1, 4.4, 2.0, 0.9, 4.4, 3.3, 2.4, 0.1)
top_i_pq(numv, 3)
$lengths
[1] 2 2 1
$values
[1] 4.2 4.4 4.5
无论是获得唯一向量、表还是（完整）排序都不是一个好主意，因为与v的长度相比，n通常很小（可能很容易大于1e6）。
到目前为止的解决方案是：
library(microbenchmark)
library(data.table)
library(DescTools)
set.seed(1789)
x <- sample(round(rnorm(1000), 3), 1e5, replace = TRUE)
n <- 5
microbenchmark(
BaseR = tail(table(x), n),
data.table = data.table(x)[, .N, keyby = x][(.N - n + 1):.N],
DescTools = Large(x, n, unique=TRUE),
Coatless = ...
)
Unit: milliseconds
expr min lq mean median uq max neval
BaseR 188.09662 190.830975 193.189422 192.306297 194.02815 253.72304 100
data.table 11.23986 11.553478 12.294456 11.768114 12.25475 15.68544 100
DescTools 4.01374 4.174854 5.796414 4.410935 6.70704 64.79134 100
嗯，DescTools仍然是最快的，但我相信Rcpp可以显著改进它（因为它是纯R）！
我想用另一种基于Rcpp的解决方案来挑战，它比DescTools方法快约7倍，比data.table方法快约13倍，使用上面的1e5长度的x和n = 5样本数据。实现过程有点长，所以我将以基准测试为例：
fn.dt <- function(v, n) {
data.table(v = v)[
,.N, keyby = v
][(.N - n + 1):.N]
}
microbenchmark(
"DescTools" = Large(x, n, unique=TRUE),
"top_n" = top_n(x, 5),
"data.table" = fn.dt(x, n),
times = 500L
)
# Unit: microseconds
# expr min lq mean median uq max neval
# DescTools 3330.527 3790.035 4832.7819 4070.573 5323.155 54921.615 500
# top_n 566.207 587.590 633.3096 593.577 640.832 3568.299 500
# data.table 6920.636 7380.786 8072.2733 7764.601 8585.472 14443.401 500
更新
如果您的编译器支持C++11，则可以利用std::priority_queue::emplace（与下面的C++98版本相比）显著提高性能。我不会发布这个版本，因为它基本上是相同的，除了对std::move和emplace的一些调用，但这里有一个链接
针对前三个函数进行测试，并使用data.table 1.9.7（比1.9.6快一点）生成
print(res2, order = "median", signif = 3)
# Unit: relative
# expr min lq mean median uq max neval cld
# top_n2 1.0 1.00 1.000000 1.00 1.00 1.00 1000 a
# top_n 1.6 1.58 1.666523 1.58 1.75 2.75 1000 b
# DescTools 10.4 10.10 8.512887 9.68 7.19 12.30 1000 c
# data.table-1.9.7 16.9 16.80 14.164139 15.50 10.50 43.70 1000 d
其中CCD_ 9是C++11版本。
top_n功能实现如下：
#include <Rcpp.h>
#include <utility>
#include <queue>
class histogram {
private:
struct paired {
typedef std::pair<double, unsigned int> pair_t;
pair_t pair;
unsigned int is_set;
paired()
: pair(pair_t()),
is_set(0)
{}
paired(double x)
: pair(std::make_pair(x, 1)),
is_set(1)
{}
bool operator==(const paired& other) const {
return pair.first == other.pair.first;
}
bool operator==(double other) const {
return is_set && (pair.first == other);
}
bool operator>(double other) const {
return is_set && (pair.first > other);
}
bool operator<(double other) const {
return is_set && (pair.first < other);
}
paired& operator++() {
++pair.second;
return *this;
}
paired operator++(int) {
paired tmp(*this);
++(*this);
return tmp;
}
};
struct greater {
bool operator()(const paired& lhs, const paired& rhs) const {
if (!lhs.is_set) return false;
if (!rhs.is_set) return true;
return lhs.pair.first > rhs.pair.first;
}
};
typedef std::priority_queue<
paired,
std::vector<paired>,
greater
> queue_t;
unsigned int sz;
queue_t queue;
void insert(double x) {
if (queue.empty()) {
queue.push(paired(x));
return;
}
if (queue.top() > x && queue.size() >= sz) return;
queue_t qtmp;
bool matched = false;
while (queue.size()) {
paired elem = queue.top();
if (elem == x) {
qtmp.push(++elem);
matched = true;
} else {
qtmp.push(elem);
}
queue.pop();
}
if (!matched) {
if (qtmp.size() >= sz) qtmp.pop();
qtmp.push(paired(x));
}
std::swap(queue, qtmp);
}
public:
histogram(unsigned int sz_)
: sz(sz_),
queue(queue_t())
{}
template <typename InputIt>
void insert(InputIt first, InputIt last) {
for ( ; first != last; ++first) {
insert(*first);
}
}
Rcpp::List get() const {
Rcpp::NumericVector values(sz);
Rcpp::IntegerVector freq(sz);
R_xlen_t i = 0;
queue_t tmp(queue);
while (tmp.size()) {
values[i] = tmp.top().pair.first;
freq[i] = tmp.top().pair.second;
++i;
tmp.pop();
}
return Rcpp::List::create(
Rcpp::Named("value") = values,
Rcpp::Named("frequency") = freq);
}
};
// [[Rcpp::export]]
Rcpp::List top_n(Rcpp::NumericVector x, int n = 5) {
histogram h(n);
h.insert(x.begin(), x.end());
return h.get();
}
在上面的histogram类中有很多内容，但只涉及一些关键点：
paired类型本质上是一个围绕std::pair<double, unsigned int>的包装类，它将值与计数相关联，提供了一些方便的功能，如用于计数的直接前/后增量的operator++()/operator++(int)，以及修改的比较运算符histogram类封装了一种"托管"优先级队列，即std::priority_queue的大小上限为特定值sz- 我没有使用
std::priority_queue的默认std::less排序，而是使用了一个大于比较器，以便可以对照std::priority_queue::top()检查候选值，以快速确定它们是否应该（a）被丢弃，（b）替换队列中的当前最小值，或（c）更新队列中现有值之一的计数。这是可能的，因为队列的大小被限制为<=CCD_ 22
我敢打赌data.table是有竞争力的：
library(data.table)
data <- data.table(v)
data[ , .N, keyby = v][(.N - n + 1):.N]
其中n是您想要获得的数字
注意：以前的版本复制了table()的功能，而不是目标。此版本已被删除，将在场外提供。
地图攻击计划
以下是使用map的解决方案。
C++98
首先，我们需要找到数字向量的"唯一"值。
为此，我们选择将被计数为key的数字存储在std::map中，并在每次观察到该数字时递增value。
使用std::map的排序结构，我们知道顶部的n数字在std::map的后面。因此，我们使用迭代器来弹出这些元素，并将它们导出到数组中。
C++11
如果可以访问C++11编译器，另一种选择是使用std::unordered_map，它的大O为O(1)用于插入和检索（如果哈希不正确，则为O(n)），而std::map的大O则为O(log(n))。
为了获得正确的顶部n，可以使用std::partial_sort()来这样做
实施
C++98
#include <Rcpp.h>
// [[Rcpp::export]]
Rcpp::List top_n_map(const Rcpp::NumericVector & v, int n)
{
// Initialize a map
std::map<double, int> Elt;
Elt.clear();
// Count each element
for (int i = 0; i != v.size(); ++i) {
Elt[ v[i] ] += 1;
}
// Find out how many unique elements exist...
int n_obs = Elt.size();
// If the top number, n, is greater than the number of observations,
// then drop it.
if(n > n_obs ) { n = n_obs; }
// Pop the last n elements as they are already sorted.
// Make an iterator to access map info
std::map<double,int>::iterator itb = Elt.end();
// Advance the end of the iterator up to 5.
std::advance(itb, -n);
// Recast for R
Rcpp::NumericVector result_vals(n);
Rcpp::NumericVector result_keys(n);
unsigned int count = 0;
// Start at the nth element and move to the last element in the map.
for( std::map<double,int>::iterator it = itb; it != Elt.end(); ++it )
{
// Move them into split vectors
result_keys(count) = it->first;
result_vals(count) = it->second;
count++;
}
return Rcpp::List::create(Rcpp::Named("lengths") = result_vals,
Rcpp::Named("values") = result_keys);
}
短期测试
让我们通过运行一些数据来验证它是否有效：
# Set seed for reproducibility
set.seed(1789)
x <- sample(round(rnorm(1000), 3), 1e5, replace = TRUE)
n <- 5
现在我们寻求获得发生信息：
# Call our function
top_n_map(a)
给我们：
$lengths
[1] 101 104 101 103 103
$values
[1] 2.468 2.638 2.819 3.099 3.509
基准
Unit: microseconds
expr min lq mean median uq max neval
BaseR 112750.403 115946.7175 119493.4501 117676.2840 120712.595 166067.530 100
data.table 6583.851 6994.3665 8311.8631 7260.9385 7972.548 47482.559 100
DescTools 3291.626 3503.5620 5047.5074 3885.4090 5057.666 43597.451 100
Coatless 6097.237 6240.1295 6421.1313 6365.7605 6528.315 7543.271 100
nrussel_c98 513.932 540.6495 571.5362 560.0115 584.628 797.315 100
nrussel_c11 489.616 512.2810 549.6581 533.2950 553.107 961.221 100
正如我们所看到的，这个实现击败了data.table，但却成为DescTools和@nrussel尝试的牺牲品。
相关文章：
- 何时在引用或唯一指针上使用移动语义
- 芬威克树(BIT）.找到具有给定累积频率的最小索引，单位为 O(logN）
- 在指针的帮助下，文本文件中单词的频率
- C++优先级队列，按对象的唯一指针的特定方法升序排列
- 计算排序向量的向量中唯一值的计数
- 如何使用Visual Studio 2017在C++中为参数化对象数组使用唯一指针
- 通过组合不同的类型来创建唯一的id
- 函数计算用户按下按钮的频率
- 使用Unique_ptr确保工厂中的对象唯一
- c++多进程编写一个唯一的文件
- 在以唯一ptr为值的C++映射中，动态内存何时会被销毁
- 在条件变量中触发错误信号的频率是多少
- 如何更改唯一指针向量的可见性
- 在C++的两个字符串中连接以逗号分隔的唯一值
- 共享指针和具有自定义删除程序的唯一指针之间的语法差异背后的任何原因
- C++尝试深度复制唯一指针时出现内存访问冲突
- 具有引用成员的结构是否具有唯一的对象表示形式
- 使用标准库对矢量的矢量进行排序，并计算唯一出现的频率
- 找出R和Rcpp中最大的n个唯一值及其频率
- 输出“std::multiset”的唯一元素及其频率，使用C++中的std::算法(无循环）
