From 70231607bbb6b5798c55174e327a45b0ec99fc82 Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Thu, 29 Oct 2026 09:51:33 +0100 Subject: [PATCH] Programma fitseg: domini di riscaldamento/raffreddamento, fit e report Co-Authored-By: Claude Sonnet 5.5 --- .gitignore | 6 +++ config.toml | 34 ++++++++++++ pyproject.toml | 19 +++++++ src/fitseg/__init__.py | 0 src/fitseg/cli.py | 112 ++++++++++++++++++++++++++++++++++++++++ src/fitseg/config.py | 35 +++++++++++++ src/fitseg/dati.py | 31 +++++++++++ src/fitseg/domini.py | 92 +++++++++++++++++++++++++++++++++ src/fitseg/fit.py | 112 ++++++++++++++++++++++++++++++++++++++++ src/fitseg/modelli.py | 43 +++++++++++++++ src/fitseg/pipeline.py | 38 ++++++++++++++ src/fitseg/riepilogo.py | 100 +++++++++++++++++++++++++++++++++++ src/fitseg/selezione.py | 25 +++++++++ tests/test_pipeline.py | 40 ++++++++++++++ 14 files changed, 687 insertions(+) create mode 100644 .gitignore create mode 100644 config.toml create mode 100644 pyproject.toml create mode 100644 src/fitseg/__init__.py create mode 100644 src/fitseg/cli.py create mode 100644 src/fitseg/config.py create mode 100644 src/fitseg/dati.py create mode 100644 src/fitseg/domini.py create mode 100644 src/fitseg/fit.py create mode 100644 src/fitseg/modelli.py create mode 100644 src/fitseg/pipeline.py create mode 100644 src/fitseg/riepilogo.py create mode 100644 src/fitseg/selezione.py create mode 100644 tests/test_pipeline.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6f9dff9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +.venv/ +lavoro/ +__pycache__/ +*.egg-info/ +.pytest_cache/ +dati/*.csv diff --git a/config.toml b/config.toml new file mode 100644 index 0000000..dc56f9f --- /dev/null +++ b/config.toml @@ -0,0 +1,34 @@ +# config.toml — UNICA fonte dei valori di fitseg. Suffissi: _s secondi, _c gradi C, _rel frazione. +# Nessun default nel codice: chiave mancante = errore (exit 3). + +[dati] +colonna_tempo = "t_ms" +unita_tempo = "ms" # "ms" | "s" +serie = ["T_C", "dT_C"] # colonne da approssimare +separatore = "," + +[domini] # riconoscimento di riposo / riscaldamento / raffreddamento (nessun LLM) +serie_riferimento = "T_C" # serie su cui si cercano i gradini di riscaldamento +soglia_gradino_c = 20.0 # salto positivo tra due campioni consecutivi che segna un gradino di riscaldamento +soglia_bordo_sigma = 6.0 # il gradino si estende ai campioni adiacenti con salto > soglia × sigma robusta dei salti +distanza_min_salti_campioni = 3 # rilevazioni più vicine di così sono lo stesso gradino +n_min_punti = 6 # campioni minimi per dominio (sotto: "indeterminato") +scarta_dopo_salto_n = 3 # campioni esclusi dal fit all'inizio di un dominio dopo un gradino (transizione del sensore) + +[fit] +modelli_abilitati = ["costante", "lineare", "quadratico", "esp_decad", "esp_doppio", "logistica"] +n_start = 20 # partenze multiple per i modelli non lineari (logistica) +perdita = "soft_l1" # "linear" | "soft_l1" | "huber" (robusta agli outlier) +scala_perdita_sigma = 3.0 # scala della perdita robusta, in unità di sigma rumore +tau_min_rel = 0.05 # limite inferiore di tau, come frazione della durata del dominio +tau_max_rel = 50.0 # limite superiore di tau, come frazione della durata del dominio +n_griglia_tau = 40 # punti della griglia log-spaziata su tau (proiezione delle variabili) +soglia_condizione = 1.0e8 # numero di condizione (jacobiano normalizzato) oltre cui i parametri sono "non identificabili" +correggi_autocorrelazione = true # allarga gli intervalli di confidenza se ACF1 dei residui > soglia_acf +soglia_acf = 0.3 +soglia_delta_aicc = 2.0 # entro questa differenza vince il modello con meno parametri +rmse_vs_rumore_max = 1.3 # RMSE <= valore × sigma rumore => modello adeguato +seed = 0 # seme delle partenze multiple + +[output] +cartella_lavoro = "lavoro" diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..ab75691 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,19 @@ +[project] +name = "fitseg" +version = "0.1.0" +description = "Approssimazione a tratti di serie temporali di riscaldamento/raffreddamento" +requires-python = ">=3.11" +dependencies = ["numpy", "scipy", "matplotlib"] + +[project.optional-dependencies] +dev = ["pytest", "ruff"] + +[project.scripts] +fitseg = "fitseg.cli:main" + +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[tool.hatch.build.targets.wheel] +packages = ["src/fitseg"] diff --git a/src/fitseg/__init__.py b/src/fitseg/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/fitseg/cli.py b/src/fitseg/cli.py new file mode 100644 index 0000000..323a530 --- /dev/null +++ b/src/fitseg/cli.py @@ -0,0 +1,112 @@ +"""CLI: fitseg domini | analizza | modelli | valuta.""" + +import argparse +import json +import shutil +import sys +from pathlib import Path + +import numpy as np + +from . import config, dati, domini as dom, riepilogo +from .modelli import MODELLI +from .pipeline import analizza_domini + +EXIT_USO, EXIT_CONFIG, EXIT_DATI, EXIT_DOMINI = 2, 3, 4, 6 + + +def errore(msg, codice): + print(f"Errore: {msg}", file=sys.stderr) + sys.exit(codice) + + +def _carica(args): + try: + cfg = config.leggi(args.config) + except config.ErroreConfig as e: + errore(e, EXIT_CONFIG) + lavoro = Path(args.lavoro or cfg["output"]["cartella_lavoro"]) + return cfg, lavoro + + +def _domini(args, cfg, lavoro, t, serie): + lavoro.mkdir(parents=True, exist_ok=True) + percorso = Path(args.domini) if getattr(args, "domini", None) else None + if percorso: + try: + return dom.leggi_toml(percorso, t) + except (OSError, ValueError, KeyError) as e: + errore(f"{percorso}: {e}", EXIT_DOMINI) + proposta = dom.proponi(t, serie, cfg) + dom.scrivi_toml(lavoro / "domini.toml", proposta, t) + return proposta + + +def cmd_domini(args): + cfg, lavoro = _carica(args) + try: + t, serie = dati.leggi_csv(args.csv, cfg) + except dati.ErroreDati as e: + errore(e, EXIT_DATI) + proposta = _domini(argparse.Namespace(domini=None), cfg, lavoro, t, serie) + for k, d in enumerate(proposta, 1): + print(f"{k:2d} {d['tipo']:15s} {t[d['i0']]:8.2f} – {t[d['i1'] - 1]:8.2f} s n={d['i1'] - d['i0']}") + print(f"scritto: {lavoro / 'domini.toml'}", file=sys.stderr) + + +def cmd_analizza(args): + cfg, lavoro = _carica(args) + try: + t, serie = dati.leggi_csv(args.csv, cfg) + except dati.ErroreDati as e: + errore(e, EXIT_DATI) + domini = _domini(args, cfg, lavoro, t, serie) + out = lavoro / "out" + out.mkdir(parents=True, exist_ok=True) + shutil.copy(args.csv, lavoro / "dati.csv") + risultati = analizza_domini(t, serie, domini, cfg) + riepilogo.scrivi_json(out, risultati) + riepilogo.scrivi_report(out, risultati, Path(args.csv).name) + riepilogo.grafici(out, t, serie, risultati, domini) + for r in risultati: + for nome, vs in r["serie"].items(): + print(f"dominio {r['k']} {r['tipo']:15s} {nome:5s} -> {vs['scelta'] or 'indeterminato'}") + print(f"scritto: {out}/ (report.md, equazioni.md, modello_finale.json, fit_*.png)", file=sys.stderr) + + +def cmd_modelli(args): + for nome, (_, par, eq, _) in MODELLI.items(): + print(f"{nome:12s} {eq:48s} parametri: {', '.join(par)}") + + +def cmd_valuta(args): + modello = json.loads(Path(args.modello).read_text()) + ts = np.array(args.t, float) + for nome in modello[0]["serie"]: + for t, y in zip(ts, riepilogo.valuta_finale(modello, nome, ts)): + print(f"{nome} t={t:g} s -> {y:.4f}") + + +def main(argv=None): + ap = argparse.ArgumentParser(prog="fitseg", description=__doc__) + ap.add_argument("--config", default="config.toml") + ap.add_argument("--lavoro", help="cartella di lavoro (default: [output] cartella_lavoro)") + sub = ap.add_subparsers(dest="cmd", required=True) + p = sub.add_parser("domini", help="propone i domini (riposo/riscaldamento/raffreddamento) e scrive domini.toml") + p.add_argument("csv") + p.set_defaults(f=cmd_domini) + p = sub.add_parser("analizza", help="domini + fit + scelta + report") + p.add_argument("csv") + p.add_argument("--domini", help="domini.toml modificato a mano, al posto della proposta automatica") + p.set_defaults(f=cmd_analizza) + sub.add_parser("modelli", help="elenca il catalogo").set_defaults(f=cmd_modelli) + p = sub.add_parser("valuta", help="valuta modello_finale.json in uno o più istanti") + p.add_argument("modello") + p.add_argument("--t", type=float, nargs="+", required=True) + p.set_defaults(f=cmd_valuta) + args = ap.parse_args(argv) + args.f(args) + + +if __name__ == "__main__": + main() diff --git a/src/fitseg/config.py b/src/fitseg/config.py new file mode 100644 index 0000000..f91767e --- /dev/null +++ b/src/fitseg/config.py @@ -0,0 +1,35 @@ +"""Lettura di config.toml: nessun default, chiave mancante = errore.""" + +import tomllib + +CHIAVI = { + "dati": ["colonna_tempo", "unita_tempo", "serie", "separatore"], + "domini": ["serie_riferimento", "soglia_gradino_c", "soglia_bordo_sigma", + "distanza_min_salti_campioni", "n_min_punti", "scarta_dopo_salto_n"], + "fit": ["modelli_abilitati", "n_start", "perdita", "scala_perdita_sigma", "tau_min_rel", "tau_max_rel", + "n_griglia_tau", "soglia_condizione", "correggi_autocorrelazione", "soglia_acf", + "soglia_delta_aicc", "rmse_vs_rumore_max", "seed"], + "output": ["cartella_lavoro"], +} + + +class ErroreConfig(Exception): + pass + + +def leggi(percorso): + try: + with open(percorso, "rb") as f: + cfg = tomllib.load(f) + except OSError as e: + raise ErroreConfig(f"impossibile leggere '{percorso}': {e.strerror}") + except tomllib.TOMLDecodeError as e: + raise ErroreConfig(f"TOML non valido in '{percorso}': {e}") + mancanti = [f"[{s}] {k}" for s, chiavi in CHIAVI.items() for k in chiavi if k not in cfg.get(s, {})] + if mancanti: + raise ErroreConfig(f"chiavi mancanti in '{percorso}': {', '.join(mancanti)}") + if cfg["dati"]["unita_tempo"] not in ("ms", "s"): + raise ErroreConfig("[dati] unita_tempo deve essere 'ms' o 's'") + if cfg["domini"]["serie_riferimento"] not in cfg["dati"]["serie"]: + raise ErroreConfig("[domini] serie_riferimento deve essere una delle [dati] serie") + return cfg diff --git a/src/fitseg/dati.py b/src/fitseg/dati.py new file mode 100644 index 0000000..70d828b --- /dev/null +++ b/src/fitseg/dati.py @@ -0,0 +1,31 @@ +"""Lettura del CSV: tempo in secondi (crescente) e serie come array.""" + +import numpy as np + + +class ErroreDati(Exception): + pass + + +def leggi_csv(percorso, cfg): + d = cfg["dati"] + try: + with open(percorso, encoding="utf-8") as f: + intestazione = f.readline().strip().split(d["separatore"]) + righe = np.genfromtxt(f, delimiter=d["separatore"]) + except OSError as e: + raise ErroreDati(f"impossibile leggere '{percorso}': {e.strerror}") + colonne = [d["colonna_tempo"], *d["serie"]] + for c in colonne: + if c not in intestazione: + raise ErroreDati(f"colonna '{c}' assente in '{percorso}' (trovate: {', '.join(intestazione)})") + righe = np.atleast_2d(righe) + if righe.size == 0 or np.isnan(righe).any(): + raise ErroreDati(f"'{percorso}': dati vuoti o non numerici") + t = righe[:, intestazione.index(d["colonna_tempo"])] + if d["unita_tempo"] == "ms": + t = t / 1000.0 + if np.any(np.diff(t) <= 0): + raise ErroreDati("il tempo non è strettamente crescente (ipotesi H1)") + serie = {s: righe[:, intestazione.index(s)] for s in d["serie"]} + return t, serie diff --git a/src/fitseg/domini.py b/src/fitseg/domini.py new file mode 100644 index 0000000..7fcc02a --- /dev/null +++ b/src/fitseg/domini.py @@ -0,0 +1,92 @@ +"""Riconoscimento dei domini: riposo, gradini di riscaldamento, raffreddamento (nessun LLM).""" + +import numpy as np + + +def sigma_robusto(d): + """σ robusto di un campione di differenze (MAD).""" + return 1.4826 * np.median(np.abs(d - np.median(d))) + + +def sigma_rumore(y): + """σ del rumore bianco dalla derivata seconda discreta: 1.4826·MAD(d²y)/√6.""" + if len(y) < 4: + return float("nan") + return sigma_robusto(np.diff(y, 2)) / np.sqrt(6) + + +def trova_gradini(y, cfg): + """Indice del primo campione di transizione di ogni gradino di riscaldamento (salto positivo > soglia).""" + c = cfg["domini"] + d = np.diff(y) + forti = np.flatnonzero(d > c["soglia_gradino_c"]) + if len(forti) == 0: + return [] + soglia_bordo = c["soglia_bordo_sigma"] * sigma_robusto(d) + gruppi, corrente = [], [forti[0]] + for i in forti[1:]: + if i - corrente[-1] < c["distanza_min_salti_campioni"]: + corrente.append(i) + else: + gruppi.append(corrente) + corrente = [i] + gruppi.append(corrente) + inizi = [] + for g in gruppi: + # il gradino parte dal primo campione già in salita: i campioni di transizione aprono il nuovo dominio + # (e sono esclusi dal fit da scarta_dopo_salto_n), così non sporcano la fine del dominio precedente + inizio = g[0] + while inizio > 0 and d[inizio - 1] > soglia_bordo: + inizio -= 1 + inizi.append(int(inizio) + 1) + return inizi + + +def proponi(t, serie, cfg): + """Lista di domini {tipo, i0, i1} (i1 escluso) ordinati nel tempo.""" + c = cfg["domini"] + y = serie[c["serie_riferimento"]] + inizi = trova_gradini(y, cfg) + tagli = [0, *inizi, len(y)] + domini = [] + for k in range(len(tagli) - 1): + if not inizi: + tipo = "riposo" + elif k == 0: + tipo = "riposo" + elif k == len(tagli) - 2: + tipo = "raffreddamento" + else: + tipo = "riscaldamento" + domini.append({"tipo": tipo, "i0": tagli[k], "i1": tagli[k + 1]}) + return domini + + +def scrivi_toml(percorso, domini, t): + """domini.toml: modificabile a mano (tipo, t0_s, t1_s); gli indici si ricalcolano dal tempo.""" + righe = ["# domini.toml — proposta automatica, modificabile: t0_s incluso, t1_s escluso (ultimo dominio: incluso)\n"] + for d in domini: + t1 = t[d["i1"]] if d["i1"] < len(t) else t[-1] + righe.append(f'[[dominio]]\ntipo = "{d["tipo"]}" # riposo | riscaldamento | raffreddamento\n' + f't0_s = {t[d["i0"]]:.4f}\nt1_s = {t1:.4f}\n') + percorso.write_text("\n".join(righe), encoding="utf-8") + + +def leggi_toml(percorso, t): + import tomllib + with open(percorso, "rb") as f: + voci = tomllib.load(f).get("dominio", []) + if not voci: + raise ValueError("domini.toml senza [[dominio]]") + domini = [] + for k, v in enumerate(voci): + if v["tipo"] not in ("riposo", "riscaldamento", "raffreddamento"): + raise ValueError(f"dominio {k + 1}: tipo '{v['tipo']}' non valido") + i0 = int(np.searchsorted(t, v["t0_s"] - 1e-6)) + i1 = len(t) if k == len(voci) - 1 else int(np.searchsorted(t, v["t1_s"] - 1e-6)) + if i1 <= i0: + raise ValueError(f"dominio {k + 1}: intervallo vuoto o in ordine errato") + domini.append({"tipo": v["tipo"], "i0": i0, "i1": i1}) + if any(domini[k]["i1"] != domini[k + 1]["i0"] for k in range(len(domini) - 1)): + raise ValueError("i domini devono essere contigui (t1_s di uno = t0_s del successivo)") + return domini diff --git a/src/fitseg/fit.py b/src/fitseg/fit.py new file mode 100644 index 0000000..4ccbb2f --- /dev/null +++ b/src/fitseg/fit.py @@ -0,0 +1,112 @@ +"""Fit di un modello su un dominio con scipy.optimize.least_squares (VARPRO su tau per gli esponenziali).""" + +import numpy as np +from scipy.optimize import least_squares + +from .modelli import MODELLI + + +def _acf1(r): + r = r - r.mean() + den = np.sum(r * r) + return float(np.sum(r[1:] * r[:-1]) / den) if den > 0 and len(r) > 2 else 0.0 + + +def _varpro(t, y, nome, cfg, L): + """Partenza per gli esponenziali: griglia log-spaziata su tau, c e ampiezze per minimi quadrati lineari.""" + f = cfg["fit"] + taus = np.geomspace(f["tau_min_rel"] * L, f["tau_max_rel"] * L, f["n_griglia_tau"]) + migliore, rss_min = None, np.inf + + def prova(ts): + nonlocal migliore, rss_min + X = np.column_stack([np.ones_like(t), *[np.exp(-t / s) for s in ts]]) + coef = np.linalg.lstsq(X, y, rcond=None)[0] + rss = np.sum((y - X @ coef) ** 2) + if rss < rss_min: + rss_min, migliore = rss, (coef, ts) + + if nome == "esp_decad": + for s in taus: + prova([s]) + c, ts = migliore + return [np.array([c[0], c[1], ts[0]])] + grossa = taus[:: max(1, len(taus) // 12)] + for i, s1 in enumerate(grossa): + for s2 in grossa[i + 1:]: + prova([s1, s2]) + c, ts = migliore + return [np.array([c[0], c[1], ts[0], c[2], ts[1]])] + + +def _partenze(t, y, nome, cfg, L, rng): + if nome in ("esp_decad", "esp_doppio"): + return _varpro(t, y, nome, cfg, L) + if nome in ("costante", "lineare", "quadratico"): + k = len(MODELLI[nome][1]) + return [np.linalg.lstsq(np.vander(t, k, increasing=True), y, rcond=None)[0]] + ampiezza = y[-1] - y[0] if abs(y[-1] - y[0]) > 0 else np.ptp(y) # logistica + return [np.array([y[0], ampiezza, rng.uniform(0.1, 0.9) * L, rng.uniform(0.02, 0.3) * L]) + for _ in range(cfg["fit"]["n_start"])] + + +def _limiti(nome, cfg, L): + f = cfg["fit"] + lo = np.full(len(MODELLI[nome][1]), -np.inf) + hi = np.full_like(lo, np.inf) + for i in MODELLI[nome][3]: + lo[i], hi[i] = f["tau_min_rel"] * L, f["tau_max_rel"] * L + if nome == "logistica": + lo[2], hi[2] = 0.0, L + return lo, hi + + +def adatta(t, y, nome, sigma, cfg, rng): + """Fit di `nome` su (t, y). Ritorna un dizionario con parametri, incertezze e metriche.""" + f = cfg["fit"] + funzione, nomi, equazione, idx_tau = MODELLI[nome] + n, k = len(t), len(nomi) + L = float(t[-1] - t[0]) or 1.0 + lo, hi = _limiti(nome, cfg, L) + scala = f["scala_perdita_sigma"] * sigma if np.isfinite(sigma) and sigma > 0 else 1.0 + res_fn = lambda p: funzione(t, p) - y + migliore = None + for p0 in _partenze(t, y, nome, cfg, L, rng): + p0 = np.clip(p0, lo, hi) + try: + r = least_squares(res_fn, p0, bounds=(lo, hi), loss=f["perdita"], f_scale=scala, x_scale="jac") + except (ValueError, FloatingPointError): + continue + if migliore is None or r.cost < migliore.cost: + migliore = r + if migliore is None: + return {"modello": nome, "convergenza": False, "n_parametri": k} + p, J = migliore.x.copy(), migliore.jac + if nome == "esp_doppio" and p[2] > p[4]: # tau1 < tau2 + perm = [0, 3, 4, 1, 2] + p, J = p[perm], J[:, perm] + r = funzione(t, p) - y + rss = float(np.sum(r**2)) + gdl = max(n - k, 1) + Jn = J / np.maximum(np.linalg.norm(J, axis=0), 1e-300) + cond = float(np.linalg.cond(Jn)) + acf1 = _acf1(r) + corr = np.sqrt((1 + acf1) / (1 - acf1)) if f["correggi_autocorrelazione"] and acf1 > f["soglia_acf"] else 1.0 + try: + cov = (rss / gdl) * np.linalg.pinv(J.T @ J) + errori = np.sqrt(np.maximum(np.diag(cov), 0)) * corr + except np.linalg.LinAlgError: + errori = np.full(k, np.nan) + a_limite = [nomi[i] for i in range(k) if np.isfinite(lo[i]) and + (abs(p[i] - lo[i]) < 1e-3 * abs(lo[i]) or abs(p[i] - hi[i]) < 1e-3 * abs(hi[i]))] + tss = float(np.sum((y - y.mean()) ** 2)) + aicc = n * np.log(max(rss, 1e-300) / n) + 2 * k + (2 * k * (k + 1) / (n - k - 1) if n - k - 1 > 0 else np.inf) + return { + "modello": nome, "equazione": equazione, "n_parametri": k, "n": n, + "parametri": dict(zip(nomi, map(float, p))), "ic95": dict(zip(nomi, map(float, 1.96 * errori))), + "rss": rss, "rmse": float(np.sqrt(rss / n)), "r2": 1 - rss / tss if tss > 0 else float("nan"), + "aicc": float(aicc), "bic": float(n * np.log(max(rss, 1e-300) / n) + k * np.log(n)), + "acf1": acf1, "condizione": cond, "parametri_a_limite": a_limite, + "non_identificabile": bool(cond > f["soglia_condizione"] or + any(errori[i] * 1.96 > abs(p[i]) for i in range(1, k))), "convergenza": bool(migliore.success), + } diff --git a/src/fitseg/modelli.py b/src/fitseg/modelli.py new file mode 100644 index 0000000..f2ad7c9 --- /dev/null +++ b/src/fitseg/modelli.py @@ -0,0 +1,43 @@ +"""Registro dei modelli (t' = t − inizio del dominio). Nuovo modello = nuova voce qui e in docs/SPEC.md.""" + +import numpy as np + + +def _costante(t, p): + return np.full_like(t, p[0]) + + +def _lineare(t, p): + return p[0] + p[1] * t + + +def _quadratico(t, p): + return p[0] + p[1] * t + p[2] * t**2 + + +def _esp_decad(t, p): + return p[0] + p[1] * np.exp(-t / p[2]) + + +def _esp_doppio(t, p): + return p[0] + p[1] * np.exp(-t / p[2]) + p[3] * np.exp(-t / p[4]) + + +def _logistica(t, p): + return p[0] + p[1] / (1 + np.exp(-(t - p[2]) / p[3])) + + +# nome: (funzione, nomi parametri, equazione testuale, indici dei parametri-tempo tau) +MODELLI = { + "costante": (_costante, ["a"], "y = a", []), + "lineare": (_lineare, ["a", "b"], "y = a + b·t'", []), + "quadratico": (_quadratico, ["a", "b", "c"], "y = a + b·t' + c·t'²", []), + "esp_decad": (_esp_decad, ["c", "A", "tau"], "y = c + A·exp(−t'/τ)", [2]), + "esp_doppio": (_esp_doppio, ["c", "A1", "tau1", "A2", "tau2"], + "y = c + A1·exp(−t'/τ1) + A2·exp(−t'/τ2)", [2, 4]), + "logistica": (_logistica, ["c", "A", "t_m", "w"], "y = c + A / (1 + exp(−(t'−t_m)/w))", [3]), +} + + +def valuta(nome, t, p): + return MODELLI[nome][0](np.asarray(t, float), p) diff --git a/src/fitseg/pipeline.py b/src/fitseg/pipeline.py new file mode 100644 index 0000000..8f9e6a0 --- /dev/null +++ b/src/fitseg/pipeline.py @@ -0,0 +1,38 @@ +"""Fit di tutti i modelli abilitati su ogni (dominio, serie) e scelta con la regola automatica.""" + +import numpy as np + +from .domini import sigma_rumore +from .fit import adatta +from .selezione import scegli + + +def analizza_domini(t, serie, domini, cfg): + """Ritorna la lista dei risultati, uno per dominio, con i fit di ogni serie e la scelta.""" + c, f = cfg["domini"], cfg["fit"] + rng = np.random.default_rng(f["seed"]) + risultati = [] + for k, d in enumerate(domini, 1): + i0 = d["i0"] + (c["scarta_dopo_salto_n"] if d["tipo"] != "riposo" else 0) + i1 = d["i1"] + t0 = t[d["i0"]] + voce = {"k": k, "tipo": d["tipo"], "t0_s": float(t[d["i0"]]), "t1_s": float(t[i1 - 1]), + "n": int(i1 - d["i0"]), "n_fit": int(i1 - i0), "serie": {}} + for nome_serie, y in serie.items(): + ty, yy = t[i0:i1] - t0, y[i0:i1] + sigma = sigma_rumore(yy) + voce_s = {"sigma": float(sigma), "fit": [], "scelta": None} + if len(yy) < c["n_min_punti"]: + voce_s["motivo"] = f"indeterminato: {len(yy)} campioni < n_min_punti" + else: + for m in f["modelli_abilitati"]: + from .modelli import MODELLI + if len(MODELLI[m][1]) > len(yy) - 2: + continue + voce_s["fit"].append(adatta(ty, yy, m, sigma, cfg, rng)) + scelto, motivo, adeguato = scegli(voce_s["fit"], sigma, cfg) + voce_s["scelta"] = scelto["modello"] if scelto else None + voce_s["motivo"], voce_s["adeguato"] = motivo, adeguato + voce["serie"][nome_serie] = voce_s + risultati.append(voce) + return risultati diff --git a/src/fitseg/riepilogo.py b/src/fitseg/riepilogo.py new file mode 100644 index 0000000..33221b7 --- /dev/null +++ b/src/fitseg/riepilogo.py @@ -0,0 +1,100 @@ +"""Report, equazioni, modello finale (JSON), valutazione e grafici.""" + +import json + +import numpy as np + +from .modelli import valuta + + +def _scelto(voce_s): + return next((x for x in voce_s["fit"] if x["modello"] == voce_s["scelta"]), None) + + +def _eq(fit): + p = ", ".join(f"{k} = {v:.4g} ± {fit['ic95'][k]:.2g}" for k, v in fit["parametri"].items()) + return f"{fit['equazione']} ({p})" + + +def modello_finale(risultati): + out = [] + for r in risultati: + s = {} + for nome, vs in r["serie"].items(): + fit = _scelto(vs) + s[nome] = ({"modello": fit["modello"], "parametri": fit["parametri"]} if fit + else {"modello": None, "motivo": vs["motivo"]}) + out.append({"tipo": r["tipo"], "t0_s": r["t0_s"], "t1_s": r["t1_s"], "serie": s}) + return out + + +def scrivi_report(cartella, risultati, nome_csv): + righe = [f"# Report fitseg — {nome_csv}\n", + "t' = t − t0 (inizio del dominio). Numeri dal motore; intervalli al 95 %.\n", + "| # | tipo | t0 [s] | t1 [s] | n |", "|---|---|---|---|---|"] + righe += [f"| {r['k']} | {r['tipo']} | {r['t0_s']:.2f} | {r['t1_s']:.2f} | {r['n']} |" for r in risultati] + eq = ["# Equazioni\n"] + for r in risultati: + righe.append(f"\n## Dominio {r['k']} — {r['tipo']} ({r['t0_s']:.2f}–{r['t1_s']:.2f} s)\n") + eq.append(f"## Dominio {r['k']} — {r['tipo']} ({r['t0_s']:.2f}–{r['t1_s']:.2f} s)\n") + for nome, vs in r["serie"].items(): + fit = _scelto(vs) + righe.append(f"### {nome} (σ rumore {vs['sigma']:.3g})\n") + if not fit: + righe.append(f"**Indeterminato** — {vs['motivo']}\n") + eq.append(f"- {nome}: indeterminato\n") + continue + righe.append(f"**Scelto: {fit['modello']}** — {_eq(fit)} \n" + f"RMSE {fit['rmse']:.3g}, R² {fit['r2']:.4f}, ACF1 {fit['acf1']:.2f}" + f"{'' if vs['adeguato'] else ' — **modello non adeguato** (RMSE > soglia·σ)'} \n" + f"Motivo: {vs['motivo']}\n") + righe += ["| modello | k | RMSE | AICc | BIC | ACF1 | valido |", "|---|---|---|---|---|---|---|"] + for x in sorted((x for x in vs["fit"] if "aicc" in x), key=lambda x: x["aicc"]): + ok = not (x["parametri_a_limite"] or x["non_identificabile"] or not x["convergenza"]) + righe.append(f"| {x['modello']} | {x['n_parametri']} | {x['rmse']:.3g} | {x['aicc']:.1f} | " + f"{x['bic']:.1f} | {x['acf1']:.2f} | {'sì' if ok else 'no'} |") + righe.append("") + eq.append(f"- {nome}: `{_eq(fit)}`\n") + (cartella / "report.md").write_text("\n".join(righe), encoding="utf-8") + (cartella / "equazioni.md").write_text("\n".join(eq), encoding="utf-8") + + +def scrivi_json(cartella, risultati): + (cartella / "risultati.json").write_text(json.dumps(risultati, indent=1, ensure_ascii=False), encoding="utf-8") + (cartella / "modello_finale.json").write_text( + json.dumps(modello_finale(risultati), indent=1, ensure_ascii=False), encoding="utf-8") + + +def valuta_finale(modello, nome_serie, t): + """Valuta il modello a tratti in t (array di secondi); NaN fuori dai domini o indeterminati.""" + t = np.asarray(t, float) + y = np.full_like(t, np.nan) + for k, d in enumerate(modello): + ultimo = k == len(modello) - 1 + m = (t >= d["t0_s"]) & ((t <= d["t1_s"]) if ultimo else (t < modello[k + 1]["t0_s"])) + s = d["serie"].get(nome_serie) + if s and s["modello"]: + y[m] = valuta(s["modello"], t[m] - d["t0_s"], list(s["parametri"].values())) + return y + + +def grafici(cartella, t, serie, risultati, domini): + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + colori = {"riposo": "#9aa0a6", "riscaldamento": "#d93025", "raffreddamento": "#1a73e8"} + modello = modello_finale(risultati) + for nome, y in serie.items(): + fig, (a, b) = plt.subplots(2, 1, figsize=(11, 7), sharex=True, height_ratios=[3, 1]) + a.plot(t, y, ".", ms=2, color="#555", label="dati") + yf = valuta_finale(modello, nome, t) + for d in domini: + sl = slice(d["i0"], d["i1"]) + a.axvspan(t[d["i0"]], t[d["i1"] - 1], color=colori[d["tipo"]], alpha=0.12) + a.plot(t[sl], yf[sl], "-", color=colori[d["tipo"]], lw=1.6) + b.plot(t, y - yf, ".", ms=2, color="#555") + b.axhline(0, color="k", lw=0.5) + a.set_ylabel(nome); b.set_ylabel("residuo"); b.set_xlabel("t [s]") + fig.tight_layout() + fig.savefig(cartella / f"fit_{nome}.png", dpi=110) + plt.close(fig) diff --git a/src/fitseg/selezione.py b/src/fitseg/selezione.py new file mode 100644 index 0000000..f10c53d --- /dev/null +++ b/src/fitseg/selezione.py @@ -0,0 +1,25 @@ +"""Regola di scelta del modello (SPEC §5.3).""" + +import math + + +def valido(fit, cfg): + return (fit.get("convergenza") and not fit["parametri_a_limite"] and not fit["non_identificabile"] + and math.isfinite(fit["aicc"])) + + +def scegli(fit_list, sigma, cfg): + """Ritorna (fit scelto, motivo, adeguato) oppure (None, motivo, False) se nessun fit è valido.""" + f = cfg["fit"] + validi = [x for x in fit_list if valido(x, cfg)] + if not validi: + return None, "indeterminato: nessun modello convergente, identificabile e lontano dai limiti", False + adeguati = [x for x in validi if x["rmse"] <= f["rmse_vs_rumore_max"] * sigma] + pool = adeguati or validi + miglior_aicc = min(x["aicc"] for x in pool) + entro = [x for x in pool if x["aicc"] <= miglior_aicc + f["soglia_delta_aicc"]] + scelto = min(entro, key=lambda x: (x["n_parametri"], x["aicc"])) + motivo = (f"{len(validi)}/{len(fit_list)} modelli validi; " + + ("adeguato (RMSE ≤ soglia·σ)" if adeguati else "nessun modello adeguato (RMSE > soglia·σ)") + + f"; meno parametri entro ΔAICc ≤ {f['soglia_delta_aicc']}") + return scelto, motivo, bool(adeguati) diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py new file mode 100644 index 0000000..7e23e97 --- /dev/null +++ b/tests/test_pipeline.py @@ -0,0 +1,40 @@ +"""Verifiche sui CSV generati da gen_data (dati/gen_.csv, config d'esempio: tau di raffreddamento nominale 15.9 s).""" + +from pathlib import Path + +import pytest + +from fitseg import config, dati, domini +from fitseg.pipeline import analizza_domini + +RADICE = Path(__file__).parent.parent +SEEDS = [1, 2, 3] + + +@pytest.fixture(scope="module") +def cfg(): + return config.leggi(RADICE / "config.toml") + + +@pytest.mark.parametrize("seed", SEEDS) +def test_domini_riposo_gradini_raffreddamento(cfg, seed): + t, serie = dati.leggi_csv(RADICE / f"dati/gen_{seed}.csv", cfg) + tipi = [d["tipo"] for d in domini.proponi(t, serie, cfg)] + assert tipi == ["riposo", "riscaldamento", "riscaldamento", "raffreddamento"] + + +@pytest.mark.parametrize("seed", SEEDS) +def test_raffreddamento_esponenziale(cfg, seed): + t, serie = dati.leggi_csv(RADICE / f"dati/gen_{seed}.csv", cfg) + ris = analizza_domini(t, serie, domini.proponi(t, serie, cfg), cfg) + for nome in serie: + vs = ris[-1]["serie"][nome] + assert vs["scelta"] == "esp_decad" + fit = next(x for x in vs["fit"] if x["modello"] == "esp_decad") + assert fit["parametri"]["tau"] == pytest.approx(15.9, rel=0.1) + + +def test_riposo_costante(cfg): + t, serie = dati.leggi_csv(RADICE / "dati/gen_1.csv", cfg) + ris = analizza_domini(t, serie, domini.proponi(t, serie, cfg), cfg) + assert ris[0]["serie"]["dT_C"]["scelta"] == "costante"