trajectory_evaluation_analysis_master.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for reading and validating data\n",
    "import emeval.input.spec_details as eisd\n",
    "import emeval.input.phone_view as eipv\n",
    "import emeval.input.eval_view as eiev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualization helpers\n",
    "import emeval.viz.phone_view as ezpv\n",
    "import emeval.viz.eval_view as ezev\n",
    "import emeval.viz.geojson as ezgj\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Metrics helpers\n",
    "import emeval.metrics.dist_calculations as emd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For computation\n",
    "import numpy as np\n",
    "import math\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import shapely as shp\n",
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATASTORE_URL = \"http://cardshark.cs.berkeley.edu\"\n",
    "AUTHOR_EMAIL = \"shankari@eecs.berkeley.edu\"\n",
    "sd_la = eisd.SpecDetails(DATASTORE_URL, AUTHOR_EMAIL, \"unimodal_trip_car_bike_mtv_la\")\n",
    "sd_sj = eisd.SpecDetails(DATASTORE_URL, AUTHOR_EMAIL, \"car_scooter_brex_san_jose\")\n",
    "sd_ucb = eisd.SpecDetails(DATASTORE_URL, AUTHOR_EMAIL, \"train_bus_ebike_mtv_ucb\")\n",
    "sd_ucb_reroute = eisd.SpecDetails(DATASTORE_URL, AUTHOR_EMAIL, \"train_bus_ebike_sm_reroute_mtv_ucb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "importlib.reload(eisd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv_la = eipv.PhoneView(sd_la)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv_sj = eipv.PhoneView(sd_sj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv_ucb = eipv.PhoneView(sd_ucb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv_ucb_reroute = eipv.PhoneView(sd_ucb_reroute)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import emeval.analysed.phone_view as eapv\n",
    "importlib.reload(eapv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "av_la = eapv.create_analysed_view(pv_la, \"http://localhost:8080\", \"analysis/recreated_location\", \"analysis/cleaned_trip\", \"analysis/cleaned_section\")\n",
    "av_sj = eapv.create_analysed_view(pv_sj, \"http://localhost:8080\", \"analysis/recreated_location\", \"analysis/cleaned_trip\", \"analysis/cleaned_section\")\n",
    "av_ucb = eapv.create_analysed_view(pv_ucb, \"http://localhost:8080\", \"analysis/recreated_location\", \"analysis/cleaned_trip\", \"analysis/cleaned_section\")\n",
    "av_ucb_reroute = eapv.create_analysed_view(pv_ucb_reroute, \"http://localhost:8080\", \"analysis/recreated_location\", \"analysis/cleaned_trip\", \"analysis/cleaned_section\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[se[\"metadata\"][\"write_fmt_time\"] for se in av_sj.map()[\"ios\"][\"ucb-sdb-ios-2\"][\"evaluation_ranges\"][0][\"sensed_trip_ranges\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Validate distance calculations\n",
    "\n",
    "Our x,y coordinates are in degrees (lon, lat). So when we calculate the distance between two points, it is also in degrees. In order for this to be meaningful, we need to convert it to a regular distance metric such as meters.\n",
    "\n",
    "This is a complicated problem in general because our distance calculation applies 2-D spatial operations to a 3-D curved space. However, as documented in the shapely documentation, since our areas of interest are small, we can use a 2-D approximation and get reasonable results.\n",
    "\n",
    "In order to get distances from degree-based calculations, we can use the following options:\n",
    "- perform the calculations in degrees and then convert them to meters. As an approximation, we can use the fact that 360 degrees represents the circumference of the earth. Therefore `dist = degree_dist * (C/360)`\n",
    "- convert degrees to x,y coordinates using utm (https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system) and then calculate the distance\n",
    "- since we calculate the distance from the ground truth linestring, calculate the closest ground truth point in (lon,lat) and then use the haversine formula (https://en.wikipedia.org/wiki/Haversine_formula) to calculate the distance between the two points\n",
    "\n",
    "Let us quickly all three calculations for three selected test cases and:\n",
    "- check whether they are largely consistent\n",
    "- compare with other distance calculators to see which are closer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results and method choice\n",
    "\n",
    "We find that the `manual_utm` and `project` methods are pretty consistent, and are significantly different from the `circumference` method. The `circumference` method appears to be consistently greater than the other two and the difference appears to be around 25%. The manual checks also appear to be closer to the `manual_utm` and `project` values. The `manual_utm` and `project` values are consistently within ~ 5% of each other, so we could really use either one.\n",
    "\n",
    "**We will use the utm approach** since it is correct, is consistent with the shapely documentation (https://shapely.readthedocs.io/en/stable/manual.html#coordinate-systems) and applicable to operations beyond distance calculation\n",
    "\n",
    "> Even though the Earth is not flat – and for that matter not exactly spherical – there are many analytic problems that can be approached by transforming Earth features to a Cartesian plane, applying tried and true algorithms, and then transforming the results back to geographic coordinates. This practice is as old as the tradition of accurate paper maps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial error calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_spatial_errors(pv):\n",
    "    spatial_error_df = pd.DataFrame()\n",
    "    \n",
    "    for phone_os, phone_map in pv.map().items():\n",
    "        for phone_label, phone_detail_map in phone_map.items():\n",
    "            for (r_idx, r) in enumerate(phone_detail_map[\"evaluation_ranges\"]):\n",
    "                run_errors = []\n",
    "                for (tr_idx, tr) in enumerate(r[\"evaluation_trip_ranges\"]):\n",
    "                    trip_errors = []\n",
    "                    for (sr_idx, sr) in enumerate(tr[\"evaluation_section_ranges\"]):\n",
    "                        # This is a Shapely LineString\n",
    "                        \n",
    "                        section_gt_leg = pv.spec_details.get_ground_truth_for_leg(tr[\"trip_id_base\"], sr[\"trip_id_base\"])\n",
    "                        section_gt_shapes = gpd.GeoSeries(eisd.SpecDetails.get_shapes_for_leg(section_gt_leg))\n",
    "                        if len(section_gt_shapes) == 1:\n",
    "                            print(\"No ground truth route for %s %s, must be polygon, skipping...\" % (tr[\"trip_id_base\"], sr[\"trip_id_base\"]))\n",
    "                            assert section_gt_leg[\"type\"] != \"TRAVEL\", \"For %s, %s, %s, %s, %s found type %s\" % (phone_os, phone_label, r_idx, tr_idx, sr_idx, section_gt_leg[\"type\"])\n",
    "                            continue\n",
    "                        if len(sr['location_df']) == 0:\n",
    "                            print(\"No sensed locations found, role = %s skipping...\" % (r[\"eval_role_base\"]))\n",
    "                            # assert r[\"eval_role_base\"] == \"power_control\", \"Found no locations for %s, %s, %s, %s, %s\" % (phone_os, phone_label, r_idx, tr_idx, sr_idx)\n",
    "                            continue\n",
    "                            \n",
    "                        print(\"Processing travel leg %s, %s, %s, %s, %s\" %\n",
    "                              (phone_os, phone_label, r[\"eval_role_base\"], tr[\"trip_id_base\"], sr[\"trip_id_base\"]))\n",
    "                        # This is a GeoDataFrame\n",
    "                        section_geo_df = emd.to_geo_df(sr[\"location_df\"])\n",
    "                        \n",
    "                        # After this point, everything is in UTM so that 2-D inside/filtering operations work\n",
    "                        utm_section_geo_df = emd.to_utm_df(section_geo_df)\n",
    "                        utm_section_gt_shapes = section_gt_shapes.apply(lambda s: shp.ops.transform(emd.to_utm_coords, s))\n",
    "                        filtered_us_gpdf = emd.filter_geo_df(utm_section_geo_df, utm_section_gt_shapes.loc[\"start_loc\":\"end_loc\"])\n",
    "                        filtered_gt_linestring = emd.filter_ground_truth_linestring(utm_section_gt_shapes)\n",
    "                        meter_dist = filtered_us_gpdf.geometry.distance(filtered_gt_linestring)\n",
    "                        ne = len(meter_dist)\n",
    "                        curr_spatial_error_df = gpd.GeoDataFrame({\"error\": meter_dist,\n",
    "                                                                  \"ts\": section_geo_df.ts,\n",
    "                                                                  \"geometry\": section_geo_df.geometry,\n",
    "                                                                  \"phone_os\": np.repeat(phone_os, ne),\n",
    "                                                                  \"phone_label\": np.repeat(phone_label, ne),\n",
    "                                                                  \"role\": np.repeat(r[\"eval_role_base\"], ne),\n",
    "                                                                  \"timeline\": np.repeat(pv.spec_details.CURR_SPEC_ID, ne), \n",
    "                                                                  \"run\": np.repeat(r_idx, ne),\n",
    "                                                                  \"trip_id\": np.repeat(tr[\"trip_id_base\"], ne),\n",
    "                                                                  \"section_id\": np.repeat(sr[\"trip_id_base\"], ne)})\n",
    "                        spatial_error_df = pd.concat([spatial_error_df, curr_spatial_error_df], axis=\"index\")\n",
    "    return spatial_error_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df = pd.DataFrame()\n",
    "spatial_errors_df = pd.concat([spatial_errors_df, get_spatial_errors(av_la)], axis=\"index\")\n",
    "spatial_errors_df = pd.concat([spatial_errors_df, get_spatial_errors(av_sj)], axis=\"index\")\n",
    "spatial_errors_df = pd.concat([spatial_errors_df, get_spatial_errors(av_ucb)], axis=\"index\")\n",
    "spatial_errors_df = pd.concat([spatial_errors_df, get_spatial_errors(av_ucb_reroute)], axis=\"index\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's merge the reroutes\n",
    "spatial_errors_df.timeline.replace(\"train_bus_ebike_sm_reroute_mtv_ucb\", \"train_bus_ebike_mtv_ucb\", inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r2q_map = {\"power_control\": 0, \"HAMFDC\": 1, \"MAHFDC\": 2, \"HAHFDC\": 3, \"accuracy_control\": 4}\n",
    "q2r_map = {0: \"power\", 1: \"HAMFDC\", 2: \"MAHFDC\", 3: \"HAHFDC\", 4: \"accuracy\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df[\"quality\"] = spatial_errors_df.role.apply(lambda r: r2q_map[r])\n",
    "spatial_errors_df[\"label\"] = spatial_errors_df.role.apply(lambda r: r.replace('_control', ''))\n",
    "timeline_list = [\"train_bus_ebike_mtv_ucb\", \"car_scooter_brex_san_jose\", \"unimodal_trip_car_bike_mtv_la\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overall stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifig, ax_array = plt.subplots(nrows=1,ncols=2,figsize=(8,2), sharey=True)\n",
    "\n",
    "spatial_errors_df.query(\"phone_os == 'android' & quality > 0\").boxplot(ax = ax_array[0], column=[\"error\"], by=[\"quality\"], showfliers=False)\n",
    "ax_array[0].set_title('android')\n",
    "spatial_errors_df.query(\"phone_os == 'ios' & quality > 0\").boxplot(ax = ax_array[1], column=[\"error\"], by=[\"quality\"], showfliers=False)\n",
    "ax_array[1].set_title(\"ios\")\n",
    "\n",
    "for i, ax in enumerate(ax_array):\n",
    "    # print([t.get_text() for t in ax.get_xticklabels()])\n",
    "    ax.set_xticklabels([q2r_map[int(t.get_text())] for t in ax.get_xticklabels()])\n",
    "    ax.set_xlabel(\"\")\n",
    "\n",
    "ax_array[0].set_ylabel(\"Spatial error (meters)\")\n",
    "# ax_array[1][0].set_ylabel(\"Spatial error (meters)\")\n",
    "ifig.suptitle(\"Spatial trajectory error v/s quality (excluding outliers)\", y = 1.1)\n",
    "# ifig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifig, ax_array = plt.subplots(nrows=1,ncols=2,figsize=(8,2), sharey=True)\n",
    "\n",
    "spatial_errors_df.query(\"phone_os == 'android' & quality > 0\").boxplot(ax = ax_array[0], column=[\"error\"], by=[\"quality\"])\n",
    "ax_array[0].set_title('android')\n",
    "spatial_errors_df.query(\"phone_os == 'ios' & quality > 0\").boxplot(ax = ax_array[1], column=[\"error\"], by=[\"quality\"])\n",
    "ax_array[1].set_title(\"ios\")\n",
    "\n",
    "for i, ax in enumerate(ax_array):\n",
    "    # print([t.get_text() for t in ax.get_xticklabels()])\n",
    "    ax.set_xticklabels([q2r_map[int(t.get_text())] for t in ax.get_xticklabels()])\n",
    "    ax.set_xlabel(\"\")\n",
    "\n",
    "ax_array[0].set_ylabel(\"Spatial error (meters)\")\n",
    "# ax_array[1][0].set_ylabel(\"Spatial error (meters)\")\n",
    "ifig.suptitle(\"Spatial trajectory error v/s quality\", y = 1.1)\n",
    "# ifig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Split out results by timeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifig, ax_array = plt.subplots(nrows=2,ncols=3,figsize=(12,6), sharex=False, sharey=False)\n",
    "timeline_list = [\"train_bus_ebike_mtv_ucb\", \"car_scooter_brex_san_jose\", \"unimodal_trip_car_bike_mtv_la\"]\n",
    "for i, tl in enumerate(timeline_list):\n",
    "    spatial_errors_df.query(\"timeline == @tl & phone_os == 'android' & quality > 0\").boxplot(ax = ax_array[0][i], column=[\"error\"], by=[\"quality\"])\n",
    "    ax_array[0][i].set_title(tl)\n",
    "    spatial_errors_df.query(\"timeline == @tl & phone_os == 'ios' & quality > 0\").boxplot(ax = ax_array[1][i], column=[\"error\"], by=[\"quality\"])\n",
    "    ax_array[1][i].set_title(\"\")\n",
    "\n",
    "for i, ax in enumerate(ax_array[0]):\n",
    "    ax.set_xticklabels([q2r_map[int(t.get_text())] for t in ax.get_xticklabels()])\n",
    "    ax.set_xlabel(\"\")\n",
    "\n",
    "for i, ax in enumerate(ax_array[1]):\n",
    "    ax.set_xticklabels([q2r_map[int(t.get_text())] for t in ax.get_xticklabels()])\n",
    "    ax.set_xlabel(\"\")\n",
    "\n",
    "ax_array[0][0].set_ylabel(\"Spatial error (android)\")\n",
    "ax_array[1][0].set_ylabel(\"Spatial error (iOS)\")\n",
    "ifig.suptitle(\"Spatial trajectory error v/s quality over multiple timelines\")\n",
    "# ifig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Split out results by section for the most complex timeline (train_bus_ebike_mtv_ucb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifig, ax_array = plt.subplots(nrows=2,ncols=4,figsize=(25,10), sharex=True, sharey=True)\n",
    "timeline_list = [\"train_bus_ebike_mtv_ucb\"]\n",
    "for i, tl in enumerate(timeline_list):\n",
    "    for q in range(1,5):\n",
    "        sel_df = spatial_errors_df.query(\"timeline == @tl & phone_os == 'android' & quality == @q\")\n",
    "        if len(sel_df) > 0:\n",
    "            sel_df.boxplot(ax = ax_array[2*i][q-1], column=[\"error\"], by=[\"section_id\"])\n",
    "        ax_array[2*i][q-1].tick_params(axis=\"x\", labelrotation=45)\n",
    "        sel_df = spatial_errors_df.query(\"timeline == @tl & phone_os == 'ios' & quality == @q\")\n",
    "        if len(sel_df) > 0:\n",
    "            sel_df.boxplot(ax = ax_array[2*i+1][q-1], column=[\"error\"], by=[\"section_id\"])\n",
    "#        ax_array[i][].set_title(\"\")\n",
    "\n",
    "def make_acronym(s):\n",
    "    ssl = s.split(\"_\")\n",
    "    # print(\"After splitting %s, we get %s\" % (s, ssl))\n",
    "    if len(ssl) == 0 or len(ssl[0]) == 0:\n",
    "        return \"\"\n",
    "    else:\n",
    "        return \"\".join([ss[0] for ss in ssl])\n",
    "\n",
    "for q in range(1,5):\n",
    "    ax_array[0][q-1].set_title(q2r_map[q])\n",
    "    curr_ticks = [t.get_text() for t in ax_array[1][q-1].get_xticklabels()]\n",
    "    new_ticks = [make_acronym(t) for t in curr_ticks]\n",
    "    ax_array[1][q-1].set_xticklabels(new_ticks)\n",
    "    \n",
    "print(list(zip(curr_ticks, new_ticks)))\n",
    "# fig.text(0,0,\"%s\"% list(zip(curr_ticks, new_ticks)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "timeline_list = [\"train_bus_ebike_mtv_ucb\"]\n",
    "for i, tl in enumerate(timeline_list):\n",
    "    unique_sections = spatial_errors_df.query(\"timeline == @tl\").section_id.unique()\n",
    "    ifig, ax_array = plt.subplots(nrows=2,ncols=len(unique_sections),figsize=(40,10), sharex=True, sharey=False)\n",
    "    for sid, s_name in enumerate(unique_sections):\n",
    "        sel_df = spatial_errors_df.query(\"timeline == @tl & phone_os == 'android' & section_id == @s_name & quality > 0\")\n",
    "        if len(sel_df) > 0:\n",
    "            sel_df.boxplot(ax = ax_array[2*i][sid], column=[\"error\"], by=[\"quality\"])\n",
    "        ax_array[2*i][sid].set_title(s_name)\n",
    "        sel_df = spatial_errors_df.query(\"timeline == @tl & phone_os == 'ios' & section_id == @s_name & quality > 0\")\n",
    "        if len(sel_df) > 0:\n",
    "            sel_df.boxplot(ax = ax_array[2*i+1][sid], column=[\"error\"], by=[\"quality\"])\n",
    "        ax_array[2*i+1][sid].set_title(\"\")\n",
    "#        ax_array[i][].set_title(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Focus only on sections where the max error is > 1000 meters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeline_list = [\"train_bus_ebike_mtv_ucb\"]\n",
    "for i, tl in enumerate(timeline_list):\n",
    "    unique_sections = pd.Series(spatial_errors_df.query(\"timeline == @tl\").section_id.unique())\n",
    "    sections_with_outliers_mask = unique_sections.apply(lambda s_name: spatial_errors_df.query(\"timeline == 'train_bus_ebike_mtv_ucb' & section_id == @s_name\").error.max() > 1000)\n",
    "    sections_with_outliers = unique_sections[sections_with_outliers_mask]   \n",
    "    ifig, ax_array = plt.subplots(nrows=2,ncols=len(sections_with_outliers),figsize=(17,4), sharex=True, sharey=False)\n",
    "    for sid, s_name in enumerate(sections_with_outliers):\n",
    "        sel_df = spatial_errors_df.query(\"timeline == @tl & phone_os == 'android' & section_id == @s_name & quality > 0\")\n",
    "        if len(sel_df) > 0:\n",
    "            sel_df.boxplot(ax = ax_array[2*i][sid], column=[\"error\"], by=[\"quality\"])\n",
    "        ax_array[2*i][sid].set_title(s_name)\n",
    "        ax_array[2*i][sid].set_xlabel(\"\")\n",
    "        sel_df = spatial_errors_df.query(\"timeline == @tl & phone_os == 'ios' & section_id == @s_name & quality > 0\")\n",
    "        if len(sel_df) > 0:\n",
    "            sel_df.boxplot(ax = ax_array[2*i+1][sid], column=[\"error\"], by=[\"quality\"])\n",
    "        ax_array[2*i+1][sid].set_title(\"\")\n",
    "        print([t.get_text() for t in ax_array[2*i+1][sid].get_xticklabels()])\n",
    "        ax_array[2*i+1][sid].set_xticklabels([q2r_map[int(t.get_text())] for t in ax_array[2*i+1][sid].get_xticklabels() if len(t.get_text()) > 0])\n",
    "        ax_array[2*i+1][sid].set_xlabel(\"\")\n",
    "    ifig.suptitle(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Validation of outliers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (express bus iOS, MAHFDC)\n",
    "\n",
    "ok, so it looks like the error is non-trivial across all runs, but run #1 is the worst and is responsible for the majority of the outliers. And this is borne out by the map, where on run #1, we end up with points in San Leandro!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'ios' & quality == 2 & section_id == 'express_bus' & error > 500\").run.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'ios' & quality == 2 & section_id == 'express_bus'\").boxplot(column=\"error\", by=\"run\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"berkeley_to_mtv_SF_express_bus\", \"express_bus\"); print(gt_leg[\"id\"])\n",
    "curr_map = ezgj.get_map_for_geojson(sd_ucb.get_geojson_for_leg(gt_leg), name=\"ground_truth\")\n",
    "ezgj.get_fg_for_loc_df(emd.linestring_to_geo_df(eisd.SpecDetails.get_shapes_for_leg(gt_leg)[\"route\"]),\n",
    "                       name=\"gt_points\", color=\"green\").add_to(curr_map)\n",
    "\n",
    "name_err_time = lambda lr: \"%d: %d, %s, %s\" % (lr[\"index\"], lr[\"df_idx\"], lr[\"error\"], sd_ucb.fmt(lr[\"ts\"], \"MM-DD HH:mm:ss\"))\n",
    "error_df = emd.to_loc_df(spatial_errors_df.query(\"phone_os == 'ios' & quality == 2 & section_id == 'express_bus' & run == 1\"))\n",
    "gt_16k = lambda lr: lr[\"error\"] == error_df.error.max()\n",
    "folium.GeoJson(ezgj.get_geojson_for_loc_df(error_df, color=\"red\"), name=\"sensed_values\").add_to(curr_map)\n",
    "ezgj.get_fg_for_loc_df(error_df, name=\"sensed_points\", color=\"red\", popupfn=name_err_time, stickyfn=gt_16k).add_to(curr_map)\n",
    "folium.LayerControl().add_to(curr_map)\n",
    "curr_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "importlib.reload(ezgj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"berkeley_to_mtv_SF_express_bus\", \"express_bus\"); print(gt_leg[\"id\"])\n",
    "curr_map = ezgj.get_map_for_geojson(sd_ucb.get_geojson_for_leg(gt_leg), name=\"ground_truth\")\n",
    "ezgj.get_fg_for_loc_df(emd.linestring_to_geo_df(eisd.SpecDetails.get_shapes_for_leg(gt_leg)[\"route\"]),\n",
    "                       name=\"gt_points\", color=\"green\").add_to(curr_map)\n",
    "\n",
    "name_err_time = lambda lr: \"%d: %d, %s, %s\" % (lr[\"index\"], lr[\"df_idx\"], lr[\"error\"], sd_ucb.fmt(lr[\"ts\"], \"MM-DD HH:mm:ss\"))\n",
    "\n",
    "colors = [\"red\", \"yellow\", \"blue\"]\n",
    "for run in range(3):\n",
    "    error_df = emd.to_loc_df(spatial_errors_df.query(\"phone_os == 'ios' & quality == 2 & section_id == 'express_bus' & run == @run\"))\n",
    "    gt_16k = lambda lr: lr[\"error\"] == error_df.error.max()\n",
    "    print(\"max error for run %d is %s\" % (run, error_df.error.max()))\n",
    "    folium.GeoJson(ezgj.get_geojson_for_loc_df(error_df, color=colors[run]), name=\"sensed_values\").add_to(curr_map)\n",
    "    ezgj.get_fg_for_loc_df(error_df, name=\"sensed_points\", color=colors[run], popupfn=name_err_time, stickyfn=gt_16k).add_to(curr_map)\n",
    "\n",
    "folium.LayerControl().add_to(curr_map)\n",
    "curr_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (commuter rail aboveground android, HAMFDC)\n",
    "\n",
    "Run 0: Multiple outliers at the start in San Jose. After that, everything is fine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & quality == 1 & section_id == 'commuter_rail_aboveground' & error > 500\").run.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & quality == 1 & section_id == 'commuter_rail_aboveground' & error > 500\").boxplot(column=\"error\", by=\"run\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"mtv_to_berkeley_sf_bart\", \"commuter_rail_aboveground\"); print(gt_leg[\"id\"])\n",
    "curr_map = ezgj.get_map_for_geojson(sd_ucb.get_geojson_for_leg(gt_leg), name=\"ground_truth\")\n",
    "ezgj.get_fg_for_loc_df(emd.linestring_to_geo_df(eisd.SpecDetails.get_shapes_for_leg(gt_leg)[\"route\"]),\n",
    "                       name=\"gt_points\", color=\"green\").add_to(curr_map)\n",
    "\n",
    "name_err_time = lambda lr: \"%d: %d, %s, %s\" % (lr[\"index\"], lr[\"df_idx\"], lr[\"error\"], sd_ucb.fmt(lr[\"ts\"], \"MM-DD HH:mm:ss\"))\n",
    "\n",
    "error_df = emd.to_loc_df(spatial_errors_df.query(\"phone_os == 'android' & quality == 1 & section_id == 'commuter_rail_aboveground' & run == 0\"))\n",
    "maxes = [error_df.error.max(), error_df[error_df.error < 10000].error.max(), error_df[error_df.error < 1000].error.max()]\n",
    "gt_16k = lambda lr: lr[\"error\"] in maxes\n",
    "folium.GeoJson(ezgj.get_geojson_for_loc_df(error_df, color=\"red\"), name=\"sensed_values\").add_to(curr_map)\n",
    "ezgj.get_fg_for_loc_df(error_df, name=\"sensed_points\", color=\"red\", popupfn=name_err_time, stickyfn=gt_16k).add_to(curr_map)\n",
    "\n",
    "folium.LayerControl().add_to(curr_map)\n",
    "curr_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & quality == 1 & section_id == 'commuter_rail_aboveground' & error > 10000\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (walk_to_bus android, HAMFDC, HAHFDC)\n",
    "\n",
    "Huge zig zag when we get out of the BART station"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & (quality == 1 | quality == 3) & section_id == 'walk_to_bus' & error > 500\").run.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & (quality == 1 | quality == 3) & section_id == 'walk_to_bus' & error > 500\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & (quality == 1 | quality == 3) & section_id == 'walk_to_bus'\").boxplot(column=\"error\", by=\"run\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & (quality == 1 | quality == 3) & section_id == 'walk_to_bus'\").error.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ucb_and_back = pv_ucb.map()[\"android\"][\"ucb-sdb-android-2\"][\"evaluation_ranges\"][0]; ucb_and_back[\"trip_id\"]\n",
    "to_trip = ucb_and_back[\"evaluation_trip_ranges\"][0]; print(to_trip[\"trip_id\"])\n",
    "wb_leg = to_trip[\"evaluation_section_ranges\"][6]; print(wb_leg[\"trip_id\"])\n",
    "gt_leg = sd_ucb.get_ground_truth_for_leg(to_trip[\"trip_id_base\"], wb_leg[\"trip_id_base\"]); gt_leg[\"id\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "importlib.reload(ezgj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"mtv_to_berkeley_sf_bart\", \"walk_to_bus\"); print(gt_leg[\"id\"])\n",
    "curr_map = ezgj.get_map_for_geojson(sd_ucb.get_geojson_for_leg(gt_leg), name=\"ground_truth\")\n",
    "ezgj.get_fg_for_loc_df(emd.linestring_to_geo_df(eisd.SpecDetails.get_shapes_for_leg(gt_leg)[\"route\"]),\n",
    "                       name=\"gt_points\", color=\"green\").add_to(curr_map)\n",
    "\n",
    "name_err_time = lambda lr: \"%d: %d, %s, %s\" % (lr[\"index\"], lr[\"df_idx\"], lr[\"error\"], sd_ucb.fmt(lr[\"ts\"], \"MM-DD HH:mm:ss\"))\n",
    "\n",
    "error_df = emd.to_loc_df(spatial_errors_df.query(\"phone_os == 'android' & quality == 1 & section_id == 'walk_to_bus'\").sort_index(axis=\"index\"))\n",
    "maxes = [error_df.error.max(), error_df[error_df.error < 16000].error.max(), error_df[error_df.error < 5000].error.max()]\n",
    "gt_16k = lambda lr: lr[\"error\"] in maxes\n",
    "print(\"Checking errors %s\" % maxes)\n",
    "folium.GeoJson(ezgj.get_geojson_for_loc_df(error_df, color=\"red\"), name=\"sensed_values\").add_to(curr_map)\n",
    "ezgj.get_fg_for_loc_df(error_df, name=\"sensed_points\", color=\"red\", popupfn=name_err_time, stickyfn=gt_16k).add_to(curr_map)\n",
    "\n",
    "folium.LayerControl().add_to(curr_map)\n",
    "curr_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (light_rail_below_above_ground, android, accuracy_control)\n",
    "\n",
    "ok, so it looks like the error is non-trivial across all runs, but run #1 is the worst and is responsible for the majority of the outliers. And this is borne out by the map, where on run #1, we end up with points in San Leandro!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & quality == 4 & section_id == 'light_rail_below_above_ground' & error > 100\").run.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# spatial_errors_df.query(\"phone_os == 'android' & (quality == 4) & section_id == 'light_rail_below_above_ground'\").boxplot(column=\"error\", by=\"run\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ucb_and_back = pv_ucb.map()[\"android\"][\"ucb-sdb-android-2\"][\"evaluation_ranges\"][0]; ucb_and_back[\"trip_id\"]\n",
    "back_trip = ucb_and_back[\"evaluation_trip_ranges\"][2]; print(back_trip[\"trip_id\"])\n",
    "lt_leg = back_trip[\"evaluation_section_ranges\"][7]; print(lt_leg[\"trip_id\"])\n",
    "gt_leg = sd_ucb.get_ground_truth_for_leg(back_trip[\"trip_id_base\"], lt_leg[\"trip_id_base\"]); gt_leg[\"id\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"berkeley_to_mtv_SF_express_bus\", \"light_rail_below_above_ground\"); print(gt_leg[\"id\"])\n",
    "curr_map = ezgj.get_map_for_geojson(sd_ucb.get_geojson_for_leg(gt_leg), name=\"ground_truth\")\n",
    "ezgj.get_fg_for_loc_df(emd.linestring_to_geo_df(eisd.SpecDetails.get_shapes_for_leg(gt_leg)[\"route\"]),\n",
    "                       name=\"gt_points\", color=\"green\").add_to(curr_map)\n",
    "\n",
    "name_err_time = lambda lr: \"%d: %d, %s, %s\" % (lr[\"index\"], lr[\"df_idx\"], lr[\"error\"], sd_ucb.fmt(lr[\"ts\"], \"MM-DD HH:mm:ss\"))\n",
    "\n",
    "colors = [\"red\", \"yellow\", \"blue\"]\n",
    "for run in range(3):\n",
    "    error_df = emd.to_loc_df(spatial_errors_df.query(\"phone_os == 'android' & quality == 2 & section_id == 'light_rail_below_above_ground' & run == @run\"))\n",
    "    gt_16k = lambda lr: lr[\"error\"] == error_df.error.max()\n",
    "    print(\"max error for run %d is %s\" % (run, error_df.error.max()))\n",
    "    folium.GeoJson(ezgj.get_geojson_for_loc_df(error_df, color=colors[run]), name=\"sensed_values\").add_to(curr_map)\n",
    "    ezgj.get_fg_for_loc_df(error_df, name=\"sensed_points\", color=colors[run], popupfn=name_err_time, stickyfn=gt_16k).add_to(curr_map)\n",
    "\n",
    "folium.LayerControl().add_to(curr_map)\n",
    "curr_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & quality == 2 & section_id == 'light_rail_below_above_ground' & run == @run\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (subway, android, HAMFDC)\n",
    "\n",
    "This is the poster child for temporal accuracy tracking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_errors_df.query(\"phone_os == 'android' & (quality == 1 | quality == 3) & section_id == 'subway_underground' & error > 8000\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bart_leg = pv_ucb.map()[\"android\"][\"ucb-sdb-android-3\"][\"evaluation_ranges\"][0][\"evaluation_trip_ranges\"][0][\"evaluation_section_ranges\"][5]\n",
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"mtv_to_berkeley_sf_bart\", \"subway_underground\"); gt_leg[\"id\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"mtv_to_berkeley_sf_bart\", \"subway_underground\"); print(gt_leg[\"id\"])\n",
    "curr_map = ezgj.get_map_for_geojson(sd_ucb.get_geojson_for_leg(gt_leg), name=\"ground_truth\")\n",
    "ezgj.get_fg_for_loc_df(emd.linestring_to_geo_df(eisd.SpecDetails.get_shapes_for_leg(gt_leg)[\"route\"]),\n",
    "                       name=\"gt_points\", color=\"green\").add_to(curr_map)\n",
    "\n",
    "name_err_time = lambda lr: \"%d: %d, %s, %s\" % (lr[\"index\"], lr[\"df_idx\"], lr[\"error\"], sd_ucb.fmt(lr[\"ts\"], \"MM-DD HH:mm:ss\"))\n",
    "\n",
    "error_df = emd.to_loc_df(spatial_errors_df.query(\"phone_os == 'android' & quality == 3 & section_id == 'subway_underground' & run == 2\").sort_index(axis=\"index\"))\n",
    "maxes = [error_df.error.max(), error_df[error_df.error < 16000].error.max(), error_df[error_df.error < 5000].error.max()]\n",
    "gt_16k = lambda lr: lr[\"error\"] in maxes\n",
    "print(\"Checking errors %s\" % maxes)\n",
    "folium.GeoJson(ezgj.get_geojson_for_loc_df(error_df, color=\"red\"), name=\"sensed_values\").add_to(curr_map)\n",
    "ezgj.get_fg_for_loc_df(error_df, name=\"sensed_points\", color=\"red\", popupfn=name_err_time, stickyfn=gt_16k).add_to(curr_map)\n",
    "\n",
    "folium.LayerControl().add_to(curr_map)\n",
    "curr_map\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_leg = sd_ucb.get_ground_truth_for_leg(\"mtv_to_berkeley_sf_bart\", \"subway_underground\"); gt_leg[\"id\"]\n",
    "eisd.SpecDetails.get_shapes_for_leg(gt_leg)[\"route\"].is_simple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pd.concat([\n",
    "    error_df.iloc[40:50],\n",
    "    error_df.iloc[55:60],\n",
    "    error_df.iloc[65:75],\n",
    "    error_df.iloc[70:75]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyproj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latlonProj = pyproj.Proj(init=\"epsg:4326\")\n",
    "xyProj = pyproj.Proj(init=\"epsg:3395\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy = pyproj.transform(latlonProj, xyProj, -122.08355963230133, 37.39091642895306); xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyproj.transform(xyProj, latlonProj, xy[0], xy[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame({\"a\": [1,2,3], \"b\": [4,5,6]}); df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.concat([pd.DataFrame([{\"a\": 10, \"b\": 14}]), df, pd.DataFrame([{\"a\": 20, \"b\": 24}])], axis='index').reset_index(drop=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}