Hands-on lab — Day 102: Linear Transformations

"""Where do the basis vectors land? Read a matrix off a picture, and back.

Run from inside examples/:

    ../.venv/bin/python3 01_columns_are_landings.py

The claim under test: if you know where (1, 0) and (0, 1) land, you know where
every vector lands, and you did not need the picture for anything else.
"""

import numpy as np

import shapes
from transforms import apply, columns_of, from_landings

TOL = shapes.TOL


def main() -> None:
    print("01_columns_are_landings.py")
    print("=" * 70)

    # -- 1. The picture, in words ---------------------------------------------
    print()
    print("1. The picture says two things, and only two things")
    print("-" * 70)
    print(f"  the arrow (1, 0) has been redrawn ending at {shapes.PICTURE_E1_LANDS_AT}")
    print(f"  the arrow (0, 1) has been redrawn ending at {shapes.PICTURE_E2_LANDS_AT}")
    print()
    print("  Write those two landing places down as COLUMNS. That is the matrix.")

    M = from_landings(shapes.PICTURE_E1_LANDS_AT, shapes.PICTURE_E2_LANDS_AT)
    print(f"  row 0: {M[0]}")
    print(f"  row 1: {M[1]}")
    assert M == shapes.PICTURE_MATRIX

    # -- 2. And back again ----------------------------------------------------
    print()
    print("2. Reading the picture back off the matrix")
    print("-" * 70)
    e1_lands, e2_lands = columns_of(M)
    print(f"  column 0 = {e1_lands}   <- where (1, 0) went")
    print(f"  column 1 = {e2_lands}   <- where (0, 1) went")
    print()
    print("  Careful: the matrix is WRITTEN as rows. Row 0 is", M[0], "and that")
    print("  is not a landing place. (3, 1) is, and it is read downwards.")
    assert e1_lands == shapes.PICTURE_E1_LANDS_AT
    assert e2_lands == shapes.PICTURE_E2_LANDS_AT

    # -- 3. Everything else follows -------------------------------------------
    print()
    print("3. Now every other vector, without looking at the picture again")
    print("-" * 70)
    v = (2.0, 1.0)
    print(f"  {v} = 2 * (1, 0) + 1 * (0, 1)")
    print("  A linear transformation keeps that combination intact, so it must")
    print("  land at 2 * (3, 1) + 1 * (-1, 2):")
    print("      2 * (3, 1) = (6, 2)")
    print("      1 * (-1, 2) = (-1, 2)")
    print("      (6, 2) + (-1, 2) = (5, 4)")
    landed = apply(M, v)
    print(f"  from-scratch apply:  {landed}")
    assert landed == shapes.PICTURE_SENDS_2_1_TO

    # -- 4. NumPy agrees ------------------------------------------------------
    print()
    print("4. NumPy, doing the same thing with the @ operator")
    print("-" * 70)
    npM = np.array(shapes.PICTURE_MATRIX)
    npv = np.array(v)
    print(f"  M @ v = {(npM @ npv).tolist()}")
    print(f"  M @ e1 = {(npM @ np.array(shapes.E1)).tolist()}   (column 0)")
    print(f"  M @ e2 = {(npM @ np.array(shapes.E2)).tolist()}   (column 1)")
    assert np.allclose(npM @ npv, np.array(landed), atol=TOL)
    assert np.allclose(npM @ np.array(shapes.E1), np.array(e1_lands), atol=TOL)
    assert np.allclose(npM @ np.array(shapes.E2), np.array(e2_lands), atol=TOL)

    # -- 5. The identity, which is the do-nothing case ------------------------
    print()
    print("5. The matrix that leaves the basis vectors exactly where they are")
    print("-" * 70)
    I = from_landings(shapes.E1, shapes.E2)
    print(f"  from_landings((1, 0), (0, 1)) = {I}")
    print("  That is the identity matrix, and it is the identity for one reason:")
    print("  nothing moved, so nothing moves.")
    for point in [(0.0, 0.0), (2.0, 1.0), (-3.5, 7.25)]:
        assert apply(I, point) == point
    print("  checked on (0, 0), (2, 1) and (-3.5, 7.25): each came back unchanged")
    assert np.allclose(np.array(I), np.eye(2), atol=TOL)

    print()
    print("01_columns_are_landings.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Scaling, reflection, shear and rotation -- each matrix derived, not given.

Run from inside examples/:

    ../.venv/bin/python3 02_building_the_transformations.py

Nothing here is memorised. Every matrix is built by asking one question --
where does (1, 0) go, and where does (0, 1) go -- and writing the two answers
down as columns.
"""

import math

import numpy as np

import shapes
from transforms import (
    apply,
    reflection_in_x_axis,
    rotation,
    scaling,
    shear_x,
    transform_polygon,
)

TOL = shapes.TOL


def show(name: str, matrix, derivation: str) -> None:
    e1 = (matrix[0][0], matrix[1][0])
    e2 = (matrix[0][1], matrix[1][1])
    print(f"  {name}")
    print(f"    {derivation}")
    print(f"    (1, 0) lands at {e1}      (0, 1) lands at {e2}")
    print(f"    matrix: [{matrix[0]}, {matrix[1]}]")


def main() -> None:
    print("02_building_the_transformations.py")
    print("=" * 70)

    # -- Scaling ---------------------------------------------------------------
    print()
    print("1. Scaling by 2 across and 3 up")
    print("-" * 70)
    S = scaling(shapes.SCALE_X, shapes.SCALE_Y)
    show(
        "scaling(2, 3)",
        S,
        "one step right becomes two steps right; one step up becomes three up",
    )
    assert S == shapes.SCALE_MATRIX
    print(f"    (1, 1) lands at {apply(S, (1.0, 1.0))}    checked by hand: (2, 3)")
    assert apply(S, (1.0, 1.0)) == (2.0, 3.0)

    # -- Reflection ------------------------------------------------------------
    print()
    print("2. Reflection in the x axis")
    print("-" * 70)
    F = reflection_in_x_axis()
    show(
        "reflection_in_x_axis()",
        F,
        "(1, 0) is ON the mirror line so it cannot move; (0, 1) mirrors to (0, -1)",
    )
    assert F == shapes.FLIP_MATRIX
    print(f"    (2, 3) lands at {apply(F, (2.0, 3.0))}   checked by hand: (2, -3)")
    assert apply(F, (2.0, 3.0)) == (2.0, -3.0)

    # -- Shear ----------------------------------------------------------------
    print()
    print("3. Shear: push sideways in proportion to height, k = 2")
    print("-" * 70)
    H = shear_x(shapes.SHEAR_K)
    show(
        "shear_x(2)",
        H,
        "(1, 0) has height 0 so nothing pushes it; (0, 1) has height 1 so it slides 2",
    )
    assert H == shapes.SHEAR_MATRIX
    print(f"    (1, 1) lands at {apply(H, (1.0, 1.0))}    checked by hand: (3, 1)")
    print(f"    (5, 0) lands at {apply(H, (5.0, 0.0))}    the x axis never moves")
    assert apply(H, (1.0, 1.0)) == (3.0, 1.0)
    assert apply(H, (5.0, 0.0)) == (5.0, 0.0)

    # -- Rotation --------------------------------------------------------------
    print()
    print("4. Rotation, derived from the unit circle")
    print("-" * 70)
    print("  Walk anticlockwise around a circle of radius 1, starting at (1, 0),")
    print("  until you have turned through theta. Where you now stand is, BY")
    print("  DEFINITION, (cos theta, sin theta). That is what those two")
    print("  functions are. So (1, 0) lands at (cos theta, sin theta), and")
    print("  (0, 1) -- already a quarter turn ahead -- lands a quarter turn")
    print("  ahead of that, at (-sin theta, cos theta).")
    print()
    for degrees in (30, 45, 90, 180):
        R = rotation(math.radians(degrees))
        print(f"  rotation({degrees} degrees):")
        print(f"    cos = {R[0][0]!r}")
        print(f"    sin = {R[1][0]!r}")
        print(f"    (1, 0) lands at ({R[0][0]!r}, {R[1][0]!r})")

    print()
    print("5. The quarter turn, and why an exact comparison would fail here")
    print("-" * 70)
    Q = rotation(math.pi / 2)
    e1_lands = apply(Q, shapes.E1)
    e2_lands = apply(Q, shapes.E2)
    print(f"  (1, 0) lands at {e1_lands!r}")
    print(f"  (0, 1) lands at {e2_lands!r}")
    print()
    print("  On paper those are (0, 1) and (-1, 0). In binary floating point")
    print(f"  cos(pi / 2) is {math.cos(math.pi / 2)!r}, which is not 0.0, because")
    print("  pi itself cannot be stored exactly and the cosine of the stored")
    print("  value is not the cosine of pi. The error is about 1e-17.")
    print(f"  So every check below uses a tolerance of {TOL!r}: five orders of")
    print("  magnitude above that rounding, and four below the smallest number")
    print("  this lab cares about. `== 0.0` would fail on a correct answer.")
    assert e1_lands != (0.0, 1.0), "exact equality really does fail here"
    assert abs(e1_lands[0] - 0.0) <= TOL and abs(e1_lands[1] - 1.0) <= TOL
    assert abs(e2_lands[0] - (-1.0)) <= TOL and abs(e2_lands[1] - 0.0) <= TOL

    # -- 6. NumPy agrees on all four ------------------------------------------
    print()
    print("6. NumPy builds the same four matrices")
    print("-" * 70)
    theta = math.pi / 2
    npQ = np.array(
        [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]
    )
    for name, mine, theirs in [
        ("scaling(2, 3)", S, np.diag([2.0, 3.0])),
        ("reflection", F, np.array([[1.0, 0.0], [0.0, -1.0]])),
        ("shear_x(2)", H, np.array([[1.0, 2.0], [0.0, 1.0]])),
        ("rotation(pi/2)", Q, npQ),
    ]:
        agree = np.allclose(np.array(mine), theirs, atol=TOL)
        print(f"  {name:<16} matches NumPy within {TOL!r}: {agree}")
        assert agree

    # -- 7. The flag, transformed ---------------------------------------------
    print()
    print("7. The flag shape under each transformation")
    print("-" * 70)
    print(f"  original      {[tuple(p) for p in shapes.FLAG]}")
    for name, M in [("scaled", S), ("reflected", F), ("sheared", H), ("turned", Q)]:
        moved = transform_polygon(M, shapes.FLAG)
        rounded = [(round(x, 6) + 0.0, round(y, 6) + 0.0) for x, y in moved]
        print(f"  {name:<13} {rounded}")
        # The corner at the origin never moves, under any of them. A linear
        # transformation always fixes the origin -- there is no matrix that
        # can move it, because M @ (0, 0) is (0, 0) whatever M holds.
        assert abs(moved[0][0]) <= TOL and abs(moved[0][1]) <= TOL
    print()
    print("  Every one of them left the first corner at the origin. That is not")
    print("  a coincidence and it is not avoidable: M @ (0, 0) is (0, 0) for")
    print("  every matrix M there has ever been. A linear transformation cannot")
    print("  move the origin, which is the first of the two limits this lab is")
    print("  really about.")

    print()
    print("02_building_the_transformations.py: every assertion held.")


if __name__ == "__main__":
    main()

"""What makes a transformation linear -- tested on one that is, and one that is not.

Run from inside examples/:

    ../.venv/bin/python3 03_linear_or_not.py

Linear means exactly two things, and nothing else:

    1. it preserves addition:  T(u + v) = T(u) + T(v)
    2. it preserves scaling:   T(s * u) = s * T(u)

Both must hold, for every u, v and s. This script checks both on a matrix,
where they hold, and on "multiply by a matrix and then add a constant", where
they do not -- and measures exactly how much they miss by, because the size of
the gap turns out to be the whole explanation.
"""

import numpy as np

import shapes
from transforms import apply, is_linear, preserves_addition, preserves_scaling

TOL = shapes.TOL

U = (1.0, 2.0)
V = (3.0, -1.0)
S = 5.0
OFFSET = (1.0, 1.0)


def main() -> None:
    print("03_linear_or_not.py")
    print("=" * 70)

    M = shapes.SCALE_MATRIX

    def linear(p):
        """v -> M @ v. Nothing else."""
        return apply(M, p)

    def affine(p):
        """v -> M @ v + b. One addition more, and no longer linear."""
        x, y = apply(M, p)
        return (x + OFFSET[0], y + OFFSET[1])

    print()
    print(f"  M = [{M[0]}, {M[1]}]   (scale by 2 across, 3 up)")
    print(f"  b = {OFFSET}")
    print(f"  u = {U}   v = {V}   s = {S}")

    # -- 1. The matrix passes both halves -------------------------------------
    print()
    print("1. T(v) = M @ v  --  preserves addition")
    print("-" * 70)
    ok, together, separately = preserves_addition(linear, U, V, TOL)
    print(f"  u + v          = {(U[0] + V[0], U[1] + V[1])}")
    print(f"  T(u + v)       = {together}")
    print(f"  T(u)           = {linear(U)}")
    print(f"  T(v)           = {linear(V)}")
    print(f"  T(u) + T(v)    = {separately}")
    print(f"  equal within {TOL!r}: {ok}")
    assert ok

    print()
    print("2. T(v) = M @ v  --  preserves scaling")
    print("-" * 70)
    ok, scaled_first, scaled_after = preserves_scaling(linear, U, S, TOL)
    print(f"  s * u          = {(S * U[0], S * U[1])}")
    print(f"  T(s * u)       = {scaled_first}")
    print(f"  s * T(u)       = {scaled_after}")
    print(f"  equal within {TOL!r}: {ok}")
    assert ok
    assert is_linear(linear, U, V, S, TOL)

    # -- 3. Adding a constant breaks it ---------------------------------------
    print()
    print("3. f(v) = M @ v + b  --  fails to preserve addition")
    print("-" * 70)
    ok_add, together, separately = preserves_addition(affine, U, V, TOL)
    print(f"  f(u + v)       = {together}")
    print(f"  f(u)           = {affine(U)}")
    print(f"  f(v)           = {affine(V)}")
    print(f"  f(u) + f(v)    = {separately}")
    gap_add = (separately[0] - together[0], separately[1] - together[1])
    print(f"  equal within {TOL!r}: {ok_add}")
    print(f"  the gap        = {gap_add}")
    print()
    print("  Look at the gap. It is exactly b. Adding the offset once on the")
    print("  left and twice on the right is the entire failure -- b sneaks in")
    print("  once per term, and adding two terms adds it twice.")
    assert not ok_add
    assert abs(gap_add[0] - OFFSET[0]) <= TOL
    assert abs(gap_add[1] - OFFSET[1]) <= TOL

    print()
    print("4. f(v) = M @ v + b  --  fails to preserve scaling too")
    print("-" * 70)
    ok_scale, scaled_first, scaled_after = preserves_scaling(affine, U, S, TOL)
    print(f"  f(s * u)       = {scaled_first}")
    print(f"  s * f(u)       = {scaled_after}")
    gap_scale = (scaled_after[0] - scaled_first[0], scaled_after[1] - scaled_first[1])
    print(f"  equal within {TOL!r}: {ok_scale}")
    print(f"  the gap        = {gap_scale}")
    print(f"  and (s - 1) * b = {((S - 1) * OFFSET[0], (S - 1) * OFFSET[1])}")
    print()
    print("  Same story: b is added once before the multiply and s times after.")
    assert not ok_scale
    assert abs(gap_scale[0] - (S - 1) * OFFSET[0]) <= TOL
    assert abs(gap_scale[1] - (S - 1) * OFFSET[1]) <= TOL
    assert not is_linear(affine, U, V, S, TOL)

    # -- 5. The one-second version of the same test ---------------------------
    print()
    print("5. The quick check: what happens to the origin?")
    print("-" * 70)
    print(f"  T((0, 0)) = {linear((0.0, 0.0))}     linear: the origin is fixed")
    print(f"  f((0, 0)) = {affine((0.0, 0.0))}     not linear: the origin moved")
    print()
    print("  Every linear transformation sends the origin to the origin, because")
    print("  M @ (0, 0) is a sum of zero lots of each column. So if a function")
    print("  moves the origin it cannot be linear, and you knew that before")
    print("  testing a single pair of vectors.")
    assert linear((0.0, 0.0)) == (0.0, 0.0)
    assert affine((0.0, 0.0)) == OFFSET

    # -- 6. Why a network keeps the bias separate -----------------------------
    print()
    print("6. Which is why a network layer is written X @ W + b")
    print("-" * 70)
    print("  The layer is deliberately NOT one operation. X @ W is the linear")
    print("  part -- it has a matrix, it has columns, it has a determinant, and")
    print("  everything in this lab applies to it. The + b is bolted on")
    print("  afterwards precisely because it cannot be folded into the matrix:")
    print("  a matrix cannot move the origin, and b exists to move the origin.")
    print("  Together they are called an AFFINE transformation: linear, plus a")
    print("  shift. Neither word is decoration.")

    # -- 7. NumPy, saying the same thing --------------------------------------
    print()
    print("7. The same four checks in NumPy")
    print("-" * 70)
    npM = np.array(M)
    npb = np.array(OFFSET)
    npu, npv = np.array(U), np.array(V)
    print(f"  M @ (u + v)             = {(npM @ (npu + npv)).tolist()}")
    print(f"  M @ u + M @ v           = {(npM @ npu + npM @ npv).tolist()}")
    print(f"  (M @ (u + v) + b)       = {(npM @ (npu + npv) + npb).tolist()}")
    print(f"  (M @ u + b) + (M @ v + b) = {(npM @ npu + npb + npM @ npv + npb).tolist()}")
    assert np.allclose(npM @ (npu + npv), npM @ npu + npM @ npv, atol=TOL)
    assert not np.allclose(
        npM @ (npu + npv) + npb, npM @ npu + npb + npM @ npv + npb, atol=TOL
    )

    print()
    print("03_linear_or_not.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Doing two transformations is one matrix -- and the order is not what it reads like.

Run from inside examples/:

    ../.venv/bin/python3 04_composition_and_order.py

Two claims, both checked against real numbers:

    1. shearing and then rotating a shape gives exactly the same answer as
       transforming it once by the product of the two matrices;
    2. that product is written ROTATE @ SHEAR, with the FIRST step on the
       RIGHT -- and writing it the other way round gives a different
       transformation, not a differently-spelled one.
"""

import math

import numpy as np

import shapes
from transforms import (
    apply,
    compose,
    determinant,
    rotation,
    shear_x,
    transform_polygon,
)

TOL = shapes.TOL


def rounded(points):
    return [(round(x, 6) + 0.0, round(y, 6) + 0.0) for x, y in points]


def main() -> None:
    print("04_composition_and_order.py")
    print("=" * 70)

    A = shear_x(shapes.SHEAR_K)          # step one
    B = rotation(math.pi / 2)            # step two

    print()
    print("  A = shear_x(2)        first")
    print(f"      [{A[0]}, {A[1]}]")
    print("  B = rotation(pi / 2)  second")
    print(f"      [[{B[0][0]:.1f}, {B[0][1]:.1f}], [{B[1][0]:.1f}, {B[1][1]:.1f}]]"
          "   (printed to one decimal; see script 02 for the raw values)")

    # -- 1. One step at a time -------------------------------------------------
    print()
    print("1. The two steps, one after the other, on the flag")
    print("-" * 70)
    step1 = transform_polygon(A, shapes.FLAG)
    step2 = transform_polygon(B, step1)
    print(f"  start        {rounded(shapes.FLAG)}")
    print(f"  after shear  {rounded(step1)}")
    print(f"  after turn   {rounded(step2)}")

    # -- 2. The single matrix that does both ----------------------------------
    print()
    print("2. Building the one matrix that does both, column by column")
    print("-" * 70)
    print("  Where does (1, 0) end up after BOTH steps?")
    after_e1 = apply(B, apply(A, shapes.E1))
    print(f"    shear sends (1, 0) to {apply(A, shapes.E1)}")
    print(f"    then the turn sends that to {rounded([after_e1])[0]}")
    print("  Where does (0, 1) end up?")
    after_e2 = apply(B, apply(A, shapes.E2))
    print(f"    shear sends (0, 1) to {apply(A, shapes.E2)}")
    print(f"    then the turn sends that to {rounded([after_e2])[0]}")
    print("  Write the two landings as columns and you have the composite.")

    C = compose(B, A)
    print(f"  compose(B, A) = [{rounded([tuple(C[0])])[0]}, {rounded([tuple(C[1])])[0]}]")
    print(f"  by hand:        [{shapes.SHEAR_THEN_ROTATE[0]}, {shapes.SHEAR_THEN_ROTATE[1]}]")
    for row_mine, row_hand in zip(C, shapes.SHEAR_THEN_ROTATE):
        for got, want in zip(row_mine, row_hand):
            assert abs(got - want) <= TOL, (got, want)

    # -- 3. One step gives the same answer as two -----------------------------
    print()
    print("3. One transformation, same landing places")
    print("-" * 70)
    at_once = transform_polygon(C, shapes.FLAG)
    print(f"  two steps    {rounded(step2)}")
    print(f"  one matrix   {rounded(at_once)}")
    for (x1, y1), (x2, y2) in zip(step2, at_once):
        assert abs(x1 - x2) <= TOL and abs(y1 - y2) <= TOL
    print(f"  every corner agrees within {TOL!r}")
    print()
    print("  This is why composition matters in practice. A stack of twenty")
    print("  transformations applied to a million points is twenty million")
    print("  operations; multiplying the twenty small matrices together first")
    print("  and applying one is a million. The answer is identical.")

    # -- 4. The order gotcha ---------------------------------------------------
    print()
    print("4. The order that trips everyone up")
    print("-" * 70)
    other = compose(A, B)   # rotate FIRST, then shear
    print(f"  compose(B, A)  shear then turn  = "
          f"[{rounded([tuple(C[0])])[0]}, {rounded([tuple(C[1])])[0]}]")
    print(f"  compose(A, B)  turn then shear  = "
          f"[{rounded([tuple(other[0])])[0]}, {rounded([tuple(other[1])])[0]}]")
    print(f"  by hand, turn then shear:         "
          f"[{shapes.ROTATE_THEN_SHEAR[0]}, {shapes.ROTATE_THEN_SHEAR[1]}]")
    for row_mine, row_hand in zip(other, shapes.ROTATE_THEN_SHEAR):
        for got, want in zip(row_mine, row_hand):
            assert abs(got - want) <= TOL, (got, want)

    differs = any(
        abs(a - b) > TOL for r1, r2 in zip(C, other) for a, b in zip(r1, r2)
    )
    print(f"  the two products differ: {differs}")
    assert differs
    probe = (1.0, 1.0)
    print(f"  and they send {probe} to different places:")
    print(f"    shear then turn -> {rounded([apply(C, probe)])[0]}")
    print(f"    turn then shear -> {rounded([apply(other, probe)])[0]}")
    assert rounded([apply(C, probe)])[0] != rounded([apply(other, probe)])[0]
    print()
    print("  Read a product RIGHT TO LEFT. B @ A means A first. It looks")
    print("  backwards until you write out B @ (A @ v): A is the one standing")
    print("  next to the vector, so A is the one that touches it first.")

    # -- 5. Determinants multiply ---------------------------------------------
    print()
    print("5. A free consequence: the area factors multiply")
    print("-" * 70)
    dA, dB, dC = determinant(A), determinant(B), determinant(C)
    print(f"  det(A) = {dA}    (a shear preserves area)")
    print(f"  det(B) = {dB}    (a turn preserves area)")
    print(f"  det(C) = {dC}")
    print(f"  det(A) * det(B) = {dA * dB}")
    assert abs(dC - dA * dB) <= TOL
    print("  Here both factors are 1, so the product is 1 and the composite")
    print("  preserves area as well. The rule is general and it is obvious once")
    print("  said out loud: if one step doubles area and the next triples it,")
    print("  area comes out six times bigger, whatever the matrices look like.")
    print("  A useful corollary: if either factor is 0, so is the product, so")
    print("  once a collapse has happened nothing downstream can undo it.")

    # -- 6. NumPy ---------------------------------------------------------------
    print()
    print("6. NumPy, with the @ operator")
    print("-" * 70)
    npA = np.array(A)
    npB = np.array(B)
    print(f"  B @ A = {np.round(npB @ npA, 6).tolist()}")
    print(f"  A @ B = {np.round(npA @ npB, 6).tolist()}")
    assert np.allclose(npB @ npA, np.array(C), atol=TOL)
    assert np.allclose(npA @ npB, np.array(other), atol=TOL)
    assert not np.allclose(npB @ npA, npA @ npB, atol=TOL)
    v = np.array([1.0, 1.0])
    assert np.allclose(npB @ (npA @ v), (npB @ npA) @ v, atol=TOL)
    print("  B @ (A @ v) equals (B @ A) @ v, which is the whole justification")
    print("  for composing matrices at all.")

    print()
    print("04_composition_and_order.py: every assertion held.")


if __name__ == "__main__":
    main()

"""The determinant as an area factor, the inverse, and what a collapse costs.

Run from inside examples/:

    ../.venv/bin/python3 05_determinant_inverse_rank.py

The determinant is introduced here the way it is actually useful: send the unit
square through the transformation and measure the area of what comes out. That
number, with its sign, IS the determinant. Everything else -- when an inverse
exists, what rank means, why a collapse is permanent -- reads off it.
"""

import numpy as np

import shapes
from transforms import (
    apply,
    compose,
    determinant,
    identity,
    inverse,
    rank,
    reflection_in_x_axis,
    scaling,
    shear_x,
    signed_area,
    transform_polygon,
    SingularMatrix,
)

TOL = shapes.TOL


def main() -> None:
    print("05_determinant_inverse_rank.py")
    print("=" * 70)

    print()
    print("1. The unit square, before anything happens to it")
    print("-" * 70)
    print(f"  corners (anticlockwise): {shapes.UNIT_SQUARE}")
    print(f"  signed area: {signed_area(shapes.UNIT_SQUARE)}")
    assert abs(signed_area(shapes.UNIT_SQUARE) - 1.0) <= TOL
    print("  Anticlockwise gives a POSITIVE area. List the same four corners")
    print("  the other way round and the shoelace formula returns -1. That sign")
    print("  is the thing to watch.")

    # -- 2. Measure, then compare against the determinant ---------------------
    print()
    print("2. Send it through four transformations and measure what comes out")
    print("-" * 70)
    cases = [
        ("scaling(2, 3)", scaling(2.0, 3.0), 6.0),
        ("shear_x(2)", shear_x(2.0), 1.0),
        ("reflection in x", reflection_in_x_axis(), -1.0),
        ("collapse", shapes.COLLAPSE_MATRIX, 0.0),
    ]
    print(f"  {'transformation':<18}{'measured area':>15}{'determinant':>14}"
          f"{'by hand':>10}")
    for name, M, by_hand in cases:
        moved = transform_polygon(M, shapes.UNIT_SQUARE)
        area = signed_area(moved)
        det = determinant(M)
        print(f"  {name:<18}{area:>15}{det:>14}{by_hand:>10}")
        assert abs(area - det) <= TOL, (name, area, det)
        assert abs(det - by_hand) <= TOL, (name, det, by_hand)

    print()
    print("  The measured area and the determinant are the same number every")
    print("  time, sign included. The determinant is not a formula that happens")
    print("  to be useful; it is the area factor, computed without drawing.")

    # -- 3. What the sign means -----------------------------------------------
    print()
    print("3. What a NEGATIVE determinant means")
    print("-" * 70)
    F = reflection_in_x_axis()
    flipped = transform_polygon(F, shapes.UNIT_SQUARE)
    print(f"  before: {shapes.UNIT_SQUARE}")
    print(f"  after : {flipped}")
    print(f"  signed area went from {signed_area(shapes.UNIT_SQUARE)} to "
          f"{signed_area(flipped)}")
    print("  The size did not change -- the square is still area 1. What")
    print("  changed is that the corners now run clockwise. The plane was")
    print("  turned over, and no amount of rotating will turn it back, the same")
    print("  way no amount of turning a left glove makes it a right one.")
    assert signed_area(flipped) < 0
    assert abs(abs(signed_area(flipped)) - 1.0) <= TOL

    # -- 4. What zero means ----------------------------------------------------
    print()
    print("4. What a ZERO determinant means")
    print("-" * 70)
    G = shapes.COLLAPSE_MATRIX
    print(f"  G = [{G[0]}, {G[1]}]")
    print("  column 0 = (1, 2)   column 1 = (2, 4)")
    print("  The second column is exactly twice the first. Both basis vectors")
    print("  land on the SAME LINE through the origin, so everything else does")
    print("  too -- there is nowhere else left to land.")
    landed = [apply(G, p) for p in [(1.0, 0.0), (0.0, 1.0), (3.0, -1.0), (7.0, 7.0)]]
    for start, end in zip([(1.0, 0.0), (0.0, 1.0), (3.0, -1.0), (7.0, 7.0)], landed):
        on_line = abs(end[1] - 2.0 * end[0]) <= TOL
        print(f"    {str(start):<12} -> {str(end):<14} on the line y = 2x: {on_line}")
        assert on_line
    squashed = transform_polygon(G, shapes.UNIT_SQUARE)
    print(f"  the unit square becomes {squashed}")
    print(f"  its area is {signed_area(squashed)}")
    assert abs(signed_area(squashed)) <= TOL
    print()
    print("  Two different starting points now share a landing place:")
    print(f"    (2, 0) -> {apply(G, (2.0, 0.0))}")
    print(f"    (0, 1) -> {apply(G, (0.0, 1.0))}")
    assert apply(G, (2.0, 0.0)) == apply(G, (0.0, 1.0))
    print("  No rule can send that shared place back to both of them. The")
    print("  information is not hidden, it is gone, and that is exactly what")
    print("  'no inverse' means.")

    # -- 5. Rank ---------------------------------------------------------------
    print()
    print("5. Rank: how many dimensions survive")
    print("-" * 70)
    for name, M in [
        ("identity", identity()),
        ("scaling(2, 3)", scaling(2.0, 3.0)),
        ("shear_x(2)", shear_x(2.0)),
        ("collapse", G),
        ("everything to the origin", [[0.0, 0.0], [0.0, 0.0]]),
    ]:
        mine = rank(M)
        theirs = int(np.linalg.matrix_rank(np.array(M)))
        print(f"  {name:<26} rank {mine}   numpy.linalg.matrix_rank: {theirs}")
        assert mine == theirs
    print()
    print("  Rank 2: the output still fills the plane. Rank 1: it is squashed")
    print("  onto a line. Rank 0: everything lands on the origin. For a square")
    print("  matrix, full rank and a non-zero determinant are the same sentence.")
    assert rank(G) == shapes.COLLAPSE_RANK

    # -- 6. The inverse --------------------------------------------------------
    print()
    print("6. The inverse, where one exists")
    print("-" * 70)
    H = shear_x(2.0)
    Hinv = inverse(H)
    print(f"  shear_x(2)          [{H[0]}, {H[1]}]")
    print(f"  its inverse         [{Hinv[0]}, {Hinv[1]}]")
    print("  which is shear_x(-2) -- push the deck of cards back the other way.")
    back = compose(Hinv, H)
    print(f"  inverse @ original = [{back[0]}, {back[1]}]   the identity")
    for row, want in zip(back, identity()):
        for got, expect in zip(row, want):
            assert abs(got - expect) <= TOL
    probe = (1.0, 1.0)
    there = apply(H, probe)
    home = apply(Hinv, there)
    print(f"  {probe} -> {there} -> {home}")
    assert abs(home[0] - probe[0]) <= TOL and abs(home[1] - probe[1]) <= TOL

    S = scaling(2.0, 3.0)
    Sinv = inverse(S)
    print(f"  scaling(2, 3) inverse: [{Sinv[0]}, {Sinv[1]}]")
    print("  which is scaling by 1/2 and 1/3, as it must be.")
    assert abs(Sinv[0][0] - 0.5) <= TOL
    assert abs(Sinv[1][1] - 1.0 / 3.0) <= TOL

    # -- 7. Asking for the impossible ------------------------------------------
    print()
    print("7. Asking for the inverse of a collapse")
    print("-" * 70)
    try:
        inverse(G)
    except SingularMatrix as exc:
        print(f"  from-scratch inverse raises SingularMatrix: {exc}")
    else:  # pragma: no cover - the assert below turns this into a failure
        raise AssertionError("the from-scratch inverse should have refused")

    try:
        np.linalg.inv(np.array(G))
    except np.linalg.LinAlgError as exc:
        print(f"  numpy.linalg.inv raises {type(exc).__name__}: {exc}")
    else:  # pragma: no cover
        raise AssertionError("numpy should have refused too")

    print()
    print("  Both refuse, and both refusals are catchable as ValueError:")
    print(f"    SingularMatrix is a ValueError:      "
          f"{issubclass(SingularMatrix, ValueError)}")
    print(f"    numpy.linalg.LinAlgError is one too: "
          f"{issubclass(np.linalg.LinAlgError, ValueError)}")
    assert issubclass(SingularMatrix, ValueError)
    assert issubclass(np.linalg.LinAlgError, ValueError)

    # -- 8. Where the two determinants disagree, and why ----------------------
    print()
    print("8. An honest difference between the two determinants")
    print("-" * 70)
    P = shapes.PICTURE_MATRIX
    mine = determinant(P)
    theirs = float(np.linalg.det(np.array(P)))
    print(f"  P = [{P[0]}, {P[1]}]")
    print(f"  by hand:            3 * 2 - (-1) * 1 = {shapes.PICTURE_DETERMINANT}")
    print(f"  from-scratch:       {mine!r}")
    print(f"  numpy.linalg.det:   {theirs!r}")
    print(f"  they differ by      {abs(mine - theirs)!r}")
    print()
    print("  This is not a bug in either one. The from-scratch version computes")
    print("  a*d - b*c directly, which for these four whole numbers is exact.")
    print("  numpy.linalg.det factorises the matrix first -- the same routine")
    print("  it uses for a 500 by 500 matrix, where the direct formula is not")
    print("  an option -- and that factorisation rounds. The general method")
    print("  pays a little accuracy on tiny inputs to stay usable on large")
    print("  ones. It is a good trade and it is worth knowing about, because")
    print("  it is why you compare determinants with a tolerance rather than")
    print("  with ==.")
    assert abs(mine - shapes.PICTURE_DETERMINANT) <= TOL
    assert abs(theirs - shapes.PICTURE_DETERMINANT) <= 1e-9
    assert mine == shapes.PICTURE_DETERMINANT

    print()
    print("05_determinant_inverse_rank.py: every assertion held.")


if __name__ == "__main__":
    main()

"""Why "linear" is a limitation, and why activation functions are not optional.

Run from inside examples/:

    ../.venv/bin/python3 06_the_limit_of_linear.py

Three facts, each demonstrated rather than asserted in prose:

    1. a linear transformation always fixes the origin;
    2. it always sends straight lines to straight lines -- the midpoint of two
       points lands on the midpoint of their landing places, every time;
    3. a stack of linear transformations, however deep, collapses to ONE
       matrix, so the stack can do nothing a single layer could not.

Put together, they say something concrete about neural networks: twenty
matrix layers with nothing between them are exactly one matrix layer, and no
matrix can draw a curved boundary. That is the reason a non-linear activation
sits between the layers, and it is the best thing this lab has to offer.
"""

import math
import random

import numpy as np

import shapes
from transforms import apply, compose, rotation, scaling, shear_x, transform_polygon

TOL = shapes.TOL


def main() -> None:
    print("06_the_limit_of_linear.py")
    print("=" * 70)

    M = compose(rotation(math.radians(30)), shear_x(1.5))

    # -- 1. The origin cannot move --------------------------------------------
    print()
    print("1. The origin never moves")
    print("-" * 70)
    print(f"  M @ (0, 0) = {apply(M, (0.0, 0.0))}")
    print("  M @ (0, 0) is 0 lots of column 0 plus 0 lots of column 1, which is")
    print("  (0, 0) whatever the columns hold. There is no matrix anywhere that")
    print("  moves the origin. If your data needs shifting, that shift has to")
    print("  arrive from somewhere else -- and in a network layer it does,")
    print("  under the name b.")
    landed = apply(M, (0.0, 0.0))
    assert abs(landed[0]) <= TOL and abs(landed[1]) <= TOL

    # -- 2. Straight stays straight -------------------------------------------
    print()
    print("2. Straight lines land as straight lines")
    print("-" * 70)
    p = (1.0, 3.0)
    q = (4.0, -2.0)
    midpoint = ((p[0] + q[0]) / 2, (p[1] + q[1]) / 2)
    print(f"  p = {p}   q = {q}   midpoint of p and q = {midpoint}")
    landed_mid = apply(M, midpoint)
    mid_of_landed = tuple(
        (a + b) / 2 for a, b in zip(apply(M, p), apply(M, q))
    )
    print(f"  M @ midpoint            = {tuple(round(c, 12) for c in landed_mid)}")
    print(f"  midpoint of M@p and M@q = {tuple(round(c, 12) for c in mid_of_landed)}")
    print(f"  the same within {TOL!r}: "
          f"{all(abs(a - b) <= TOL for a, b in zip(landed_mid, mid_of_landed))}")
    assert all(abs(a - b) <= TOL for a, b in zip(landed_mid, mid_of_landed))

    print()
    print("  And not just the midpoint. Eleven points spaced evenly along the")
    print("  line from p to q stay evenly spaced along a line after the transform:")
    worst = 0.0
    for i in range(11):
        t = i / 10
        on_line = (p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1]))
        expected = tuple(
            a + t * (b - a) for a, b in zip(apply(M, p), apply(M, q))
        )
        got = apply(M, on_line)
        worst = max(worst, max(abs(a - b) for a, b in zip(got, expected)))
    print(f"  worst disagreement across all eleven points: {worst!r}")
    assert worst <= TOL
    print("  Evenly spaced in, evenly spaced out. A transformation that cannot")
    print("  bend a line cannot draw a curve, and a boundary that needs a curve")
    print("  is therefore out of reach -- not hard, out of reach.")

    # -- 3. A stack collapses --------------------------------------------------
    print()
    print("3. A stack of linear layers is one linear layer")
    print("-" * 70)
    random.seed(102)
    stack = [
        [[random.uniform(-2, 2) for _ in range(2)] for _ in range(2)]
        for _ in range(20)
    ]
    print("  Twenty 2 by 2 matrices with pseudo-random entries (seed 102),")
    print("  applied one after another to a point.")

    point = (0.7, -0.4)
    one_at_a_time = point
    for layer in stack:
        one_at_a_time = apply(layer, one_at_a_time)

    combined = stack[0]
    for layer in stack[1:]:
        combined = compose(layer, combined)
    all_at_once = apply(combined, point)

    print(f"  twenty applications:      {one_at_a_time}")
    print(f"  one combined matrix:      {all_at_once}")
    print(f"  the combined matrix is    [{combined[0]}, {combined[1]}]")
    relative = max(
        abs(a - b) / max(1.0, abs(a)) for a, b in zip(one_at_a_time, all_at_once)
    )
    print(f"  largest relative difference: {relative!r}")
    print()
    print("  A wider tolerance is used here on purpose, and saying why is more")
    print("  useful than hiding it: these entries are not small whole numbers,")
    print("  the two routes multiply them in a different ORDER, and floating")
    print("  point addition is not associative, so twenty layers of rounding")
    print("  accumulate. On this run the two routes agreed to within a")
    print("  relative 1e-15 -- close, but not to the last bit, and demanding")
    print("  the last bit would be demanding something arithmetic does not")
    print("  promise. The check allows 1e-9, which is loose enough to stay")
    print("  true on another machine and tight enough to catch a real error.")
    assert relative <= 1e-9

    print()
    print("  The point stands regardless of the last few digits: twenty layers")
    print("  did nothing that one 2 by 2 matrix could not do. Depth bought")
    print("  nothing. Stack a thousand and it is still one matrix.")

    # -- 4. The consequence, made concrete -------------------------------------
    print()
    print("4. The data that no stack of these can separate")
    print("-" * 70)
    print("  Four points, the classic exclusive-or arrangement:")
    print("    (0, 0) -> class A      (1, 1) -> class A")
    print("    (1, 0) -> class B      (0, 1) -> class B")
    print()
    print("  Class A sits on one diagonal and class B on the other. No straight")
    print("  line separates them -- try it on paper, it takes about ten seconds")
    print("  to convince yourself. A linear transformation followed by a")
    print("  threshold can only ever cut the plane with a straight line, and a")
    print("  stack of them still cuts with a straight line, because point 3")
    print("  says the stack IS one transformation.")
    print()
    print("  Notice that (0, 0) is one of the four. It cannot be moved at all")
    print("  by any matrix, so it is not even a matter of finding good weights.")
    a_points = [(0.0, 0.0), (1.0, 1.0)]
    b_points = [(1.0, 0.0), (0.0, 1.0)]
    for M_try in [scaling(3.0, -2.0), shear_x(4.0), rotation(1.0), combined]:
        moved_a = [apply(M_try, p) for p in a_points]
        assert abs(moved_a[0][0]) <= TOL and abs(moved_a[0][1]) <= TOL
    print("  Checked against four different matrices: (0, 0) stayed at (0, 0)")
    print("  in every one.")

    print()
    print("5. Which is what the activation function is for")
    print("-" * 70)
    print("  Put a function that is NOT linear between the layers -- ReLU, which")
    print("  replaces every negative number with zero, is the usual one -- and")
    print("  the collapse in point 3 stops working, because you can no longer")
    print("  slide the matrices together past it. Depth starts to buy something.")
    print()
    relu_out = [max(0.0, c) for c in apply(M, (1.0, -1.0))]
    print(f"  M @ (1, -1)            = {tuple(round(c, 6) for c in apply(M, (1.0, -1.0)))}")
    print(f"  ReLU of that           = {tuple(round(c, 6) for c in relu_out)}")
    u, v = (1.0, -1.0), (0.5, 2.0)
    relu = lambda pt: tuple(max(0.0, c) for c in apply(M, pt))
    together = relu((u[0] + v[0], u[1] + v[1]))
    separately = tuple(a + b for a, b in zip(relu(u), relu(v)))
    print(f"  relu(M @ (u + v))      = {tuple(round(c, 6) for c in together)}")
    print(f"  relu(M@u) + relu(M@v)  = {tuple(round(c, 6) for c in separately)}")
    print("  Not equal -- so this composite is not linear, so it is not a")
    print("  matrix, so no amount of algebra folds the layers together. The")
    print("  non-linearity is doing the one job the matrix cannot.")
    assert any(abs(a - b) > TOL for a, b in zip(together, separately))

    # -- 6. NumPy -------------------------------------------------------------
    print()
    print("6. NumPy says the same about the stack")
    print("-" * 70)
    npstack = [np.array(layer) for layer in stack]
    npcombined = npstack[0]
    for layer in npstack[1:]:
        npcombined = layer @ npcombined
    print(f"  the product of all twenty: {np.round(npcombined, 6).tolist()}")
    assert np.allclose(npcombined, np.array(combined), rtol=1e-9, atol=0.0)
    print("  agrees with the from-scratch product to a relative 1e-9")

    square = transform_polygon(combined, shapes.UNIT_SQUARE)
    print(f"  the unit square under the whole stack: "
          f"{[(round(x, 4), round(y, 4)) for x, y in square]}")
    print("  still a parallelogram -- four straight edges, opposite sides")
    print("  parallel, one corner nailed to the origin. Twenty layers deep.")

    print()
    print("06_the_limit_of_linear.py: every assertion held.")


