Fix emeval.metrics.reference_trajectory for trajectory evaluation notebooks

[singhish]

@shankari as I had mentioned earlier today, trajectory_evaluation_timestamp_metric_selection and trajectory_evaluation_spatio_temporal are currently failing due to issues with the emeval.metrics.reference_trajectory module.

The error in the former notebook occurs at:

ref_tree = {}
coverage_list = []
coverage_list.extend(get_reference_coverage(pv_la, ref_tree))
coverage_list.extend(get_reference_coverage(pv_sj, ref_tree))
coverage_list.extend(get_reference_coverage(pv_ucb, ref_tree))

coverage_df = pd.DataFrame(coverage_list)

and the stack trace is as follows:

ValueError                                Traceback (most recent call last)
<ipython-input-33-2252b9dd9ae9> in <module>
      3 coverage_list.extend(get_reference_coverage(pv_la, ref_tree))
      4 coverage_list.extend(get_reference_coverage(pv_sj, ref_tree))
----> 5 coverage_list.extend(get_reference_coverage(pv_ucb, ref_tree))
      6 
      7 coverage_df = pd.DataFrame(coverage_list)

<ipython-input-31-95ac980592cc> in get_reference_coverage(pv, ref_tree_root)
     74     ref_tree = get_reference_trajectory_tree(pv)
     75     ref_tree_root[pv.spec_details.curr_spec["id"]] = ref_tree
---> 76     [fill_ref_tree_entry(e) for e in ref_tree.values()]
     77     return get_coverage(pv, ref_tree)

<ipython-input-31-95ac980592cc> in <listcomp>(.0)
     74     ref_tree = get_reference_trajectory_tree(pv)
     75     ref_tree_root[pv.spec_details.curr_spec["id"]] = ref_tree
---> 76     [fill_ref_tree_entry(e) for e in ref_tree.values()]
     77     return get_coverage(pv, ref_tree)

<ipython-input-31-95ac980592cc> in fill_ref_tree_entry(e)
     49     # So the "latitude" and "longitude" here is actually in meters and the points can't be displayed directly
     50     # TODO: figure out the best way to support both
---> 51     e["travel_forward_25"] = emr.ref_travel_forward(e, 25, tz=curr_tz)
     52     e["travel_forward_50"] = emr.ref_travel_forward(e, 50, tz=curr_tz)
     53     e["ct_midpoint_25"] = emr.ref_ct_general(e, emr.b_merge_midpoint, 25, tz=curr_tz)

~/mobilitynet-analysis-scripts/emeval/metrics/reference_trajectory.py in ref_travel_forward(e, dist_threshold, tz)
    362         emd.to_geo_df(e["temporal_control"]["ios"]["location_df"]),
    363         section_gt_shapes.filter(["start_loc","end_loc"]))
--> 364     new_location_df_a = get_int_aligned_trajectory(filtered_utm_loc_df_a, tz)
    365     new_location_df_i = get_int_aligned_trajectory(filtered_utm_loc_df_b, tz)
    366 

~/mobilitynet-analysis-scripts/emeval/metrics/reference_trajectory.py in get_int_aligned_trajectory(location_df, tz)
     41 
     42 def get_int_aligned_trajectory(location_df, tz="UTC"):
---> 43     lat_fn = sci.interp1d(x=location_df.ts, y=location_df.latitude)
     44     lon_fn = sci.interp1d(x=location_df.ts, y=location_df.longitude)
     45     # In order to avoid extrapolation, we use ceil for the first int and floor

~/miniconda3/envs/emissioneval/lib/python3.8/site-packages/scipy/interpolate/interpolate.py in __init__(self, x, y, kind, axis, copy, bounds_error, fill_value, assume_sorted)
    466         # Interpolation goes internally along the first axis
    467         self.y = y
--> 468         self._y = self._reshape_yi(self.y)
    469         self.x = x
    470         del y, x  # clean up namespace to prevent misuse; use attributes

~/miniconda3/envs/emissioneval/lib/python3.8/site-packages/scipy/interpolate/polyint.py in _reshape_yi(self, yi, check)
    109                                            self._y_extra_shape[:-self._y_axis])
    110             raise ValueError("Data must be of shape %s" % ok_shape)
--> 111         return yi.reshape((yi.shape[0], -1))
    112 
    113     def _set_yi(self, yi, xi=None, axis=None):

ValueError: cannot reshape array of size 0 into shape (0,newaxis)

Here, it is clear that the error occurs in reference_trajectory's get_int_aligned_trajectory function (line 42).

For the latter notebook, the error occurs at:

