计算具有不同属性的范围内的点

Calcualte points in range with with different properties

我有以下数据框,其中包含来自不同样本的点。每个点都有一个类型。 对于属于给定类型的给定样本的每个点(例如“Sample_1”类型“A”),我需要计算在给定截止点周围有多少其他类型的点。 我当前的实现使用“future.apply”,我想知道是否有更有效的方法来解决这个问题。这里的例子是有限的,应该 运行 很快,真正的问题是由几千行组成的,而且速度要慢得多。 最后,我将结果存储在一个列表中。 对于“starting_point”中具有“type”的每个元素,此列表具有阈值为 40 的“target_point”类型元素的数量。

library(future)
library(future.apply)

a_test=data.frame(ID=sample(c("Sample_1", "Sample_2", "Sample_3"), 100, replace=TRUE), type=sample(c("A", "B", "C", "D"), 100, replace = TRUE), xpos=sample(1:200, 100, replace=TRUE), ypos=sample(1:200, 100, replace=TRUE))

starting_point=c("A", "B")
target_point=c("C", "D")
threshold=40
result_per_pair=list()

for(sp in starting_point){
  ## Here I select a data frame of "Starting points" without looking
  ## from which ID they came from
  sp_tdf=a_test[a_test$type==sp, ]
  for(tp in target_point){
    ## Here I select a data frame of "Target points" without looking
    ## from which ID they came from
    tp_tdf=a_test[a_test$type==tp, ]
    ## I use future_sapply here, parallelizing on each line of "sp_tdf"
    plan(multisession)
    elements_around=future_sapply(1:nrow(sp_tdf), function(x, sp_tdf, tp_tdf, treshold2){
      xc=sp_tdf$xpos[x]
      yc=sp_tdf$ypos[x]
      ###  NOTE HERE:  At this point I select the points that are in the same
      ###  ID as the current line of sp_tdf
      tp_tdf2=tp_tdf[tp_tdf$ID == sp_tdf$ID[x],]
      ares=tp_tdf2[ (tp_tdf2$xpos-xc)^2 + (tp_tdf2$ypos-yc)^2 <threshold2, ]
      return(nrow(ares))
    },sp_tdf=sp_tdf, tp_tdf=tp_tdf, threshold2=threshold*threshold)
    
    a_newcol=paste0(tp, "_around_", sp)
    ## we need to create a copy of sp_tdf otherwise we add columns to the
    ## initial sp_tdf and we memorize them in the wrong place in the list
    sp_tdf_temp=sp_tdf
    sp_tdf_temp[,  a_newcol]=elements_around
    
    result_per_pair[[ paste0(tp, "_around_", sp ) ]]=rbind(result_per_pair[[ paste0(tp, "_around_", sp ) ]], sp_tdf_temp)
  }
}

可以看到table的类型我这里是:

head(result_per_pair[[1]])

$C_around_A
          ID type xpos ypos C_around_A
1   Sample_2    A   26   74          1
2   Sample_3    A   64    8          1
3   Sample_3    A  121    2          1
5   Sample_2    A   62   94          0

您可以尝试使用 RANN::nn2 函数:

id_list <- split(a_test, a_test$ID) 

res <- id_list %>%
  map(~select(.x, xpos, ypos)) %>%
  map(~RANN::nn2(.x, .x, k = nrow(.), searchtype = "radius", radius = threshold)) %>%
  map(1) %>%
  map2(
    id_list, 
    function(x, y){ 
      seq_len(nrow(x)) %>% 
        map(~x[.x,] %>% .[. > 0]) %>% 
        map(~y[.x,]) %>% 
        map("type") %>%
        map_dfr(table) %>%
        mutate(across(everything(), as.integer))
    }
  ) %>%
  map2_dfr(id_list, ~bind_cols(.y, .x))

替换 tidyverse 函数可能会进行一些时间改进(很难说它在您的示例中有多快)。结果:

res %>% head()

ID type xpos ypos A B C D
Sample_1    C   48  157 0 0 3 1
Sample_1    D  177   97 1 1 1 3
Sample_1    C   10   71 0 0 3 0
Sample_1    C   71  168 1 1 2 0
Sample_1    D   82   48 1 0 1 2
Sample_1    C  165   71 3 3 1 1

其中 A-D 列表示同一 ID 中的类型数。我使用种子 123 生成 a_test。您可以调整算法以与 starting_pointtarget_point 一起使用,将每个 id_list 分成两部分 - 由 starting_pointtarget_point 定义的部分并调整 data & query RANN::nn2 中的参数。

编辑

基于楼上评论思路的功能:

f <- function(df, threshold, start = levels(df$type), target = levels(df$type)){
  
  my_lists <- df %>%
    filter(type %in% c(start, target)) %>%
    split(.$ID) %>%
    map(
      function(x){ 
        map(
          list(start, target), 
          ~filter(x, type %in% .x) %>% mutate(type = droplevels(type))
        )
      }
    ) %>%
    discard(~any(map_int(.x, nrow) == 0))  
    
  indices <- my_lists %>%
    map(
      ~RANN::nn2(
        data = select(.x[[2]], xpos, ypos), 
        query = select(.x[[1]], xpos, ypos), 
        k = nrow(.x[[2]]), 
        searchtype = "radius", 
        radius = threshold
      )
    ) %>%
    map(1) %>%
    map(function(x) seq_len(nrow(x)) %>% map(~x[.x,] %>% .[. > 0]))
  
  my_lists %>%
    map(2) %>%
    map2(indices, function(x, y) map_dfr(y, ~summary(x[.x,]$type))) %>%
    {map2_dfr(map(my_lists, 1), ., bind_cols)}
}

以半径 40 围绕 A 计算 C:

f(a_test, 40, "A", "C")