Riscaldamento e raffreddamento esponenziali con c e tau condivisi; grafici per dominio nel report
- modelli per tipo di dominio ([fit.modelli_per_tipo]) - nei domini di riscaldamento c e tau fissati ai valori del raffreddamento - un grafico dati/modelli per dominio e serie, incorporato in report.md Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
This commit is contained in:
1 parent
d4fbd02a51
commit
b8c8aeb60e
11 files changed
+167
-58
No files matched your search
@@ -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.
|
||||
@@ -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_<serie>.png` | dati, fit a tratti (un colore per tipo di dominio) e residui |
|
||||
| `out/fit_<serie>.png` | panoramica: dati, fit a tratti (un colore per tipo di dominio) e residui |
|
||||
| `out/dominio<k>_<serie>.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.
|
||||
+6
-1
@@ -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"
|
||||
+13
-4
@@ -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_<serie>.png` | dati, fit a tratti (colore per tipo di dominio) e residui |
|
||||
| `out/fit_<serie>.png` | panoramica: dati, fit a tratti (colore per tipo di dominio) e residui |
|
||||
| `out/dominio<k>_<serie>.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_<seed>.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.
|
||||
+2
-2
@@ -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`.
|
||||
+2
-1
@@ -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):
|
||||
|
||||
@@ -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"]:
|
||||
|
||||
+34
-14
@@ -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),
|
||||
}
|
||||
+35
-24
@@ -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
|
||||
+44
-3
@@ -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"\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)
|
||||
@@ -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
|
||||
Reference in new issue
Block a user