ref_tree = {}
get_reference_trajectory_tree(pv_la, ref_tree)
get_reference_trajectory_tree(pv_sj, ref_tree)
get_reference_trajectory_tree(pv_ucb, ref_tree)

with the following stack trace:

AttributeError                            Traceback (most recent call last)
<ipython-input-12-10e416badc4c> in <module>
      2 get_reference_trajectory_tree(pv_la, ref_tree)
      3 get_reference_trajectory_tree(pv_sj, ref_tree)
----> 4 get_reference_trajectory_tree(pv_ucb, ref_tree)

<ipython-input-11-00a6964e029f> in get_reference_trajectory_tree(pv, ref_tree_root)
     55     curr_ref_tree = get_reference_trajectory_input_tree(pv)
     56     ref_tree_root[pv.spec_details.CURR_SPEC_ID] = curr_ref_tree
---> 57     [fill_ref_tree_entry(pv, e) for e in curr_ref_tree.values()]

<ipython-input-11-00a6964e029f> in <listcomp>(.0)
     55     curr_ref_tree = get_reference_trajectory_input_tree(pv)
     56     ref_tree_root[pv.spec_details.CURR_SPEC_ID] = curr_ref_tree
---> 57     [fill_ref_tree_entry(pv, e) for e in curr_ref_tree.values()]

<ipython-input-11-00a6964e029f> in fill_ref_tree_entry(pv, e)
     50     curr_tz = pv.spec_details.eval_tz
     51     assert "android" in e["temporal_control"] and "ios" in e["temporal_control"]
---> 52     (e["reference_algo"], e["reference_df"]) = emr.final_ref_ensemble(e, 25, tz=curr_tz)
     53 
     54 def get_reference_trajectory_tree(pv, ref_tree_root):

~/mobilitynet-analysis-scripts/emeval/metrics/reference_trajectory.py in final_ref_ensemble(e, dist_threshold, tz)
    432     ct_stats = {
    433         "coverage_density": coverage_density(ct_ref_df, e),
--> 434         "coverage_time": coverage_time(ct_ref_df, e),
    435         "coverage_max_gap": coverage_max_gap(ct_ref_df, e)
    436     }

~/mobilitynet-analysis-scripts/emeval/metrics/reference_trajectory.py in <lambda>(df, sr)
    417 
    418 coverage_density = lambda df, sr: len(df)/(sr["end_ts"] - sr["start_ts"])
--> 419 coverage_time = lambda df, sr: (df.ts.iloc[-1] - df.ts.iloc[0])/(sr["end_ts"] - sr["start_ts"])
    420 coverage_max_gap = lambda df, sr: df.ts.diff().max()/(sr["end_ts"] - sr["start_ts"])
    421 

~/miniconda3/envs/emissioneval/lib/python3.8/site-packages/pandas/core/generic.py in __getattr__(self, name)
   5272             if self._info_axis._can_hold_identifiers_and_holds_name(name):
   5273                 return self[name]
-> 5274             return object.__getattribute__(self, name)
   5275 
   5276     def __setattr__(self, name: str, value) -> None:

AttributeError: 'GeoDataFrame' object has no attribute 'ts'

Here, final_ref_ensemble is the problematic function.

These errors need to be fixed for the above notebooks to work properly.

[shankari]

@singhish the first error just indicates that there is no data for the interpolation

ValueError: cannot reshape array of size 0 into shape (0,newaxis)

You need to verify that you are in fact storing and loading the data correctly from the files.

Note from #19 that both the trajectory notebooks were working properly when we were reading data from the server.

[shankari]

wait a minute, I don't see any confirmation in #19 that you got the metric selection notebooks to work.
I thought we got all the notebooks to work before starting on the refactor for the file specs.
Can you confirm?

[shankari]

I triple checked, and trajectory_evaluation_spatio_temporal is listed as working in #19 (#19 (comment)). Please try running it on master to see whether there is an error with the reference trajectory code or a regression in your changes.

[shankari]

And for trajectory_evaluation_timestamp_metric_selection, as I said in #25 (comment) the error is that there are no location points.

The error report does not indicate the run which caused the error, but on reproducing, I found that this was

Considering entry mtv_to_berkeley_sf_bart walk_to_bus 11

Given the shortness of the walk_to_bus leg, it is very likely that there were no trajectory points outside the start and end polygons. Would be good to display that and validate. The reason this wasn't encountered earlier is because we developed the notebooks in Fall 2019 before the last round of data collection that generated run 11

[shankari]

Yes, looking at the 11 evaluation range and the walk_to_bus section using

>>> third_repetition = pv.map()["ios"]["ucb-sdb-ios-1"]["evaluation_ranges"][11]; third_repetition["trip_id"]

