"""Check of e:grass-gue (Section 3, p:grass-gue):

    E det(w - H)/det(z - H)
      = (N/2pi) int int e^{-N(x^2+y^2)/2} (w-iy)^{N-1} ((z-x)(w-iy) - 1)/(z-x)^{N+1} dx dy

for the GUE with E|H_ij|^2 = 1/N.  For N = 1, 2, 3 the left side is computed by Gauss-Hermite
quadrature over the eigenvalue density (prop. to Vandermonde^2 exp(-N sum lambda^2/2));
for N = 4 by Monte Carlo.  The case w = z, where both sides equal one, is included.
"""
import numpy as np
from numpy.polynomial.hermite import hermgauss


def direct(N, z, w, n=60):
    t, wt = hermgauss(n)
    lam = t * np.sqrt(2.0 / N)                  # e^{-N lam^2/2} = e^{-t^2}
    grids = np.meshgrid(*([lam] * N), indexing='ij')
    W = np.ones_like(grids[0])
    for g in np.meshgrid(*([wt] * N), indexing='ij'):
        W = W * g
    V = np.ones_like(grids[0])
    for i in range(N):
        for j in range(i + 1, N):
            V = V * (grids[i] - grids[j]) ** 2
    f = np.ones_like(grids[0], dtype=complex)
    for i in range(N):
        f = f * (w - grids[i]) / (z - grids[i])
    return (W * V * f).sum() / (W * V).sum()


def susy(N, z, w, n=200):
    t, wt = hermgauss(n)
    s = t * np.sqrt(2.0 / N)
    X, Y = np.meshgrid(s, s, indexing='ij')
    WW = np.outer(wt, wt) * (2.0 / N)           # dx dy = (2/N) dt dt
    integrand = (w - 1j * Y) ** (N - 1) / (z - X) ** (N + 1) * ((z - X) * (w - 1j * Y) - 1)
    return N / (2 * np.pi) * (WW * integrand).sum()


if __name__ == '__main__':
    for (z, w) in [(0.3 + 0.7j, -0.2 + 0.1j), (-0.5 + 0.9j, 1.1 - 0.4j), (0.1 + 0.8j, 0.1 + 0.8j)]:
        for N in [1, 2, 3]:
            a, b = direct(N, z, w), susy(N, z, w)
            print(f"N={N} z={z} w={w}  density: {a:.10f}  formula: {b:.10f}  |diff| = {abs(a - b):.1e}")
    rng = np.random.default_rng(1)
    N, z, w, M = 4, 0.3 + 0.7j, -0.2 + 0.1j, 200000
    acc = 0
    for _ in range(M // 1000):
        A = (rng.normal(size=(1000, N, N)) + 1j * rng.normal(size=(1000, N, N))) / np.sqrt(2)
        H = (A + np.conj(np.transpose(A, (0, 2, 1)))) / np.sqrt(2 * N)   # E|H_ij|^2 = 1/N
        ev = np.linalg.eigvalsh(H)
        acc += np.prod((w - ev) / (z - ev), axis=1).sum()
    print(f"N=4 Monte Carlo ({M} samples): {acc / M:.4f}  formula: {susy(N, z, w):.4f}")