Hands-on lab — Day 106: Eigenvalues and Eigenvectors, Intuitively

"""Exercise 1 — most vectors get knocked off their line. A few do not.

Run from inside examples/:

    ../.venv/bin/python3 01_the_fan_of_vectors.py

No symbols, no equation, no theory. Take twenty-four directions spread evenly
around the circle, apply the matrix to each one, and measure how far each
output has swung away from its input. Then look at the column of numbers and
notice that two of them are zero.
"""

from __future__ import annotations

import numpy as np

from dataset import A, A_EIGEN_ANGLES_DEG
from eigen import deviation_degrees, direction_degrees, eigen_lines_by_sweep, sweep_deviations

SCRIPT = "01_the_fan_of_vectors.py"


def main() -> None:
    print(f"{SCRIPT}")
    print("=" * 72)
    print()
    print("The matrix, which is just a transformation that moves the grid:")
    print()
    print("    A = [[4, 1],")
    print("         [2, 3]]")
    print()

    # ---------------------------------------------------------------- 1
    print("1. Twenty-four directions, one every 15 degrees, each one a unit vector.")
    print("   For each: where it points, where its output points, and how far")
    print("   apart those two directions are.")
    print()
    print("   in (deg)   out (deg)   swung by (deg)   length in -> out")
    print("   " + "-" * 60)

    zero_deviation = []
    for angle in range(0, 360, 15):
        radians = np.radians(angle)
        v = np.array([np.cos(radians), np.sin(radians)])
        out = A @ v
        swing = deviation_degrees(v, out)
        length = float(np.linalg.norm(out))
        marker = "   <-- kept its direction" if swing < 1e-9 else ""
        print(
            f"   {angle:6d}     {direction_degrees(out):7.3f}     {swing:9.4f}"
            f"        1.000 -> {length:5.3f}{marker}"
        )
        if swing < 1e-9:
            zero_deviation.append(angle)

    print()
    print(f"   Directions that came back on their own line: {zero_deviation}")
    print("   45 and 225 are the same line, pointing opposite ways along it.")
    print("   So out of twenty-four directions, ONE line survived.")
    print()
    assert zero_deviation == [45, 225], zero_deviation

    # ---------------------------------------------------------------- 2
    print("2. But look again at the column. It dips towards zero twice, not once.")
    print("   Around 45 degrees it reaches zero exactly. Around 120 degrees it")
    print("   gets down to 5.36 and climbs again — so the true minimum is")
    print("   somewhere between the sample points, and the coarse fan stepped")
    print("   right over it.")
    print()
    print("   Sweeping every thousandth of a degree from 0 to 180 instead:")
    print()

    found = eigen_lines_by_sweep(A)
    for angle in found["lines"]:
        deviation = sweep_deviations(A, [angle])[0][0]
        print(f"     a surviving line near {angle:11.6f} degrees  (deviation {deviation:.3e})")
    print()
    print(f"   Directions swept: 180,000. Directions that kept their line to")
    print(f"   within a hundredth of a degree: {found['fraction'] * 180000:.0f}, in {len(found['lines'])} separate bands.")
    print()
    print("   Two lines, not one. Only 180 degrees was swept because a line and")
    print("   its reverse are the same line, so the other half is a repeat.")
    print()
    assert found["verdict"] == "some"
    assert len(found["lines"]) == 2, found
    assert np.allclose(found["lines"], [45.0, 116.565], atol=1e-2), found

    # ---------------------------------------------------------------- 3
    print("3. Those two lines have names you can write down exactly.")
    print()
    for angle in A_EIGEN_ANGLES_DEG:
        radians = np.radians(angle)
        v = np.array([np.cos(radians), np.sin(radians)])
        out = A @ v
        stretch = float(np.linalg.norm(out)) / float(np.linalg.norm(v))
        print(f"     direction {angle:18.14f} degrees")
        print(f"       in  = [{v[0]: .6f}, {v[1]: .6f}]")
        print(f"       out = [{out[0]: .6f}, {out[1]: .6f}]")
        print(f"       swung by {deviation_degrees(v, out):.3e} degrees, stretched by {stretch:.6f}")
        print()
        assert deviation_degrees(v, out) < 1e-9

    print("   45 degrees is the direction of (1, 1). The output is 5 times longer.")
    print("   116.565... degrees is the direction of (1, -2). The output is 2 times longer.")
    print()
    print("   Check both on paper, with no decimals at all:")
    print()
    for vector in ((1.0, 1.0), (1.0, -2.0)):
        v = np.array(vector)
        out = A @ v
        factor = out[0] / v[0]
        print(f"     A @ ({v[0]:.0f}, {v[1]:.0f}) = ({out[0]:.0f}, {out[1]:.0f}) = {factor:.0f} * ({v[0]:.0f}, {v[1]:.0f})")
        assert np.allclose(out, factor * v)
    print()
    print("   That is the entire definition, arrived at by measurement:")
    print("   a vector that the matrix only STRETCHES is an eigenvector, and")
    print("   the stretch factor is its eigenvalue.")
    print()

    # ---------------------------------------------------------------- 4
    print("4. Why the zero vector is excluded, even though it fits the equation.")
    print()
    zero = np.array([0.0, 0.0])
    print(f"     A @ (0, 0) = ({(A @ zero)[0]:.0f}, {(A @ zero)[1]:.0f})")
    print("     and that equals lambda * (0, 0) for lambda = 5, for lambda = 2,")
    print("     for lambda = 1000, and for every other number there is.")
    print("     It satisfies the equation for EVERY value, so it distinguishes")
    print("     nothing, which is another way of saying it tells you nothing.")
    print("     It also has no direction to keep. So it is ruled out by")
    print("     definition, and that exclusion is doing real work rather than")
    print("     being bookkeeping.")
    print()
    assert np.allclose(A @ zero, 0.0)
    for lam in (5.0, 2.0, 1000.0, -3.5):
        assert np.allclose(A @ zero, lam * zero)

    print(f"{SCRIPT}: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 2 — solving a 2x2 with a pencil, then checking against NumPy.

Run from inside examples/:

    ../.venv/bin/python3 02_by_hand_2x2.py

Exercise 1 found the eigendirections by brute-force measurement, which works
and does not scale. This script does the same job with algebra: derive the
characteristic equation, solve the quadratic, and read the eigenvectors out of
the squashed matrix. Then hand the same matrix to numpy.linalg.eig and check.

The check is where the interesting part is, and it is not the part anyone
expects.
"""

from __future__ import annotations

import warnings

import numpy as np

from dataset import A, A_EIGENVALUES, A_EIGENVECTORS
from eigen import abs_cosine, characteristic_coefficients, eigenvalues_2x2, eigenvector_2x2

SCRIPT = "02_by_hand_2x2.py"

TOL = 1e-9


def main() -> None:
    print(f"{SCRIPT}")
    print("=" * 72)
    print()

    # ---------------------------------------------------------------- 1
    print("1. The characteristic equation, derived rather than quoted.")
    print()
    print("   Ask: for which numbers lambda does A v = lambda v have a solution")
    print("   other than v = 0? Rearrange it so one side is zero:")
    print()
    print("       A v = lambda v")
    print("       A v - lambda v = 0")
    print("       (A - lambda I) v = 0")
    print()
    print("   The I is there because you cannot subtract a number from a matrix.")
    print("   lambda I is the number lambda spread down the diagonal, which does")
    print("   to v exactly what multiplying by lambda does.")
    print()
    print("   Now (A - lambda I) is one matrix, and it sends some non-zero v to")
    print("   the origin. Day 102 named the matrices that do that: the ones with")
    print("   determinant zero, the ones that squash the plane onto a line. So:")
    print()
    print("       det(A - lambda I) = 0")
    print()

    trace, determinant = characteristic_coefficients(A)
    print(f"   For A = [[4, 1], [2, 3]]:  trace = {trace:.0f}, determinant = {determinant:.0f}")
    print()
    print("       A - lambda I = [[4 - lambda,          1],")
    print("                       [         2, 3 - lambda]]")
    print()
    print("       det = (4 - lambda)(3 - lambda) - 1*2")
    print("           = 12 - 4*lambda - 3*lambda + lambda^2 - 2")
    print("           = lambda^2 - 7*lambda + 10")
    print()
    print("   and the 7 is the trace while the 10 is the determinant, which is")
    print("   not a coincidence and holds for every 2x2.")
    print()
    assert trace == 7.0 and determinant == 10.0

    # ---------------------------------------------------------------- 2
    print("2. Solve the quadratic.")
    print()
    discriminant = trace * trace - 4.0 * determinant
    print(f"       discriminant = trace^2 - 4*det = {trace:.0f}^2 - 4*{determinant:.0f} = {discriminant:.0f}")
    print(f"       sqrt({discriminant:.0f}) = {np.sqrt(discriminant):.0f}")
    print(f"       lambda = (7 +/- 3) / 2  ->  5 and 2")
    print()
    print("   Or factorise it and skip the formula:  (lambda - 5)(lambda - 2) = 0")
    print()

    by_hand = eigenvalues_2x2(A)
    print(f"   eigenvalues_2x2(A) = {tuple(round(value.real, 12) for value in by_hand)}")
    print(f"   imaginary parts    = {tuple(round(value.imag, 12) for value in by_hand)}  (both zero: two real eigenvalues)")
    print()
    assert abs(by_hand[0].real - 5.0) < TOL
    assert abs(by_hand[1].real - 2.0) < TOL
    assert all(abs(value.imag) < TOL for value in by_hand)

    # ---------------------------------------------------------------- 3
    print("3. Read each eigenvector out of the squashed matrix.")
    print()
    for eigenvalue, expected in zip(A_EIGENVALUES, A_EIGENVECTORS):
        shifted = A - eigenvalue * np.eye(2)
        v = eigenvector_2x2(A, eigenvalue)
        print(f"   lambda = {eigenvalue:.0f}:")
        print(f"       A - {eigenvalue:.0f}I = [[{shifted[0, 0]:5.0f}, {shifted[0, 1]:5.0f}],")
        print(f"{'':17}[{shifted[1, 0]:5.0f}, {shifted[1, 1]:5.0f}]]")
        print(f"       determinant of that = {np.linalg.det(shifted):.1e}  (zero, as it must be)")
        print(f"       row 0 says {shifted[0, 0]:.0f}*x + {shifted[0, 1]:.0f}*y = 0")
        print(f"       one solution: ({expected[0]:.0f}, {expected[1]:.0f})")
        print(f"       eigenvector_2x2 returned: [{v[0]: .6f}, {v[1]: .6f}]")
        print(f"       same line? abs_cosine = {abs_cosine(v, expected):.15f}")
        if float(np.dot(v, expected)) < 0.0:
            print("       (it came back pointing the other way along that line.")
            print("        Not a bug. The row [p, q] gives (-q, p) and the sign of")
            print("        that depends on which row happened to be non-zero.)")
        out = A @ np.array(expected)
        print(f"       check: A @ ({expected[0]:.0f}, {expected[1]:.0f}) = ({out[0]:.0f}, {out[1]:.0f}) = {eigenvalue:.0f} * ({expected[0]:.0f}, {expected[1]:.0f})")
        print()
        assert abs(np.linalg.det(shifted)) < 1e-12
        assert abs_cosine(v, expected) > 1.0 - TOL
        assert np.allclose(out, eigenvalue * np.array(expected))

    print("   Notice what is NOT determined here. Row 0 of (A - 5I) says")
    print("   -x + y = 0. That is one equation for two unknowns, so it names a")
    print("   whole LINE of solutions: (1, 1), (2, 2), (-1, -1), (0.3, 0.3).")
    print("   All of them are eigenvectors. There is no such thing as THE")
    print("   eigenvector for an eigenvalue, only the eigen-LINE, and every")
    print("   library that returns one has made an arbitrary choice on your")
    print("   behalf. Section 5 shows NumPy making exactly that choice.")
    print()

    # ---------------------------------------------------------------- 4
    print("4. Now numpy.linalg.eig on the same matrix.")
    print()
    values, vectors = np.linalg.eig(A)
    print(f"   eigenvalues  = {values}")
    print(f"   dtype        = {values.dtype}")
    print(f"   eigenvectors =")
    for row in vectors:
        print(f"                  [{row[0]!s:>26}  {row[1]!s:>26}]")
    print(f"   dtype        = {vectors.dtype}")
    print()
    print("   The eigenvalues are 5 and 2, matching the hand calculation. But")
    print("   they came back as complex128 with zero imaginary parts, on a real")
    print("   matrix with two real eigenvalues.")
    print()
    print("   That is worth pausing on, because the docstring shipped with this")
    print("   very version of NumPy says otherwise. Quoting it exactly:")
    print()
    print('       "The resulting array will be of complex type, unless the')
    print('        imaginary part is zero in which case it will be cast to a')
    print('        real type."')
    print()
    print(f"   Observed on numpy {np.__version__} on the authoring machine: the")
    print("   imaginary part IS zero, and the array was NOT cast to a real type.")
    print("   Every case tried came back complex128 — including numpy.eye(2),")
    print("   whose eigenvalues are both exactly 1.")
    print()
    for name, matrix in (
        ("numpy.eye(2)", np.eye(2)),
        ("numpy.diag([1., 2., 3.])", np.diag([1.0, 2.0, 3.0])),
        ("an integer matrix [[2,0],[0,3]]", np.array([[2, 0], [0, 3]])),
    ):
        dtype = np.linalg.eig(matrix)[0].dtype
        print(f"       {name:<32} -> {dtype}")
        assert dtype == np.complex128
    print()
    assert values.dtype == np.complex128
    assert vectors.dtype == np.complex128
    assert np.all(values.imag == 0.0)

    # ---------------------------------------------------------------- 5
    print("5. Why that costs you something, and what to do about it.")
    print()
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        as_float = values.astype(float)
        categories = [item.category.__name__ for item in caught]
    print(f"   values.astype(float) = {as_float}   and it warned: {categories}")
    print()
    print("   The fix is one attribute, and it should be a reflex:")
    print()
    real_values = values.real
    real_vectors = vectors.real
    print(f"       values.real  = {real_values}   dtype {real_values.dtype}")
    print()
    print("   But .real is a claim, so check it before you make it:")
    print()
    print("       if numpy.all(numpy.abs(values.imag) < 1e-12):")
    print("           values = values.real")
    print("       else:")
    print("           ... this matrix has no real eigenvalues, handle that case")
    print()
    print("   Taking .real without the check on a rotation matrix silently")
    print("   throws away the entire answer. Exercise 3 shows that happening.")
    print()
    assert as_float.dtype == np.float64
    assert categories == ["ComplexWarning"]
    assert real_values.dtype == np.float64

    # ---------------------------------------------------------------- 6
    print("6. Do the hand answer and the NumPy answer agree?")
    print()
    order = np.argsort(real_values)[::-1]
    print("   eigenvalues, largest first:")
    for hand, numpy_value in zip(A_EIGENVALUES, real_values[order]):
        print(f"       by hand {hand:.0f}   numpy {numpy_value:.12f}   difference {abs(hand - numpy_value):.3e}")
        assert abs(hand - numpy_value) < TOL
    print()
    print("   eigenvectors — and here the naive comparison FAILS:")
    print()
    for eigenvalue, hand, column in zip(A_EIGENVALUES, A_EIGENVECTORS, order):
        theirs = real_vectors[:, column]
        mine = np.array(hand) / np.linalg.norm(hand)
        print(f"       lambda = {eigenvalue:.0f}")
        print(f"         mine   [{mine[0]: .6f}, {mine[1]: .6f}]")
        print(f"         numpy  [{theirs[0]: .6f}, {theirs[1]: .6f}]")
        print(f"         numpy.allclose(mine, theirs)      -> {np.allclose(mine, theirs)}")
        print(f"         abs_cosine(mine, theirs)          -> {abs_cosine(mine, theirs):.15f}")
        print()
        assert abs_cosine(mine, theirs) > 1.0 - TOL

    print("   For lambda = 2 the two answers are exact negatives of each other,")
    print("   and numpy.allclose says False. Nothing is wrong. (1, -2) and")
    print("   (-1, 2) name the same line, and the equation A v = lambda v holds")
    print("   for both — multiply both sides by -1 and it is the same statement.")
    print()
    print("   So: never compare eigenvectors component by component. Compare")
    print("   the absolute cosine, which asks the only question that has a")
    print("   determinate answer: do these two lie on the same LINE?")
    print()

    # ---------------------------------------------------------------- 7
    print("7. The equation itself, checked numerically for every returned pair.")
    print()
    for index in range(len(real_values)):
        lam = real_values[index]
        v = real_vectors[:, index]
        left = A @ v
        right = lam * v
        residual = float(np.linalg.norm(left - right))
        print(f"       lambda = {lam:.6f}")
        print(f"         A @ v      = [{left[0]: .12f}, {left[1]: .12f}]")
        print(f"         lambda * v = [{right[0]: .12f}, {right[1]: .12f}]")
        print(f"         residual   = {residual:.3e}   (tolerance 1e-12)")
        assert residual < 1e-12
    print()
    print("   That residual check is the one to keep. It does not care about")
    print("   sign, scale, ordering or dtype: it asks whether the matrix really")
    print("   does to v what multiplying by a single number does.")
    print()

    print(f"{SCRIPT}: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 3 — the eigenvectors of the transformations you already know.

Run from inside examples/:

    ../.venv/bin/python3 03_standard_transformations.py

Day 102 built a small vocabulary of transformations: scale, reflect, shear,
rotate, project. Each one is a matrix, so each one has eigenvalues, and the
answers are not all of the same kind. One has every direction. One has two.
One has exactly one. One has none at all.

The "none at all" case is the one to spend time on, because it is the moment
the geometry and the algebra agree loudly.
"""

from __future__ import annotations

import numpy as np

from dataset import (
    ROTATION_60,
    ROTATION_90,
    SHEAR,
    STANDARD_TRANSFORMATIONS,
    SYMMETRIC,
    SYMMETRIC_3X3,
)
from eigen import abs_cosine, direction_degrees, eigen_lines_by_sweep, sweep_deviations

SCRIPT = "03_standard_transformations.py"

TOL = 1e-9


def describe(matrix) -> tuple[np.ndarray, np.ndarray, bool]:
    """eigenvalues, eigenvectors, and whether the eigenvalues are all real."""
    values, vectors = np.linalg.eig(np.asarray(matrix, dtype=float))
    is_real = bool(np.all(np.abs(values.imag) < 1e-12))
    return values, vectors, is_real


def main() -> None:
    print(f"{SCRIPT}")
    print("=" * 72)
    print()

    # ---------------------------------------------------------------- 1
    print("1. Every transformation from Day 102, and what it does to directions.")
    print()
    print("   Two independent answers per matrix: what numpy.linalg.eig says,")
    print("   and what a brute-force sweep of 180,000 directions measures. They")
    print("   have to agree, and where they seem not to, the disagreement is")
    print("   itself the lesson.")
    print()
    for name, (matrix, note) in STANDARD_TRANSFORMATIONS.items():
        values, vectors, is_real = describe(matrix)
        determinant = float(np.linalg.det(matrix))
        found = eigen_lines_by_sweep(matrix)
        print(f"   {name}")
        print(f"     matrix       [[{matrix[0, 0]: .4f}, {matrix[0, 1]: .4f}], [{matrix[1, 0]: .4f}, {matrix[1, 1]: .4f}]]")
        print(f"     determinant  {determinant: .6f}")
        if is_real:
            real_values = values.real
            real_vectors = vectors.real
            print(f"     eigenvalues  {np.array2string(real_values, precision=6, suppress_small=True)}   (real)")
            directions = sorted({round(direction_degrees(real_vectors[:, i]), 6) for i in range(2)})
            same_line = abs_cosine(real_vectors[:, 0], real_vectors[:, 1]) > 1.0 - 1e-8
            print(f"     eig columns  {directions} degrees" + ("  <-- BOTH THE SAME LINE" if same_line else ""))
        else:
            print(f"     eigenvalues  {np.array2string(values, precision=6, suppress_small=True)}   (COMPLEX: no real eigenvector)")
        summary = found["verdict"] if found["verdict"] != "some" else f"{len(found['lines'])} line(s) near {found['lines']} degrees"
        print(f"     measured     {summary}")
        if found["collapsed"]:
            print(f"     collapsed    {found['collapsed']} degrees sent to the origin (eigenvalue 0)")
        print(f"     {note}")
        print()

    print("   Three rows in that table are worth arguing with.")
    print()
    print("   The identity and the uniform scaling: eig returned the columns 0")
    print("   and 90 degrees, which reads like two special directions. The sweep")
    print("   says 'every direction', and the sweep is right. When every")
    print("   direction is an eigenvector, eig still has to return exactly two")
    print("   columns, so it returns a basis and the arbitrariness is invisible.")
    print()
    print("   The projection: the sweep found ONE surviving line and separately")
    print("   reported that 90 degrees was collapsed to the origin. The y-axis")
    print("   IS an eigenvector, with eigenvalue 0, but 'did it keep its")
    print("   direction?' cannot be answered about a vector that no longer has")
    print("   one. Measuring angles cannot see eigenvalue 0; the algebra can.")
    print()
    print("   The shear: eig returned two columns, the sweep found one line.")
    print("   Section 2 is about that.")
    print()

    # ---------------------------------------------------------------- 2
    print("2. The shear: one eigenvalue, repeated, and only one eigen-line.")
    print()
    values, vectors, _ = describe(SHEAR)
    print("   S = [[1, 1],")
    print("        [0, 1]]  — slides the top of the unit square to the right.")
    print()
    print(f"   eigenvalues: {np.array2string(values.real, precision=6, suppress_small=True)}  — 1 twice. Its ALGEBRAIC multiplicity is 2.")
    print(f"   NumPy returned two eigenvectors anyway:")
    for i in range(2):
        column = vectors.real[:, i]
        print(f"       column {i}: [{column[0]: .17f}, {column[1]: .17f}]")
    similarity = abs_cosine(vectors.real[:, 0], vectors.real[:, 1])
    print()
    print(f"   abs_cosine between them = {similarity:.15f}")
    print("   They are the same line. NumPy has to return a square array of")
    print("   eigenvectors, so when there are not enough distinct directions to")
    print("   fill it, it fills the space anyway. Counting COLUMNS would tell you")
    print("   there are two eigendirections. There is one. The number of")
    print("   independent directions is the GEOMETRIC multiplicity, and here it")
    print("   is 1 while the algebraic multiplicity is 2 — the gap is exactly")
    print("   what makes this matrix impossible to diagonalise.")
    print()
    print("   Geometrically it is obvious. A shear leaves the x-axis alone and")
    print("   tilts everything else. Only the horizontal line survives.")
    print()
    assert np.allclose(values.real, [1.0, 1.0], atol=TOL)
    assert similarity > 1.0 - 1e-8

    # Confirm it by brute force too: sweep the circle and count the bands.
    found = eigen_lines_by_sweep(SHEAR)
    print(f"   Brute-force check: sweeping 180,000 directions finds")
    print(f"   {len(found['lines'])} surviving line, near {found['lines'][0]} degrees — the x-axis, to")
    print(f"   the accuracy a sampled sweep can offer. One line, found by")
    print(f"   measurement, agreeing with the one line found by algebra.")
    print()
    assert found["verdict"] == "some"
    assert len(found["lines"]) == 1
    assert min(found["lines"][0], 180.0 - found["lines"][0]) < 0.05

    # ---------------------------------------------------------------- 3
    print("3. The rotation: no real eigenvector at all, and you can SEE why.")
    print()
    for name, matrix, degrees in (("90 degrees", ROTATION_90, 90.0), ("60 degrees", ROTATION_60, 60.0)):
        values, _, is_real = describe(matrix)
        print(f"   rotation by {name}:")
        print(f"       eigenvalues  {np.array2string(values, precision=6, suppress_small=True)}")
        print(f"       all real?    {is_real}")
        print(f"       magnitudes   {np.round(np.abs(values), 12).tolist()}  (rotation changes no lengths, so both are 1)")
        print(f"       real parts   {np.round(values.real, 6).tolist()}  = cos({degrees:.0f} degrees) = {np.cos(np.radians(degrees)):.6f}")
        print(f"       imag parts   {np.round(values.imag, 6).tolist()}  = +/- sin({degrees:.0f} degrees) = {np.sin(np.radians(degrees)):.6f}")
        print()
        assert not is_real
        assert np.allclose(np.abs(values), 1.0, atol=TOL)

    print("   The geometry says it first: a rotation turns EVERY vector by the")
    print("   same angle. If a vector kept its line it would have to be turned")
    print("   by 0 or by 180 degrees, and this rotation turns by neither. So no")
    print("   direction survives, and there is nothing for the algebra to find.")
    print()
    print("   The algebra says the same thing in its own words. For the 90-degree")
    print("   rotation, trace = 0 and determinant = 1, so:")
    print("       lambda^2 - 0*lambda + 1 = 0   ->   lambda^2 = -1")
    print("   and no real number squares to -1. The negative discriminant is not")
    print("   an error message; it is the algebra reporting the geometry.")
    print()
    trace = float(ROTATION_90[0, 0] + ROTATION_90[1, 1])
    determinant = float(np.linalg.det(ROTATION_90))
    print(f"       measured: trace = {trace:.0f}, determinant = {determinant:.0f}, discriminant = {trace * trace - 4 * determinant:.0f}")
    print()
    assert trace == 0.0 and abs(determinant - 1.0) < TOL

    grid = np.arange(0.0, 180.0, 0.001)
    deviations, _collapsed = sweep_deviations(ROTATION_90, grid)
    smallest = float(np.nanmin(deviations))
    found = eigen_lines_by_sweep(ROTATION_90)
    print(f"   Brute-force check: over 180,000 directions the SMALLEST swing was")
    print(f"   {smallest:.6f} degrees, and the sweep reports verdict '{found['verdict']}'.")
    print(f"   Nothing came close to keeping its line. A 90-degree rotation turns")
    print(f"   every single direction by 90 degrees, which is the largest swing")
    print(f"   there is once you measure lines rather than arrows.")
    print()
    assert smallest > 89.0
    assert found["verdict"] == "none" and found["lines"] == []

    print("   And here is the trap from exercise 2, sprung. Take .real without")
    print("   checking, and the answer becomes:")
    print()
    values, _, _ = describe(ROTATION_90)
    print(f"       values        = {values}")
    print(f"       values.real   = {values.real}   <-- both zero. The answer is gone.")
    print("   Two eigenvalues of magnitude 1 have been silently reported as 0.")
    print()
    assert np.allclose(values.real, 0.0, atol=TOL)

    # ---------------------------------------------------------------- 4
    print("4. Eigenvalue 0 means the matrix collapsed a direction, and that is")
    print("   the same fact as determinant 0 from Day 102.")
    print()
    projection = STANDARD_TRANSFORMATIONS["projection onto x-axis"][0]
    values, vectors, _ = describe(projection)
    real_values = values.real
    print("   P = [[1, 0],")
    print("        [0, 0]]  — flattens the whole plane onto the x-axis.")
    print()
    print(f"       eigenvalues  {np.array2string(real_values, precision=6, suppress_small=True)}")
    print(f"       determinant  {np.linalg.det(projection):.6f}")
    print(f"       product of the eigenvalues = {float(np.prod(real_values)):.6f}")
    print(f"       sum of the eigenvalues     = {float(np.sum(real_values)):.6f},  trace = {float(np.trace(projection)):.6f}")
    print()
    print("   Two identities worth carrying, both checkable on every matrix in")
    print("   section 1: the eigenvalues multiply to the determinant and add to")
    print("   the trace. So a zero eigenvalue and a zero determinant are the same")
    print("   news arriving twice — some direction was squashed to nothing, and")
    print("   the transformation cannot be undone.")
    print()
    for name, (matrix, _note) in STANDARD_TRANSFORMATIONS.items():
        values, _, _ = describe(matrix)
        product = complex(np.prod(values))
        total = complex(np.sum(values))
        assert abs(product.real - np.linalg.det(matrix)) < 1e-9, name
        assert abs(total.real - np.trace(matrix)) < 1e-9, name
        print(f"       {name:<24} product {product.real: .6f} = det,   sum {total.real: .6f} = trace")
    print()

    # ---------------------------------------------------------------- 5
    print("5. Symmetric matrices: two guarantees, checked rather than proved.")
    print()
    print("   A symmetric matrix — one equal to its own transpose — always has")
    print("   REAL eigenvalues, and eigenvectors at RIGHT ANGLES to each other.")
    print("   The proof is standard and is in any linear algebra text; this lab")
    print("   only checks that the guarantee holds on the matrices it has.")
    print()
    for name, matrix in (("SYMMETRIC (2x2)", SYMMETRIC), ("SYMMETRIC_3X3", SYMMETRIC_3X3)):
        assert np.allclose(matrix, matrix.T)
        values, vectors = np.linalg.eigh(matrix)
        print(f"   {name}: symmetric? {np.array_equal(matrix, matrix.T)}")
        print(f"       numpy.linalg.eigh eigenvalues {np.array2string(values, precision=6)}   dtype {values.dtype}")
        gram = vectors.T @ vectors
        off_diagonal = float(np.abs(gram - np.eye(len(values))).max())
        print(f"       eigenvectors mutually perpendicular: largest deviation from")
        print(f"       the identity in V.T @ V is {off_diagonal:.3e}")
        for i in range(len(values)):
            residual = float(np.linalg.norm(matrix @ vectors[:, i] - values[i] * vectors[:, i]))
            assert residual < 1e-12
        print(f"       every A v = lambda v residual below 1e-12: True")
        print()
        assert values.dtype == np.float64
        assert off_diagonal < 1e-12

    print("   Note the dtype. numpy.linalg.eigh returned float64 — no complex,")
    print("   no .real needed — and its values came back sorted ascending. eig")
    print("   guarantees neither. When your matrix is symmetric, and a covariance")
    print("   matrix always is, eigh is the right call.")
    print()

    print(f"{SCRIPT}: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 4 — apply the matrix over and over and watch a direction win.

Run from inside examples/:

    ../.venv/bin/python3 04_power_method.py

Multiply a random vector by A. Then multiply the answer by A. Then again.
Nothing here is clever, and after twenty-five rounds you have found the
dominant eigenvector to ten decimal places.

That algorithm has a name — the power method — and it is not a toy. It is
where PageRank came from, it is what an implicitly restarted Arnoldi
iteration is a sophisticated version of, and it is how you find the leading
eigenvector of a matrix far too large to factorise, because it never needs the
matrix itself, only the ability to multiply by it.
"""

from __future__ import annotations

import numpy as np

from dataset import (
    A,
    A_EIGEN_ANGLES_DEG,
    A_EIGENVALUES,
    A_EIGENVECTORS,
    SYMMETRIC,
    power_method_start,
)
from eigen import abs_cosine, direction_degrees, power_method, rayleigh_quotient

SCRIPT = "04_power_method.py"

TOL = 1e-9


def main() -> None:
    print(f"{SCRIPT}")
    print("=" * 72)
    print()

    start = power_method_start()
    dominant = np.array(A_EIGENVECTORS[0])

    # ---------------------------------------------------------------- 1
    print("1. Why repeated multiplication should work at all.")
    print()
    print("   Write the starting vector as a mixture of the two eigenvectors:")
    print()
    basis = np.column_stack([np.array(A_EIGENVECTORS[0]), np.array(A_EIGENVECTORS[1])])
    weights = np.linalg.solve(basis, start)
    print(f"       v0 = [{start[0]: .6f}, {start[1]: .6f}]")
    print(f"          = {weights[0]: .6f} * (1, 1)  +  {weights[1]: .6f} * (1, -2)")
    reconstructed = weights[0] * basis[:, 0] + weights[1] * basis[:, 1]
    print(f"       check: that mixture is [{reconstructed[0]: .6f}, {reconstructed[1]: .6f}]")
    print()
    assert np.allclose(reconstructed, start, atol=1e-12)

    print("   Now apply A. It does not mix the ingredients: it multiplies each")
    print("   one by its own eigenvalue, because that is what an eigenvector IS.")
    print()
    print("       A^k v0 = c1 * 5^k * (1, 1)  +  c2 * 2^k * (1, -2)")
    print()
    print("   So the (1, 1) ingredient grows by 5 each round and the (1, -2)")
    print("   ingredient grows by 2. After k rounds the second is smaller than")
    print("   the first by a factor of (2/5)^k, which goes to nothing:")
    print()
    for k in (1, 5, 10, 20, 25):
        print(f"       k = {k:2d}:  (2/5)^k = {(0.4 ** k):.3e}")
    print()
    print("   Nothing needed to be eliminated. The dominant direction simply")
    print("   outgrew the other one, and that is the whole idea.")
    print()

    # ---------------------------------------------------------------- 2
    print("2. Watch it happen, one iteration at a time.")
    print()
    print("   The starting vector is drawn from a seeded random generator, on")
    print("   purpose: the claim is that ALMOST ANY start converges, so a start")
    print("   chosen to work would prove nothing.")
    print()
    print(f"       v0 = [{start[0]: .6f}, {start[1]: .6f}], pointing at {direction_degrees(start):.6f} degrees")
    print(f"       target: the eigen-line at {direction_degrees(dominant):.6f} degrees, eigenvalue {A_EIGENVALUES[0]:.0f}")
    print()
    print("    k    direction (deg)   Rayleigh quotient   off target (deg)")
    print("   " + "-" * 62)
    v = start.copy()
    for k in range(0, 13):
        quotient = rayleigh_quotient(A, v)
        off = float(np.degrees(np.arccos(min(1.0, abs_cosine(v, dominant)))))
        print(f"   {k:2d}     {direction_degrees(v):11.6f}      {quotient:12.8f}      {off:10.6f}")
        w = A @ v
        v = w / np.linalg.norm(w)
    print()
    print("   The direction crawls in from 20 degrees and settles on 45, and the")
    print("   Rayleigh quotient — (v . A v) / (v . v), the best single number to")
    print("   call the eigenvalue for a given v — settles on 5 alongside it.")
    print()
    print("   The textbook claim about the Rayleigh quotient is that its error is")
    print("   the SQUARE of the vector's, so it converges twice as fast. Measure")
    print("   it rather than repeat it, because on this matrix it is not true.")
    print()
    print("      k   angle error (rad)   quotient error   ratio to angle   to angle^2")
    print("   " + "-" * 72)
    v = start.copy()
    for k in range(1, 9):
        w = A @ v
        v = w / np.linalg.norm(w)
        angle = float(np.arccos(min(1.0, abs_cosine(v, dominant))))
        error = abs(rayleigh_quotient(A, v) - 5.0)
        print(f"   {k:3d}      {angle:.6e}     {error:.6e}     {error / angle:10.6f}   {error / angle**2:10.2f}")
    print()
    print(f"   The ratio to the ANGLE settles on 1.0. The ratio to the angle")
    print(f"   SQUARED runs away. On A the Rayleigh quotient converges linearly,")
    print(f"   at exactly the same rate as the vector, and buys nothing.")
    print()
    print("   Now the same measurement on the symmetric matrix [[2, 1], [1, 2]],")
    print("   whose dominant eigenvalue is 3:")
    print()
    print("      k   angle error (rad)   quotient error   ratio to angle   to angle^2")
    print("   " + "-" * 72)
    v = start.copy()
    symmetric_target = np.array([1.0, 1.0])
    for k in range(1, 9):
        w = SYMMETRIC @ v
        v = w / np.linalg.norm(w)
        angle = float(np.arccos(min(1.0, abs_cosine(v, symmetric_target))))
        error = abs(rayleigh_quotient(SYMMETRIC, v) - 3.0)
        print(f"   {k:3d}      {angle:.6e}     {error:.6e}     {error / angle:10.6f}   {error / angle**2:10.4f}")
    print()
    print("   There the ratio to the angle squared locks onto 2.0000 and stays,")
    print("   which is textbook quadratic convergence.")
    print()
    print("   So the textbook claim is right, and it has a condition attached")
    print("   that is easy to drop: the quadratic result needs the eigenvectors")
    print("   to be at right angles, which SYMMETRY guarantees and A does not")
    print("   have. A is not symmetric — 1 in one corner and 2 in the other —")
    print("   so its eigen-lines at 45 and 116.6 degrees meet at 71.6 degrees,")
    print("   not 90, and the speed-up does not apply.")
    print()
    angle_between = A_EIGEN_ANGLES_DEG[1] - A_EIGEN_ANGLES_DEG[0]
    print(f"       angle between A's two eigen-lines: {angle_between:.4f} degrees")
    print(f"       angle between the symmetric matrix's:  90.0000 degrees")
    print()
    assert abs(angle_between - 71.565051) < 1e-4

    # ---------------------------------------------------------------- 3
    print("3. Run it to convergence and report the count.")
    print()
    result = power_method(A, start, tol=1e-10)
    print(f"       tolerance          1e-10 on the distance between successive unit vectors")
    print(f"       iterations         {result['iterations']}")
    print(f"       converged          {result['converged']}")
    print(f"       final change       {result['change']:.6e}")
    print(f"       vector             [{result['vector'][0]: .12f}, {result['vector'][1]: .12f}]")
    print(f"       direction          {direction_degrees(result['vector']):.12f} degrees")
    print(f"       eigenvalue         {result['eigenvalue']:.12f}")
    print(f"       abs_cosine with (1, 1)  {abs_cosine(result['vector'], dominant):.15f}")
    print()
    assert result["converged"]
    assert result["iterations"] == 25
    assert abs_cosine(result["vector"], dominant) > 1.0 - 1e-12
    assert abs(result["eigenvalue"] - 5.0) < 1e-9

    # ---------------------------------------------------------------- 4
    print("4. The convergence RATE is the eigenvalue ratio, and you can measure it.")
    print()
    print("   Theory says the error should shrink by |lambda2 / lambda1| = 2/5 =")
    print("   0.4 each round. Divide each step's change by the previous one:")
    print()
    print("      k    change        ratio to previous")
    print("   " + "-" * 44)
    history = result["history"]
    for index in range(4, 14):
        ratio = history[index] / history[index - 1]
        print(f"     {index + 1:3d}   {history[index]:.6e}     {ratio:.6f}")
    print()
    final_ratio = history[13] / history[12]
    print(f"   Measured ratio at step 14: {final_ratio:.6f}. Predicted: {2 / 5:.6f}.")
    print("   The algorithm is telling you the SECOND eigenvalue through the")
    print("   speed at which it finds the first one.")
    print()
    assert abs(final_ratio - 0.4) < 1e-3

    print("   That also tells you when the power method is a bad idea. If the")
    print("   two largest eigenvalues are close, the ratio is near 1 and")
    print("   convergence is glacial. Here is the same code on a matrix whose")
    print("   eigenvalues are 5 and 4.9:")
    print()
    slow = np.array([[5.0, 0.0], [0.0, 4.9]])
    slow_result = power_method(slow, np.array([0.6, 0.8]), tol=1e-10)
    print(f"       eigenvalue ratio 4.9 / 5 = {4.9 / 5:.4f}")
    print(f"       iterations to the same 1e-10 tolerance: {slow_result['iterations']}")
    print(f"       compared with {result['iterations']} for the ratio 0.4 matrix")
    print()
    assert slow_result["iterations"] > 200
    print("   And if they are exactly equal in magnitude the method does not")
    print("   converge at all, because there is no single dominant direction to")
    print("   converge to.")
    print()

    # ---------------------------------------------------------------- 5
    print("5. Why normalise at all? Because of what happens if you do not.")
    print()
    v = start.copy()
    lengths: dict[int, float] = {}
    print("      k    length of A^k v0")
    print("   " + "-" * 34)
    with np.errstate(over="ignore"):
        for k in range(0, 401):
            if k in (0, 10, 50, 100, 200, 300, 400):
                lengths[k] = float(np.linalg.norm(v))
                print(f"     {k:3d}    {lengths[k]:.6e}")
            v = A @ v
    print()
    print("   The direction was right after twenty-five rounds. The LENGTH keeps")
    print("   multiplying by five, and float64 stops at about 1.8e308:")
    print()
    print(f"       largest float64:      {np.finfo(np.float64).max:.6e}")
    print(f"       length after 200:     {lengths[200]:.6e}  (still fine)")
    print(f"       length after 300:     {lengths[300]:.6e}")
    print()
    assert np.isfinite(lengths[200]) and not np.isfinite(lengths[300])
    overflow = start.copy()
    with np.errstate(over="ignore", invalid="ignore"):
        for _ in range(600):
            overflow = A @ overflow
    print(f"       after 600 un-normalised rounds:  {overflow}")
    assert not np.all(np.isfinite(overflow))

    with np.errstate(invalid="ignore"):
        recovered = overflow / np.linalg.norm(overflow)
    print(f"       and normalising it now:          {recovered}")
    assert np.all(np.isnan(recovered))
    print()
    print("       inf divided by inf is nan. The direction was correct at round")
    print("       twenty-five and is now unrecoverable, destroyed by a magnitude")
    print("       nobody asked for. Normalising each round changes no direction")
    print("       and costs one division.")
    print()

    # ---------------------------------------------------------------- 6
    print("6. Against numpy.linalg.eig, which solves the whole problem at once.")
    print()
    values, vectors = np.linalg.eig(A)
    real_values = values.real
    top = int(np.argmax(np.abs(real_values)))
    print(f"       numpy dominant eigenvalue   {real_values[top]:.12f}")
    print(f"       power method eigenvalue     {result['eigenvalue']:.12f}")
    print(f"       difference                  {abs(real_values[top] - result['eigenvalue']):.3e}")
    print()
    print(f"       numpy dominant eigenvector  [{vectors.real[0, top]: .12f}, {vectors.real[1, top]: .12f}]")
    print(f"       power method eigenvector    [{result['vector'][0]: .12f}, {result['vector'][1]: .12f}]")
    print(f"       abs_cosine between them     {abs_cosine(vectors.real[:, top], result['vector']):.15f}")
    print()
    assert abs(real_values[top] - result["eigenvalue"]) < 1e-9
    assert abs_cosine(vectors.real[:, top], result["vector"]) > 1.0 - 1e-12

    print("   Two honest points about that comparison.")
    print()
    print("   The power method found ONE eigenvector; eig found all of them, and")
    print("   for a 2x2 there is no reason at all to use anything else.")
    print()
    print("   But eig needs the matrix as an array and works on all of it. The")
    print("   power method needs only a function that computes A @ v. On a graph")
    print("   with a hundred million nodes the adjacency matrix does not fit in")
    print("   memory as an array, while multiplying by it is just a walk over the")
    print("   edges — which is why the method that looks naive here is the one")
    print("   that survives at scale.")
    print()

    print(f"{SCRIPT}: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 5 — principal component analysis is eigenvectors of a covariance matrix.

Run from inside examples/:

    ../.venv/bin/python3 05_pca_from_covariance.py

Everything so far has been about matrices that were handed to you. This is the
step that makes eigenvectors matter: build a matrix OUT OF DATA, take its
eigenvectors, and the top one is the direction the data is most spread along.

The dataset is invented, and it is stretched along a direction chosen in
advance — 30 degrees from the x-axis. Nothing in the array of coordinates
records that number. PCA has to find it, from 400 pairs of numbers, and the
whole of PCA is five lines long.
"""

from __future__ import annotations

import numpy as np

from dataset import (
    CENTRE,
    ELONGATION_DEG,
    N_POINTS,
    SPREAD_ACROSS,
    SPREAD_ALONG,
    elongation_direction,
    make_cloud,
)
from eigen import abs_cosine, covariance_matrix, direction_degrees, principal_components

SCRIPT = "05_pca_from_covariance.py"


def main() -> None:
    print(f"{SCRIPT}")
    print("=" * 72)
    print()

    cloud = make_cloud()
    truth = elongation_direction()

    # ---------------------------------------------------------------- 1
    print("1. The invented data. 400 points, deliberately cigar-shaped.")
    print()
    print(f"       shape                 {cloud.shape}")
    print(f"       first three points    {np.round(cloud[:3], 6).tolist()}")
    print(f"       column means          {np.round(cloud.mean(axis=0), 6).tolist()}")
    print(f"       built around centre   {CENTRE.tolist()}")
    print()
    print(f"   It was built by drawing a wide spread (sd {SPREAD_ALONG}) along the")
    print(f"   direction {ELONGATION_DEG} degrees, and a narrow spread (sd {SPREAD_ACROSS}) across it.")
    print("   That direction is the answer. It appears nowhere in the array.")
    print()
    print(f"       the answer, kept aside:  [{truth[0]:.6f}, {truth[1]:.6f}]  at {ELONGATION_DEG} degrees")
    print()
    assert cloud.shape == (N_POINTS, 2)

    # ---------------------------------------------------------------- 2
    print("2. Centre the data. This step is not optional and it is the one")
    print("   people skip.")
    print()
    centred = cloud - cloud.mean(axis=0)
    print(f"       means before centring  {np.round(cloud.mean(axis=0), 6).tolist()}")
    print(f"       means after centring   {np.round(centred.mean(axis=0), 12).tolist()}")
    print()
    print("   Covariance is about how the points vary AROUND THEIR OWN MEAN. Skip")
    print("   the subtraction and every product picks up the offset of the cloud")
    print("   from the origin, which here is (5, -2) and has nothing to do with")
    print("   the shape. Section 6 shows exactly how wrong the answer goes.")
    print()
    assert np.allclose(centred.mean(axis=0), 0.0, atol=1e-12)

    # ---------------------------------------------------------------- 3
    print("3. The covariance matrix, from scratch and then from NumPy.")
    print()
    print("       C = (Xc.T @ Xc) / (n - 1)")
    print()
    covariance = covariance_matrix(cloud)
    print("   from scratch:")
    for row in covariance:
        print(f"       [{row[0]: .8f}  {row[1]: .8f}]")
    print()
    print("   numpy.cov(cloud, rowvar=False):")
    reference = np.cov(cloud, rowvar=False)
    for row in reference:
        print(f"       [{row[0]: .8f}  {row[1]: .8f}]")
    print()
    print(f"   identical to 1e-12? {np.allclose(covariance, reference, atol=1e-12)}")
    print()
    assert np.allclose(covariance, reference, atol=1e-12)

    print("   Read the entries. C[0][0] is how much x varies, C[1][1] is how much")
    print("   y varies, and C[0][1] is how much they vary TOGETHER — positive")
    print("   here, meaning that points to the right also tend to be higher up,")
    print("   which is what a cloud tilted upwards at 30 degrees looks like.")
    print()
    print(f"       variance of x        {covariance[0, 0]:.6f}")
    print(f"       variance of y        {covariance[1, 1]:.6f}")
    print(f"       covariance of x, y   {covariance[0, 1]:.6f}")
    print(f"       symmetric?           {np.allclose(covariance, covariance.T, atol=1e-15)}")
    print()
    print("   The symmetry is not a coincidence: entry (0,1) and entry (1,0) are")
    print("   the same sum of products written in the other order. And symmetry")
    print("   is exactly what guarantees the eigenvalues are real and the")
    print("   eigenvectors at right angles, which is why PCA always works and")
    print("   never comes back with a complex answer you have to interpret.")
    print()
    assert np.allclose(covariance, covariance.T, atol=1e-15)

    # ---------------------------------------------------------------- 4
    print("4. Take the eigenvectors. That is PCA. There is no step five.")
    print()
    variances, directions = principal_components(cloud)
    print(f"       numpy.linalg.eigh eigenvalues, sorted large to small:")
    print(f"           {np.round(variances, 8).tolist()}")
    print(f"       eigenvectors, as columns:")
    for row in directions:
        print(f"           [{row[0]: .8f}  {row[1]: .8f}]")
    print()
    top = directions[:, 0]
    second = directions[:, 1]
    print(f"       top component      [{top[0]: .8f}, {top[1]: .8f}]")
    print(f"       its direction      {direction_degrees(top):.6f} degrees")
    print(f"       the truth          {ELONGATION_DEG} degrees")
    print(f"       error              {abs(direction_degrees(top) - ELONGATION_DEG):.6f} degrees")
    print()
    similarity = abs_cosine(top, truth)
    print(f"       abs_cosine(top component, true direction) = {similarity:.10f}")
    print()
    assert similarity > 0.999
    assert abs(direction_degrees(top) - ELONGATION_DEG) < 1.0

    print("   Found, from 400 pairs of coordinates and nothing else.")
    print()
    print("   And note the SIGN of what came back:")
    print(f"       top component  [{top[0]: .6f}, {top[1]: .6f}]")
    print(f"       true direction [{truth[0]: .6f}, {truth[1]: .6f}]")
    if float(np.dot(top, truth)) < 0:
        print("       They point OPPOSITE ways along the same line. numpy.allclose")
        print(f"       says {np.allclose(top, truth)}, and the answer is nonetheless exactly right.")
        print("       A principal component names an AXIS, not an arrow, and any")
        print("       code that cares which end is which has a bug waiting.")
    print()
    assert np.dot(top, truth) < 0  # observed: eigh returned the reversed sign here

    # ---------------------------------------------------------------- 5
    print("5. What the eigenVALUES mean here: variance along each axis.")
    print()
    print(f"       eigenvalue 1   {variances[0]:.6f}   sqrt = {np.sqrt(variances[0]):.6f}   (built with sd {SPREAD_ALONG})")
    print(f"       eigenvalue 2   {variances[1]:.6f}   sqrt = {np.sqrt(variances[1]):.6f}   (built with sd {SPREAD_ACROSS})")
    print()
    print("   The eigenvalue IS the variance along its own eigenvector, so its")
    print("   square root is the standard deviation — and both come back close")
    print("   to the numbers the cloud was built with. Not exactly equal: 400")
    print("   samples estimate a spread, they do not reproduce it.")
    print()
    proportion = variances / variances.sum()
    print(f"       proportion of variance explained: {np.round(proportion, 8).tolist()}")
    print(f"       the first component alone carries {100 * proportion[0]:.4f}% of it")
    print()
    print("   That is the sentence behind every 'we reduced 768 dimensions to 50'")
    print("   claim you will read: sort the eigenvalues, keep the ones that add")
    print("   up to enough of the total, and throw the rest away.")
    print()
    assert proportion[0] > 0.97
    assert abs(np.sqrt(variances[0]) - SPREAD_ALONG) < 0.2
    assert abs(np.sqrt(variances[1]) - SPREAD_ACROSS) < 0.1

    print("   Check it directly: project every point onto each component and")
    print("   measure the spread of what comes out.")
    print()
    projected = centred @ directions
    print(f"       spread along component 1  {projected[:, 0].std(ddof=1):.6f}")
    print(f"       spread along component 2  {projected[:, 1].std(ddof=1):.6f}")
    print(f"       correlation between them  {np.corrcoef(projected.T)[0, 1]:.3e}")
    print()
    print("   The two projections are uncorrelated to within rounding, which is")
    print("   the other half of what PCA buys you: not just the best directions,")
    print("   but directions along which the data carries no shared information.")
    print()
    assert abs(projected[:, 0].std(ddof=1) - np.sqrt(variances[0])) < 1e-9
    assert abs(float(np.corrcoef(projected.T)[0, 1])) < 1e-12

    # ---------------------------------------------------------------- 6
    print("6. What forgetting to centre actually costs.")
    print()
    uncentred = (cloud.T @ cloud) / (N_POINTS - 1)
    bad_variances, bad_directions = np.linalg.eigh(uncentred)
    order = np.argsort(bad_variances)[::-1]
    bad_top = bad_directions[:, order][:, 0]
    print("       Xc.T @ Xc / (n-1) with NO centring:")
    for row in uncentred:
        print(f"           [{row[0]: .8f}  {row[1]: .8f}]")
    print(f"       its top eigenvector points at {direction_degrees(bad_top):.6f} degrees")
    print(f"       the correct answer is         {ELONGATION_DEG} degrees")
    print(f"       error                         {abs(direction_degrees(bad_top) - ELONGATION_DEG):.6f} degrees")
    print()
    centre_angle = float(np.degrees(np.arctan2(CENTRE[1], CENTRE[0])) % 180.0)
    bad_angle = direction_degrees(bad_top)
    print(f"       the cloud's CENTRE lies at {centre_angle:.6f} degrees from the origin")
    print(f"       distance from the uncentred answer to that:  {abs(bad_angle - centre_angle):.6f} degrees")
    print(f"       distance from the uncentred answer to the truth: {abs(bad_angle - ELONGATION_DEG):.6f} degrees")
    print()
    print("       So the uncentred answer sits close to the direction of the")
    print("       cloud's OFFSET from the origin and nowhere near its shape. That")
    print("       is what the calculation measured, because without centring the")
    print("       squared offset (5^2 + 2^2 = 29) swamps the actual spread (8.4).")
    print("       No exception, no warning, a confident wrong answer.")
    print()
    assert abs(bad_angle - ELONGATION_DEG) > 5.0
    assert abs(bad_angle - centre_angle) < abs(bad_angle - ELONGATION_DEG)

    # ---------------------------------------------------------------- 7
    print("7. Against the two NumPy routines.")
    print()
    values_eig, vectors_eig = np.linalg.eig(covariance)
    print(f"       numpy.linalg.eig  values {values_eig}  dtype {values_eig.dtype}")
    values_eigh, vectors_eigh = np.linalg.eigh(covariance)
    print(f"       numpy.linalg.eigh values {values_eigh}  dtype {values_eigh.dtype}")
    print()
    print("   Same numbers, different packaging. eigh knows the input is")
    print("   symmetric, so it returns float64 in ascending order; eig does not")
    print("   assume it, so it returns complex128 unordered. For a covariance")
    print("   matrix — always symmetric, by construction — eigh is the right")
    print("   call every time.")
    print()
    eig_top = vectors_eig.real[:, int(np.argmax(values_eig.real))]
    print(f"       eig  top component, abs_cosine with the truth  {abs_cosine(eig_top, truth):.10f}")
    print(f"       eigh top component, abs_cosine with the truth  {abs_cosine(top, truth):.10f}")
    print()
    assert abs_cosine(eig_top, truth) > 0.999
    assert np.allclose(np.sort(values_eig.real), np.sort(values_eigh), atol=1e-12)

    print("   The AI connection, in one line: replace this 2-column cloud with a")
    print("   matrix of 768-dimensional sentence embeddings and nothing about")
    print("   the method changes. The covariance matrix becomes 768 by 768, its")
    print("   top eigenvectors are the directions those embeddings actually vary")
    print("   along, and the eigenvalues tell you how many of the 768 dimensions")
    print("   are carrying real information rather than noise.")
    print()

    print(f"{SCRIPT}: every assertion held.")


if __name__ == "__main__":
    main()

"""Exercise 6 — eig or eigh, and what the choice is actually worth.

Run from inside examples/:

    ../.venv/bin/python3 06_eig_against_eigh.py

NumPy offers four routines for this job and they are not interchangeable.
This script runs all four on the same inputs and measures the one difference
that is usually quoted without a number attached: how much faster eigh is.

Timings are measured here, on one machine, on one day. They are printed rather
than asserted, because a test that asserts milliseconds fails on somebody
else's laptop for no good reason. The SHAPES and DTYPES are asserted, because
those are properties of the routines rather than of the hardware.
"""

from __future__ import annotations

import platform
import sys
import time

import numpy as np

from dataset import A, SYMMETRIC

SCRIPT = "06_eig_against_eigh.py"

SIZE = 400
REPEATS = 5


def best_of(function, matrix, repeats: int = REPEATS) -> float:
    """Fastest of several runs, in seconds.

    The fastest run rather than the mean, because everything that makes a run
    slower — another process waking up, a page fault, the CPU changing its
    clock speed — is noise added on top. There is no source of noise that
    makes a run faster than the machine can do it, so the minimum is the
    cleanest estimate of the real cost.
    """
    timings = []
    for _ in range(repeats):
        start = time.perf_counter()
        function(matrix)
        timings.append(time.perf_counter() - start)
    return min(timings)


def main() -> None:
    print(f"{SCRIPT}")
    print("=" * 72)
    print()
    print(f"   python   {platform.python_version()}")
    print(f"   numpy    {np.__version__}")
    print(f"   platform {platform.platform()}")
    print(f"   exe      {sys.executable.rsplit('/', 3)[-1]}")
    print()

    # ---------------------------------------------------------------- 1
    print("1. The four routines, and what each one is for.")
    print()
    print("   numpy.linalg.eig      any square matrix; values AND vectors")
    print("   numpy.linalg.eigvals  any square matrix; values only")
    print("   numpy.linalg.eigh     symmetric input only; values AND vectors")
    print("   numpy.linalg.eigvalsh symmetric input only; values only")
    print()
    print("   On the symmetric matrix [[2, 1], [1, 2]]:")
    print()
    for name, function in (
        ("eig", np.linalg.eig),
        ("eigh", np.linalg.eigh),
    ):
        values, vectors = function(SYMMETRIC)
        print(f"       {name:<6} values {values}  dtype {values.dtype}")
        print(f"       {'':<6} sorted ascending? {bool(np.all(np.diff(values.real) >= 0))}")
    for name, function in (
        ("eigvals", np.linalg.eigvals),
        ("eigvalsh", np.linalg.eigvalsh),
    ):
        values = function(SYMMETRIC)
        print(f"       {name:<8} values {values}  dtype {values.dtype}")
    print()
    assert np.linalg.eig(SYMMETRIC)[0].dtype == np.complex128
    assert np.linalg.eigh(SYMMETRIC)[0].dtype == np.float64
    assert np.linalg.eigvalsh(SYMMETRIC).dtype == np.float64
    assert np.all(np.diff(np.linalg.eigh(SYMMETRIC)[0]) >= 0)

    print("   Three differences that matter in practice:")
    print()
    print("     * eigh returns float64; eig returns complex128 even when every")
    print("       eigenvalue is real. Exercise 2 covers that at length.")
    print("     * eigh returns its values sorted ascending. eig makes no")
    print("       ordering promise at all, so 'the largest eigenvalue' needs an")
    print("       argmax rather than an index.")
    print("     * eigh reads only ONE triangle of the input and assumes the")
    print("       other matches. Feed it a non-symmetric matrix and it does not")
    print("       complain — it answers a question about a different matrix.")
    print()

    # ---------------------------------------------------------------- 2
    print("2. eigh on non-symmetric input: silently the wrong answer.")
    print()
    print("   A = [[4, 1], [2, 3]] is NOT symmetric. Its real eigenvalues are 5 and 2.")
    print()
    honest = np.sort(np.linalg.eig(A)[0].real)
    wrong = np.linalg.eigh(A)[0]
    print(f"       eig  on A -> {honest}   (correct)")
    print(f"       eigh on A -> {wrong}   (no error, no warning)")
    print()
    print("   eigh took the lower triangle, [[4, ...], [2, 3]], assumed the")
    print("   upper matched it, and solved [[4, 2], [2, 3]] instead — a")
    print("   different matrix with different eigenvalues.")
    print()
    substitute = np.array([[4.0, 2.0], [2.0, 3.0]])
    print(f"       eigvalsh([[4, 2], [2, 3]]) -> {np.linalg.eigvalsh(substitute)}")
    print(f"       matches what eigh returned for A? {np.allclose(wrong, np.linalg.eigvalsh(substitute))}")
    print()
    print("   So check symmetry before reaching for eigh, or know for structural")
    print("   reasons that it holds — as it does for every covariance matrix.")
    print()
    assert not np.allclose(np.sort(wrong), honest)
    assert np.allclose(wrong, np.linalg.eigvalsh(substitute))

    # ---------------------------------------------------------------- 3
    print(f"3. What eigh is worth, measured on a {SIZE} by {SIZE} symmetric float64 matrix.")
    print()
    rng = np.random.default_rng(7)
    noise = rng.normal(size=(SIZE, SIZE))
    symmetric = (noise + noise.T) / 2.0
    assert np.allclose(symmetric, symmetric.T)

    print(f"       shape {symmetric.shape}, dtype {symmetric.dtype}, symmetric: True")
    print(f"       best of {REPEATS} runs each:")
    print()
    results = {}
    for name, function in (
        ("numpy.linalg.eig", np.linalg.eig),
        ("numpy.linalg.eigh", np.linalg.eigh),
        ("numpy.linalg.eigvals", np.linalg.eigvals),
        ("numpy.linalg.eigvalsh", np.linalg.eigvalsh),
    ):
        seconds = best_of(function, symmetric)
        results[name] = seconds
        print(f"         {name:<24} {seconds * 1000:8.2f} ms")
    print()
    ratio = results["numpy.linalg.eig"] / results["numpy.linalg.eigh"]
    print(f"       eig / eigh = {ratio:.2f}x on this run, on this machine.")
    print()
    print("   float64 on purpose. A float32 or an integer array would be a")
    print("   different measurement, and mixing them would make the number")
    print("   meaningless. Nothing here is asserted — the ratio is real but it")
    print("   is one machine on one day, and your number will differ.")
    print()
    print("   Both answers agree, which IS asserted:")
    values_eig = np.sort(np.linalg.eig(symmetric)[0].real)
    values_eigh = np.linalg.eigh(symmetric)[0]
    difference = float(np.abs(values_eig - values_eigh).max())
    print(f"       largest disagreement across all {SIZE} eigenvalues: {difference:.3e}")
    print()
    assert difference < 1e-10

    # ---------------------------------------------------------------- 4
    print("4. Everything else in this area, described honestly and NOT run here.")
    print()
    print("   None of the following is installed in this lab, and no output from")
    print("   any of them is reproduced anywhere in this day. They are described")
    print("   from their own documentation, and that is all.")
    print()
    for name, note in (
        (
            "scipy.linalg.eig / eigh",
            "the same jobs with more knobs — a generalized problem A v = lambda B v, "
            "the option to ask for only a range of eigenvalues, and a choice of "
            "LAPACK driver. Reach for it when NumPy's version does not take the "
            "argument you need. Free and open source, BSD 3-Clause.",
        ),
        (
            "scipy.sparse.linalg.eigs / eigsh",
            "for matrices too large to hold densely. Asks for the k largest or "
            "smallest eigenvalues rather than all of them, and needs only the "
            "ability to multiply by the matrix — the same requirement the power "
            "method has, which is not a coincidence: this is an industrial-grade "
            "relative of it. Free and open source, BSD 3-Clause.",
        ),
        (
            "torch.linalg.eig / eigh",
            "the same interface on tensors, so it runs on a GPU and takes part in "
            "automatic differentiation. Choose it when the eigendecomposition sits "
            "inside a model rather than beside it. Free and open source, BSD-style.",
        ),
        (
            "sklearn.decomposition.PCA",
            "PCA as a fitted object, with the centring, the sorting, the "
            "variance ratios and the transform handled for you — and with a "
            "singular value decomposition underneath rather than an explicit "
            "covariance matrix, which is more accurate on ill-conditioned data. "
            "Use it for real work; the five lines in exercise 5 are for "
            "understanding what it does. Free and open source, BSD 3-Clause.",
        ),
    ):
        print(f"     {name}")
        for line in _wrap(note, 66):
            print(f"         {line}")
        print()

    print(f"{SCRIPT}: every assertion held.")


def _wrap(text: str, width: int) -> list[str]:
    words = text.split()
    lines: list[str] = []
    current = ""
    for word in words:
        candidate = f"{current} {word}".strip()
        if len(candidate) > width and current:
            lines.append(current)
            current = word
        else:
            current = candidate
    if current:
        lines.append(current)
    return lines


if __name__ == "__main__":
    main()

"""Make this directory's own eigen.py the one its tests import.

Both `examples/` and `starter/` contain modules called `eigen` and `dataset`,
and pytest imports test files by putting their directory on `sys.path`.
Without this file, running a bare `pytest` across both directories at once
would import whichever `eigen` was seen first and then reuse it for the other
suite — so the starter tests would silently pass against the reference
solution instead of skipping. That is a wrong answer with a green tick on it,
which is the worst kind.

So: put this directory first on the import path, and drop any already-imported
`eigen` or `dataset` that came from somewhere else.

Section 4 of `tests/run_tests.sh` checks that this still works, by comparing
the skip count from `pytest starter` against the skip count from a bare
`pytest` over the whole lab. If this guard ever stops working, that check goes
red rather than the lab quietly lying to a learner.
"""

import sys
from pathlib import Path

HERE = str(Path(__file__).parent.resolve())

if HERE in sys.path:
    sys.path.remove(HERE)
sys.path.insert(0, HERE)

for name in ("eigen", "dataset"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""The matrices this lab works with, and the answers worked out by hand.

Everything here is invented and every number is small enough to check on
paper. Nothing is read from disk and nothing is downloaded.

The star of the lab is `A`. It was chosen so that the characteristic equation
factorises over the integers, which means the whole eigenvalue calculation can
be done by hand in about a minute and then checked against NumPy:

    A = [[4, 1],
         [2, 3]]

    trace  = 4 + 3 = 7
    det    = 4*3 - 1*2 = 10
    characteristic equation:  lambda^2 - 7*lambda + 10 = 0
                              (lambda - 5)(lambda - 2) = 0
    eigenvalues:              5 and 2

    for lambda = 5:  A - 5I = [[-1, 1], [2, -2]]  ->  -x + y = 0  ->  (1, 1)
    for lambda = 2:  A - 2I = [[ 2, 1], [2,  1]]  ->  2x + y = 0  ->  (1, -2)

Check both by hand:

    A @ (1,  1) = (4 + 1, 2 + 3) = ( 5,  5) = 5 * (1,  1)
    A @ (1, -2) = (4 - 2, 2 - 6) = ( 2, -4) = 2 * (1, -2)
"""

from __future__ import annotations

import numpy as np

# --------------------------------------------------------------------------
# The worked 2x2
# --------------------------------------------------------------------------

A = np.array([[4.0, 1.0], [2.0, 3.0]])

#: The eigenvalues of `A`, worked out by hand above, largest first.
A_EIGENVALUES = (5.0, 2.0)

#: One eigenvector per eigenvalue, in the same order. NOT normalised, and
#: deliberately so: any non-zero multiple of these is equally correct, which
#: is the point exercise 4 is built around.
A_EIGENVECTORS = ((1.0, 1.0), (1.0, -2.0))

#: The direction of each eigenvector as an angle in degrees, measured from the
#: positive x-axis and folded into [0, 180) because a direction and its
#: reverse lie on the same line.
#:   atan2( 1,  1) =  45 degrees
#:   atan2(-2,  1) = -63.4349... degrees, which is 116.5650... modulo 180
A_EIGEN_ANGLES_DEG = (45.0, 116.56505117707799)

# --------------------------------------------------------------------------
# The standard transformations from Day 102, and what each one does to a line
# --------------------------------------------------------------------------

IDENTITY = np.eye(2)
UNIFORM_SCALE = 2.0 * np.eye(2)
NON_UNIFORM_SCALE = np.diag([3.0, 0.5])
REFLECTION_IN_X = np.array([[1.0, 0.0], [0.0, -1.0]])
SHEAR = np.array([[1.0, 1.0], [0.0, 1.0]])
ROTATION_90 = np.array([[0.0, -1.0], [1.0, 0.0]])
ROTATION_60 = np.array(
    [
        [0.5, -np.sqrt(3.0) / 2.0],
        [np.sqrt(3.0) / 2.0, 0.5],
    ]
)
PROJECTION_ONTO_X = np.array([[1.0, 0.0], [0.0, 0.0]])

#: A symmetric matrix, for the two guarantees symmetry buys: real eigenvalues
#: and orthogonal eigenvectors. trace 4, det 3, so lambda^2 - 4L + 3 = 0 and
#: the eigenvalues are 3 and 1.
SYMMETRIC = np.array([[2.0, 1.0], [1.0, 2.0]])

#: A larger symmetric matrix, so the orthogonality claim is checked on
#: something that is not a 2x2 special case.
SYMMETRIC_3X3 = np.array(
    [
        [4.0, 1.0, 2.0],
        [1.0, 3.0, 0.0],
        [2.0, 0.0, 5.0],
    ]
)

#: name -> (matrix, one-line description of what it does to directions)
STANDARD_TRANSFORMATIONS = {
    "identity": (IDENTITY, "leaves everything alone; every direction is an eigenvector"),
    "uniform scale 2x": (UNIFORM_SCALE, "stretches everything equally; every direction survives"),
    "non-uniform scale": (NON_UNIFORM_SCALE, "stretches x by 3 and squashes y by half"),
    "reflection in x-axis": (REFLECTION_IN_X, "flips the sign of y; one eigenvalue is negative"),
    "shear": (SHEAR, "slides the top of the square sideways; only ONE eigendirection"),
    "rotation 90 degrees": (ROTATION_90, "turns every vector; NO real eigenvector at all"),
    "rotation 60 degrees": (ROTATION_60, "the same story at a different angle"),
    "projection onto x-axis": (PROJECTION_ONTO_X, "flattens the plane onto a line; det 0, eigenvalue 0"),
}

# --------------------------------------------------------------------------
# The tiny 2-D dataset for the PCA demonstration
# --------------------------------------------------------------------------

#: The seed for every random draw in this lab. Fixed so that every number in
#: expected-output/ is reproducible on any machine.
SEED = 2106

#: The direction the invented cloud is deliberately stretched along, in
#: degrees from the positive x-axis. PCA has to rediscover this number from
#: the data alone.
ELONGATION_DEG = 30.0

#: Standard deviation of the cloud along the elongation direction, and across
#: it. The ratio 3.0 : 0.4 is what makes the cloud visibly cigar-shaped.
SPREAD_ALONG = 3.0
SPREAD_ACROSS = 0.4

#: How many points, and where the cloud is centred. The centre is deliberately
#: away from the origin so that forgetting to subtract the mean is a mistake
#: with visible consequences.
N_POINTS = 400
CENTRE = np.array([5.0, -2.0])


def elongation_direction() -> np.ndarray:
    """The unit vector the cloud is stretched along. The answer PCA must find."""
    radians = np.radians(ELONGATION_DEG)
    return np.array([np.cos(radians), np.sin(radians)])


def make_cloud() -> np.ndarray:
    """Build the invented 2-D dataset: shape (N_POINTS, 2), one point per row.

    The construction is the whole trick, and it is worth reading rather than
    running: take a unit vector `along` pointing at ELONGATION_DEG, and the
    unit vector `across` at right angles to it. Draw a wide spread of numbers
    to travel `along` and a narrow spread to travel `across`, add them, and
    shift the whole cloud to CENTRE. The result is a cigar-shaped cloud whose
    long axis is known exactly, because it was put there on purpose.

    Nothing about the direction is stored in the array that comes back. The
    covariance matrix has to recover it from 400 pairs of coordinates.
    """
    rng = np.random.default_rng(SEED)
    along = elongation_direction()
    across = np.array([-along[1], along[0]])

    travel_along = rng.normal(0.0, SPREAD_ALONG, size=N_POINTS)
    travel_across = rng.normal(0.0, SPREAD_ACROSS, size=N_POINTS)

    return travel_along[:, None] * along + travel_across[:, None] * across + CENTRE


#: The starting vector for the power method. Drawn from the same seeded
#: generator so the iteration count in expected-output/ is reproducible.
def power_method_start() -> np.ndarray:
    """A random unit vector to start the power method from.

    Random on purpose: the whole claim is that ALMOST ANY starting vector
    converges to the dominant eigenvector, so starting from something chosen
    to work would prove nothing.
    """
    rng = np.random.default_rng(106)
    v = rng.normal(size=2)
    return v / np.linalg.norm(v)

"""Eigenvalues and eigenvectors, written out by hand.

The reference implementation. Every function here does something NumPy will
do for you in one call; the point of writing them is that afterwards you know
what that one call is doing, and you know which of its answers are forced by
mathematics and which are arbitrary choices the library made for you.

Nothing in this file imports anything beyond NumPy, and no function here
touches the file system, the clock or the network.
"""

from __future__ import annotations

import numpy as np

# --------------------------------------------------------------------------
# Measuring whether a vector kept its direction
# --------------------------------------------------------------------------


def abs_cosine(u, v) -> float:
    """The absolute cosine of the angle between two vectors: 1.0 when they lie
    on the same LINE, 0.0 when they are at right angles.

    The absolute value is the whole point, and it is the single most important
    habit in this lab. An eigenvector is defined only up to sign and scale: if
    v satisfies A v = lambda v, then so does -v, and so does 3.7 v. NumPy
    returns *a* unit-length eigenvector, and which of the two possible signs
    it hands back is a detail of the LAPACK routine underneath, not a fact
    about the matrix. So a test that compares components will fail on a
    correct answer roughly half the time. Compare DIRECTIONS instead, and
    "direction" means "line", which is what taking the absolute value does.
    """
    u = np.asarray(u, dtype=float).ravel()
    v = np.asarray(v, dtype=float).ravel()
    nu = float(np.linalg.norm(u))
    nv = float(np.linalg.norm(v))
    if nu == 0.0 or nv == 0.0:
        raise ValueError("the zero vector has no direction, so no angle is defined")
    return abs(float(np.dot(u, v)) / (nu * nv))


def deviation_degrees(u, v) -> float:
    """How far apart two directions are, in degrees, ignoring sign.

    Returns a number in [0, 90]. Zero means "these lie on the same line",
    which is exactly the test for "this vector kept its direction".
    """
    c = min(1.0, abs_cosine(u, v))
    return float(np.degrees(np.arccos(c)))


def direction_degrees(v) -> float:
    """The direction of a 2-D vector as an angle in [0, 180) degrees.

    Folded modulo 180 because a vector and its reverse describe the same line,
    and a line is what an eigenvector really names.
    """
    v = np.asarray(v, dtype=float).ravel()
    return float(np.degrees(np.arctan2(v[1], v[0])) % 180.0)


def fan_of_directions(n: int) -> np.ndarray:
    """`n` unit vectors spread evenly around the circle, as rows of an (n, 2) array."""
    angles = np.radians(np.linspace(0.0, 360.0, n, endpoint=False))
    return np.column_stack([np.cos(angles), np.sin(angles)])


def sweep_deviations(matrix, angles_deg) -> tuple[np.ndarray, np.ndarray]:
    """For each angle, how far the matrix knocks that direction off its line.

    Returns (deviations, collapsed) where `deviations` is in degrees and
    `collapsed` is a boolean mask marking directions the matrix sends to the
    origin. Those get `numpy.nan` for their deviation, and that is the honest
    answer rather than a defect: the zero vector has no direction, so "did it
    keep its direction?" has nothing to compare against. A collapsed direction
    is an eigenvector with eigenvalue 0 — it is squashed rather than turned —
    and the caller has to decide what to do about it, because this measurement
    cannot.

    Vectorised on purpose: this is called with 180,000 angles in exercise 1,
    and a Python loop over that is the difference between instant and tedious.
    """
    radians = np.radians(np.asarray(angles_deg, dtype=float))
    vectors = np.column_stack([np.cos(radians), np.sin(radians)])
    outputs = vectors @ np.asarray(matrix, dtype=float).T

    input_norms = np.linalg.norm(vectors, axis=1)
    output_norms = np.linalg.norm(outputs, axis=1)
    collapsed = output_norms <= 1e-12 * max(1.0, float(np.abs(matrix).max()))

    dots = np.einsum("ij,ij->i", vectors, outputs)
    deviations = np.full(len(radians), np.nan)
    live = ~collapsed
    cosines = np.clip(
        np.abs(dots[live] / (input_norms[live] * output_norms[live])), 0.0, 1.0
    )
    deviations[live] = np.degrees(np.arccos(cosines))
    return deviations, collapsed


def eigen_lines_by_sweep(
    matrix,
    step: float = 0.001,
    keep_below: float = 0.01,
    gap: float = 1.0,
) -> dict:
    """Find the eigendirections of a 2x2 by brute-force measurement.

    Sweep every direction from 0 to 180 degrees, keep the ones the matrix
    barely moves, group the survivors into clusters, and report the very best
    angle in each cluster.

    Returns a dict with:
        verdict   one of "none", "every direction", "some"
        lines     one representative angle in degrees per distinct eigen-line
        collapsed angles the matrix sends to the origin (eigenvalue 0)
        fraction  what proportion of the swept directions kept their line

    `len(result["lines"])` is the number of distinct eigen-lines, which the
    algebra calls the geometric multiplicity count. That COUNT is the reliable
    output. The angles are approximate — good to roughly the width of the
    surviving band, which is a few hundredths of a degree here — because they
    come from sampling a smooth curve rather than from solving anything. The
    shear's single line comes back as 0.005 rather than 0.000 for exactly that
    reason: its deviation curve is not symmetric about the eigendirection, so
    the surviving band is not centred on it. Use this to find out HOW MANY
    directions survive and roughly where; use the algebra to get them exactly.

    Four details that are not obvious and that a first attempt gets wrong:

      * Only 0 to 180 is swept, because a direction and its reverse are the
        same line and sweeping the full circle would double-count every answer.
      * The clustering wraps around, because 179.999 degrees and 0.001 degrees
        are neighbours on a line-through-the-origin, not opposite ends of a
        range. A shear's single eigen-line sits exactly on that seam, so a
        version without the wrap reports two lines where there is one.
      * `keep_below` cannot be made arbitrarily small. Near an eigendirection
        the deviation curve is smooth, so the sampled grid usually straddles
        the true minimum rather than landing on it — the second eigen-line of
        this lab's matrix sits at 116.56505 degrees and the nearest sample at
        0.001-degree spacing still deviates by 7.7e-05. A threshold of 1e-06
        would find nothing there and report one eigen-line instead of two.
        0.01 degrees is far below the several-degree deviations of every
        non-eigendirection and far above the sampling error.
      * A matrix that keeps EVERY direction — the identity, or any uniform
        scaling — makes every sample a survivor, and they all merge into one
        giant cluster. Reporting "1 eigen-line" there would be exactly wrong,
        so that case is detected by the survivor fraction and named.
    """
    angles = np.arange(0.0, 180.0, step)
    deviations, collapsed = sweep_deviations(matrix, angles)
    keep = deviations < keep_below  # nan compares False, so collapsed is excluded
    fraction = float(np.count_nonzero(keep)) / len(angles)
    collapsed_angles = [float(a) for a in angles[collapsed]]

    if fraction > 0.99:
        return {
            "verdict": "every direction",
            "lines": [],
            "collapsed": collapsed_angles,
            "fraction": fraction,
        }
    if not np.any(keep):
        return {
            "verdict": "none",
            "lines": [],
            "collapsed": collapsed_angles,
            "fraction": fraction,
        }

    surviving_angles = angles[keep]
    surviving_deviations = deviations[keep]

    clusters: list[list[int]] = [[0]]
    for index in range(1, len(surviving_angles)):
        if surviving_angles[index] - surviving_angles[index - 1] <= gap:
            clusters[-1].append(index)
        else:
            clusters.append([index])

    # The seam: the last cluster may be the same line as the first one.
    if len(clusters) > 1:
        wrapped = 180.0 - surviving_angles[clusters[-1][0]] + surviving_angles[clusters[0][0]]
        if wrapped <= gap:
            clusters[0] = clusters[-1] + clusters[0]
            clusters.pop()

    # Report each cluster by its CENTRE, not by its lowest sample. Near an
    # eigendirection the deviation is so small that arccos rounds a whole run
    # of samples to exactly 0.0 — fifteen of them for the shear — so "the
    # sample with the smallest deviation" is decided by a floating-point tie
    # and lands wherever the tie happens to break. The centre of the surviving
    # band is stable, and the band is symmetric about the true eigendirection.
    #
    # Angles are folded RELATIVE TO THE CLUSTER'S OWN FIRST MEMBER, so that a
    # cluster straddling the 0/180 seam becomes contiguous. Folding at a fixed
    # boundary instead — say, everything above 90 minus 180 — splits the
    # perfectly ordinary cluster sitting at 90 degrees and reports it as two.
    best = []
    for cluster in clusters:
        raw = surviving_angles[np.array(cluster)]
        anchor = float(raw[0])
        relative = raw - anchor
        relative = np.where(relative > 90.0, relative - 180.0, relative)
        relative = np.where(relative < -90.0, relative + 180.0, relative)
        centre = anchor + (float(relative.min()) + float(relative.max())) / 2.0
        best.append(round(centre % 180.0, 6))
    return {
        "verdict": "some",
        "lines": sorted(best),
        "collapsed": collapsed_angles,
        "fraction": fraction,
    }


# --------------------------------------------------------------------------
# Solving a 2x2 by hand
# --------------------------------------------------------------------------


def characteristic_coefficients(matrix) -> tuple[float, float]:
    """The two numbers that define the characteristic equation of a 2x2.

    Returns (trace, determinant), which is everything you need, because for a
    2x2 the determinant of (A - lambda*I) always works out to

        lambda^2 - (trace)*lambda + (determinant)

    Derive it once and you never have to again. Writing A as [[a, b], [c, d]]:

        A - lambda*I = [[a - lambda, b], [c, d - lambda]]
        det          = (a - lambda)(d - lambda) - b*c
                     = lambda^2 - (a + d)*lambda + (a*d - b*c)

    and (a + d) is the trace while (a*d - b*c) is the determinant of A.
    """
    m = np.asarray(matrix, dtype=float)
    if m.shape != (2, 2):
        raise ValueError(f"this hand method is for 2x2 matrices only, got {m.shape}")
    trace = float(m[0, 0] + m[1, 1])
    determinant = float(m[0, 0] * m[1, 1] - m[0, 1] * m[1, 0])
    return trace, determinant


def eigenvalues_2x2(matrix) -> tuple[complex, complex]:
    """Solve the characteristic equation with the school quadratic formula.

    lambda = (trace +/- sqrt(trace^2 - 4*det)) / 2

    The discriminant decides the whole character of the matrix:

      * positive  -> two different real eigenvalues, two separate eigendirections
      * zero      -> one repeated real eigenvalue; there may be one
                     eigendirection or every direction, and the eigenvalue
                     alone cannot tell you which
      * negative  -> no real eigenvalues at all, which geometrically means the
                     matrix knocks EVERY direction off its line. A rotation in
                     the plane is the example to remember.

    Returns complex numbers always, so that the negative-discriminant case
    needs no special handling by the caller. Take `.real` when you have
    already checked that the imaginary part is zero.
    """
    trace, determinant = characteristic_coefficients(matrix)
    discriminant = trace * trace - 4.0 * determinant
    root = np.emath.sqrt(discriminant)  # returns a complex root when negative
    first = (trace + root) / 2.0
    second = (trace - root) / 2.0
    return complex(first), complex(second)


def eigenvector_2x2(matrix, eigenvalue: float) -> np.ndarray:
    """A non-zero solution v of (A - lambda*I) v = 0, for a real eigenvalue.

    The reasoning, which is the part worth carrying away. Once lambda is a
    root of the characteristic equation, det(A - lambda*I) is zero, and Day
    102 taught what a zero determinant means: the matrix squashes the plane
    onto a line (or onto a point). Every vector on the line that gets squashed
    to the origin is an eigenvector, so a whole line of solutions exists and
    we only have to name one point on it.

    Concretely: writing B = A - lambda*I as [[p, q], [r, s]], the row [p, q]
    says p*x + q*y = 0. The vector (-q, p) satisfies that for free. If that
    row happens to be all zeros it tells us nothing, so we try the other row.
    If BOTH rows are zero then B is the zero matrix, every direction works,
    and we return (1, 0) with a note in the docstring rather than pretending
    the answer is unique.

    Returned as a unit vector, matching what NumPy does, and with the sign
    left exactly as the arithmetic produced it — because the sign is not
    determined by anything, and pretending otherwise is the trap this lab is
    built to teach.
    """
    m = np.asarray(matrix, dtype=float)
    if m.shape != (2, 2):
        raise ValueError(f"this hand method is for 2x2 matrices only, got {m.shape}")
    shifted = m - eigenvalue * np.eye(2)

    scale = max(1.0, float(np.abs(m).max()))
    tiny = 1e-12 * scale

    for row in (shifted[0], shifted[1]):
        candidate = np.array([-row[1], row[0]])
        if float(np.linalg.norm(candidate)) > tiny:
            return candidate / float(np.linalg.norm(candidate))

    # Both rows vanished: A was lambda*I, so EVERY direction is an
    # eigenvector. One representative is as good as another.
    return np.array([1.0, 0.0])


def solve_2x2(matrix) -> tuple[tuple[complex, complex], list[np.ndarray] | None]:
    """The whole hand method in one call: eigenvalues, and eigenvectors if real.

    Returns (eigenvalues, eigenvectors) where eigenvectors is None when the
    eigenvalues are not real — because in that case there is no real vector
    that keeps its direction, and returning something anyway would be a lie
    dressed as an answer.
    """
    values = eigenvalues_2x2(matrix)
    if any(abs(value.imag) > 1e-12 for value in values):
        return values, None
    vectors = [eigenvector_2x2(matrix, value.real) for value in values]
    return values, vectors


# --------------------------------------------------------------------------
# The power method
# --------------------------------------------------------------------------


def power_method(
    matrix,
    start=None,
    tol: float = 1e-10,
    max_iter: int = 1000,
) -> dict:
    """Find the dominant eigenvector by multiplying, over and over.

    The whole algorithm is three lines, and the reason it works is one
    sentence: write the starting vector as a mixture of the eigenvectors, and
    every application of A multiplies each ingredient by its own eigenvalue,
    so the ingredient with the largest-magnitude eigenvalue outgrows all the
    others and eventually the mixture is nothing but that one direction.

    Normalising after each step is not part of the mathematics — it is there
    to stop the numbers running away. On this lab's matrix the vector would
    grow by a factor of 5 per step, so after 500 steps it would overflow to
    infinity and the direction, which is the thing we actually want, would be
    lost inside a NaN.

    Convergence is measured as the distance between successive UNIT vectors,
    after aligning their signs — because the iteration can flip the sign at
    every step when the dominant eigenvalue is negative, and a sign flip is
    not a failure to converge.

    Returns a dict with the vector, the eigenvalue estimate (the Rayleigh
    quotient), the iteration count, the final change, and the change at every
    step, so that the RATE of convergence can be inspected as well as the
    result.
    """
    m = np.asarray(matrix, dtype=float)
    if start is None:
        v = np.ones(m.shape[0], dtype=float)
    else:
        v = np.asarray(start, dtype=float).ravel().copy()
    norm = float(np.linalg.norm(v))
    if norm == 0.0:
        raise ValueError(
            "the power method cannot start from the zero vector: "
            "A @ 0 is 0 forever, which is why the zero vector is excluded "
            "from the definition of an eigenvector in the first place"
        )
    v = v / norm

    history: list[float] = []
    for iteration in range(1, max_iter + 1):
        w = m @ v
        length = float(np.linalg.norm(w))
        if length == 0.0:
            raise ValueError(
                "the iteration collapsed to the zero vector: the starting "
                "vector lay entirely in the part of space this matrix sends "
                "to the origin"
            )
        w = w / length
        if float(np.dot(w, v)) < 0.0:
            w = -w  # align signs so a flip is not mistaken for wandering
        change = float(np.linalg.norm(w - v))
        history.append(change)
        v = w
        if change < tol:
            return {
                "vector": v,
                "eigenvalue": rayleigh_quotient(m, v),
                "iterations": iteration,
                "change": change,
                "history": history,
                "converged": True,
            }

    return {
        "vector": v,
        "eigenvalue": rayleigh_quotient(m, v),
        "iterations": max_iter,
        "change": history[-1],
        "history": history,
        "converged": False,
    }


def rayleigh_quotient(matrix, v) -> float:
    """The best scalar estimate of the eigenvalue for a given vector.

    If v really is an eigenvector then A v = lambda v exactly, so

        (v . A v) / (v . v) = (v . lambda v) / (v . v) = lambda

    and if v is merely close to an eigenvector this is the closest thing to an
    eigenvalue that v admits.

    It is usually introduced with the claim that its error is the SQUARE of
    the vector's, so it converges twice as fast. That claim has a condition
    attached which is easy to lose: it needs the eigenvectors to be at right
    angles, which symmetry guarantees. Measured in exercise 4, the ratio of
    quotient error to squared angle locks onto 2.0 for the symmetric matrix
    in this lab and runs away for the non-symmetric one, where the quotient
    converges merely linearly and buys nothing over the vector.

    So: on a symmetric matrix — a covariance matrix, a Gram matrix, a graph
    Laplacian — take the Rayleigh quotient. On a general matrix, take it
    because it is one line and is never worse, but do not expect the speed-up.
    """
    m = np.asarray(matrix, dtype=float)
    v = np.asarray(v, dtype=float).ravel()
    return float(np.dot(v, m @ v) / np.dot(v, v))


# --------------------------------------------------------------------------
# Covariance, which is where PCA starts
# --------------------------------------------------------------------------


def covariance_matrix(data) -> np.ndarray:
    """The covariance matrix of an (n_points, n_features) array, from scratch.

    Two steps, and the first is the one people forget:

      1. Subtract the mean of each column, so the cloud is centred on the
         origin. Skip this and you measure how far the cloud is from the
         origin instead of how it is shaped, and the top eigenvector points
         at the cloud rather than along it.
      2. Take Xc.T @ Xc and divide by (n - 1).

    The (n - 1) rather than n is Bessel's correction, and it is what NumPy's
    numpy.cov uses by default. It does not change the eigenVECTORS at all —
    scaling a matrix by a constant scales its eigenvalues and leaves its
    eigenvectors exactly where they were — so for PCA's directions the choice
    is irrelevant. It matters only if you quote the eigenvalues as variances.

    The result is always symmetric, because entry (i, j) and entry (j, i) are
    the same sum of products written in the other order. That symmetry is not
    decoration: it is the guarantee that the eigenvalues are real and the
    eigenvectors are at right angles, which is what makes PCA well behaved.
    """
    x = np.asarray(data, dtype=float)
    if x.ndim != 2:
        raise ValueError(f"expected a 2-D (n_points, n_features) array, got shape {x.shape}")
    n_points = x.shape[0]
    if n_points < 2:
        raise ValueError("covariance needs at least two points")
    centred = x - x.mean(axis=0)
    return (centred.T @ centred) / (n_points - 1)


def principal_components(data) -> tuple[np.ndarray, np.ndarray]:
    """PCA in five lines: centre, covariance, eigh, sort, return.

    Returns (variances, directions) with the largest variance first, and
    directions given as COLUMNS so that directions[:, 0] is the top principal
    component — matching NumPy's own convention, which is worth matching
    exactly rather than improving on.

    numpy.linalg.eigh rather than numpy.linalg.eig, because a covariance
    matrix is symmetric and eigh is the routine written for that case: it
    returns real values in ascending order rather than unsorted complex ones,
    and on the authoring machine it was an order of magnitude faster on a
    400 by 400 matrix.
    """
    covariance = covariance_matrix(data)
    variances, directions = np.linalg.eigh(covariance)
    order = np.argsort(variances)[::-1]
    return variances[order], directions[:, order]

"""The reference suite: every claim this lab makes, checked against real values.

Run from the lab directory:

    .venv/bin/pytest examples -q -p no:cacheprovider

Every float comparison here names its tolerance in the assertion rather than
relying on a default, and every eigenvector comparison goes through
`abs_cosine` rather than comparing components — because an eigenvector is
defined only up to sign and scale, so a component-wise comparison fails on a
correct answer roughly half the time.
"""

from __future__ import annotations

import warnings

import numpy as np
import pytest

from dataset import (
    A,
    A_EIGEN_ANGLES_DEG,
    A_EIGENVALUES,
    A_EIGENVECTORS,
    CENTRE,
    ELONGATION_DEG,
    N_POINTS,
    PROJECTION_ONTO_X,
    REFLECTION_IN_X,
    ROTATION_60,
    ROTATION_90,
    SHEAR,
    SPREAD_ACROSS,
    SPREAD_ALONG,
    STANDARD_TRANSFORMATIONS,
    SYMMETRIC,
    SYMMETRIC_3X3,
    elongation_direction,
    make_cloud,
    power_method_start,
)
from eigen import (
    abs_cosine,
    characteristic_coefficients,
    covariance_matrix,
    deviation_degrees,
    direction_degrees,
    eigen_lines_by_sweep,
    eigenvalues_2x2,
    eigenvector_2x2,
    power_method,
    principal_components,
    rayleigh_quotient,
    solve_2x2,
    sweep_deviations,
)

#: Tolerances, named once so every test says which one it used.
EXACT = 1e-12
TIGHT = 1e-9
ANGLE = 1e-6
SWEEP = 1e-2  # a sampled sweep cannot beat its own grid spacing


# ==========================================================================
# The fan of vectors: which directions survive
# ==========================================================================


def test_only_two_directions_out_of_twentyfour_keep_their_line():
    kept = []
    for angle in range(0, 360, 15):
        radians = np.radians(angle)
        v = np.array([np.cos(radians), np.sin(radians)])
        if deviation_degrees(v, A @ v) < TIGHT:
            kept.append(angle)
    assert kept == [45, 225]


def test_45_and_225_are_the_same_line():
    up = np.array([np.cos(np.radians(45)), np.sin(np.radians(45))])
    down = np.array([np.cos(np.radians(225)), np.sin(np.radians(225))])
    assert abs_cosine(up, down) == pytest.approx(1.0, abs=EXACT)
    assert np.allclose(down, -up, atol=EXACT)


def test_a_full_sweep_finds_exactly_two_eigen_lines():
    found = eigen_lines_by_sweep(A)
    assert found["verdict"] == "some"
    assert len(found["lines"]) == 2
    assert found["lines"] == pytest.approx(list(A_EIGEN_ANGLES_DEG), abs=SWEEP)


def test_the_eigen_directions_deviate_by_nothing_at_all():
    for angle in A_EIGEN_ANGLES_DEG:
        radians = np.radians(angle)
        v = np.array([np.cos(radians), np.sin(radians)])
        assert deviation_degrees(v, A @ v) < TIGHT


def test_the_stretch_factors_are_the_eigenvalues():
    for angle, expected in zip(A_EIGEN_ANGLES_DEG, A_EIGENVALUES):
        radians = np.radians(angle)
        v = np.array([np.cos(radians), np.sin(radians)])
        stretch = float(np.linalg.norm(A @ v)) / float(np.linalg.norm(v))
        assert stretch == pytest.approx(expected, abs=TIGHT)


def test_a_non_eigen_direction_really_is_knocked_off_its_line():
    v = np.array([1.0, 0.0])
    assert deviation_degrees(v, A @ v) == pytest.approx(26.565051, abs=1e-6)


@pytest.mark.parametrize("vector,eigenvalue", list(zip(A_EIGENVECTORS, A_EIGENVALUES)))
def test_the_integer_eigenvectors_check_out_on_paper(vector, eigenvalue):
    v = np.array(vector)
    assert np.allclose(A @ v, eigenvalue * v, atol=EXACT)


def test_the_zero_vector_satisfies_the_equation_for_every_lambda():
    zero = np.zeros(2)
    for candidate in (5.0, 2.0, 0.0, -3.5, 1000.0):
        assert np.allclose(A @ zero, candidate * zero, atol=EXACT)


def test_abs_cosine_refuses_the_zero_vector():
    with pytest.raises(ValueError, match="no direction"):
        abs_cosine(np.zeros(2), np.array([1.0, 1.0]))


# ==========================================================================
# The hand solution
# ==========================================================================


def test_trace_and_determinant_are_the_characteristic_coefficients():
    trace, determinant = characteristic_coefficients(A)
    assert trace == 7.0
    assert determinant == 10.0


def test_the_characteristic_polynomial_really_vanishes_at_both_eigenvalues():
    trace, determinant = characteristic_coefficients(A)
    for eigenvalue in A_EIGENVALUES:
        value = eigenvalue**2 - trace * eigenvalue + determinant
        assert value == pytest.approx(0.0, abs=EXACT)


def test_hand_eigenvalues_match_the_worked_answer():
    values = eigenvalues_2x2(A)
    assert sorted(value.real for value in values) == pytest.approx([2.0, 5.0], abs=TIGHT)
    assert all(abs(value.imag) < EXACT for value in values)


def test_hand_eigenvalues_match_numpy_sorted_with_a_stated_tolerance():
    hand = sorted(value.real for value in eigenvalues_2x2(A))
    theirs = sorted(np.linalg.eig(A)[0].real)
    assert hand == pytest.approx(theirs, abs=TIGHT)


def test_a_minus_lambda_i_has_determinant_zero_at_each_eigenvalue():
    for eigenvalue in A_EIGENVALUES:
        shifted = A - eigenvalue * np.eye(2)
        assert float(np.linalg.det(shifted)) == pytest.approx(0.0, abs=EXACT)


@pytest.mark.parametrize("eigenvalue,expected", list(zip(A_EIGENVALUES, A_EIGENVECTORS)))
def test_hand_eigenvector_lies_on_the_right_line(eigenvalue, expected):
    v = eigenvector_2x2(A, eigenvalue)
    assert abs_cosine(v, expected) == pytest.approx(1.0, abs=TIGHT)


@pytest.mark.parametrize("eigenvalue", A_EIGENVALUES)
def test_hand_eigenvector_satisfies_the_defining_equation(eigenvalue):
    v = eigenvector_2x2(A, eigenvalue)
    assert np.allclose(A @ v, eigenvalue * v, atol=EXACT)


def test_hand_eigenvectors_come_back_as_unit_vectors():
    for eigenvalue in A_EIGENVALUES:
        assert float(np.linalg.norm(eigenvector_2x2(A, eigenvalue))) == pytest.approx(1.0, abs=EXACT)


def test_any_multiple_of_an_eigenvector_is_also_an_eigenvector():
    v = np.array(A_EIGENVECTORS[0])
    for scale in (-7.5, -1.0, 0.25, 3.0, 1000.0):
        scaled = scale * v
        assert np.allclose(A @ scaled, A_EIGENVALUES[0] * scaled, atol=1e-9)


def test_solve_2x2_returns_no_eigenvectors_for_a_rotation():
    values, vectors = solve_2x2(ROTATION_90)
    assert vectors is None
    assert all(abs(value.imag) > 0.5 for value in values)


def test_the_hand_method_rejects_a_non_2x2():
    with pytest.raises(ValueError, match="2x2"):
        characteristic_coefficients(np.eye(3))


# ==========================================================================
# What NumPy actually returns
# ==========================================================================


def test_numpy_eig_returns_complex_even_when_every_eigenvalue_is_real():
    """Observed on numpy 2.5.2, and it contradicts numpy's own docstring.

    The docstring for numpy.linalg.eig says the result "will be of complex
    type, unless the imaginary part is zero in which case it will be cast to a
    real type". On this version the imaginary part IS zero and the cast does
    NOT happen. The measurement is what this test records; if a future version
    changes the behaviour, this test going red is the correct outcome and the
    lesson text needs updating with it.
    """
    values, vectors = np.linalg.eig(A)
    assert values.dtype == np.complex128
    assert vectors.dtype == np.complex128
    assert np.all(values.imag == 0.0)


@pytest.mark.parametrize(
    "matrix",
    [np.eye(2), np.diag([1.0, 2.0, 3.0]), np.array([[2, 0], [0, 3]]), SYMMETRIC],
)
def test_eig_is_complex_for_every_real_eigenvalued_matrix_tried(matrix):
    assert np.linalg.eig(matrix)[0].dtype == np.complex128


def test_casting_those_eigenvalues_to_float_warns():
    values = np.linalg.eig(A)[0]
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        values.astype(float)
    assert [item.category.__name__ for item in caught] == ["ComplexWarning"]


def test_taking_real_without_checking_destroys_a_rotations_answer():
    values = np.linalg.eig(ROTATION_90)[0]
    assert np.allclose(np.abs(values), 1.0, atol=TIGHT)
    assert np.allclose(values.real, 0.0, atol=TIGHT)


def test_eigh_returns_real_sorted_values_on_symmetric_input():
    values = np.linalg.eigh(SYMMETRIC)[0]
    assert values.dtype == np.float64
    assert np.all(np.diff(values) >= 0)
    assert values == pytest.approx([1.0, 3.0], abs=TIGHT)


def test_eigh_silently_answers_a_different_question_on_non_symmetric_input():
    wrong = np.linalg.eigh(A)[0]
    lower_triangle_mirrored = np.array([[4.0, 2.0], [2.0, 3.0]])
    assert np.allclose(wrong, np.linalg.eigvalsh(lower_triangle_mirrored), atol=EXACT)
    assert not np.allclose(np.sort(wrong), sorted(A_EIGENVALUES), atol=1e-6)


def test_every_returned_pair_satisfies_a_v_equals_lambda_v():
    for matrix in [A, SYMMETRIC, SHEAR, REFLECTION_IN_X, PROJECTION_ONTO_X]:
        values, vectors = np.linalg.eig(matrix)
        for index in range(len(values)):
            left = matrix @ vectors[:, index]
            right = values[index] * vectors[:, index]
            assert float(np.abs(left - right).max()) < EXACT


def test_eigenvectors_from_numpy_are_unit_length():
    for matrix in [A, SYMMETRIC, SYMMETRIC_3X3]:
        vectors = np.linalg.eig(matrix)[1]
        assert np.allclose(np.linalg.norm(vectors, axis=0), 1.0, atol=EXACT)


def test_the_sign_ambiguity_is_real_and_component_comparison_fails():
    """The trap, pinned down. Do not 'fix' this test by flipping a sign."""
    values, vectors = np.linalg.eig(A)
    index = int(np.argmin(values.real))  # the lambda = 2 column
    mine = np.array(A_EIGENVECTORS[1], dtype=float)
    mine = mine / np.linalg.norm(mine)
    theirs = vectors.real[:, index]
    assert not np.allclose(mine, theirs, atol=1e-6)
    assert np.allclose(mine, -theirs, atol=TIGHT)
    assert abs_cosine(mine, theirs) == pytest.approx(1.0, abs=EXACT)


# ==========================================================================
# The standard transformations
# ==========================================================================


def test_a_shear_has_exactly_one_eigen_line():
    found = eigen_lines_by_sweep(SHEAR)
    assert len(found["lines"]) == 1
    line = found["lines"][0]
    assert min(line, 180.0 - line) < 0.05


def test_a_shear_returns_two_columns_that_are_the_same_direction():
    values, vectors = np.linalg.eig(SHEAR)
    assert np.allclose(values.real, [1.0, 1.0], atol=TIGHT)
    assert abs_cosine(vectors.real[:, 0], vectors.real[:, 1]) == pytest.approx(1.0, abs=1e-8)


@pytest.mark.parametrize("rotation", [ROTATION_90, ROTATION_60])
def test_a_plane_rotation_has_no_real_eigenvalues(rotation):
    values = np.linalg.eig(rotation)[0]
    assert np.all(np.abs(values.imag) > 0.5)
    assert np.allclose(np.abs(values), 1.0, atol=TIGHT)


@pytest.mark.parametrize("rotation", [ROTATION_90, ROTATION_60])
def test_a_plane_rotation_leaves_no_direction_on_its_own_line(rotation):
    found = eigen_lines_by_sweep(rotation)
    assert found["verdict"] == "none"
    assert found["lines"] == []


def test_the_rotations_negative_discriminant_is_the_geometry_speaking():
    trace, determinant = characteristic_coefficients(ROTATION_90)
    assert trace == 0.0
    assert determinant == pytest.approx(1.0, abs=EXACT)
    assert trace * trace - 4.0 * determinant < 0.0


@pytest.mark.parametrize("matrix", [np.eye(2), 2.0 * np.eye(2)])
def test_a_uniform_scaling_keeps_every_direction(matrix):
    found = eigen_lines_by_sweep(matrix)
    assert found["verdict"] == "every direction"
    assert found["fraction"] == pytest.approx(1.0, abs=EXACT)


def test_a_reflection_has_a_negative_eigenvalue_and_still_keeps_both_lines():
    values = np.sort(np.linalg.eig(REFLECTION_IN_X)[0].real)
    assert values == pytest.approx([-1.0, 1.0], abs=TIGHT)
    found = eigen_lines_by_sweep(REFLECTION_IN_X)
    assert len(found["lines"]) == 2


def test_a_projection_has_eigenvalue_zero_and_determinant_zero():
    values = np.sort(np.linalg.eig(PROJECTION_ONTO_X)[0].real)
    assert values == pytest.approx([0.0, 1.0], abs=TIGHT)
    assert float(np.linalg.det(PROJECTION_ONTO_X)) == pytest.approx(0.0, abs=EXACT)


def test_the_collapsed_direction_has_no_angle_to_measure():
    deviations, collapsed = sweep_deviations(PROJECTION_ONTO_X, [0.0, 45.0, 90.0])
    assert collapsed.tolist() == [False, False, True]
    assert np.isnan(deviations[2])
    assert deviations[0] == pytest.approx(0.0, abs=ANGLE)


@pytest.mark.parametrize("name", list(STANDARD_TRANSFORMATIONS))
def test_eigenvalues_multiply_to_the_determinant(name):
    matrix = STANDARD_TRANSFORMATIONS[name][0]
    product = complex(np.prod(np.linalg.eig(matrix)[0]))
    assert product.real == pytest.approx(float(np.linalg.det(matrix)), abs=TIGHT)
    assert abs(product.imag) < TIGHT


@pytest.mark.parametrize("name", list(STANDARD_TRANSFORMATIONS))
def test_eigenvalues_add_to_the_trace(name):
    matrix = STANDARD_TRANSFORMATIONS[name][0]
    total = complex(np.sum(np.linalg.eig(matrix)[0]))
    assert total.real == pytest.approx(float(np.trace(matrix)), abs=TIGHT)
    assert abs(total.imag) < TIGHT


@pytest.mark.parametrize("matrix", [SYMMETRIC, SYMMETRIC_3X3])
def test_a_symmetric_matrix_has_real_eigenvalues_and_orthogonal_eigenvectors(matrix):
    assert np.allclose(matrix, matrix.T, atol=EXACT)
    values, vectors = np.linalg.eigh(matrix)
    assert values.dtype == np.float64
    assert float(np.abs(vectors.T @ vectors - np.eye(len(values))).max()) < EXACT


def test_diagonalisation_reconstructs_the_original_matrix():
    """A = V D V-inverse: change basis, scale, change back."""
    values, vectors = np.linalg.eig(A)
    rebuilt = vectors @ np.diag(values) @ np.linalg.inv(vectors)
    assert np.allclose(rebuilt.real, A, atol=TIGHT)
    assert float(np.abs(rebuilt.imag).max()) < TIGHT


def test_the_shear_cannot_be_diagonalised_and_fails_silently():
    """The shear has one eigen-line, so its eigenvector matrix is singular and
    V D V-inverse cannot rebuild it.

    What is worth recording is HOW it fails. numpy.linalg.inv does not raise:
    the determinant is 2.2e-16 rather than exactly 0, so LAPACK inverts it and
    returns entries around 4.5e15. The reconstruction then comes back as the
    identity matrix — a clean, plausible, completely wrong answer, with no
    exception and no warning anywhere.

    The reliable check is the condition number, not an exception.
    """
    values, vectors = np.linalg.eig(SHEAR)
    vectors = vectors.real
    assert abs(float(np.linalg.det(vectors))) < 1e-15
    assert float(np.linalg.cond(vectors)) > 1e15

    inverse = np.linalg.inv(vectors)  # does NOT raise
    rebuilt = vectors @ np.diag(values.real) @ inverse
    assert not np.allclose(rebuilt, SHEAR, atol=1e-6)
    assert np.allclose(rebuilt, np.eye(2), atol=1e-9)


# ==========================================================================
# The power method
# ==========================================================================


def test_power_method_converges_in_25_iterations_to_the_stated_tolerance():
    result = power_method(A, power_method_start(), tol=1e-10)
    assert result["converged"] is True
    assert result["iterations"] == 25
    assert result["change"] < 1e-10


def test_power_method_finds_the_dominant_eigen_line():
    result = power_method(A, power_method_start(), tol=1e-10)
    assert abs_cosine(result["vector"], A_EIGENVECTORS[0]) == pytest.approx(1.0, abs=EXACT)
    assert direction_degrees(result["vector"]) == pytest.approx(45.0, abs=ANGLE)


def test_power_method_finds_the_dominant_eigenvalue():
    result = power_method(A, power_method_start(), tol=1e-10)
    assert result["eigenvalue"] == pytest.approx(5.0, abs=TIGHT)


def test_power_method_agrees_with_numpy():
    result = power_method(A, power_method_start(), tol=1e-10)
    values, vectors = np.linalg.eig(A)
    top = int(np.argmax(np.abs(values.real)))
    assert result["eigenvalue"] == pytest.approx(float(values.real[top]), abs=TIGHT)
    assert abs_cosine(result["vector"], vectors.real[:, top]) == pytest.approx(1.0, abs=EXACT)


def test_the_convergence_rate_equals_the_eigenvalue_ratio():
    history = power_method(A, power_method_start(), tol=1e-10)["history"]
    ratio = history[13] / history[12]
    assert ratio == pytest.approx(A_EIGENVALUES[1] / A_EIGENVALUES[0], abs=1e-3)


def test_close_eigenvalues_make_the_power_method_crawl():
    fast = power_method(A, power_method_start(), tol=1e-10)
    slow = power_method(np.diag([5.0, 4.9]), np.array([0.6, 0.8]), tol=1e-10)
    assert slow["iterations"] > 10 * fast["iterations"]


def test_power_method_works_from_several_different_starts():
    rng = np.random.default_rng(99)
    for _ in range(8):
        start = rng.normal(size=2)
        result = power_method(A, start, tol=1e-10)
        assert result["converged"] is True
        assert abs_cosine(result["vector"], A_EIGENVECTORS[0]) == pytest.approx(1.0, abs=1e-10)


def test_power_method_refuses_the_zero_vector():
    with pytest.raises(ValueError, match="zero vector"):
        power_method(A, np.zeros(2))


def test_power_method_reports_failure_rather_than_lying():
    result = power_method(A, power_method_start(), tol=1e-16, max_iter=5)
    assert result["converged"] is False
    assert result["iterations"] == 5


def test_the_rayleigh_quotient_is_exact_on_a_true_eigenvector():
    for vector, eigenvalue in zip(A_EIGENVECTORS, A_EIGENVALUES):
        assert rayleigh_quotient(A, vector) == pytest.approx(eigenvalue, abs=EXACT)


def _quotient_error_ratios(matrix, target, eigenvalue, steps=8):
    """(angle error, quotient error) after each of `steps` power iterations."""
    v = power_method_start()
    out = []
    for _ in range(steps):
        w = matrix @ v
        v = w / np.linalg.norm(w)
        angle = float(np.arccos(min(1.0, abs_cosine(v, target))))
        out.append((angle, abs(rayleigh_quotient(matrix, v) - eigenvalue)))
    return out


def test_the_rayleigh_quotient_is_quadratic_only_for_a_symmetric_matrix():
    """The textbook says the Rayleigh quotient converges twice as fast as the
    vector. Measured here, that is true for a SYMMETRIC matrix and false for
    this lab's non-symmetric one.

    On SYMMETRIC the quotient error divided by the squared angle settles on
    2.0 — dead-on quadratic. On A the quotient error divided by the angle
    itself settles on 1.0 — merely linear, the same rate as the vector.

    The quadratic result depends on the eigenvectors being orthogonal, which
    symmetry guarantees and A does not have.
    """
    symmetric = _quotient_error_ratios(SYMMETRIC, np.array([1.0, 1.0]), 3.0)
    angle, quotient = symmetric[-1]
    assert quotient / angle**2 == pytest.approx(2.0, abs=1e-3)

    unsymmetric = _quotient_error_ratios(A, np.array(A_EIGENVECTORS[0]), 5.0)
    angle, quotient = unsymmetric[-1]
    assert quotient / angle == pytest.approx(1.0, abs=1e-2)
    assert quotient / angle**2 > 100.0  # nowhere near quadratic


def test_unnormalised_iteration_overflows_to_infinity():
    v = power_method_start()
    with np.errstate(over="ignore"):
        for _ in range(600):
            v = A @ v
    assert not np.all(np.isfinite(v))


# ==========================================================================
# PCA
# ==========================================================================


def test_the_cloud_is_reproducible_from_the_seed():
    first = make_cloud()
    second = make_cloud()
    assert np.array_equal(first, second)
    assert first.shape == (N_POINTS, 2)


def test_the_cloud_sits_where_it_was_built_to_sit():
    assert make_cloud().mean(axis=0) == pytest.approx(CENTRE, abs=0.2)


def test_covariance_from_scratch_matches_numpy_cov():
    cloud = make_cloud()
    assert np.allclose(covariance_matrix(cloud), np.cov(cloud, rowvar=False), atol=EXACT)


def test_a_covariance_matrix_is_always_symmetric():
    covariance = covariance_matrix(make_cloud())
    assert np.allclose(covariance, covariance.T, atol=1e-15)


def test_covariance_refuses_a_1d_array():
    with pytest.raises(ValueError, match="2-D"):
        covariance_matrix(np.array([1.0, 2.0, 3.0]))


def test_the_top_principal_component_lies_along_the_known_elongation():
    _variances, directions = principal_components(make_cloud())
    similarity = abs_cosine(directions[:, 0], elongation_direction())
    assert similarity > 0.999
    assert similarity == pytest.approx(0.9999984422, abs=1e-9)


def test_the_top_component_is_within_a_fifth_of_a_degree_of_the_truth():
    _variances, directions = principal_components(make_cloud())
    error = abs(direction_degrees(directions[:, 0]) - ELONGATION_DEG)
    assert error < 0.2


def test_the_returned_top_component_points_the_other_way_along_that_line():
    """The sign ambiguity, in the artifact that matters. Not a defect."""
    _variances, directions = principal_components(make_cloud())
    truth = elongation_direction()
    assert float(np.dot(directions[:, 0], truth)) < 0.0
    assert not np.allclose(directions[:, 0], truth, atol=1e-3)
    assert abs_cosine(directions[:, 0], truth) > 0.999


def test_the_eigenvalues_recover_the_spreads_the_cloud_was_built_with():
    variances, _directions = principal_components(make_cloud())
    assert float(np.sqrt(variances[0])) == pytest.approx(SPREAD_ALONG, abs=0.2)
    assert float(np.sqrt(variances[1])) == pytest.approx(SPREAD_ACROSS, abs=0.1)


def test_the_first_component_carries_almost_all_the_variance():
    variances, _directions = principal_components(make_cloud())
    assert float(variances[0] / variances.sum()) > 0.97


def test_the_components_come_back_largest_first():
    variances, _directions = principal_components(make_cloud())
    assert np.all(np.diff(variances) <= 0)


def test_the_two_principal_components_are_perpendicular():
    _variances, directions = principal_components(make_cloud())
    assert float(np.dot(directions[:, 0], directions[:, 1])) == pytest.approx(0.0, abs=EXACT)


def test_projections_onto_the_components_are_uncorrelated():
    _variances, directions = principal_components(make_cloud())
    cloud = make_cloud()
    projected = (cloud - cloud.mean(axis=0)) @ directions
    assert abs(float(np.corrcoef(projected.T)[0, 1])) < EXACT


def test_skipping_the_centring_gives_a_confidently_wrong_answer():
    cloud = make_cloud()
    uncentred = (cloud.T @ cloud) / (N_POINTS - 1)
    values, vectors = np.linalg.eigh(uncentred)
    top = vectors[:, int(np.argmax(values))]
    assert abs(direction_degrees(top) - ELONGATION_DEG) > 5.0


def test_eig_and_eigh_agree_on_the_covariance_matrix():
    covariance = covariance_matrix(make_cloud())
    assert np.allclose(
        np.sort(np.linalg.eig(covariance)[0].real),
        np.sort(np.linalg.eigh(covariance)[0]),
        atol=EXACT,
    )

lesson_id: D106
day: 106
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-106-eigenvalues-and-eigenvectors-intuitively
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - .venv/bin/python3 -c "import numpy; print(numpy.__version__)"
run_commands:
  - 'cd examples && ../.venv/bin/python3 01_the_fan_of_vectors.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_by_hand_2x2.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_standard_transformations.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_power_method.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_pca_from_covariance.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_eig_against_eigh.py && cd ..'
  - .venv/bin/pytest examples -q -p no:cacheprovider
  - .venv/bin/pytest starter -q -p no:cacheprovider
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - "find . -type d -name '__pycache__' -prune -exec rm -rf -- {} +"
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 35
last_executed: '2026-08-17'
executed_on: 'macOS 26.5.2 (Apple Silicon, arm64), Python 3.14.0, numpy 2.5.2, pytest 9.1.1, bash 3.2.57 — bash tests/run_tests.sh -> 110 checks, 0 failure(s), exit 0; pytest examples -> 94 passed; pytest starter -> 1 passed, 52 skipped on an untouched checkout. All six reference scripts exit 0 with every internal assertion holding. The harness was additionally run with the lab-local .venv removed entirely and PYTEST pointed at an external interpreter: 110 checks, 0 failure(s), exit 0 in that configuration too. Section 7 prunes .venv before every find, which is load-bearing rather than cosmetic — the installed .venv contains 113 __pycache__ directories and 6 data files, so without the prune the lab would fail the reader for following its own installation instructions. Network is needed once to install numpy and pytest; nothing else in the lab opens a socket, and the 400-point dataset is GENERATED from numpy.random.default_rng(2106) rather than downloaded — section 7 greps the sources for network calls and also asserts that no data file exists anywhere under the lab. Section 6 re-runs the harness with one expectation deliberately swapped for the naive belief that a shear has two eigendirections because numpy.linalg.eig returns two columns, and asserts that the re-run exits non-zero and reports exactly one failure, so the suite is demonstrated to be capable of failing rather than merely claimed to be. Four measured results are asserted rather than smoothed over. (1) numpy.linalg.eig returns complex128 for a real matrix whose eigenvalues are both real, which CONTRADICTS the docstring shipped with this very version — it says the array "will be of complex type, unless the imaginary part is zero in which case it will be cast to a real type", and on numpy 2.5.2 the imaginary part is zero and the cast does not happen, for A, for numpy.eye(2), for numpy.diag([1., 2., 3.]) and for an integer matrix. The measurement is recorded; if a future version casts, that test going red is the correct outcome. (2) A shear has ONE eigen-line, confirmed independently by algebra and by a 180,000-direction brute-force sweep, while eig returns TWO columns whose absolute cosine is 1.0 — geometric multiplicity 1 against algebraic multiplicity 2. Attempting V D V-inverse on it does not raise: numpy.linalg.inv inverts a matrix of determinant 2.2e-16 and returns entries around 4.5e15, and the reconstruction comes back as a clean, plausible, completely wrong identity matrix with no exception and no warning. The reliable check is the condition number. (3) The textbook claim that the Rayleigh quotient converges quadratically is TRUE for the symmetric matrix here (quotient error over squared angle locks onto 2.0000) and FALSE for this lab non-symmetric A (quotient error over angle settles on 1.0, merely linear). The condition the textbook attaches and that is easy to drop is orthogonal eigenvectors, which symmetry guarantees and A does not have — its eigen-lines meet at 71.5651 degrees, not 90. (4) PCA recovers 30.101134 degrees from 400 points built along 30.0 degrees, and returns the component pointing the OPPOSITE way along that line, so numpy.allclose says False on an answer that is exactly right (absolute cosine 0.9999984422). That sign and scale ambiguity is treated as the central teaching point rather than papered over: every comparison in the lab goes through absolute cosine, and expected-output/FIELDS.md states that a flipped sign in a reader output is not a difference.'

# Dependencies for the Day 106 lab

Two packages. Both free and open source, both installed from the Python Package
Index with `pip`, both running entirely on your own machine.

| Package | Pinned version | Why this lab needs it |
| --- | --- | --- |
| `numpy` | `2.5.2` | Holds the vectors and matrices, and supplies the independent answer every hand calculation is checked against: `numpy.linalg.eig`, `numpy.linalg.eigh`, `numpy.linalg.eigvals`, `numpy.linalg.eigvalsh`, `numpy.linalg.det`, `numpy.cov` and the seeded generator behind `numpy.random.default_rng`. |
| `pytest` | `9.1.1` | The test runner from Days 071–074. Nothing new here except what it is pointed at. |

That is the whole list. This lab needs no plotting library, no image library and
no machine-learning library, and that is deliberate — see below.

## Why the from-scratch code does not call NumPy's eigensolvers

`examples/eigen.py` and `starter/eigen.py` compute eigenvalues from the
characteristic equation with the school quadratic formula, find eigenvectors by
reading a row of `A - lambda*I`, and find the dominant eigenvector by repeated
multiplication. None of those functions calls `numpy.linalg.eig`.

If they did, checking them against `numpy.linalg.eig` would be checking NumPy
against itself, and would prove exactly nothing.

NumPy holds the **arrays** and does the **arithmetic** — dot products, norms,
matrix-vector products — because writing those by hand teaches nothing you did
not learn on Day 104. The **eigen-mathematics** is ours. That split is what
makes `examples/06_eig_against_eigh.py` mean something: two implementations
that share no eigen-code, agreeing on this lab's matrix to within 4.5e-11 on
the eigenvalue and to fifteen digits on the direction.

## Why `scikit-learn` is not here, even though this lab does PCA

Exercise 5 builds PCA from a covariance matrix and its eigenvectors in about
fifteen lines. `sklearn.decomposition.PCA` does the same job better — it uses a
singular value decomposition rather than an explicit covariance matrix, which
is more accurate on ill-conditioned data, and it hands you the centring, the
sorting, the variance ratios and a `transform` method for free.

**Use scikit-learn for real work.** The fifteen lines here exist so that when
you later call `PCA(n_components=50)` you know precisely what it did and why the
answer sometimes comes back with the sign flipped. Installing it would let you
run the tool without ever seeing the mechanism, which is the opposite of the
point.

`examples/06_eig_against_eigh.py` describes scikit-learn, SciPy and PyTorch from
their own documentation and reproduces **no output from any of them**, because
none of them is installed here.

## Why the versions are pinned

They are *checked* rather than assumed. Section 1 of `tests/run_tests.sh` reads
the installed versions and compares them against this file, so a mismatch is
reported at the top of the run rather than surfacing later as a confusing diff.

One place the version genuinely matters, and it is handled by measurement
rather than by trusting the pin:

**`numpy.linalg.eig` returns `complex128` on a real matrix with real
eigenvalues.** The docstring shipped with numpy 2.5.2 says the result "will be
of complex type, unless the imaginary part is zero in which case it will be
cast to a real type". On this version, on the authoring machine, the imaginary
part *is* zero and the cast does *not* happen — for `A`, for `numpy.eye(2)`,
and for every other real-eigenvalued matrix tried.

The lab measures that every run instead of asserting it from memory.
`test_numpy_eig_returns_complex_even_when_every_eigenvalue_is_real` will fail
loudly if a future version changes the behaviour, and that failure is the
**correct** outcome: it means the lesson text needs updating, not the test.

## Licences

NumPy is distributed under the BSD 3-Clause licence and pytest under the MIT
licence, each stated on that project's own documentation site. Both are
maintained in the open, cost nothing, and need no account, no key and no
signup — personally or commercially.

## One-time install

From the lab directory:

```bash
python3 -m venv .venv
.venv/bin/pip install -r requirements/requirements.txt
.venv/bin/python3 -c "import numpy; print(numpy.__version__)"
```

Expect `2.5.2`. Day 43 covered `python3 -m venv` in full; this is the same
pattern. The environment lives in `.venv/` inside the lab, is already excluded
from version control, and can be deleted at any time with `rm -rf .venv`.

## Network

Installing needs the network, once. **Nothing else in this lab does.**

In particular, the dataset is *generated in code*, not downloaded. A 400-point
cloud drawn from a seeded generator is reproducible on every machine, needs no
licence check, and cannot break on a train. `examples/dataset.py` builds it
from `numpy.random.default_rng(2106)` and every number in `expected-output/`
follows from that seed. Section 7 of `tests/run_tests.sh` greps every file under
`examples/` and `starter/` for the patterns that would indicate a socket being
opened.

## Running without a lab-local environment

If NumPy and pytest are already available in an environment you have activated,
the harness will find `pytest` on your `PATH`. You can also point it at a
specific binary:

```bash
PYTEST=/path/to/pytest bash tests/run_tests.sh
```

The harness uses the `python3` that sits beside that `pytest`, because that is
the interpreter the packages are installed into. If NumPy is not importable
from it, the harness says so and stops rather than skipping checks quietly.

## What you would give up with an older NumPy

Very little, and the lab will tell you rather than guess.

`numpy.random.default_rng` needs NumPy 1.17 or later, `numpy.emath.sqrt` has
been there far longer, and every eigensolver used here is ancient by library
standards. Section 1 of the harness checks only that NumPy's major version is 2
or later.

What would change on an older version is the seeded cloud. `default_rng` is
guaranteed reproducible for a given NumPy generation, not across all of them,
so on NumPy 1.x the 400 points would differ and the PCA numbers in
`expected-output/05-pca-from-covariance.txt` would not match to the last digit.
The *claims* would still hold — the top component still lands within a fifth of
a degree of 30 — because the lab asserts tolerances around the construction
rather than the exact digits. `expected-output/FIELDS.md` says which is which.

numpy==2.5.2
pytest==9.1.1

# The Vectors That Keep Their Direction — your brief

Five exercises, in order. Exercise 1 is code; exercises 2 to 5 are predictions
you write down before you run anything.

Everything you write goes in exactly two files:

- `starter/eigen.py` — six functions, each currently `return NotImplemented`
- `starter/answers.py` — twenty-six predictions, each currently `None`

Nothing else in `starter/` needs editing, and `starter/dataset.py` is read-only
data you should read rather than change.

## Before you start

From the **lab directory** (the one with `README.md` in it), not from
`starter/`:

```bash
python3 -m venv .venv
.venv/bin/pip install -r requirements/requirements.txt
.venv/bin/pytest starter -q
```

That last command should print:

```
1 passed, 52 skipped
```

The one pass is `test_the_environment_is_ready`. The fifty-two skips are the
exercises. **A skip means "not attempted", not "broken".** When it says
`53 passed`, you are finished.

Run that command as often as you like. A failure prints your answer beside the
real one, so a wrong guess still teaches you something.

## The one habit that matters more than any other

An eigenvector is defined **only up to sign and scale**.

If `A v = lambda v`, then the same equation holds for `-v`, for `3.7 v`, and
for every other non-zero multiple. There is no such thing as *the* eigenvector
for an eigenvalue — there is an eigen-**line**, and every library that hands
you one vector has made an arbitrary choice on your behalf.

So `(1, -2)`, `(-1, 2)` and `(0.447, -0.894)` are all the same answer, and
`numpy.allclose` will tell you two of them are wrong.

Every test in this lab compares **directions**, using the absolute cosine.
Write your code the same way and this trap never catches you. Fight it — try to
"fix" a sign somewhere — and you will spend an hour on a bug that was never
there. Exercise 5f is built to spring this trap on purpose.

---

## Exercise 1 — write the six functions (`starter/eigen.py`)

Each function's docstring gives the steps and, where it matters, says why the
obvious implementation is wrong. Read the docstring before writing the body.

| Function | What it does | Watch out for |
| --- | --- | --- |
| `abs_cosine(u, v)` | 1.0 when `u` and `v` lie on the same line | Take the **absolute** value. Raise `ValueError` mentioning "no direction" on a zero vector. |
| `characteristic_coefficients(matrix)` | `(trace, determinant)` of a 2x2 | Raise `ValueError` mentioning "2x2" on anything else. |
| `eigenvalues_2x2(matrix)` | Solve the characteristic equation | Use `numpy.emath.sqrt`, **not** `numpy.sqrt`. Return complex always. |
| `eigenvector_2x2(matrix, eigenvalue)` | A unit eigenvector for a real eigenvalue | Try both rows. Do not try to fix the sign. |
| `power_method(matrix, start, ...)` | The dominant eigenvector by repeated multiplication | Normalise every round. Align the signs. Report non-convergence rather than raising. |
| `covariance_matrix(data)` | Covariance of an `(n_points, n_features)` array | **Centre it first.** Divide by `n - 1`. |

Check your progress at any point:

```bash
.venv/bin/pytest starter -q -k "1a or 1b or 1c or 1d or 1e or 1f"
```

Two of these have a failure mode that produces a plausible wrong answer rather
than an error, which is why they get their own warning here:

- **`eigenvalues_2x2` with `numpy.sqrt`.** On a rotation the discriminant is
  negative. `numpy.sqrt(-4.0)` returns `nan` and emits a `RuntimeWarning`;
  `numpy.emath.sqrt(-4.0)` returns `2j`, which is the actual answer. Test 1c
  checks the rotation case for exactly this reason.
- **`power_method` without the sign alignment.** Give it a matrix whose
  dominant eigenvalue is *negative* and the iterate flips direction every
  single step. The answer converged on round three; your `change` measurement
  never drops below `tol` and the loop runs to `max_iter` reporting failure.
  Test 1e uses `numpy.diag([-5.0, 2.0])` to catch this.

## Exercise 2 — solve the 2x2 by hand (`starter/answers.py`, 2a–2g)

```
A = [[4, 1],
     [2, 3]]
```

Pencil and paper. Seven answers, none of them needing a calculator:

1. The trace, then the determinant.
2. The characteristic equation is `lambda^2 - b*lambda + c = 0`. Those two
   numbers are the trace and the determinant, in that order — derive that once
   from `det(A - lambda*I)` and you never have to again.
3. The discriminant `b^2 - 4c`. Its **sign** is the interesting part.
4. Both eigenvalues, largest first. The quadratic factorises over the integers.
5. An eigenvector for each. For `lambda = 5`, the first row of `A - 5I` is
   `[-1, 1]`, which says `-x + y = 0`. For `lambda = 2`, the first row of
   `A - 2I` is `[2, 1]`.

Any non-zero multiple of the right eigenvector passes. Check yourself before
running: `A @ (1, 1)` should come out as exactly `5 * (1, 1)`.

## Exercise 3 — the standard transformations (`answers.py`, 3a–3g)

The matrices from Day 102, and what each one does to directions. Predict from
the **geometry** first, then run `examples/03_standard_transformations.py` and
see whether you were right.

- **3a** — the shear `[[1, 1], [0, 1]]`. Count eigen-**lines**, not columns
  returned by NumPy. Those are different numbers here, and the difference is
  the question.
- **3b, 3c** — the 90-degree rotation. Picture an arrow being turned before you
  reach for any algebra. Then remember that a rotation changes no lengths.
- **3d** — what dtype does `numpy.linalg.eig` return for `A`, whose eigenvalues
  are 5 and 2, both real? Predict from the documentation, then measure. If your
  prediction and the machine disagree, **the machine is right**, and that
  disagreement is one of the things this lab exists to show you.
- **3e** — the projection `[[1, 0], [0, 0]]`. Its determinant and its smaller
  eigenvalue are the same number, and Day 102 explains why.
- **3f** — the reflection. One eigenvalue is negative. Say out loud what a
  negative eigenvalue means before you write it down.
- **3g** — for a symmetric matrix, at what angle do the eigenvectors meet?

## Exercise 4 — the power method (`answers.py`, 4a–4e)

Multiply, normalise, repeat. Predict what it converges to and how fast.

- **4a, 4b** — which eigenvalue and which direction. The clue is in the name
  "dominant".
- **4c** — each round the error shrinks by a constant factor. It is a ratio of
  the two eigenvalues, and the **smaller one is on top**. Reason about which
  ingredient of the mixture is dying out relative to which.
- **4d** — eigenvalues of 5 and 4.9 instead of 5 and 2: more iterations or
  fewer?
- **4e** — what happens to the raw length of `A^k v0` if you never normalise?

`examples/04_power_method.py` shows every one of these happening, with the
iteration table printed step by step. Predict first, then read it.

## Exercise 5 — PCA (`answers.py`, 5a–5g)

The cloud in `dataset.py` is 400 points deliberately stretched along **30
degrees**, with a standard deviation of 3.0 along that direction and 0.4 across
it, centred at `(5, -2)`. That direction appears nowhere in the array. PCA has
to rediscover it from 400 pairs of coordinates.

- **5a, 5b** — the shape of the covariance matrix of a `(400, 2)` dataset, and
  whether it is symmetric. It is neither `(400, 400)` nor `(400, 2)`.
- **5c, 5d** — where the top eigenvector points, and what the square root of
  the top eigenvalue comes out near.
- **5e** — `eig` or `eigh` for a covariance matrix?
- **5f** — the sign trap, sprung. The top component is **correct** and
  `numpy.allclose` against the true direction returns... what?
- **5g** — forget to subtract the mean. Is the answer still within 5 degrees?

## When you are done

```bash
.venv/bin/pytest starter -q          # expect: 53 passed
bash tests/run_tests.sh              # expect: 0 failure(s), exit 0
```

Then read `examples/` end to end. Every function you wrote has a reference
version there with a docstring explaining the choices, and the six numbered
scripts print the whole story with real numbers. Reading them *after* writing
your own is worth several times reading them before.

"""Exercises 2 to 5 — predictions. Work them out, then let the tests check you.

Every answer below is `None`, which makes its test SKIP. Replace a None with
your answer and the skip becomes a pass or a failure — and a failure prints
both your answer and the real one, so a wrong guess still teaches you
something.

The rule for this file: PREDICT FIRST, then run. Every one of these can be
reasoned out from the lesson. Looking them up by running NumPy first turns a
thinking exercise into a typing exercise.
"""

# ==========================================================================
# EXERCISE 2 — solving the 2x2 by hand
#
#     A = [[4, 1],
#          [2, 3]]
# ==========================================================================

#: 2a. The trace of A (the sum of the diagonal entries). An integer.
TRACE_OF_A = None

#: 2b. The determinant of A. An integer.
DETERMINANT_OF_A = None

#: 2c. The characteristic equation is lambda^2 - b*lambda + c = 0.
#:     Give (b, c) as a tuple of two integers.
CHARACTERISTIC_COEFFICIENTS = None

#: 2d. The discriminant, b^2 - 4c. An integer.
#:     Its SIGN is the interesting part: positive means two real eigenvalues,
#:     zero means one repeated, negative means none that are real.
DISCRIMINANT = None

#: 2e. The two eigenvalues, as a tuple of two integers, LARGEST FIRST.
EIGENVALUES_LARGEST_FIRST = None

#: 2f. An eigenvector for the LARGER eigenvalue, as a tuple of two integers
#:     with no common factor. There are infinitely many correct answers that
#:     differ by scale and sign; the test compares directions, so any
#:     non-zero multiple of the right one passes.
#:     Hint: solve (A - 5I) v = 0. The first row of A - 5I is [-1, 1].
EIGENVECTOR_FOR_LARGER = None

#: 2g. An eigenvector for the SMALLER eigenvalue, same rules.
#:     Hint: the first row of A - 2I is [2, 1].
EIGENVECTOR_FOR_SMALLER = None


# ==========================================================================
# EXERCISE 3 — the standard transformations
# ==========================================================================

#: 3a. How many DISTINCT eigen-lines does the shear [[1, 1], [0, 1]] have?
#:     An integer. Count lines, not columns returned by NumPy — those are not
#:     the same number here, and the difference is the point of the question.
SHEAR_EIGEN_LINE_COUNT = None

#: 3b. Are the eigenvalues of the 90-degree rotation [[0, -1], [1, 0]] real?
#:     True or False. Picture what a rotation does to an arrow before you
#:     reach for the algebra.
ROTATION_EIGENVALUES_ARE_REAL = None

#: 3c. The MAGNITUDE (absolute value) of each eigenvalue of that rotation.
#:     A single number, since both are the same. Hint: a rotation changes no
#:     lengths at all.
ROTATION_EIGENVALUE_MAGNITUDE = None

#: 3d. numpy.linalg.eig is given the real matrix A, whose eigenvalues are 5
#:     and 2 — both real. What dtype does it return the eigenvalues in?
#:     The string 'float64' or the string 'complex128'.
#:     Predict from the documentation, then run it. If your prediction and the
#:     machine disagree, the machine is right and that disagreement is one of
#:     the things this lab exists to show you.
EIG_DTYPE_ON_A = None

#: 3e. The determinant of the projection [[1, 0], [0, 0]], and the smaller of
#:     its two eigenvalues, as a tuple of two numbers.
#:     They are the same number, and that is not a coincidence.
PROJECTION_DET_AND_SMALLEST_EIGENVALUE = None

#: 3f. The eigenvalues of the reflection [[1, 0], [0, -1]], as a tuple of two
#:     integers, largest first. One of them is negative — say what a negative
#:     eigenvalue means to yourself before you write it down.
REFLECTION_EIGENVALUES = None

#: 3g. For a symmetric matrix, at what angle (in degrees) do the eigenvectors
#:     meet? A single number.
SYMMETRIC_EIGENVECTOR_ANGLE_DEG = None


# ==========================================================================
# EXERCISE 4 — the power method on A, started from a seeded random vector
# ==========================================================================

#: 4a. Which eigenvalue does the power method converge towards? A number.
POWER_METHOD_FINDS_EIGENVALUE = None

#: 4b. The direction it converges to, in degrees in [0, 180). A number.
#:     Hint: it is the direction of the eigenvector for 4a.
POWER_METHOD_FINDS_DIRECTION_DEG = None

#: 4c. Each iteration, the remaining error shrinks by roughly a constant
#:     factor. What factor? A number between 0 and 1.
#:     Hint: it is a ratio of the two eigenvalues, and the smaller one is on
#:     top. Reason about which ingredient is dying out relative to which.
CONVERGENCE_RATIO = None

#: 4d. If the two eigenvalues were 5 and 4.9 instead of 5 and 2, would the
#:     power method need MORE or FEWER iterations to reach the same
#:     tolerance? The string 'more' or the string 'fewer'.
CLOSE_EIGENVALUES_NEED = None

#: 4e. If you never normalise, what does the length of A^k v0 do as k grows
#:     past a few hundred? One of the strings 'overflows to inf',
#:     'shrinks to zero', 'stays at 1'.
UNNORMALISED_LENGTH_BEHAVIOUR = None


# ==========================================================================
# EXERCISE 5 — PCA on the invented cloud
#
# The cloud is 400 points, deliberately stretched along 30 degrees, with a
# standard deviation of 3.0 along that direction and 0.4 across it, centred
# at (5, -2).
# ==========================================================================

#: 5a. What is the shape of the covariance matrix of a (400, 2) dataset?
#:     A tuple. It is NOT (400, 400) and it is NOT (400, 2).
COVARIANCE_SHAPE = None

#: 5b. Is the covariance matrix symmetric? True or False.
COVARIANCE_IS_SYMMETRIC = None

#: 5c. Roughly what direction, in degrees, does the top eigenvector point
#:     along? A number. The test allows one degree of slack, because 400
#:     samples estimate a direction rather than reproduce it.
TOP_COMPONENT_DIRECTION_DEG = None

#: 5d. The square root of the LARGEST eigenvalue should come out close to one
#:     of the numbers the cloud was built with. Which one? A number.
#:     Half a unit of slack is allowed.
SQRT_OF_TOP_EIGENVALUE = None

#: 5e. Which NumPy routine is the right one for a covariance matrix?
#:     The string 'eig' or the string 'eigh'.
RIGHT_ROUTINE_FOR_COVARIANCE = None

#: 5f. numpy.allclose is used to compare the top eigenvector against the true
#:     direction (0.866, 0.5) that the cloud was built along. The answer is
#:     correct. Does numpy.allclose return True or False?
#:     Think about what is and is not determined about an eigenvector.
ALLCLOSE_ON_THE_CORRECT_COMPONENT = None

#: 5g. If you forget to subtract the mean before computing the covariance,
#:     is the top eigenvector still within 5 degrees of the truth?
#:     True or False.
UNCENTRED_STILL_CORRECT = None

"""Make this directory's own eigen.py the one its tests import.

Both `examples/` and `starter/` contain modules called `eigen` and `dataset`,
and pytest imports test files by putting their directory on `sys.path`.
Without this file, running a bare `pytest` across both directories at once
would import whichever `eigen` was seen first and then reuse it for the other
suite — so these starter tests would silently pass against the reference
solution instead of skipping. That is a wrong answer with a green tick on it,
which is the worst kind.

So: put this directory first on the import path, and drop any already-imported
`eigen` or `dataset` that came from somewhere else.

Section 4 of `tests/run_tests.sh` checks that this still works, by comparing
the skip count from `pytest starter` against the skip count from a bare
`pytest` over the whole lab. If this guard ever stops working, that check goes
red rather than the lab quietly lying to you.
"""

import sys
from pathlib import Path

HERE = str(Path(__file__).parent.resolve())

if HERE in sys.path:
    sys.path.remove(HERE)
sys.path.insert(0, HERE)

for name in ("eigen", "dataset"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""The data for the exercises. Read this file; you do not need to change it.

The star of the lab is `A`, chosen so the whole eigenvalue calculation can be
done with a pencil in about a minute:

    A = [[4, 1],
         [2, 3]]

Work it out before you write any code. You will need the answer in exercise 2.
"""

from __future__ import annotations

import numpy as np

A = np.array([[4.0, 1.0], [2.0, 3.0]])

SHEAR = np.array([[1.0, 1.0], [0.0, 1.0]])
ROTATION_90 = np.array([[0.0, -1.0], [1.0, 0.0]])
ROTATION_60 = np.array(
    [
        [0.5, -np.sqrt(3.0) / 2.0],
        [np.sqrt(3.0) / 2.0, 0.5],
    ]
)
SYMMETRIC = np.array([[2.0, 1.0], [1.0, 2.0]])
PROJECTION_ONTO_X = np.array([[1.0, 0.0], [0.0, 0.0]])
REFLECTION_IN_X = np.array([[1.0, 0.0], [0.0, -1.0]])

# --------------------------------------------------------------------------
# The invented 2-D dataset for exercise 5
# --------------------------------------------------------------------------

SEED = 2106
ELONGATION_DEG = 30.0
SPREAD_ALONG = 3.0
SPREAD_ACROSS = 0.4
N_POINTS = 400
CENTRE = np.array([5.0, -2.0])


def elongation_direction() -> np.ndarray:
    """The unit vector the cloud is stretched along — the answer PCA must find."""
    radians = np.radians(ELONGATION_DEG)
    return np.array([np.cos(radians), np.sin(radians)])


def make_cloud() -> np.ndarray:
    """400 points, shape (400, 2), stretched along ELONGATION_DEG degrees.

    Seeded, so your numbers will match the captured output exactly.
    """
    rng = np.random.default_rng(SEED)
    along = elongation_direction()
    across = np.array([-along[1], along[0]])
    travel_along = rng.normal(0.0, SPREAD_ALONG, size=N_POINTS)
    travel_across = rng.normal(0.0, SPREAD_ACROSS, size=N_POINTS)
    return travel_along[:, None] * along + travel_across[:, None] * across + CENTRE


def power_method_start() -> np.ndarray:
    """A seeded random unit vector, for exercise 4."""
    rng = np.random.default_rng(106)
    v = rng.normal(size=2)
    return v / np.linalg.norm(v)

"""Exercise 1 — six functions to write. Your work goes here.

Every function currently `return NotImplemented`, which is what makes the
matching test SKIP rather than fail. Replace each one and the skip turns into
a pass. Check yourself at any point with, from the lab directory:

    .venv/bin/pytest starter -q

Read `00_brief.md` first. It gives the exercises in order with the exact
commands.

The docstrings below tell you what each function must do and, where it
matters, WHY the obvious implementation is wrong. Read them before writing.
"""

from __future__ import annotations

import numpy as np


def abs_cosine(u, v):
    """EXERCISE 1a — return the absolute cosine of the angle between u and v.

    1.0 when the two lie on the same LINE, 0.0 when they are at right angles.

    Steps:
      1. Convert both to float arrays with numpy.asarray(..., dtype=float)
         and flatten them with .ravel().
      2. Compute each one's length with numpy.linalg.norm.
      3. If either length is 0.0, raise ValueError with a message containing
         the words "no direction" — the zero vector has no direction, so no
         angle exists. Returning 0.0 or nan instead would be a lie.
      4. Otherwise return abs(numpy.dot(u, v)) divided by the product of the
         two lengths, as a plain Python float.

    Take the ABSOLUTE value. That is the whole point of this function and the
    single most important habit in this lab: an eigenvector is defined only up
    to sign and scale, so (1, -2) and (-1, 2) are the same answer. Without the
    abs, half of your correct answers will look wrong.
    """
    return NotImplemented


def characteristic_coefficients(matrix):
    """EXERCISE 1b — return (trace, determinant) of a 2x2 matrix, as floats.

    These two numbers ARE the characteristic equation, because for any 2x2

        det(A - lambda*I) = lambda^2 - (trace)*lambda + (determinant)

    Steps:
      1. numpy.asarray(matrix, dtype=float).
      2. If its shape is not (2, 2), raise ValueError with a message
         containing "2x2".
      3. trace       = m[0,0] + m[1,1]
         determinant = m[0,0]*m[1,1] - m[0,1]*m[1,0]
      4. Return them as a tuple of two plain floats.
    """
    return NotImplemented


def eigenvalues_2x2(matrix):
    """EXERCISE 1c — solve the characteristic equation. Return two complex numbers.

        lambda = (trace +/- sqrt(trace^2 - 4*determinant)) / 2

    Steps:
      1. Get trace and determinant from characteristic_coefficients.
      2. discriminant = trace*trace - 4*determinant
      3. Take its square root with numpy.emath.sqrt, NOT numpy.sqrt.
         numpy.sqrt of a negative float returns nan and warns; numpy.emath.sqrt
         returns a complex number, which is the correct answer and is what
         makes the rotation case in exercise 3 work without special-casing.
      4. Return (complex((trace + root) / 2), complex((trace - root) / 2)).

    Always complex, even when the imaginary part is zero, so the caller has one
    code path instead of two.
    """
    return NotImplemented


def eigenvector_2x2(matrix, eigenvalue):
    """EXERCISE 1d — return a unit-length eigenvector for a real eigenvalue.

    Solve (A - lambda*I) v = 0 for a non-zero v.

    Why a solution exists: lambda is a root of the characteristic equation, so
    det(A - lambda*I) is 0, so that matrix squashes the plane onto a line
    (Day 102). A whole line of vectors gets sent to the origin, and any one of
    them is an eigenvector.

    Steps:
      1. shifted = matrix - eigenvalue * numpy.eye(2)
      2. For each ROW [p, q] of `shifted`, the candidate (-q, p) satisfies
         p*x + q*y = 0 automatically. Check it: p*(-q) + q*(p) = 0.
      3. Take the first row whose candidate is not the zero vector — compare
         its norm against a small tolerance, not against 0.0 exactly, because
         the entries are floats.
      4. Return that candidate divided by its own norm.
      5. If BOTH rows give the zero vector, the matrix was lambda*I and EVERY
         direction is an eigenvector; return numpy.array([1.0, 0.0]).

    Do not try to fix the sign. Whichever sign falls out is correct, and the
    test compares directions with abs_cosine for exactly that reason.
    """
    return NotImplemented


def power_method(matrix, start, tol=1e-10, max_iter=1000):
    """EXERCISE 1e — find the dominant eigenvector by repeated multiplication.

    Return a dict with keys "vector", "eigenvalue", "iterations", "converged".

    The algorithm:
      1. v = start as a float array, divided by its own norm.
         If start has norm 0.0, raise ValueError mentioning "zero vector".
      2. Loop `iteration` from 1 to max_iter:
           w = matrix @ v
           divide w by its norm
           IF numpy.dot(w, v) < 0, negate w
               — this aligns the signs. Without it, a matrix with a negative
                 dominant eigenvalue flips direction every single step and
                 your convergence test never fires, even though the ANSWER
                 converged on iteration three.
           change = numpy.linalg.norm(w - v)
           v = w
           if change < tol: return the dict with converged=True
      3. If the loop finishes, return the dict with converged=False and
         iterations=max_iter. Report the failure; do not raise and do not
         pretend it converged.

    For "eigenvalue" use the Rayleigh quotient: (v . matrix @ v) / (v . v).

    Normalising each round is not optional. Without it the vector's length
    multiplies by the eigenvalue every step and overflows to inf after a few
    hundred rounds, destroying a direction that was already correct.
    """
    return NotImplemented


def covariance_matrix(data):
    """EXERCISE 1f — the covariance matrix of an (n_points, n_features) array.

    Steps:
      1. numpy.asarray(data, dtype=float).
      2. If it is not 2-D, raise ValueError with a message containing "2-D".
      3. CENTRE IT: subtract data.mean(axis=0). Do not skip this. Without it
         you measure where the cloud sits rather than how it is shaped, and
         exercise 5 shows the answer coming out 136 degrees wrong with no
         error raised.
      4. Return (centred.T @ centred) / (n_points - 1).

    The result is always symmetric, and that symmetry is what guarantees PCA
    gets real eigenvalues and perpendicular eigenvectors.
    """
    return NotImplemented

"""Your running score. Run from the lab directory:

    .venv/bin/pytest starter -q

A SKIP means "not attempted yet". A FAILURE means "attempted and wrong", and
it prints your answer beside the real one. On an untouched checkout this is
1 passed, 52 skipped. When it says 53 passed, you are finished.

Every float comparison here names its tolerance. Every eigenvector comparison
goes through abs_cosine, because an eigenvector is defined only up to sign and
scale — so if your (1, -2) comes back as (-1, 2), or as (0.447, -0.894), all
three are correct and the tests treat them as equal.
"""

from __future__ import annotations

import numpy as np
import pytest

import answers
import eigen
from dataset import (
    A,
    CENTRE,
    ELONGATION_DEG,
    N_POINTS,
    PROJECTION_ONTO_X,
    REFLECTION_IN_X,
    ROTATION_90,
    SHEAR,
    SPREAD_ALONG,
    SYMMETRIC,
    elongation_direction,
    make_cloud,
    power_method_start,
)

EXACT = 1e-12
TIGHT = 1e-9

#: The known-correct facts, used to check your predictions. Reading this list
#: is of course possible, and equally of course misses the entire point.
TRUE_EIGENVALUES = (5.0, 2.0)
TRUE_EIGENVECTORS = ((1.0, 1.0), (1.0, -2.0))


def written(function):
    """Skip cleanly if the exercise has not been attempted yet."""
    try:
        result = function()
    except NotImplementedError:
        result = NotImplemented
    if result is NotImplemented:
        pytest.skip("not attempted yet — write this function in starter/eigen.py")
    return result


def answered(name):
    """Skip cleanly if this prediction is still None."""
    value = getattr(answers, name)
    if value is None:
        pytest.skip(f"not attempted yet — set {name} in starter/answers.py")
    return value


def test_the_environment_is_ready():
    """The one test that passes on an untouched checkout."""
    assert np.__version__.split(".")[0] >= "2"
    assert A.shape == (2, 2)
    assert make_cloud().shape == (N_POINTS, 2)


# ==========================================================================
# Exercise 1 — the six functions
# ==========================================================================


def test_1a_abs_cosine_on_the_same_line():
    fn = written(lambda: eigen.abs_cosine([1.0, 1.0], [1.0, 1.0]))
    assert fn == pytest.approx(1.0, abs=EXACT)


def test_1a_abs_cosine_ignores_a_sign_flip():
    value = written(lambda: eigen.abs_cosine([1.0, -2.0], [-1.0, 2.0]))
    assert value == pytest.approx(1.0, abs=EXACT), (
        "the sign must not matter: (1, -2) and (-1, 2) name the same line. "
        "Did you take the ABSOLUTE value of the dot product?"
    )


def test_1a_abs_cosine_ignores_scale():
    value = written(lambda: eigen.abs_cosine([1.0, 1.0], [37.5, 37.5]))
    assert value == pytest.approx(1.0, abs=EXACT)


def test_1a_abs_cosine_at_right_angles():
    value = written(lambda: eigen.abs_cosine([1.0, 0.0], [0.0, 1.0]))
    assert value == pytest.approx(0.0, abs=EXACT)


def test_1a_abs_cosine_refuses_the_zero_vector():
    written(lambda: eigen.abs_cosine([1.0, 1.0], [1.0, 1.0]))
    with pytest.raises(ValueError, match="no direction"):
        eigen.abs_cosine([0.0, 0.0], [1.0, 1.0])


def test_1b_characteristic_coefficients():
    value = written(lambda: eigen.characteristic_coefficients(A))
    assert tuple(value) == pytest.approx((7.0, 10.0), abs=EXACT)


def test_1b_characteristic_coefficients_rejects_a_3x3():
    written(lambda: eigen.characteristic_coefficients(A))
    with pytest.raises(ValueError, match="2x2"):
        eigen.characteristic_coefficients(np.eye(3))


def test_1c_eigenvalues_of_a():
    values = written(lambda: eigen.eigenvalues_2x2(A))
    assert sorted(complex(v).real for v in values) == pytest.approx([2.0, 5.0], abs=TIGHT)


def test_1c_eigenvalues_of_a_have_no_imaginary_part():
    values = written(lambda: eigen.eigenvalues_2x2(A))
    assert all(abs(complex(v).imag) < EXACT for v in values)


def test_1c_a_rotation_gives_complex_eigenvalues():
    values = written(lambda: eigen.eigenvalues_2x2(ROTATION_90))
    assert all(abs(complex(v).imag) > 0.5 for v in values), (
        "a 90-degree rotation has NO real eigenvalues. Did you use "
        "numpy.emath.sqrt rather than numpy.sqrt?"
    )


def test_1c_matches_numpy():
    values = written(lambda: eigen.eigenvalues_2x2(A))
    mine = sorted(complex(v).real for v in values)
    theirs = sorted(np.linalg.eig(A)[0].real)
    assert mine == pytest.approx(theirs, abs=TIGHT)


@pytest.mark.parametrize("eigenvalue,expected", list(zip(TRUE_EIGENVALUES, TRUE_EIGENVECTORS)))
def test_1d_eigenvector_lies_on_the_right_line(eigenvalue, expected):
    v = written(lambda: eigen.eigenvector_2x2(A, eigenvalue))
    cosine = abs(float(np.dot(v, expected))) / (
        float(np.linalg.norm(v)) * float(np.linalg.norm(expected))
    )
    assert cosine == pytest.approx(1.0, abs=TIGHT), (
        f"for lambda = {eigenvalue} you returned {np.asarray(v).tolist()}, which is not "
        f"on the same line as {list(expected)}"
    )


@pytest.mark.parametrize("eigenvalue", TRUE_EIGENVALUES)
def test_1d_eigenvector_satisfies_the_equation(eigenvalue):
    v = written(lambda: eigen.eigenvector_2x2(A, eigenvalue))
    assert np.allclose(A @ np.asarray(v), eigenvalue * np.asarray(v), atol=TIGHT)


def test_1d_eigenvector_is_unit_length():
    v = written(lambda: eigen.eigenvector_2x2(A, 5.0))
    assert float(np.linalg.norm(v)) == pytest.approx(1.0, abs=TIGHT)


def test_1e_power_method_converges():
    result = written(lambda: eigen.power_method(A, power_method_start(), tol=1e-10))
    assert result["converged"] is True
    assert result["iterations"] == 25, (
        f"expected 25 iterations to reach 1e-10, got {result['iterations']}. "
        "Check that you normalise BEFORE measuring the change, and that you "
        "align the signs."
    )


def test_1e_power_method_finds_the_dominant_line():
    result = written(lambda: eigen.power_method(A, power_method_start(), tol=1e-10))
    v = np.asarray(result["vector"], dtype=float)
    cosine = abs(float(np.dot(v, [1.0, 1.0]))) / (float(np.linalg.norm(v)) * np.sqrt(2.0))
    assert cosine == pytest.approx(1.0, abs=1e-10)


def test_1e_power_method_finds_the_dominant_eigenvalue():
    result = written(lambda: eigen.power_method(A, power_method_start(), tol=1e-10))
    assert float(result["eigenvalue"]) == pytest.approx(5.0, abs=TIGHT)


def test_1e_power_method_works_from_many_starts():
    written(lambda: eigen.power_method(A, power_method_start(), tol=1e-10))
    rng = np.random.default_rng(99)
    for _ in range(6):
        result = eigen.power_method(A, rng.normal(size=2), tol=1e-10)
        v = np.asarray(result["vector"], dtype=float)
        cosine = abs(float(np.dot(v, [1.0, 1.0]))) / (float(np.linalg.norm(v)) * np.sqrt(2.0))
        assert cosine == pytest.approx(1.0, abs=1e-10)


def test_1e_power_method_handles_a_negative_dominant_eigenvalue():
    """[[-5, 0], [0, 2]] has dominant eigenvalue -5, so the raw iteration
    flips sign every step. Without the sign alignment this never converges."""
    written(lambda: eigen.power_method(A, power_method_start(), tol=1e-10))
    result = eigen.power_method(np.diag([-5.0, 2.0]), np.array([0.6, 0.8]), tol=1e-10)
    assert result["converged"] is True, (
        "the dominant eigenvalue is negative, so w flips sign every step. "
        "Negate w when numpy.dot(w, v) < 0."
    )
    assert float(result["eigenvalue"]) == pytest.approx(-5.0, abs=TIGHT)


def test_1e_power_method_refuses_the_zero_vector():
    written(lambda: eigen.power_method(A, power_method_start(), tol=1e-10))
    with pytest.raises(ValueError, match="zero vector"):
        eigen.power_method(A, np.zeros(2))


def test_1e_power_method_reports_failure_rather_than_lying():
    written(lambda: eigen.power_method(A, power_method_start(), tol=1e-10))
    result = eigen.power_method(A, power_method_start(), tol=1e-16, max_iter=5)
    assert result["converged"] is False


def test_1f_covariance_matches_numpy():
    cloud = make_cloud()
    mine = written(lambda: eigen.covariance_matrix(cloud))
    assert np.allclose(mine, np.cov(cloud, rowvar=False), atol=EXACT), (
        "if this is close but not equal, check that you divided by (n - 1) "
        "and not by n; if it is far out, check that you subtracted the mean"
    )


def test_1f_covariance_is_symmetric():
    mine = written(lambda: eigen.covariance_matrix(make_cloud()))
    assert np.allclose(mine, np.asarray(mine).T, atol=1e-15)


def test_1f_covariance_rejects_a_1d_array():
    written(lambda: eigen.covariance_matrix(make_cloud()))
    with pytest.raises(ValueError, match="2-D"):
        eigen.covariance_matrix(np.array([1.0, 2.0, 3.0]))


# ==========================================================================
# Exercise 2 — the hand solution
# ==========================================================================


def test_2a_trace():
    assert answered("TRACE_OF_A") == 7


def test_2b_determinant():
    assert answered("DETERMINANT_OF_A") == 10


def test_2c_characteristic_coefficients():
    assert tuple(answered("CHARACTERISTIC_COEFFICIENTS")) == (7, 10)


def test_2d_discriminant():
    assert answered("DISCRIMINANT") == 9


def test_2e_eigenvalues():
    assert tuple(answered("EIGENVALUES_LARGEST_FIRST")) == (5, 2)


def test_2f_eigenvector_for_the_larger_eigenvalue():
    v = np.asarray(answered("EIGENVECTOR_FOR_LARGER"), dtype=float)
    assert np.allclose(A @ v, 5.0 * v, atol=TIGHT), (
        f"A @ {v.tolist()} is {(A @ v).tolist()}, which is not 5 times {v.tolist()}"
    )


def test_2g_eigenvector_for_the_smaller_eigenvalue():
    v = np.asarray(answered("EIGENVECTOR_FOR_SMALLER"), dtype=float)
    assert np.allclose(A @ v, 2.0 * v, atol=TIGHT), (
        f"A @ {v.tolist()} is {(A @ v).tolist()}, which is not 2 times {v.tolist()}"
    )


# ==========================================================================
# Exercise 3 — the standard transformations
# ==========================================================================


def test_3a_shear_has_one_eigen_line():
    assert answered("SHEAR_EIGEN_LINE_COUNT") == 1, (
        "numpy.linalg.eig returns two COLUMNS for the shear, but check whether "
        "they point along different lines"
    )
    values, vectors = np.linalg.eig(SHEAR)
    cosine = abs(float(np.dot(vectors.real[:, 0], vectors.real[:, 1])))
    assert cosine == pytest.approx(1.0, abs=1e-8)


def test_3b_rotation_eigenvalues_are_not_real():
    assert answered("ROTATION_EIGENVALUES_ARE_REAL") is False
    assert np.all(np.abs(np.linalg.eig(ROTATION_90)[0].imag) > 0.5)


def test_3c_rotation_eigenvalue_magnitude():
    assert float(answered("ROTATION_EIGENVALUE_MAGNITUDE")) == pytest.approx(1.0, abs=TIGHT)


def test_3d_eig_dtype_on_a_real_matrix():
    predicted = answered("EIG_DTYPE_ON_A")
    observed = str(np.linalg.eig(A)[0].dtype)
    assert predicted == observed, (
        f"you predicted {predicted!r}; numpy {np.__version__} on this machine "
        f"returned {observed!r}. Both eigenvalues are real and it still handed "
        "back complex. Its own docstring says otherwise. When documentation "
        "and measurement disagree, the measurement wins."
    )


def test_3e_projection_determinant_and_smallest_eigenvalue():
    pair = tuple(float(x) for x in answered("PROJECTION_DET_AND_SMALLEST_EIGENVALUE"))
    assert pair == pytest.approx((0.0, 0.0), abs=TIGHT)
    assert float(np.linalg.det(PROJECTION_ONTO_X)) == pytest.approx(0.0, abs=EXACT)
    assert float(np.min(np.linalg.eig(PROJECTION_ONTO_X)[0].real)) == pytest.approx(0.0, abs=TIGHT)


def test_3f_reflection_eigenvalues():
    assert tuple(answered("REFLECTION_EIGENVALUES")) == (1, -1)
    assert np.sort(np.linalg.eig(REFLECTION_IN_X)[0].real) == pytest.approx([-1.0, 1.0], abs=TIGHT)


def test_3g_symmetric_eigenvectors_are_perpendicular():
    assert float(answered("SYMMETRIC_EIGENVECTOR_ANGLE_DEG")) == pytest.approx(90.0, abs=TIGHT)
    vectors = np.linalg.eigh(SYMMETRIC)[1]
    assert float(np.dot(vectors[:, 0], vectors[:, 1])) == pytest.approx(0.0, abs=EXACT)


# ==========================================================================
# Exercise 4 — the power method
# ==========================================================================


def test_4a_which_eigenvalue():
    assert float(answered("POWER_METHOD_FINDS_EIGENVALUE")) == pytest.approx(5.0, abs=TIGHT)


def test_4b_which_direction():
    assert float(answered("POWER_METHOD_FINDS_DIRECTION_DEG")) == pytest.approx(45.0, abs=1e-6)


def test_4c_convergence_ratio():
    assert float(answered("CONVERGENCE_RATIO")) == pytest.approx(0.4, abs=1e-6), (
        "the surviving error is the second eigenvector's share, which shrinks "
        "relative to the first by |lambda2 / lambda1| every round"
    )


def test_4d_close_eigenvalues_are_slower():
    assert answered("CLOSE_EIGENVALUES_NEED") == "more"


def test_4e_unnormalised_length():
    assert answered("UNNORMALISED_LENGTH_BEHAVIOUR") == "overflows to inf"
    v = power_method_start()
    with np.errstate(over="ignore"):
        for _ in range(600):
            v = A @ v
    assert not np.all(np.isfinite(v))


# ==========================================================================
# Exercise 5 — PCA
# ==========================================================================


def test_5a_covariance_shape():
    assert tuple(answered("COVARIANCE_SHAPE")) == (2, 2)


def test_5b_covariance_is_symmetric():
    assert answered("COVARIANCE_IS_SYMMETRIC") is True


def test_5c_top_component_direction():
    predicted = float(answered("TOP_COMPONENT_DIRECTION_DEG"))
    assert predicted == pytest.approx(ELONGATION_DEG, abs=1.0)


def test_5d_sqrt_of_the_top_eigenvalue():
    predicted = float(answered("SQRT_OF_TOP_EIGENVALUE"))
    assert predicted == pytest.approx(SPREAD_ALONG, abs=0.5)


def test_5e_right_routine_for_a_covariance_matrix():
    assert answered("RIGHT_ROUTINE_FOR_COVARIANCE") == "eigh", (
        "a covariance matrix is symmetric by construction, and eigh is the "
        "routine written for symmetric input: real dtype, sorted values, and "
        "measurably faster"
    )


def test_5f_allclose_says_false_on_a_correct_answer():
    assert answered("ALLCLOSE_ON_THE_CORRECT_COMPONENT") is False, (
        "the answer is correct and numpy.allclose still says False, because "
        "eigh returned the reversed sign. An eigenvector names an AXIS, not "
        "an arrow. This is the trap the whole lab is built around."
    )
    cloud = make_cloud()
    centred = cloud - cloud.mean(axis=0)
    covariance = (centred.T @ centred) / (N_POINTS - 1)
    values, vectors = np.linalg.eigh(covariance)
    top = vectors[:, int(np.argmax(values))]
    truth = elongation_direction()
    assert not np.allclose(top, truth, atol=1e-3)
    assert abs(float(np.dot(top, truth))) > 0.999


def test_5g_forgetting_to_centre_is_not_survivable():
    assert answered("UNCENTRED_STILL_CORRECT") is False
    cloud = make_cloud()
    uncentred = (cloud.T @ cloud) / (N_POINTS - 1)
    values, vectors = np.linalg.eigh(uncentred)
    top = vectors[:, int(np.argmax(values))]
    angle = float(np.degrees(np.arctan2(top[1], top[0])) % 180.0)
    assert abs(angle - ELONGATION_DEG) > 5.0
    assert np.allclose(cloud.mean(axis=0), CENTRE, atol=0.2)

#!/usr/bin/env bash
# Tests for the Day 106 lab. Run from the lab directory:
#   bash tests/run_tests.sh
#
# The harness proves the lesson's claims by running code and reading real
# values, never by reading source:
#
#   * out of twenty-four directions spread around the circle, exactly two come
#     back on their own line -- and those two are the SAME line pointing
#     opposite ways, so one line survived, not two;
#   * a sweep of 180,000 directions finds the second eigen-line the coarse fan
#     stepped straight over, at 116.565 degrees;
#   * the hand solution -- trace 7, determinant 10, discriminant 9, eigenvalues
#     5 and 2 -- agrees with numpy.linalg.eig to 1e-9;
#   * numpy.linalg.eig returns complex128 for a real matrix with two real
#     eigenvalues, which contradicts its own docstring and is recorded as
#     measured rather than smoothed over;
#   * a shear has ONE eigen-line while eig returns TWO columns, and both
#     columns lie on that one line;
#   * a plane rotation has no real eigenvalues at all, and taking .real without
#     checking silently reports two eigenvalues of magnitude 1 as 0;
#   * the eigenvalues of every matrix here multiply to the determinant and add
#     to the trace;
#   * the power method converges in 25 iterations to 1e-10, its error shrinks
#     by the eigenvalue ratio 0.4 every round, and the same code needs 962
#     iterations when that ratio is 0.98;
#   * un-normalised iteration overflows to inf and then to nan, destroying a
#     direction that was already correct;
#   * PCA on a cloud built along 30 degrees recovers 30.101 degrees from the
#     coordinates alone -- and returns it with the sign flipped, which is
#     correct and is the trap the whole lab is built around;
#   * forgetting to centre gives a confident answer 136.58 degrees wrong;
#   * nothing is downloaded, and nothing is left behind on disk.
#
# Everything runs offline. Nothing binds a port, nothing writes outside the
# lab or a temporary directory, nothing needs a key. Deterministic,
# non-interactive, exits 0 only if every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1

lab_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"

# Bytecode left by an EARLIER command is not this run's litter. The README
# documents `pytest starter -q`, and running it writes .pyc files that would
# then fail the cleanliness check at the end of this script -- failing the
# reader for following the instructions. Clearing them here makes that final
# check measure what it claims to: what THIS run left behind. `.venv` is
# untouched, because the packages' own bytecode is theirs, not ours.
find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -exec rm -rf {} + 2>/dev/null || true
find "${lab_dir}" -name '.venv' -prune -o -type d -name '.pytest_cache' -exec rm -rf {} + 2>/dev/null || true

failures=0
checks=0

check() {
  local label="$1" ok="$2"
  checks=$((checks + 1))
  if [ "${ok}" = "yes" ]; then
    echo "  ok: ${label}"
  else
    echo "  FAIL: ${label}"
    failures=$((failures + 1))
  fi
}

check_eq() {
  # check_eq <label> <expected> <actual>
  if [ "$2" = "$3" ]; then
    check "$1" "yes"
  else
    check "$1 (expected [$2], got [$3])" "no"
  fi
}

# Resolve pytest: an explicit override, then this lab's .venv, then PATH.
# Fails loudly with instructions rather than silently skipping checks.
resolve_tool() {
  local tool="$1" override="$2"
  if [ -n "${override}" ] && [ -x "${override}" ]; then echo "${override}"; return 0; fi
  if [ -x "${lab_dir}/.venv/bin/${tool}" ]; then echo "${lab_dir}/.venv/bin/${tool}"; return 0; fi
  if command -v "${tool}" >/dev/null 2>&1; then command -v "${tool}"; return 0; fi
  return 1
}

pytest_bin="$(resolve_tool pytest "${PYTEST:-}")" || {
  echo "FAIL: pytest not found." >&2
  echo "  Install the lab's dependencies with:" >&2
  echo "    python3 -m venv .venv" >&2
  echo "    .venv/bin/pip install -r requirements/requirements.txt" >&2
  echo "  Or point this suite at an existing pytest:" >&2
  echo "    PYTEST=/path/to/pytest bash tests/run_tests.sh" >&2
  exit 1
}

# The Python that owns that pytest is the one with numpy installed.
python_bin="$(dirname "${pytest_bin}")/python3"
if [ ! -x "${python_bin}" ]; then
  python_bin="$(command -v python3 || true)"
fi
if [ -z "${python_bin}" ]; then
  echo "FAIL: python3 not found on PATH." >&2
  exit 1
fi

if ! "${python_bin}" -c "import numpy" >/dev/null 2>&1; then
  echo "FAIL: numpy is not importable from ${python_bin}." >&2
  echo "  Install the lab's dependencies with:" >&2
  echo "    python3 -m venv .venv" >&2
  echo "    .venv/bin/pip install -r requirements/requirements.txt" >&2
  exit 1
fi

echo "Day 106 — The Vectors That Keep Their Direction"
echo

# --------------------------------------------------------------------------
echo "1. The tools and the versions this lab was written against"
# --------------------------------------------------------------------------

versions="$("${python_bin}" - <<'PY'
import platform
import sys
from importlib.metadata import version

print(f"python   {platform.python_version()}")
for name in ("numpy", "pytest"):
    print(f"{name:<8} {version(name)}")
print(f"platform {platform.platform()}")
print(f"exe      {sys.executable.rsplit('/', 3)[-1]}")
PY
)"
echo "${versions}" | sed 's/^/  /'

for package in numpy pytest; do
  pinned="$(grep -iE "^${package}==" "${lab_dir}/requirements/requirements.txt" | cut -d= -f3)"
  installed="$("${python_bin}" -c "from importlib.metadata import version; print(version('${package}'))")"
  check_eq "installed ${package} matches requirements.txt" "${pinned}" "${installed}"
done

major="$("${python_bin}" -c "import numpy; print(numpy.__version__.split('.')[0])")"
check_eq "numpy is version 2 or later" "2" "${major}"

# --------------------------------------------------------------------------
echo
echo "2. Every reference script runs and every assertion inside it holds"
# --------------------------------------------------------------------------

for script in 01_the_fan_of_vectors 02_by_hand_2x2 03_standard_transformations \
              04_power_method 05_pca_from_covariance 06_eig_against_eigh; do
  out="$(cd "${lab_dir}/examples" && "${python_bin}" "${script}.py" 2>&1)"
  status=$?
  if [ "${status}" -ne 0 ]; then
    check "${script}.py exits 0" "no"
    echo "${out}" | tail -5 | sed 's/^/      /'
  else
    check "${script}.py exits 0" "yes"
  fi
  case "${out}" in
    *"${script}.py: every assertion held."*)
      check "${script}.py reports every assertion held" "yes" ;;
    *) check "${script}.py reports every assertion held" "no" ;;
  esac
