gpt4 book ai didi

r - 按组划分的地理距离 - 在每对行上应用一个函数

转载 作者:行者123 更新时间:2023-12-03 15:11:53 26 4
gpt4 key购买 nike

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

假设我有以下数据。

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

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

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

最佳答案

我最初的想法是查看 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助手远非最佳,
这个版本现在使用给定的公式 here .

编辑2 : 我刚刚发现 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;
}

此代码不会分配任何中间矩阵,
它将简单地计算每对下三角形的距离,并最终累积平均值。
here对于 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++ 代码更快的东西,
所以你要么必须等待它,
或者在问题上投入更多的核心,
使用一台机器或多台机器一起工作。
否则你无法解决这个问题。

关于r - 按组划分的地理距离 - 在每对行上应用一个函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/55621827/

26 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com