7

我正在处理在不同列中包含原始 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

由reprex 包（v0.3.0）于 2019 年 7 月 27 日创建

这种方法的问题是对于大型数据集来说速度很慢。

我已经研究过从矩阵中的每一列中获取 Min 的最快方法？这表明这pmin是在多列（不是 2 列）中获取最小值的最有效方法，我们已经在使用它。

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

比这更快的解决方案是szudzik_pairing函数，它由我的同事 Malcolm Morgan 创建，基于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

由reprex 包（v0.3.0）于 2019 年 7 月 27 日创建

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
4

1 回答 1

4

以下是从元组 (x, y) 生成唯一键的四种可能方法，其中顺序无关紧要：od_id_order_base并且szudzik_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
于 2019-07-27T22:56:09.567 回答