Files
EquationMiner/docs/algoritmi.md
T

16 KiB
Raw Blame History

Algoritmi

Tutti i passi, nell'ordine in cui il programma li esegue. Le assunzioni richiamate (A1, A2, …) sono in assunzioni.md; le chiavi di configurazione (in_questo_formato) in configurazione.md; il modulo che implementa ogni passo in moduli.md.

Visione d'insieme

flowchart TD
    CSV(["CSV in input/"]) --> P1

    subgraph P1["PASSO 1 · Domini"]
        direction LR
        L["<b>1 · Lettura</b><br/>validazione"] --> S["<b>2 · Rumore</b><br/>stima di σ"] --> G["<b>3 · Gradini</b><br/>salite sopra soglia"] --> D["<b>4 · Domini</b><br/>divisione del tempo"]
    end

    P1 --> RIP["<b>Riposo</b>"]
    P1 --> RIS["<b>Riscaldamento</b><br/>2–4 campioni"]
    P1 --> INT["<b>Raffreddamenti intermedi</b><br/>tra due gradini, corti"]
    P1 --> FIN["<b>Raffreddamento finale</b><br/>l'ultimo, il più lungo"]

    subgraph P2["PASSO 2 · Equazioni"]
        direction TB
        F0["<b>5 · Fit</b><br/>costante, lineare,<br/>esponenziale"]
        F1["<b>5 · Fit esponenziale</b><br/>c, A, τ liberi"]
        F2["<b>5b · Fit esponenziale</b><br/>c e τ fissi, solo A libera"]
        SALTO["<b>7 · Salto Δ</b><br/>descritto, non fittato"]
        SC["<b>6 · Scelta del modello</b><br/>validi → adeguati → meno parametri"]
        F0 --> SC
        F1 --> SC
        F2 --> SC
    end

    RIP --> F0
    FIN --> F1
    INT --> F2
    RIS --> SALTO
    F1 -. "passa c, τ" .-> F2

    SC --> REP
    SALTO --> REP
    REP(["<b>PASSO 3</b> · Sovrapposizione e report<br/>dati grezzi contro equazioni<br/>report.md · grafici · risultati.json"])

    classDef dato fill:#eceff1,stroke:#607d8b,color:#263238
    classDef calcolo fill:#e3f2fd,stroke:#1a73e8,color:#0d47a1
    classDef dominio fill:#fff8e1,stroke:#f9a825,color:#5d4037
    classDef uscita fill:#e8f5e9,stroke:#2e7d32,color:#1b5e20
    class CSV dato
    class L,S,G,D,F0,F1,F2,SC,SALTO calcolo
    class RIP,RIS,INT,FIN dominio
    class REP uscita

    linkStyle 15 stroke:#e65100,stroke-width:2px,color:#e65100

Come si legge: dall'alto in basso. I numeri nei riquadri sono i paragrafi di questo documento. Colori: grigio dati, azzurro calcoli, giallo tipi di dominio, verde risultati. Il Passo 1 produce i domini; ogni tipo ha il suo trattamento nel Passo 2; la freccia arancione tratteggiata indica che il raffreddamento finale passa c e \tau ai raffreddamenti intermedi (§5b), per questo si elabora per primo. Il riscaldamento non si fitta: se ne riporta il salto (§7).

Simboli e chiavi di configurazione

Nelle formule si usano questi simboli; a destra la chiave di configurazione.md con il valore.

Simbolo Chiave Significato
\theta_g soglia_gradino_c salita complessiva minima che segna un gradino [°C]
\theta_b soglia_bordo_sigma soglia, in multipli di \sigma_d, sopra la quale un campione è «in salita»
\theta_\rho soglia_acf soglia di \rho_1 oltre la quale si allargano gli intervalli
\kappa scala_perdita_sigma scala della perdita robusta, in multipli di \sigma
\alpha_{\min},\ \alpha_{\max} tau_min_rel, tau_max_rel limiti di \tau come frazione di L
\lambda scala_residui_vs_rumore_max un modello è «adeguato» se \hat\sigma_\text{res} \le \lambda\,\sigma (scala robusta dei residui, §5.5)
\delta soglia_delta_aicc tolleranza su \Delta\mathrm{AICc} per preferire il modello più semplice
L — durata del dominio, t_\text{ultimo} - t_\text{primo}
n, k — campioni usati nel fit; parametri liberi del modello

