look into using geodist for distance calculations

[njtierney]

https://cran.r-project.org/web/packages/geosphere/index.html

[mpadge]

Oh yeah, sorry, just wrote what i did on #51, and hadn't yet seen this one ...

[njtierney]

geodist is now on CRAN (:tada:), so another good reason to get looking at applying it!

I think that I could probably just get away with use geodist instead of: distance_matrix_cpp, which then means I can remove spherical_distance_cpp and binary_distance_matrix.

This then leaves nearest_facility_dist, but perhaps this could just take the output of geodist on the facilities and users.

[mpadge]

some motivation for you here:

devtools::load_all ("<...>/maxcovr", export_all = TRUE)
#> Loading maxcovr
suppressMessages (library (dplyr))
york_selected <- york %>% filter(grade == "I")
york_unselected <- york %>% filter(grade != "I")
distance_cutoff <- 100
user <- tibble::rowid_to_column(york_crime, var = "user_id")
user_not_covered <- find_users_not_covered(york_selected,
                                           user,
                                           distance_cutoff)
proposed_facility <- york_unselected

knitr::kable (rbenchmark::benchmark (
    A <- binary_distance_matrix(facility = proposed_facility,
                                user = user_not_covered,
                                distance_cutoff = distance_cutoff),
    d1 <- which (geodist::geodist (user_not_covered, proposed_facility, measure = "cheap") < distance_cutoff),
    d2 <- which (geodist::geodist (user_not_covered, proposed_facility, measure = "haversine") < distance_cutoff),
    d3 <- which (geodist::geodist (user_not_covered, proposed_facility, measure = "geodesic") < distance_cutoff),
replications = 2, order = "relative"))

	test	replications	elapsed	relative	user.self	sys.self	user.child	sys.child
2	d1 <- which(geodist::geodist(user_not_covered, proposed_facility, measure = “cheap”) < distance_cutoff)	2	0.313	1.000	0.252	0.059	0	0
3	d2 <- which(geodist::geodist(user_not_covered, proposed_facility, measure = “haversine”) < distance_cutoff)	2	0.440	1.406	0.379	0.060	0	0
1	A <- binary_distance_matrix(facility = proposed_facility, user = user_not_covered, distance_cutoff = distance_cutoff)	2	0.644	2.058	0.584	0.060	0	0
4	d3 <- which(geodist::geodist(user_not_covered, proposed_facility, measure = “geodesic”) < distance_cutoff)	2	6.716	21.457	6.657	0.053	0	0

^{Created on 2019-06-05 by the reprex package (v0.3.0)}

Your C++ code is not too bad, but geodist is still faster, and offers the distinct advantage of multiple measures, including nanometre-accuracy geodesic, which some folk may consider important. This would simply involve replacing your binary_distance_matrix function with the code used above of:

which (geodist::geodist (...) < distance_cutoff)

Shall i PR that too?

[mpadge]

Note also that your radius is not the "official" WSG84 standard as in geodist, and most other places, and so geodist(..., measure = "haversine") does not give exactly the same results as your code. But just switching to geodist will resolve that anyway ...

[njtierney]

Hi @mpadge !

That is indeed great motivation!

I would love a PR to implement this - I think that the unit tests should help make sure not too much else breaks, I think that one problem may be that the format of geodist() and binary_distance_matrix() , so maybe there might need to be some syntactic sugar around that?

[mpadge]

I'll gladly get on to it, but wait for this PR first, noting especially that travis is currently failing because neither gurobi nor Rglpk can be installed. Once you've fixed that, then I'll get moving. :+1:

[njtierney]

Yeah so I'm not sure what to do about gurobi - it's an internal package from gurobi that isn't on CRAN or GitHub. I'll make another issue for that.