if __name__ == "__main__":
    main()

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

Both `examples/` and `starter/` contain modules called `transforms` and
`shapes`, and pytest imports test files by putting their directory on
`sys.path`. Without this file, running `pytest` across both directories at once
would import whichever `transforms` 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
`transforms` or `shapes` that came from somewhere else.
"""

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 ("transforms", "shapes", "answers"):
    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 invented data this lab works on, and the answers worked out by hand.

Everything here is made up. There is no dataset, no download and no file to
read: a 2 by 2 matrix has four numbers in it, and the whole point of the day is
that you can check every one of them on paper.

Points are `(x, y)` pairs of plain Python numbers. Polygons are lists of points
in counter-clockwise order, because the SIGN of a polygon's area depends on
which way round you list its corners, and that sign is what tells you whether a
transformation flipped the plane over.
"""

# -- The standard basis -------------------------------------------------------

E1 = (1.0, 0.0)
E2 = (0.0, 1.0)

# -- The unit square, listed counter-clockwise from the origin ----------------
#
# Corners: origin, one step right, the far corner, one step up. Its area is 1
# and its signed area is +1. After a transformation, the signed area of the
# image is exactly the determinant of the matrix -- sign included.

UNIT_SQUARE = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]

# -- The flag ------------------------------------------------------------------
#
# A deliberately lopsided shape, so you can see at a glance whether it has been
# turned, stretched, sheared or mirrored. It is an L lying on its back: a long
# foot along the x axis and a short mast going up at the left.

