import numpy as np import matplotlib.pyplot as plt

1. 1 for GOE, 2 for GUE, 4 for GSE

betas = 1, 2, 4 Ns = [1, 2, 4, 8, 16, 32] # matrix sizes

1. Choose number of samples

Nmatr = 10000 repeats = 10

1. The following condition selects the desired ensemble: a number Nmatr
2. of matrices are diagonalized, and the eigenvalues are collected in the vector E

Es = {} for _ in range(repeats):

   for N in Ns:
     for beta in betas:
       if beta == 1:  # Gaussian Orthogonal Ensemble
           M = np.random.randn(Nmatr, N, N)
           M = (M + M.transpose((0, 2, 1))) / 2
           E = np.linalg.eigvals(M.reshape(Nmatr, N, N)).flatten()
       elif beta == 2:  # Gaussian Unitary Ensemble
           M_real = np.random.randn(Nmatr, N, N)
           M_imag = np.random.randn(Nmatr, N, N)
           M = (M_real + 1j * M_imag + M_real.transpose((0, 2, 1)) - 1j * M_imag.transpose((0, 2, 1))) / 2
           E = np.linalg.eigvals(M.reshape(Nmatr, N, N)).flatten()
       elif beta == 4:  # Gaussian Symplectic Ensemble
           A = np.random.randn(Nmatr, N, N) + 1j * np.random.randn(Nmatr, N, N)
           B = np.random.randn(Nmatr, N, N) + 1j * np.random.randn(Nmatr, N, N)

           M_top = np.block(A, B)
           M_bottom = np.block(-np.conj(B), np.conj(A))
           M = np.block([[M_top], [M_bottom]])

           M = (M + np.conj(M.transpose((0, 2, 1)))) / 2
           E = np.linalg.eigvals(M.reshape(Nmatr, 2 * N, 2 * N)).flatten()
       
       if (N, beta) in Es: 
           Es[(N, beta)]= np.append(Es[(N, beta)], E)
       else: 
           Es[(N, beta)] = E

fig, axs = plt.subplots(2, 3, figsize=(18, 9))

legends = {1: "GOE", 2: "GUE", 4: "GSE"} colors = {1: "blue", 2: "red", 4: "green"}

for i, N in enumerate(Ns):

   row = i // 3
   col = i % 3
   ax = axs[row, col]

   for beta in betas:
       color = colors[beta]
       E = Es[(N, beta)]
       xs = np.real(E) / np.sqrt(2 * beta * N)

       bin_heights, bin_borders, _ = ax.hist(xs, bins=500, density=True, color=color, alpha=0.1)
       bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2

       # Compute sliding window average
       window_size = 5
       window = np.ones(window_size) / window_size
       smoothed_heights = np.convolve(bin_heights, window, mode='same')

       # Plot sliding window average
       ax.plot(bin_centers, smoothed_heights, label=legends[beta], color=color)

   # Add plot labels and title
   ax.set_xlabel('x', fontsize=14)
   ax.set_ylabel('ρ(x)', fontsize=14)
   ax.grid(True)
   ax.legend()

plt.tight_layout() fig.suptitle(r'Eigenvalues $/\sqrtTemplate:2N\beta$, with N = {} to {}'.format(Ns[0], Ns[-1]), fontsize=18, y=1.04) plt.show()