done

# --------------------------------------------------------------------------
echo
echo "3. The reference pytest suite: real values, stated tolerances"
# --------------------------------------------------------------------------

ref_out="$(cd "${lab_dir}" && "${pytest_bin}" examples -q -p no:cacheprovider 2>&1)"
ref_status=$?
echo "${ref_out}" | tail -3 | sed 's/^/  /'
if [ "${ref_status}" -eq 0 ]; then
  check "pytest examples exits 0" "yes"
else
  check "pytest examples exits 0" "no"
fi
case "${ref_out}" in
  *" failed"*) check "no test in the reference suite failed" "no" ;;
  *)           check "no test in the reference suite failed" "yes" ;;
esac
ref_passed="$(printf '%s\n' "${ref_out}" | grep -o '[0-9][0-9]* passed' | head -1 | cut -d' ' -f1)"
if [ "${ref_passed:-0}" -ge 90 ]; then
  check "the reference suite ran at least 90 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 90 tests (ran ${ref_passed:-0})" "no"
fi

# --------------------------------------------------------------------------
echo
echo "4. The starter suite skips unattempted work instead of failing it"
# --------------------------------------------------------------------------

start_out="$(cd "${lab_dir}" && "${pytest_bin}" starter -q -p no:cacheprovider 2>&1)"
start_status=$?
echo "${start_out}" | tail -3 | sed 's/^/  /'
if [ "${start_status}" -eq 0 ]; then
  check "pytest starter exits 0 on an untouched checkout" "yes"
