current gmm state

PhilippSchmelter · PhilippSchmelter · commit f47870fdf88a · 2025-06-28T02:07:29.000+02:00
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - `CODEOWNERS` file and dependabot automation [#7](https://github.com/ie3-institute/simonaMarkovLoad/issues/7)
 - Added full Markov-model pipeline [#10](https://github.com/ie3-institute/simonaMarkovLoad/issues/10)
+- Added GMM feature to project [#12](https://github.com/ie3-institute/simonaMarkovLoad/issues/12)
 
 ### Changed
 - Compute instantaneous kW from cumulative kWh via 15-minute differencing [#1](https://github.com/ie3-institute/simonaMarkovLoad/issues/1)
diff --git a/src/main.py b/src/main.py
@@ -1,31 +1,133 @@
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
+from matplotlib.colors import Normalize
 
+from src.markov.buckets import assign_buckets, bucket_id
+from src.markov.gmm import fit_gmms, sample_value
 from src.markov.transition_counts import build_transition_counts
 from src.markov.transitions import build_transition_matrices
 from src.preprocessing.loader import load_timeseries
 
-df = load_timeseries(normalize=True, discretize=True)
+SIM_DAYS = 3
 
 
-counts = build_transition_counts(df)
-probs = build_transition_matrices(df, alpha=1.0)
+def _detect_value_col(df: pd.DataFrame) -> str:
+    candidate_cols = ["x", "value", "load", "power", "p_norm", "load_norm"]
+    for c in candidate_cols:
+        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.")
 
-print("counts shape :", counts.shape)
-print("probs  shape :", probs.shape)
 
+def _simulate_series(
+    probs: np.ndarray,
+    gmms,
+    start_ts: pd.Timestamp,
+    start_state: int,
+    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")
+    states = np.empty(periods, dtype=int)
+    xs = np.empty(periods, dtype=float)
 
-active_buckets = np.where(counts.sum(axis=(1, 2)) > 0)[0]
-bucket = int(active_buckets[0]) if active_buckets.size else 0
-print(f"\nUsing bucket {bucket}")
+    s = start_state
+    for i, ts in enumerate(timestamps):
+        b = bucket_id(ts)
+        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})
 
 
-print("row sums :", probs[bucket].sum(axis=1))
+def main() -> None:
+    df = load_timeseries(normalize=True, discretize=True)
+    if "bucket" not in df.columns:
+        df = assign_buckets(df)
 
-plt.imshow(probs[bucket], aspect="auto")
-plt.title(f"Bucket {bucket} – transition probabilities")
-plt.xlabel("state t+1")
-plt.ylabel("state t")
-plt.colorbar()
-plt.tight_layout()
-plt.show()
+    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)
+
+
+def _plot_first_25_buckets(counts: np.ndarray, probs: np.ndarray) -> None:
+    """Heat‑map grid for buckets 0‑24."""
+    buckets = list(range(25))
+    fig, axes = plt.subplots(5, 5, figsize=(15, 15), sharex=True, sharey=True)
+    vmax = probs[buckets].max()
+    norm = Normalize(vmin=0, vmax=vmax)
+
+    for idx, b in enumerate(buckets):
+        ax = axes.flat[idx]
+        if counts[b].sum() == 0:
+            ax.axis("off")
+            continue
+        im = ax.imshow(probs[b], aspect="auto", origin="lower", norm=norm)
+        ax.set_title(f"Bucket {b}", fontsize=8)
+        ax.set_xticks([])
+        ax.set_yticks([])
+
+    for ax in axes.flat[len(buckets) :]:
+        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])
+    plt.show()
+
+
+def _plot_simulation_diagnostics(
+    df: pd.DataFrame, sim: pd.DataFrame, value_col: str
+) -> None:
+    first_day = sim.iloc[:96]
+    plt.figure(figsize=(10, 3))
+    plt.plot(first_day["timestamp"], first_day["x_sim"], marker=".")
+    plt.title("Simulated power – first day")
+    plt.ylabel("normalised load x")
+    plt.tight_layout()
+    plt.show()
+
+    plt.figure(figsize=(6, 4))
+    plt.hist(df[value_col], bins=50, alpha=0.6, density=True, label="original")
+    plt.hist(sim["x_sim"], bins=50, alpha=0.6, density=True, label="simulated")
+    plt.title("Original vs simulated load distribution")
+    plt.xlabel("normalised load x")
+    plt.ylabel("density")
+    plt.legend()
+    plt.tight_layout()
+    plt.show()
+
+    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()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/markov/gmm.py b/src/markov/gmm.py
@@ -0,0 +1,129 @@
+import math
+from typing import Dict, List, 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
+except ImportError:
+    tqdm = None  # type: ignore
+
+from .buckets import NUM_BUCKETS
+from .transition_counts import N_STATES
+
+__all__ = [
+    "GaussianBucketModels",
+    "fit_gmms",
+    "sample_value",
+]
+
+
+GmmTuple = Tuple[np.ndarray, np.ndarray, np.ndarray]
+GaussianBucketModels = List[List[GmmTuple]]
+
+
+def _fit_single(
+    x: np.ndarray,
+    *,
+    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."""
+    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]))
+
+    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",
+            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)
+
+
+def fit_gmms(
+    df: pd.DataFrame,
+    *,
+    value_col: str = "x",
+    bucket_col: str = "bucket",
+    state_col: str = "state",
+    min_samples: int = 30,
+    k_candidates: Tuple[int, ...] = (1, 2, 3),
+    n_jobs: int = -1,
+    random_state: int | None = None,
+    verbose: int = 0,
+) -> GaussianBucketModels:
+    """Return list [bucket][state] -> (weights, means, variances)."""
+
+    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])
+
+    tasks = [
+        ((b, s), samples.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")
+
+    def _worker(item: Tuple[Tuple[int, int], List[float]]):
+        (b, s), x_list = item
+        x = np.asarray(x_list, dtype=float)
+        return (
+            b,
+            s,
+            _fit_single(
+                x,
+                min_samples=min_samples,
+                k_candidates=k_candidates,
+                random_state=random_state,
+            ),
+        )
+
+    results = joblib.Parallel(n_jobs=n_jobs, verbose=verbose)(
+        joblib.delayed(_worker)(task) for task in iterable
+    )
+
+    gmms: GaussianBucketModels = [
+        [None for _ in range(N_STATES)] for _ in range(NUM_BUCKETS)
+    ]
+    for b, s, gmm_tuple in results:
+        gmms[b][s] = gmm_tuple
+    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]
+    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])))
