"""
svd.py — numerical companion to the glossary entry
'singular value decomposition (SVD)'.

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 SVD A = V Lambda U^T of a rectangular matrix: the factors
         returned by np.linalg.svd reconstruct A to machine precision,
         V and U are orthonormal (V^T V = I, U^T U = I), and — in
         contrast to an EVD — the factorization exists for every matrix,
         including the rectangular 5 x 3 example and the defective
         matrix [[0, 1], [0, 0]] that admits no EVD.
[P-sing] Lambda is nonzero only on its main diagonal and its diagonal
         entries, the singular values, are nonnegative.

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

Data generated by pythondemos/svd.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 U^T with orthonormal V, U — any matrix")
A = rng.normal(size=(5, 3))                     # rectangular
Vf, s, Ut = np.linalg.svd(A)                    # full SVD
Lam = np.zeros((5, 3)); Lam[:3, :3] = np.diag(s)
check("reconstruction V Lambda U^T = A (err < 1e-12)",
      np.max(np.abs(Vf @ Lam @ Ut - A)) < 1e-12)
check("V orthonormal: V^T V = I", np.allclose(Vf.T @ Vf, np.eye(5)))
check("U orthonormal: U^T U = I", np.allclose(Ut @ Ut.T, np.eye(3)))
# exists even for the defective matrix that has no EVD
D = np.array([[0.0, 1.0], [0.0, 0.0]])
Vd, sd, Utd = np.linalg.svd(D)
check("defective [[0,1],[0,0]] (no EVD) still has an exact SVD",
      np.max(np.abs(Vd @ np.diag(sd) @ Utd - D)) < 1e-14)

# ---------------------------------------------------------- [P-sing]
print("[P-sing] Lambda diagonal, singular values nonnegative")
off_diag = Lam.copy(); np.fill_diagonal(off_diag, 0.0)
check("Lambda vanishes off the main diagonal", np.all(off_diag == 0))
check("singular values are nonnegative", np.all(s >= 0))
check("np.linalg.svd returns them sorted (descending)",
      np.all(np.diff(s) <= 0))

# ------------------------------------------------------------ preview
fig, ax = plt.subplots(1, 3, figsize=(9, 2.8))
for a, M, t in ((ax[0], A, "A (5 x 3)"), (ax[1], Lam, "Lambda"),
                (ax[2], Vf @ Lam @ Ut - A, "reconstruction error")):
    im = a.imshow(M, cmap="gray"); a.set_title(t)
    fig.colorbar(im, ax=a, shrink=0.75)
fig.suptitle("[P-def] SVD of a rectangular matrix")
fig.tight_layout()
fig.savefig("svd.png", dpi=110)
print(f"\n{sum(ok for _, ok in report)}/{len(report)} checks passed")
assert all(ok for _, ok in report)