Programma fitseg: domini di riscaldamento/raffreddamento, fit e report
Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
This commit is contained in:
1 parent
55722f0e13
commit
70231607bb
14 files changed
+687
No files matched your search
@@ -0,0 +1,6 @@
|
|||||||
|
.venv/
|
||||||
|
lavoro/
|
||||||
|
__pycache__/
|
||||||
|
*.egg-info/
|
||||||
|
.pytest_cache/
|
||||||
|
dati/*.csv
|
||||||
+34
@@ -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"
|
||||||
@@ -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"]
|
||||||
Whitespace-only changes.
@@ -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()
|
||||||
@@ -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
|
||||||
@@ -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
|
||||||
@@ -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
|
||||||
@@ -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),
|
||||||
|
}
|
||||||
@@ -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)
|
||||||
@@ -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
|
||||||
@@ -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)
|
||||||
@@ -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)
|
||||||
@@ -0,0 +1,40 @@
|
|||||||
|
"""Verifiche sui CSV generati da gen_data (dati/gen_<seed>.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"
|
||||||
Reference in new issue
Block a user