Я пытаюсь скопировать следующий код в Rcpp (исходный код панды по следующей ссылке: https://engineering.upside.com/a-beginners-guide-to-optimizing-pandas-code-for-speed-c09ef2c6a4d6:
library(data.table)
library(microbenchmark)
deg2rad <- function(deg) {(deg * pi) / (180)}
haversine = function(lat1, lon1, lat2, lon2) {
MILES = 3959
lat1 = deg2rad(lat1)
lon1 = deg2rad(lon1)
lat2 = deg2rad(lat2)
lon2 = deg2rad(lon2)
dlat = lat2 - lat1
dlon = lon2 - lon1
a = sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlon/2)^2
c = 2 * asin(sqrt(a))
total_miles = MILES * c
return(total_miles)
}
# get data from here
download.file("https://raw.githubusercontent.com/sversh/pycon2017-
optimizing-pandas/master/new_york_hotels.csv","new_york_hotels.csv")
nyc_hotels = fread("new_york_hotels.csv", na.strings = c("NA", "N/A",
"NULL"))
summary(microbenchmark({
nyc_hotels[, greater_circle := haversine(40.671, -73.985, latitude,
longitude)]
},times=1000))[,-1]
# min lq mean median uq max neval
# 290.161 318.559 366.6786 329.491 345.0295 4365.697 1000
##########
#version 2 - invoke update differently, no change to function
summary(microbenchmark({
set(nyc_hotels,j="greater_circle",value=haversine(40.671, -73.985,
nyc_hotels[['latitude']], nyc_hotels[['longitude']]))
},times=1000))[,-1]
# min lq mean median uq max neval
# 81.395 89.5985 123.2211 96.1635 103.476 3670.193 1000
Я создал
haversine.cpp
файл в моем домашнем каталоге следующим образом:
#include <Rcpp.h>
#include <iostream>
using namespace Rcpp;// [[Rcpp::export]]
NumericVector haversine_cpp_fun(double lat1_cpp,double lon1_cpp,NumericVector lat2_cpp,NumericVector lon2_cpp){
double Miles = 3959.0;
int n = lat2_cpp.size();
NumericVector dlat_cpp;
NumericVector dlon_cpp;
NumericVector a_cpp;
NumericVector c_cpp;
NumericVector total_mile_cpp;
lat1_cpp = (lat1_cpp*3.14159)/180.0;
lon1_cpp = (lon1_cpp*3.14159)/180.0;
for (int i=0 ; i<n ; ++i){
lat2_cpp[i] = (lat2_cpp[i]*3.14159)/180.0;
lon2_cpp[i] = (lon2_cpp[i]*3.14159)/180.0;
dlat_cpp[i] = lat2_cpp[i] - lat1_cpp;
dlon_cpp[i] = lon2_cpp[i] - lon1_cpp;
a_cpp[i] = pow(sin(dlat_cpp[i]/2.0),2.0) + cos(lat1_cpp) * cos(lat2_cpp[i]) * pow(sin(dlon_cpp[i]/2.0),2.0);
c_cpp[i] = 2 * asin(sqrt(a_cpp[i]));
total_mile_cpp[i] = Miles * c_cpp[i];
}
return total_mile_cpp;
}
/***R
# Approach 1: Trying to use the set statement from data.table--- fails without giving error. The session just crashes
summary(microbenchmark({
set(nyc_hotels,j="greater_circle",value=haversine_cpp_fun(40.671, -73.985,
nyc_hotels[['latitude']], nyc_hotels[['longitude']]))
},times=1000))[,-1]
# Approach 2: Without using the set statement from data.table and doing thing in a simple way by a simple function call--- again fails without giving error. The R session just crashes again.
microbenchmark({
nyc_hotels[, greater_circle := haversine_cpp_fun(40.671, -73.985, latitude,
longitude)]
})
*/
и вызвал его, используя sourceCpp как
sourceCpp('./haversine.cpp')
На мой взгляд, что-то не так с циклом for, который вызывает его сбой, но я просто не могу понять, что это такое. Я говорю об этой причине, когда я выполнил пробный прогон без цикла и только одного элемента вектора с индексом 0, функция rcpp запустилась. Единственная ссылка, которую я нашел полезной, это где цикл for был написан неправильно (Сбой функции RCPP), но каким-то образом я перепробовал все, что он сказал, и до сих пор не могу выяснить причину сбоя.
Пожалуйста помоги!
Ваш сеанс завершается сбоем, потому что вы создаете объекты NumericVector нулевой длины, а затем пытаетесь присвоить им значения с помощью небезопасной скобки ([i]
) обозначение. Если вы инициализируете NumericVectors с правильной длиной, ваш код выполняется (хотя я не проверял его точность):
#include <Rcpp.h>
#include <iostream>
using namespace Rcpp;// [[Rcpp::export]]
NumericVector haversine_cpp_fun(double lat1_cpp, double lon1_cpp,
NumericVector lat2_cpp, NumericVector lon2_cpp){
double Miles = 3959.0;
int n = lat2_cpp.size();
NumericVector dlat_cpp(n);
NumericVector dlon_cpp(n);
NumericVector a_cpp(n);
NumericVector c_cpp(n);
NumericVector total_mile_cpp(n);
lat1_cpp = (lat1_cpp*3.14159)/180.0;
lon1_cpp = (lon1_cpp*3.14159)/180.0;
for (int i=0 ; i<n ; ++i){
lat2_cpp[i] = (lat2_cpp[i]*3.14159)/180.0;
lon2_cpp[i] = (lon2_cpp[i]*3.14159)/180.0;
dlat_cpp[i] = lat2_cpp[i] - lat1_cpp;
dlon_cpp[i] = lon2_cpp[i] - lon1_cpp;
a_cpp[i] = pow(sin(dlat_cpp[i]/2.0),2.0) + cos(lat1_cpp) * cos(lat2_cpp[i]) * pow(sin(dlon_cpp[i]/2.0),2.0);
c_cpp[i] = 2 * asin(sqrt(a_cpp[i]));
total_mile_cpp[i] = Miles * c_cpp[i];
}
return total_mile_cpp;
}
На более общее примечание: использование безопаснее .at(i)
метод делает ваш код более изящным.
Других решений пока нет …