<?php

namespace App\Support;

/**
 * Regresja wielomianowa - wierny port legacy polynomial/PolynomialRegression.php
 * (Andrew Que, GPL). Ten sam algorytm (równania normalne + eliminacja
 * wierszowa), arytmetyka natywna (float) zamiast bcmath - wynik identyczny dla
 * zakresów danych pomp, bez zależności od ext-bcmath.
 *
 * Użycie jak w legacy: new PolynomialRegression(4) → 4 współczynniki (wielomian
 * stopnia 3); addData(x,y) ×N; getCoefficients(); interpolate(coeffs, x).
 */
class PolynomialRegression
{
    /** @var float[] */
    private array $xPowers = [];
    /** @var float[] */
    private array $xyPowers = [];

    public function __construct(private int $degree = 3)
    {
        $this->reset();
    }

    public function reset(): void
    {
        $squares = ($this->degree - 1) * 2;
        for ($i = 0; $i <= $squares; $i++) {
            $this->xPowers[$i] = 0.0;
            $this->xyPowers[$i] = 0.0;
        }
    }

    public function addData(float $x, float $y): void
    {
        $squares = ($this->degree - 1) * 2;
        for ($i = 0; $i <= $squares; $i++) {
            $this->xPowers[$i] += $x ** $i;
            $this->xyPowers[$i] += $y * ($x ** $i);
        }
    }

    /**
     * Współczynniki [c0, c1, ...] (c0 = wyraz wolny).
     *
     * @return float[]
     */
    public function getCoefficients(?int $degree = null): array
    {
        $degree ??= $this->degree;
        $eps = 1e-12;

        // Macierz równań normalnych + kolumna augmentacji.
        $matrix = [];
        for ($row = 0; $row < $degree; $row++) {
            for ($col = 0; $col < $degree; $col++) {
                $matrix[$row][$col] = $this->xPowers[$row + $col];
            }
            $matrix[$row][$degree] = $this->xyPowers[$row];
        }

        $rows = count($matrix);
        $isDone = array_fill(0, $rows, false);
        $order = [];

        // Eliminacja wierszowa do macierzy górnotrójkątnej z 1 na przekątnej.
        for ($column = 0; $column < $rows; $column++) {
            $activeRow = 0;
            while ($activeRow < $rows
                && (abs($matrix[$activeRow][$column]) < $eps || $isDone[$activeRow])) {
                $activeRow++;
            }

            if ($activeRow < $rows) {
                $order[$column] = $activeRow;

                $firstTerm = $matrix[$activeRow][$column];
                for ($sub = $column; $sub <= $rows; $sub++) {
                    $matrix[$activeRow][$sub] /= $firstTerm;
                }
                $isDone[$activeRow] = true;

                for ($row = 0; $row < $rows; $row++) {
                    if (! $isDone[$row] && abs($matrix[$row][$column]) >= $eps) {
                        $factor = $matrix[$row][$column];
                        for ($sub = $column; $sub <= $rows; $sub++) {
                            $matrix[$row][$sub] -= $factor * $matrix[$activeRow][$sub];
                        }
                    }
                }
            }
        }

        // Podstawianie wsteczne.
        $isDone = array_fill(0, $rows, false);
        $coefficients = [];
        for ($column = $rows - 1; $column >= 0; $column--) {
            $activeRow = $order[$column];
            $isDone[$activeRow] = true;

            for ($row = 0; $row < $rows; $row++) {
                if (! $isDone[$row]) {
                    $factor = $matrix[$row][$column];
                    for ($sub = $column; $sub <= $rows; $sub++) {
                        $matrix[$row][$sub] -= $factor * $matrix[$activeRow][$sub];
                    }
                }
            }
            $coefficients[$column] = $matrix[$activeRow][$rows];
        }

        ksort($coefficients);

        return $coefficients;
    }

    /** @param float[] $coefficients */
    public function interpolate(array $coefficients, float $x): float
    {
        $y = 0.0;
        foreach ($coefficients as $i => $c) {
            $y += $c * ($x ** $i);
        }

        return $y;
    }
}