#!/usr/bin/env python
# coding: utf-8
# _Since we announced [our collaboration with the World Bank and more partners to create the Open Traffic platform](https://mapzen.com/blog/announcing-open-traffic/), we’ve been busy. We’ve shared [two](https://mapzen.com/blog/open-traffic-osmlr-technical-preview/) [technical](https://mapzen.com/blog/osmlr-2nd-technical-preview/) previews of the OSMLR linear referencing system. Now we’re ready to share more about how we’re using [Mapzen Map Matching](https://mapzen.com/blog/map-matching/) to “snap” GPS-derived locations to OSMLR segments, and how we’re using a data-driven approach to evaluate and improve the algorithms._
# # A "data-driven" approach to improving map-matching - Part II:
# ## _PARAMETER TUNING_
# ============================================================================================
# In the [last](https://mapzen.com/blog/map-matching-validation/) blog post on Mapzen's map-matching service, we showed you how Mapzen uses synthetic GPS data to validate the results of our map-matching algorithm. This time around, we'll dive a _bit_ deeper into the internals of the algorithm itself to see how we can use our validation metrics to fine-tune the map-matching parameters.
# ## 0. Setup test environment
# In[3]:
import os
import sys; sys.path.insert(0, os.path.abspath('..'));
import validator.validator as val
import numpy as np
import glob
import pandas as pd
import pickle
import seaborn as sns
from matplotlib import pyplot as plt
from IPython.display import Image
from IPython.core.display import HTML
get_ipython().run_line_magic('matplotlib', 'inline')
# #### User vars
# In[2]:
mapzenKey = os.environ.get('MAPZEN_API')
gmapsKey = os.environ.get('GOOGLE_MAPS')
cityName = 'San Francisco'
numRoutes = 200
# #### Load our old routes
# Check out the last [notebook](http://nbviewer.jupyter.org/github/opentraffic/reporter-quality-testing-rig/blob/33868a9703f9c00d42dd520339eac3ec04fb4ea0/notebooks/map_matching_part_I.ipynb) in this series if you don't have routes yet. Or just read along!
# In[3]:
routeList = pickle.load(open('{0}_{1}_routes.pkl'.format(cityName, numRoutes), 'rb'))
# ## 1. Parameter Definitions
# The Open Traffic Reporter map-matching service is based on the Hidden Markov Model (HMM) design of [Newson and Krumm (2009)](http://research.microsoft.com/en-us/um/people/jckrumm/Publications%202009/map%20matching%20ACM%20GIS%20camera%20ready.pdf). HMM's define a class of models that utilize a directed graph structure to represent the probability of observing a particular sequence of events -- or in our case, a particular sequence of road segments that defines a route. For our purposes, it is enough to know that HMM's are, in general, governed by two probability distributions which we've parameterized using $\sigma_z$ and $\beta$, respectively. If, however, you are interested in a more detailed explanation of how HMM's work in the context of map-matching, please see the excellent map-matching primer [here](https://github.com/valhalla/meili/blob/master/docs/meili/algorithms.md) written by Mapzen's own routing experts.
# 
# *The graph structure of an HMM for an example route (from Mapzen's map-matching algorithms [primer](https://github.com/valhalla/meili/blob/master/docs/meili/algorithms.md))*
# ### $\sigma_z$
#
# The parameter $\sigma_z$ appears in the equation below which expresses the likelihood of recording a GPS measurement $z_t$ from road segment $r_i$ as the following Gaussian probability density function (PDF):
#
# ### $$ p(z_t|r_i) = \frac{1}{\sqrt{2 \pi \sigma_z}} e^{-0.5(\frac{||z_t - x_{t,i}||_{\text{great circle}}}{\sigma_z})^2}$$
# where $x_{t,i}$ is the point on road segment $r_i$ nearest the measurement $z_t$ at time $t$.
# In practice, $\sigma_z$ can be thought of as an estimate of the **standard deviation of GPS noise**. The more we trust our measurements, the lower the value of our $\sigma_z$ should be. Newson and Krumm (2009) derive $\sigma_z$ from the median absolute deviation over their dataset, arriving at a value of 4.07, which we have adopted as our default parameter value.
# ### $\beta$
#
# The $\beta$ parameter comes from the following exponential PDF:
# ### $$p(d_t) = \frac{1}{\beta}e^{\frac{-d_t}{\beta}}$$
# where $d_t$ is the difference between the great circle distance and route-traveled distance between time $t$ and $t+1$. We use an exponential PDF to define this probability because Newson and Krumm (2009) found that for two successive GPS measurements, the differences in length between the "great circle" distance (i.e. "as-the-crow-flies" distances) and the corresponding route-traveled distance are exponentially distributed.
# 
# From Newson and Krumm (2009)
# In this context, $\beta$ can be thought of as the expected difference between great circle distances and route distance traveled, or in other words, the **degree of circuitousness** of GPS routes. We have adopted a $\beta$ of 3 as our default parameter value.
# ## 2. Grid search
# Although we can empirically derive our parameter values using the equations above, we can also just search the space of feasible parameter values to find the ones that give the best results. This is a common machine learning approach for algorithm optimization, typically referred to as a **grid search** or parameter sweep. Grid search requires a reliable scoring mechanism in order to compare results. Lucky for us, we already implemented a number of these in our exploration of [validation](http://nbviewer.jupyter.org/github/opentraffic/reporter-quality-testing-rig/blob/33868a9703f9c00d42dd520339eac3ec04fb4ea0/notebooks/map_matching_part_I.ipynb) metrics.
# Since we have to pick one to use in the grid search, we'll stick with Newson and Krumm's distance-based error metric for now, since it appears most frequently in the literature. Let's take a quick look back at what those errors looked like over our set of validation data:
# 
# Clearly the quality of our data, by way of GPS noise and sample rate, has a huge impact on match error. Taking that into consideration, in addition to searching over the feasible space of $\beta$'s and $\sigma_z$'s, we'll want to iterate over different noise levels and sample rates so that we can optimize our map-matching regardless of the quality of the data we receive.
# ### a) Define the search space
# Because this grid search takes place in such high dimensional space - $\text{noise levels} \times \text{sample rates} \times \beta \ \text{values} \times \sigma\ \text{values} \times \text{# routes}$ - we use a smaller number of discrete noise levels for this section of the analysis:
# In[4]:
noiseLevels = np.linspace(0, 100, 6) # in meters
sampleRates = [1, 5, 10, 20, 30] # in seconds
betas = np.linspace(0.5, 10, 20) # no units
sigmas = np.linspace(0.5, 10, 20) # no units
# ### b) Search and score
# We've wrapped the entire grid search code into the `grid_search_hmm_params()` function since the internals are pretty messy. Behind the scenes, though, `grid_search_hmm_params()` proceeds as follows:
# ```python
# for route in routes:
# for beta in betas:
# for sigma in sigmas:
# for noise in noiseLevels:
# for rate in sampleRates:
# i) simulate gps
# ii) match routes
# iii) score the route
# write scores to a route-specific file
# ```
# Depending on the dimensions of your search space, grid search can be a long and arduous process. In our case, on a 2015 MacBook Pro, it seems to take ~1 min. per route for a given noise level and sample rate. With the configuration specified above, that works out to roughly 30 minutes per route. While this runtime does indeed seem a bit absurd, the beauty of grid search is that if you properly define your search space, you should really only need to do it once.
# In[5]:
val.grid_search_hmm_params(cityName, routeList, sampleRates, noiseLevels, betas, sigmas)
# ## 3. Get optimal params
# ### a) Combine the scores for all routes
# Because of the slow runtime, `grid_search_hmm_params()` saves its results incrementally in route-specific files. We must first recombine them into a single master dataframe so that we can query, plot, and otherwise manipulate the data all together
# In[6]:
frame = val.mergeCsvsToDataFrame(dataDir='../data/sf200/', matchString='*.csv')
# The first five rows of the combined dataframe:
# In[8]:
frame.head()
# ### b) Visualize scores
# Now we visually can assess the parameter scores to check for global optima or tease out any other patterns that might persist across the various dimensions of our grid search. The `get_param_scores()` function will also return a dictionary of optimal parameter values corresponding to each one of the subplots below:
# In[9]:
paramDict = val.get_param_scores(frame, sampleRates, noiseLevels, betas, sigmas, plot=True)
# ### c) Making sense of the results
# In our case there does appear to be a "global" trend of decreasing error toward the upper right quadrants of our plots, which correspond to higher $\sigma_z$'s and lower $\beta$'s. If we wanted to optimize our parameters across all sample rates and noise levels at once, the results suggest that our optimal parameters would lie at the maximum $\sigma_z$ and minimum $\beta$ values of our search space. This is an interesting finding, and not entirely satistfying, since it is typically unwise to pick parameter values at the limits of your search space -- who knows what may lie on the other side? Although this could be evidence that we didn't cast a wide enough net in our parameter sweep, there is reason to believe that it might not be worthwhile to search any further:
#
# 1. Focusing specifically at the plots in the 0 - 40 m range of noise, it's pretty clear that expanding the range of feasible parameter values wouldn't do any good. **You can't get much better than 0 error, and our plots in this range are already dominated by dark green**. We can obtain optimal results by picking parameters anywhere in this dark green zone.
#
#
#