From 9fe9735286c24f0f5a0e5c2567e7cbbdfd0b8750 Mon Sep 17 00:00:00 2001
From: tiffanychu90 <tiffany.ku@dot.ca.gov>
Date: Wed, 18 Dec 2024 01:26:34 +0000
Subject: [PATCH 1/5] (remove): no longer need filtering nearest 10 to 2 vp
 script

---
 rt_segment_speeds/scripts/vp_around_stops.py | 269 -------------------
 1 file changed, 269 deletions(-)
 delete mode 100644 rt_segment_speeds/scripts/vp_around_stops.py

diff --git a/rt_segment_speeds/scripts/vp_around_stops.py b/rt_segment_speeds/scripts/vp_around_stops.py
deleted file mode 100644
index 9f76378e0..000000000
--- a/rt_segment_speeds/scripts/vp_around_stops.py
+++ /dev/null
@@ -1,269 +0,0 @@
-"""
-Filter the nearest 10 neighbors down to the 
-nearest 2 neighbors for each stop position.
-Attach the projected stop position against shape,
-projected vp position against shape, and timestamps.
-"""
-import datetime
-import geopandas as gpd
-import pandas as pd
-import sys
-
-from dask import delayed, compute
-from loguru import logger
-from pathlib import Path
-from typing import Literal, Optional
-
-from segment_speed_utils import helpers
-from shared_utils import geo_utils
-from update_vars import SEGMENT_GCS, GTFS_DATA_DICT
-from segment_speed_utils.project_vars import SEGMENT_TYPES, PROJECT_CRS
-
-def stops_projected_against_shape(
-    input_file: str,
-    analysis_date: str, 
-    trip_stop_cols: list,
-) -> pd.DataFrame:
-    """
-    From nearest 10 vp points, project the stop geometry
-    onto shape geometry and get 
-    stop_meters.
-    """
-    stop_position = gpd.read_parquet(
-        f"{SEGMENT_GCS}{input_file}_{analysis_date}.parquet",
-        columns = trip_stop_cols + [
-            "shape_array_key", "stop_geometry"],
-    ).to_crs(PROJECT_CRS)
-    
-    shapes = helpers.import_scheduled_shapes(
-        analysis_date,
-        columns = ["shape_array_key", "geometry"],
-        crs = PROJECT_CRS,
-        get_pandas = True
-    )
-    
-    gdf = pd.merge(
-        stop_position,
-        shapes.rename(columns = {"geometry": "shape_geometry"}),
-        on = "shape_array_key",
-        how = "inner"
-    )
-    
-    gdf = gdf.assign(
-        stop_meters = gdf.shape_geometry.project(gdf.stop_geometry),
-    )[trip_stop_cols + ["stop_meters"]]
-    
-    del shapes, stop_position
-    
-    return gdf
-
-
-def explode_vp_nearest(
-    input_file: str,
-    analysis_date: str, 
-    trip_stop_cols: list,
-) -> pd.DataFrame:
-    """
-    Take nearest 10 vp, which holds vp_idx as an array,
-    and explode it so it becomes long.
-    """
-    vp_nearest = pd.read_parquet(
-        f"{SEGMENT_GCS}{input_file}_{analysis_date}.parquet",
-        columns = trip_stop_cols + [
-            "shape_array_key",
-            "nearest_vp_arr"],
-    ).explode(
-        "nearest_vp_arr"
-    ).drop_duplicates().reset_index(
-        drop=True
-    ).rename(
-        columns = {"nearest_vp_arr": "vp_idx"}
-    ).astype({"vp_idx": "int64"})
-    
-    return vp_nearest
-
-
-def get_vp_projected_against_shape(
-    input_file: str,
-    analysis_date: str, 
-    **kwargs
-) -> pd.DataFrame:
-    """
-    Put in subset of vp_idx (using the kwargs)
-    and turn the x, y into vp point geometry.
-    Merge in shapes and project the vp position
-    against shape geometry, and save out
-    shape_meters.
-    """    
-    # Get crosswalk of trip to shapes
-    trips_to_shapes = helpers.import_scheduled_trips(
-        analysis_date,
-        columns = ["trip_instance_key", "shape_array_key"],
-        get_pandas = True
-    )
-    
-    # Get shapes
-    shapes = helpers.import_scheduled_shapes(
-        analysis_date,
-        columns = ["shape_array_key", "geometry"],
-        crs = PROJECT_CRS,
-        get_pandas = True
-    )
-    
-    # Subset usable vp with only the ones present in exploded vp
-    # and turn those into vp geometry
-    vp = pd.read_parquet(
-        f"{SEGMENT_GCS}{input_file}_{analysis_date}",
-        columns = ["trip_instance_key", "vp_idx", "x", "y"],
-        **kwargs
-    ).pipe(geo_utils.vp_as_gdf, crs = PROJECT_CRS)
-    
-    # Merge all together so we can project vp point goem
-    # against shape line geom
-    gdf = pd.merge(
-        vp.rename(columns = {"geometry": "vp_geometry"}),
-        trips_to_shapes,
-        on = "trip_instance_key",
-        how = "inner"
-    ).merge(
-        shapes.rename(columns = {"geometry": "shape_geometry"}),
-        on = "shape_array_key",
-        how = "inner"
-    ).set_geometry("vp_geometry")
-
-    del trips_to_shapes, shapes, vp
-    
-    gdf = gdf.assign(
-        shape_meters = gdf.shape_geometry.project(gdf.vp_geometry),
-    )[["vp_idx", "shape_meters"]]
-            
-    return gdf
-
-
-def find_two_closest_vp(
-    df: pd.DataFrame, 
-    group_cols: list
-) -> pd.DataFrame:
-    """
-    Based on the distances calculated between vp and stop, 
-    keep the 2 observations that are closest. Find the smallest
-    positive distance and negative distance.
-    
-    This filters down the nearest 10 into nearest 2.
-    """
-    positive_distances_df = df.loc[df.stop_vp_distance_meters >= 0]
-    negative_distances_df = df.loc[df.stop_vp_distance_meters < 0]
-    
-    #https://github.com/pandas-dev/pandas/issues/45089
-    # add dropna=False or else too many combos are lost
-    min_pos_distance = (
-        positive_distances_df
-        .groupby(group_cols, 
-                 observed=True, group_keys=False, dropna=False)
-        .agg({"stop_vp_distance_meters": "min"})
-        .reset_index()
-    )
-    
-    min_neg_distance = (
-        negative_distances_df
-        .groupby(group_cols, 
-                 observed=True, group_keys=False, dropna=False)
-        .agg({"stop_vp_distance_meters": "max"})
-        .reset_index()
-    )
-
-    two_vp = pd.concat(
-        [min_pos_distance, min_neg_distance], 
-        axis=0, ignore_index=True
-    )
-    
-    return two_vp
-
-    
-def filter_to_nearest_two_vp(
-    analysis_date: str,
-    segment_type: Literal[SEGMENT_TYPES],
-    config_path: Optional[Path] = GTFS_DATA_DICT
-):
-    dict_inputs = config_path[segment_type]
-    trip_stop_cols = [*dict_inputs["trip_stop_cols"]]
-    USABLE_VP_FILE = dict_inputs["stage1"]
-    INPUT_FILE = dict_inputs["stage2"]
-    EXPORT_FILE = dict_inputs["stage2b"]
-
-    start = datetime.datetime.now()
-    
-    stop_meters_df = delayed(stops_projected_against_shape)(
-        INPUT_FILE, analysis_date, trip_stop_cols)
-    
-    vp_nearest = delayed(explode_vp_nearest)(
-        INPUT_FILE, analysis_date, trip_stop_cols)
-
-    subset_vp = vp_nearest.vp_idx.unique()
-        
-    vp_meters_df = delayed(get_vp_projected_against_shape)(
-        USABLE_VP_FILE,
-        analysis_date, 
-        filters = [[("vp_idx", "in", subset_vp)]]
-    )
-    
-    gdf = delayed(pd.merge)(
-        vp_nearest,
-        stop_meters_df,
-        on = trip_stop_cols,
-        how = "inner"
-    ).merge(
-        vp_meters_df,
-        on = "vp_idx",
-        how = "inner"
-    )
-    
-    # Calculate the distance between the stop and vp position
-    # This is used to find the minimum positive and minimum negative
-    # distance (get at vp before and after stop)
-    gdf = gdf.assign(
-        stop_meters = gdf.stop_meters.round(3),
-        shape_meters = gdf.shape_meters.round(3),
-        stop_vp_distance_meters = (gdf.stop_meters - gdf.shape_meters).round(2)
-    )
-    
-    gdf_filtered = delayed(find_two_closest_vp)(gdf, trip_stop_cols)
-    
-    gdf2 = delayed(pd.merge)(
-        gdf,
-        gdf_filtered,
-        on = trip_stop_cols + ["stop_vp_distance_meters"],
-        how = "inner"
-    )
-    
-    gdf2 = compute(gdf2)[0]
-    
-    del gdf, gdf_filtered, vp_nearest, stop_meters_df, vp_meters_df
-        
-    gdf2.to_parquet(
-        f"{SEGMENT_GCS}{EXPORT_FILE}_{analysis_date}.parquet",
-    )
-        
-    end = datetime.datetime.now()
-    logger.info(f"nearest 2 vp for {segment_type} "
-                f"{analysis_date}: {end - start}")
-    
-    del gdf2
-    
-    return
-    
-'''
-if __name__ == "__main__":
-    
-    from segment_speed_utils.project_vars import analysis_date_list
-    
-    delayed_dfs = [
-        delayed(filter_to_nearest_two_vp)(
-            analysis_date = analysis_date,
-            segment_type = segment_type,
-            config_path = GTFS_DATA_DICT
-        ) for analysis_date in analysis_date_list
-    ]
-
-    [compute(i)[0] for i in delayed_dfs]
-'''
\ No newline at end of file

