directional headway consistency

notebooks/4_validate_gtfs_io.ipynb

@@ -11,7 +11,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "id": "a691e35d",
    "metadata": {},
    "outputs": [],
@@ -34,7 +34,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 2,
    "id": "cdc591c4",
    "metadata": {},
    "outputs": [
@@ -90,8 +90,7 @@
       "Route 11963: Round-trip 259.9min exceeds limit (240.0min), filtered out\n",
       "Filtered out 39 routes (excessive round-trip time)\n",
       "Failed to process 1 routes (no valid data)\n",
-      "Found 12 cells with >10min mapping difference\n",
-      "Large discrepancy between raw and discretized fleet calculations!\n"
+      "Found 12 cells with >10min mapping difference\n"
      ]
     },
     {
@@ -128,7 +127,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "id": "35900e19",
    "metadata": {},
    "outputs": [
@@ -170,7 +169,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "id": "079bea2b",
    "metadata": {},
    "outputs": [
@@ -311,7 +310,7 @@
        "60   12627             338                  240    -98         98"
       ]
      },
-     "execution_count": 5,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -355,7 +354,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 5,
    "id": "6638e32d",
    "metadata": {},
    "outputs": [
@@ -437,7 +436,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 6,
    "id": "8db79bae",
    "metadata": {},
    "outputs": [
@@ -525,7 +524,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 7,
    "id": "bde2e350",
    "metadata": {},
    "outputs": [

src/transit_opt/gtfs/solution_manager.py

@@ -86,6 +86,7 @@ def _calculate_solution_fleet_rows(self, solution_id: str, solution_data: dict[s
             no_service_threshold=params.get("no_service_threshold", 480),
             allowed_headways=np.array(allowed_headways),
             no_service_index=no_service_idx,
+            n_directions=self.optimization_data.get("routes", {}).get("n_directions", None),
         )
 
         rows = []

src/transit_opt/optimisation/problems/base.py

@@ -140,12 +140,13 @@ def _calculate_fleet_from_solution(self, solution_matrix: np.ndarray) -> np.ndar
 
         # Use shared calculation logic
         fleet_results = calculate_fleet_requirements(
-            headways_matrix=solution_matrix,
+            headways_matrix=pt_matrix,
             round_trip_times=self.round_trip_times,
             operational_buffer=operational_buffer,
             no_service_threshold=no_service_threshold,
             allowed_headways=self.allowed_headways,
             no_service_index=self.no_service_index,
+            n_directions=self.opt_data.get("routes", {}).get("n_directions", None),
         )
 
         return fleet_results["fleet_per_interval"]

src/transit_opt/optimisation/utils/fleet_calculations.py

@@ -5,6 +5,7 @@
 GTFSDataPreparator and optimization constraint handlers to ensure consistency.
 """
 
+import numbers
 from typing import Any
 
 import numpy as np
@@ -17,6 +18,7 @@ def calculate_fleet_requirements(
     no_service_threshold: float = 480,
     allowed_headways: np.ndarray | None = None,
     no_service_index: int | None = None,
+    n_directions: np.ndarray | None = None,
 ) -> dict[str, Any]:
     """
     Unified fleet calculation for both baseline analysis and optimization constraints.
@@ -45,9 +47,10 @@ def calculate_fleet_requirements(
 
     # Validate inputs
     if len(round_trip_times) != n_routes:
-        raise ValueError(
-            f"round_trip_times length ({len(round_trip_times)}) must match n_routes ({n_routes})"
-        )
+        raise ValueError(f"round_trip_times length ({len(round_trip_times)}) must match n_routes ({n_routes})")
+
+    if n_directions is not None and len(n_directions) != n_routes:
+        raise ValueError(f"n_directions length ({len(n_directions)}) must match n_routes ({n_routes})")
 
     # Initialize output arrays
     route_fleet_matrix = np.zeros((n_routes, n_intervals), dtype=int)
@@ -57,15 +60,20 @@ def calculate_fleet_requirements(
     for route_idx in range(n_routes):
         round_trip_time = round_trip_times[route_idx]
 
+        raw_dirs = n_directions[route_idx] if n_directions is not None else 2.0
+        try:
+            dirs = float(raw_dirs)
+            dirs = 2.0 if np.isnan(dirs) or np.isinf(dirs) else max(1.0, dirs)
+        except (ValueError, TypeError):
+            dirs = 2.0
+
         for interval_idx in range(n_intervals):
             headway_value = headways_matrix[route_idx, interval_idx]
 
             # Decode headway value based on context
             if allowed_headways is not None and no_service_index is not None:
                 # Optimization context: decode choice index to headway value
-                if isinstance(headway_value, (int, np.integer)) and headway_value < len(
-                    allowed_headways
-                ):
+                if isinstance(headway_value, numbers.Integral) and headway_value < len(allowed_headways):
                     if headway_value == no_service_index:
                         actual_headway = np.inf  # No service
                     else:
@@ -77,15 +85,11 @@ def calculate_fleet_requirements(
                 actual_headway = headway_value
 
             # Calculate vehicles needed using standardized logic
-            if (
-                not np.isnan(actual_headway)
-                and not np.isinf(actual_headway)
-                and actual_headway < no_service_threshold
-            ):
-                # Valid service headway - apply same formula as GTFSDataPreparator
-                vehicles_needed = np.ceil(
-                    (round_trip_time * operational_buffer) / actual_headway
-                )
+            if not np.isnan(actual_headway) and not np.isinf(actual_headway) and actual_headway < no_service_threshold:
+                # Using Little's Law: We derive directional headway so it correctly
+                # scales against round_trip_time without double counting across forks.
+                directional_headway = actual_headway * dirs
+                vehicles_needed = np.ceil((round_trip_time * operational_buffer) / directional_headway)
                 vehicles_needed = max(1, int(vehicles_needed))  # At least 1 vehicle
             else:
                 # No service or invalid headway

src/transit_opt/preprocessing/prepare_gtfs.py

@@ -370,6 +370,7 @@ def extract_optimization_data(self, allowed_headways: list[float]) -> dict[str,
                 "ids": route_ids,
                 "round_trip_times": round_trip_times,
                 "current_headways": current_headways,
+                "n_directions": np.array([r["n_directions"] for r in route_data]),
             },
             "constraints": {
                 "fleet_analysis": fleet_analysis,
@@ -453,8 +454,8 @@ def _extract_route_essentials(self) -> list[dict[str, Any]]:
                 failed_count += 1
                 continue
 
-            # Calculate round-trip time
-            round_trip_time = self._calculate_round_trip_time(route_id, route_trips)
+            # Calculate round-trip time and number of directions
+            round_trip_time, n_directions = self._calculate_round_trip_time(route_id, route_trips)
 
             # Track default usage
             if round_trip_time == self.default_round_trip_time:
@@ -481,6 +482,7 @@ def _extract_route_essentials(self) -> list[dict[str, Any]]:
                     "route_id": route_id,
                     "headways_by_interval": headways_by_interval,
                     "round_trip_time": round_trip_time,
+                    "n_directions": n_directions,
                 }
             )
 
@@ -585,7 +587,7 @@ def _calculate_route_headways(self, route_id: str, route_trips: pd.DataFrame) ->
             logger.debug(f"Route {route_id}: Exception in headway calculation: {e}")
             return headways
 
-    def _calculate_round_trip_time(self, route_id: str, route_trips: pd.DataFrame) -> float:
+    def _calculate_round_trip_time(self, route_id: str, route_trips: pd.DataFrame) -> tuple[float, int]:
         """
         Calculate round-trip time with turnaround buffer for fleet sizing.
 
@@ -622,7 +624,7 @@ def _calculate_round_trip_time(self, route_id: str, route_trips: pd.DataFrame) -
 
             if len(route_stop_times) == 0:
                 logger.debug(f"Route {route_id}: No stop times, using default {self.default_round_trip_time}min")
-                return self.default_round_trip_time
+                return self.default_round_trip_time, 2
 
             trip_durations = []
             for trip_id, trip_stops in route_stop_times.groupby("trip_id"):
@@ -654,6 +656,10 @@ def _calculate_round_trip_time(self, route_id: str, route_trips: pd.DataFrame) -
                 else:
                     n_directions = 2  # Fallback to standard 2-way assumption if no info
 
+                # Normalize directions to avoid division-by-zero or massive fleet inflation
+                # from fragmented headsigns. Assume it's either a 1-way loop or a 2-way route.
+                n_directions = 1 if n_directions == 1 else 2
+
                 # If we have exactly 1 direction, we treat it as a loop (1x).
                 # All other cases (2, 8, etc.) we treat as needing return trips (2x).
                 multiplier = 1.0 if n_directions == 1 else 2.0
@@ -665,17 +671,17 @@ def _calculate_round_trip_time(self, route_id: str, route_trips: pd.DataFrame) -
                     f"directions: {n_directions}, multiplier: {multiplier}, "
                     f"buffer: {self.turnaround_buffer})"
                 )
-                return round_trip
+                return round_trip, n_directions
             else:
                 logger.debug(f"Route {route_id}: No valid durations, using default {self.default_round_trip_time}min")
-                return self.default_round_trip_time
+                return self.default_round_trip_time, 2
 
         except Exception as e:
             logger.debug(
                 f"Route {route_id}: Exception calculating round-trip time: {e}, "
                 f"using default {self.default_round_trip_time}min"
             )
-            return self.default_round_trip_time
+            return self.default_round_trip_time, 2
 
     def _create_initial_solution(self, current_headways: np.ndarray, headway_to_index: dict[float, int]) -> np.ndarray:
         """
@@ -859,12 +865,14 @@ def _analyze_current_fleet(self, route_data: list[dict[str, Any]]) -> dict[str,
                 * fleet_distribution: Count of routes by fleet size category
 
         **Fleet Calculation Formula**:
-            vehicles_needed = ceil((round_trip_time * operational_buffer) / headway)
+            directional_headway = headway * n_directions
+            vehicles_needed = ceil((round_trip_time * operational_buffer) / directional_headway)
 
             Where:
             - round_trip_time: Total time for vehicle to complete route and return (minutes)
             - operational_buffer: Extra time factor for maintenance, delays, crew relief (1.15 = 15%)
-            - headway: Time between consecutive departures (minutes)
+            - headway: Time between consecutive departures across all directions (aggregate headway, minutes)
+            - n_directions: Number of directions served (1 or 2, limits double counting)
             - ceil(): Round up to next integer (can't have fractional vehicles)
 
         **Calculation Examples**:
@@ -899,6 +907,7 @@ def _analyze_current_fleet(self, route_data: list[dict[str, Any]]) -> dict[str,
         # Extract data for calculation
         round_trip_times = np.array([r["round_trip_time"] for r in route_data])
         raw_headways_matrix = np.array([r["headways_by_interval"] for r in route_data])
+        n_directions = np.array([r["n_directions"] for r in route_data])
 
         # CALCULATION 1: Raw GTFS headways (original baseline)
         logger.debug("Calculating fleet with raw GTFS headways...")
@@ -908,6 +917,7 @@ def _analyze_current_fleet(self, route_data: list[dict[str, Any]]) -> dict[str,
             round_trip_times=round_trip_times,
             operational_buffer=operational_buffer,
             no_service_threshold=self.no_service_threshold_minutes,
+            n_directions=n_directions,
         )
 
         # CALCULATION 2: Discretized headways (constraint-consistent baseline)
@@ -932,6 +942,7 @@ def _analyze_current_fleet(self, route_data: list[dict[str, Any]]) -> dict[str,
             round_trip_times=round_trip_times,
             operational_buffer=operational_buffer,
             no_service_threshold=self.no_service_threshold_minutes,
+            n_directions=n_directions,
         )
 
         # Extract results (use discretized for optimization, keep raw for reporting)

tests/fixtures/__init__.py

@@ -139,8 +139,7 @@ def precalculated_fleet_data(sample_optimization_data):
     print(f"      Fleet by interval: {baseline_data['current_fleet_by_interval']}")
 
     # Get parameters for fleet calculations (same as GTFSDataPreparator used)
-    from transit_opt.optimisation.utils.fleet_calculations import \
-        calculate_fleet_requirements
+    from transit_opt.optimisation.utils.fleet_calculations import calculate_fleet_requirements
 
     allowed_headways = sample_optimization_data["allowed_headways"]
     round_trip_times = sample_optimization_data["routes"]["round_trip_times"]
@@ -173,10 +172,7 @@ def precalculated_fleet_data(sample_optimization_data):
         # Convert solution indices to actual headway minutes
         # This is the key step - solution_matrix contains indices, we need actual headway values
         headways_matrix = np.array(
-            [
-                [allowed_headways[solution_matrix[i, j]] for j in range(n_intervals)]
-                for i in range(n_routes)
-            ]
+            [[allowed_headways[solution_matrix[i, j]] for j in range(n_intervals)] for i in range(n_routes)]
         )
 
         print(f"        Solution indices: {np.unique(solution_matrix)}")
@@ -190,6 +186,7 @@ def precalculated_fleet_data(sample_optimization_data):
             no_service_threshold=no_service_threshold,
             allowed_headways=allowed_headways,
             no_service_index=no_service_index,
+            n_directions=sample_optimization_data.get("routes", {}).get("n_directions", None),
         )
 
         # Store calculated results
@@ -212,6 +209,7 @@ def precalculated_fleet_data(sample_optimization_data):
     print("\n✅ PRECALCULATED FLEET DATA READY")
     return result
 
+
 @pytest.fixture
 def usa_population_path():
     """Path to real USA WorldPop data."""

tests/test_constraints.py

@@ -831,7 +831,7 @@ def test_evaluate_pt_only_mode(self):
             "no_service_index": 3,
             "routes": {
                 "round_trip_times": np.array([60.0, 120.0]),  # Floats for safety
-                "route_ids": ["r1", "r2"],
+                "route_ids": ["r1", "r2"], "n_directions": np.array([1.0, 1.0]),
             },
             "constraints": {
                 "fleet_analysis": {
@@ -878,7 +878,7 @@ def test_evaluate_pt_drt_mode(self):
             "n_intervals": 3,
             "allowed_headways": [10, 20, 30],
             "no_service_index": 3,
-            "routes": {"round_trip_times": np.array([60.0, 120.0]), "route_ids": ["r1", "r2"]},
+            "routes": {"round_trip_times": np.array([60.0, 120.0]), "route_ids": ["r1", "r2"], "n_directions": np.array([1.0, 1.0])},
             "constraints": {
                 "fleet_analysis": {
                     "operational_buffer": 1.0,