0

I am trying to replicate the following code in Rcpp (original pandas source from the following link- 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

I created a

haversine.cpp

file in my home directory as follows:

#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)]
})
*/

and called it using sourceCpp as

 sourceCpp('./haversine.cpp')

In my opinion, there is something wrong with the for loop which is causing it to crash, but I just cant seem to find out what it is. I say this cause when I did a dry run without a loop and only a single element of the vector at index 0, the rcpp function ran. The only link I have found useful is where the for loop was not written correctly (Rcpp function crashes), but somehow I have tried everything it said and still cant find out the crash reason. Please help!

Gompu
  • 415
  • 1
  • 6
  • 21
  • 1
    maybe [my implementation](https://github.com/SymbolixAU/spatialdatatable/blob/master/src/DistanceCalculations.cpp#L315) can help you find your error – SymbolixAU Nov 28 '17 at 22:39
  • I also [wrote a blog](https://www.symbolix.com.au/blog-main/j26ynx6awfl32brcaxdjxnwalwmmbc) which I think outlines the approach you're taking. Maybe this could be of assistance too. – SymbolixAU Nov 29 '17 at 00:20
  • Thanks SymbolixAU, I will take a look at your work and add anything interesting I might find in your solutions. The solution posted by @hunzikp works for me. – Gompu Nov 29 '17 at 15:30

1 Answers1

1

Your session crashes because you create NumericVector objects of length zero and then try to assign them values using the unsafe bracket ([i]) notation. If you initialize the NumericVectors with correct length your code runs (I haven't checked its accuracy, though):

#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;
}

On a more general note: Using the safer .at(i) method makes your code fail more gracefully.

hunzikp
  • 131
  • 4