From 873a97cc80e6aacade880757160a74cb21c6c59e Mon Sep 17 00:00:00 2001
From: tiffanychu90 <tiffany.ku@dot.ca.gov>
Date: Wed, 18 Dec 2024 01:31:04 +0000
Subject: [PATCH 2/5] (remove): array_utils, move into segment_calcs

---
 .../segment_speed_utils/array_utils.py        | 37 -------------------
 1 file changed, 37 deletions(-)
 delete mode 100644 rt_segment_speeds/segment_speed_utils/array_utils.py

diff --git a/rt_segment_speeds/segment_speed_utils/array_utils.py b/rt_segment_speeds/segment_speed_utils/array_utils.py
deleted file mode 100644
index 295aae47e..000000000
--- a/rt_segment_speeds/segment_speed_utils/array_utils.py
+++ /dev/null
@@ -1,37 +0,0 @@
-"""
-Functions for working with numpy arrays.
-"""
-import numpy as np
-import pandas as pd
-
-from numba import jit
-
-def rolling_window_make_array(
-    df: pd.DataFrame, 
-    window: int, 
-    rolling_col: str
-) -> pd.DataFrame:
-    # https://stackoverflow.com/questions/47482009/pandas-rolling-window-to-return-an-array
-    df[f"rolling_{rolling_col}"] = [
-        np.asarray(window) for window in 
-        df.groupby("trip_instance_key")[rolling_col].rolling(
-            window = window, center=True)
-    ]
-    
-    is_monotonic_series = np.vectorize(monotonic_check)(df[f"rolling_{rolling_col}"])
-    df[f"{rolling_col}_monotonic"] = is_monotonic_series
-    
-    return df
-    
-@jit(nopython=True)
-def monotonic_check(arr: np.ndarray) -> bool:
-    """
-    For an array, check if it's monotonically increasing. 
-    https://stackoverflow.com/questions/4983258/check-list-monotonicity
-    """
-    diff_arr = np.diff(arr)
-    
-    if np.all(diff_arr > 0):
-        return True
-    else:
-        return False
\ No newline at end of file