else
  check "pytest starter exits 0 on an untouched checkout" "no"
fi
case "${start_out}" in
  *" failed"*) check "the starter suite reports no failures" "no" ;;
  *)           check "the starter suite reports no failures" "yes" ;;
esac
case "${start_out}" in
  *skipped*) check "unwritten exercises are reported as skipped, not passed" "yes" ;;
  *) check "unwritten exercises are reported as skipped, not passed" "no" ;;
esac

# The import guard. Both directories contain modules called `eigen` and
# `dataset`, and pytest imports test files by putting their directory on
# sys.path -- so collecting both suites at once would otherwise let the starter
# tests import the REFERENCE solution and report unwritten exercises as
# passing. Each directory's conftest.py prevents that. This check proves it
# still does: across both suites, the skip count must be unchanged.
both_out="$(cd "${lab_dir}" && "${pytest_bin}" -q -p no:cacheprovider 2>&1)"
start_skipped="$(printf '%s\n' "${start_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
both_skipped="$(printf '%s\n' "${both_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
check_eq "collecting both suites at once does not turn skips into passes" \
  "${start_skipped:-none}" "${both_skipped:-none}"

# --------------------------------------------------------------------------
echo
echo "5. The lesson's claims, checked one value at a time"
# --------------------------------------------------------------------------

facts="$(cd "${lab_dir}/examples" && "${python_bin}" - <<'PY'
import numpy as np

import dataset
import eigen
from dataset import (
    A, A_EIGENVALUES, A_EIGENVECTORS, ELONGATION_DEG, PROJECTION_ONTO_X,
    REFLECTION_IN_X, ROTATION_60, ROTATION_90, SHEAR, SPREAD_ALONG,
    STANDARD_TRANSFORMATIONS, SYMMETRIC, SYMMETRIC_3X3, elongation_direction,
    make_cloud, power_method_start,
)

# -- the fan of vectors
kept = [
    a for a in range(0, 360, 15)
    if eigen.deviation_degrees(
        [np.cos(np.radians(a)), np.sin(np.radians(a))],
        A @ np.array([np.cos(np.radians(a)), np.sin(np.radians(a))]),
    ) < 1e-9
]
print("fan_kept", kept)
up = np.array([np.cos(np.radians(45)), np.sin(np.radians(45))])
down = np.array([np.cos(np.radians(225)), np.sin(np.radians(225))])
print("fan_45_and_225_same_line", bool(np.allclose(down, -up, atol=1e-12)))
print("fan_x_axis_swing", round(eigen.deviation_degrees([1.0, 0.0], A @ np.array([1.0, 0.0])), 6))

swept = eigen.eigen_lines_by_sweep(A)
print("sweep_verdict", swept["verdict"])
print("sweep_line_count", len(swept["lines"]))
print("sweep_lines_match_exact_angles",
      bool(np.allclose(swept["lines"], dataset.A_EIGEN_ANGLES_DEG, atol=1e-2)))

for angle, expected in zip(dataset.A_EIGEN_ANGLES_DEG, A_EIGENVALUES):
    v = np.array([np.cos(np.radians(angle)), np.sin(np.radians(angle))])
    print(f"stretch_at_{expected:.0f}", round(float(np.linalg.norm(A @ v)), 12))

print("integer_checks_hold", bool(all(
    np.allclose(A @ np.array(v), lam * np.array(v), atol=1e-12)
    for v, lam in zip(A_EIGENVECTORS, A_EIGENVALUES))))
print("zero_vector_fits_every_lambda", bool(all(
    np.allclose(A @ np.zeros(2), c * np.zeros(2), atol=1e-12)
    for c in (5.0, 2.0, 0.0, -3.5, 1000.0))))

# -- the hand solution
trace, det = eigen.characteristic_coefficients(A)
print("trace", trace)
print("determinant", det)
print("discriminant", trace * trace - 4.0 * det)
hand = sorted(v.real for v in eigen.eigenvalues_2x2(A))
print("hand_eigenvalues", [round(v, 12) for v in hand])
print("hand_matches_numpy",
      bool(np.allclose(hand, sorted(np.linalg.eig(A)[0].real), atol=1e-9)))
print("det_of_a_minus_lambda_i_is_zero", bool(all(
    abs(float(np.linalg.det(A - lam * np.eye(2)))) < 1e-12 for lam in A_EIGENVALUES)))

# -- what numpy actually returns
values, vectors = np.linalg.eig(A)
print("eig_dtype", str(values.dtype))
print("eig_imag_all_zero", bool(np.all(values.imag == 0.0)))
print("eig_dtype_on_identity", str(np.linalg.eig(np.eye(2))[0].dtype))
smaller = int(np.argmin(values.real))
mine = np.array(A_EIGENVECTORS[1], dtype=float)
mine = mine / np.linalg.norm(mine)
theirs = vectors.real[:, smaller]
print("sign_flip_allclose_says_false", bool(not np.allclose(mine, theirs, atol=1e-6)))
print("sign_flip_is_exact_negative", bool(np.allclose(mine, -theirs, atol=1e-9)))
print("sign_flip_abs_cosine", round(eigen.abs_cosine(mine, theirs), 12))

worst = 0.0
for m in (A, SYMMETRIC, SHEAR, REFLECTION_IN_X, PROJECTION_ONTO_X):
    vals, vecs = np.linalg.eig(m)
    for i in range(len(vals)):
        worst = max(worst, float(np.abs(m @ vecs[:, i] - vals[i] * vecs[:, i]).max()))
print("worst_a_v_minus_lambda_v_residual_below_1e12", bool(worst < 1e-12))

# -- the standard transformations
shear_found = eigen.eigen_lines_by_sweep(SHEAR)
sv, svec = np.linalg.eig(SHEAR)
print("shear_eigen_lines", len(shear_found["lines"]))
print("shear_eig_columns", svec.shape[1])
print("shear_columns_same_line",
      round(eigen.abs_cosine(svec.real[:, 0], svec.real[:, 1]), 8))
print("shear_eigenvalues", [round(float(v), 12) for v in sv.real])
print("shear_eigenvector_matrix_is_singular",
      bool(abs(float(np.linalg.det(svec.real))) < 1e-15))

for name, rot in (("90", ROTATION_90), ("60", ROTATION_60)):
    rv = np.linalg.eig(rot)[0]
    print(f"rotation{name}_all_complex", bool(np.all(np.abs(rv.imag) > 0.5)))
    print(f"rotation{name}_magnitudes_are_one", bool(np.allclose(np.abs(rv), 1.0, atol=1e-9)))
    print(f"rotation{name}_verdict", eigen.eigen_lines_by_sweep(rot)["verdict"])
print("rotation90_real_parts_destroy_the_answer",
      bool(np.allclose(np.linalg.eig(ROTATION_90)[0].real, 0.0, atol=1e-9)))
rt, rd = eigen.characteristic_coefficients(ROTATION_90)
print("rotation90_discriminant", rt * rt - 4.0 * rd)

print("identity_verdict", eigen.eigen_lines_by_sweep(np.eye(2))["verdict"])
print("uniform_scale_verdict", eigen.eigen_lines_by_sweep(2.0 * np.eye(2))["verdict"])
print("reflection_eigenvalues",
      [round(float(v), 12) for v in np.sort(np.linalg.eig(REFLECTION_IN_X)[0].real)])
print("projection_eigenvalues",
      [round(float(v), 12) for v in np.sort(np.linalg.eig(PROJECTION_ONTO_X)[0].real)])
print("projection_determinant", round(float(np.linalg.det(PROJECTION_ONTO_X)), 12))
_devs, collapsed = eigen.sweep_deviations(PROJECTION_ONTO_X, [0.0, 45.0, 90.0])
print("projection_collapses_the_y_axis", collapsed.tolist())

prod_ok = trace_ok = True
for name in STANDARD_TRANSFORMATIONS:
    m = STANDARD_TRANSFORMATIONS[name][0]
    vals = np.linalg.eig(m)[0]
    prod_ok &= abs(complex(np.prod(vals)).real - float(np.linalg.det(m))) < 1e-9
    trace_ok &= abs(complex(np.sum(vals)).real - float(np.trace(m))) < 1e-9
print("eigenvalues_multiply_to_determinant_on_all_eight", bool(prod_ok))
print("eigenvalues_add_to_trace_on_all_eight", bool(trace_ok))

for label, m in (("2x2", SYMMETRIC), ("3x3", SYMMETRIC_3X3)):
    vals, vecs = np.linalg.eigh(m)
    print(f"symmetric_{label}_dtype", str(vals.dtype))
    print(f"symmetric_{label}_orthogonal",
          bool(float(np.abs(vecs.T @ vecs - np.eye(len(vals))).max()) < 1e-12))
print("eigh_sorted_ascending", bool(np.all(np.diff(np.linalg.eigh(SYMMETRIC)[0]) >= 0)))
print("eigh_on_non_symmetric_is_silently_wrong", bool(
    np.allclose(np.linalg.eigh(A)[0], np.linalg.eigvalsh([[4.0, 2.0], [2.0, 3.0]]), atol=1e-12)
    and not np.allclose(np.sort(np.linalg.eigh(A)[0]), sorted(A_EIGENVALUES), atol=1e-6)))

vals, vecs = np.linalg.eig(A)
rebuilt = vecs @ np.diag(vals) @ np.linalg.inv(vecs)
print("diagonalisation_rebuilds_a", bool(np.allclose(rebuilt.real, A, atol=1e-9)))

# -- the power method
result = eigen.power_method(A, power_method_start(), tol=1e-10)
print("power_iterations", result["iterations"])
print("power_converged", result["converged"])
print("power_change_below_tolerance", bool(result["change"] < 1e-10))
print("power_direction", round(eigen.direction_degrees(result["vector"]), 6))
print("power_eigenvalue", round(result["eigenvalue"], 9))
print("power_abs_cosine_with_1_1", round(eigen.abs_cosine(result["vector"], (1.0, 1.0)), 12))
top = int(np.argmax(np.abs(np.linalg.eig(A)[0].real)))
print("power_agrees_with_numpy_eigenvalue", bool(
    abs(result["eigenvalue"] - float(np.linalg.eig(A)[0].real[top])) < 1e-9))
hist = result["history"]
print("power_ratio_at_step_14", round(hist[13] / hist[12], 6))
slow = eigen.power_method(np.diag([5.0, 4.9]), np.array([0.6, 0.8]), tol=1e-10)
print("power_close_eigenvalues_iterations", slow["iterations"])
print("power_negative_dominant_converges", eigen.power_method(
    np.diag([-5.0, 2.0]), np.array([0.6, 0.8]), tol=1e-10)["converged"])
print("power_rayleigh_exact_on_true_eigenvector", bool(all(
    abs(eigen.rayleigh_quotient(A, v) - lam) < 1e-12
    for v, lam in zip(A_EIGENVECTORS, A_EIGENVALUES))))
print("power_reports_failure_rather_than_lying", eigen.power_method(
    A, power_method_start(), tol=1e-16, max_iter=5)["converged"])

v = power_method_start()
with np.errstate(over="ignore"):
    for _ in range(600):
        v = A @ v
print("unnormalised_overflows", bool(not np.all(np.isfinite(v))))
with np.errstate(invalid="ignore"):
    print("unnormalised_then_normalised_is_nan",
          bool(np.all(np.isnan(v / np.linalg.norm(v)))))

# -- PCA
cloud = make_cloud()
print("cloud_shape", cloud.shape)
print("cloud_is_reproducible", bool(np.array_equal(make_cloud(), make_cloud())))
cov = eigen.covariance_matrix(cloud)
print("covariance_matches_numpy_cov",
      bool(np.allclose(cov, np.cov(cloud, rowvar=False), atol=1e-12)))
print("covariance_is_symmetric", bool(np.allclose(cov, cov.T, atol=1e-15)))
print("covariance_shape", cov.shape)
variances, directions = eigen.principal_components(cloud)
truth = elongation_direction()
print("pca_top_direction", round(eigen.direction_degrees(directions[:, 0]), 6))
print("pca_true_direction", ELONGATION_DEG)
print("pca_abs_cosine", round(eigen.abs_cosine(directions[:, 0], truth), 10))
print("pca_top_component_points_the_other_way",
      bool(float(np.dot(directions[:, 0], truth)) < 0.0))
print("pca_allclose_says_false_on_a_correct_answer",
      bool(not np.allclose(directions[:, 0], truth, atol=1e-3)))
print("pca_sqrt_top_eigenvalue", round(float(np.sqrt(variances[0])), 6))
print("pca_sqrt_second_eigenvalue", round(float(np.sqrt(variances[1])), 6))
print("pca_variance_explained", round(float(variances[0] / variances.sum()), 6))
print("pca_components_perpendicular",
      bool(abs(float(np.dot(directions[:, 0], directions[:, 1]))) < 1e-12))
projected = (cloud - cloud.mean(axis=0)) @ directions
print("pca_projections_uncorrelated",
      bool(abs(float(np.corrcoef(projected.T)[0, 1])) < 1e-12))

uncentred = (cloud.T @ cloud) / (cloud.shape[0] - 1)
uv, uvec = np.linalg.eigh(uncentred)
uncentred_deg = eigen.direction_degrees(uvec[:, int(np.argmax(uv))])
print("pca_uncentred_direction", round(uncentred_deg, 6))
print("pca_uncentred_error", round(abs(uncentred_deg - ELONGATION_DEG), 6))
print("pca_eig_and_eigh_agree_on_covariance", bool(np.allclose(
    np.sort(np.linalg.eig(cov)[0].real), np.sort(np.linalg.eigh(cov)[0]), atol=1e-12)))
PY
)"

get() { printf '%s\n' "${facts}" | grep "^$1 " | cut -d' ' -f2-; }

# -- the fan
check_eq "of 24 directions, exactly 45 and 225 degrees keep their line" \
  "[45, 225]" "$(get fan_kept)"
check_eq "and 45 and 225 are the SAME line, so one line survived, not two" \
  "True" "$(get fan_45_and_225_same_line)"
check_eq "the x-axis is knocked 26.565051 degrees off its line" \
  "26.565051" "$(get fan_x_axis_swing)"
check_eq "a 180,000-direction sweep finds surviving lines" \
  "some" "$(get sweep_verdict)"
check_eq "and finds TWO of them: the coarse fan stepped over the second" \
  "2" "$(get sweep_line_count)"
check_eq "both agree with the exact angles 45 and 116.565051 to 0.01 degrees" \
  "True" "$(get sweep_lines_match_exact_angles)"
check_eq "the 45-degree direction is stretched by exactly 5" \
  "5.0" "$(get stretch_at_5)"
check_eq "the 116.565-degree direction is stretched by exactly 2" \
  "2.0" "$(get stretch_at_2)"
check_eq "the integer eigenvectors (1,1) and (1,-2) check out on paper" \
  "True" "$(get integer_checks_hold)"
check_eq "the zero vector satisfies A v = lambda v for EVERY lambda, so it is excluded" \
  "True" "$(get zero_vector_fits_every_lambda)"

# -- the hand solution
check_eq "trace of A is 7" "7.0" "$(get trace)"
check_eq "determinant of A is 10" "10.0" "$(get determinant)"
check_eq "discriminant is 9, so two distinct real eigenvalues" "9.0" "$(get discriminant)"
check_eq "the hand quadratic gives exactly 2 and 5" \
  "[2.0, 5.0]" "$(get hand_eigenvalues)"
check_eq "and matches numpy.linalg.eig, sorted, to 1e-9" \
  "True" "$(get hand_matches_numpy)"
check_eq "det(A - lambda I) is zero at each eigenvalue, which is the whole derivation" \
  "True" "$(get det_of_a_minus_lambda_i_is_zero)"

# -- what numpy returns
check_eq "numpy.linalg.eig returns complex128 for a real matrix" \
  "complex128" "$(get eig_dtype)"
check_eq "with every imaginary part exactly zero" "True" "$(get eig_imag_all_zero)"
check_eq "and does the same for numpy.eye(2), whose eigenvalues are both 1" \
  "complex128" "$(get eig_dtype_on_identity)"
check_eq "numpy.allclose says False on a CORRECT eigenvector" \
  "True" "$(get sign_flip_allclose_says_false)"
check_eq "because the two answers are exact negatives of each other" \
  "True" "$(get sign_flip_is_exact_negative)"
check_eq "and the absolute cosine, which asks the right question, says 1" \
  "1.0" "$(get sign_flip_abs_cosine)"
check_eq "every pair numpy returned satisfies A v = lambda v below 1e-12" \
  "True" "$(get worst_a_v_minus_lambda_v_residual_below_1e12)"

# -- the standard transformations
#
# Section 6 re-runs this script with D106_SELF_TEST=1, which swaps ONE
# expectation below for the naive belief that counting eig's columns counts
# eigendirections. That is how the harness proves it can fail rather than
# merely asserting that it could.
expected_shear_lines="1"
if [ -n "${D106_SELF_TEST:-}" ]; then
  expected_shear_lines="2"   # the naive belief, deliberately wrong here
fi
check_eq "a shear has exactly ONE eigen-line" \
  "${expected_shear_lines}" "$(get shear_eigen_lines)"
check_eq "while numpy.linalg.eig returns TWO columns for it" \
  "2" "$(get shear_eig_columns)"
check_eq "and those two columns lie on the same line" \
  "1.0" "$(get shear_columns_same_line)"
check_eq "its eigenvalue is 1, repeated" "[1.0, 1.0]" "$(get shear_eigenvalues)"
check_eq "so its eigenvector matrix is singular and it cannot be diagonalised" \
  "True" "$(get shear_eigenvector_matrix_is_singular)"

check_eq "a 90-degree rotation has no real eigenvalues" \
  "True" "$(get rotation90_all_complex)"
check_eq "both of magnitude 1, because a rotation changes no lengths" \
  "True" "$(get rotation90_magnitudes_are_one)"
check_eq "and a 180,000-direction sweep finds nothing that kept its line" \
  "none" "$(get rotation90_verdict)"
check_eq "the same holds for a 60-degree rotation" \
  "True" "$(get rotation60_all_complex)"
check_eq "with the same verdict from measurement" "none" "$(get rotation60_verdict)"
check_eq "taking .real without checking reports both eigenvalues as 0" \
  "True" "$(get rotation90_real_parts_destroy_the_answer)"
check_eq "the negative discriminant is the algebra reporting the geometry" \
  "-4.0" "$(get rotation90_discriminant)"

check_eq "the identity keeps EVERY direction" \
  "every direction" "$(get identity_verdict)"
check_eq "and so does a uniform scaling" \
  "every direction" "$(get uniform_scale_verdict)"
check_eq "a reflection has eigenvalues 1 and -1: one direction reversed" \
  "[-1.0, 1.0]" "$(get reflection_eigenvalues)"
check_eq "a projection has eigenvalue 0" "[0.0, 1.0]" "$(get projection_eigenvalues)"
check_eq "and determinant 0, which is Day 102's news arriving twice" \
  "0.0" "$(get projection_determinant)"
check_eq "the collapsed y-axis has no direction left to measure" \
  "[False, False, True]" "$(get projection_collapses_the_y_axis)"
check_eq "on all eight transformations the eigenvalues multiply to the determinant" \
  "True" "$(get eigenvalues_multiply_to_determinant_on_all_eight)"
check_eq "and add to the trace" "True" "$(get eigenvalues_add_to_trace_on_all_eight)"

check_eq "a symmetric 2x2 gives eigh real float64 values" \
  "float64" "$(get symmetric_2x2_dtype)"
check_eq "with eigenvectors at right angles" "True" "$(get symmetric_2x2_orthogonal)"
check_eq "and the same holds for a symmetric 3x3, so it is not a 2x2 accident" \
  "True" "$(get symmetric_3x3_orthogonal)"
check_eq "eigh returns its values sorted ascending; eig promises no order" \
  "True" "$(get eigh_sorted_ascending)"
check_eq "eigh on NON-symmetric input answers a different question, silently" \
  "True" "$(get eigh_on_non_symmetric_is_silently_wrong)"
check_eq "diagonalisation V D V-inverse rebuilds A exactly" \
  "True" "$(get diagonalisation_rebuilds_a)"

# -- the power method
check_eq "the power method converges in 25 iterations to 1e-10" \
  "25" "$(get power_iterations)"
check_eq "and says so rather than being assumed" "True" "$(get power_converged)"
check_eq "with the final change below the stated tolerance" \
  "True" "$(get power_change_below_tolerance)"
check_eq "it lands on the 45-degree eigen-line" "45.0" "$(get power_direction)"
check_eq "with abs_cosine 1 against (1, 1)" "1.0" "$(get power_abs_cosine_with_1_1)"
check_eq "and a Rayleigh quotient of 5 to nine decimal places" \
  "5.0" "$(get power_eigenvalue)"
check_eq "agreeing with numpy.linalg.eig to 1e-9" \
  "True" "$(get power_agrees_with_numpy_eigenvalue)"
check_eq "its error shrinks by the eigenvalue ratio 2/5 every round" \
  "0.399999" "$(get power_ratio_at_step_14)"
check_eq "eigenvalues of 5 and 4.9 need 962 iterations instead of 25" \
  "962" "$(get power_close_eigenvalues_iterations)"
check_eq "a NEGATIVE dominant eigenvalue still converges, given sign alignment" \
  "True" "$(get power_negative_dominant_converges)"
check_eq "the Rayleigh quotient is exact on a true eigenvector" \
  "True" "$(get power_rayleigh_exact_on_true_eigenvector)"
check_eq "and non-convergence is REPORTED, not hidden" \
  "False" "$(get power_reports_failure_rather_than_lying)"
check_eq "600 un-normalised rounds overflow to infinity" \
  "True" "$(get unnormalised_overflows)"
check_eq "and normalising afterwards gives nan: the direction is unrecoverable" \
  "True" "$(get unnormalised_then_normalised_is_nan)"

# -- PCA
check_eq "the cloud is 400 points in 2 dimensions" "(400, 2)" "$(get cloud_shape)"
check_eq "generated from a seed, so it is identical on every run" \
  "True" "$(get cloud_is_reproducible)"
check_eq "the from-scratch covariance matches numpy.cov to 1e-12" \
  "True" "$(get covariance_matches_numpy_cov)"
check_eq "a covariance matrix is always symmetric" \
  "True" "$(get covariance_is_symmetric)"
check_eq "and is 2 by 2 for a (400, 2) dataset, not 400 by 400" \
  "(2, 2)" "$(get covariance_shape)"
check_eq "PCA recovers 30.101134 degrees from the coordinates alone" \
  "30.101134" "$(get pca_top_direction)"
check_eq "against a true elongation of 30.0 degrees it was never told" \
  "30.0" "$(get pca_true_direction)"
check_eq "abs_cosine with the truth is 0.9999984422" \
  "0.9999984422" "$(get pca_abs_cosine)"
check_eq "and the component came back pointing the OTHER way along that line" \
  "True" "$(get pca_top_component_points_the_other_way)"
check_eq "so numpy.allclose says False on an answer that is exactly right" \
  "True" "$(get pca_allclose_says_false_on_a_correct_answer)"
check_eq "the top eigenvalue's square root recovers the spread 3.0 it was built with" \
  "2.902251" "$(get pca_sqrt_top_eigenvalue)"
check_eq "and the second recovers the across-spread 0.4" \
  "0.420504" "$(get pca_sqrt_second_eigenvalue)"
check_eq "the first component alone carries 97.9439% of the variance" \
  "0.979439" "$(get pca_variance_explained)"
check_eq "the two components are perpendicular" \
  "True" "$(get pca_components_perpendicular)"
check_eq "and the projections onto them are uncorrelated" \
  "True" "$(get pca_projections_uncorrelated)"
check_eq "forgetting to centre points the answer at 166.583965 degrees" \
  "166.583965" "$(get pca_uncentred_direction)"
check_eq "which is 136.583965 degrees wrong, with no error and no warning" \
  "136.583965" "$(get pca_uncentred_error)"
check_eq "eig and eigh agree on the covariance matrix" \
  "True" "$(get pca_eig_and_eigh_agree_on_covariance)"

# --------------------------------------------------------------------------
echo
echo "6. The harness can actually fail"
# --------------------------------------------------------------------------

# A green test suite proves nothing until you have watched it go red. This
# section re-runs the whole script with one expectation deliberately swapped
# for the naive belief that a shear has two eigendirections because eig returns
# two columns, and asserts that the re-run reports the failure and exits
# non-zero. If this section passes, section 5 is not decorative.
if [ -z "${D106_SELF_TEST:-}" ]; then
  self_out="$(D106_SELF_TEST=1 bash "${BASH_SOURCE[0]}" 2>&1)"
  self_status=$?
  if [ "${self_status}" -ne 0 ]; then
    check "a deliberately wrong expectation makes the harness exit non-zero (${self_status})" "yes"
  else
    check "a deliberately wrong expectation makes the harness exit non-zero" "no"
  fi
  case "${self_out}" in
    *"FAIL: a shear has exactly ONE eigen-line"*)
      check "the failing check is named in the output with both values" "yes" ;;
    *) check "the failing check is named in the output with both values" "no" ;;
  esac
  case "${self_out}" in
    *", 1 failure(s)."*)
      check "the summary line counts exactly one failure" "yes" ;;
    *) check "the summary line counts exactly one failure" "no" ;;
  esac
