r"""
PT-Symmetry Breaking in Pseudo-Hermitian Ensembles
==================================================

Bender and Boettcher observed that certain non-Hermitian Hamiltonians
with PT (parity-time) symmetry can nonetheless have entirely real
spectra. This ensemble realizes the equivalent condition of
pseudo-Hermiticity, :math:`P H P^{-1} = H^\dagger` for a fixed
Hermitian involution :math:`P = \mathrm{diag}(I_p, -I_q)`
(:math:`n=p+q`). Writing :math:`H` in the matching :math:`2\times 2`
block form

.. math::

    H = \begin{pmatrix} A & B \\ -B^\dagger & D \end{pmatrix},

pseudo-Hermiticity forces :math:`A = A^\dagger` (:math:`p \times p`),
:math:`D = D^\dagger` (:math:`q \times q`) -- independent GOE/GUE-type
Hermitian blocks -- while the off-diagonal (coupling) block :math:`B`
is an independent Ginibre-type :math:`p \times q` block scaled by the
non-Hermiticity strength :math:`g`. At :math:`g=0`, :math:`H` is
block-diagonal and trivially has a real spectrum (the "unbroken"
PT-symmetric phase); since :math:`H` is always similar to
:math:`H^\dagger`, its eigenvalues remain closed under complex
conjugation for any :math:`g`, either real or in exact
complex-conjugate pairs. As :math:`g` grows, an increasing fraction of
eigenvalues collide (an "exceptional point", where eigenvectors also
coalesce) and split off the real axis as conjugate pairs -- the
PT-symmetry-breaking transition.

This example reproduces the PT-symmetry-breaking transition: the
fraction of real eigenvalues of a pseudo-Hermitian ensemble drops from
1 (fully real, "unbroken" PT symmetry) toward 0 as the
non-Hermiticity/coupling strength g increases.

Reference: C. M. Bender, S. Boettcher, Phys. Rev. Lett. 80 (1998) 5243.

Run:
    python examples/paper_replications/pt_symmetric_demo.py
"""

import matplotlib.pyplot as plt
import numpy as np

import physicskit.rmt as rmt

P = Q = 20
N_SAMPLES = 40
SEED = 2026

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))

g_values = np.linspace(0.0, 3.0, 25)
for beta, color in [(1, "steelblue"), (2, "indianred")]:
    fractions = []
    for g in g_values:
        ensemble = rmt.ensembles.PTSymmetricEnsemble(p=P, q=Q, g=g, beta=beta, seed=SEED)
        spectrum = ensemble.sample(n_samples=N_SAMPLES)
        fractions.append(rmt.stats.real_eigenvalue_fraction(spectrum.eigenvalues).mean())
    axes[0].plot(g_values, fractions, "o-", ms=3, color=color, label=f"beta={beta}")

axes[0].set_xlabel("coupling strength g")
axes[0].set_ylabel("mean fraction of real eigenvalues")
axes[0].set_title("PT-symmetry-breaking transition")
axes[0].legend()

# Eigenvalue scatter at a representative g, showing the real-axis
# survivors and the complex-conjugate pairs that have already split off.
ensemble = rmt.ensembles.PTSymmetricEnsemble(p=P, q=Q, g=1.0, beta=2, seed=SEED)
spectrum = ensemble.sample(n_samples=1)
eigs = spectrum.eigenvalues[0]
axes[1].scatter(eigs.real, eigs.imag, s=20, color="indianred")
axes[1].axhline(0.0, color="k", lw=1)
axes[1].set_xlabel("Re")
axes[1].set_ylabel("Im")
axes[1].set_title("Eigenvalues at g=1.0 (beta=2)\nreal survivors + conjugate pairs")

# --- Panel 3: phase diagram over (g, system size) ---
# Panel 1 fixes n=p+q=40 and only varies g; the transition itself shifts
# with system size (a larger ensemble has more real-eigenvalue pairs
# that can collide), so a 2D map of the same real_eigenvalue_fraction
# statistic over BOTH g and n resolves that size-dependence directly,
# rather than a single 1D slice through it.
g_grid_2d = np.concatenate([[0.0], np.geomspace(0.02, 3.0, 15)])
n_half_values = [5, 10, 20, 40]
fraction_grid = np.empty((len(n_half_values), len(g_grid_2d)))
for row, half in enumerate(n_half_values):
    for col, g in enumerate(g_grid_2d):
        # beta=1 (GOE-type blocks): the slower of the two transitions in
        # panel 1, giving better resolution across this g range than
        # beta=2's much sharper drop.
        ensemble = rmt.ensembles.PTSymmetricEnsemble(p=half, q=half, g=float(g), beta=1, seed=SEED)
        spectrum = ensemble.sample(n_samples=N_SAMPLES)
        fraction_grid[row, col] = rmt.stats.real_eigenvalue_fraction(spectrum.eigenvalues).mean()

ax = axes[2]
im = ax.imshow(fraction_grid, aspect="auto", origin="lower", cmap="viridis", vmin=0.0, vmax=1.0)
ax.set_yticks(range(len(n_half_values)))
ax.set_yticklabels([f"{2 * h}" for h in n_half_values])
tick_idx = [0, 3, 7, 11, 15]
ax.set_xticks(tick_idx)
ax.set_xticklabels([f"{g_grid_2d[i]:.2g}" for i in tick_idx])
ax.set_xlabel("coupling strength g (log-spaced grid)")
ax.set_ylabel("n = p+q")
ax.set_title("Phase diagram (beta=1):\nmean real-eigenvalue fraction")
fig.colorbar(im, ax=ax, label="mean fraction real")

fig.suptitle(f"PT-symmetric ensemble -- p=q={P}")
fig.tight_layout()
out_path = "pt_symmetric_replication.png"
fig.savefig(out_path, dpi=150)
print(f"Saved {out_path}")