Given an augment stop times DataFrame as output by the function
    :func:`build_augmented_stop_times`, a distance threshold (float)
    in the same units as the ``'shape_dist_traveled'`` column of
    ``augmented_stop_times``, if that column is present, and a delay
    threshold (integer number of seconds), alter the delay values
    of the augmented stop times as follows.

    Drop all delays with absolute value more than ``delay_threshold``
    seconds.
    For each trip and for each delay type (arrival delay or departure
    delay) do the following.
    If the trip has all null values for the delay type,
    then leave the values as is.
    Otherwise:

    - If the first delay is more than ``dist_threshold`` distance units
      from the first stop, then set the first stop delay to zero
      (charitably);
      otherwise set the first stop delay to the first delay.
    - If the last delay is more than ``dist_threshold`` distance units
      from the last stop, then set the last stop delay to zero
      (charitably); otherwise set the last stop delay to the last delay.
    - Linearly interpolate the remaining stop delays by distance.

    Return the resulting DataFrame.

    If a list of delay column names is given in ``delay_cols``,
    then alter those columns instead of the ``arrival_delay`` and
    ``departure_delay`` columns.
    This is useful if the given stop times have extra delay columns.

def interpolate_delays(augmented_stop_times, dist_threshold,
  delay_threshold=3600, delay_cols=None):
    """
    Given an augment stop times DataFrame as output by the function
    :func:`build_augmented_stop_times`, a distance threshold (float)
    in the same units as the ``'shape_dist_traveled'`` column of
    ``augmented_stop_times``, if that column is present, and a delay
    threshold (integer number of seconds), alter the delay values
    of the augmented stop times as follows.

    Drop all delays with absolute value more than ``delay_threshold``
    seconds.
    For each trip and for each delay type (arrival delay or departure
    delay) do the following.
    If the trip has all null values for the delay type,
    then leave the values as is.
    Otherwise:

    - If the first delay is more than ``dist_threshold`` distance units
      from the first stop, then set the first stop delay to zero
      (charitably);
      otherwise set the first stop delay to the first delay.
    - If the last delay is more than ``dist_threshold`` distance units
      from the last stop, then set the last stop delay to zero
      (charitably); otherwise set the last stop delay to the last delay.
    - Linearly interpolate the remaining stop delays by distance.

    Return the resulting DataFrame.

    If a list of delay column names is given in ``delay_cols``,
    then alter those columns instead of the ``arrival_delay`` and
    ``departure_delay`` columns.
    This is useful if the given stop times have extra delay columns.
    """
    f = augmented_stop_times.copy()

    if delay_cols is None or not set(delay_cols) <= set(f.columns):
        delay_cols = ['arrival_delay', 'departure_delay']

    # Return f if nothing to do
    if 'shape_dist_traveled' not in f.columns or\
      not f['shape_dist_traveled'].notnull().any() or\
      all([f[col].count() == f[col].shape[0] for col in delay_cols]):
        return f

    # Nullify fishy delays
    for col in delay_cols:
        f.loc[abs(f[col]) > delay_threshold, col] = np.nan

    # Fill null delays
    def fill(group):
        # Only columns that have at least one nonnull value.
        fill_cols = []
        for col in delay_cols:
            if group[col].count() >= 1:
                fill_cols.append(col)

        for col in fill_cols:
            # Set first and last delays
            for i in [0, -1]:
                j = group[col].dropna().index[i]
                dist_diff = (abs(group['shape_dist_traveled'].iat[i] -
                  group['shape_dist_traveled'].ix[j]))
                if dist_diff > dist_threshold:
                    group[col].iat[i] = 0
                else:
                    group[col].iat[i] = group[col].ix[j]

            # Interpolate remaining delays
            ind = np.where(group[col].notnull())[0]
            group[col] = np.interp(group['shape_dist_traveled'],
              group.iloc[ind]['shape_dist_traveled'],
              group.iloc[ind][col])

        return group

    f = f.groupby('trip_id').apply(fill)

    # Round
    f[delay_cols] = f[delay_cols].round(0)

    return f