"""
Gauss-Seidel AC Power Flow Solver (independent crosscheck benchmark)
"""

import numpy as np

def run_gauss_seidel(Ybus, bus_types, P_spec, Q_spec, V0=None,
                      tol=1e-6, max_iter=1000, acceleration=1.6):
    n = len(bus_types)
    V = np.ones(n, dtype=complex) if V0 is None else V0.copy()

    slack_idx = [i for i, t in enumerate(bus_types) if t == "slack"][0]
    pv_idx = [i for i, t in enumerate(bus_types) if t == "PV"]
    pq_idx = [i for i, t in enumerate(bus_types) if t == "PQ"]

    v_mag_set = np.abs(V).copy()

    history = []
    converged = False
    it = 0
    for it in range(1, max_iter + 1):
        V_prev = V.copy()
        for i in range(n):
            if i == slack_idx:
                continue
            if i in pv_idx:
                # recompute Q at PV bus from current voltage estimate
                Q_i = -np.imag(np.conj(V[i]) * (Ybus[i, :] @ V))
                S_i = complex(P_spec[i], Q_i)
            else:
                S_i = complex(P_spec[i], Q_spec[i])

            sum_yv = Ybus[i, :] @ V - Ybus[i, i] * V[i]
            V_new = (np.conj(S_i) / np.conj(V[i]) - sum_yv) / Ybus[i, i]
            V[i] = V[i] + acceleration * (V_new - V[i])

            if i in pv_idx:
                V[i] = v_mag_set[i] * V[i] / abs(V[i])

        max_change = np.max(np.abs(V - V_prev))
        history.append(max_change)
        if max_change < tol:
            converged = True
            break

    return {
        "converged": converged,
        "iterations": it,
        "V": np.abs(V),
        "theta": np.angle(V),
        "history": history,
    }