FLAG = [
    (0.0, 0.0),
    (2.0, 0.0),
    (2.0, 0.5),
    (0.5, 0.5),
    (0.5, 2.0),
    (0.0, 2.0),
]

# -- The picture you read a matrix off of --------------------------------------
#
# Exercise 1 describes a drawing in words rather than showing it, because the
# skill being trained is going from "where did the basis vectors land" to "what
# is the matrix", and a picture would let you skip the step. The drawing shows:
#
#     the arrow (1, 0) has been redrawn ending at (3, 1)
#     the arrow (0, 1) has been redrawn ending at (-1, 2)
#
# Those two landing places, written as COLUMNS, are the matrix.

PICTURE_E1_LANDS_AT = (3.0, 1.0)
PICTURE_E2_LANDS_AT = (-1.0, 2.0)

PICTURE_MATRIX = [
    [3.0, -1.0],
    [1.0, 2.0],
]

# Worked by hand, so the test has something honest to compare against.
#
#   (2, 1) = 2 * (1, 0) + 1 * (0, 1)
#   so it must land at 2 * (3, 1) + 1 * (-1, 2)
#                    = (6, 2) + (-1, 2)
#                    = (5, 4)
#
# and the determinant is 3 * 2 - (-1) * 1 = 6 + 1 = 7.
PICTURE_SENDS_2_1_TO = (5.0, 4.0)
PICTURE_DETERMINANT = 7.0

# -- The four standard transformations, with their hand-worked matrices --------

SCALE_X, SCALE_Y = 2.0, 3.0
SCALE_MATRIX = [[2.0, 0.0], [0.0, 3.0]]      # e1 -> (2, 0), e2 -> (0, 3)

FLIP_MATRIX = [[1.0, 0.0], [0.0, -1.0]]      # reflection in the x axis

SHEAR_K = 2.0
SHEAR_MATRIX = [[1.0, 2.0], [0.0, 1.0]]      # e1 stays, e2 -> (2, 1)

# Rotation by a quarter turn. Written exactly here; computed from cosine and
# sine in transforms.py, where it comes out as 6.123233995736766e-17 rather
# than 0. That difference is the reason every float check in this lab states a
# tolerance.
QUARTER_TURN_MATRIX = [[0.0, -1.0], [1.0, 0.0]]

# -- Composition ---------------------------------------------------------------
#
# Shear first, then rotate a quarter turn. In matrix form that is ROT @ SHEAR,
# with the FIRST transformation written on the RIGHT, because it is the one
# standing next to the vector.
#
#   ROT @ SHEAR = [[0, -1], [1, 0]] @ [[1, 2], [0, 1]]
#               = [[0*1 + -1*0, 0*2 + -1*1],
#                  [1*1 +  0*0, 1*2 +  0*1]]
#               = [[0, -1], [1, 2]]
SHEAR_THEN_ROTATE = [[0.0, -1.0], [1.0, 2.0]]

# The other order, to show that it is a different transformation entirely.
#
#   SHEAR @ ROT = [[1, 2], [0, 1]] @ [[0, -1], [1, 0]]
#               = [[1*0 + 2*1, 1*-1 + 2*0],
#                  [0*0 + 1*1, 0*-1 + 1*0]]
#               = [[2, -1], [1, 0]]
ROTATE_THEN_SHEAR = [[2.0, -1.0], [1.0, 0.0]]

# -- The collapse --------------------------------------------------------------
#
# Both columns point along the same line: (2, 4) is exactly twice (1, 2). So
# every vector in the plane lands somewhere on the line through (1, 2), the
# whole plane is squashed onto a line, and the area of anything you send
# through is zero. Nothing can undo it, because everything on that line came
# from a whole line's worth of starting points.
COLLAPSE_MATRIX = [[1.0, 2.0], [2.0, 4.0]]
COLLAPSE_DETERMINANT = 0.0                    # 1 * 4 - 2 * 2
COLLAPSE_RANK = 1

# -- Tolerance -----------------------------------------------------------------
#
# Why 1e-12 and not equality: cos(pi / 2) is 6.123233995736766e-17 in binary
# floating point, not 0.0, and sin(pi / 6) is 0.49999999999999994, not 0.5.
# Both are about 1e-17 away from the exact answer. 1e-12 sits five orders of
# magnitude above that error and about four below the smallest quantity this
# lab cares about (0.5), so it accepts the rounding and would still catch a
# genuinely wrong answer.
TOL = 1e-12

"""The reference test suite: real values, real shapes, real exceptions.

Run from the LAB DIRECTORY:

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

Every float comparison here states a tolerance and the module docstring of
shapes.py says why that number was chosen. Nothing is compared with == unless
the arithmetic that produced it is exact -- and where it is exact, the test
says == on purpose, because that is a stronger claim.
"""

import math

import numpy as np
import pytest

import shapes
from transforms import (
    SingularMatrix,
    apply,
    columns_of,
    compose,
    determinant,
    from_landings,
    identity,
    inverse,
    is_linear,
    preserves_addition,
    preserves_scaling,
    rank,
    reflection_in_x_axis,
    reflection_in_y_axis,
    rotation,
    scaling,
    shear_x,
    shear_y,
    signed_area,
    transform_polygon,
)

TOL = shapes.TOL


# -- Reading a matrix off its landings ---------------------------------------


def test_from_landings_builds_the_picture_matrix():
    M = from_landings(shapes.PICTURE_E1_LANDS_AT, shapes.PICTURE_E2_LANDS_AT)
    assert M == shapes.PICTURE_MATRIX


def test_columns_are_the_landings_not_the_rows():
    e1, e2 = columns_of(shapes.PICTURE_MATRIX)
    assert e1 == (3.0, 1.0)
    assert e2 == (-1.0, 2.0)
    # The trap this test exists for: row 0 is (3, -1) and is NOT a landing.
    assert tuple(shapes.PICTURE_MATRIX[0]) != e1


def test_the_matrix_sends_2_1_where_hand_arithmetic_says():
    assert apply(shapes.PICTURE_MATRIX, (2.0, 1.0)) == shapes.PICTURE_SENDS_2_1_TO


def test_numpy_agrees_on_the_picture_matrix():
    M = np.array(shapes.PICTURE_MATRIX)
    assert np.allclose(M @ np.array([2.0, 1.0]), np.array([5.0, 4.0]), atol=TOL)
    assert np.allclose(M @ np.array(shapes.E1), np.array([3.0, 1.0]), atol=TOL)
    assert np.allclose(M @ np.array(shapes.E2), np.array([-1.0, 2.0]), atol=TOL)


def test_applying_to_a_basis_vector_just_reads_a_column():
    for M in (shapes.PICTURE_MATRIX, shapes.SCALE_MATRIX, shapes.SHEAR_MATRIX):
        e1, e2 = columns_of(M)
        assert apply(M, shapes.E1) == e1
        assert apply(M, shapes.E2) == e2


# -- The identity -------------------------------------------------------------


def test_identity_leaves_everything_alone():
    I = identity()
    for point in [(0.0, 0.0), (1.0, 0.0), (-3.5, 7.25), (1e6, -1e-6)]:
        assert apply(I, point) == point


def test_identity_matches_numpy_eye():
    assert np.allclose(np.array(identity()), np.eye(2), atol=TOL)


def test_identity_has_determinant_one_and_full_rank():
    assert determinant(identity()) == 1.0
    assert rank(identity()) == 2


# -- The four standard transformations ---------------------------------------


