""" features.py — Feature Extraction for Hybrid Warehouse Scheduler Implements a stateful FeatureExtractor that computes 39 features split into: - 32 scenario-level features describing system-wide state (including 4 disruption-aware + 10 composition-adaptive novel features) - 7 job-level features for per-job priority prediction NEW in DAHS_2: - get_feature_ranges() method: returns {feature_name: (min, max)} from training data """ from __future__ import annotations import json import logging import math from collections import deque from pathlib import Path from typing import Any, Dict, List, Optional, Tuple, Union import numpy as np logger = logging.getLogger(__name__) # ------------------------------------------------------------------------- # Feature name lists (used for DataFrame column labeling) # ------------------------------------------------------------------------- SCENARIO_FEATURE_NAMES: List[str] = [ "n_orders_in_system", "n_express_orders_pct", "avg_due_date_tightness", "fraction_already_late", "zone_utilization_avg", "zone_utilization_std", "bottleneck_zone", "avg_remaining_proc_time", "std_remaining_proc_time", "throughput_last_30min", "breakdown_flag", "n_broken_stations", "lunch_break_flag", "surge_multiplier", "batch_pending_flag", "avg_priority_weight", "max_tardiness_so_far", "sla_breach_rate_current", # Disruption-aware features (novel contribution) "disruption_intensity", "queue_imbalance", "job_mix_entropy", "time_pressure_ratio", # Composition-adaptive features (novel contribution, DAHS 2.1) "pct_type_A", "pct_type_B", "pct_type_C", "pct_type_D", "pct_type_E", "count_type_A", "count_type_B", "count_type_C", "count_type_D", "count_type_E", ] JOB_FEATURE_NAMES: List[str] = [ "job_type_encoded", "proc_time_next_station", "remaining_proc_time", "time_to_due", "time_in_system", "critical_ratio", "station_queue_at_next", ] FEATURE_DESCRIPTIONS = { "n_orders_in_system": "Total jobs currently in the system (waiting + processing)", "n_express_orders_pct": "Fraction of waiting jobs that are express (type E)", "avg_due_date_tightness": "Average (due_date - now) for waiting jobs", "fraction_already_late": "Fraction of waiting jobs past their due date", "zone_utilization_avg": "Average utilization across all 8 zones", "zone_utilization_std": "Std deviation of zone utilization (imbalance indicator)", "bottleneck_zone": "Utilization of the most-loaded zone", "avg_remaining_proc_time": "Average remaining processing time for waiting jobs", "std_remaining_proc_time": "Std deviation of remaining processing times", "throughput_last_30min": "Jobs completed per minute in the last 30 minutes", "breakdown_flag": "1 if any station is currently broken, else 0", "n_broken_stations": "Number of stations currently under repair", "lunch_break_flag": "1 if shift is currently in lunch break (t=300-360), else 0", "surge_multiplier": "Current time-of-day arrival rate multiplier", "batch_pending_flag": "1 if a truck batch arrival is imminent", "avg_priority_weight": "Average priority weight of waiting jobs", "max_tardiness_so_far": "Maximum job tardiness observed so far", "sla_breach_rate_current": "Fraction of completed jobs that breached SLA", "disruption_intensity": "[NOVEL] Composite disruption score: breakdowns + lunch + surge", "queue_imbalance": "[NOVEL] Coefficient of variation of queue sizes across zones", "job_mix_entropy": "[NOVEL] Shannon entropy of job-type distribution in queue", "time_pressure_ratio": "[NOVEL] Fraction of waiting jobs with Critical Ratio < 1", "pct_type_A": "[NOVEL] Fraction of waiting jobs of type A (standard)", "pct_type_B": "[NOVEL] Fraction of waiting jobs of type B (picking-intensive)", "pct_type_C": "[NOVEL] Fraction of waiting jobs of type C (value-add)", "pct_type_D": "[NOVEL] Fraction of waiting jobs of type D (complex/bulk)", "pct_type_E": "[NOVEL] Fraction of waiting jobs of type E (express)", "count_type_A": "[NOVEL] Absolute count of waiting type-A jobs", "count_type_B": "[NOVEL] Absolute count of waiting type-B jobs", "count_type_C": "[NOVEL] Absolute count of waiting type-C jobs", "count_type_D": "[NOVEL] Absolute count of waiting type-D jobs", "count_type_E": "[NOVEL] Absolute count of waiting type-E jobs", } # Job type → integer encoding _JOB_TYPE_ENC: Dict[str, int] = {"A": 0, "B": 1, "C": 2, "D": 3, "E": 4} # Job type → priority weight (mirrors simulator definitions) _JOB_PRIORITY_WEIGHT: Dict[str, float] = { "A": 2.0, "B": 1.5, "C": 1.0, "D": 0.8, "E": 3.0 } class FeatureExtractor: """Stateful extractor that maintains running statistics across events. Call ``update(event_type, data)`` as events occur during simulation, then call ``extract_scenario_features`` or ``extract_job_features`` to obtain the feature vectors. NEW in DAHS_2: - get_feature_ranges(): returns {feature_name: (min, max)} from a training DataFrame """ # Window size for throughput tracking (minutes) THROUGHPUT_WINDOW = 30.0 def __init__(self) -> None: # Circular buffer of (timestamp, job_id) for throughput window self._completion_times: deque = deque() # Batch pending flag set externally when a truck batch is imminent self.batch_pending: bool = False # Stored feature ranges for OOD detection (set after training) self._feature_ranges: Optional[Dict[str, Tuple[float, float]]] = None # Metadata loaded alongside the ranges (run hash etc.) — used by the # selector loader to detect stale artifacts. self._feature_ranges_meta: Dict[str, Any] = {} # ------------------------------------------------------------------ # Event update # ------------------------------------------------------------------ def update(self, event_type: str, data: Dict[str, Any]) -> None: """Update running statistics on job events.""" if event_type == "job_complete": self._completion_times.append(data.get("timestamp", 0.0)) # ------------------------------------------------------------------ # Scenario-level features (22) # ------------------------------------------------------------------ def extract_scenario_features(self, sim_state: Dict[str, Any]) -> np.ndarray: """Extract 32 scenario-level features from a system state snapshot. 22 system-state features (F1-F22, including 4 disruption-aware novel) + 10 composition-adaptive features (F23-F32, novel in DAHS 2.1). Parameters ---------- sim_state : dict Output of ``WarehouseSimulator.get_state_snapshot()``. Returns ------- np.ndarray of shape (32,) """ now: float = sim_state.get("current_time", 0.0) waiting_jobs: List[Any] = sim_state.get("waiting_jobs", []) completed_jobs: List[Any] = sim_state.get("completed_jobs", []) queue_sizes: Dict[int, int] = sim_state.get("queue_sizes", {}) zone_util: Dict[int, float] = sim_state.get("zone_utilization", {}) n_broken: int = sim_state.get("n_broken_stations", 0) lunch: bool = sim_state.get("lunch_active", False) surge: float = sim_state.get("surge_multiplier", 1.0) # F1: n_orders_in_system n_in_system = float(sim_state.get("n_orders_in_system", 0)) # F2: n_express_orders_pct n_express = sum(1 for j in waiting_jobs if j.job_type == "E") n_express_pct = n_express / max(1.0, n_in_system) # F3: avg_due_date_tightness = avg(due_date - now) for waiting jobs if waiting_jobs: tightness = float(np.mean([j.due_date - now for j in waiting_jobs])) else: tightness = 999.0 # F4: fraction_already_late if waiting_jobs: frac_late = sum(1 for j in waiting_jobs if j.due_date < now) / len(waiting_jobs) else: frac_late = 0.0 # F5/F6: zone utilization avg and std util_vals = list(zone_util.values()) util_avg = float(np.mean(util_vals)) if util_vals else 0.0 util_std = float(np.std(util_vals)) if util_vals else 0.0 # F7: bottleneck_zone (utilization value of the most-loaded zone) # Bug fix from DAHS_1: use max(zone_util.values()) NOT zone_id if zone_util: bottleneck = float(max(zone_util.values())) else: bottleneck = 0.0 # F8/F9: avg and std remaining proc time for waiting jobs rem_times = [j.remaining_proc_time() for j in waiting_jobs] avg_rem = float(np.mean(rem_times)) if rem_times else 0.0 std_rem = float(np.std(rem_times)) if rem_times else 0.0 # F10: throughput in last 30 min (completions per minute) cutoff = now - self.THROUGHPUT_WINDOW while self._completion_times and self._completion_times[0] < cutoff: self._completion_times.popleft() throughput_30 = len(self._completion_times) / self.THROUGHPUT_WINDOW # F11: breakdown_flag breakdown_flag = 1.0 if n_broken > 0 else 0.0 # F12: n_broken_stations n_broken_f = float(n_broken) # F13: lunch_break_flag lunch_flag = 1.0 if lunch else 0.0 # F14: surge_multiplier surge_f = float(surge) # F15: batch_pending_flag batch_flag = 1.0 if self.batch_pending else 0.0 # F16: avg_priority_weight if waiting_jobs: avg_prio_w = float(np.mean([ _JOB_PRIORITY_WEIGHT.get(j.job_type, 1.0) for j in waiting_jobs ])) else: avg_prio_w = 1.0 # F17: max_tardiness_so_far if completed_jobs: max_tard = float(max( max(0.0, j.completion_time - j.due_date) for j in completed_jobs )) else: max_tard = 0.0 # F18: sla_breach_rate_current if completed_jobs: breach_rate = sum( 1 for j in completed_jobs if j.completion_time > j.due_date ) / len(completed_jobs) else: breach_rate = 0.0 # F19: disruption_intensity — composite disruption score [0, 1] breakdown_severity = min(1.0, n_broken / 5.0) lunch_severity = 1.0 if lunch else 0.0 surge_deviation = abs(surge - 1.0) disruption_intensity = 0.5 * breakdown_severity + 0.25 * lunch_severity + 0.25 * surge_deviation # F20: queue_imbalance — coefficient of variation of queue sizes # Bug fix: guard with mean > 1e-6 (not > 0) q_vals = list(queue_sizes.values()) if q_vals and np.mean(q_vals) > 1e-6: queue_imbalance = float(min(np.std(q_vals) / np.mean(q_vals), 10.0)) else: queue_imbalance = 0.0 # F21: job_mix_entropy — Shannon entropy of job type distribution in queue if waiting_jobs: type_counts: Dict[str, int] = {} for j in waiting_jobs: type_counts[j.job_type] = type_counts.get(j.job_type, 0) + 1 total_w = len(waiting_jobs) job_mix_entropy = 0.0 for cnt in type_counts.values(): p = cnt / total_w if p > 0: job_mix_entropy -= p * math.log2(p) else: job_mix_entropy = 0.0 # F22: time_pressure_ratio — fraction of waiting jobs with CR < 1 if waiting_jobs: n_under_pressure = 0 for j in waiting_jobs: rem = j.remaining_proc_time() ttd = j.due_date - now cr = ttd / max(rem, 0.001) if rem > 0 else 999.0 if cr < 1.0: n_under_pressure += 1 time_pressure_ratio = n_under_pressure / len(waiting_jobs) else: time_pressure_ratio = 0.0 # F23-F32: composition-adaptive features (per-type % and absolute counts) # These give the selector explicit, non-lossy signal about the current # batch composition — crucial for heuristic adaptation. type_counts: Dict[str, int] = {"A": 0, "B": 0, "C": 0, "D": 0, "E": 0} for j in waiting_jobs: if j.job_type in type_counts: type_counts[j.job_type] += 1 total_w = max(len(waiting_jobs), 1) pct_A = type_counts["A"] / total_w if waiting_jobs else 0.0 pct_B = type_counts["B"] / total_w if waiting_jobs else 0.0 pct_C = type_counts["C"] / total_w if waiting_jobs else 0.0 pct_D = type_counts["D"] / total_w if waiting_jobs else 0.0 pct_E = type_counts["E"] / total_w if waiting_jobs else 0.0 features = np.array([ n_in_system, # F1 n_express_pct, # F2 tightness, # F3 frac_late, # F4 util_avg, # F5 util_std, # F6 bottleneck, # F7 avg_rem, # F8 std_rem, # F9 throughput_30, # F10 breakdown_flag, # F11 n_broken_f, # F12 lunch_flag, # F13 surge_f, # F14 batch_flag, # F15 avg_prio_w, # F16 max_tard, # F17 breach_rate, # F18 disruption_intensity, # F19 (novel) queue_imbalance, # F20 (novel) job_mix_entropy, # F21 (novel) time_pressure_ratio, # F22 (novel) pct_A, # F23 (novel) pct_B, # F24 (novel) pct_C, # F25 (novel) pct_D, # F26 (novel) pct_E, # F27 (novel) float(type_counts["A"]),# F28 (novel) float(type_counts["B"]),# F29 (novel) float(type_counts["C"]),# F30 (novel) float(type_counts["D"]),# F31 (novel) float(type_counts["E"]),# F32 (novel) ], dtype=np.float64) # Sanitize: replace NaN/inf with safe values (training pipeline bug fix) features = np.nan_to_num(features, nan=0.0, posinf=999.0, neginf=-999.0) return features.astype(np.float32) # ------------------------------------------------------------------ # Job-level features (7) # ------------------------------------------------------------------ def extract_job_features(self, job: Any, sim_state: Dict[str, Any]) -> np.ndarray: """Extract 7 job-level features for priority prediction.""" now: float = sim_state.get("current_time", 0.0) queue_sizes: Dict[int, int] = sim_state.get("queue_sizes", {}) jt_enc = float(_JOB_TYPE_ENC.get(job.job_type, 0)) if not job.is_complete: next_op = job.operations[job.current_op_idx] proc_next = float(next_op.nominal_proc_time) else: proc_next = 0.0 rem_proc = float(job.remaining_proc_time()) time_to_due = float(job.due_date - now) time_in_sys = float(now - job.arrival_time) if rem_proc > 0: cr = time_to_due / rem_proc else: cr = 999.0 # large finite value, safe for ML models if not job.is_complete: next_zone = job.operations[job.current_op_idx].zone_id queue_at_next = float(queue_sizes.get(next_zone, 0)) else: queue_at_next = 0.0 features = np.array([ jt_enc, proc_next, rem_proc, time_to_due, time_in_sys, cr, queue_at_next, ], dtype=np.float32) return features # ------------------------------------------------------------------ # Feature names # ------------------------------------------------------------------ def get_feature_names(self, level: str = "scenario") -> List[str]: """Return the ordered list of feature names.""" if level == "scenario": return SCENARIO_FEATURE_NAMES elif level == "job": return JOB_FEATURE_NAMES elif level == "all": return SCENARIO_FEATURE_NAMES + JOB_FEATURE_NAMES else: raise ValueError(f"Unknown level: {level!r}. Use 'scenario', 'job', or 'all'.") # ------------------------------------------------------------------ # NEW in DAHS_2: Feature ranges for OOD detection # ------------------------------------------------------------------ def get_feature_ranges( self, X_train: Optional[np.ndarray] = None, feature_names: Optional[List[str]] = None, ) -> Dict[str, Tuple[float, float]]: """Compute {feature_name: (min, max)} from training data. If X_train is None, returns stored ranges (set by set_feature_ranges()). Parameters ---------- X_train : np.ndarray of shape (n_samples, 22) Training feature matrix. If None, returns cached ranges. feature_names : list of str, optional Column names. Defaults to SCENARIO_FEATURE_NAMES. Returns ------- dict mapping feature_name -> (min_val, max_val) """ if X_train is None: if self._feature_ranges is None: raise ValueError("No training data provided and no cached feature ranges.") return self._feature_ranges names = feature_names or SCENARIO_FEATURE_NAMES ranges = {} for i, name in enumerate(names): if i < X_train.shape[1]: ranges[name] = (float(X_train[:, i].min()), float(X_train[:, i].max())) self._feature_ranges = ranges return ranges def set_feature_ranges(self, ranges: Dict[str, Tuple[float, float]]) -> None: """Set feature ranges for OOD detection (loaded from JSON artifact).""" self._feature_ranges = ranges def load_feature_ranges(self, json_path: "Union[Path, str]") -> Dict[str, Tuple[float, float]]: """Load feature ranges from a JSON file saved by train_selector.py. Accepts both the legacy flat format ({feature_name: [min, max]}) and the wrapped format ({"_meta": {...}, "ranges": {feature_name: [...]}}). Stores any meta payload on `self._feature_ranges_meta` so callers can verify the artifact was produced in the same training run as the model. """ with open(json_path, "r") as f: data = json.load(f) if isinstance(data, dict) and "ranges" in data: self._feature_ranges_meta = data.get("_meta", {}) raw = data["ranges"] else: self._feature_ranges_meta = {} raw = data ranges = {k: (v[0], v[1]) for k, v in raw.items()} self._feature_ranges = ranges return ranges def is_out_of_distribution( self, features: np.ndarray, tolerance: float = 0.10, ) -> bool: """Check if any feature falls outside training range ±10%. Parameters ---------- features : np.ndarray of shape (22,) Scenario features to check. tolerance : float Fractional tolerance beyond training range (default 10%). Returns ------- bool: True if OOD """ if self._feature_ranges is None: return False # no ranges loaded → assume in-distribution for i, name in enumerate(SCENARIO_FEATURE_NAMES): if name not in self._feature_ranges: continue lo, hi = self._feature_ranges[name] val = float(features[i]) span = max(hi - lo, 1e-6) if val < lo - tolerance * span or val > hi + tolerance * span: return True return False