"""optimize_phi_sym.py -- multi-start optimization of Phi_sym over biseparable pure
states (1|3 and 2|2 cuts) and over ALL pure 4-qubit states (unconstrained).

Key finding: all three searches converge to the SAME value, 6/sqrt(7) -- i.e. Phi_sym's
biseparable supremum equals its global supremum over the entire state space, achieved
both by connected 4-qubit graph states AND by a trivial biseparable state (two Bell
pairs). This shows Phi_sym cannot certify genuine multipartite entanglement at that
threshold.
"""
import numpy as np
from scipy.optimize import minimize
from core import state_1_3, state_2_2
from core2 import phi_sym


def neg_phi_1_3(params, source):
    psi = state_1_3(params, source)
    val, _ = phi_sym(psi)
    return -val


def neg_phi_2_2(params, pair1, pair2):
    psi = state_2_2(params, pair1, pair2)
    val, _ = phi_sym(psi)
    return -val


def neg_phi_general(params):
    v = params[:16] + 1j * params[16:]
    v = v / np.linalg.norm(v)
    val, _ = phi_sym(v)
    return -val


def run_multistart(objective, nparams, args, n_restarts, seed0, label):
    best = -np.inf
    bx = None
    for i in range(n_restarts):
        rng = np.random.default_rng(seed0 + i)
        x0 = rng.normal(size=nparams)
        res = minimize(objective, x0, args=args, method='Powell',
                        options={'maxiter': 2000, 'xtol': 1e-9, 'ftol': 1e-11})
        v = -res.fun
        if v > best:
            best = v
            bx = res.x
    print(f'{label}: best Phi_sym = {best:.10f}   ({n_restarts} restarts)')
    return best, bx


if __name__ == "__main__":
    # biseparable across 2|2 cut AB|CD -> exact optimum 6/sqrt(7), achieved by
    # Bell_AB (x) Bell_CD (verified in closed form: two maximal Bell pairs)
    v22, x22 = run_multistart(neg_phi_2_2, 16, ((0, 1), (2, 3)), 15, 100,
                               "biseparable 2|2 (AB|CD)")
    # biseparable across 1|3 cut A|BCD -> exact optimum (1+18/sqrt(7))/4,
    # achieved by (any single qubit) (x) GHZ_3(B,C,D)
    v13, x13 = run_multistart(neg_phi_1_3, 18, (0,), 15, 200,
                               "biseparable 1|3 (A|BCD)")
    # fully unconstrained over ALL pure 4-qubit states -> same ceiling 6/sqrt(7)
    vgen, xgen = run_multistart(neg_phi_general, 32, (), 20, 500,
                                 "unconstrained (all 4-qubit states)")

    print()
    print("6/sqrt(7)          =", 6 / np.sqrt(7))
    print("(1+18/sqrt(7))/4   =", (1 + 18 / np.sqrt(7)) / 4)