from __future__ import annotations from dataclasses import dataclass import numpy as np import pandas as pd from sklearn.ensemble import ( HistGradientBoostingClassifier, HistGradientBoostingRegressor, RandomForestClassifier, ) from sklearn.impute import SimpleImputer from sklearn.linear_model import LogisticRegression from sklearn.model_selection import GroupKFold from scipy.stats import rankdata from .competition_features import ( CONTEXT_COLUMNS, MIN_TRAJECTORY_COVERAGE, TRAJECTORY_COLUMNS, TRAJECTORY_REFERENCE_OFFSET, TRAJECTORY_COVERAGE_COLUMN, TRAJECTORY_LAGS, add_context_columns, numeric_feature_columns, ) def _logit(values: np.ndarray) -> np.ndarray: clipped = np.clip(values, 1e-5, 1.0 - 1e-5) return np.log(clipped / (1.0 - clipped)) @dataclass class FrozenLogisticCalibrator: """Minimal binary logistic inference independent of sklearn pickle internals.""" coefficients: np.ndarray intercept: float @classmethod def from_estimator(cls, estimator: LogisticRegression) -> "FrozenLogisticCalibrator": return cls( coefficients=np.asarray(estimator.coef_[0], dtype=float).copy(), intercept=float(estimator.intercept_[0]), ) def predict_proba(self, values: np.ndarray) -> np.ndarray: scores = np.asarray(values, dtype=float) @ self.coefficients + self.intercept positive = np.empty_like(scores, dtype=float) nonnegative = scores >= 0.0 positive[nonnegative] = 1.0 / (1.0 + np.exp(-scores[nonnegative])) exp_scores = np.exp(scores[~nonnegative]) positive[~nonnegative] = exp_scores / (1.0 + exp_scores) return np.column_stack([1.0 - positive, positive]) def _clean(frame: pd.DataFrame, columns: list[str]) -> pd.DataFrame: return frame.reindex(columns=columns).replace([np.inf, -np.inf], np.nan).astype(float) def _derived_features(base: pd.DataFrame, kind: str) -> pd.DataFrame: derived: dict[str, pd.Series] = {} if kind in {"physics", "all"}: for level in ("voltage_mean_3d", "voltage_mean_7d", "voltage_max_3d"): for window in (30, 60, 90, 180): rate = np.maximum(-base[f"voltage_slope_{window}d"], 1e-5) derived[f"knee_{level}_{window}"] = np.clip( (base[level] - 2.40) / rate, -100.0, 1000.0, ) if kind == "all": for short, long in ((3, 14), (7, 30), (14, 60), (30, 90), (30, 180), (60, 365)): derived[f"slope_accel_{short}_{long}"] = ( base[f"voltage_slope_{short}d"] - base[f"voltage_slope_{long}d"] ) derived[f"mean_drop_{short}_{long}"] = ( base[f"voltage_mean_{short}d"] - base[f"voltage_mean_{long}d"] ) for window in (3, 7, 14, 30, 60, 90): derived[f"volt_temp_{window}"] = base[f"voltage_mean_{window}d"] / ( base[f"temperature_mean_{window}d"] + 273.15 ) return pd.concat([base, pd.DataFrame(derived, index=base.index)], axis=1) HAZARD_EDGES = (0, 14, 28, 42, 70, 120, 200, 365) FORECAST_HORIZONS = (30, 60, 90) FORECAST_THRESHOLDS = (2.50, 2.45) TRAJECTORY_BLEND_WEIGHT = 0.25 def _trajectory_features(bins: np.ndarray) -> pd.DataFrame: """Shape of the degradation trajectory: levels, drops, steps and summaries.""" current = bins[:, 0] highest = np.nanmax(bins, axis=1) lowest = np.nanmin(bins, axis=1) spread = np.nanstd(bins, axis=1) columns: dict[str, np.ndarray] = {} for position, lag in enumerate(TRAJECTORY_LAGS): columns[f"v_lag{lag}"] = bins[:, position] for position, lag in enumerate(TRAJECTORY_LAGS[1:], start=1): columns[f"drop{lag}"] = current - bins[:, position] columns[f"ndrop{lag}"] = (current - bins[:, position]) / np.where( spread > 1e-6, spread, np.nan ) steps = np.diff(bins, axis=1) for position in range(steps.shape[1]): columns[f"step{position}"] = -steps[:, position] columns.update( own_max=highest, own_min=lowest, own_std=spread, cur_minus_max=current - highest, cur_over_max=current / highest, pos=(current - lowest) / np.where(highest - lowest > 1e-6, highest - lowest, np.nan), fb50=np.nanmean(bins < 2.50, axis=1), fb45=np.nanmean(bins < 2.45, axis=1), mono=np.nanmean(steps < 0, axis=1), ) return pd.DataFrame(columns) def _snapshot_bins(snapshot: pd.DataFrame) -> np.ndarray | None: if not set(TRAJECTORY_COLUMNS).issubset(snapshot.columns): return None return snapshot.reindex(columns=list(TRAJECTORY_COLUMNS)).to_numpy(dtype=float) def _dense_transitions(matrix: dict, buildings: dict[str, str]): """Every observed degradation transition, not just the terminal failures.""" grid = matrix["grid"] first = matrix["first"] devices = list(matrix["index"]) rows, cuts = [], [] for name, position in matrix["index"].items(): begin = int(first[position]) if begin >= grid.shape[1]: continue end = int(np.flatnonzero(np.isfinite(grid[position]))[-1]) start = begin + 180 for cut in range(start, end + 1, 14): rows.append(position) cuts.append(cut) rows = np.asarray(rows, dtype=np.int64) cuts = np.asarray(cuts, dtype=np.int64) lags = np.asarray(TRAJECTORY_LAGS, dtype=np.int64) + TRAJECTORY_REFERENCE_OFFSET columns = cuts[:, None] - lags[None, :] safe = np.clip(columns, 0, grid.shape[1] - 1) values = grid[rows[:, None], safe] usable = (columns >= 0) & (columns >= first[rows][:, None]) bins = np.where(usable, values, np.nan) names = np.array([devices[position] for position in rows]) where = np.array([buildings.get(name, "") for name in names]) return rows, cuts, bins, names, where def _forecast_targets(matrix: dict, rows, cuts, eol_day: dict[str, float], names): """Threshold crossings and future voltage, read only from observed history.""" grid = matrix["grid"] last = matrix["observed_last"] ends_at = np.array([eol_day.get(name, np.nan) for name in names], dtype=float) targets: dict[tuple, np.ndarray] = {} for horizon in FORECAST_HORIZONS: finish = cuts + horizon reachable = finish <= last[rows] died = np.isfinite(ends_at) & (ends_at <= finish) lowest = np.full(len(rows), np.nan) for position in range(len(rows)): begin = cuts[position] + 1 stop = min(int(finish[position]), int(last[rows[position]])) if stop >= begin: lowest[position] = np.nanmin(grid[rows[position], begin : stop + 1]) for threshold in FORECAST_THRESHOLDS: label = np.where( died, 1.0, np.where(reachable, (lowest <= threshold).astype(float), np.nan) ) targets[(horizon, threshold)] = np.where(~reachable & ~died, np.nan, label) future = np.where(reachable, grid[rows, np.clip(finish, 0, grid.shape[1] - 1)], np.nan) targets[(horizon, "voltage")] = np.where(died, 2.30, future) return targets def _forecaster(random_state: int, regression: bool): common = dict(learning_rate=0.06, max_iter=200, max_leaf_nodes=15, min_samples_leaf=40, l2_regularization=5.0, early_stopping=True, validation_fraction=0.1, n_iter_no_change=15, random_state=random_state) if regression: return HistGradientBoostingRegressor(**common) return HistGradientBoostingClassifier(loss="log_loss", **common) def _fit_forecasters(matrix, buildings, eol_day, random_state, snapshot_features, snapshot_buildings): """Learn degradation dynamics from dense transitions; score snapshots by held-out building.""" rows, cuts, bins, names, where = _dense_transitions(matrix, buildings) if len(rows) < 500: return None, None features = _trajectory_features(bins) targets = _forecast_targets(matrix, rows, cuts, eol_day, names) unique = np.array(sorted(set(where) - {""})) if unique.size < 3: return None, None assignment = {name: position % 3 for position, name in enumerate(unique)} dense_fold = np.array([assignment.get(name, -1) for name in where]) snap_fold = np.array([assignment.get(str(name), -1) for name in snapshot_buildings]) fitted: dict[str, object] = {} holdout: dict[str, np.ndarray] = {} for (horizon, kind), values in targets.items(): regression = kind == "voltage" name = f"v{horizon}" if regression else f"p{horizon}_{kind}" usable = np.isfinite(values) if usable.sum() < 200: continue estimator = _forecaster(random_state, regression) estimator.fit(features[usable], values[usable]) fitted[name] = estimator out = np.full(len(snapshot_features), np.nan) for group in range(3): train = usable & (dense_fold != group) target_rows = np.flatnonzero(snap_fold == group) if train.sum() < 100 or target_rows.size == 0: continue fold_model = _forecaster(random_state + group + 1, regression) fold_model.fit(features[train], values[train]) block = snapshot_features.iloc[target_rows] out[target_rows] = ( fold_model.predict(block) if regression else fold_model.predict_proba(block)[:, 1] ) holdout[name] = out if not fitted: return None, None return fitted, pd.DataFrame(holdout) def _forecast_features(forecasters, features: pd.DataFrame) -> pd.DataFrame: out: dict[str, np.ndarray] = {} for name, estimator in forecasters.items(): if name.startswith("v"): out[name] = estimator.predict(features) else: out[name] = estimator.predict_proba(features)[:, 1] return pd.DataFrame(out, index=features.index) def _context_features(snapshot: pd.DataFrame) -> pd.DataFrame: """Scenario-level context, already attached by scenario_snapshot.""" missing = [name for name in CONTEXT_COLUMNS if name not in snapshot.columns] if missing: snapshot = add_context_columns(snapshot.copy()) return snapshot.reindex(columns=list(CONTEXT_COLUMNS)).astype(float) def _context_estimators(random_state: int) -> dict[str, object]: estimators: dict[str, object] = {} for seed in range(3): estimators[f"hist_{seed}"] = HistGradientBoostingClassifier( loss="log_loss", learning_rate=0.05, max_iter=300, max_leaf_nodes=15, min_samples_leaf=30, l2_regularization=5.0, early_stopping=True, validation_fraction=0.12, n_iter_no_change=25, class_weight={0: 1.0, 1: 5.0}, random_state=random_state + seed, ) estimators["forest"] = RandomForestClassifier( n_estimators=350, min_samples_leaf=5, max_features=0.7, class_weight="balanced_subsample", n_jobs=-1, random_state=random_state, ) return estimators def _fit_context(features: pd.DataFrame, target: np.ndarray, random_state: int): estimators = _context_estimators(random_state) imputer = SimpleImputer(strategy="median") imputed = imputer.fit_transform(features) for name, estimator in estimators.items(): estimator.fit(imputed if name == "forest" else features, target) return estimators, imputer def _context_risk(estimators, imputer, features: pd.DataFrame) -> np.ndarray: imputed = imputer.transform(features) parts = [ estimator.predict_proba(imputed if name == "forest" else features)[:, 1] for name, estimator in estimators.items() ] return np.mean(parts, axis=0) def _fit_hazard(base: pd.DataFrame, target_rul: pd.Series, observed: pd.Series, random_state: int): """Person-period model: one row per battery per interval it is still alive.""" rul = target_rul.to_numpy(dtype=float) events = observed.to_numpy(dtype=bool) rows, bins, labels = [], [], [] for index in range(len(HAZARD_EDGES) - 1): low, high = HAZARD_EDGES[index], HAZARD_EDGES[index + 1] alive = np.flatnonzero(rul >= low) if alive.size == 0: continue rows.append(alive) bins.append(np.full(alive.shape, index)) labels.append((events[alive] & (rul[alive] < high)).astype(int)) if not rows: return None positions = np.concatenate(rows) features = base.iloc[positions].copy() features["hazard_interval"] = np.concatenate(bins) target = np.concatenate(labels) if target.sum() < 20: return None estimator = HistGradientBoostingClassifier( loss="log_loss", learning_rate=0.05, max_iter=300, max_leaf_nodes=15, min_samples_leaf=40, l2_regularization=5.0, early_stopping=True, validation_fraction=0.12, n_iter_no_change=25, random_state=random_state, ) estimator.fit(features, target) return estimator def _hazard_risk(estimator, base: pd.DataFrame, horizon_days: int) -> np.ndarray: survival = np.ones(len(base), dtype=float) for index in range(len(HAZARD_EDGES) - 1): if HAZARD_EDGES[index + 1] > horizon_days: break features = base.copy() features["hazard_interval"] = index survival *= 1.0 - estimator.predict_proba(features)[:, 1] return 1.0 - survival def _scenario_ranks(values: np.ndarray, groups: np.ndarray | None) -> np.ndarray: if groups is None: return rankdata(values) / max(len(values), 1) ranks = np.empty(len(values), dtype=float) for key in pd.unique(groups): mask = groups == key ranks[mask] = rankdata(values[mask]) / max(int(mask.sum()), 1) return ranks def _event_estimators(random_state: int) -> dict[str, object]: common = dict( loss="log_loss", learning_rate=0.05, max_iter=240, early_stopping=True, validation_fraction=0.12, n_iter_no_change=25, random_state=random_state, ) return { "weighted_hist": HistGradientBoostingClassifier( max_leaf_nodes=15, min_samples_leaf=35, l2_regularization=4.0, class_weight={0: 1.0, 1: 5.0}, **common, ), "forest": RandomForestClassifier( n_estimators=350, min_samples_leaf=5, max_features=0.7, class_weight="balanced_subsample", n_jobs=-1, random_state=random_state, ), "all_hist": HistGradientBoostingClassifier( max_leaf_nodes=31, min_samples_leaf=25, l2_regularization=5.0, **common, ), "physics_hist": HistGradientBoostingClassifier( max_leaf_nodes=31, min_samples_leaf=25, l2_regularization=5.0, **common, ), } def _fit_event_estimators( estimators: dict[str, object], base: pd.DataFrame, all_features: pd.DataFrame, physics_features: pd.DataFrame, target: np.ndarray, sample_weight: np.ndarray | None = None, ) -> tuple[dict[str, object], SimpleImputer]: imputer = SimpleImputer(strategy="median") forest_features = imputer.fit_transform(base) estimators["weighted_hist"].fit(base, target, sample_weight=sample_weight) estimators["forest"].fit(forest_features, target, sample_weight=sample_weight) estimators["all_hist"].fit(all_features, target, sample_weight=sample_weight) estimators["physics_hist"].fit(physics_features, target, sample_weight=sample_weight) return estimators, imputer def _device_weights(devices: pd.Series) -> np.ndarray: """One vote per battery, so devices seen in many scenarios do not dominate.""" counts = devices.groupby(devices).transform("size").to_numpy(dtype=float) weights = 1.0 / np.maximum(counts, 1.0) return weights * (len(weights) / weights.sum()) def _event_components( estimators: dict[str, object], imputer: SimpleImputer, base: pd.DataFrame, all_features: pd.DataFrame, physics_features: pd.DataFrame, ) -> np.ndarray: return np.column_stack( [ estimators["weighted_hist"].predict_proba(base)[:, 1], estimators["forest"].predict_proba(imputer.transform(base))[:, 1], estimators["all_hist"].predict_proba(all_features)[:, 1], estimators["physics_hist"].predict_proba(physics_features)[:, 1], ] ) @dataclass class EventTimeModel: event_estimators: dict[str, object] event_imputer: SimpleImputer event_calibrator: FrozenLogisticCalibrator rul_estimator: HistGradientBoostingRegressor survivor_rul_estimator: HistGradientBoostingRegressor | None feature_columns: list[str] all_feature_columns: list[str] physics_feature_columns: list[str] quantile: float horizon_days: int dataset_revision: str balanced: "EventTimeModel | None" = None context_estimators: dict | None = None context_imputer: SimpleImputer | None = None context_columns: list[str] | None = None hazard_estimator: object | None = None blend_weights: tuple[float, float, float, float] = (1.0, 0.0, 0.0, 0.0) forecasters: dict | None = None residual_model: object | None = None residual_columns: list[str] | None = None trajectory_weight: float = 0.0 def _feature_sets( self, snapshot: pd.DataFrame ) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]: base = _clean(snapshot, self.feature_columns) all_features = _clean(_derived_features(base, "all"), self.all_feature_columns) physics_features = _clean( _derived_features(base, "physics"), self.physics_feature_columns ) return base, all_features, physics_features def _calibrated_risk(self, snapshot: pd.DataFrame) -> np.ndarray: base, all_features, physics_features = self._feature_sets(snapshot) components = _event_components( self.event_estimators, self.event_imputer, base, all_features, physics_features, ) stacked = _logit(components) return np.asarray(self.event_calibrator.predict_proba(stacked)[:, 1]) def predict_risk_components(self, snapshot: pd.DataFrame) -> dict[str, np.ndarray]: components = {"primary": self._calibrated_risk(snapshot)} if self.balanced is not None: components["balanced"] = self.balanced._calibrated_risk(snapshot) if self.context_estimators is not None: features = _clean( pd.concat( [_clean(snapshot, self.feature_columns), _context_features(snapshot)], axis=1, ), self.context_columns, ) components["context"] = _context_risk( self.context_estimators, self.context_imputer, features ) if self.hazard_estimator is not None: components["hazard"] = _hazard_risk( self.hazard_estimator, _clean(snapshot, self.feature_columns), self.horizon_days ) return components def _residual_risk(self, snapshot: pd.DataFrame): """Degradation-dynamics view; None when the trajectory history is unusable.""" if self.residual_model is None or self.forecasters is None: return None, None bins = _snapshot_bins(snapshot) if bins is None: return None, None trajectory = _trajectory_features(bins) forecast = _forecast_features(self.forecasters, trajectory) frame = pd.concat([forecast, trajectory], axis=1) frame = frame.reindex(columns=self.residual_columns).astype(float) risk = self.residual_model.predict_proba(frame)[:, 1] if TRAJECTORY_COVERAGE_COLUMN in snapshot.columns: coverage = pd.to_numeric( snapshot[TRAJECTORY_COVERAGE_COLUMN], errors="coerce" ).to_numpy(dtype=float) else: coverage = np.isfinite(bins).mean(axis=1) coverage = np.where(np.isfinite(coverage), coverage, 0.0) # A battery without enough history keeps its blend position rather than moving on noise. weight = np.where(coverage >= MIN_TRAJECTORY_COVERAGE, self.trajectory_weight, 0.0) return risk, weight def predict_event_risk(self, snapshot: pd.DataFrame) -> np.ndarray: groups = ( snapshot["scenario"].astype(str).to_numpy() if "scenario" in snapshot.columns else None ) primary = blend_risk_components( self.predict_risk_components(snapshot), self.blend_weights, groups ) risk, weight = self._residual_risk(snapshot) if risk is None or not np.any(weight > 0.0): return primary return blend_trajectory(primary, risk, weight, groups) def predict_rul(self, snapshot: pd.DataFrame) -> np.ndarray: base = _clean(snapshot, self.feature_columns) return np.asarray(self.rul_estimator.predict(base), dtype=float) def predict_survivor_rul(self, snapshot: pd.DataFrame) -> np.ndarray | None: """Expected remaining life of batteries that outlive the planning window.""" if self.survivor_rul_estimator is None: return None base = _clean(snapshot, self.feature_columns) return np.asarray(self.survivor_rul_estimator.predict(base), dtype=float) def blend_risk_components( components: dict[str, np.ndarray], weights: tuple[float, float, float, float], groups: np.ndarray | None, ) -> np.ndarray: """Rank each view inside a scenario, then map back onto the calibrated scale.""" order = ("primary", "balanced", "context", "hazard") primary = np.asarray(components["primary"], dtype=float) active = [ (np.asarray(components[name], dtype=float), float(weights[index])) for index, name in enumerate(order) if name in components and float(weights[index]) > 0.0 ] if len(active) <= 1: return primary blended = np.zeros(len(primary), dtype=float) total = 0.0 for values, weight in active: blended += weight * _scenario_ranks(values, groups) total += weight blended /= max(total, 1e-9) return _map_to_scale(blended, primary, groups) def blend_trajectory(primary, residual, weight, groups): """Nudge the ranking by the degradation view, keeping the calibrated scale.""" primary = np.asarray(primary, dtype=float) residual = np.asarray(residual, dtype=float) weight = np.asarray(weight, dtype=float) residual = np.where(np.isfinite(residual), residual, np.nanmin(residual[np.isfinite(residual)]) if np.isfinite(residual).any() else 0.0) blended = np.empty(len(primary), dtype=float) keys = [None] if groups is None else list(pd.unique(groups)) for key in keys: mask = slice(None) if key is None else (groups == key) size = len(primary[mask]) first = rankdata(primary[mask]) / max(size, 1) second = rankdata(residual[mask]) / max(size, 1) local = weight[mask] blended[mask] = (1.0 - local) * first + local * second return _map_to_scale(blended, primary, groups) def _map_to_scale(blended: np.ndarray, reference: np.ndarray, groups) -> np.ndarray: """Keep the calibrated marginal distribution and take only the new ordering.""" mapped = np.empty(len(blended), dtype=float) keys = [None] if groups is None else list(pd.unique(groups)) for key in keys: mask = slice(None) if key is None else (groups == key) scores = blended[mask] scale = np.sort(np.asarray(reference)[mask]) order = rankdata(scores, method="ordinal").astype(int) - 1 mapped[mask] = scale[order] return mapped def _fit_trajectory_stack(work, trajectory, eol_times, event_target, random_state): """Dense degradation forecasting plus the temporal-bin view, as one residual ranker.""" bins = _snapshot_bins(work) if bins is None: return None, None, None, 0.0 snapshot_features = _trajectory_features(bins) buildings = dict(zip(work["battery"].astype(str), work["building"].astype(str))) eol_day: dict[str, float] = {} if eol_times is not None: origin = trajectory["origin"] for name, when in eol_times.items(): if pd.notna(when): eol_day[str(name)] = float((pd.Timestamp(when).normalize() - origin).days) forecasters, holdout = _fit_forecasters( trajectory, buildings, eol_day, random_state, snapshot_features, work["building"].astype(str).to_numpy(), ) if forecasters is None: return None, None, None, 0.0 frame = pd.concat([holdout.reset_index(drop=True), snapshot_features.reset_index(drop=True)], axis=1) columns = list(frame.columns) residual = HistGradientBoostingClassifier( loss="log_loss", learning_rate=0.05, max_iter=300, max_leaf_nodes=15, min_samples_leaf=30, l2_regularization=5.0, early_stopping=True, validation_fraction=0.12, n_iter_no_change=25, class_weight={0: 1.0, 1: 5.0}, random_state=random_state, ) residual.fit(frame.astype(float), event_target) return forecasters, residual, columns, TRAJECTORY_BLEND_WEIGHT def _fit_survivor_estimator( base: pd.DataFrame, target_rul: pd.Series, observed: pd.Series, horizon_days: int, maximum_days: float, random_state: int, ) -> HistGradientBoostingRegressor | None: """Median remaining life among batteries observed to outlive the window.""" mask = observed & target_rul.gt(float(horizon_days)) & target_rul.le(float(maximum_days)) if int(mask.sum()) < 50: return None estimator = HistGradientBoostingRegressor( loss="quantile", quantile=0.5, learning_rate=0.05, max_iter=200, max_leaf_nodes=15, min_samples_leaf=40, l2_regularization=4.0, early_stopping=True, validation_fraction=0.15, n_iter_no_change=20, random_state=random_state, ) estimator.fit(base.loc[mask], target_rul.loc[mask]) return estimator def fit_event_time_model( training_snapshots: pd.DataFrame, *, quantile: float = 0.05, horizon_days: int = 42, rul_training_max_days: float = 90.0, survivor_training_max_days: float = 400.0, dataset_revision: str = "unknown", random_state: int = 2026, calibration_folds: int = 5, device_balanced: bool = False, with_auxiliary: bool = True, trajectory: dict | None = None, eol_times: pd.Series | None = None, blend_weights: tuple[float, float, float, float] = (0.55, 0.15, 0.15, 0.15), max_iter: int | None = None, ) -> EventTimeModel: del max_iter # retained for compatibility with earlier validation commands work = training_snapshots.dropna(subset=["target_rul_days"]).copy() columns = numeric_feature_columns(work) base = _clean(work, columns) all_features = _derived_features(base, "all") physics_features = _derived_features(base, "physics") all_columns = list(all_features.columns) physics_columns = list(physics_features.columns) target_rul = pd.to_numeric(work["target_rul_days"], errors="coerce").astype(float) event_target = ( work["event_observed"].astype(bool) & target_rul.between(0.0, float(horizon_days)) ).astype("int8").to_numpy() devices = work["battery"].astype(str) groups = devices.to_numpy() weights = _device_weights(devices) if device_balanced else None folds = min(int(calibration_folds), len(np.unique(groups))) if folds < 2: raise ValueError("At least two device groups are required for event calibration") oof_components = np.full((len(work), 4), np.nan) splitter = GroupKFold(n_splits=folds) for fold_number, (train_index, valid_index) in enumerate( splitter.split(base, groups=groups), start=1 ): fold_estimators = _event_estimators(random_state + fold_number) fold_estimators, fold_imputer = _fit_event_estimators( fold_estimators, base.iloc[train_index], all_features.iloc[train_index], physics_features.iloc[train_index], event_target[train_index], None if weights is None else weights[train_index], ) oof_components[valid_index] = _event_components( fold_estimators, fold_imputer, base.iloc[valid_index], all_features.iloc[valid_index], physics_features.iloc[valid_index], ) fitted_calibrator = LogisticRegression(C=0.03, max_iter=1000) fitted_calibrator.fit(_logit(oof_components), event_target, sample_weight=weights) calibrator = FrozenLogisticCalibrator.from_estimator(fitted_calibrator) estimators = _event_estimators(random_state) estimators, imputer = _fit_event_estimators( estimators, base, all_features, physics_features, event_target, weights, ) observed_near_horizon = ( work["event_observed"].astype(bool) & target_rul.gt(0.0) & target_rul.le(float(rul_training_max_days)) ) rul_estimator = HistGradientBoostingRegressor( loss="quantile", quantile=quantile, learning_rate=0.055, max_iter=260, max_leaf_nodes=19, min_samples_leaf=30, l2_regularization=3.0, early_stopping=True, validation_fraction=0.12, n_iter_no_change=25, random_state=random_state, ) rul_estimator.fit(base.loc[observed_near_horizon], target_rul.loc[observed_near_horizon]) survivor_rul_estimator = _fit_survivor_estimator( base, target_rul, work["event_observed"].astype(bool), horizon_days, survivor_training_max_days, random_state, ) balanced = None forecasters = None residual_model = None residual_columns = None trajectory_weight = 0.0 context_estimators = None context_imputer = None context_columns = None hazard_estimator = None weights_used = (1.0, 0.0, 0.0, 0.0) if with_auxiliary: # Diverse views of the same 82 failures: reweighted, context-aware, and time-to-event. balanced = fit_event_time_model( training_snapshots, quantile=quantile, horizon_days=horizon_days, rul_training_max_days=rul_training_max_days, survivor_training_max_days=survivor_training_max_days, dataset_revision=dataset_revision, random_state=random_state, calibration_folds=calibration_folds, device_balanced=True, with_auxiliary=False, ) context_frame = pd.concat([base, _context_features(work)], axis=1) context_columns = list(context_frame.columns) context_estimators, context_imputer = _fit_context( context_frame, event_target, random_state ) hazard_estimator = _fit_hazard( base, target_rul, work["event_observed"].astype(bool), random_state ) weights_used = tuple(float(value) for value in blend_weights) if trajectory is not None: forecasters, residual_model, residual_columns, trajectory_weight = _fit_trajectory_stack( work, trajectory, eol_times, event_target, random_state ) return EventTimeModel( estimators, imputer, calibrator, rul_estimator, survivor_rul_estimator, columns, all_columns, physics_columns, quantile, horizon_days, dataset_revision, balanced, context_estimators, context_imputer, context_columns, hazard_estimator, weights_used, forecasters, residual_model, residual_columns, trajectory_weight, )