import sympy as sp
from dataclasses import dataclass, field
from typing import Dict, Tuple

@dataclass
class StochastischesModellApriori:

    sigma_obs: Iterable[float]           # σ_i
    group_ids: Iterable[int]            # Gruppenzugehörigkeit der i-ten Beobachtung
    sigma0_sq_groups: Dict[int, float] = field(default_factory=dict)

    def __post_init__(self):
        # In sympy-Objekte konvertieren
        self.sigma_obs = sp.Matrix(list(self.sigma_obs))  # Spaltenvektor
        self.group_ids = sp.Matrix(list(self.group_ids))  # Spaltenvektor

        if self.sigma_obs.rows != self.group_ids.rows:
            raise ValueError("sigma_obs und group_ids müssen gleich viele Einträge haben.")

        # Fehlende Gruppen mit σ_0j^2 = 1.0 initialisieren
        unique_groups = sorted({int(g) for g in self.group_ids})
        for g in unique_groups:
            if g not in self.sigma0_sq_groups:
                self.sigma0_sq_groups[g] = 1.0

    @property

    def n_obs(self) -> int:
        return int(self.sigma_obs.rows)


    def build_Qll_P(self) -> Tuple[sp.Matrix, sp.Matrix]:

        n = self.n_obs
        Q_ll = sp.zeros(n, n)
        P = sp.zeros(n, n)

        for i in range(n):
            sigma_i = self.sigma_obs[i, 0]
            g = int(self.group_ids[i, 0])
            sigma0_sq = self.sigma0_sq_groups[g]

            q_ii = sigma_i**2
            Q_ll[i, i] = q_ii

            P[i, i] = 1 / (sigma0_sq * q_ii)

        return Q_ll, P

    @staticmethod
    def _redundanz_pro_beobachtung(A: sp.Matrix, P: sp.Matrix) -> sp.Matrix:

        n_obs = P.rows
        n_param = A.cols

        # P^(1/2) aufbauen (diagonal, sqrt der Diagonale)
        sqrtP = sp.zeros(n_obs, n_obs)
        for i in range(n_obs):
            sqrtP[i, i] = sp.sqrt(P[i, i])

        A_tilde = sqrtP * A  # Ã

        # M = (Ãᵀ Ã)^(-1)
        M = (A_tilde.T * A_tilde).inv()

        r_vec = sp.zeros(n_obs, 1)

        for i in range(n_obs):
            a_i = A_tilde.row(i)          # 1 × n_param
            a_i_row = sp.Matrix([a_i])    # explizit 1×n-Matrix
            r_i = 1 - (a_i_row * M * a_i_row.T)[0, 0]
            r_vec[i, 0] = r_i

        return r_vec


    def varianzkomponenten_schaetzung(
        self,
        v: sp.Matrix,   # Residuenvektor (n × 1)
        A: sp.Matrix,   # Designmatrix
    ) -> Dict[int, float]:

        if v.rows != self.n_obs:
            raise ValueError("Länge von v passt nicht zur Anzahl Beobachtungen im Modell.")

        # Aktuelle Gewichte
        Q_ll, P = self.build_Qll_P()

        # Redundanzzahlen pro Beobachtung
        r_vec = self._redundanz_pro_beobachtung(A, P)

        new_sigma0_sq: Dict[int, float] = {}

        # Für jede Gruppe j:
        unique_groups = sorted({int(g) for g in self.group_ids})

        for g in unique_groups:
            # Indizes der Beobachtungen in dieser Gruppe
            idx = [i for i in range(self.n_obs) if int(self.group_ids[i, 0]) == g]
            if not idx:
                continue

            # v_j, P_j, r_j extrahieren
            v_j = sp.Matrix([v[i, 0] for i in idx])               # (m_j × 1)
            P_j = sp.zeros(len(idx), len(idx))
            r_j = 0
            for ii, i in enumerate(idx):
                P_j[ii, ii] = P[i, i]
                r_j += r_vec[i, 0]

            # σ̂_j^2 = (v_jᵀ P_j v_j) / r_j
            sigma_hat_j_sq = (v_j.T * P_j * v_j)[0, 0] / r_j

            # als float rausgeben, kann man aber auch symbolisch lassen
            new_sigma0_sq[g] = float(sigma_hat_j_sq)

        return new_sigma0_sq

    def update_sigma0(self, new_sigma0_sq: Dict[int, float]) -> None:

        for g, val in new_sigma0_sq.items():
            self.sigma0_sq_groups[int(g)] = float(val)