def test_scaling_matrix_is_derived_correctly():
    assert scaling(2.0, 3.0) == shapes.SCALE_MATRIX
    assert apply(scaling(2.0, 3.0), (1.0, 1.0)) == (2.0, 3.0)
    assert apply(scaling(2.0, 3.0), (-4.0, 0.5)) == (-8.0, 1.5)


def test_reflection_in_x_axis_fixes_the_axis_and_flips_the_rest():
    F = reflection_in_x_axis()
    assert F == shapes.FLIP_MATRIX
    assert apply(F, (5.0, 0.0)) == (5.0, 0.0)
    assert apply(F, (2.0, 3.0)) == (2.0, -3.0)


def test_reflection_in_y_axis():
    F = reflection_in_y_axis()
    assert apply(F, (0.0, 5.0)) == (0.0, 5.0)
    assert apply(F, (2.0, 3.0)) == (-2.0, 3.0)


def test_reflecting_twice_is_the_identity():
    F = reflection_in_x_axis()
    assert compose(F, F) == identity()


def test_shear_leaves_the_x_axis_alone_and_slides_the_rest():
    H = shear_x(2.0)
    assert H == shapes.SHEAR_MATRIX
    assert apply(H, (5.0, 0.0)) == (5.0, 0.0)
    assert apply(H, (1.0, 1.0)) == (3.0, 1.0)
    assert apply(H, (0.0, 1.0)) == (2.0, 1.0)


def test_shear_y_is_the_same_idea_the_other_way_up():
    V = shear_y(3.0)
    assert apply(V, (0.0, 5.0)) == (0.0, 5.0)
    assert apply(V, (1.0, 0.0)) == (1.0, 3.0)


@pytest.mark.parametrize(
    "degrees, expected_cos, expected_sin",
    [
        (0, 1.0, 0.0),
        (30, math.sqrt(3) / 2, 0.5),
        (45, math.sqrt(2) / 2, math.sqrt(2) / 2),
        (90, 0.0, 1.0),
        (180, -1.0, 0.0),
    ],
)
def test_rotation_columns_are_the_unit_circle_coordinates(
    degrees, expected_cos, expected_sin
):
    R = rotation(math.radians(degrees))
    e1, e2 = columns_of(R)
    assert abs(e1[0] - expected_cos) <= TOL
    assert abs(e1[1] - expected_sin) <= TOL
    assert abs(e2[0] + expected_sin) <= TOL
    assert abs(e2[1] - expected_cos) <= TOL


def test_a_quarter_turn_sends_the_basis_where_the_picture_says():
    Q = rotation(math.pi / 2)
    x, y = apply(Q, shapes.E1)
    assert abs(x - 0.0) <= TOL and abs(y - 1.0) <= TOL
    x, y = apply(Q, shapes.E2)
    assert abs(x + 1.0) <= TOL and abs(y - 0.0) <= TOL


def test_the_quarter_turn_is_NOT_exactly_zero_which_is_why_tolerance_exists():
    """The honest reason this lab never uses == on a rotation.

    cos(pi / 2) is 6.123233995736766e-17 rather than 0.0, because pi cannot be
    stored exactly in binary and the cosine of the stored value is not the
    cosine of pi. This test asserts the inexactness itself, so that if a future
    NumPy or libm ever made it exact, the suite would say so rather than
    silently keeping a comment that had stopped being true.
    """
    assert math.cos(math.pi / 2) != 0.0
    assert 0.0 < abs(math.cos(math.pi / 2)) < 1e-15
    assert apply(rotation(math.pi / 2), shapes.E1) != (0.0, 1.0)


def test_sin_of_thirty_degrees_is_also_not_exactly_a_half():
    assert math.sin(math.radians(30)) != 0.5
    assert abs(math.sin(math.radians(30)) - 0.5) <= TOL


def test_four_quarter_turns_return_to_the_start():
    Q = rotation(math.pi / 2)
    back = compose(Q, compose(Q, compose(Q, Q)))
    for row, want in zip(back, identity()):
        for got, expect in zip(row, want):
            assert abs(got - expect) <= TOL


def test_rotation_preserves_length():
    R = rotation(0.9)
    for point in [(3.0, 4.0), (1.0, 0.0), (-2.0, 7.0)]:
        before = math.hypot(*point)
        after = math.hypot(*apply(R, point))
        assert abs(before - after) <= TOL


def test_all_four_match_numpy():
    theta = math.pi / 2
    pairs = [
        (scaling(2.0, 3.0), np.diag([2.0, 3.0])),
        (reflection_in_x_axis(), np.array([[1.0, 0.0], [0.0, -1.0]])),
        (shear_x(2.0), np.array([[1.0, 2.0], [0.0, 1.0]])),
        (
            rotation(theta),
            np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]),
        ),
    ]
    for mine, theirs in pairs:
        assert np.allclose(np.array(mine), theirs, atol=TOL)


# -- Linearity ----------------------------------------------------------------

U = (1.0, 2.0)
V = (3.0, -1.0)
S = 5.0
B = (1.0, 1.0)


def _linear(point):
    return apply(shapes.SCALE_MATRIX, point)


def _affine(point):
    x, y = apply(shapes.SCALE_MATRIX, point)
    return (x + B[0], y + B[1])


def test_a_matrix_preserves_addition():
    ok, together, separately = preserves_addition(_linear, U, V, TOL)
    assert ok
    assert together == (8.0, 3.0)
    assert separately == (8.0, 3.0)


def test_a_matrix_preserves_scaling():
    ok, first, after = preserves_scaling(_linear, U, S, TOL)
    assert ok
    assert first == (10.0, 30.0)
    assert after == (10.0, 30.0)


def test_a_matrix_is_linear():
    assert is_linear(_linear, U, V, S, TOL)


def test_adding_a_constant_breaks_addition_by_exactly_b():
    ok, together, separately = preserves_addition(_affine, U, V, TOL)
    assert not ok
    assert together == (9.0, 4.0)
    assert separately == (10.0, 5.0)
    gap = (separately[0] - together[0], separately[1] - together[1])
    assert gap == B


def test_adding_a_constant_breaks_scaling_by_exactly_s_minus_one_times_b():
    ok, first, after = preserves_scaling(_affine, U, S, TOL)
    assert not ok
    assert first == (11.0, 31.0)
    assert after == (15.0, 35.0)
    gap = (after[0] - first[0], after[1] - first[1])
    assert gap == ((S - 1) * B[0], (S - 1) * B[1])


def test_an_affine_function_is_not_linear():
    assert not is_linear(_affine, U, V, S, TOL)


def test_linear_fixes_the_origin_and_affine_does_not():
    assert _linear((0.0, 0.0)) == (0.0, 0.0)
    assert _affine((0.0, 0.0)) == B


def test_every_matrix_in_this_lab_fixes_the_origin():
    matrices = [
        identity(),
        scaling(2.0, 3.0),
        shear_x(2.0),
        reflection_in_x_axis(),
        rotation(1.234),
        shapes.PICTURE_MATRIX,
        shapes.COLLAPSE_MATRIX,
    ]
    for M in matrices:
        x, y = apply(M, (0.0, 0.0))
        assert abs(x) <= TOL and abs(y) <= TOL


def test_squaring_a_coordinate_is_not_linear():
    """A second non-linear example, so the point is not just about b."""

    def squarer(point):
        return (point[0] ** 2, point[1])

    assert squarer((0.0, 0.0)) == (0.0, 0.0)  # fixes the origin, yet still not linear
    assert not is_linear(squarer, U, V, S, TOL)


def test_a_linear_map_sends_midpoints_to_midpoints():
    M = compose(rotation(math.radians(30)), shear_x(1.5))
    p, q = (1.0, 3.0), (4.0, -2.0)
    mid = ((p[0] + q[0]) / 2, (p[1] + q[1]) / 2)
    landed_mid = apply(M, mid)
    mid_of_landed = tuple((a + b) / 2 for a, b in zip(apply(M, p), apply(M, q)))
    assert all(abs(a - b) <= TOL for a, b in zip(landed_mid, mid_of_landed))


# -- Composition --------------------------------------------------------------


def test_composition_matches_the_hand_worked_product():
    C = compose(rotation(math.pi / 2), shear_x(2.0))
    for row, want in zip(C, shapes.SHEAR_THEN_ROTATE):
        for got, expect in zip(row, want):
            assert abs(got - expect) <= TOL


def test_the_other_order_is_a_different_transformation():
    A, Bm = shear_x(2.0), rotation(math.pi / 2)
    for row, want in zip(compose(A, Bm), shapes.ROTATE_THEN_SHEAR):
        for got, expect in zip(row, want):
            assert abs(got - expect) <= TOL
    probe = (1.0, 1.0)
    one = apply(compose(Bm, A), probe)
    other = apply(compose(A, Bm), probe)
    assert any(abs(a - b) > TOL for a, b in zip(one, other))


def test_one_matrix_reproduces_the_two_step_sequence_on_every_corner():
    A, Bm = shear_x(2.0), rotation(math.pi / 2)
    two_steps = transform_polygon(Bm, transform_polygon(A, shapes.FLAG))
    at_once = transform_polygon(compose(Bm, A), shapes.FLAG)
    for (x1, y1), (x2, y2) in zip(two_steps, at_once):
        assert abs(x1 - x2) <= TOL and abs(y1 - y2) <= TOL


def test_composition_matches_numpy_matmul_in_both_orders():
    A = np.array(shear_x(2.0))
    Bm = np.array(rotation(math.pi / 2))
    assert np.allclose(Bm @ A, np.array(compose(rotation(math.pi / 2), shear_x(2.0))), atol=TOL)
    assert np.allclose(A @ Bm, np.array(compose(shear_x(2.0), rotation(math.pi / 2))), atol=TOL)
    assert not np.allclose(A @ Bm, Bm @ A, atol=TOL)


def test_composing_with_the_identity_changes_nothing():
    H = shear_x(2.0)
    assert compose(identity(), H) == H
    assert compose(H, identity()) == H


def test_determinants_multiply_under_composition():
    A, Bm = scaling(2.0, 3.0), shear_x(4.0)
    assert abs(determinant(compose(Bm, A)) - determinant(A) * determinant(Bm)) <= TOL


# -- Determinant, area, orientation ------------------------------------------


def test_the_unit_square_starts_with_signed_area_one():
    assert abs(signed_area(shapes.UNIT_SQUARE) - 1.0) <= TOL


def test_listing_the_corners_clockwise_flips_the_sign():
    assert abs(signed_area(list(reversed(shapes.UNIT_SQUARE))) + 1.0) <= TOL


@pytest.mark.parametrize(
    "matrix, expected",
    [
        (scaling(2.0, 3.0), 6.0),
        (scaling(0.5, 0.5), 0.25),
        (shear_x(2.0), 1.0),
        (reflection_in_x_axis(), -1.0),
        (reflection_in_y_axis(), -1.0),
        (identity(), 1.0),
        (shapes.PICTURE_MATRIX, 7.0),
        (shapes.COLLAPSE_MATRIX, 0.0),
    ],
)
def test_transformed_area_equals_the_determinant(matrix, expected):
    area = signed_area(transform_polygon(matrix, shapes.UNIT_SQUARE))
    assert abs(area - expected) <= TOL
    assert abs(determinant(matrix) - expected) <= TOL


def test_a_negative_determinant_means_the_orientation_flipped():
    flipped = transform_polygon(reflection_in_x_axis(), shapes.UNIT_SQUARE)
    assert signed_area(flipped) < 0
    assert abs(abs(signed_area(flipped)) - 1.0) <= TOL


def test_a_rotation_never_flips_orientation():
    for degrees in (0, 30, 90, 180, 270, 359):
        assert determinant(rotation(math.radians(degrees))) > 0


def test_rotation_determinant_is_one_which_is_the_pythagorean_identity():
    for degrees in (17, 45, 123, 300):
        theta = math.radians(degrees)
        assert abs(determinant(rotation(theta)) - 1.0) <= TOL
        # det = cos*cos - (-sin)*sin = cos^2 + sin^2, which is 1 by Pythagoras
        assert abs(math.cos(theta) ** 2 + math.sin(theta) ** 2 - 1.0) <= TOL


def test_the_from_scratch_determinant_is_exact_on_whole_numbers():
    assert determinant(shapes.PICTURE_MATRIX) == 7.0
    assert determinant(shapes.COLLAPSE_MATRIX) == 0.0


def test_numpy_determinant_is_close_but_not_always_equal():
    """The honest difference, asserted rather than described.

    numpy.linalg.det factorises the matrix -- the general method that also
    works at 500 by 500 -- and that rounds. The direct a*d - b*c does not. On
    this matrix, on this machine, they differ in the last bit.
    """
    mine = determinant(shapes.PICTURE_MATRIX)
    theirs = float(np.linalg.det(np.array(shapes.PICTURE_MATRIX)))
    assert mine == 7.0
    assert theirs != 7.0
    assert abs(theirs - 7.0) < 1e-14


# -- Rank and collapse --------------------------------------------------------


def test_the_collapse_puts_everything_on_one_line():
    G = shapes.COLLAPSE_MATRIX
    for point in [(1.0, 0.0), (0.0, 1.0), (3.0, -1.0), (7.0, 7.0), (-2.5, 0.25)]:
        x, y = apply(G, point)
        assert abs(y - 2.0 * x) <= TOL


def test_the_collapse_sends_two_different_points_to_the_same_place():
    G = shapes.COLLAPSE_MATRIX
    assert apply(G, (2.0, 0.0)) == apply(G, (0.0, 1.0))


@pytest.mark.parametrize(
    "matrix, expected_rank",
    [
        (identity(), 2),
        (scaling(2.0, 3.0), 2),
        (shear_x(2.0), 2),
        (shapes.PICTURE_MATRIX, 2),
        (shapes.COLLAPSE_MATRIX, 1),
        ([[0.0, 0.0], [0.0, 0.0]], 0),
        ([[1.0, 0.0], [0.0, 0.0]], 1),
    ],
)
def test_rank_matches_numpy(matrix, expected_rank):
    assert rank(matrix) == expected_rank
    assert int(np.linalg.matrix_rank(np.array(matrix))) == expected_rank


def test_scaling_by_zero_in_one_direction_loses_a_dimension():
    flat = scaling(3.0, 0.0)
    assert determinant(flat) == 0.0
    assert rank(flat) == 1
    for point in [(1.0, 5.0), (-2.0, 100.0)]:
        assert apply(flat, point)[1] == 0.0


# -- Inverses ------------------------------------------------------------------


def test_inverse_of_a_shear_is_the_opposite_shear():
    assert inverse(shear_x(2.0)) == shear_x(-2.0)


def test_inverse_of_a_scaling_divides():
    inv = inverse(scaling(2.0, 4.0))
    assert inv == [[0.5, 0.0], [0.0, 0.25]]


def test_inverse_composed_with_the_original_is_the_identity():
    for M in [shear_x(2.0), scaling(2.0, 3.0), shapes.PICTURE_MATRIX, rotation(0.7)]:
        back = compose(inverse(M), M)
        for row, want in zip(back, identity()):
            for got, expect in zip(row, want):
                assert abs(got - expect) <= TOL


def test_a_round_trip_returns_the_original_point():
    M = shapes.PICTURE_MATRIX
    for point in [(1.0, 1.0), (-3.0, 2.5), (0.0, 0.0)]:
        there = apply(M, point)
        home = apply(inverse(M), there)
        assert abs(home[0] - point[0]) <= TOL
        assert abs(home[1] - point[1]) <= TOL


def test_the_inverse_determinant_is_the_reciprocal():
    M = scaling(2.0, 3.0)
    assert abs(determinant(inverse(M)) - 1.0 / determinant(M)) <= TOL


def test_inverting_a_collapse_raises():
    with pytest.raises(SingularMatrix):
        inverse(shapes.COLLAPSE_MATRIX)


def test_the_from_scratch_refusal_is_catchable_as_a_ValueError():
    with pytest.raises(ValueError):
        inverse(shapes.COLLAPSE_MATRIX)


def test_numpy_raises_LinAlgError_for_the_same_matrix():
    with pytest.raises(np.linalg.LinAlgError) as caught:
        np.linalg.inv(np.array(shapes.COLLAPSE_MATRIX))
    assert "Singular matrix" in str(caught.value)


def test_numpys_error_is_also_a_ValueError():
    assert issubclass(np.linalg.LinAlgError, ValueError)
    with pytest.raises(ValueError):
        np.linalg.inv(np.array(shapes.COLLAPSE_MATRIX))


def test_numpy_inverse_agrees_with_the_from_scratch_one():
    for M in [shear_x(2.0), scaling(2.0, 3.0), shapes.PICTURE_MATRIX]:
        assert np.allclose(np.linalg.inv(np.array(M)), np.array(inverse(M)), atol=TOL)


# -- The limit of linear -------------------------------------------------------


def test_a_stack_of_twenty_layers_is_one_matrix():
    import random

    random.seed(102)
    stack = [
        [[random.uniform(-2, 2) for _ in range(2)] for _ in range(2)]
        for _ in range(20)
    ]
    point = (0.7, -0.4)
    stepwise = point
    for layer in stack:
        stepwise = apply(layer, stepwise)

    combined = stack[0]
    for layer in stack[1:]:
        combined = compose(layer, combined)
    at_once = apply(combined, point)

    # A relative tolerance, because twenty layers of rounding on entries that
    # are not small whole numbers is a real accumulation and pretending
    # otherwise would be dishonest.
    for a, b in zip(stepwise, at_once):
        assert abs(a - b) <= 1e-9 * max(1.0, abs(a))


