update

Author: PhilippSchmelter

src/main.py

@@ -9,15 +9,15 @@
 from src.markov.transitions import build_transition_matrices
 from src.preprocessing.loader import load_timeseries
 
-SIM_DAYS = 3
+SIM_DAYS = 10
+PER_DAY = 96
 
 
 def _detect_value_col(df: pd.DataFrame) -> str:
-    candidate_cols = ["x", "value", "load", "power", "p_norm", "load_norm"]
-    for c in candidate_cols:
+    for c in ["x", "value", "load", "power", "p_norm", "load_norm"]:
         if c in df.columns and np.issubdtype(df[c].dtype, np.number):
             return c
-    raise KeyError("No numeric load column found – please inspect the dataframe.")
+    raise KeyError("numeric load column missing")
 
 
 def _simulate_series(
@@ -28,54 +28,25 @@ def _simulate_series(
     periods: int,
     rng: np.random.Generator | None = None,
 ) -> pd.DataFrame:
-    """Generate synthetic 15‑min series (timestamp, state, x)."""
     rng = np.random.default_rng() if rng is None else rng
-    timestamps = pd.date_range(start_ts, periods=periods, freq="15min")
+    ts = pd.date_range(start_ts, periods=periods, freq="15min")
     states = np.empty(periods, dtype=int)
     xs = np.empty(periods, dtype=float)
 
     s = start_state
-    for i, ts in enumerate(timestamps):
-        b = bucket_id(ts)
+    for i, t in enumerate(ts):
+        b = bucket_id(t)
         s = rng.choice(probs.shape[1], p=probs[b, s])
         states[i] = s
         xs[i] = sample_value(gmms, b, s, rng=rng)
-    return pd.DataFrame({"timestamp": timestamps, "state": states, "x_sim": xs})
 
-
-def main() -> None:
-    df = load_timeseries(normalize=True, discretize=True)
-    if "bucket" not in df.columns:
-        df = assign_buckets(df)
-
-    value_col = _detect_value_col(df)
-    print("Using load column:", value_col)
-
-    counts = build_transition_counts(df)
-    probs = build_transition_matrices(df)
-
-    _plot_first_25_buckets(counts, probs)
-
-    print("Fitting GMMs … (this may take a moment)")
-    gmms = fit_gmms(df, value_col=value_col)
-
-    periods = SIM_DAYS * 96
-    sim_df = _simulate_series(
-        probs,
-        gmms,
-        start_ts=df["timestamp"].min().normalize(),
-        start_state=int(df["state"].iloc[0]),
-        periods=periods,
-    )
-
-    _plot_simulation_diagnostics(df, sim_df, value_col)
+    return pd.DataFrame({"timestamp": ts, "state": states, "x_sim": xs})
 
 
 def _plot_first_25_buckets(counts: np.ndarray, probs: np.ndarray) -> None:
-    """Heat‑map grid for buckets 0‑24."""
-    buckets = list(range(25))
+    buckets = range(25)
     fig, axes = plt.subplots(5, 5, figsize=(15, 15), sharex=True, sharey=True)
-    vmax = probs[buckets].max()
+    vmax = probs[list(buckets)].max()
     norm = Normalize(vmin=0, vmax=vmax)
 
     for idx, b in enumerate(buckets):
@@ -92,15 +63,15 @@ def _plot_first_25_buckets(counts: np.ndarray, probs: np.ndarray) -> None:
         ax.axis("off")
 
     fig.colorbar(im, ax=axes.ravel().tolist(), shrink=0.6, label="p")
-    fig.suptitle("Transition probabilities – buckets 0‑24", fontsize=14)
-    fig.tight_layout(rect=[0, 0, 0.97, 0.96])
+    fig.suptitle("Transition probabilities – buckets 0–24", fontsize=14)
+    plt.tight_layout(rect=[0, 0, 0.97, 0.96])
     plt.show()
 
 
 def _plot_simulation_diagnostics(
     df: pd.DataFrame, sim: pd.DataFrame, value_col: str
 ) -> None:
-    first_day = sim.iloc[:96]
+    first_day = sim.iloc[:PER_DAY]
     plt.figure(figsize=(10, 3))
     plt.plot(first_day["timestamp"], first_day["x_sim"], marker=".")
     plt.title("Simulated power – first day")
@@ -121,13 +92,43 @@ def _plot_simulation_diagnostics(
     sim["hour"] = sim["timestamp"].dt.hour
     plt.figure(figsize=(10, 4))
     sim.boxplot(column="x_sim", by="hour", grid=False)
-    plt.suptitle("")
     plt.title("Simulated power by hour of day")
     plt.xlabel("hour of day")
     plt.ylabel("normalised load x")
     plt.tight_layout()
     plt.show()
 
 
+def main() -> None:
+    df = load_timeseries(normalize=True, discretize=True)
+    if "bucket" not in df.columns:
+        df = assign_buckets(df)
+
+    val_col = _detect_value_col(df)
+
+    counts = build_transition_counts(df)
+    probs = build_transition_matrices(df)
+
+    _plot_first_25_buckets(counts, probs)
+
+    gmms = fit_gmms(
+        df,
+        value_col=val_col,
+        verbose=1,
+        heartbeat_seconds=60,
+    )
+
+    periods = SIM_DAYS * PER_DAY
+    sim = _simulate_series(
+        probs,
+        gmms,
+        start_ts=df["timestamp"].min().normalize(),
+        start_state=int(df["state"].iloc[0]),
+        periods=periods,
+    )
+
+    _plot_simulation_diagnostics(df, sim, val_col)
+
+
 if __name__ == "__main__":
     main()

src/markov/gmm.py

@@ -1,29 +1,25 @@
 import math
-from typing import Dict, List, Tuple
+from typing import List, Optional, Tuple
 
 import joblib
 import numpy as np
 import pandas as pd
 from sklearn.mixture import GaussianMixture
 
 try:
-    # rich progress bar for terminals / notebooks; optional dependency
-    from tqdm.auto import tqdm  # type: ignore
+    from tqdm.auto import tqdm
 except ImportError:
-    tqdm = None  # type: ignore
+    tqdm = None
 
 from .buckets import NUM_BUCKETS
 from .transition_counts import N_STATES
 
-__all__ = [
-    "GaussianBucketModels",
-    "fit_gmms",
-    "sample_value",
-]
-
+__all__ = ["GaussianBucketModels", "fit_gmms", "sample_value"]
 
 GmmTuple = Tuple[np.ndarray, np.ndarray, np.ndarray]
-GaussianBucketModels = List[List[GmmTuple]]
+GaussianBucketModels = List[List[Optional[GmmTuple]]]
+
+_rng = np.random.default_rng()
 
 
 def _fit_single(
@@ -32,31 +28,35 @@ def _fit_single(
     min_samples: int = 30,
     k_candidates: Tuple[int, ...] = (1, 2, 3),
     random_state: int | None = None,
-) -> GmmTuple:
-    """Fit 1‑3 component spherical GMM; fallback to Normal if too few samples."""
+) -> Optional[GmmTuple]:
+    if x.size == 0:
+        return None
     if len(x) < min_samples:
         mean = float(np.mean(x))
-        var = float(np.var(x) + 1e-6)
-        return (np.array([1.0]), np.array([mean]), np.array([var]))
-
+        var = float(max(np.var(x), 1e-5))
+        return (
+            np.array([1.0], dtype=float),
+            np.array([mean], dtype=float),
+            np.array([var], dtype=float),
+        )
     best_bic = math.inf
     best_gmm: GaussianMixture | None = None
     for k in k_candidates:
         gmm = GaussianMixture(
             n_components=k,
             covariance_type="spherical",
-            n_init="auto",
+            n_init=1,
             random_state=random_state,
         ).fit(x.reshape(-1, 1))
         bic = gmm.bic(x.reshape(-1, 1))
         if bic < best_bic:
             best_bic = bic
             best_gmm = gmm
-
-    weights = best_gmm.weights_
-    means = best_gmm.means_.ravel()
-    variances = best_gmm.covariances_.ravel()
-    return (weights, means, variances)
+    return (
+        best_gmm.weights_,
+        best_gmm.means_.ravel(),
+        best_gmm.covariances_.ravel(),
+    )
 
 
 def fit_gmms(
@@ -65,29 +65,37 @@ def fit_gmms(
     value_col: str = "x",
     bucket_col: str = "bucket",
     state_col: str = "state",
-    min_samples: int = 30,
+    min_samples: int = 5,
     k_candidates: Tuple[int, ...] = (1, 2, 3),
     n_jobs: int = -1,
     random_state: int | None = None,
     verbose: int = 0,
+    heartbeat_seconds: int | None = None,
 ) -> GaussianBucketModels:
-    """Return list [bucket][state] -> (weights, means, variances)."""
+    if heartbeat_seconds:
+        import faulthandler
+
+        faulthandler.enable()
+        faulthandler.dump_traceback_later(heartbeat_seconds, repeat=True)
 
-    samples: Dict[Tuple[int, int], List[float]] = {}
-    for _, row in df[[bucket_col, state_col, value_col]].iterrows():
-        samples.setdefault((row[bucket_col], row[state_col]), []).append(row[value_col])
+    grouped = (
+        df[[bucket_col, state_col, value_col]]
+        .groupby([bucket_col, state_col])[value_col]
+        .apply(list)
+        .to_dict()
+    )
 
     tasks = [
-        ((b, s), samples.get((b, s), []))
+        ((b, s), grouped.get((b, s), []))
         for b in range(NUM_BUCKETS)
         for s in range(N_STATES)
     ]
 
-    iterable = tasks
-    if verbose > 0 and tqdm is not None:
-        iterable = tqdm(tasks, desc="Fitting GMMs", unit="model")
+    iterable = (
+        tqdm(tasks, desc="Fitting GMMs", unit="model") if verbose and tqdm else tasks
+    )
 
-    def _worker(item: Tuple[Tuple[int, int], List[float]]):
+    def _worker(item):
         (b, s), x_list = item
         x = np.asarray(x_list, dtype=float)
         return (
@@ -113,17 +121,16 @@ def _worker(item: Tuple[Tuple[int, int], List[float]]):
     return gmms
 
 
-_rng = np.random.default_rng()
-
-
 def sample_value(
     gmms: GaussianBucketModels,
     bucket: int,
     state: int,
     rng: np.random.Generator | None = None,
 ) -> float:
-    """Draw a normalised load value from the GMM for (bucket, state)."""
-    weights, means, vars_ = gmms[bucket][state]
+    gmm = gmms[bucket][state]
+    if gmm is None:
+        raise ValueError(f"No GMM trained for bucket {bucket}, state {state}")
+    weights, means, vars_ = gmm
     rng = _rng if rng is None else rng
     comp = rng.choice(len(weights), p=weights)
     return float(rng.normal(means[comp], math.sqrt(vars_[comp])))

src/markov/transition_counts.py

@@ -15,11 +15,6 @@ def build_transition_counts(
     bucket_col: str = "bucket",
     dtype=np.uint32,
 ) -> np.ndarray:
-    """
-    Absolute transition counts:
-        C[b, i, j] = # of times state_t=i  → state_{t+1}=j in bucket b
-    Shape = (2 304, 10, 10).
-    """
     df = df.sort_values("timestamp")
 
     s_t = df[state_col].to_numpy(dtype=int)[:-1]