如何避免计算竞争指数的循环

Question

我必须为几个实验计算所谓的竞争指数。我知道物体的位置及其大小。我想计算某个半径内的大小总和以及到该半径内对象的距离总和。数据示例在这里：

set.seed(13181938)
    df <- data.frame(exp = rep(LETTERS[1:20], each = 100), x = rnorm(1000, 100, 50), 
                     y = rnorm(1000, 100, 50), di = rnorm(5, 2, 2))
df$comp1 <- 0
df$dist <- 0

我使用循环进行计算，但完成 1000 个对象的计算需要很长时间。在真实数据集中，我有超过 10000 个对象。

fori <- function(x) {
  for (i in 1:nrow(x)){
    for (j in 1:nrow(x)){
      dist = sqrt((x$x[j] - x$x[i])^2 + (x$y[j] - x$y[i])^2)
        #print(paste(x$exp[i], x$exp[j], dist))
        if(dist < 2 & x$exp[i] == x$exp[j]){
        x$comp1[i] = x$comp1[i] + x$di[j]
        x$dist[i] = x$dist[i] + dist
      }
    }
  }
  df <- data.frame(x)
  return(df)
}

abc <- fori(df)

为示例运行此循环需要很长时间，这意味着整个数据集将花费更多时间。你能建议任何其他方式吗？我尝试了apply 和DT，但没有成功。

Answer 1

像这样的循环是使用 Rcpp 加速的完美候选。逻辑转换不变：

library(Rcpp)

cppFunction('
List
computeIndex(const NumericVector x,
             const NumericVector y, 
             const NumericVector di,
             const CharacterVector ex)
{
    int n = x.size();
    NumericVector comp1(n), dist(n);

    for(int i = 0; i < n; ++i)
    {
        for(int j = 0; j < n; ++j)
        {
            double dx = x[j] - x[i], dy = y[j] - y[i];
            double d = std::sqrt(dx*dx + dy*dy);

            if((d < 2) && (ex[i] == ex[j]))
            {
                comp1[i] += di[j];
                dist[i] +=  d;
            }
        }
    }

    return List::create(Named("comp1") = comp1,
                        Named("dist") = dist);
}
')

res <- data.frame(computeIndex(df$x, df$y, df$di, df$exp))

这不仅比等效的纯 R 代码更快，而且避免了 分配任何 O(N^2) 对象。您还可以将其与 dplyr 结合使用，以避免在具有不同 exp 值的行之间进行不必要的比较：

df %>%
    group_by(exp) %>%
    do({
        res <- computeIndex(.$x, .$y, .$di, .$exp)
        data.frame(., res)
    })

Answer 2

我使用dplyr 并加入exp。然后summarise 对应每个（生成的）id。

res <- df %>% mutate(id = row_number()) %>%
  merge(df, by='exp') %>% 
  mutate(dist = sqrt((x.x - x.y)^2 + (y.x - y.y)^2)) %>% 
  filter(dist < 2 ) %>%
  group_by(id,x.x,y.x,di.x) %>%
  summarise(comp1 = sum(di.y),
                      dist = sum(dist))

结果：

Source: local data frame [2,000 x 6]
Groups: id, x.x, y.x [?]

      id       x.x       y.x       di.x      comp1     dist
   <int>     <dbl>     <dbl>      <dbl>      <dbl>    <dbl>
1      1 127.36166  89.64637 -0.2508979 -0.2508979 0.000000
2      2  90.98491 153.17911  1.4561061  1.4561061 0.000000
3      3  58.96620 144.72710  2.7909274  2.7909274 0.000000
4      4 162.44443 132.35379  3.0175213  3.0175213 0.000000
5      5 184.52673  47.12997  1.1127618  1.1127618 0.000000
6      6  57.07334 126.03554 -0.2508979 -0.2508979 0.000000
7      7  22.28946 110.69319  1.4561061  2.5688679 1.267998
8      8  40.54007 123.32645  2.7909274  2.7909274 0.000000
9      9 179.37667  61.45213  3.0175213  3.0175213 0.000000
10    10  73.82714  67.86194  1.1127618  1.1127618 0.000000
# ... with 1,990 more rows

PS：查看标准if(dist < 2 & x$exp[i] == x$exp[j])意味着只有几行符合dist