def test_a_stack_still_fixes_the_origin_and_still_makes_a_parallelogram():
    M = compose(compose(rotation(0.4), shear_x(3.0)), scaling(2.0, -1.0))
    corners = transform_polygon(M, shapes.UNIT_SQUARE)
    assert abs(corners[0][0]) <= TOL and abs(corners[0][1]) <= TOL
    # Opposite sides of the image are still parallel and equal, which is what
    # "still a parallelogram" means numerically.
    side_a = (corners[1][0] - corners[0][0], corners[1][1] - corners[0][1])
    side_c = (corners[2][0] - corners[3][0], corners[2][1] - corners[3][1])
    assert abs(side_a[0] - side_c[0]) <= TOL
    assert abs(side_a[1] - side_c[1]) <= TOL


def test_relu_after_a_matrix_is_not_linear():
    M = compose(rotation(math.radians(30)), shear_x(1.5))

    def relu_layer(point):
        return tuple(max(0.0, c) for c in apply(M, point))

    assert not is_linear(relu_layer, (1.0, -1.0), (0.5, 2.0), 3.0, TOL)


def test_the_area_factor_of_a_stack_is_the_product_of_its_factors():
    layers = [scaling(2.0, 3.0), shear_x(4.0), reflection_in_x_axis()]
    combined = layers[0]
    for layer in layers[1:]:
        combined = compose(layer, combined)
    expected = 1.0
    for layer in layers:
        expected *= determinant(layer)
    assert abs(determinant(combined) - expected) <= TOL
    assert abs(expected + 6.0) <= TOL  # 6 * 1 * -1, and the sign says it flipped

"""Linear transformations of the plane, built from nothing but arithmetic.

The reference implementation. Everything here works on plain Python lists and
tuples, uses only `math` from the standard library, and never imports NumPy --
so you can read every line and see exactly what a transformation is doing. The
scripts beside this file then check all of it against NumPy, which is the only
way to know that "I wrote it myself" and "it is right" are both true.

The one idea the whole module is built on:

    A matrix IS a function. Its columns are where the basis vectors land.

Column 0 is where (1, 0) goes. Column 1 is where (0, 1) goes. Everything else
follows, because every vector (x, y) is x * (1, 0) + y * (0, 1), and a linear
transformation is exactly one that keeps that combination intact.
"""

from __future__ import annotations

import math

Point = tuple[float, float]
Matrix = list[list[float]]


class SingularMatrix(ValueError):
    """Raised when a matrix has no inverse.

    Subclasses ValueError deliberately, to match NumPy: `numpy.linalg.inv` on a
    matrix with no inverse raises `numpy.linalg.LinAlgError`, and that class is
    itself a subclass of ValueError. So `except ValueError` catches both, and
    code written against one behaves the same against the other.
    """


# -- Reading a matrix, and reading it back -------------------------------------


def columns_of(matrix: Matrix) -> tuple[Point, Point]:
    """Return the two columns as points: (where e1 lands, where e2 lands).

    A 2 by 2 matrix is written as rows -- [[a, b], [c, d]] -- but it MEANS its
    columns. Column 0 is (a, c) and column 1 is (b, d). Half of all confusion
    about transformation matrices is reading (a, b) as a landing place when it
    is a row.
    """
    (a, b), (c, d) = matrix
    return (a, c), (b, d)


def from_landings(e1_lands_at: Point, e2_lands_at: Point) -> Matrix:
    """Build the matrix that sends (1, 0) and (0, 1) to the two given points.

    This is the whole day in four lines. If you can see where the basis vectors
    land, you can write the matrix down, and the matrix then tells you where
    every other vector lands without you having to look at the picture again.
    """
    (a, c), (b, d) = e1_lands_at, e2_lands_at
    return [[a, b], [c, d]]


# -- Applying and composing ----------------------------------------------------


def apply(matrix: Matrix, point: Point) -> Point:
    """Send one point through the transformation.

    Written the way the idea is stated rather than the way a textbook writes
    it: the answer is x lots of the first column plus y lots of the second.
    That is the same arithmetic as the usual row-by-row rule, but it says out
    loud why the columns are the landing places.
    """
    x, y = point
    (e1x, e1y), (e2x, e2y) = columns_of(matrix)
    return (x * e1x + y * e2x, x * e1y + y * e2y)


def compose(second: Matrix, first: Matrix) -> Matrix:
    """Return the single matrix that does `first` and then `second`.

    Note the argument order, and note that it is not an accident. Written out,
    applying `first` and then `second` to a vector v is

        second @ (first @ v)

    and the matrix that does both in one step is `second @ first` -- the one
    that happens FIRST is written on the RIGHT, because it is the one standing
    next to the vector. Reading a product right to left is not a quirk to
    memorise; it is what the notation means.

    The product is built one column at a time, which is again the day's idea:
    column 0 of the answer is wherever (1, 0) ends up after both steps.
    """
    e1, e2 = columns_of(first)
    return from_landings(apply(second, e1), apply(second, e2))


def identity() -> Matrix:
    """The do-nothing transformation: e1 stays at (1, 0), e2 stays at (0, 1)."""
    return [[1.0, 0.0], [0.0, 1.0]]


# -- The four standard transformations, each DERIVED from its landings ---------


def scaling(sx: float, sy: float) -> Matrix:
    """Stretch by sx horizontally and sy vertically.

    Derivation: (1, 0) is one step right, and stretching horizontally by sx
    makes it sx steps right, so it lands at (sx, 0). (0, 1) is one step up and
    lands at (0, sy). Write those two down as columns and you have the matrix.
    """
    return from_landings((sx, 0.0), (0.0, sy))


def reflection_in_x_axis() -> Matrix:
    """Mirror the plane in the horizontal axis: up becomes down.

    Derivation: (1, 0) already lies on the mirror line, so it does not move.
    (0, 1) is one step up, and its mirror image is one step down, at (0, -1).
    """
    return from_landings((1.0, 0.0), (0.0, -1.0))


def reflection_in_y_axis() -> Matrix:
    """Mirror the plane in the vertical axis: right becomes left."""
    return from_landings((-1.0, 0.0), (0.0, 1.0))


def shear_x(k: float) -> Matrix:
    """Push the plane sideways by k times its height.

    Derivation: a point's sideways push is proportional to how high it is.
    (1, 0) has height 0, so it is pushed by nothing and does not move. (0, 1)
    has height 1, so it is pushed k to the right and lands at (k, 1).

    This is the transformation to picture when someone says a deck of cards
    slid sideways: the bottom card stays put, every card above it slides
    further, and the deck's volume never changes.
    """
    return from_landings((1.0, 0.0), (k, 1.0))


def shear_y(k: float) -> Matrix:
    """Push the plane upwards by k times its horizontal distance."""
    return from_landings((1.0, k), (0.0, 1.0))


def rotation(theta: float) -> Matrix:
    """Turn the whole plane anticlockwise by theta RADIANS about the origin.

    Derivation, which needs the unit circle and nothing else. Draw a circle of
    radius 1 around the origin. Start at (1, 0) and walk anticlockwise around
    the rim until you have turned through an angle theta. The two numbers that
    name where you now stand are, by definition, the cosine and the sine of
    theta: cos(theta) across, sin(theta) up. That is what those two functions
    ARE -- the coordinates of a point on the unit circle -- and every identity
    about them is a fact about that picture.

    So (1, 0) lands at (cos(theta), sin(theta)).

    (0, 1) is (1, 0) already turned a quarter of a turn anticlockwise, and
    turning it a further theta puts it a quarter turn ahead of the first
    landing place. A quarter turn anticlockwise sends any point (x, y) to
    (-y, x) -- push it round and the across-ness becomes up-ness. Applying
    that to (cos(theta), sin(theta)) gives (-sin(theta), cos(theta)).

    Write the two landings as columns:

        [[cos(theta), -sin(theta)],
         [sin(theta),  cos(theta)]]

    which is the rotation matrix, derived rather than remembered.

    Radians, briefly: an angle measured as the distance you walked around the
    rim of that unit circle. A full turn is the whole circumference, 2 * pi. So
    a quarter turn is pi / 2, and `math.radians(90)` converts if you would
    rather think in degrees.
    """
    cos_t, sin_t = math.cos(theta), math.sin(theta)
    return from_landings((cos_t, sin_t), (-sin_t, cos_t))


# -- Determinant, inverse, rank ------------------------------------------------


def determinant(matrix: Matrix) -> float:
    """The factor by which the transformation multiplies area, sign included.

    For a 2 by 2 matrix [[a, b], [c, d]] the answer is a*d - b*c, and this
    function computes that exactly -- one multiply, one multiply, one subtract.
    No rearrangement, so no rounding beyond the inputs themselves.

    What the number means:

      * its SIZE is the area factor. A unit square of area 1 comes out with
        area |determinant|.
      * its SIGN is the orientation. Positive means the plane was not flipped
        over; negative means it was, and a shape listed anticlockwise comes out
        listed clockwise.
      * ZERO means the plane was flattened onto a line or onto a point. Area
        1 became area 0, information was destroyed, and nothing can undo it.
    """
    (a, b), (c, d) = matrix
    return a * d - b * c


def inverse(matrix: Matrix) -> Matrix:
    """The transformation that undoes this one.

    Derived by asking the only question that matters: which matrix, composed
    with this one, leaves everything where it started? For 2 by 2 the answer is
    the standard formula, one over the determinant times [[d, -b], [-c, a]].

    It exists precisely when the determinant is not zero, which is the same
    sentence as "precisely when no area was destroyed". If two different
    starting points landed on the same place, no rule could send that place
    back to both of them, so there is nothing to return.
    """
    det = determinant(matrix)
    if det == 0.0:
        raise SingularMatrix(
            "Singular matrix: the determinant is 0, so this transformation "
            "collapses the plane and cannot be undone"
        )
    (a, b), (c, d) = matrix
    return [[d / det, -b / det], [-c / det, a / det]]


def rank(matrix: Matrix, tol: float = 1e-12) -> int:
    """How many dimensions survive the transformation.

    Plain language: feed the whole plane in, and look at what comes out. If the
    output fills the plane, the rank is 2. If it is squashed onto a line, the
    rank is 1. If everything lands on the origin, the rank is 0.

    For a 2 by 2 matrix that reads straight off the columns. If the determinant
    is not zero the two columns point in genuinely different directions and
    between them they reach everywhere, so the rank is 2. If the determinant is
    zero but at least one column is not the zero vector, everything lands on
    the line through that column, so the rank is 1. If both columns are zero,
    everything lands on the origin and the rank is 0.

    `tol` is here because the determinant of a matrix built from cosines will
    rarely be exactly 0.0, and asking `== 0` of a computed float is how you get
    a rank of 2 for a matrix that has plainly collapsed.
    """
    if abs(determinant(matrix)) > tol:
        return 2
    if any(abs(entry) > tol for row in matrix for entry in row):
        return 1
    return 0


# -- Polygons ------------------------------------------------------------------


def transform_polygon(matrix: Matrix, polygon: list[Point]) -> list[Point]:
    """Send every corner of a polygon through the transformation.

    Corners are enough. A linear transformation sends straight lines to
    straight lines, so the edges take care of themselves -- which is exactly
    the property that makes these transformations cheap, and exactly the
    property that stops a stack of them ever drawing a curve.
    """
    return [apply(matrix, point) for point in polygon]


def signed_area(polygon: list[Point]) -> float:
    """The area of a polygon, negative if its corners run clockwise.

    The shoelace formula: walk the corners in order, and for each edge add
    x_here * y_next - x_next * y_here. Halve the total. It is called the
    shoelace formula because the cross-multiplied pairs criss-cross like the
    lacing on a shoe.

    The sign is the part this lab uses. List the unit square anticlockwise and
    the signed area is +1; transform it by a matrix with a negative determinant
    and the signed area comes out negative, because the corners now run the
    other way round. That is what "the plane was flipped over" means when it is
    measured rather than described.
    """
    total = 0.0
    count = len(polygon)
    for i in range(count):
        x_here, y_here = polygon[i]
        x_next, y_next = polygon[(i + 1) % count]
        total += x_here * y_next - x_next * y_here
    return total / 2.0


# -- Linearity -----------------------------------------------------------------


def preserves_addition(
    func, u: Point, v: Point, tol: float = 1e-12
) -> tuple[bool, Point, Point]:
    """Test whether func(u + v) equals func(u) + func(v).

    The first half of the definition of linear. Returns the verdict and both
    sides, because when a function fails this test the interesting part is not
    that it failed -- it is by how much, and whether the gap is the same every
    time.
    """
    together = func((u[0] + v[0], u[1] + v[1]))
    separately = tuple(a + b for a, b in zip(func(u), func(v)))
    ok = all(abs(a - b) <= tol for a, b in zip(together, separately))
    return ok, together, separately  # type: ignore[return-value]


def preserves_scaling(
    func, u: Point, s: float, tol: float = 1e-12
) -> tuple[bool, Point, Point]:
    """Test whether func(s * u) equals s * func(u).

    The second half of the definition. A function needs BOTH halves to be
    linear, and a function that fails either one cannot be written as a matrix,
    no matter how simple it looks.
    """
    scaled_first = func((s * u[0], s * u[1]))
    scaled_after = tuple(s * component for component in func(u))
    ok = all(abs(a - b) <= tol for a, b in zip(scaled_first, scaled_after))
    return ok, scaled_first, scaled_after  # type: ignore[return-value]


def is_linear(func, u: Point, v: Point, s: float, tol: float = 1e-12) -> bool:
    """Both halves at once, on one chosen pair of vectors and one scalar.

    An honest caveat, and it matters: passing this on one example does not
    prove a function is linear. The definition quantifies over EVERY pair of
    vectors and every scalar, and no finite number of examples can settle that.
    What a test like this does well is the other direction -- a single failure
    is a complete disproof, and that is how it is used here.
    """
    return (
        preserves_addition(func, u, v, tol)[0]
        and preserves_scaling(func, u, s, tol)[0]
    )

lesson_id: D102
day: 102
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-102-linear-transformations
  - 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_columns_are_landings.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_building_the_transformations.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_linear_or_not.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_composition_and_order.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_determinant_inverse_rank.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_the_limit_of_linear.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: 30
last_executed: '2026-08-16'
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 -> 64 checks, 0 failure(s), exit 0; pytest examples -> 80 passed; pytest starter -> 1 passed, 53 skipped on an untouched checkout, and 54 passed against a fully solved copy of starter/ kept outside the lab. All six reference scripts exit 0 with every internal assertion holding. Network is needed once to install numpy and pytest; nothing else in the lab opens a socket, and section 7 of the harness greps the sources to prove it. Section 6 re-runs the harness with one expectation deliberately swapped for the naive belief that cos(pi/2) is exactly 0.0, 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. Two numbers in the captured output are real floating-point artefacts and are asserted rather than hidden: cos(pi/2) is 6.123233995736766e-17 rather than 0.0, and numpy.linalg.det on [[3, -1], [1, 2]] returns 7.000000000000001 where the direct a*d-b*c formula returns exactly 7.0.'

# Dependencies for the Day 102 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` | The library your from-scratch transformations are checked against: the `@` operator for matrix products, `numpy.linalg.det`, `numpy.linalg.inv`, `numpy.linalg.matrix_rank`, and the `LinAlgError` a singular matrix raises. |
| `pytest` | `9.1.1` | The test runner from Days 071–074. Nothing new here except what it is pointed at. |

Nothing else is required. Both suites and all six reference scripts run on
those two packages plus the standard library — and the from-scratch module
itself imports nothing but `math`.

## Why the from-scratch code deliberately does not use NumPy

`examples/transforms.py` and `starter/transforms.py` are written with plain
lists and tuples on purpose. If your `rotation` function returned a NumPy array
built by a NumPy helper, then checking it against NumPy would be checking NumPy
against itself. Writing the arithmetic yourself and *then* having a mature
library agree with you is worth more than either half alone.

## Why numpy is pinned

The version is *checked* rather than assumed. Section 1 of `tests/run_tests.sh`
reads the installed version and compares it against `requirements.txt`, so a
mismatch is reported at the top of the run instead of surfacing later as a
confusing diff.

There is one place where the version could matter. `numpy.linalg.det` on the
matrix `[[3, -1], [1, 2]]` returns `7.000000000000001` on this machine with
this version — one bit away from the exact 7 the direct formula gives. That is
a property of the factorisation it uses and of the underlying linear algebra
library, and it is not guaranteed to be identical everywhere. The lab handles
this honestly: it asserts that NumPy's answer is *within 1e-14 of 7* and that
the from-scratch answer is *exactly 7*, rather than pinning NumPy's last digit.

The version was read from the installed package rather than guessed:

```bash
.venv/bin/python3 -c "from importlib.metadata import version; print(version('numpy'))"
```

On the authoring machine, on 16 August 2026, that printed `2.5.2`.

## 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.** No
script opens a socket, reads a URL or contacts a service, and section 7 of
`tests/run_tests.sh` greps every file under `examples/` and `starter/` for the
patterns that would indicate otherwise.

## 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 NumPy is 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 without NumPy

More than you might expect, and less than you might fear. Exercise 1 — writing
all ten transformation functions — needs nothing but `math` from the standard
library, and you can complete the whole of it on a bare `python3`. What you
lose is every cross-check: the tests that confirm your `compose` matches `@`,
your `inverse` matches `numpy.linalg.inv`, and your `rank` matches
`numpy.linalg.matrix_rank`, plus the exercise on the exception a singular
matrix raises. Those parts cannot be faked and the lab does not pretend
otherwise.

numpy==2.5.2
pytest==9.1.1

# Day 102 lab — Where Do the Basis Vectors Land?

One idea holds this whole lab together:

> **A matrix IS a function, and its columns are where the basis vectors land.**

