"""
evd.py — numerical companion to the glossary entry
'eigenvalue decomposition (EVD)'.

One block per paragraph of the entry (marked [P...]): each block verifies
numerically what the corresponding statement asserts. Self-contained
(numpy/matplotlib only), fixed seed.

Blocks
------
[P-def]  The EVD A = V Lambda V^{-1}: reconstructing A from the factors
         returned by np.linalg.eig gives A back to machine precision; the
         columns of V are eigenvectors and Lambda is diagonal with the
         matching eigenvalues.
[P-diag] Matrices that admit an EVD are the diagonalizable ones: for the
         defective matrix [[0, 1], [0, 0]] the eigenvector matrix is
         singular (rank 1), so V^{-1} does not exist and no EVD is
         possible, while the diagonalizable matrix A above passes.

Outputs
-------
evd.png : preview figure (checking only).

Data generated by pythondemos/evd.py.
"""

import numpy as np
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
report = []


def check(name, ok):
    report.append((name, bool(ok)))
    print(f"  [{'ok' if ok else 'FAIL'}] {name}")


# ----------------------------------------------------------- [P-def]
print("[P-def] A = V Lambda V^{-1} reconstructs A")
A = np.array([[2.0, 1.0, 0.0], [1.0, 3.0, 0.5], [0.0, 0.5, 1.5]])
lam, V = np.linalg.eig(A)
Lam = np.diag(lam)
A_rec = V @ Lam @ np.linalg.inv(V)
check("reconstruction error < 1e-12",
      np.max(np.abs(A_rec - A)) < 1e-12)
check("columns of V are eigenvectors (A v = lambda v)",
      all(np.linalg.norm(A @ V[:, i] - lam[i] * V[:, i]) < 1e-12
          for i in range(3)))
check("Lambda is diagonal", np.allclose(Lam, np.diag(np.diag(Lam))))

# ---------------------------------------------------------- [P-diag]
print("[P-diag] EVD exists iff the matrix is diagonalizable")
D = np.array([[0.0, 1.0], [0.0, 0.0]])         # defective (not diagonalizable)
lamD, VD = np.linalg.eig(D)
rank_VD = np.linalg.matrix_rank(VD)
check("defective matrix: eigenvector matrix is singular (rank 1)",
      rank_VD == 1)
sym = A                                         # symmetric example above
check("diagonalizable matrix: eigenvector matrix invertible (rank 3)",
      np.linalg.matrix_rank(V) == 3)

# ------------------------------------------------------------ preview
fig, ax = plt.subplots(1, 3, figsize=(9, 2.8))
for a, M, t in ((ax[0], A, "A"), (ax[1], Lam.real, "Lambda"),
                (ax[2], (A_rec - A).real, "V Lambda V^{-1} - A")):
    im = a.imshow(M, cmap="gray"); a.set_title(t)
    fig.colorbar(im, ax=a, shrink=0.75)
fig.suptitle("[P-def] EVD factors and reconstruction error")
fig.tight_layout()
fig.savefig("evd.png", dpi=110)
print(f"\n{sum(ok for _, ok in report)}/{len(report)} checks passed")
assert all(ok for _, ok in report)