From 8255c381d4c4eca68af3404540923cbaa0026edd Mon Sep 17 00:00:00 2001
From: tiffanychu90 <tiffany.ku@dot.ca.gov>
Date: Wed, 18 Dec 2024 01:31:37 +0000
Subject: [PATCH 3/5] (segment_speed_utils): remove unused functions

---
 .../segment_speed_utils/segment_calcs.py      | 31 ++++++++++++++++++-
 .../segment_speed_utils/vp_transform.py       | 17 +---------
 2 files changed, 31 insertions(+), 17 deletions(-)

diff --git a/rt_segment_speeds/segment_speed_utils/segment_calcs.py b/rt_segment_speeds/segment_speed_utils/segment_calcs.py
index 9bfac4613..004490ece 100644
--- a/rt_segment_speeds/segment_speed_utils/segment_calcs.py
+++ b/rt_segment_speeds/segment_speed_utils/segment_calcs.py
@@ -1,3 +1,6 @@
+"""
+Functions related to calculating segment speeds.
+"""
 import dask.dataframe as dd
 import dask_geopandas as dg
 import geopandas as gpd
@@ -149,4 +152,30 @@ def interpolate_stop_arrival_time(
 
     return np.interp(
         stop_position, np.asarray(shape_meters_arr), timestamp_arr
-    ).astype("datetime64[s]")
\ No newline at end of file
+    ).astype("datetime64[s]")
+
+
+def rolling_window_make_array(
+    df: pd.DataFrame, 
+    window: int, 
+    rolling_col: str
+) -> pd.DataFrame:
+    """
+    Interpolated stop arrival times are checked
+    to see if they are monotonically increasing.
+    If it isn't, it gets recalculated based on 
+    stop_meters (the stop's position) relative to
+    other stop arrival times.
+    
+    https://stackoverflow.com/questions/47482009/pandas-rolling-window-to-return-an-array
+    """
+    df[f"rolling_{rolling_col}"] = [
+        np.asarray(window) for window in 
+        df.groupby("trip_instance_key")[rolling_col].rolling(
+            window = window, center=True)
+    ]
+    
+    is_monotonic_series = np.vectorize(monotonic_check)(df[f"rolling_{rolling_col}"])
+    df[f"{rolling_col}_monotonic"] = is_monotonic_series
+    
+    return df
\ No newline at end of file
diff --git a/rt_segment_speeds/segment_speed_utils/vp_transform.py b/rt_segment_speeds/segment_speed_utils/vp_transform.py
index 48694b585..bc41a9a90 100644
--- a/rt_segment_speeds/segment_speed_utils/vp_transform.py
+++ b/rt_segment_speeds/segment_speed_utils/vp_transform.py
@@ -53,19 +53,4 @@ def condense_point_geom_to_line(
            .reset_index()
           )
     
