DAHS / src /features.py
Vittal-M's picture
Upload 39 files
f3fc7bb verified
"""
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