Geographical Distances Benchmarks (Crystal)

001.txt

# Benchmark run with Crystal 0.17.3 (LLVM 3.6) on "Intel Core i7-4712HQ CPU @ 2.30GHz" CPU

Brooklin (block)
haversine: 0.144098
     sqrt: 0.19007 (+31.904%)
 spheroid: 0.144098 (+1.59146e-09%)
ellipsoid: 0.188279 (+30.6604%)

Paris -> Marseille
haversine: 660.815
     sqrt: 953.726 (+44.3257%)
 spheroid: 660.958 (+0.021547%)
ellipsoid: 641.647 (+-2.90069%)

Paris -> New York
haversine: 5836.69
     sqrt: 8585.3 (+47.0919%)
 spheroid: 6093.28 (+4.39621%)
ellipsoid: 5993.13 (+2.68024%)

Paris -> Tokyo
haversine: 9712.67
     sqrt: 15400.5 (+58.5613%)
 spheroid: 11395 (+17.3213%)
ellipsoid: 2615.34 (+-73.0729%)

# M stands for Millions of Iterations per Second:
     sqrt 354.86M (± 4.72%)  1.12× slower
haversine  11.08M (± 2.69%) 35.93× slower
 spheroid  398.2M (± 5.17%)       fastest
ellipsoid 354.25M (± 5.41%)  1.12× slower

bench.cr

# All formulaes come from:
#   - https://en.wikipedia.org/wiki/Geographical_distance
#   - https://en.wikipedia.org/wiki/Haversine_formula
#
# TODO:
#   - Vincenty's Formulae: https://en.wikipedia.org/wiki/Vincenty%27s_formulae
#   - Compare distances over different regions over the globe
#   - Compare distances involving the N/S polars
# 
# Run me with `crystal run --release bench.cr`
require "benchmark"

RAD_PER_DEG = Math::PI / 180
GREAT_CIRCLE_RADIUS_METERS = 6_371_009.0

macro to_radians(p)
  {{p}}.map { |n| n * RAD_PER_DEG }
end

def world_mercator(p)
  lon, lat = to_radians(p)
  x = GREAT_CIRCLE_RADIUS_METERS * lon
  y = GREAT_CIRCLE_RADIUS_METERS * Math.log((Math.sin(lat) + 1) / Math.cos(lat))
  {x, y}
end

# Pythagore (Square Root) with World Mercator projection.
def sqrt(p, q)
  x1, y1 = world_mercator(p)
  x2, y2 = world_mercator(q)
  Math.sqrt((x1 - x2) ** 2 + (y2 - y1) ** 2)
end

macro hav(delta)
  Math.sin({{delta}} / 2) ** 2
end

# Haversine (Spherical projection plane).
def haversine(p, q)
  lon1, lat1 = to_radians(p)
  lon2, lat2 = to_radians(q)

  dlon = lon2 - lon1
  dlat = lat2 - lat1

  a = hav(dlat) + Math.cos(lat1) * Math.cos(lat2) * hav(dlon)
  d = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a))
  d * GREAT_CIRCLE_RADIUS_METERS
end

# Great Circle (Spheroid projection plane).
def spheroid(p, q)
  lon1, lat1 = to_radians(p)
  lon2, lat2 = to_radians(q)

  dlon = lon2 - lon1
  dlat = lat2 - lat1
  mlat = (lat1 + lat2) / 2

  d = Math.sqrt(dlat ** 2 + (Math.cos(mlat) * dlon) ** 2)
  d * GREAT_CIRCLE_RADIUS_METERS
end

# Ellipsoid projection plane (FCC).
def ellipsoid(p, q)
  lon1, lat1 = p
  lon2, lat2 = q

  dlon = lon2 - lon1
  dlat = lat2 - lat1
  mlat = (lat1 + lat2) / 2

  klat = 111.13209 - 0.56605 * Math.cos(2 * mlat) + 0.00120 * Math.cos(4 * mlat)
  klon = 111.41513 * Math.cos(mlat) - 0.09455 * Math.cos(3 * mlat) + 0.00012 * Math.cos(5 * mlat)

  d = Math.sqrt((klat * dlat) ** 2 + (klon * dlon) ** 2)
  d * 1000
end

def compare(p, q)
  hdist = haversine(p, q) / 1000
  rdist = sqrt(p, q) / 1000
  sdist = spheroid(p, q) / 1000
  edist = ellipsoid(p, q) / 1000
  puts "haversine: #{hdist}"
  puts "     sqrt: #{rdist} (+#{100 / hdist * rdist - 100}%)"
  puts " spheroid: #{sdist} (+#{100 / hdist * sdist - 100}%)"
  puts "ellipsoid: #{edist} (+#{100 / hdist * edist - 100}%)"
end

Brooklin1 = {-73.993683, 40.700667}
Brooklin2 = {-73.995389, 40.700586}
Marseille = {5.38,43.29}
NewYork = {-74.0, 40.717}
Paris = {2.35, 48.85}
Tokyo = {139.69,35.69}

puts "\nBrooklin (one block)"
compare(Brooklin1, Brooklin2)

puts "\nParis -> Marseille"
compare(Paris, Marseille)

puts "\nParis -> New York"
compare(Paris, NewYork)

puts "\nParis -> Tokyo"
compare(Paris, Tokyo)

puts
Benchmark.ips do |x|
  x.report("sqrt") { sqrt(Paris, NewYork) }
  x.report("haversine") { haversine(Paris, NewYork) }
  x.report("spheroid") { spheroid(Paris, NewYork) }
  x.report("ellipsoid") { ellipsoid(Paris, NewYork) }
end