hexagon.py

"""
NETWORK SURFACE: none. Pure math, fully local.
SPDX-License-Identifier: LicenseRef-Sovereign-MIT-1.0
See LICENSE for terms.

Matematyka heksagonu 6 cnót na sferze Blocha.
Q = pure-state pairwise fidelity (inspired by SWSSB framework,
    Lessa et al. PRX Quantum 6, 010344, 2025; arXiv:2405.03639).
    qsim implementuje pure-state proxy, nie mixed-state SWSSB order parameter.
Berry phase: Berry, Proc. R. Soc. A 392, 45 (1984);
    Wilczek-Zee, Phys. Rev. Lett. 52, 2111 (1984);
    discrete formulation: Resta, Phys. Rev. Lett. 80, 1800 (1998).
"""

from math import sqrt, radians, cos, acos
from cmath import exp as cexp
import numpy as np

HEXAGON = [
    "odwaga",        # 0°
    "pokora",        # 60°
    "przebaczenie",  # 120°
    "współczucie",   # 180°
    "wdzięczność",   # 240°
    "rozumienie",    # 300°
]

VIRTUE_COLORS = {
    "odwaga":       "#DC143C",
    "pokora":       "#6A0DAD",
    "przebaczenie": "#FFD700",
    "współczucie":  "#00CED1",
    "wdzięczność":  "#C0C0C0",
    "rozumienie":   "#2ECC71",
}

Q_STOPS = [
    (0.00, "#4B0082"),
    (0.40, "#1E90FF"),
    (0.60, "#00CED1"),
    (0.75, "#FFB800"),
    (0.84, "#FF6B00"),
    (0.90, "#FFFFFF"),
]

Q_ATTRACTOR = 0.83929
_x = (1 + sqrt(4 * Q_ATTRACTOR - 3)) / 2   # 0.79882
THETA_IDEAL = 2 * acos(sqrt(_x))             # ~0.932 rad


def _hex_color(t: float, r1: str, r2: str) -> str:
    """Interpolacja liniowa między dwoma hex kolorami, t ∈ [0,1]."""
    def parse(h): return tuple(int(h[i:i+2], 16) for i in (1, 3, 5))
    a, b = parse(r1), parse(r2)
    rgb = tuple(round(a[i] + (b[i] - a[i]) * t) for i in range(3))
    return "#{:02X}{:02X}{:02X}".format(*rgb)


def q_to_chroma(q: float) -> str:
    """Q ∈ [0,1] → hex kolor z Q_STOPS."""
    q = max(0.0, min(1.0, q))
    for i in range(len(Q_STOPS) - 1):
        q0, c0 = Q_STOPS[i]
        q1, c1 = Q_STOPS[i + 1]
        if q <= q1:
            t = (q - q0) / (q1 - q0) if q1 > q0 else 0.0
            return _hex_color(t, c0, c1)
    return Q_STOPS[-1][1]


def virtue_state(phi_deg: float, strength: float = 1.0) -> np.ndarray:
    """Kubit na sferze Blocha. strength ∈ [0,1]."""
    phi = radians(phi_deg)
    theta = THETA_IDEAL * float(np.clip(strength, 0, 1))
    return np.array([
        cos(theta / 2),
        cexp(1j * phi) * sin(theta / 2)
    ], dtype=complex)


def _fidelity(s1: np.ndarray, s2: np.ndarray) -> float:
    return float(abs(np.dot(s1.conj(), s2)) ** 2)


def compute_q(strengths: list[float]) -> dict:
    """
    strengths: lista 6 floatów ∈ [0,1] w kolejności HEXAGON.

    Zwraca dict:
      q        — fidelity correlator ∈ [0,1]
      gamma    — Berry phase [rad]
      aniso    — std(fidelities)
      q_renyi  — mean(fids^2)
      dominant — HEXAGON[argmin(fids)]  (najsłabszy sprzęg)
      fids     — list[float] len=6
      chroma   — hex kolor Q
    """
    states = [
        virtue_state(i * 60, s)
        for i, s in enumerate(strengths)
    ]
    # sąsiednie pary (cyklicznie)
    fids = [
        _fidelity(states[i], states[(i + 1) % 6])
        for i in range(6)
    ]
    q = float(np.mean(fids))

    # Berry phase — iloczyn amplitud sąsiednich kubitów
    berry = np.prod([np.dot(states[i].conj(), states[(i + 1) % 6])
                     for i in range(6)])
    gamma = float(np.angle(berry))

    aniso = float(np.std(fids))
    q_renyi = float(np.mean(np.array(fids) ** 2))
    dominant = HEXAGON[int(np.argmin(fids))]
    chroma = q_to_chroma(q)

    return dict(q=q, gamma=gamma, aniso=aniso, q_renyi=q_renyi,
                dominant=dominant, fids=fids, chroma=chroma)


def is_crystallized(q_result: dict) -> bool:
    """HARD: aniso < 0.001 AND |Q - Q_ATTRACTOR| < 0.002"""
    return (q_result["aniso"] < 0.001 and
            abs(q_result["q"] - Q_ATTRACTOR) < 0.002)


def anneal(init_strengths: list[float],
           T_start: float = 1.0,
           T_end: float = 0.001,
           n_steps: int = 2000,
           sigma: float = 0.05) -> list[float]:
    """
    Simulated annealing → s_crystal bliskie Q_ATTRACTOR.
    Zweryfikowano: dist<0.001 w ~150ms.
    """
    s = np.array(init_strengths, dtype=float)
    qr = compute_q(s.tolist())

    # cold start skip
    if qr["aniso"] < 0.001 and abs(qr["q"] - Q_ATTRACTOR) < 0.001:
        return s.tolist()

    def cost(sv):
        r = compute_q(sv.tolist())
        return abs(r["q"] - Q_ATTRACTOR) + r["aniso"]

    current_cost = cost(s)
    best_s = s.copy()
    best_cost = current_cost

    rng = np.random.default_rng(42)
    for step in range(n_steps):
        T = T_start * (T_end / T_start) ** (step / n_steps)
        delta = rng.normal(0, sigma, size=6)
        s_new = np.clip(s + delta, 0.0, 1.0)
        new_cost = cost(s_new)
        diff = new_cost - current_cost
        if diff < 0 or rng.random() < np.exp(-diff / T):
            s = s_new
            current_cost = new_cost
            if current_cost < best_cost:
                best_s = s.copy()
                best_cost = current_cost

    return best_s.tolist()


# math.sin nie importuje się razem z from math — uzupełnienie
from math import sin  # noqa: E402