If you know where `(1, 0)` and `(0, 1)` go, you know where everything goes —
because every vector `(x, y)` is `x * (1, 0) + y * (0, 1)`, and a linear
transformation is precisely one that keeps that combination intact.

Everything below is a consequence of that sentence. Work in order; each
exercise uses the one before it.

Check yourself at any point, from the **lab directory** (the one above this
file):

```bash
.venv/bin/pytest starter -q
```

Unattempted work is **skipped**, not failed. On an untouched checkout you will
see `1 passed, 53 skipped`. When it says `54 passed`, you are finished.

---

## Exercise 1 — `transforms.py` (ten functions)

Write the ten functions marked `raise NotImplementedError`. Each docstring
gives the derivation and a worked example you can check on paper. Use only
`math` from the standard library in this file — NumPy appears in the tests,
where it checks your work, which is the right way round.

| Step | Function | What it must do |
| --- | --- | --- |
| 1.1 | `from_landings` | Two landing places in, one matrix out. The landings are the **columns**. |
| 1.2 | `columns_of` | The exact inverse of 1.1. |
| 1.3 | `apply` | `x` lots of the first column plus `y` lots of the second. |
| 1.4 | `scaling` | Derive it: where does one step right go when the plane stretches? |
| 1.5 | `reflection_in_x_axis` | Derive it: what happens to a point that is already **on** the mirror line? |
| 1.6 | `shear_x` | Derive it: the sideways push is proportional to the **height**. |
| 1.7 | `rotation` | Derive it from the unit circle. The docstring defines cosine and sine from scratch if you have not met them. |
| 1.8 | `compose` | Column 0 is where `(1, 0)` ends up after **both** steps. Mind the argument order. |
| 1.9 | `determinant` | `a*d - b*c`, computed directly so whole numbers stay exact. |
| 1.10 | `inverse` | The undo. Raise `SingularMatrix` when the determinant is 0. |

Three helpers are written for you at the bottom of the file — `identity`,
`transform_polygon`, `signed_area` and `rank`. Read them; the tests use them.

**The gotcha in 1.8.** `compose(second, first)` returns the matrix that does
`first` and then `second`, so the step that happens **first** is written on the
**right**. That is not an arbitrary convention: applying `first` and then
`second` to a vector `v` is `second @ (first @ v)`, and `first` is the one
standing next to the vector.

---

## Exercise 2 — read a matrix off a picture (`answers.py`)

There is no picture file. The drawing is described in words in `shapes.py`,
because the skill being trained is going from *where did the basis vectors
land* to *what is the matrix*, and looking at a picture would let you skip the
step.

The drawing shows two things and only two things:

```
the arrow (1, 0) redrawn ending at ( 3, 1)
the arrow (0, 1) redrawn ending at (-1, 2)
```

Write the matrix. Then say where `(2, 1)` lands, without drawing anything.

Exercise 2.3 asks which **row** of that matrix is a landing place. Read the
question carefully before answering; it is the single most common mistake with
transformation matrices.

---

## Exercise 3 — the four standard transformations

Predict where specific points land under scaling, reflection, shear and a
quarter turn. Each answer is one line of arithmetic.

Exercise 3.6 is the one worth slowing down for. On paper, a quarter turn sends
`(1, 0)` to exactly `(0, 1)`. In binary floating point, `cos(pi / 2)` comes out
as `6.123233995736766e-17`, because `pi` cannot be stored exactly and the
cosine of the stored value is not the cosine of pi. **That is why every float
comparison in this lab uses a tolerance of `1e-12` and none of them use `==`.**
The number is not arbitrary: it sits about five orders of magnitude above that
rounding error and about four below the smallest quantity the lab cares about.

Exercise 3.7 asks for the one point no matrix anywhere can move. Once you see
why, half of exercise 4 answers itself.

---

## Exercise 4 — linear, and not linear

*Linear* means exactly two things:

```
T(u + v) = T(u) + T(v)        it preserves addition
T(s * u) = s * T(u)           it preserves scalar multiplication
```

Both, for every `u`, `v` and `s`. Nothing else.

With `M = [[2, 0], [0, 3]]`, `b = (1, 1)`, `u = (1, 2)`, `v = (3, -1)`, `s = 5`,
you check both properties for `T(v) = M @ v` and for `f(v) = M @ v + b`.

`f` fails both. Compute the gap in 4.5 before you read any further explanation
— the size of the gap is a recognisable quantity, and recognising it is the
exercise.

This is not a curiosity. It is the reason a neural network layer is written
`X @ W + b` with the bias kept separate rather than folded into the matrix, and
the reason the word *affine* exists.

---

## Exercise 5 — composition and order

`A = shear_x(2)`, `B = rotation(pi / 2)`. You shear first, then rotate.

Which product is that? Write it out. Then decide whether the other order gives
the same matrix. It does not, and 5.3 wants you to have checked rather than
assumed.

5.4 is a freebie with a real consequence behind it: determinants multiply under
composition, so if any step in a pipeline has determinant 0, the whole pipeline
does, and nothing downstream can recover what that step destroyed.

---

## Exercise 6 — determinant, area, orientation, rank and the inverse

The determinant is introduced here the way it is actually useful: send the unit
square through the transformation and measure the area of what comes out. That
number, **with its sign**, is the determinant.

- 6.1–6.4: positive determinants, and one negative one. What does the sign say?
- 6.5–6.7: `[[1, 2], [2, 4]]`. Look at its two columns before predicting
  anything. Where does everything end up?
- 6.8: name the exception class `numpy.linalg.inv` raises on it. Give the class
  itself, not a string. `numpy` is already imported at the top of `answers.py`.
- 6.9–6.10: two inverses you can write down without the formula, if you think
  about what would undo each transformation.

---

## When you are done

Read the reference. Each script prints its working and asserts every claim it
makes, so nothing in it is decoration:

```bash
cd examples
../.venv/bin/python3 01_columns_are_landings.py
../.venv/bin/python3 02_building_the_transformations.py
../.venv/bin/python3 03_linear_or_not.py
../.venv/bin/python3 04_composition_and_order.py
../.venv/bin/python3 05_determinant_inverse_rank.py
../.venv/bin/python3 06_the_limit_of_linear.py
cd ..
```

Script 06 is the payoff, and it is worth reading even if you stop everything
else. It shows that a linear transformation always fixes the origin, always
sends straight lines to straight lines, and that a stack of twenty of them
collapses into a single 2 by 2 matrix. Depth bought nothing. That is the
concrete, measured reason a non-linear activation function sits between the
layers of a network — not a rule to memorise, a limitation you can watch
happening.

"""Exercises 2 to 6 -- your predictions. Work them out BEFORE running anything.

Every one of these can be done on paper in under a minute. That is the point:
a lab about transformations whose answers you cannot check by hand is a lab
that teaches you to trust output.

Replace each `None` with your answer. Anything still `None` is SKIPPED by the
test suite rather than failed, so your score only ever counts work you actually
attempted.

Check yourself from the LAB DIRECTORY:

    .venv/bin/pytest starter -q
"""

# Imported for you: exercise 6.8 asks for an exception CLASS, and the one it is
# looking for lives at numpy.linalg.LinAlgError.
import numpy

# =============================================================================
# Exercise 2 -- reading a matrix off a picture
# =============================================================================
#
# The picture shows the arrow (1, 0) redrawn ending at (3, 1), and the arrow
# (0, 1) redrawn ending at (-1, 2). Nothing else.

# 2.1 Write the matrix, as a list of two ROWS. Careful: the landing places are
#     the COLUMNS, so they are read downwards, not across.
#     Example of the format: [[1.0, 0.0], [0.0, 1.0]]
PICTURE_MATRIX = None

# 2.2 Where does (2, 1) land? Work it out as 2 lots of the first landing place
#     plus 1 lot of the second, and give an (x, y) tuple of floats.
PICTURE_SENDS_2_1_TO = None

# 2.3 Which of the two ROWS of that matrix is a landing place for a basis
#     vector -- 0, 1, or neither? Answer with the integer 0, the integer 1, or
#     the string "neither".
WHICH_ROW_IS_A_LANDING = None


# =============================================================================
# Exercise 3 -- the four standard transformations
# =============================================================================

# 3.1 scaling(2, 3) applied to (1, 1). An (x, y) tuple.
SCALE_SENDS_1_1_TO = None

# 3.2 reflection_in_x_axis() applied to (2, 3). An (x, y) tuple.
FLIP_SENDS_2_3_TO = None

# 3.3 shear_x(2) applied to (1, 1). An (x, y) tuple.
#     Remember: the sideways push is proportional to the HEIGHT.
SHEAR_SENDS_1_1_TO = None

# 3.4 shear_x(2) applied to (5, 0). An (x, y) tuple. Think before you compute.
SHEAR_SENDS_5_0_TO = None

# 3.5 A quarter turn anticlockwise, rotation(pi / 2), applied to (1, 0).
#     Give the answer you would write on paper, as an (x, y) tuple. The test
#     compares with a tolerance of 1e-12 rather than with ==, and exercise 3.6
#     is about why.
QUARTER_TURN_SENDS_1_0_TO = None

# 3.6 In binary floating point, is math.cos(math.pi / 2) exactly 0.0?
#     True or False.
COS_OF_QUARTER_TURN_IS_EXACTLY_ZERO = None

# 3.7 Every matrix in this lab sends one particular point to itself, no matter
#     what the four entries are. Which point? An (x, y) tuple.
THE_POINT_NO_MATRIX_CAN_MOVE = None


# =============================================================================
# Exercise 4 -- linearity, and the function that fails it
# =============================================================================
#
# T(v) = M @ v          with M = [[2, 0], [0, 3]]
# f(v) = M @ v + b      with b = (1, 1)
# u = (1, 2), v = (3, -1), s = 5

# 4.1 T(u + v). An (x, y) tuple.
T_OF_U_PLUS_V = None

# 4.2 T(u) + T(v). An (x, y) tuple.
T_OF_U_PLUS_T_OF_V = None

# 4.3 f(u + v). An (x, y) tuple.
F_OF_U_PLUS_V = None

# 4.4 f(u) + f(v). An (x, y) tuple.
F_OF_U_PLUS_F_OF_V = None

# 4.5 Subtract 4.3 from 4.4. The gap is one recognisable quantity -- which?
#     An (x, y) tuple.
THE_GAP_BETWEEN_THEM = None

# 4.6 Is f linear? True or False.
F_IS_LINEAR = None

# 4.7 f((0, 0)). An (x, y) tuple -- and notice what it tells you about 4.6
#     without needing 4.1 to 4.5 at all.
F_OF_THE_ORIGIN = None


# =============================================================================
# Exercise 5 -- composition and order
# =============================================================================
#
# A = shear_x(2), B = rotation(pi / 2). You shear FIRST and then rotate.

# 5.1 Which expression is the single matrix that does shear-then-rotate?
#     Answer with the string "compose(B, A)" or the string "compose(A, B)".
SHEAR_THEN_ROTATE_IS = None

# 5.2 Write that matrix out, as two rows of floats. Work it out from where the
#     two basis vectors end up after both steps -- it is easier than the
#     row-by-column rule and it is the same answer.
SHEAR_THEN_ROTATE_MATRIX = None

# 5.3 Do the two orders give the same matrix? True or False.
BOTH_ORDERS_AGREE = None

# 5.4 det(A) is 1 and det(B) is 1. What is the determinant of the composite?
#     A float.
DET_OF_THE_COMPOSITE = None


# =============================================================================
# Exercise 6 -- determinant, area, orientation, rank and the inverse
# =============================================================================

# 6.1 The unit square has area 1. What area does scaling(2, 3) give it?
#     A float.
AREA_AFTER_SCALING = None

# 6.2 What is the SIGNED area of the unit square after reflection in the x
#     axis? A float. Mind the sign; that is the whole question.
SIGNED_AREA_AFTER_REFLECTION = None

# 6.3 What does a negative determinant tell you? Answer with one of the
#     strings: "the shape got smaller", "the plane was flipped over",
#     "the transformation cannot be undone".
A_NEGATIVE_DETERMINANT_MEANS = None

# 6.4 What is the determinant of shear_x(2)? A float, and it should surprise
#     you slightly less once you have drawn the sheared square.
DET_OF_SHEAR = None

# 6.5 COLLAPSE_MATRIX is [[1, 2], [2, 4]]. Look at its two columns. What is its
#     determinant? A float.
DET_OF_COLLAPSE = None

# 6.6 What is its rank -- how many dimensions survive? The integer 0, 1 or 2.
RANK_OF_COLLAPSE = None

# 6.7 Every vector it touches lands on one line through the origin. The line is
#     y = m * x. What is m? A float.
COLLAPSE_LANDS_EVERYTHING_ON_THE_LINE_Y_EQUALS = None

# 6.8 numpy.linalg.inv of that matrix raises an exception. Which class?
#     Give the class itself, not a string, for example: ValueError
#     `numpy` is already imported for you below.
COLLAPSE_INVERSE_EXCEPTION = None

# 6.9 The inverse of shear_x(2) is another shear. With what k? A float.
INVERSE_OF_SHEAR_IS_SHEAR_WITH_K = None

# 6.10 The inverse of scaling(2, 4), as two rows of floats.
INVERSE_OF_SCALING_2_4 = None

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

Both `examples/` and `starter/` contain modules called `transforms` and
`shapes`, and pytest imports test files by putting their directory on
`sys.path`. Without this file, running `pytest` across both directories at once
would import whichever `transforms` 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
`transforms` or `shapes` that came from somewhere else.
"""

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 ("transforms", "shapes", "answers"):
    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 this lab works on. Read it; you do not need to change it.

Everything here is invented. A 2 by 2 matrix has four numbers in it, and the
whole point of the day is that you can check every one of them on paper.

Points are `(x, y)` pairs. Polygons are lists of corners in counter-clockwise
order, because the SIGN of a polygon's area depends on which way round its
corners are listed, and that sign is what tells you whether a transformation
turned the plane over.
"""

# -- The standard basis -------------------------------------------------------

E1 = (1.0, 0.0)
E2 = (0.0, 1.0)

# -- The unit square, listed counter-clockwise from the origin ----------------

UNIT_SQUARE = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]

# -- The flag ------------------------------------------------------------------
#
# A deliberately lopsided shape, so you can tell at a glance whether it has
# been turned, stretched, sheared or mirrored: an L lying on its back.

FLAG = [
    (0.0, 0.0),
    (2.0, 0.0),
    (2.0, 0.5),
    (0.5, 0.5),
    (0.5, 2.0),
    (0.0, 2.0),
]

# -- The picture you read a matrix off of --------------------------------------
#
# Exercise 2 describes a drawing in words rather than showing it, because the
# skill is going from "where did the basis vectors land" to "what is the
# matrix", and a picture would let you skip the step. The drawing shows:
#
#     the arrow (1, 0) redrawn ending at (3, 1)
#     the arrow (0, 1) redrawn ending at (-1, 2)

PICTURE_E1_LANDS_AT = (3.0, 1.0)
PICTURE_E2_LANDS_AT = (-1.0, 2.0)

# -- The transformations the exercises use ------------------------------------

SCALE_X, SCALE_Y = 2.0, 3.0
SHEAR_K = 2.0

# -- The collapse --------------------------------------------------------------
#
# Look at the two columns before you predict anything about this one.

COLLAPSE_MATRIX = [[1.0, 2.0], [2.0, 4.0]]

# -- Tolerance -----------------------------------------------------------------
#
# Why 1e-12 and not equality: cos(pi / 2) is 6.123233995736766e-17 in binary
# floating point, not 0.0, and sin(pi / 6) is 0.49999999999999994, not 0.5.
# Both are about 1e-17 from the exact answer. 1e-12 sits five orders of
# magnitude above that error and about four below the smallest quantity this
# lab cares about, so it accepts the rounding and would still catch a genuinely
# wrong answer.
TOL = 1e-12

"""Your running score. Run from the LAB DIRECTORY:

    .venv/bin/pytest starter -q

Anything you have not written yet is SKIPPED, not failed. A skip means "not
attempted"; a failure means "attempted and wrong", and the failure prints both
your answer and the real one.

Float comparisons use a tolerance of TOL, stated below, and `shapes.py`
explains why that number and not equality.
"""

import math

import numpy as np
import pytest

import answers
import shapes
from transforms import (
    SingularMatrix,
    apply,
    columns_of,
    compose,
    determinant,
    from_landings,
    identity,
    inverse,
    reflection_in_x_axis,
    rotation,
    scaling,
    shear_x,
    signed_area,
    transform_polygon,
)

TOL = shapes.TOL


def written(fn, *args, **kwargs):
    """Run part of your work, or skip the test if it is not written yet."""
    try:
        return fn(*args, **kwargs)
    except NotImplementedError as exc:
        pytest.skip(f"not written yet: {exc}")


def predicted(name):
    """Read one prediction from answers.py, or skip if it is still None."""
    value = getattr(answers, name)
    if value is None:
        pytest.skip(f"answers.{name} is still unanswered")
    return value


def close(a, b, tol=TOL):
    """Elementwise closeness for points and for 2 by 2 matrices."""
    if isinstance(a[0], (list, tuple)):
        return all(close(ra, rb, tol) for ra, rb in zip(a, b))
    return all(abs(x - y) <= tol for x, y in zip(a, b))


# -- Exercise 0: the environment ---------------------------------------------


def test_0_the_environment_is_ready():
    """Always passes once the install worked. Everything below is your work."""
    assert np.__version__, "numpy is importable"
    assert identity() == [[1.0, 0.0], [0.0, 1.0]], "the written-for-you helpers load"
    assert abs(signed_area(shapes.UNIT_SQUARE) - 1.0) <= TOL


