Titano/Titano/Core/LocalRegression.cs

namespace Titano.Core;

/// <summary>
/// Regressione lineare locale pesata, con seconda passata robusta agli scarti anomali.
///
/// È il nucleo numerico condiviso da tutto ciò che, nel motore, deve distinguere una
/// variazione reale e lenta da un disturbo veloce: l'esposizione che sfarfalla sopra una
/// rampa di tramonto, la deriva del bilanciamento del bianco, il tremolio del treppiede
/// sopra una panoramica voluta. In tutti questi casi il segnale utile è la componente
/// liscia e il disturbo è quella ad alta frequenza.
///
/// Rispetto a una media mobile, la componente lineare della regressione segue senza ritardo
/// le rampe: una media mobile su una salita costante restituisce sempre un valore vecchio di
/// mezza finestra, e la correzione risultante spegnerebbe proprio il fenomeno da preservare.
/// </summary>
public static class LocalRegression
{
    /// <summary>Costante del peso di Tukey: oltre 4.685 deviazioni il campione non pesa più nulla.</summary>
    private const double TukeyConstant = 4.685;

    /// <summary>
    /// Lisciatura della serie su una finestra mobile di <paramref name="windowFrames"/> campioni.
    /// Con <paramref name="rejectOutliers"/> si esegue una seconda passata in cui i campioni
    /// lontani dalla prima stima vengono neutralizzati.
    /// </summary>
    public static double[] Smooth(double[] series, int windowFrames, bool rejectOutliers)
    {
        int n = series.Length;
        var result = new double[n];
        if (n == 0) return result;
        if (n <= 2)
        {
            Array.Copy(series, result, n);
            return result;
        }

        int radius = Math.Clamp((Math.Max(3, windowFrames) - 1) / 2, 1, Math.Max(1, n - 1));
        double sigma = Math.Max(radius / 2.0, 0.5);

        var robust = new double[n];
        Array.Fill(robust, 1.0);

        FitAll(series, result, radius, sigma, robust);

        if (rejectOutliers)
        {
            UpdateRobustWeights(series, result, robust);
            FitAll(series, result, radius, sigma, robust);
        }

        return result;
    }

    private static void FitAll(double[] y, double[] output, int radius, double sigma, double[] robust)
    {
        int n = y.Length;
        double twoSigmaSq = 2.0 * sigma * sigma;

        for (int i = 0; i < n; i++)
        {
            int from = Math.Max(0, i - radius);
            int to = Math.Min(n - 1, i + radius);

            // Sistema normale della retta pesata y = a + b·t, con t = j - i.
            double sw = 0, swt = 0, swt2 = 0, swy = 0, swty = 0;
            for (int j = from; j <= to; j++)
            {
                double t = j - i;
                double w = Math.Exp(-(t * t) / twoSigmaSq) * robust[j];
                if (w <= 1e-9) continue;
                sw += w;
                swt += w * t;
                swt2 += w * t * t;
                swy += w * y[j];
                swty += w * t * y[j];
            }

            if (sw <= 1e-9)
            {
                output[i] = y[i];
                continue;
            }

            double det = sw * swt2 - swt * swt;
            if (Math.Abs(det) < 1e-12)
            {
                output[i] = swy / sw;                       // finestra degenere: media pesata
                continue;
            }

            double a = (swt2 * swy - swt * swty) / det;     // intercetta = valore stimato in t = 0
            output[i] = a;
        }
    }

    private static void UpdateRobustWeights(double[] y, double[] fit, double[] robust)
    {
        int n = y.Length;
        var residuals = new double[n];
        for (int i = 0; i < n; i++) residuals[i] = Math.Abs(y[i] - fit[i]);

        var sorted = (double[])residuals.Clone();
        Array.Sort(sorted);
        double mad = sorted[n / 2];
        double scale = 1.4826 * mad;

        // Se la sequenza è già pulita non si scarta nulla: evita di amplificare il rumore numerico.
        if (scale < 1e-4)
        {
            Array.Fill(robust, 1.0);
            return;
        }

        for (int i = 0; i < n; i++)
        {
            double u = residuals[i] / (TukeyConstant * scale);
            robust[i] = u >= 1.0 ? 0.0 : Math.Pow(1.0 - u * u, 2.0);
        }
    }
}