1. Lettura dei dati (dati.leggi_csv)

Legge l'intestazione, controlla che ci siano la colonna del tempo e tutte le serie (A2), converte il tempo in secondi (A3) e verifica che sia strettamente crescente (A1). Restituisce t (array in secondi) e un dizionario serie {nome: array}.

2. Stima del rumore (domini.sigma_rumore)

Per un rumore bianco \varepsilon_i di deviazione \sigma, la derivata seconda discreta \Delta^2 y_i = y_{i+2} - 2y_{i+1} + y_i elimina le componenti lisce (costante, lineare, e in buona parte l'esponenziale) e ha varianza 6\sigma^2. Per renderla robusta agli outlier si usa la deviazione mediana assoluta:

\sigma = \frac{1.4826\,\mathrm{MAD}(\Delta^2 y)}{\sqrt{6}}, \qquad \mathrm{MAD}(x) = \mathrm{mediana}\,|x_i - \mathrm{mediana}(x)|

Il fattore 1.4826 rende la MAD uguale a \sigma per dati gaussiani. La stima si calcola per dominio (A13) sui campioni usati nel fit e serve come scala della perdita robusta e come riferimento per l'«adeguatezza» (§6). Ignora la deriva lenta (A14).

Servono almeno 8 campioni (domini.N_MIN_SIGMA): con meno, la MAD di pochi \Delta^2 y può annullarsi, perché i valori sono quantizzati e due \Delta^2 y su tre possono coincidere. Con \sigma \approx 0 la scala della perdita robusta diventa quasi nulla, ogni residuo è trattato come outlier, lo jacobiano pesato tende a zero e gli intervalli di confidenza esplodono (fit scartato come «non identificabile» pur essendo corretto). Sotto 8 campioni \sigma è «non stimabile»: sigma vale null, la scala della perdita è 1 e l'adeguatezza non è valutabile.

3. Rilevamento dei gradini (domini.trova_gradini)

Sulla serie_riferimento, con d_i = y_{i+1} - y_i e \sigma_d = 1.4826\,\mathrm{MAD}(d) calcolata su tutti i salti:

  1. Campioni in salita: gli indici con d_i > \theta_b\,\sigma_d (solo salite, A12). Il rumore bianco e la pendenza del raffreddamento restano sotto questa soglia.
  2. Raggruppamento: i campioni in salita a meno di distanza_min_salti_campioni indici l'uno dall'altro formano una salita (un gradino si distribuisce su 2–4 campioni, A6).
  3. Gradino: una salita è un gradino se la sua salita complessiva y_{\text{ultimo}+1} - y_{\text{primo}} supera \theta_g (A5). Si somma su tutta la salita perché ogni singolo salto può restare sotto soglia: un gradino di 38 °C diviso in 19 + 19 è un gradino, e uno scarto casuale (≤ 4 °C) no.
  4. Ogni gradino diventa la coppia (i_0, i_1): i_0 è il primo campione già in salita e i_1 il campione successivo all'ultimo in salita (estremo escluso). I campioni di salita sono quelli in [i_0, i_1).

4. Costruzione dei domini (domini.proponi)

Dati N gradini, in ordine di tempo:

  • riposo: dall'inizio al primo gradino (o tutto il segnale se non ci sono gradini);
  • per ogni gradino: un dominio riscaldamento sui suoi campioni di salita, seguito da un dominio raffreddamento fino al gradino successivo (o alla fine). Un raffreddamento senza campioni (gradini adiacenti) non viene creato.

Con N gradini i domini sono quindi 1 + 2N. I domini sono contigui e coprono tutti i campioni. L'utente può sostituire la proposta con domini.toml (A21), che viene validato: tipi ammessi, intervalli non vuoti, contiguità.

5. Fit (fit.adatta)

5.1 Modelli (modelli.MODELLI, t' = t - t_0)

Modello Equazione Parametri Dove si prova
costante y = a a riposo
lineare y = a + b\,t' a, b riposo
esp_decad y = c + A\,e^{-t'/\tau} c, A, \tau riposo, raffreddamento

Nel raffreddamento c è l'asintoto (≈ ambiente per T_C, ≈ 0 per dT_C), A l'ampiezza a t'=0 e \tau = 1/k la costante di tempo.

5.2 Partenza

  • Polinomi (costante, lineare): soluzione dei minimi quadrati lineari.
  • esp_decad: proiezione delle variabili. Per \tau fissato il modello è lineare in (c, A). Si provano n_griglia_tau valori di \tau su una griglia logaritmica tra \alpha_{\min} L e \alpha_{\max} L, per ognuno si risolvono (c, A) in modo lineare e si tiene il \tau con \mathrm{RSS} minima. Poi il raffinamento non lineare parte da lì: nessuna partenza casuale (A22).
  • Se c e \tau sono fissi e solo A è libero: A = \sum (y - c)\,e^{-t'/\tau} \big/ \sum e^{-2t'/\tau}.

5.3 Minimizzazione

scipy.optimize.least_squares sui soli parametri liberi, con x_scale="jac", limiti \tau \in [\alpha_{\min} L,\ \alpha_{\max} L] e perdita robusta (perdita, default soft_l1):

\min_p \sum_i \rho\!\left(\frac{r_i^2}{s^2}\right), \quad \rho(z) = 2\left(\sqrt{1+z}-1\right), \quad s = \kappa\,\sigma, \quad r_i = f(t'_i; p) - y_i

Per |r_i| \ll s si comporta come i minimi quadrati; per |r_i| \gg s cresce linearmente, così gli scarti casuali pesano molto meno (A14). Se \sigma non è finita o è nulla la scala è 1 (A20).

5.4 Incertezze

s^2 = \frac{\mathrm{RSS}}{n-k}, \qquad \mathrm{cov}(p) = s^2\,(J^\top J)^{-1}, \qquad \mathrm{IC}_{95\%} = \pm 1.96\,\sqrt{\mathrm{diag}(\mathrm{cov})}\cdot f_\rho

k è il numero di parametri liberi, J lo jacobiano dei residui nella soluzione (con la pseudo-inversa di J^\top J), \mathrm{RSS} = \sum r_i^2 (A16). Il fattore di autocorrelazione vale f_\rho = \sqrt{(1+\rho_1)/(1-\rho_1)} se correggi_autocorrelazione è attivo e \rho_1 > \theta_\rho, altrimenti 1 (A17). Per i parametri fissi l'incertezza è 0 (A18).

5.5 Metriche

Metrica Definizione
\mathrm{RMSE} \sqrt{\mathrm{RSS}/n}
R^2 1 - \mathrm{RSS}/\mathrm{TSS}, con \mathrm{TSS} = \sum (y_i - \bar y)^2
\hat\sigma_\text{res} scala robusta dei residui: 1.4826\,\mathrm{MAD}(r), con r i residui. Come la perdita del fit, ignora gli scarti casuali (la RMSE no)
\rho_1 autocorrelazione dei residui a lag 1: \sum_{i} \tilde r_{i+1}\tilde r_i \big/ \sum_i \tilde r_i^2, con \tilde r = r - \bar r
\mathrm{AICc} n\ln(\mathrm{RSS}/n) + 2k + \dfrac{2k(k+1)}{n-k-1} (infinito se n-k-1 \le 0)
\mathrm{BIC} n\ln(\mathrm{RSS}/n) + k\ln n
condizione numero di condizione di J con le colonne normalizzate in norma euclidea

AICc e BIC sono calcolati da \mathrm{RSS}, quindi senza peso robusto (A19).

5.6 Validità e identificabilità

Un fit è valido se converge, nessun parametro libero è sul limite (distanza dal limite inferiore a 10^{-3} volte il limite stesso: \tau al bordo della griglia), il numero di condizione non supera soglia_condizione e l'AICc è finito. Un fit è non identificabile se la condizione supera la soglia oppure un parametro libero diverso da a e c ha 1.96\,\sigma_p > |p| (intervallo che include lo zero: il parametro non è distinguibile da 0).

5b. Costanti condivise (analisi.analizza_domini)

Per A8 i raffreddamenti hanno lo stesso \tau e lo stesso c. L'ordine di elaborazione è quindi:

  1. il raffreddamento finale (ultimo dominio di tipo raffreddamento) si fitta per primo con c, A, \tau liberi;
  2. se la scelta è valida, i valori di c e \tau (parametri in raffreddamento_intermedio_fissa_da_finale) sono salvati per serie;
  3. ogni raffreddamento intermedio si fitta con quei valori fissi, e solo A libera (k = 1).

Motivo: un intermedio ha pochi campioni (da 4 a qualche decina) con scarti casuali, e \tau da solo non è identificabile (su dati puliti lo sarebbe: i test lo verificano). Con un solo parametro libero bastano pochi campioni: il minimo per questi fit è n_min_punti_parametri_fissi (default 4) invece di n_min_punti (6), così anche i tratti tra gradini ravvicinati (1–2 s) si fittano. Se il finale non dà un fit valido o la lista è vuota gli intermedi restano liberi.

6. Scelta del modello (analisi.scegli)

Per ogni (dominio, serie) con più modelli candidati:

  1. Scarto: via i fit non validi (§5.6). Se non ne resta nessuno: indeterminato.
  2. Adeguatezza: sono adeguati i fit con \hat\sigma_\text{res} \le \lambda\,\sigma (scala robusta dei residui, §5.5). Se ce n'è almeno uno si sceglie tra questi, altrimenti tra tutti i validi (e il dominio è segnato non adeguato). Se \sigma non è stimabile (meno di 8 campioni) l'adeguatezza non è valutabile (adeguato = null) e si sceglie tra tutti i validi.
  3. Parsimonia: tra i candidati con \mathrm{AICc} \le \mathrm{AICc}_{\min} + \delta vince quello con meno parametri (a parità, AICc minore).

Un dominio con meno di n_min_punti campioni (n_min_punti_parametri_fissi se c e \tau sono fissi) non viene fittato (indeterminato); un modello con più di n - 2 parametri liberi non viene provato. Con un solo modello candidato (raffreddamento) la scelta è lui, se valido: il flag adeguato resta informativo.

Perché la scala robusta e non la RMSE: gli scarti casuali del sensore gonfiano la RMSE (\mathrm{RMSE}/\sigma da 1.07 a 5.26 per il modello giusto) ma non la scala robusta (\hat\sigma_\text{res}/\sigma da 0.97 a 1.40 su 120 raffreddamenti finali). Con un modello sbagliato (lineare o costante sullo stesso tratto) la scala robusta vale da 11 a 100 volte \sigma, quindi discrimina bene. Con la RMSE solo il 27 % dei fit corretti risultava «adeguato»; con la scala robusta il 100 %.

7. Riscaldamento come salto (analisi._gradino, analisi._affina_salti)

Il dominio di riscaldamento ha 2–4 campioni: nessun fit. Per ogni serie si riporta il salto \Delta con due stime.

Misura sui campioni (delta_misurato), con i_0, i_1 gli estremi del dominio:

\Delta_\text{misurato} = y_{i_1-1} - y_{i_0-1}, \qquad \text{durata} = t_{i_1-1} - t_{i_0-1}

cioè dal campione prima della salita all'ultimo campione di salita. È una stima per difetto: l'ultimo campione di salita ha già cominciato a decadere (0.3–0.45 s dopo il gradino) e l'errore cresce con la temperatura, fino a ≈ 10 % (A23).

Stima dai modelli adiacenti (delta, usata quando disponibile). L'istante del gradino t_s è il punto medio tra l'ultimo campione prima della salita e il primo in salita. Con f_\text{prima} e f_\text{dopo} i modelli scelti per il dominio precedente (riposo o raffreddamento) e per il successivo (raffreddamento), valutati in t_s con t' = t - t_0 del rispettivo dominio:

\Delta = f_\text{dopo}(t_s) - f_\text{prima}(t_s)

Per due raffreddamenti, che condividono c e \tau, equivale a A_\text{dopo}\,e^{(t_{0,\text{dopo}}-t_s)/\tau} - A_\text{prima}\,e^{-(t_s-t_{0,\text{prima}})/\tau}. I campi da e a sono allora i valori dei due modelli in t_s.

Misura corretta (se manca un modello adiacente, p.es. un tratto indeterminato, ma c'è il raffreddamento finale). Si riportano i due campioni all'istante del gradino con la legge di decadimento del finale (c e \tau condivisi):

y_\text{prima}^* = c + (y_{i_0-1} - c)\,e^{-(t_s - t_{i_0-1})/\tau}, \qquad y_\text{dopo}^* = c + (y_{i_1-1} - c)\,e^{(t_{i_1-1} - t_s)/\tau}, \qquad \Delta = y_\text{dopo}^* - y_\text{prima}^*

Ripiego finale: la misura sui campioni, se non c'è nemmeno il raffreddamento finale (gradino senza fit dopo), o se la stima si discosta più del 50 % dalla misura (guardia contro fit anomali). Il campo stima vale "modelli", "misura corretta" o "misura"; il report mostra sempre anche la misura. sigma vale null e fit è vuoto.

Accuratezza (validazione su 60 seed di gen_data, test.md): con la stima dai modelli l'errore sull'ampiezza vera del gradino ha media +0.33 %, scarto quadratico medio 0.86 % e massimo 3.16 % su 124 salti (con scarta_dopo_salto_n = 0 nessun salto ha richiesto il ripiego). Con la sola misura l'errore medio era −5.5 % e arrivava a −20 %; la misura corretta, usata da sola, ha errore mediano 1.85 % e massimo 11.9 % (è un ripiego, non un sostituto dei modelli).

Esempio di calcolo completo (seed 7 di gen_data, T_C, raffreddamento finale)

  1. n = 620 campioni da t_0 = 26.79 s (scarta_dopo_salto_n = 0: tutti usati nel fit); \sigma = 0.115 °C.
  2. Griglia su \tau \in [0.05 L,\ 50 L] con L \approx 93 s → \tau_0 vicino a 12 s, c, A lineari.
  3. Raffinamento soft_l1 con scala 3\sigma = 0.35 °C: c = 23.92 \pm 0.06, A = 163.2 \pm 0.33, \tau = 12.12 \pm 0.04 s (valore vero \tau = 12.07 s).
  4. Metriche: \mathrm{RMSE} = 0.247, \hat\sigma_\text{res} = 0.159, \rho_1 = 0.61 (> 0.3: intervalli moltiplicati per \sqrt{1.61/0.39} = 2.03).
  5. \hat\sigma_\text{res} = 0.159 \le 1.5 \cdot 0.115 = 0.173: adeguato, anche se \mathrm{RMSE} = 0.247 > 1.5\,\sigma per gli scarti casuali (con il vecchio criterio sarebbe stato scartato).
  6. Gli intermedi 3 e 5 usano c = 23.92, \tau = 12.12 e stimano solo A (76.7 e 124.5).
  7. Il salto del terzo gradino: misurato +80.9, stimato dai modelli +90.2 (ampiezza vera +89.4).