gen_data/genera.py

# /// script
# requires-python = ">=3.11"
# dependencies = ["numpy"]
# ///
"""Genera un CSV con T e ΔT = T − T_amb con rumore realistico (vedi docs/SPEC.md)."""

import argparse
import secrets
import sys
import tomllib
from pathlib import Path

import numpy as np

CARTELLA_OUTPUT = Path("output")
GRIGLIA_MS = 5.0  # passo della griglia fine per il filtro del sensore

CHIAVI = [
    "durata_s", "dt_ms", "jitter_ms",
    "amb_inf_c", "amb_A_c", "amb_tau_s", "amb_sigma_c", "amb_passo_c",
    "sigma_c", "sigma_caldo_c", "x_rif_c", "wander_sigma_c", "wander_tau_s", "passo_T_c",
    "t_impulsi_s", "dT_impulsi_c", "k_raffreddamento_hz", "tau_sensore_s",
    "impulsi_jitter_t_s", "impulsi_jitter_amp_rel", "k_jitter_rel",
    "salti_rate_hz", "salti_media_c", "salti_sigma_c", "salti_tau_s", "salti_x_min_c",
    "salti_max_c", "salti_risc_rate_hz", "salti_risc_media_c", "salti_risc_sigma_c", "salti_risc_tau_s",
]


def errore(msg):
    print(f"Errore: {msg}", file=sys.stderr)
    sys.exit(1)


def leggi_config(percorso):
    try:
        with open(percorso, "rb") as f:
            cfg = tomllib.load(f)
    except OSError as e:
        errore(f"impossibile leggere '{percorso}': {e.strerror}")
    except tomllib.TOMLDecodeError as e:
        errore(f"TOML non valido in '{percorso}': {e}")
    mancanti = [k for k in CHIAVI if k not in cfg]
    if mancanti:
        errore(f"chiavi mancanti in '{percorso}': {', '.join(mancanti)}")
    if len(cfg["t_impulsi_s"]) != len(cfg["dT_impulsi_c"]):
        errore("'t_impulsi_s' e 'dT_impulsi_c' hanno lunghezze diverse")
    return cfg


def eccesso(t_s, cfg):
    """Eccesso termico x(t): un gradino per impulso, poi raffreddamento esponenziale."""
    x = np.zeros_like(t_s)
    for tj, dj in zip(cfg["t_impulsi_s"], cfg["dT_impulsi_c"]):
        dopo = t_s >= tj
        x[dopo] += dj * np.exp(-cfg["k_raffreddamento_hz"] * (t_s[dopo] - tj))
    return x


def perturba(cfg, rng):
    """Copia di cfg con istanti, ampiezze e k degli impulsi variati a caso (ogni corsa è diversa)."""
    n = len(cfg["t_impulsi_s"])
    t = np.sort(np.asarray(cfg["t_impulsi_s"], float) + cfg["impulsi_jitter_t_s"] * rng.standard_normal(n))
    d = np.asarray(cfg["dT_impulsi_c"], float) * (1 + cfg["impulsi_jitter_amp_rel"] * rng.standard_normal(n))
    k = cfg["k_raffreddamento_hz"] * (1 + cfg["k_jitter_rel"] * rng.standard_normal())
    return {**cfg, "t_impulsi_s": t, "dT_impulsi_c": d, "k_raffreddamento_hz": k}


def eventi(g_s, cfg, rng, t_ini, t_fin, rate, media, sigma, tau, x_min):
    """Scarti casuali (Poisson) in [t_ini, t_fin], con rientro esponenziale, se x ≥ x_min."""
    s = np.zeros_like(g_s)
    n = rng.poisson(rate * (t_fin - t_ini))
    ts = np.sort(rng.uniform(t_ini, t_fin, n))
    ampiezze = np.clip(rng.normal(media, sigma, n), -cfg["salti_max_c"], cfg["salti_max_c"])
    for tj, aj in zip(ts, ampiezze):
        if eccesso(np.array([tj]), cfg)[0] < x_min:
            continue
        dopo = g_s >= tj
        s[dopo] += aj * np.exp(-(g_s[dopo] - tj) / tau)
    return s


def salti(g_s, cfg, rng):
    """Disturbi a salto: durante il riscaldamento (tra primo e ultimo gradino) e il raffreddamento."""
    if len(cfg["t_impulsi_s"]) == 0:
        return np.zeros_like(g_s)
    t_primo, t_ultimo = cfg["t_impulsi_s"][0], cfg["t_impulsi_s"][-1]
    risc = eventi(g_s, cfg, rng, t_primo, t_ultimo, cfg["salti_risc_rate_hz"],
                  cfg["salti_risc_media_c"], cfg["salti_risc_sigma_c"], cfg["salti_risc_tau_s"], 0.0)
    raff = eventi(g_s, cfg, rng, t_ultimo, cfg["durata_s"], cfg["salti_rate_hz"],
                  cfg["salti_media_c"], cfg["salti_sigma_c"], cfg["salti_tau_s"], cfg["salti_x_min_c"])
    return risc + raff


