diff --git a/CLAUDE.md b/CLAUDE.md index 1da48a0..776b053 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -29,13 +29,14 @@ Flusso: `domini.proponi` → `pipeline.analizza_domini` (per ogni dominio e seri - [domini.py](src/fitseg/domini.py): i domini derivano **solo dai gradini positivi** della `serie_riferimento` (`diff > soglia_gradino_c`), con arretramento ai campioni di transizione; non c'è più PELT/fusione BIC (presente nelle prime versioni della spec). `domini.toml` è modificabile e rientra con `--domini`; deve essere contiguo. - [modelli.py](src/fitseg/modelli.py): registro `MODELLI` (funzione, parametri, equazione testuale, indici dei τ). Nuovo modello = voce qui + riga in SPEC §4. - [fit.py](src/fitseg/fit.py): `least_squares` con perdita robusta; gli esponenziali partono da una griglia su τ (VARPRO). "Non identificabile" = condizione > soglia **o** IC 95 % > |valore| per un parametro diverso dal primo; la `logistica` ha `t_m` limitato a [0, L] altrimenti degenera in esponenziale. +- [pipeline.py](src/fitseg/pipeline.py): i domini di raffreddamento si fittano **per primi** e i loro `c`/`tau` vengono fissati nei domini di riscaldamento (`riscaldamento_fissa_da_raffreddamento`, `adatta(..., fissi=)`); i modelli provati dipendono dal tipo (`[fit.modelli_per_tipo]`). - [selezione.py](src/fitseg/selezione.py): regola SPEC §5.3 (validi → adeguati se RMSE ≤ soglia·σ → entro ΔAICc vince chi ha meno parametri). - σ del rumore è stimata **per dominio** da `MAD(d²y)`; il tempo del fit è traslato a `t' = t − t0` del dominio. ## Gotcha - `scarta_dopo_salto_n` deve coprire i campioni di transizione del sensore (τ_sensore 0.07 s, dt 0.15 s ⇒ ≥ 2–3 campioni); se troppo basso, 1–2 punti a ~100 °C dal fit sporcano la fine del dominio. -- Domini di riscaldamento di 14–26 campioni: τ non identificabile; il modello scelto (`lineare`, a volte `logistica`) è solo descrittivo. +- Domini di riscaldamento di 14–26 campioni: τ non identificabile da solo, quindi è fisso al valore del raffreddamento; gli IC di `A` sono condizionati a quei fissi. Se il raffreddamento non è valido, il riscaldamento ricade sul fit libero e può risultare indeterminato. - Cambiando `config.toml` del generatore (gradini, rumore) ricontrolla `soglia_gradino_c` e i test in [tests/test_pipeline.py](tests/test_pipeline.py), che usano `dati/gen_{1,2,3}.csv`. - Il parser CSV richiede la colonna tempo e tutte le `serie` in config; il tempo deve essere strettamente crescente (exit 4 altrimenti). - Exit code: 0 ok · 2 uso · 3 config · 4 CSV · 6 domini non validi. diff --git a/README.md b/README.md index f8eb638..5623676 100644 --- a/README.md +++ b/README.md @@ -77,11 +77,12 @@ Tutto in `lavoro/` (ignorata da git): |---|---| | `domini.toml` | proposta dei domini (o quelli che hai usato) | | `dati.csv` | copia del CSV analizzato | -| `out/report.md` | per dominio e serie: modello scelto, parametri ± IC 95 %, metriche, tabella di confronto di tutti i modelli | +| `out/report.md` | con i grafici incorporati (si vede bene in un visualizzatore Markdown); per dominio e serie: modello scelto, parametri ± IC 95 %, metriche, tabella di confronto di tutti i modelli | | `out/equazioni.md` | solo le equazioni scelte | | `out/modello_finale.json` | modello e parametri per dominio e serie (input di `valuta`) | | `out/risultati.json` | tutti i fit con tutte le metriche | -| `out/fit_.png` | dati, fit a tratti (un colore per tipo di dominio) e residui | +| `out/fit_.png` | panoramica: dati, fit a tratti (un colore per tipo di dominio) e residui | +| `out/dominio_.png` | un grafico per dominio e serie: dati grezzi contro il modello scelto (rosso) e gli altri modelli validi (tratteggiati), più residui; i campioni esclusi dal fit sono in grigio | ### Come leggere il report @@ -93,7 +94,7 @@ Tutto in `lavoro/` (ignorata da git): ## Come funziona 1. **Domini**: un salto positivo tra due campioni consecutivi sopra `soglia_gradino_c` è un gradino di riscaldamento. Prima del primo gradino il dominio è `riposo`, tra due gradini `riscaldamento`, dopo l'ultimo `raffreddamento`. -2. **Fit**: in ogni dominio, con il tempo traslato a `t' = t − t0`, si fittano tutti i modelli di `modelli_abilitati` con `scipy.optimize.least_squares` e perdita robusta. +2. **Fit**: in ogni dominio, con il tempo traslato a `t' = t − t0`, si fittano i modelli del suo tipo con `scipy.optimize.least_squares` e perdita robusta. Riscaldamento e raffreddamento sono esponenziali; nel riscaldamento `c` e `τ` sono fissati ai valori del raffreddamento (i tratti sono troppo corti per stimarli) e si stima solo l'ampiezza. 3. **Scelta**: regola automatica (SPEC §5.3), senza intervento umano. Dettagli, modelli e limiti: [docs/SPEC.md](docs/SPEC.md). Uso esteso: [docs/USO.md](docs/USO.md). @@ -109,7 +110,7 @@ Dettagli, modelli e limiti: [docs/SPEC.md](docs/SPEC.md). Uso esteso: [docs/USO. | `esp_doppio` | `y = c + A1·exp(−t'/τ1) + A2·exp(−t'/τ2)` | | `logistica` | `y = c + A / (1 + exp(−(t'−t_m)/w))` | -Per aggiungerne uno: nuova voce in `MODELLI` di [src/fitseg/modelli.py](src/fitseg/modelli.py), nuova riga in SPEC §4 e nome in `modelli_abilitati`. +Per aggiungerne uno: nuova voce in `MODELLI` di [src/fitseg/modelli.py](src/fitseg/modelli.py), nuova riga in SPEC §4 e nome nelle liste di `[fit.modelli_per_tipo]` dei tipi di dominio in cui va provato. ## Configurazione @@ -120,7 +121,8 @@ Tutti i valori sono in [config.toml](config.toml); non c'è nessun default nel c | `[domini] soglia_gradino_c` | salto minimo che segna un gradino (20 °C); da adattare se i gradini reali sono più piccoli | | `[domini] scarta_dopo_salto_n` | campioni esclusi all'inizio di un dominio dopo un gradino (transizione del sensore) | | `[domini] n_min_punti` | campioni minimi per dominio | -| `[fit] modelli_abilitati` | modelli provati in ogni dominio | +| `[fit.modelli_per_tipo]` | modelli provati per tipo di dominio (riposo: costante/lineare/esponenziale; riscaldamento e raffreddamento: esponenziale) | +| `[fit] riscaldamento_fissa_da_raffreddamento` | parametri dei domini di riscaldamento fissati al valore del raffreddamento (`["c", "tau"]`; `[]` = liberi) | | `[fit] soglia_delta_aicc` | entro questa differenza vince il modello con meno parametri | | `[fit] rmse_vs_rumore_max` | soglia di «adeguato» (RMSE ≤ valore × σ rumore) | | `[fit] seed` | seme delle partenze multiple: stessi dati e stessa config danno gli stessi risultati | @@ -131,7 +133,7 @@ Tutti i valori sono in [config.toml](config.toml); non c'è nessun default nel c .venv/bin/pytest # tutti .venv/bin/pytest -k raffreddamento # uno solo ``` -I test verificano, sui seed 1–3, che i domini siano `riposo, riscaldamento, riscaldamento, raffreddamento`, che nel raffreddamento vinca `esp_decad` con τ entro il 10 % di 15.9 s e che il riposo sia `costante`. Richiedono `dati/gen_1.csv`, `gen_2.csv`, `gen_3.csv` (vedi sopra). +I test verificano, sui seed 1–3, che i domini siano `riposo, riscaldamento, riscaldamento, raffreddamento`, che nel raffreddamento vinca `esp_decad` con τ entro il 10 % di 15.9 s, che nel riscaldamento sia `esp_decad` con `c` e `τ` del raffreddamento e che il riposo sia `costante`. Richiedono `dati/gen_1.csv`, `gen_2.csv`, `gen_3.csv` (vedi sopra). ## Codici di uscita @@ -140,6 +142,6 @@ I test verificano, sui seed 1–3, che i domini siano `riposo, riscaldamento, ri ## Limiti noti - Nel raffreddamento i residui sono autocorrelati (ACF1 ≈ 0.7) per scarti casuali e deriva lenta del sensore: i parametri sono buoni, ma gli intervalli di confidenza sono troppo stretti e il dominio risulta «non adeguato». -- I domini di riscaldamento sono corti (2–4 s): il decadimento è quasi lineare e τ non è identificabile; vince un modello semplice, con valore descrittivo e non fisico. +- I domini di riscaldamento sono corti (2–4 s): τ non è identificabile da solo e viene preso dal raffreddamento; gli IC di `A` sono condizionati a quei valori fissi. - Ogni dominio è fittato in modo indipendente: nessun vincolo di continuità tra domini. - Provato su 3 seed; non su file senza gradini o con un solo gradino. diff --git a/config.toml b/config.toml index dc56f9f..e3c3eeb 100644 --- a/config.toml +++ b/config.toml @@ -16,7 +16,7 @@ n_min_punti = 6 # campioni minimi per dominio (sotto: "indete 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"] +riscaldamento_fissa_da_raffreddamento = ["c", "tau"] # parametri dei domini di riscaldamento fissati ai valori del fit del raffreddamento ([] = liberi) 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 @@ -30,5 +30,10 @@ soglia_delta_aicc = 2.0 # entro questa differenza vince il modello co rmse_vs_rumore_max = 1.3 # RMSE <= valore × sigma rumore => modello adeguato seed = 0 # seme delle partenze multiple +[fit.modelli_per_tipo] # modelli provati in ogni tipo di dominio (riscaldamento e raffreddamento: esponenziali) +riposo = ["costante", "lineare", "esp_decad"] +riscaldamento = ["esp_decad"] +raffreddamento = ["esp_decad"] + [output] cartella_lavoro = "lavoro" diff --git a/docs/SPEC.md b/docs/SPEC.md index a973d49..d3be3e6 100644 --- a/docs/SPEC.md +++ b/docs/SPEC.md @@ -59,7 +59,7 @@ I domini devono essere contigui; in caso contrario exit 6. | `esp_doppio` | `y = c + A1·exp(−t'/τ1) + A2·exp(−t'/τ2)`, τ1 < τ2 | c, A1, τ1, A2, τ2 | | `logistica` | `y = c + A / (1 + exp(−(t'−t_m)/w))`, 0 ≤ t_m ≤ L | c, A, t_m, w | -Nuovo modello = nuova voce nel registro `MODELLI` e nuova riga qui. Nel raffreddamento `c` è la temperatura asintotica (≈ ambiente) e +Nuovo modello = nuova voce nel registro `MODELLI` e nuova riga qui. I modelli provati dipendono dal tipo di dominio (`[fit.modelli_per_tipo]`): `riposo` → `costante`, `lineare`, `esp_decad`; `riscaldamento` e `raffreddamento` → solo `esp_decad` (sono decadimenti esponenziali). Nel raffreddamento `c` è la temperatura asintotica (≈ ambiente) e `τ = 1/k`; nel riposo `costante`/`lineare` descrivono la baseline. ## 5. Fit @@ -72,6 +72,13 @@ Nuovo modello = nuova voce nel registro `MODELLI` e nuova riga qui. Nel raffredd 5. Metriche: RSS, RMSE, R², AICc, BIC, ACF1 dei residui, parametri a limite, numero di condizione (jacobiano normalizzato). **Non identificabile** se condizione > `soglia_condizione` oppure un parametro (escluso il primo) ha IC 95 % > |valore|. +### 5.2 Costanti condivise tra riscaldamento e raffreddamento +Dopo ogni gradino la temperatura decade con la stessa legge del raffreddamento finale (stessa costante di tempo, stesso asintoto ≈ ambiente). +I domini di raffreddamento si fittano per primi; nei domini di riscaldamento i parametri di `riscaldamento_fissa_da_raffreddamento` +(default `c` e `tau`) sono **fissati** ai valori del raffreddamento e resta libera solo l'ampiezza `A`. Motivo: in 2–4 s di dati τ non si +identifica da solo. Gli IC 95 % di `A` sono condizionati ai valori fissi (l'incertezza di `c` e `τ` non si propaga); i parametri fissi sono +marcati "(fisso)" nel report. Se manca un raffreddamento valido, o la lista è vuota, i parametri restano liberi. Con `[]` si torna al fit libero. + ### 5.3 Regola di scelta (`selezione.py`) Scarta non convergenti, con parametri a limite o non identificabili. Se esistono modelli con `RMSE ≤ rmse_vs_rumore_max·σ` ("adeguati") si sceglie tra questi, altrimenti tra tutti i validi; tra i modelli entro `soglia_delta_aicc` dal miglior AICc vince quello con meno @@ -99,7 +106,8 @@ Exit code: 0 ok · 2 uso errato · 3 config non valida · 4 CSV illeggibile · 6 | `out/risultati.json` | tutti i fit per (dominio, serie), con metriche e scelta | | `out/modello_finale.json` | per dominio e serie: modello scelto e parametri | | `out/report.md`, `out/equazioni.md` | tabelle di confronto, equazioni con IC 95 % | -| `out/fit_.png` | dati, fit a tratti (colore per tipo di dominio) e residui | +| `out/fit_.png` | panoramica: dati, fit a tratti (colore per tipo di dominio) e residui | +| `out/dominio_.png` | per dominio e serie: dati grezzi vs modello scelto e altri modelli validi, residui; incorporati nel `report.md` | ## 8. Test (`uv run pytest`, dati in `dati/gen_.csv`) @@ -107,6 +115,7 @@ Exit code: 0 ok · 2 uso errato · 3 config non valida · 4 CSV illeggibile · 6 |---|---|---| | T1 | Domini sui seed 1–3 | `riposo, riscaldamento, riscaldamento, raffreddamento` ✔ | | T2 | Raffreddamento | `esp_decad` scelto per T e ΔT, τ entro il 10 % di 15.9 s ✔ (seed 1: τ = 16.08 ± 0.07 s) | +| T2b | Riscaldamento | `esp_decad` con `c` e `τ` fissi uguali a quelli del raffreddamento, `A` > 0 ✔ | | T3 | Riposo | `costante` per ΔT ✔ | | — | da fare | più seed, robustezza a config diverse (nessun gradino, un solo gradino), `domini.toml` manuale, riproducibilità | @@ -114,7 +123,7 @@ Exit code: 0 ok · 2 uso errato · 3 config non valida · 4 CSV illeggibile · 6 - **Residui autocorrelati nel raffreddamento** (ACF1 ≈ 0.7, RMSE ≈ 0.36 contro σ ≈ 0.15): scarti casuali e deriva lenta non sono nel modello; i parametri sono buoni ma gli IC 95 % sono troppo stretti (nonostante la correzione) e il dominio risulta "non adeguato". -- **Domini di riscaldamento corti** (14–26 campioni, 2–4 s): il decadimento è quasi lineare, τ non identificabile; vince un modello - semplice (`lineare`, talvolta `logistica`, che qui è solo una curva flessibile e non ha significato fisico). +- **Domini di riscaldamento corti** (14–26 campioni, 2–4 s): τ non è identificabile da solo, per questo è fissato al valore del raffreddamento (§5.2). + Il secondo tratto di riscaldamento ha ACF1 ≈ 0.7 e RMSE ≈ 2 σ (scarti casuali): risulta "non adeguato" anche se la forma è giusta. - La soglia `soglia_gradino_c` va adattata se i gradini reali sono più piccoli di 20 °C o gli scarti più grandi. - Il vincolo di continuità tra domini è solo ignorato: ogni dominio è fittato in modo indipendente. diff --git a/docs/USO.md b/docs/USO.md index b74cbbe..e210afb 100644 --- a/docs/USO.md +++ b/docs/USO.md @@ -19,9 +19,9 @@ cd ../gen_data && .venv/bin/python genera.py config.toml --seed 4 -o ../Equation .venv/bin/fitseg valuta lavoro/out/modello_finale.json --t 10 20 50 .venv/bin/pytest # test (singolo: -k raffreddamento) ``` -Il report è in `lavoro/out/report.md`: per ogni dominio e serie, modello scelto con parametri ± IC 95 %, metriche e tabella di +Il report è in `lavoro/out/report.md` (con i grafici dati-contro-equazioni incorporati, uno per dominio e serie): per ogni dominio e serie, modello scelto con parametri ± IC 95 %, metriche e tabella di confronto di tutti i modelli. Leggi sempre la riga "non adeguato" e l'ACF1: indicano residui che il modello non spiega. ## 4. Configurazione Tutto in [config.toml](../config.toml) (nessun default nel codice). Chiavi di rilievo: `[domini] soglia_gradino_c`, -`scarta_dopo_salto_n`, `n_min_punti`; `[fit] modelli_abilitati`, `soglia_delta_aicc`, `rmse_vs_rumore_max`, `seed`. +`scarta_dopo_salto_n`, `n_min_punti`; `[fit.modelli_per_tipo]`, `riscaldamento_fissa_da_raffreddamento`, `soglia_delta_aicc`, `rmse_vs_rumore_max`, `seed`. diff --git a/src/fitseg/cli.py b/src/fitseg/cli.py index 323a530..66279f5 100644 --- a/src/fitseg/cli.py +++ b/src/fitseg/cli.py @@ -68,10 +68,11 @@ def cmd_analizza(args): riepilogo.scrivi_json(out, risultati) riepilogo.scrivi_report(out, risultati, Path(args.csv).name) riepilogo.grafici(out, t, serie, risultati, domini) + riepilogo.grafici_domini(out, t, serie, risultati, domini, cfg) 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) + print(f"scritto: {out}/ (report.md, equazioni.md, modello_finale.json, fit_*.png, dominio*_*.png)", file=sys.stderr) def cmd_modelli(args): diff --git a/src/fitseg/config.py b/src/fitseg/config.py index f91767e..bf87fc1 100644 --- a/src/fitseg/config.py +++ b/src/fitseg/config.py @@ -6,7 +6,7 @@ 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", + "fit": ["modelli_per_tipo", "riscaldamento_fissa_da_raffreddamento", "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"], @@ -28,6 +28,10 @@ def leggi(percorso): 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)}") + mpt = cfg["fit"]["modelli_per_tipo"] + for tipo in ("riposo", "riscaldamento", "raffreddamento"): + if not mpt.get(tipo): + raise ErroreConfig(f"[fit.modelli_per_tipo] manca la lista di modelli per '{tipo}'") 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"]: diff --git a/src/fitseg/fit.py b/src/fitseg/fit.py index 4ccbb2f..365311c 100644 --- a/src/fitseg/fit.py +++ b/src/fitseg/fit.py @@ -61,28 +61,45 @@ def _limiti(nome, cfg, 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.""" +def adatta(t, y, nome, sigma, cfg, rng, fissi=None): + """Fit di `nome` su (t, y). `fissi` = {parametro: valore} tenuti costanti (non contano nei gradi di libertà).""" f = cfg["fit"] funzione, nomi, equazione, idx_tau = MODELLI[nome] - n, k = len(t), len(nomi) + fissi = {k: v for k, v in (fissi or {}).items() if k in nomi} + n, k_tot = len(t), len(nomi) + libero = [i for i, nm in enumerate(nomi) if nm not in fissi] + k = len(libero) 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 + + def completo(x): + p = np.empty(k_tot) + p[libero] = x + for i, nm in enumerate(nomi): + if nm in fissi: + p[i] = fissi[nm] + return p + + res_fn = lambda x: funzione(t, completo(x)) - y migliore = None for p0 in _partenze(t, y, nome, cfg, L, rng): - p0 = np.clip(p0, lo, hi) + p0 = completo(np.asarray(p0)[libero]) + if nome == "esp_decad" and libero == [1]: # c e tau fissi: l'ampiezza è lineare + e = np.exp(-t / p0[2]) + p0[1] = np.sum((y - p0[0]) * e) / np.sum(e * e) + x0 = np.clip(p0[libero], lo[libero], hi[libero]) try: - r = least_squares(res_fn, p0, bounds=(lo, hi), loss=f["perdita"], f_scale=scala, x_scale="jac") + r = least_squares(res_fn, x0, bounds=(lo[libero], hi[libero]), 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 + p, J = completo(migliore.x), migliore.jac + if nome == "esp_doppio" and p[2] > p[4] and not fissi: # tau1 < tau2 perm = [0, 3, 4, 1, 2] p, J = p[perm], J[:, perm] r = funzione(t, p) - y @@ -94,19 +111,22 @@ def adatta(t, y, nome, sigma, cfg, rng): 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 + err_liberi = 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 + err_liberi = np.full(k, np.nan) + errori = np.zeros(k_tot) + errori[libero] = err_liberi + a_limite = [nomi[i] for i in libero 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]))] + non_ident = cond > f["soglia_condizione"] or any( + errori[i] * 1.96 > abs(p[i]) for i in libero if nomi[i] not in ("a", "c")) 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, + "modello": nome, "equazione": equazione, "n_parametri": k, "n": n, "fissi": sorted(fissi), "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), + "non_identificabile": bool(non_ident), "convergenza": bool(migliore.success), } diff --git a/src/fitseg/pipeline.py b/src/fitseg/pipeline.py index 8f9e6a0..d769fc7 100644 --- a/src/fitseg/pipeline.py +++ b/src/fitseg/pipeline.py @@ -1,38 +1,49 @@ -"""Fit di tutti i modelli abilitati su ogni (dominio, serie) e scelta con la regola automatica.""" +"""Fit dei modelli di ogni tipo di dominio su ogni (dominio, serie) e scelta con la regola automatica.""" import numpy as np from .domini import sigma_rumore from .fit import adatta +from .modelli import MODELLI from .selezione import scegli +def _analizza_serie(ty, yy, sigma, tipo, fissi, cfg, rng): + c, f = cfg["domini"], cfg["fit"] + voce = {"sigma": float(sigma), "fit": [], "scelta": None} + if len(yy) < c["n_min_punti"]: + voce["motivo"] = f"indeterminato: {len(yy)} campioni < n_min_punti" + return voce + for m in f["modelli_per_tipo"][tipo]: + if len(MODELLI[m][1]) > len(yy) - 2: + continue + voce["fit"].append(adatta(ty, yy, m, sigma, cfg, rng, fissi if tipo == "riscaldamento" else None)) + scelto, motivo, adeguato = scegli(voce["fit"], sigma, cfg) + voce["scelta"] = scelto["modello"] if scelto else None + voce["motivo"], voce["adeguato"] = motivo, adeguato + return voce + + def analizza_domini(t, serie, domini, cfg): - """Ritorna la lista dei risultati, uno per dominio, con i fit di ogni serie e la scelta.""" + """Ritorna un risultato per dominio. I domini di raffreddamento si fittano per primi: i loro parametri + (`riscaldamento_fissa_da_raffreddamento`) vengono fissati nei domini di riscaldamento (stesse costanti fisiche).""" c, f = cfg["domini"], cfg["fit"] rng = np.random.default_rng(f["seed"]) - risultati = [] - for k, d in enumerate(domini, 1): + risultati = [None] * len(domini) + fissi = {nome: {} for nome in serie} # nome serie -> parametri del raffreddamento + ordine = sorted(range(len(domini)), key=lambda i: domini[i]["tipo"] != "raffreddamento") + for idx in ordine: + d = domini[idx] 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]), + i1, t0 = d["i1"], t[d["i0"]] + voce = {"k": idx + 1, "tipo": d["tipo"], "t0_s": float(t0), "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) + for nome, y in serie.items(): + yy = y[i0:i1] + vs = _analizza_serie(t[i0:i1] - t0, yy, sigma_rumore(yy), d["tipo"], fissi[nome], cfg, rng) + voce["serie"][nome] = vs + if d["tipo"] == "raffreddamento" and vs["scelta"] and not fissi[nome]: + par = next(x for x in vs["fit"] if x["modello"] == vs["scelta"])["parametri"] + fissi[nome] = {k: par[k] for k in f["riscaldamento_fissa_da_raffreddamento"] if k in par} + risultati[idx] = voce return risultati diff --git a/src/fitseg/riepilogo.py b/src/fitseg/riepilogo.py index 33221b7..8f76731 100644 --- a/src/fitseg/riepilogo.py +++ b/src/fitseg/riepilogo.py @@ -12,7 +12,8 @@ def _scelto(voce_s): def _eq(fit): - p = ", ".join(f"{k} = {v:.4g} ± {fit['ic95'][k]:.2g}" for k, v in fit["parametri"].items()) + p = ", ".join(f"{k} = {v:.4g} (fisso)" if k in fit.get("fissi", []) else f"{k} = {v:.4g} ± {fit['ic95'][k]:.2g}" + for k, v in fit["parametri"].items()) return f"{fit['equazione']} ({p})" @@ -30,8 +31,9 @@ def modello_finale(risultati): 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 |", "|---|---|---|---|---|"] + "t' = t − t0 (inizio del dominio). Numeri dal motore; intervalli al 95 %.\n"] + righe += [f"![{nome}](fit_{nome}.png)\n" for nome in risultati[0]["serie"]] + righe += ["| # | 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: @@ -40,6 +42,7 @@ def scrivi_report(cartella, risultati, nome_csv): for nome, vs in r["serie"].items(): fit = _scelto(vs) righe.append(f"### {nome} (σ rumore {vs['sigma']:.3g})\n") + righe.append(f"![dominio {r['k']} {nome}](dominio{r['k']}_{nome}.png)\n") if not fit: righe.append(f"**Indeterminato** — {vs['motivo']}\n") eq.append(f"- {nome}: indeterminato\n") @@ -98,3 +101,41 @@ def grafici(cartella, t, serie, risultati, domini): fig.tight_layout() fig.savefig(cartella / f"fit_{nome}.png", dpi=110) plt.close(fig) + + +def grafici_domini(cartella, t, serie, risultati, domini, cfg): + """Un PNG per (dominio, serie): dati grezzi contro il modello scelto e gli altri modelli validi, più residui.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + from .selezione import valido + scarta = cfg["domini"]["scarta_dopo_salto_n"] + for r, d in zip(risultati, domini): + i0 = d["i0"] + (scarta if d["tipo"] != "riposo" else 0) + sl_tutti, sl_fit = slice(d["i0"], d["i1"]), slice(i0, d["i1"]) + tt = np.linspace(t[i0], t[d["i1"] - 1], 400) + for nome, vs in r["serie"].items(): + y = serie[nome] + fig, (a, b) = plt.subplots(2, 1, figsize=(9, 6), sharex=True, height_ratios=[3, 1]) + a.plot(t[sl_tutti], y[sl_tutti], ".", ms=4, color="#aaa", label="dati (esclusi dal fit)") + a.plot(t[sl_fit], y[sl_fit], ".", ms=5, color="#222", label="dati grezzi") + scelto = _scelto(vs) + for x in vs["fit"]: + if "parametri" not in x or x is scelto or not valido(x, cfg): + continue + a.plot(tt, valuta(x["modello"], tt - r["t0_s"], list(x["parametri"].values())), "--", lw=1, + alpha=0.8, label=f"{x['modello']} (RMSE {x['rmse']:.3g})") + if scelto: + p = list(scelto["parametri"].values()) + a.plot(tt, valuta(scelto["modello"], tt - r["t0_s"], p), "-", lw=2.2, color="#d93025", + label=f"scelto: {scelto['modello']} (RMSE {scelto['rmse']:.3g})") + b.plot(t[sl_fit], y[sl_fit] - valuta(scelto["modello"], t[sl_fit] - r["t0_s"], p), ".", ms=4, color="#222") + else: + a.set_title("indeterminato: nessun modello valido") + b.axhline(0, color="k", lw=0.5) + a.set_ylabel(nome); b.set_ylabel("residuo"); b.set_xlabel("t [s]") + a.set_title(f"Dominio {r['k']} — {r['tipo']} — {nome}", fontsize=10) + a.legend(fontsize=8) + fig.tight_layout() + fig.savefig(cartella / f"dominio{r['k']}_{nome}.png", dpi=110) + plt.close(fig) diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py index 7e23e97..e6e7be8 100644 --- a/tests/test_pipeline.py +++ b/tests/test_pipeline.py @@ -38,3 +38,18 @@ 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" + + +@pytest.mark.parametrize("seed", SEEDS) +def test_riscaldamento_esponenziale_con_costanti_del_raffreddamento(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: + raff = next(x for x in ris[-1]["serie"][nome]["fit"] if x["modello"] == "esp_decad")["parametri"] + for r in ris[1:3]: + vs = r["serie"][nome] + assert vs["scelta"] == "esp_decad" + fit = next(x for x in vs["fit"] if x["modello"] == "esp_decad") + assert fit["fissi"] == ["c", "tau"] + assert fit["parametri"]["tau"] == raff["tau"] and fit["parametri"]["c"] == raff["c"] + assert fit["parametri"]["A"] > 0