else
  echo "  (self-test run: section 6 does not recurse)"
fi

# --------------------------------------------------------------------------
echo
echo "7. Nothing was downloaded, and nothing was left behind"
# --------------------------------------------------------------------------

# Every find below PRUNES .venv first, and it is not optional. The README tells
# you to create a lab-local virtual environment, so `.venv` is the documented
# setup rather than litter -- and NumPy ships its own compiled bytecode inside
# it. Without the prune, this section would fail the lab for following its own
# installation instructions.
if find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -print -quit 2>/dev/null | grep -q .; then
  check "no __pycache__ directory left under the lab (ignoring .venv)" "no"
else
  check "no __pycache__ directory left under the lab (ignoring .venv)" "yes"
fi

if find "${lab_dir}" -name '.venv' -prune -o -type d -name '.pytest_cache' -print -quit 2>/dev/null | grep -q .; then
  check "no .pytest_cache directory left under the lab (ignoring .venv)" "no"
else
  check "no .pytest_cache directory left under the lab (ignoring .venv)" "yes"
fi

# The dataset is GENERATED from a seed, not downloaded and not committed. If a
# data file ever appears in the lab's own tree, either something was committed
# by mistake or a script wrote one and failed to clean up. NumPy ships plenty
# of its own data inside site-packages, so .venv is pruned here too.
data_files="$(find "${lab_dir}" -name '.venv' -prune -o -type f \
  \( -name '*.csv' -o -name '*.npy' -o -name '*.npz' -o -name '*.json' \
     -o -name '*.parquet' -o -name '*.pkl' \) -print 2>/dev/null \
  | wc -l | tr -d ' ')"
check_eq "no data file in the lab's own tree: the cloud is generated from a seed" \
  "0" "${data_files}"

if grep -rqE 'urlopen|requests\.|socket\.|http://|https://' \
     "${lab_dir}/examples" "${lab_dir}/starter" 2>/dev/null; then
  check "no lab source opens a network connection" "no"
else
  check "no lab source opens a network connection" "yes"
fi

echo
echo "${checks} checks, ${failures} failure(s)."
[ "${failures}" -eq 0 ]