def ambiente_vero(t_s, cfg):
    return cfg["amb_inf_c"] + cfg["amb_A_c"] * np.exp(-t_s / cfg["amb_tau_s"])


def quantizza(v, passo):
    return np.round(v / passo) * passo


def genera(cfg, rng):
    cfg = perturba(cfg, rng)
    # 1. tempi
    n_max = int(cfg["durata_s"] * 1000 / (cfg["dt_ms"] - cfg["jitter_ms"])) + 2
    dt = cfg["dt_ms"] + rng.uniform(-cfg["jitter_ms"], cfg["jitter_ms"], n_max)
    t_ms = np.cumsum(dt)
    t_ms = t_ms[t_ms <= cfg["durata_s"] * 1000]
    dt = np.diff(t_ms, prepend=0.0)
    t_s = t_ms / 1000

    # 2. T vera su griglia fine = ambiente vero + eccesso termico + salti
    g_s = np.arange(0.0, cfg["durata_s"] + GRIGLIA_MS / 1000, GRIGLIA_MS / 1000)
    t_vera = ambiente_vero(g_s, cfg) + eccesso(g_s, cfg) + salti(g_s, cfg, rng)

    # 3. risposta del sensore: filtro del 1° ordine, interpolato ai t_i
    a = np.exp(-GRIGLIA_MS / 1000 / cfg["tau_sensore_s"])
    y = np.empty_like(t_vera)
    y[0] = t_vera[0]
    for n in range(1, len(y)):
        y[n] = a * y[n - 1] + (1 - a) * t_vera[n]
    y_i = np.interp(t_s, g_s, y)

    # 4. rumore di T: bianco con σ(x) + deriva Ornstein-Uhlenbeck, poi quantizzazione
    x_i = eccesso(t_s, cfg)
    sigma = cfg["sigma_c"] + (cfg["sigma_caldo_c"] - cfg["sigma_c"]) * np.minimum(1.0, x_i / cfg["x_rif_c"])
    phi = np.exp(-(dt / 1000) / cfg["wander_tau_s"])
    inn = rng.standard_normal(len(t_s))
    w = np.empty_like(t_s)
    w[0] = cfg["wander_sigma_c"] * inn[0]
    for i in range(1, len(w)):
        w[i] = phi[i] * w[i - 1] + cfg["wander_sigma_c"] * np.sqrt(1 - phi[i] ** 2) * inn[i]
    T = quantizza(y_i + w + sigma * rng.standard_normal(len(t_s)), cfg["passo_T_c"])

    # 5. ambiente misurato e ΔT (dopo rumore e quantizzazione di entrambi i canali)
    amb = quantizza(
        ambiente_vero(t_s, cfg) + cfg["amb_sigma_c"] * rng.standard_normal(len(t_s)),
        cfg["amb_passo_c"],
    )
    dT = quantizza(T - amb, cfg["passo_T_c"])
    return np.round(t_ms).astype(int), T, dT


def scrivi(f, t_ms, T, dT):
    f.write("t_ms,T_C,dT_C\n")
    for t, a, b in zip(t_ms, T, dT):
        # "+ 0.0" evita "-0.00"
        f.write(f"{t},{round(a, 2) + 0.0:.2f},{round(b, 2) + 0.0:.2f}\n")


def main():
    ap = argparse.ArgumentParser(description="Genera un CSV con T e ΔT con rumore realistico.")
    ap.add_argument("config", help="file di configurazione TOML")
    ap.add_argument("-o", "--output", help="file CSV di uscita ('-' = stdout; default: output/gen_<seed>.csv)")
    ap.add_argument("--seed", type=int, help="seme del generatore (default: casuale, stampato su stderr)")
    args = ap.parse_args()

    cfg = leggi_config(args.config)
    seed = args.seed if args.seed is not None else secrets.randbelow(2**32)
    if args.seed is None:
        print(f"seed: {seed}", file=sys.stderr)

    t_ms, T, dT = genera(cfg, np.random.default_rng(seed))

    if args.output == "-":
        scrivi(sys.stdout, t_ms, T, dT)
        return
    percorso = Path(args.output) if args.output else CARTELLA_OUTPUT / f"gen_{seed}.csv"
    try:
        percorso.parent.mkdir(parents=True, exist_ok=True)
        with open(percorso, "w", newline="\n") as f:
            scrivi(f, t_ms, T, dT)
    except OSError as e:
        errore(f"impossibile scrivere '{percorso}': {e.strerror}")
    print(f"scritto: {percorso}", file=sys.stderr)


if __name__ == "__main__":
    main()