# -- Exercise 1: your transforms.py ------------------------------------------
#
# Every test below runs its whole body inside `written(...)`, so a test is
# skipped if ANY function it needs is still unwritten -- not just the first
# one. Arguments in Python are evaluated before the call, so gating on one
# function while calling another inside the arguments would let a
# NotImplementedError escape and be reported as a failure. It would say
# "attempted and wrong" about work you had not attempted, which is precisely
# the lie this suite exists to avoid.


def test_1_1_from_landings():
    M = written(lambda: from_landings((3.0, 1.0), (-1.0, 2.0)))
    assert M == [[3.0, -1.0], [1.0, 2.0]], (
        "the landing places are the COLUMNS, so they read downwards"
    )


def test_1_2_columns_of():
    e1, e2 = written(lambda: columns_of([[3.0, -1.0], [1.0, 2.0]]))
    assert tuple(e1) == (3.0, 1.0)
    assert tuple(e2) == (-1.0, 2.0)


def test_1_2_columns_of_undoes_from_landings():
    def check():
        for a, b in [((2.0, -5.0), (0.5, 7.0)), ((0.0, 1.0), (1.0, 0.0))]:
            got_a, got_b = columns_of(from_landings(a, b))
            assert tuple(got_a) == a and tuple(got_b) == b

    written(check)


def test_1_3_apply():
    assert written(lambda: apply([[3.0, -1.0], [1.0, 2.0]], (2.0, 1.0))) == (5.0, 4.0)


def test_1_3_apply_to_a_basis_vector_reads_a_column():
    M = [[3.0, -1.0], [1.0, 2.0]]
    assert tuple(written(lambda: apply(M, shapes.E1))) == (3.0, 1.0)
    assert tuple(apply(M, shapes.E2)) == (-1.0, 2.0)


def test_1_3_apply_matches_numpy():
    def check():
        M = [[3.0, -1.0], [1.0, 2.0]]
        for point in [(2.0, 1.0), (-1.5, 4.0), (0.0, 0.0)]:
            assert close(apply(M, point), (np.array(M) @ np.array(point)).tolist())

    written(check)


def test_1_4_scaling():
    assert written(lambda: scaling(2.0, 3.0)) == [[2.0, 0.0], [0.0, 3.0]]
    assert tuple(written(lambda: apply(scaling(2.0, 3.0), (1.0, 1.0)))) == (2.0, 3.0)


def test_1_5_reflection():
    assert written(reflection_in_x_axis) == [[1.0, 0.0], [0.0, -1.0]]
    assert tuple(written(lambda: apply(reflection_in_x_axis(), (5.0, 0.0)))) == (
        5.0,
        0.0,
    ), "a point on the mirror line cannot move"


def test_1_6_shear():
    assert written(lambda: shear_x(2.0)) == [[1.0, 2.0], [0.0, 1.0]]
    assert tuple(written(lambda: apply(shear_x(2.0), (5.0, 0.0)))) == (5.0, 0.0), (
        "height 0 means no sideways push at all"
    )
    assert tuple(apply(shear_x(2.0), (1.0, 1.0))) == (3.0, 1.0)


def test_1_7_rotation_quarter_turn():
    Q = written(lambda: rotation(math.pi / 2))
    assert close(Q, [[0.0, -1.0], [1.0, 0.0]]), (
        "compared with a tolerance, because cos(pi / 2) is not exactly 0.0"
    )


def test_1_7_rotation_thirty_degrees():
    R = written(lambda: rotation(math.radians(30)))
    assert close(R, [[math.sqrt(3) / 2, -0.5], [0.5, math.sqrt(3) / 2]])


def test_1_7_rotation_preserves_length():
    def check():
        for point in [(3.0, 4.0), (-2.0, 7.0)]:
            turned = apply(rotation(0.9), point)
            assert abs(math.hypot(*point) - math.hypot(*turned)) <= TOL

    written(check)


def test_1_8_compose_matches_the_hand_worked_product():
    C = written(lambda: compose(rotation(math.pi / 2), shear_x(2.0)))
    assert close(C, [[0.0, -1.0], [1.0, 2.0]])


def test_1_8_compose_reproduces_two_separate_steps():
    def check():
        A, B = shear_x(2.0), rotation(math.pi / 2)
        two_steps = transform_polygon(B, transform_polygon(A, shapes.FLAG))
        at_once = transform_polygon(compose(B, A), shapes.FLAG)
        for p, q in zip(two_steps, at_once):
            assert close(p, q)

    written(check)


def test_1_8_compose_matches_numpy_matmul():
    def check():
        A, B = shear_x(2.0), rotation(math.pi / 2)
        assert close(compose(B, A), (np.array(B) @ np.array(A)).tolist())

    written(check)


def test_1_9_determinant():
    assert written(lambda: determinant([[3.0, -1.0], [1.0, 2.0]])) == 7.0, (
        "compute a*d - b*c directly; on whole numbers it should be exact"
    )
    assert determinant([[1.0, 2.0], [2.0, 4.0]]) == 0.0
    assert determinant(identity()) == 1.0


def test_1_9_determinant_is_the_area_factor():
    def check():
        for M, expected in [(scaling(2.0, 3.0), 6.0), (reflection_in_x_axis(), -1.0)]:
            area = signed_area(transform_polygon(M, shapes.UNIT_SQUARE))
            assert abs(area - expected) <= TOL
            assert abs(determinant(M) - expected) <= TOL

    written(check)


def test_1_10_inverse():
    assert written(lambda: inverse([[1.0, 2.0], [0.0, 1.0]])) == [[1.0, -2.0], [0.0, 1.0]]


def test_1_10_inverse_undoes_the_original():
    def check():
        for M in [shear_x(2.0), scaling(2.0, 3.0), [[3.0, -1.0], [1.0, 2.0]]]:
            assert close(compose(inverse(M), M), identity())

    written(check)


def test_1_10_inverse_refuses_a_collapse():
    written(lambda: inverse(identity()))
    with pytest.raises(SingularMatrix):
        inverse(shapes.COLLAPSE_MATRIX)


def test_1_10_that_refusal_is_catchable_as_a_ValueError():
    written(lambda: inverse(identity()))
    with pytest.raises(ValueError):
        inverse(shapes.COLLAPSE_MATRIX)


def test_1_10_inverse_matches_numpy():
    def check():
        for M in [shear_x(2.0), scaling(2.0, 3.0), [[3.0, -1.0], [1.0, 2.0]]]:
            assert close(inverse(M), np.linalg.inv(np.array(M)).tolist())

    written(check)


# -- Exercise 2: reading a matrix off a picture ------------------------------


def test_2_1_picture_matrix():
    assert predicted("PICTURE_MATRIX") == [[3.0, -1.0], [1.0, 2.0]]


def test_2_2_where_2_1_lands():
    assert tuple(predicted("PICTURE_SENDS_2_1_TO")) == (5.0, 4.0)


def test_2_3_neither_row_is_a_landing():
    assert predicted("WHICH_ROW_IS_A_LANDING") == "neither"


# -- Exercise 3: the four standard transformations ---------------------------


def test_3_1_scaling():
    assert tuple(predicted("SCALE_SENDS_1_1_TO")) == (2.0, 3.0)


def test_3_2_reflection():
    assert tuple(predicted("FLIP_SENDS_2_3_TO")) == (2.0, -3.0)


def test_3_3_shear():
    assert tuple(predicted("SHEAR_SENDS_1_1_TO")) == (3.0, 1.0)


def test_3_4_shear_on_the_axis():
    assert tuple(predicted("SHEAR_SENDS_5_0_TO")) == (5.0, 0.0)


def test_3_5_quarter_turn():
    assert close(tuple(predicted("QUARTER_TURN_SENDS_1_0_TO")), (0.0, 1.0))


def test_3_6_cosine_is_not_exactly_zero():
    assert predicted("COS_OF_QUARTER_TURN_IS_EXACTLY_ZERO") is False
    assert math.cos(math.pi / 2) != 0.0


def test_3_7_the_origin_is_fixed_by_every_matrix():
    assert tuple(predicted("THE_POINT_NO_MATRIX_CAN_MOVE")) == (0.0, 0.0)


# -- Exercise 4: linearity ----------------------------------------------------


def test_4_1_T_of_u_plus_v():
    assert tuple(predicted("T_OF_U_PLUS_V")) == (8.0, 3.0)


def test_4_2_T_of_u_plus_T_of_v():
    assert tuple(predicted("T_OF_U_PLUS_T_OF_V")) == (8.0, 3.0)


def test_4_3_f_of_u_plus_v():
    assert tuple(predicted("F_OF_U_PLUS_V")) == (9.0, 4.0)


def test_4_4_f_of_u_plus_f_of_v():
    assert tuple(predicted("F_OF_U_PLUS_F_OF_V")) == (10.0, 5.0)


def test_4_5_the_gap_is_exactly_b():
    assert tuple(predicted("THE_GAP_BETWEEN_THEM")) == (1.0, 1.0)


def test_4_6_f_is_not_linear():
    assert predicted("F_IS_LINEAR") is False


def test_4_7_f_moves_the_origin():
    assert tuple(predicted("F_OF_THE_ORIGIN")) == (1.0, 1.0)


# -- Exercise 5: composition and order ---------------------------------------


def test_5_1_which_order():
    assert predicted("SHEAR_THEN_ROTATE_IS") == "compose(B, A)"


def test_5_2_the_composite_matrix():
    assert close(predicted("SHEAR_THEN_ROTATE_MATRIX"), [[0.0, -1.0], [1.0, 2.0]])


def test_5_3_order_matters():
    assert predicted("BOTH_ORDERS_AGREE") is False


def test_5_4_determinants_multiply():
    assert abs(predicted("DET_OF_THE_COMPOSITE") - 1.0) <= TOL


# -- Exercise 6: determinant, rank and the inverse ---------------------------


def test_6_1_area_after_scaling():
    assert abs(predicted("AREA_AFTER_SCALING") - 6.0) <= TOL


def test_6_2_signed_area_after_reflection():
    assert abs(predicted("SIGNED_AREA_AFTER_REFLECTION") + 1.0) <= TOL


def test_6_3_what_a_negative_determinant_means():
    assert predicted("A_NEGATIVE_DETERMINANT_MEANS") == "the plane was flipped over"


def test_6_4_shear_preserves_area():
    assert abs(predicted("DET_OF_SHEAR") - 1.0) <= TOL


def test_6_5_determinant_of_the_collapse():
    assert abs(predicted("DET_OF_COLLAPSE")) <= TOL


def test_6_6_rank_of_the_collapse():
    assert predicted("RANK_OF_COLLAPSE") == 1
    assert int(np.linalg.matrix_rank(np.array(shapes.COLLAPSE_MATRIX))) == 1


def test_6_7_the_line_everything_lands_on():
    m = predicted("COLLAPSE_LANDS_EVERYTHING_ON_THE_LINE_Y_EQUALS")
    M = np.array(shapes.COLLAPSE_MATRIX)
    for point in [(1.0, 0.0), (0.0, 1.0), (3.0, -1.0)]:
        x, y = (M @ np.array(point)).tolist()
        assert abs(y - m * x) <= TOL


def test_6_8_the_exception_numpy_raises():
    cls = predicted("COLLAPSE_INVERSE_EXCEPTION")
    assert cls is np.linalg.LinAlgError
    with pytest.raises(cls):
        np.linalg.inv(np.array(shapes.COLLAPSE_MATRIX))


def test_6_9_inverse_of_a_shear():
    assert abs(predicted("INVERSE_OF_SHEAR_IS_SHEAR_WITH_K") + 2.0) <= TOL


def test_6_10_inverse_of_a_scaling():
    assert close(predicted("INVERSE_OF_SCALING_2_4"), [[0.5, 0.0], [0.0, 0.25]])

"""Exercise 1 -- your linear transformations, built from arithmetic alone.

Ten functions to write. Each one has a docstring saying exactly what it must
do, a worked example you can check on paper, and a `raise NotImplementedError`
to delete when you write it.

Check yourself as you go, from the LAB DIRECTORY:

    .venv/bin/pytest starter -q

Anything you have not written yet is SKIPPED rather than failed. A skip means
"not attempted"; a failure means "attempted and wrong", and it prints your
answer beside the real one.

Use only `math` from the standard library in this file. NumPy appears in the
tests, where it checks your work -- which is the right way round. Writing it
yourself and then having a mature library agree is worth far more than either
half on its own.

The one idea everything here is built on:

    A matrix IS a function. Its columns are where the basis vectors land.
"""

from __future__ import annotations

import math

Point = tuple[float, float]
Matrix = list[list[float]]


class SingularMatrix(ValueError):
    """Raised when a matrix has no inverse.

    Written for you. It subclasses ValueError to match NumPy, whose
    `numpy.linalg.LinAlgError` is itself a ValueError -- so `except ValueError`
    catches your version and NumPy's alike.
    """


# -- Exercise 1.1 -------------------------------------------------------------


def from_landings(e1_lands_at: Point, e2_lands_at: Point) -> Matrix:
    """Build the 2 by 2 matrix that sends (1, 0) and (0, 1) to these two points.

    A matrix is written as a list of ROWS, but it MEANS its columns. So if
    (1, 0) lands at (3, 1) and (0, 1) lands at (-1, 2), the matrix is

        [[3, -1],
         [1,  2]]

    -- the first landing place read downwards in the left-hand column, the
    second read downwards in the right-hand column.

    >>> from_landings((3.0, 1.0), (-1.0, 2.0))
    [[3.0, -1.0], [1.0, 2.0]]
    """
    raise NotImplementedError("exercise 1.1: from_landings")


# -- Exercise 1.2 -------------------------------------------------------------


def columns_of(matrix: Matrix) -> tuple[Point, Point]:
    """Return the two columns as points: (where e1 lands, where e2 lands).

    The exact inverse of exercise 1.1. For [[a, b], [c, d]] the answer is
    ((a, c), (b, d)).

    >>> columns_of([[3.0, -1.0], [1.0, 2.0]])
    ((3.0, 1.0), (-1.0, 2.0))
    """
    raise NotImplementedError("exercise 1.2: columns_of")


# -- Exercise 1.3 -------------------------------------------------------------


def apply(matrix: Matrix, point: Point) -> Point:
    """Send one point through the transformation.

    Write it the way the idea is stated, not the way a textbook writes it: the
    answer is x lots of the first column, plus y lots of the second.

        apply(M, (x, y)) = x * (where e1 landed) + y * (where e2 landed)

    That is the same arithmetic as the row-by-row rule and it says out loud why
    the columns are the landing places. Hint: `columns_of` is already written
    by the time you get here.

    >>> apply([[3.0, -1.0], [1.0, 2.0]], (2.0, 1.0))
    (5.0, 4.0)
    """
    raise NotImplementedError("exercise 1.3: apply")


# -- Exercise 1.4 -------------------------------------------------------------


def scaling(sx: float, sy: float) -> Matrix:
    """Stretch by sx across and sy up.

    Derive it, do not look it up. Where does (1, 0) -- one step right -- go
    when the plane is stretched sx times horizontally? Where does (0, 1) go?
    Write the two answers down as columns with `from_landings`.

    >>> scaling(2.0, 3.0)
    [[2.0, 0.0], [0.0, 3.0]]
    """
    raise NotImplementedError("exercise 1.4: scaling")


# -- Exercise 1.5 -------------------------------------------------------------


def reflection_in_x_axis() -> Matrix:
    """Mirror the plane in the horizontal axis: up becomes down.

    Derive it. (1, 0) lies ON the mirror line, so ask yourself whether it can
    move at all. (0, 1) is one step up -- where is its reflection?

    >>> reflection_in_x_axis()
    [[1.0, 0.0], [0.0, -1.0]]
    """
    raise NotImplementedError("exercise 1.5: reflection_in_x_axis")


# -- Exercise 1.6 -------------------------------------------------------------


def shear_x(k: float) -> Matrix:
    """Push the plane sideways by k times its height.

    A shear slides each point sideways in proportion to how high it is -- the
    deck of cards pushed over, where the bottom card does not move and the top
    one moves furthest.

    Derive it. (1, 0) has height 0, so how far is it pushed? (0, 1) has height
    1, so how far is it pushed, and where does it end up?

    >>> shear_x(2.0)
    [[1.0, 2.0], [0.0, 1.0]]
    """
    raise NotImplementedError("exercise 1.6: shear_x")


# -- Exercise 1.7 -------------------------------------------------------------


def rotation(theta: float) -> Matrix:
    """Turn the plane anticlockwise by theta RADIANS about the origin.

    The derivation, which needs the unit circle and nothing else:

      * draw a circle of radius 1 about the origin;
      * start at (1, 0) and walk anticlockwise around the rim until you have
        turned through theta;
      * the coordinates of where you now stand are, BY DEFINITION,
        (cos(theta), sin(theta)). That is what cosine and sine ARE.

    So (1, 0) lands at (cos(theta), sin(theta)).

    For (0, 1): it is (1, 0) already turned a quarter turn, so after turning by
    theta it sits a quarter turn ahead of the first landing place. A quarter
    turn anticlockwise sends any (x, y) to (-y, x). Apply that to
    (cos(theta), sin(theta)) and you have the second column.

    Radians: an angle measured as distance walked around the rim of that unit
    circle. A full turn is 2 * pi, so a quarter turn is pi / 2. Use
    `math.cos`, `math.sin`, and `math.radians` if you would rather think in
    degrees.

    >>> [round(v, 10) for row in rotation(math.pi / 2) for v in row]
    [0.0, -1.0, 1.0, 0.0]
    """
    raise NotImplementedError("exercise 1.7: rotation")


