r - 如何有效地识别重复的有序对



>我正在使用在单独的列中包含起点 ID 和目的地 ID 的起点-目的地 (OD) 数据。有时,聚合相同的 OD 对很重要,只是起点和目的地是交换的。

OD 数据如下所示:

orign             dest      value
E02002361         E02002361 109
E02002361         E02002363  38
E02002361         E02002367  10
E02002361         E02002371  44
E02002363         E02002361  34

在上面的示例中,第一行和最后一行可以被视为同一对,只是方向相反。挑战在于如何有效地识别它们是重复的。

我创建了一个包 stplanr,它可以回答这个问题,如下面的可重现示例所示:

x = read.csv(stringsAsFactors = FALSE, text = "orign,dest,value
E02002361,E02002361,109
E02002361,E02002363,38
E02002361,E02002367,10
E02002361,E02002371,44
E02002363,E02002361,34")
duplicated(stplanr::od_id_order(x)[[3]])
#> Registered S3 method overwritten by 'R.oo':
#>   method        from       
#>   throw.default R.methodsS3
#> [1] FALSE FALSE FALSE FALSE  TRUE

创建于 2019-07-27 由 reprex 软件包 (v0.3.0)

这种方法的问题在于它对于大型数据集来说很慢。

我已经研究了从矩阵中的每一列中获取 Min 的最快方法?这表明pmin是跨多列(而不是 2)获取最小值的最有效方法,我们已经在使用它。

与删除重复组合(无论顺序如何)不同,这个问题是关于仅 2 列的重复识别和效率。删除重复组合(无论顺序如何)中发布的解决方案似乎比以下计时中显示的最慢解决方案慢。

比这更快的解决方案是szudzik_pairing函数,由我的同事马尔科姆·摩根创建,基于 Matthew Szudzik 开发的方法。

我们已经尝试了每种方法,Szudzik 方法确实看起来更快,但我想知道:有没有更有效的方法(在任何语言中,但最好可以在 R 中实现)?

以下是我们所做的一个快速可重现的示例,包括一个显示计时的简单基准测试:

od_id_order_base <- function(x, id1 = names(x)[1], id2 = names(x)[2]) {
data.frame(
stringsAsFactors = FALSE,
stplanr.id1 = x[[id1]],
stplanr.id1 = x[[id2]],
stplanr.key = paste(
pmin(x[[id1]], x[[id2]]),
pmax(x[[id1]], x[[id2]])
)
)
}
od_id_order_rfast <- function(x, id1 = names(x)[1], id2 = names(x)[2]) {
data.frame(
stringsAsFactors = FALSE,
stplanr.id1 = x[[id1]],
stplanr.id1 = x[[id2]],
stplanr.key = paste(
Rfast::colPmin(as.numeric(as.factor(x[[id1]])), as.numeric(as.factor(x[[id1]]))),
pmax(x[[id1]], x[[id2]])
)
)
}
od_id_order_dplyr <- function(x, id1 = names(x)[1], id2 = names(x)[2]) {
dplyr::transmute_(x,
stplanr.id1 = as.name(id1),
stplanr.id2 = as.name(id2),
stplanr.key = ~paste(pmin(stplanr.id1, stplanr.id2), pmax(stplanr.id1, stplanr.id2))
)
}
szudzik_pairing <- function(val1, val2, ordermatters = FALSE) {
if(length(val1) != length(val2)){
stop("val1 and val2 are not of equal length")
}
if(class(val1) == "factor"){
val1 <- as.character(val1)
}
if(class(val2) == "factor"){
val2 <- as.character(val2)
}
lvls <- unique(c(val1, val2))
val1 <- as.integer(factor(val1, levels = lvls))
val2 <- as.integer(factor(val2, levels = lvls))
if(ordermatters){
ismax <- val1 > val2
stplanr.key <- (ismax * 1) * (val1^2 + val1 + val2) + ((!ismax) * 1) * (val2^2 + val1)
}else{
a <- ifelse(val1 > val2, val2, val1)
b <- ifelse(val1 > val2, val1, val2)
stplanr.key <- b^2 + a
}
return(stplanr.key)
}

n = 1000
ids <- as.character(runif(n, 1e4, 1e7 - 1))
x <- data.frame(id1 = rep(ids, times = n),
id2 = rep(ids, each = n),
val = 1,
stringsAsFactors = FALSE)

head(od_id_order_base(x))
#>        stplanr.id1    stplanr.id1.1                       stplanr.key
#> 1 8515501.50763425 8515501.50763425 8515501.50763425 8515501.50763425
#> 2 2454738.52108038 8515501.50763425 2454738.52108038 8515501.50763425
#> 3  223811.25236322 8515501.50763425  223811.25236322 8515501.50763425
#> 4 4882305.41496906 8515501.50763425 4882305.41496906 8515501.50763425
#> 5  4663684.5752892 8515501.50763425  4663684.5752892 8515501.50763425
#> 6 725621.968830239 8515501.50763425 725621.968830239 8515501.50763425
head(od_id_order_rfast(x))
#>        stplanr.id1    stplanr.id1.1          stplanr.key
#> 1 8515501.50763425 8515501.50763425 830 8515501.50763425
#> 2 2454738.52108038 8515501.50763425 163 8515501.50763425
#> 3  223811.25236322 8515501.50763425 135 8515501.50763425
#> 4 4882305.41496906 8515501.50763425 435 8515501.50763425
#> 5  4663684.5752892 8515501.50763425 408 8515501.50763425
#> 6 725621.968830239 8515501.50763425 689 8515501.50763425
head(od_id_order_dplyr(x))
#> Warning: transmute_() is deprecated. 
#> Please use transmute() instead
#> 
#> The 'programming' vignette or the tidyeval book can help you
#> to program with transmute() : https://tidyeval.tidyverse.org
#> This warning is displayed once per session.
#>        stplanr.id1      stplanr.id2                       stplanr.key
#> 1 8515501.50763425 8515501.50763425 8515501.50763425 8515501.50763425
#> 2 2454738.52108038 8515501.50763425 2454738.52108038 8515501.50763425
#> 3  223811.25236322 8515501.50763425  223811.25236322 8515501.50763425
#> 4 4882305.41496906 8515501.50763425 4882305.41496906 8515501.50763425
#> 5  4663684.5752892 8515501.50763425  4663684.5752892 8515501.50763425
#> 6 725621.968830239 8515501.50763425 725621.968830239 8515501.50763425
head(szudzik_pairing(x$id1, x$id2))
#> [1]  2  5 10 17 26 37

system.time(od_id_order_base(x))
#>    user  system elapsed 
#>   0.467   0.000   0.467
system.time(od_id_order_rfast(x))
#>    user  system elapsed 
#>   1.063   0.001   1.064
system.time(od_id_order_dplyr(x))
#>    user  system elapsed 
#>   0.493   0.000   0.493
system.time(szudzik_pairing(x$id1, x$id2))
#>    user  system elapsed 
#>   0.100   0.000   0.101

创建于 2019-07-27 由 reprex 软件包 (v0.3.0)

devtools::session_info()
#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.6.0 (2019-04-26)
#>  os       Debian GNU/Linux 9 (stretch)
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       Etc/UTC                     
#>  date     2019-07-27                  
#> 
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package      * version date       lib source        
#>  assertthat     0.2.1   2019-03-21 [1] CRAN (R 3.6.0)
#>  backports      1.1.4   2019-04-10 [1] CRAN (R 3.6.0)
#>  callr          3.3.0   2019-07-04 [1] CRAN (R 3.6.0)
#>  cli            1.1.0   2019-03-19 [1] CRAN (R 3.6.0)
#>  crayon         1.3.4   2017-09-16 [1] CRAN (R 3.6.0)
#>  desc           1.2.0   2018-05-01 [1] CRAN (R 3.6.0)
#>  devtools       2.1.0   2019-07-06 [1] CRAN (R 3.6.0)
#>  digest         0.6.20  2019-07-04 [1] CRAN (R 3.6.0)
#>  dplyr          0.8.3   2019-07-04 [1] CRAN (R 3.6.0)
#>  evaluate       0.14    2019-05-28 [1] CRAN (R 3.6.0)
#>  fs             1.3.1   2019-05-06 [1] CRAN (R 3.6.0)
#>  glue           1.3.1   2019-03-12 [1] CRAN (R 3.6.0)
#>  highr          0.8     2019-03-20 [1] CRAN (R 3.6.0)
#>  htmltools      0.3.6   2017-04-28 [1] CRAN (R 3.6.0)
#>  knitr          1.23    2019-05-18 [1] CRAN (R 3.6.0)
#>  magrittr       1.5     2014-11-22 [1] CRAN (R 3.6.0)
#>  memoise        1.1.0   2017-04-21 [1] CRAN (R 3.6.0)
#>  pillar         1.4.2   2019-06-29 [1] CRAN (R 3.6.0)
#>  pkgbuild       1.0.3   2019-03-20 [1] CRAN (R 3.6.0)
#>  pkgconfig      2.0.2   2018-08-16 [1] CRAN (R 3.6.0)
#>  pkgload        1.0.2   2018-10-29 [1] CRAN (R 3.6.0)
#>  prettyunits    1.0.2   2015-07-13 [1] CRAN (R 3.6.0)
#>  processx       3.4.0   2019-07-03 [1] CRAN (R 3.6.0)
#>  ps             1.3.0   2018-12-21 [1] CRAN (R 3.6.0)
#>  purrr          0.3.2   2019-03-15 [1] CRAN (R 3.6.0)
#>  R6             2.4.0   2019-02-14 [1] CRAN (R 3.6.0)
#>  Rcpp           1.0.1   2019-03-17 [1] CRAN (R 3.6.0)
#>  RcppZiggurat   0.1.5   2018-06-10 [1] CRAN (R 3.6.0)
#>  remotes        2.1.0   2019-06-24 [1] CRAN (R 3.6.0)
#>  Rfast          1.9.4   2019-05-25 [1] CRAN (R 3.6.0)
#>  rlang          0.4.0   2019-06-25 [1] CRAN (R 3.6.0)
#>  rmarkdown      1.13    2019-05-22 [1] CRAN (R 3.6.0)
#>  rprojroot      1.3-2   2018-01-03 [1] CRAN (R 3.6.0)
#>  sessioninfo    1.1.1   2018-11-05 [1] CRAN (R 3.6.0)
#>  stringi        1.4.3   2019-03-12 [1] CRAN (R 3.6.0)
#>  stringr        1.4.0   2019-02-10 [1] CRAN (R 3.6.0)
#>  testthat       2.1.1   2019-04-23 [1] CRAN (R 3.6.0)
#>  tibble         2.1.3   2019-06-06 [1] CRAN (R 3.6.0)
#>  tidyselect     0.2.5   2018-10-11 [1] CRAN (R 3.6.0)
#>  usethis        1.5.1   2019-07-04 [1] CRAN (R 3.6.0)
#>  withr          2.1.2   2018-03-15 [1] CRAN (R 3.6.0)
#>  xfun           0.8     2019-06-25 [1] CRAN (R 3.6.0)
#>  yaml           2.2.0   2018-07-25 [1] CRAN (R 3.6.0)
#> 
#> [1] /usr/local/lib/R/site-library
#> [2] /usr/local/lib/R/library

以下是从元组(x,y)生成唯一键的四种可能方法,其中顺序无关紧要:od_id_order_baseszudzik_pairing,根据OP的问题; 修改后的Szudzik方法szudzik_pairing_alt; 以及使用此答案中显示的公式的方法max_min

convert_to_numeric <- function(x, y) {
if (length(x) != length(y)) stop("x and y are not of equal length")
if (class(x) == "factor") x <- as.character(x)
if (class(y) == "factor") y <- as.character(y)
lvls <- unique(c(x, y))
x <- as.integer(factor(x, levels = lvls))
y <- as.integer(factor(y, levels = lvls))
list(x = x, y = y)  
}
od_id_order_base <- function(x, y) {
d <- convert_to_numeric(x, y)
x <- d$x
y <- d$y
paste(pmin(x, y), pmax(x, y))
}
szudzik_pairing <- function(x, y) {
d <- convert_to_numeric(x, y)
x <- d$x
y <- d$y
a <- ifelse(x > y, y, x)
b <- ifelse(x > y, x, y)
b^2 + a
}
szudzik_pairing_alt <- function(x, y) {
d <- convert_to_numeric(x, y)
x <- d$x
y <- d$y
z <- y^2 + x
ifelse(y < x, x^2 + y, z)
}
max_min <- function(x, y) {
d <- convert_to_numeric(x, y)
x <- d$x
y <- d$y
a <- pmax(x, y)
b <- pmin(x, y)
a * (a + 1) / 2 + b
}

从一些示例数据生成密钥,并验证我们是否得到相同的结果:

check_dupe <- function(f, x, y) duplicated(f(x, y))
set.seed(123)
n <- 1000^2
x <- ceiling(runif(n) * 1000)
y <- ceiling(runif(n) * 1000)
p <- lapply(list(od_id_order_base, szudzik_pairing, szudzik_pairing_alt, 
max_min), check_dupe, x, y)
all(sapply(p[-1], function(x) identical(p[[1]], as.vector(x))))
# [1] TRUE

标杆:

bmk <- microbenchmark::microbenchmark(
p1 = check_dupe(od_id_order_base, x, y),
p2 = check_dupe(szudzik_pairing, x, y),
p3 = check_dupe(szudzik_pairing_alt, x, y),
p4 = check_dupe(max_min, x, y),
times = 100
)
# Unit: seconds
#  expr      min       lq     mean   median       uq      max neval cld
#    p1 2.958512 3.089615 3.336621 3.336915 3.474973 4.721378   100   b
#    p2 1.934742 2.058185 2.191331 2.190588 2.306203 2.983729   100  a 
#    p3 1.889201 1.990306 2.173845 2.138995 2.259218 5.186751   100  a 
#    p4 1.870261 1.980756 2.143026 2.145458 2.234580 3.111324   100  a 

相关内容

  • 没有找到相关文章

最新更新