按组的地理距离-在每对行上应用函数

Question

我想计算每个省的房屋数量之间的平均地理距离。

假设我有以下数据。

df1 <- data.frame(province = c(1, 1, 1, 2, 2, 2),
              house = c(1, 2, 3, 4, 5, 6),
              lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7), 
              lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2))

使用geosphere图书馆，我可以找到两所房子之间的距离。例如：

library(geosphere)
distm(c(df1$lon[1], df1$lat[1]), c(df1$lon[2], df1$lat[2]), fun = distHaversine)

#11429.1

如何计算该省所有房屋之间的距离并收集每个省的平均距离？

原始数据集每个省都有数百万个观测值，因此此处的性能也是一个问题。

Answer 1

我最初的想法是查看的源代码，distHaversine并将其复制到要与一起使用的函数中proxy。可以这样工作（请注意应该lon是第一列）：

library(geosphere)
library(dplyr)
library(proxy)

df1 <- data.frame(province = as.integer(c(1, 1, 1, 2, 2, 2)),
                  house = as.integer(c(1, 2, 3, 4, 5, 6)),
                  lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7), 
                  lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2))

custom_haversine <- function(x, y) {
  toRad <- pi / 180

  diff <- (y - x) * toRad
  dLon <- diff[1L]
  dLat <- diff[2L]

  a <- sin(dLat / 2) ^ 2 + cos(x[2L] * toRad) * cos(y[2L] * toRad) * sin(dLon / 2) ^ 2
  a <- min(a, 1)
  # return
  2 * atan2(sqrt(a), sqrt(1 - a)) * 6378137
}

pr_DB$set_entry(FUN=custom_haversine, names="haversine", loop=TRUE, distance=TRUE)

average_dist <- df1 %>%
  select(-house) %>%
  group_by(province) %>%
  group_map(~ data.frame(avg=mean(proxy::dist(.x[ , c("lon", "lat")], method="haversine"))))

但是，如果您期望每个省有数百万行，则 proxy可能无法分配中间（矩阵的下三角）矩阵。因此，我将代码移植到C ++并添加了多线程作为奖励：

编辑：事实证明，该s2d帮助程序远非最佳，此版本现在使用此处给出的公式。

EDIT2：我刚刚发现了RcppThread，它可以用来检测用户中断。

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel,RcppThread)]]

#include <cstddef> // size_t
#include <math.h> // sin, cos, sqrt, atan2, pow
#include <vector>

#include <RcppThread.h>
#include <Rcpp.h>
#include <RcppParallel.h>

using namespace std;
using namespace Rcpp;
using namespace RcppParallel;

// single to double indices for lower triangular of matrices without diagonal
void s2d(const size_t id, const size_t nrow, size_t& i, size_t& j) {
  j = nrow - 2 - static_cast<size_t>(sqrt(-8 * id + 4 * nrow * (nrow - 1) - 7) / 2 - 0.5);
  i = id + j + 1 - nrow * (nrow - 1) / 2 + (nrow - j) * ((nrow - j) - 1) / 2;
}

class HaversineCalculator : public Worker
{
public:
  HaversineCalculator(const NumericVector& lon,
                      const NumericVector& lat,
                      double& avg,
                      const int n)
    : lon_(lon)
    , lat_(lat)
    , avg_(avg)
    , n_(n)
    , cos_lat_(lon.length())
  {
    // terms for distance calculation
    for (size_t i = 0; i < cos_lat_.size(); i++) {
      cos_lat_[i] = cos(lat_[i] * 3.1415926535897 / 180);
    }
  }

  void operator()(size_t begin, size_t end) {
    // for Kahan summation
    double sum = 0;
    double c = 0;

    double to_rad = 3.1415926535897 / 180;

    size_t i, j;
    for (size_t ind = begin; ind < end; ind++) {
      if (RcppThread::isInterrupted(ind % static_cast<int>(1e5) == 0)) return;

      s2d(ind, lon_.length(), i, j);

      // haversine distance
      double d_lon = (lon_[j] - lon_[i]) * to_rad;
      double d_lat = (lat_[j] - lat_[i]) * to_rad;
      double d_hav = pow(sin(d_lat / 2), 2) + cos_lat_[i] * cos_lat_[j] * pow(sin(d_lon / 2), 2);
      if (d_hav > 1) d_hav = 1;
      d_hav = 2 * atan2(sqrt(d_hav), sqrt(1 - d_hav)) * 6378137;

      // the average part
      d_hav /= n_;

      // Kahan sum step
      double y = d_hav - c;
      double t = sum + y;
      c = (t - sum) - y;
      sum = t;
    }

    mutex_.lock();
    avg_ += sum;
    mutex_.unlock();
  }

private:
  const RVector<double> lon_;
  const RVector<double> lat_;
  double& avg_;
  const int n_;
  tthread::mutex mutex_;
  vector<double> cos_lat_;
};

// [[Rcpp::export]]
double avg_haversine(const DataFrame& input, const int nthreads) {
  NumericVector lon = input["lon"];
  NumericVector lat = input["lat"];

  double avg = 0;
  int size = lon.length() * (lon.length() - 1) / 2;
  HaversineCalculator hc(lon, lat, avg, size);

  int grain = size / nthreads / 10;
  RcppParallel::parallelFor(0, size, hc, grain);
  RcppThread::checkUserInterrupt();

  return avg;
}