# -- Exercise 1.8 -------------------------------------------------------------


def compose(second: Matrix, first: Matrix) -> Matrix:
    """Return the single matrix that does `first` and then `second`.

    Mind the argument order, and mind that it is not arbitrary. Applying
    `first` and then `second` to a vector v is

        second @ (first @ v)

    so the single matrix that does both is `second @ first`: the step that
    happens FIRST is written on the RIGHT, because it is the one standing next
    to the vector.

    Build it the day's way rather than with the row-by-column rule: column 0 of
    the answer is wherever (1, 0) ends up after BOTH steps, and column 1 is
    wherever (0, 1) ends up. Two calls to `apply` and one to `from_landings`.

    >>> compose([[0.0, -1.0], [1.0, 0.0]], [[1.0, 2.0], [0.0, 1.0]])
    [[0.0, -1.0], [1.0, 2.0]]
    """
    raise NotImplementedError("exercise 1.8: compose")


# -- Exercise 1.9 -------------------------------------------------------------


def determinant(matrix: Matrix) -> float:
    """The factor by which this transformation multiplies area, sign included.

    For [[a, b], [c, d]] it is a*d - b*c. Compute exactly that -- one multiply,
    one multiply, one subtract -- so that whole-number inputs give an exact
    whole-number answer.

    What the number means, which matters more than the formula:

      * its size is the area factor: a unit square of area 1 comes out with
        area |determinant|;
      * negative means the plane was flipped over;
      * zero means the plane was flattened onto a line, and nothing can undo it.

    >>> determinant([[3.0, -1.0], [1.0, 2.0]])
    7.0
    >>> determinant([[1.0, 2.0], [2.0, 4.0]])
    0.0
    """
    raise NotImplementedError("exercise 1.9: determinant")


# -- Exercise 1.10 ------------------------------------------------------------


def inverse(matrix: Matrix) -> Matrix:
    """The transformation that undoes this one.

    For 2 by 2, the inverse of [[a, b], [c, d]] is

        1 / determinant  *  [[d, -b], [-c, a]]

    and it exists exactly when the determinant is not zero -- which is the same
    sentence as "exactly when no area was destroyed". If two starting points
    landed on the same place, no rule could send that place back to both, so
    there is nothing to return.

    So: compute the determinant first. If it is zero, `raise SingularMatrix`
    with a message that says why. Otherwise divide each of the four rearranged
    entries by it.

    >>> inverse([[1.0, 2.0], [0.0, 1.0]])
    [[1.0, -2.0], [0.0, 1.0]]
    """
    raise NotImplementedError("exercise 1.10: inverse")


# =============================================================================
# Written for you below this line -- read them, they are used by the tests.
# =============================================================================


def identity() -> Matrix:
    """The do-nothing transformation."""
    return [[1.0, 0.0], [0.0, 1.0]]


def transform_polygon(matrix: Matrix, polygon: list[Point]) -> list[Point]:
    """Send every corner of a polygon through the transformation.

    Corners are enough, because a linear transformation sends straight lines to
    straight lines -- which is exactly what makes it cheap, and exactly what
    stops any stack of them ever drawing a curve.
    """
    return [apply(matrix, point) for point in polygon]


def signed_area(polygon: list[Point]) -> float:
    """The area of a polygon, negative if its corners run clockwise.

    The shoelace formula, given to you because measuring polygons is not what
    today is about: walk the corners in order and for each edge add
    x_here * y_next - x_next * y_here, then halve the total.

    The SIGN is the part this lab uses. Corners listed anticlockwise give a
    positive area; a transformation with a negative determinant turns the shape
    over and the same corners now run clockwise, so the answer comes out
    negative. That is "the plane was flipped" measured rather than asserted.
    """
    total = 0.0
    count = len(polygon)
    for i in range(count):
        x_here, y_here = polygon[i]
        x_next, y_next = polygon[(i + 1) % count]
        total += x_here * y_next - x_next * y_here
    return total / 2.0


def rank(matrix: Matrix, tol: float = 1e-12) -> int:
    """How many dimensions survive: 2 fills the plane, 1 a line, 0 the origin."""
    if abs(determinant(matrix)) > tol:
        return 2
    if any(abs(entry) > tol for row in matrix for entry in row):
        return 1
    return 0

#!/usr/bin/env bash
# Tests for the Day 102 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:
#
#   * a matrix's COLUMNS are where the basis vectors land, and the two landings
#     alone determine where every other vector goes;
#   * scaling, reflection, shear and rotation each come out of that one
#     question, and each agrees with NumPy;
#   * a matrix preserves addition and scalar multiplication, and "matrix plus a
#     constant" fails BOTH -- by exactly the constant, and by exactly (s - 1)
#     times the constant, which is asserted rather than described;
#   * a quarter turn's cosine is 6.123233995736766e-17 and not 0.0, so every
#     float check here states a tolerance and the harness asserts the
#     inexactness itself;
#   * composing two transformations is one matrix product, BA means A first,
#     and AB is a different transformation;
#   * the determinant IS the signed area of the transformed unit square --
#     positive, negative and zero cases all measured;
#   * a singular matrix's inverse raises numpy.linalg.LinAlgError with the
#     message "Singular matrix", and the exact class is asserted;
#   * twenty stacked linear layers collapse to one 2 by 2 matrix;
#   * nothing is left behind on disk.
#
# Everything runs offline. Nothing binds a port, nothing writes outside the
# lab, 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 102 — Where Do the Basis Vectors Land?"
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/^/  /'

pinned_numpy="$(grep -E '^numpy==' "${lab_dir}/requirements/requirements.txt" | cut -d= -f3)"
installed_numpy="$("${python_bin}" -c "from importlib.metadata import version; print(version('numpy'))")"
check_eq "installed numpy matches requirements.txt" "${pinned_numpy}" "${installed_numpy}"

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_columns_are_landings 02_building_the_transformations \
              03_linear_or_not 04_composition_and_order \
              05_determinant_inverse_rank 06_the_limit_of_linear; 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, real exceptions"
# --------------------------------------------------------------------------

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 75 ]; then
  check "the reference suite ran at least 75 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 75 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 `transforms` and
# `shapes`, 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 math
import random

import numpy as np

import shapes
from transforms import (
    apply,
    columns_of,
    compose,
    determinant,
    inverse,
    is_linear,
    preserves_addition,
    preserves_scaling,
    rank,
    reflection_in_x_axis,
    rotation,
    scaling,
    shear_x,
    signed_area,
    transform_polygon,
    SingularMatrix,
)

TOL = shapes.TOL

M = shapes.PICTURE_MATRIX
print("columns", columns_of(M))
print("row0_is_not_a_landing", tuple(M[0]) != columns_of(M)[0])
print("sends_2_1", apply(M, (2.0, 1.0)))
print("det_picture_scratch", determinant(M))
print("det_picture_numpy_exact", float(np.linalg.det(np.array(M))) == 7.0)
print("det_picture_numpy_close", abs(float(np.linalg.det(np.array(M))) - 7.0) < 1e-14)

print("scale", scaling(2.0, 3.0))
print("flip", reflection_in_x_axis())
print("shear", shear_x(2.0))
print("shear_on_axis", apply(shear_x(2.0), (5.0, 0.0)))

Q = rotation(math.pi / 2)
print("cos_quarter_exact_zero", math.cos(math.pi / 2) == 0.0)
print("quarter_within_tol", all(
    abs(a - b) <= TOL for a, b in zip(apply(Q, shapes.E1), (0.0, 1.0))
))
print("sin_thirty_exact_half", math.sin(math.radians(30)) == 0.5)

S = shapes.SCALE_MATRIX
b = (1.0, 1.0)
u, v, s = (1.0, 2.0), (3.0, -1.0), 5.0
linear = lambda p: apply(S, p)
affine = lambda p: tuple(c + o for c, o in zip(apply(S, p), b))
print("linear_is_linear", is_linear(linear, u, v, s, TOL))
print("affine_is_linear", is_linear(affine, u, v, s, TOL))
ok, together, separately = preserves_addition(affine, u, v, TOL)
print("affine_add_gap", tuple(round(x - y, 12) for x, y in zip(separately, together)))
ok, first, after = preserves_scaling(affine, u, s, TOL)
print("affine_scale_gap", tuple(round(x - y, 12) for x, y in zip(after, first)))
print("linear_fixes_origin", linear((0.0, 0.0)) == (0.0, 0.0))
print("affine_moves_origin", affine((0.0, 0.0)) == b)

A, B = shear_x(2.0), Q
BA, AB = compose(B, A), compose(A, B)
print("BA_matches_hand", all(
    abs(x - y) <= TOL
    for r1, r2 in zip(BA, shapes.SHEAR_THEN_ROTATE) for x, y in zip(r1, r2)
))
print("AB_matches_hand", all(
    abs(x - y) <= TOL
    for r1, r2 in zip(AB, shapes.ROTATE_THEN_SHEAR) for x, y in zip(r1, r2)
))
print("orders_differ", any(
    abs(x - y) > TOL for r1, r2 in zip(BA, AB) for x, y in zip(r1, r2)
))
two_steps = transform_polygon(B, transform_polygon(A, shapes.FLAG))
at_once = transform_polygon(BA, shapes.FLAG)
print("one_matrix_equals_two_steps", all(
    abs(x - y) <= TOL for p, q in zip(two_steps, at_once) for x, y in zip(p, q)
))

for name, matrix, expected in (
    ("area_scale", scaling(2.0, 3.0), 6.0),
    ("area_shear", shear_x(2.0), 1.0),
    ("area_flip", reflection_in_x_axis(), -1.0),
    ("area_collapse", shapes.COLLAPSE_MATRIX, 0.0),
):
    area = signed_area(transform_polygon(matrix, shapes.UNIT_SQUARE))
    print(name, round(area, 12), round(determinant(matrix), 12), expected)

G = shapes.COLLAPSE_MATRIX
print("collapse_on_one_line", all(
    abs(apply(G, p)[1] - 2.0 * apply(G, p)[0]) <= TOL
    for p in [(1.0, 0.0), (0.0, 1.0), (3.0, -1.0)]
))
print("collapse_two_points_one_landing", apply(G, (2.0, 0.0)) == apply(G, (0.0, 1.0)))
print("rank_collapse", rank(G), int(np.linalg.matrix_rank(np.array(G))))
print("rank_zero_matrix", rank([[0.0, 0.0], [0.0, 0.0]]))

print("inverse_of_shear", inverse(shear_x(2.0)) == shear_x(-2.0))
try:
    inverse(G)
except SingularMatrix:
    print("scratch_inverse_refuses", "SingularMatrix")
else:
    print("scratch_inverse_refuses", "NOTHING_RAISED")
try:
    np.linalg.inv(np.array(G))
except Exception as exc:  # deliberately broad: the TYPE is what is asserted
    print("numpy_inverse_refuses", type(exc).__name__, str(exc))
else:
    print("numpy_inverse_refuses", "NOTHING_RAISED", "")
print("linalgerror_is_valueerror", issubclass(np.linalg.LinAlgError, ValueError))

random.seed(102)
stack = [
    [[random.uniform(-2, 2) for _ in range(2)] for _ in range(2)] for _ in range(20)
]
point = (0.7, -0.4)
stepwise = point
for layer in stack:
    stepwise = apply(layer, stepwise)
combined = stack[0]
for layer in stack[1:]:
    combined = compose(layer, combined)
at_once_pt = apply(combined, point)
print("stack_collapses", all(
    abs(a - c) <= 1e-9 * max(1.0, abs(a)) for a, c in zip(stepwise, at_once_pt)
))
print("stack_still_fixes_origin", all(
    abs(c) <= TOL for c in apply(combined, (0.0, 0.0))
))
PY
)"

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

check_eq "the columns are the two landing places" \
  "((3.0, 1.0), (-1.0, 2.0))" "$(get columns)"
check_eq "row 0 is NOT a landing place" "True" "$(get row0_is_not_a_landing)"
check_eq "the matrix sends (2, 1) to the hand-worked (5, 4)" \
  "(5.0, 4.0)" "$(get sends_2_1)"
check_eq "the from-scratch determinant of the picture matrix is exactly 7" \
  "7.0" "$(get det_picture_scratch)"
check_eq "numpy.linalg.det is NOT exactly 7 on the same matrix" \
  "False" "$(get det_picture_numpy_exact)"
check_eq "numpy.linalg.det is within 1e-14 of 7" \
  "True" "$(get det_picture_numpy_close)"

check_eq "scaling(2, 3) is derived correctly" \
  "[[2.0, 0.0], [0.0, 3.0]]" "$(get scale)"
check_eq "reflection in the x axis is derived correctly" \
  "[[1.0, 0.0], [0.0, -1.0]]" "$(get flip)"
check_eq "shear_x(2) is derived correctly" \
  "[[1.0, 2.0], [0.0, 1.0]]" "$(get shear)"
check_eq "a shear leaves a point at height 0 exactly where it was" \
  "(5.0, 0.0)" "$(get shear_on_axis)"

# Section 6 re-runs this script with D102_SELF_TEST=1, which swaps ONE
# expectation below for a deliberately wrong one. That is how the harness
# proves it can fail rather than merely asserting that it could.
expected_cos_exact="False"
if [ -n "${D102_SELF_TEST:-}" ]; then
  expected_cos_exact="True"   # the naive belief, deliberately wrong here
fi
check_eq "cos(pi / 2) is not exactly 0.0, which is why a tolerance is required" \
  "${expected_cos_exact}" "$(get cos_quarter_exact_zero)"
check_eq "the quarter turn still lands on (0, 1) within the stated tolerance" \
  "True" "$(get quarter_within_tol)"
check_eq "sin(30 degrees) is not exactly 0.5 either" \
  "False" "$(get sin_thirty_exact_half)"

check_eq "a matrix is linear on the tested pair" "True" "$(get linear_is_linear)"
check_eq "matrix-plus-a-constant is not linear" "False" "$(get affine_is_linear)"
check_eq "the addition failure is exactly b" "(1.0, 1.0)" "$(get affine_add_gap)"
check_eq "the scaling failure is exactly (s - 1) times b" \
  "(4.0, 4.0)" "$(get affine_scale_gap)"
check_eq "a linear map fixes the origin" "True" "$(get linear_fixes_origin)"
check_eq "an affine map moves it" "True" "$(get affine_moves_origin)"

check_eq "compose(B, A) matches the hand-worked shear-then-rotate" \
  "True" "$(get BA_matches_hand)"
check_eq "compose(A, B) matches the hand-worked rotate-then-shear" \
  "True" "$(get AB_matches_hand)"
check_eq "the two orders are different transformations" "True" "$(get orders_differ)"
check_eq "one composed matrix reproduces the two-step sequence on every corner" \
  "True" "$(get one_matrix_equals_two_steps)"

check_eq "scaling(2, 3) multiplies the unit square's area by 6" \
  "6.0 6.0 6.0" "$(get area_scale)"
check_eq "a shear preserves area exactly" "1.0 1.0 1.0" "$(get area_shear)"
check_eq "a reflection gives a NEGATIVE area of the same size" \
  "-1.0 -1.0 -1.0" "$(get area_flip)"
check_eq "a collapse gives area 0" "0.0 0.0 0.0" "$(get area_collapse)"

check_eq "the collapse puts every vector on the line y = 2x" \
  "True" "$(get collapse_on_one_line)"
check_eq "two different points land on the same place, so nothing can undo it" \
  "True" "$(get collapse_two_points_one_landing)"
check_eq "the collapse has rank 1, and numpy agrees" "1 1" "$(get rank_collapse)"
check_eq "the all-zero matrix has rank 0" "0" "$(get rank_zero_matrix)"

check_eq "the inverse of shear_x(2) is shear_x(-2)" "True" "$(get inverse_of_shear)"
check_eq "the from-scratch inverse refuses a singular matrix" \
  "SingularMatrix" "$(get scratch_inverse_refuses)"
check_eq "numpy.linalg.inv raises LinAlgError with the message 'Singular matrix'" \
  "LinAlgError Singular matrix" "$(get numpy_inverse_refuses)"
check_eq "numpy.linalg.LinAlgError is catchable as a ValueError" \
  "True" "$(get linalgerror_is_valueerror)"

check_eq "twenty stacked linear layers collapse to one 2 by 2 matrix" \
  "True" "$(get stack_collapses)"
check_eq "and the stack still cannot move the origin" \
  "True" "$(get stack_still_fixes_origin)"

# --------------------------------------------------------------------------
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 cos(pi / 2) is exactly 0.0, and asserts that the
# re-run reports the failure and exits non-zero. If this section passes,
# section 5 is not decorative.
if [ -z "${D102_SELF_TEST:-}" ]; then
  self_out="$(D102_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: cos(pi / 2) is not exactly 0.0"*)
      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 left behind"
# --------------------------------------------------------------------------

# `.venv` is pruned from the searches below. A virtual environment ships the
# installed packages' own precompiled bytecode -- hundreds of __pycache__
# directories that came with NumPy or pytest and have nothing to do with
# whether THIS lab tidied up after itself. Without the prune, following the
# README's own setup instructions makes this check fail, which reports a
# problem the reader cannot fix and did not cause.
if find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -print -quit 2>/dev/null | grep -q .; then
  check "no __pycache__ directory anywhere under the lab after a full run" "no"
else
  check "no __pycache__ directory anywhere under the lab after a full run" "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" "no"
else
  check "no .pytest_cache directory left under the lab" "yes"
fi

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 ]