在 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
获得唯一向量、table 或(完整)排序都不是一个好主意,因为与 v 的长度(可能很容易 >1e6)相比,n 通常较小。
目前的解决方案是:
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)!
我敢打赌 data.table
很有竞争力:
library(data.table)
data <- data.table(v)
data[ , .N, keyby = v][(.N - n + 1):.N]
其中 n
是您要获取的号码
注意:以前的版本为 table()
而不是目标复制了功能。此版本已被删除,将在场外提供。
绘制攻击计划
下面是使用 map
.
的解决方案
C++98
首先,我们需要找到数字向量的 "unique" 值。
为此,我们选择将被计数为 key
的数字存储在 std::map
中,并在每次观察到该数字时递增 value
。
使用 std::map
的排序结构,我们知道前 n
个数字在 std::map
的后面。因此,我们使用迭代器弹出这些元素并将它们导出到数组中。
C++11
如果可以访问 C++11 编译器,另一种方法是使用 std::unordered_map
,它具有 O(1)
的大 O 用于插入和检索(O(n)
如果坏散列)与 std::map
相比,后者具有 O(log(n))
的大 O。
要获得正确的顶部 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 尝试的牺牲品。
我想用另一个基于 Rcpp 的解决方案来证明自己的能力,它比 DescTools
方法快 ~7 倍,比 data.table
方法快 ~13 倍,使用上面的 1e5-length 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 版本相比)。我不会 post 这个版本,因为它大部分是相同的,除了对 std::move
和 emplace
的几次调用,但是 here's a link to it。
针对前三个函数进行测试,并使用 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
其中 top_n2
是 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
class 中发生了很多事情,但只是触及一些关键点:
paired
类型本质上是 class 围绕 std::pair<double, unsigned int>
的包装器,它将值与计数相关联,提供一些便利功能,例如 operator++()
/ operator++(int)
用于直接 pre-/post-increment 的计数,并修改了比较运算符。
-
histogram
class 包装了一种 "managed" 优先级队列,在某种意义上 std::priority_queue
的大小上限为特定值 sz
.
- 我没有使用
std::priority_queue
的默认 std::less
排序,而是使用大于比较器,以便可以根据 std::priority_queue::top()
检查候选值以快速确定它们是否应该 (a) 被丢弃,(b) 替换队列中的当前最小值,或 (c) 更新队列中现有值之一的计数。这是唯一可能的,因为队列的大小被限制为 <= sz
。
我有一个数值向量 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
获得唯一向量、table 或(完整)排序都不是一个好主意,因为与 v 的长度(可能很容易 >1e6)相比,n 通常较小。
目前的解决方案是:
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)!
我敢打赌 data.table
很有竞争力:
library(data.table)
data <- data.table(v)
data[ , .N, keyby = v][(.N - n + 1):.N]
其中 n
是您要获取的号码
注意:以前的版本为 table()
而不是目标复制了功能。此版本已被删除,将在场外提供。
绘制攻击计划
下面是使用 map
.
C++98
首先,我们需要找到数字向量的 "unique" 值。
为此,我们选择将被计数为 key
的数字存储在 std::map
中,并在每次观察到该数字时递增 value
。
使用 std::map
的排序结构,我们知道前 n
个数字在 std::map
的后面。因此,我们使用迭代器弹出这些元素并将它们导出到数组中。
C++11
如果可以访问 C++11 编译器,另一种方法是使用 std::unordered_map
,它具有 O(1)
的大 O 用于插入和检索(O(n)
如果坏散列)与 std::map
相比,后者具有 O(log(n))
的大 O。
要获得正确的顶部 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 尝试的牺牲品。
我想用另一个基于 Rcpp 的解决方案来证明自己的能力,它比 DescTools
方法快 ~7 倍,比 data.table
方法快 ~13 倍,使用上面的 1e5-length 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 版本相比)。我不会 post 这个版本,因为它大部分是相同的,除了对 std::move
和 emplace
的几次调用,但是 here's a link to it。
针对前三个函数进行测试,并使用 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
其中 top_n2
是 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
class 中发生了很多事情,但只是触及一些关键点:
paired
类型本质上是 class 围绕std::pair<double, unsigned int>
的包装器,它将值与计数相关联,提供一些便利功能,例如operator++()
/operator++(int)
用于直接 pre-/post-increment 的计数,并修改了比较运算符。-
histogram
class 包装了一种 "managed" 优先级队列,在某种意义上std::priority_queue
的大小上限为特定值sz
. - 我没有使用
std::priority_queue
的默认std::less
排序,而是使用大于比较器,以便可以根据std::priority_queue::top()
检查候选值以快速确定它们是否应该 (a) 被丢弃,(b) 替换队列中的当前最小值,或 (c) 更新队列中现有值之一的计数。这是唯一可能的,因为队列的大小被限制为 <=sz
。