这段代码不会分配任何中间矩阵，它只会简单地计算出每对下三角形的距离，并最终累加平均值。有关Kahan求和部分，请参见此处。

如果您将该代码保存在中haversine.cpp，则可以执行以下操作：

library(dplyr)
library(Rcpp)
library(RcppParallel)
library(RcppThread)

sourceCpp("haversine.cpp")

df1 %>%
  group_by(province) %>%
  group_map(~ data.frame(avg=avg_haversine(.x, parallel::detectCores())))
# A tibble: 2 x 2
# Groups:   province [2]
  province     avg
     <int>   <dbl>
1        1  15379.
2        2 793612.

这也是一个健全性检查：

pr_DB$set_entry(FUN=geosphere::distHaversine, names="distHaversine", loop=TRUE, distance=TRUE)

df1 %>%
  select(-house) %>%
  group_by(province) %>%
  group_map(~ data.frame(avg=mean(proxy::dist(.x[ , c("lon", "lat")], method="distHaversine"))))

不过请注意：

df <- data.frame(lon=runif(1e3, -90, 90), lat=runif(1e3, -90, 90))

system.time(proxy::dist(df, method="distHaversine"))
   user  system elapsed 
 34.353   0.005  34.394

system.time(proxy::dist(df, method="haversine"))
   user  system elapsed 
  0.789   0.020   0.809

system.time(avg_haversine(df, 4L))
   user  system elapsed 
  0.054   0.000   0.014

df <- data.frame(lon=runif(1e5, -90, 90), lat=runif(1e5, -90, 90))

system.time(avg_haversine(df, 4L))
   user  system elapsed 
 73.861   0.238  19.670

如果您有数百万行，则可能需要等待一段时间。

我还应该提到，不可能在通过创建的线程中检测到用户中断RcppParallel，因此，如果您开始计算，则应该等到计算完成或完全重新启动R / RStudio。 请参阅上面的EDIT2。

关于复杂性

根据您的实际数据和计算机有多少个内核，您很可能最终需要等待数天才能完成计算。这个问题具有二次复杂性（可以说每个省）。这行：

int size = lon.length() * (lon.length() - 1) / 2;

表示必须执行的（haversine）距离计算量。因此，粗略地说，如果行数增加了倍n，则计算数量增加了倍n^2 / 2。

没有办法对此进行优化。您必须N先实际计算每个数字，然后才能计算数字的平均值，并且要比多线程C ++代码更快地找到一些东西，所以您要么必须等待它，要么扔掉更多的内核这个问题，无论是单台机器还是多台机器一起工作。否则，您将无法解决此问题。

Answer 2

鉴于您的数据具有数百万行，这听起来像是“ XY”问题。即，您真正需要的答案不是您所提出问题的答案。

让我举个比喻：如果您想知道森林中树木的平均高度，则不必测量每棵树木。您只需测量足够大的样本，以确保您的估计有足够高的概率接近所需的真实平均值。

使用每个房屋到其他房屋的距离进行暴力计算，不仅会占用过多的资源（即使使用优化的代码），而且还会提供比您可能需要的小数位数更多的数据，或者由数据准确性来证明（GPS坐标通常最多只能校正到几米以内）。

因此，我建议您对样本量进行计算，样本量仅与问题所需的准确性水平所需的大小一样大。例如，以下内容将在短短几秒钟内提供200万行的估计值，该估计值对4个有效数字有好处。您可以通过增加样本大小来提高准确性，但是鉴于GPS坐标本身存在不确定性，我怀疑这样做是否必要。

sample.size=1e6    
lapply(split(df1[3:4], df1$province), 
  function(x) {
    s1 = x[sample(nrow(x), sample.size, T), ]
    s2 = x[sample(nrow(x), sample.size, T), ]
    mean(distHaversine(s1, s2))
  })

一些大数据可以测试：

N=1e6
df1 <- data.frame(
  province = c(rep(1,N),rep(2,N)),
  house = 1:(2*N),
  lat = c(rnorm(N,-76), rnorm(N,-85)), 
  lon = c(rnorm(N,39), rnorm(N,-55,2)))

为了了解这种方法的准确性，我们可以使用自举。对于以下演示，我仅使用100,000行数据，以便我们可以在短时间内执行1000次引导程序迭代：

N=1e5
df1 <- data.frame(lat = rnorm(N,-76,0.1), lon = rnorm(N,39,0.1))

dist.f = function(i) {
    s1 = df1[sample(N, replace = T), ]
    s2 = df1[sample(N, replace = T), ]
    mean(distHaversine(s1, s2))
    }

boot.dist = sapply(1:1000, dist.f)
mean(boot.dist)
# [1] 17580.63
sd(boot.dist)
# [1] 29.39302

hist(boot.dist, 20) 

即，对于这些测试数据，平均距离为17,580 +/- 29 m。那是0.1％的变异系数，对于大多数目的而言，它可能足够准确。正如我所说，如果确实需要，可以通过增加样本数量来提高准确性。