>>> to_ucb_trip = third_repetition["evaluation_trip_ranges"][0]; print(to_ucb_trip["trip_id"])
>>> walk_leg = to_ucb_trip["evaluation_section_ranges"][6]; walk_leg["trip_id"] 
>>> gt_leg = sd.get_ground_truth_for_leg(to_ucb_trip["trip_id_base"], walk_leg["trip_id_base"]); gt_leg["id"]

We see that the entire trajectory is within the start (green) polygon

[shankari]

Ran into this again while generating reference trajectories to allow us to generate a distance based confusion matrix.
Gah! Should have checked issues before investigating.

There were actually two separate issues.

1. In get_ground_truth_for_leg, we pick the correct entry from the ground truth list, but then overwrite the list in the spec with the single matching value. This ensures that downstream analyses continue to work, but it means that after the first call to get_ground_truth_for_leg, the spec only has the first set of values. So when we call get_ground_truth_for_leg at a time when we should get the second entry in the list, it is no longer there!
   https://github.com/MobilityNet/mobilitynet-analysis-scripts/blame/eac11f8356dae7825257a40495183f846077d746/emeval/input/spec_details.py#L107

This is why the map in the comment above has the end_loc polygon way off from the actual trip trajectory. The real map should be as in the screenshot below (left: android, right: iOS). The fix is simple, we need to modify a copy() of the original spec instead of modifying the spec directly.

2. Even after fixing that, the entire sensed trajectory is within the green start_loc polygon. This is because during the last data collection on the UCB timeline, the bus stop for the city bus was changed. So the start_loc and end_loc polygon now overlap. This means that when we filter for points outside the start and end polygon, the length is zero. So when we try to match against spatial ground truth, the loc_df is empty and we get an exception

There are a couple of ways to fix this:

1. We can change the start_loc ground truth polygon as well so that the start and end loc don't overlap. There is only one trip with this new ground truth (taken after the thesis was complete) so we can narrow the start_loc to fit fairly tightly around the start location of the trip. This is likely to give better results because we have more of the trip to match to ground truth.
2. We can make the ref_travel_forward to be more robust to zero entries outside the start and end loc. This is likely better as we collect more timeline data

Or we could do both!

[shankari]

If we are ever going to do (1) we should do it now, so that the results are consistent and reproducible later. This has not been published yet - it was too late to make my thesis (#25 (comment)), and we haven't updated the spatio-temporal results later. So we should fix before the publish.

We should think about what to do for (2). If there are no points outside the start and end loc, we don't have a ground truth trajectory that we can project to. Presumably we will then use a non_gt approach for the trajectory inside the start and end loc? We should fix that when we clean up and understand the reference trajectory implementations.

Fixing (1), including in the mobilitynet mongodump, now.

[shankari]

For the record, before I re-run my notebook, I originally thought that this might be because some part of the trajectory was outside the end_loc, but that happens in run 1 as well without any errors

map_list = []
for phone_os in ["android", "ios"]:
    breaking_e = ref_tree['train_bus_ebike_mtv_ucb']['mtv_to_berkeley_sf_bart/walk_to_bus_1']
    curr_map = folium.Map()
    gt_leg_gj = pv_ucb.spec_details.get_geojson_for_leg(breaking_e["ground_truth"]["leg"])
    sensed_section_gj = ezgj.get_geojson_for_loc_df(breaking_e["temporal_control"][phone_os]["location_df"])
    gt_leg_gj_feature = folium.GeoJson(gt_leg_gj, name="ground_truth")
    sensed_leg_gj_feature = folium.GeoJson(sensed_section_gj, name="sensed_values")
    section_gt_shapes = breaking_e["ground_truth"]["gt_shapes"]
    filtered_loc = emd.filter_geo_df(breaking_e["temporal_control"][phone_os]["location_df"],
                      section_gt_shapes.filter(["start_loc","end_loc"]))
    filtered_section_gj = ezgj.get_geojson_for_loc_df(filtered_loc)
    filtered_section_gj_feature = folium.GeoJson(filtered_section_gj, name="filtered_values")
    curr_map.add_child(gt_leg_gj_feature)
    curr_map.add_child(sensed_leg_gj_feature)
    curr_map.add_child(filtered_section_gj_feature)
    curr_map.fit_bounds(sensed_leg_gj_feature.get_bounds())
    folium.LayerControl().add_to(curr_map)
    map_list.append(curr_map)

import branca.element as bre
fig = bre.Figure()
for i, curr_map in enumerate(map_list):
    fig.add_subplot(1, 2, i+1).add_child(curr_map)
fig