-    return df3
-
-
-def sort_by_vp_idx_order(
-    vp_idx_array: np.ndarray, 
-    geometry_array: np.ndarray,
-    timestamp_array: np.ndarray,
-) -> tuple[np.ndarray]:    
-    
-    sort_order = np.argsort(vp_idx_array, axis=0)
-    
-    vp_sorted = np.take_along_axis(vp_idx_array, sort_order, axis=0)
-    geom_sorted = np.take_along_axis(geometry_array, sort_order, axis=0)
-    timestamp_sorted = np.take_along_axis(timestamp_array, sort_order, axis=0)
-    
-    return vp_sorted, geom_sorted, timestamp_sorted
\ No newline at end of file
+    return df3
\ No newline at end of file

From d17d8b309507f3447c7746d20e4af8ca30ac89d1 Mon Sep 17 00:00:00 2001
From: tiffanychu90 <tiffany.ku@dot.ca.gov>
Date: Wed, 18 Dec 2024 01:31:55 +0000
Subject: [PATCH 4/5] update init after array_utils removed

---
 rt_segment_speeds/segment_speed_utils/__init__.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/rt_segment_speeds/segment_speed_utils/__init__.py b/rt_segment_speeds/segment_speed_utils/__init__.py
index 4d10f7d0e..d7e66563e 100644
--- a/rt_segment_speeds/segment_speed_utils/__init__.py
+++ b/rt_segment_speeds/segment_speed_utils/__init__.py
@@ -1,5 +1,4 @@
 from . import (
-    array_utils,
     gtfs_schedule_wrangling,
     helpers,
     metrics,
@@ -12,7 +11,6 @@
 )
 
 __all__ = [
-    "array_utils",
     "gtfs_schedule_wrangling",
     "helpers",
     "metrics",

From bc6469835495caade0d52dea926bfdf0ab74ca58 Mon Sep 17 00:00:00 2001
From: tiffanychu90 <tiffany.ku@dot.ca.gov>
Date: Wed, 18 Dec 2024 16:17:55 +0000
Subject: [PATCH 5/5] rename functions for clarity, test sep-nov2024 dates for
 nearest vp step

---
 .../scripts/interpolate_stop_arrival.py       |  5 +-
 .../scripts/nearest_vp_to_stop.py             |  3 +-
 .../scripts/pipeline_rt_stop_times.py         |  5 +-
 .../scripts/pipeline_segment_speeds.py        |  7 +--
 .../scripts/pipeline_speedmap.py              |  7 +--
 .../segment_speed_utils/neighbor.py           | 54 +++++++++----------
 .../segment_speed_utils/segment_calcs.py      | 14 +++++
 7 files changed, 54 insertions(+), 41 deletions(-)

diff --git a/rt_segment_speeds/scripts/interpolate_stop_arrival.py b/rt_segment_speeds/scripts/interpolate_stop_arrival.py
index fafb920b8..0495893e9 100644
--- a/rt_segment_speeds/scripts/interpolate_stop_arrival.py
+++ b/rt_segment_speeds/scripts/interpolate_stop_arrival.py
@@ -17,8 +17,7 @@
 from pathlib import Path
 from typing import Literal, Optional
 
-from segment_speed_utils import (array_utils, helpers, 
-                                 segment_calcs)
+from segment_speed_utils import helpers, segment_calcs
 from update_vars import SEGMENT_GCS, GTFS_DATA_DICT
 from segment_speed_utils.project_vars import PROJECT_CRS, SEGMENT_TYPES
 from shared_utils import rt_dates
@@ -166,7 +165,7 @@ def enforce_monotonicity_and_interpolate_across_stops(
     df = segment_calcs.convert_timestamp_to_seconds(
         df, ["arrival_time"])
 
-    df = array_utils.rolling_window_make_array(
+    df = segment_calcs.rolling_window_make_array(
         df, 
         window = 3, rolling_col = "arrival_time_sec"
     )
diff --git a/rt_segment_speeds/scripts/nearest_vp_to_stop.py b/rt_segment_speeds/scripts/nearest_vp_to_stop.py
index 2ea13aae0..6706637b4 100644
--- a/rt_segment_speeds/scripts/nearest_vp_to_stop.py
+++ b/rt_segment_speeds/scripts/nearest_vp_to_stop.py
@@ -4,6 +4,7 @@
 """
 import datetime
 import geopandas as gpd
+import numpy as np
 import pandas as pd
 import sys
 
@@ -158,7 +159,7 @@ def nearest_neighbor_for_stop(
     gdf = neighbor.merge_stop_vp_for_nearest_neighbor(stop_times, analysis_date)
     
     vp_before, vp_after, vp_before_meters, vp_after_meters = np.vectorize(
-        neighbor.subset_arrays_to_valid_directions
+        neighbor.two_nearest_neighbor_near_stop
     )(
         gdf.vp_primary_direction, 
         gdf.vp_geometry, 
diff --git a/rt_segment_speeds/scripts/pipeline_rt_stop_times.py b/rt_segment_speeds/scripts/pipeline_rt_stop_times.py
index a37199c96..39cc7567e 100644
--- a/rt_segment_speeds/scripts/pipeline_rt_stop_times.py
+++ b/rt_segment_speeds/scripts/pipeline_rt_stop_times.py
@@ -38,7 +38,7 @@
         
     logger.remove()
 
-
+    
     LOG_FILE = "../logs/interpolate_stop_arrival.log"
     logger.add(LOG_FILE, retention="3 months")
     logger.add(sys.stderr, 
@@ -57,7 +57,7 @@
 
     logger.remove()
 
-
+    
     LOG_FILE = "../logs/speeds_by_segment_trip.log"
     logger.add(LOG_FILE, retention="3 months")
     logger.add(sys.stderr, 
@@ -75,3 +75,4 @@
     [compute(i)[0] for i in delayed_dfs]
 
     logger.remove()
+    
\ No newline at end of file
diff --git a/rt_segment_speeds/scripts/pipeline_segment_speeds.py b/rt_segment_speeds/scripts/pipeline_segment_speeds.py
index fe8084eba..fa7922c14 100644
--- a/rt_segment_speeds/scripts/pipeline_segment_speeds.py
+++ b/rt_segment_speeds/scripts/pipeline_segment_speeds.py
@@ -17,7 +17,7 @@
 if __name__ == "__main__":
     
     from segment_speed_utils.project_vars import analysis_date_list    
-
+    
     segment_type = "stop_segments"
     print(f"segment_type: {segment_type}")
     
@@ -59,7 +59,7 @@
 
     logger.remove()
 
-
+    
     LOG_FILE = "../logs/speeds_by_segment_trip.log"
     logger.add(LOG_FILE, retention="3 months")
     logger.add(sys.stderr, 
@@ -76,4 +76,5 @@
 
     [compute(i)[0] for i in delayed_dfs]
 
-    logger.remove()
\ No newline at end of file
+    logger.remove()
+    
\ No newline at end of file
diff --git a/rt_segment_speeds/scripts/pipeline_speedmap.py b/rt_segment_speeds/scripts/pipeline_speedmap.py
index 293db2e44..86997ebfe 100644
--- a/rt_segment_speeds/scripts/pipeline_speedmap.py
+++ b/rt_segment_speeds/scripts/pipeline_speedmap.py
@@ -90,7 +90,7 @@ def concatenate_speedmap_proxy_arrivals_with_remaining(
 
     logger.remove()
 
-
+    
     LOG_FILE = "../logs/interpolate_stop_arrival.log"
     logger.add(LOG_FILE, retention="3 months")
     logger.add(sys.stderr, 
@@ -109,7 +109,7 @@ def concatenate_speedmap_proxy_arrivals_with_remaining(
 
     logger.remove()
     
-        
+    
     t0 = datetime.datetime.now()
     delayed_dfs = [
         delayed(concatenate_speedmap_proxy_arrivals_with_remaining)(
@@ -139,4 +139,5 @@ def concatenate_speedmap_proxy_arrivals_with_remaining(
 
     [compute(i)[0] for i in delayed_dfs]
     
-    logger.remove()
\ No newline at end of file
+    logger.remove()
+    
\ No newline at end of file
diff --git a/rt_segment_speeds/segment_speed_utils/neighbor.py b/rt_segment_speeds/segment_speed_utils/neighbor.py
index f5d21ce84..2e3f9a307 100644
--- a/rt_segment_speeds/segment_speed_utils/neighbor.py
+++ b/rt_segment_speeds/segment_speed_utils/neighbor.py
@@ -15,7 +15,7 @@ def merge_stop_vp_for_nearest_neighbor(
     stop_times: gpd.GeoDataFrame,
     analysis_date: str,
     **kwargs
-):
+) -> gpd.GeoDataFrame:
     """
     Merge stop times file with vp.
     vp gdf has been condensed so that all the vp coords
@@ -64,7 +64,7 @@ def merge_stop_vp_for_nearest_neighbor(
 
 
 def find_nearest_points(
-    vp_coords_line: np.ndarray, 
+    vp_coords_array: np.ndarray, 
     target_stop: shapely.Point, 
     vp_idx_array: np.ndarray,
 ) -> np.ndarray:
@@ -79,7 +79,7 @@ def find_nearest_points(
     (ex: nearest 5 vp to each stop).
     """
     indices = geo_utils.nearest_snap(
-        vp_coords_line, 
+        vp_coords_array, 
         target_stop, 
         k_neighbors = 5
     )
@@ -89,17 +89,17 @@ def find_nearest_points(
     # if we want 10 nearest neighbors and 8th, 9th, 10th are all
     # the same result, the 8th will have a result, then 9th and 10th will
     # return the length of the array (which is out-of-bounds)
+    # using vp_coords_array keeps too many points (is this because coords can be dupes?)
     indices2 = indices[indices < vp_idx_array.size]
     
     return indices2
 
 
 def filter_to_nearest2_vp(
-    vp_coords_line: np.ndarray,
+    nearest_vp_coords_array: np.ndarray,
     shape_geometry: shapely.LineString,
-    vp_idx_array: np.ndarray,
+    nearest_vp_idx_array: np.ndarray,
     stop_meters: float,
-    indices_of_nearest: np.ndarray,
 ) -> tuple[np.ndarray]:
     """
     Take the indices that are the nearest.
@@ -109,16 +109,11 @@ def filter_to_nearest2_vp(
     Filter down to the nearest 2 vp before and after a stop.
     If there isn't one before or after, a value of -1 is returned.
     """
-    # Subset the array of vp coords and vp_idx_array with 
-    # the indices that show the nearest k neighbors.
-    nearest_vp = vp_coords_line[indices_of_nearest]
-    nearest_vp_idx = vp_idx_array[indices_of_nearest]
-    
     # Project these vp coords to shape geometry and see how far it is
     # from the stop's position on the shape
     nearest_vp_projected = np.asarray(
         [shape_geometry.project(shapely.Point(i)) 
-         for i in nearest_vp]
+         for i in nearest_vp_coords_array]
     )
 
     # Negative values are before the stop
@@ -126,26 +121,28 @@ def filter_to_nearest2_vp(
     before_indices = np.where(nearest_vp_projected - stop_meters < 0)[0]
     after_indices = np.where(nearest_vp_projected - stop_meters > 0)[0]
     
+    # Set missing values when we're not able to find a nearest neighbor result
+    # use -1 as vp_idx (since this is not present in vp_usable)
+    # and zeroes for meters
+    before_idx = -1
+    after_idx = -1
+    before_vp_meters = 0
+    after_vp_meters = 0
+    
     # Grab the closest vp before a stop (-1 means array was empty)
     if before_indices.size > 0:
-        before_idx = nearest_vp_idx[before_indices][-1] 
+        before_idx = nearest_vp_idx_array[before_indices][-1] 
         before_vp_meters = nearest_vp_projected[before_indices][-1]
-    else:
-        before_idx = -1
-        before_vp_meters = 0
-        
+   
     # Grab the closest vp after a stop (-1 means array was empty)
     if after_indices.size > 0:
-        after_idx = nearest_vp_idx[after_indices][0]
+        after_idx = nearest_vp_idx_array[after_indices][0]
         after_vp_meters = nearest_vp_projected[after_indices][0]
-    else:
-        after_idx = -1
-        after_vp_meters = 0
     
     return before_idx, after_idx, before_vp_meters, after_vp_meters
 
 
-def subset_arrays_to_valid_directions(
+def two_nearest_neighbor_near_stop(
     vp_direction_array: np.ndarray,
     vp_geometry: shapely.LineString,
     vp_idx_array: np.ndarray,
@@ -170,23 +167,22 @@ def subset_arrays_to_valid_directions(
     valid_indices = (vp_direction_array != opposite_direction).nonzero()   
 
     # These are vp coords where index values of opposite direction is excluded
-    valid_vp_coords_line = np.array(vp_geometry.coords)[valid_indices]
+    valid_vp_coords_array = np.array(vp_geometry.coords)[valid_indices]
     
     # These are the subset of vp_idx values where opposite direction is excluded
-    valid_vp_idx_arr = np.asarray(vp_idx_array)[valid_indices]  
+    valid_vp_idx_array = np.asarray(vp_idx_array)[valid_indices]  
     
     nearest_indices = find_nearest_points(
-        valid_vp_coords_line, 
+        valid_vp_coords_array, 
         stop_geometry, 
-        valid_vp_idx_arr,
+        valid_vp_idx_array,
     )
  
     before_vp, after_vp, before_meters, after_meters = filter_to_nearest2_vp(
-        valid_vp_coords_line,
+        valid_vp_idx_array[nearest_indices], # subset of coords in nn
         shape_geometry,
-        valid_vp_idx_arr,
+        valid_vp_idx_array[nearest_indices], # subset of vp_idx in nn
         stop_meters,
-        nearest_indices,
     )
     
     return before_vp, after_vp, before_meters, after_meters
\ No newline at end of file
diff --git a/rt_segment_speeds/segment_speed_utils/segment_calcs.py b/rt_segment_speeds/segment_speed_utils/segment_calcs.py
index 004490ece..4fb21c5b8 100644
--- a/rt_segment_speeds/segment_speed_utils/segment_calcs.py
+++ b/rt_segment_speeds/segment_speed_utils/segment_calcs.py
@@ -7,6 +7,7 @@
 import numpy as np
 import pandas as pd
 
+from numba import jit
 from typing import Union
 
 from shared_utils.rt_utils import MPH_PER_MPS
@@ -155,6 +156,19 @@ def interpolate_stop_arrival_time(
     ).astype("datetime64[s]")
 
 
+@jit(nopython=True)
+def monotonic_check(arr: np.ndarray) -> bool:
+    """
+    For an array, check if it's monotonically increasing. 
+    https://stackoverflow.com/questions/4983258/check-list-monotonicity
+    """
+    diff_arr = np.diff(arr)
+    
+    if np.all(diff_arr > 0):
+        return True
+    else:
+        return False
+    
 def rolling_window_make_array(
     df: pd.DataFrame, 
     window: int,