Hands-on lab — Day 105: Transforming Images with Matrices

"""01 — An image is a matrix, and the ordering trap that comes with it.

Run from the examples directory:

    ../.venv/bin/python3 01_an_image_is_a_matrix.py

Nothing is downloaded. The picture is built from arithmetic, its exact pixel
values are asserted, and the (row, column) versus (x, y) mismatch is measured
rather than described.
"""

import numpy as np

import pattern
import warp

SCRIPT = "01_an_image_is_a_matrix.py"


def main():
    img = pattern.make_pattern()

    print("An image IS a matrix. Here is the whole of this lab's test picture,")
    print("printed twice: once as characters, once as the numbers it really is.")
    print()
    print(pattern.as_text(img))
    print()
    print("The same array, as numbers:")
    print()
    for r, row in enumerate(img):
        print(f"  row {r}: " + " ".join(f"{int(v):3d}" for v in row))
    print()

    # ----------------------------------------------------------------------
    print("Shape and dtype")
    print("-" * 60)
    print(f"  img.shape       = {img.shape}   <- (height, width) = (rows, columns)")
    print(f"  img.dtype       = {img.dtype}      <- one byte per pixel, 0 to 255")
    print(f"  img.size        = {img.size}       <- total pixels")
    print(f"  img.nbytes      = {img.nbytes}       <- one byte each, so size == nbytes")
    print()
    assert img.shape == pattern.EXPECTED_SHAPE
    assert img.dtype == np.uint8
    assert img.size == 81
    assert img.nbytes == 81

    # ----------------------------------------------------------------------
    print("The ordering trap, measured")
    print("-" * 60)
    print("  All week a point was written (x, y). A NumPy array is indexed")
    print("  rows first. Those are the same two numbers in the OPPOSITE order.")
    print()
    # The corner mark is at row 8, column 8, which is symmetric and therefore
    # useless for showing the swap. Use a cell where row and column differ.
    row, col = 4, 3  # inside the middle bar of the F
    print(f"  img[{row}, {col}] = {int(img[row, col])}   (row {row}, column {col})")
    print(f"  img[{col}, {row}] = {int(img[col, row])}   (row {col}, column {row})")
    print("  Different pixels. Swapping the two numbers silently reads the")
    print("  wrong place -- it does not raise, it just returns a wrong answer.")
    print()
    assert int(img[row, col]) == pattern.INK
    assert int(img[col, row]) == pattern.PAPER
    assert img[row, col] != img[col, row]

    print("  As a point in the language of Week 15, that same pixel is")
    print(f"    (x, y) = ({col}, {row})   because x is the COLUMN and y is the ROW")
    print("  and to read it back you must swap again: img[y, x].")
    print()
    x, y = col, row
    assert int(img[y, x]) == pattern.INK

    # ----------------------------------------------------------------------
    print("Why the origin is top-left here and bottom-left on Day 102")
    print("-" * 60)
    print("  Row 0 is the first row of memory, and screens have always been")
    print("  drawn top row first. So y = 0 is the TOP and y grows DOWNWARD.")
    print("  Day 102's graphs had y growing upward. The matrices are identical;")
    print("  the picture is flipped relative to the graph paper, which is why a")
    print("  counter-clockwise rotation matrix turns an image clockwise.")
    print()
    top_row_ink = int((img[0] == pattern.INK).sum())
    bottom_row_ink = int((img[8] == pattern.INK).sum())
    print(f"  ink pixels in row 0 (the TOP bar of the F): {top_row_ink}")
    print(f"  ink pixels in row 8 (the bottom of the stem): {bottom_row_ink}")
    print()
    assert top_row_ink == 6, top_row_ink
    assert bottom_row_ink == 2, bottom_row_ink

    # ----------------------------------------------------------------------
    print("The asymmetry, and why the pattern was chosen this way")
    print("-" * 60)
    flipped_lr = np.fliplr(img)
    flipped_ud = np.flipud(img)
    transposed = img.T
    print(f"  equal to its left-right mirror?  {np.array_equal(img, flipped_lr)}")
    print(f"  equal to its up-down mirror?     {np.array_equal(img, flipped_ud)}")
    print(f"  equal to its transpose?          {np.array_equal(img, transposed)}")
    print("  Three Falses. A square would have given three Trues, and a broken")
    print("  flip would then have passed its test. That is the whole reason the")
    print("  test image is an F.")
    print()
    assert not np.array_equal(img, flipped_lr)
    assert not np.array_equal(img, flipped_ud)
    assert not np.array_equal(img, transposed)

    # ----------------------------------------------------------------------
    print("Colour: the same thing, three times")
    print("-" * 60)
    colour = pattern.make_colour_pattern()
    print(f"  colour.shape = {colour.shape}   <- (height, width, 3)")
    print("  The last axis is the channel. Three stacked planes, each one a")
    print("  matrix of exactly the kind printed above.")
    print()
    for c, name in enumerate(("red", "green", "blue")):
        plane = colour[:, :, c]
        print(f"  {name:<6} plane: shape {plane.shape}, "
              f"{int((plane == pattern.INK).sum())} pixels at {pattern.INK}, "
              f"mean {plane.mean():.2f}")
    print()
    assert colour.shape == pattern.EXPECTED_COLOUR_SHAPE
    assert np.array_equal(colour[:, :, 0], img)
    assert np.array_equal(colour[:, :, 1], np.fliplr(img))
    assert int(colour[:, :, 2].min()) == int(colour[:, :, 2].max()) == pattern.MARK

    print("  One pixel of the colour image is three numbers:")
    print(f"    colour[0, 1] = {tuple(int(v) for v in colour[0, 1])}  (red, green, blue)")
    print()
    assert tuple(int(v) for v in colour[0, 1]) == (pattern.INK, pattern.PAPER, pattern.MARK)

    # ----------------------------------------------------------------------
    print("A pixel is a point, and a point can be transformed")
    print("-" * 60)
    print("  That is the whole of the rest of this lab. The corner mark sits at")
    print(f"  row {pattern.MARK_CELL[0]}, column {pattern.MARK_CELL[1]}, so as a")
    print(f"  point it is (x, y) = ({pattern.MARK_CELL[1]}, {pattern.MARK_CELL[0]}).")
    moved = warp.apply_point(warp.translation(2, -3), (8.0, 8.0))
    print(f"  Move it by (+2, -3):  {moved}")
    print("  Nothing was done to any pixel VALUE. The coordinate moved.")
    print()
    assert moved == (10.0, 5.0)

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


if __name__ == "__main__":
    main()

"""02 — Forward mapping, done wrong on purpose, and the holes counted.

Run from the examples directory:

    ../.venv/bin/python3 02_forward_mapping_leaves_holes.py

The obvious way to rotate an image is to walk the input, work out where each
pixel goes, and put it there. This script does exactly that and then counts the
damage. The holes are not a bug in the code below; they are a property of the
method, and no amount of care in the loop removes them.
"""

import math

import numpy as np

import pattern
import warp

SCRIPT = "02_forward_mapping_leaves_holes.py"


def report(label, matrix, img):
    forward, holes = warp.warp_forward(img, matrix, fill=pattern.FILL)
    inverse = warp.warp_nearest(img, matrix, fill=pattern.FILL)
    n_holes = int(holes.sum())
    print(f"{label}")
    print("-" * 60)
    print(f"  forward mapping: {n_holes} of {holes.size} output pixels were never "
          f"written ({100.0 * n_holes / holes.size:.1f}%)")
    print()
    print("  forward mapping (holes show as ~)      inverse mapping")
    left = pattern.as_text(forward).split("\n")
    right = pattern.as_text(inverse).split("\n")
    for a, b in zip(left, right):
        print(f"    {a}                          {b}")
    print()
    return forward, holes, inverse, n_holes


def main():
    img = pattern.make_pattern()
    height, width = img.shape

    print("The idea that does not work")
    print("=" * 60)
    print("  for each INPUT pixel:")
    print("      send its centre through the matrix")
    print("      round to the nearest output pixel")
    print("      write the value there")
    print()
    print("  It reads like the definition of the transformation, and it is.")
    print("  The problem is not the arithmetic. The problem is that the loop is")
    print("  over the WRONG array: it iterates the input, so nothing guarantees")
    print("  every output pixel gets visited.")
    print()

    # ----------------------------------------------------------------------
    # A rotation: pixels spread apart, and gaps open between them.
    # ----------------------------------------------------------------------
    rot30 = warp.about_centre(warp.rotation(math.radians(30)), width, height)
    _, holes30, _, n30 = report("A 30 degree rotation about the centre", rot30, img)
    assert n30 == 22, n30
    assert holes30.sum() > 0

    print("  Look at where the holes are. They are not only round the edges,")
    print("  where the picture genuinely ran out of source. They are INSIDE the")
    print("  glyph -- single missing pixels punched through solid ink. A rotation")
    print("  is area-preserving (its determinant is 1), so it cannot create more")
    print("  room; what it does is land the input pixels on non-integer positions")
    print("  that round unevenly, so some output pixels collect two input pixels")
    print("  and their neighbours collect none.")
    print()
    assert abs(warp.determinant(rot30) - 1.0) <= warp.TOL

    # ----------------------------------------------------------------------
    # Enlarging: the worst case, because there are simply more output pixels
    # than input pixels to fill them.
    # ----------------------------------------------------------------------
    print("Enlarging: the same failure, but arithmetically unavoidable")
    print("=" * 60)
    grow = warp.scaling(2.0, 2.0)
    forward_big, holes_big = warp.warp_forward(
        img, grow, out_shape=(height * 2, width * 2), fill=pattern.FILL
    )
    n_big = int(holes_big.sum())
    print(f"  Scale by 2. Input has {img.size} pixels; output has "
          f"{holes_big.size}.")
    print(f"  At most {img.size} output pixels can ever be written, so at least")
    print(f"  {holes_big.size - img.size} MUST be holes. Measured: {n_big}.")
    print()
    print(pattern.as_text(forward_big))
    print()
    assert n_big >= holes_big.size - img.size
    assert n_big == 243, n_big
    assert holes_big.size == 324

    print("  Three quarters of the output is missing, in a regular lattice. This")
    print("  is a counting argument, not a rounding accident: 81 input pixels")
    print("  cannot fill 324 output pixels however carefully you place them.")
    print()

    # ----------------------------------------------------------------------
    # Shrinking: no holes, but a different loss.
    # ----------------------------------------------------------------------
    print("Shrinking: no holes, and still not right")
    print("=" * 60)
    shrink = warp.scaling(0.5, 0.5)
    forward_small, holes_small = warp.warp_forward(
        img, shrink, out_shape=(5, 5), fill=pattern.FILL
    )
    n_small = int(holes_small.sum())
    print(f"  Scale by a half into a 5 by 5 output: {n_small} holes.")
    print("  No holes at all -- but now several input pixels land on the SAME")
    print("  output pixel and overwrite each other, so which one survives is")
    print("  decided by the order of the loop rather than by the picture.")
    print()
    print(pattern.as_text(forward_small))
    print()
    assert n_small == 0, n_small

    # Prove the overwriting rather than describing it: count how many input
    # pixels land on each output pixel.
    landings = {}
    for y in range(height):
        for x in range(width):
            fx, fy = warp.apply_point(shrink, (x + 0.5, y + 0.5))
            key = (math.floor(fy), math.floor(fx))
            landings.setdefault(key, []).append((y, x))
    collisions = {k: v for k, v in landings.items() if len(v) > 1}
    worst = max(len(v) for v in landings.values())
    print(f"  {len(collisions)} of the 25 output pixels were written more than")
    print(f"  once; the busiest received {worst} input pixels.")
    print(f"  Example: output pixel {sorted(collisions)[0]} was written by input "
          f"pixels {collisions[sorted(collisions)[0]]}.")
    print()
    # 24, not 25. The 9 input rows map to output rows 0,0,1,1,2,2,3,3,4 -- row 8
    # is odd one out and lands alone. The same happens in x, so output pixel
    # (4, 4) is the single one that receives exactly one input pixel.
    assert len(collisions) == 24, len(collisions)
    assert worst == 4, worst
    assert len(landings[(4, 4)]) == 1, landings[(4, 4)]

    # ----------------------------------------------------------------------
    print("Why patching the holes is the wrong instinct")
    print("=" * 60)
    print("  The tempting fix is to find the holes and fill each one from its")
    print("  neighbours. That is more code, it is slower, it needs a second pass")
    print("  over the image, and it still guesses. Turning the loop inside out")
    print("  costs nothing and removes the problem entirely, because a loop over")
    print("  the OUTPUT visits every output pixel exactly once by construction.")
    print()
    print("  That is the next script. The count to beat is:")
    print(f"    30 degree rotation, forward mapping: {n30} holes")
    print()

    inverse30 = warp.warp_nearest(img, rot30, fill=pattern.FILL)
    unset = int((inverse30 == pattern.FILL).sum())
    print(f"  Inverse mapping on the same rotation leaves {unset} pixels at the")
    print("  fill value -- and every one of those is a CORNER whose source lies")
    print("  outside the input image, which is clipping, not a hole. Proof: each")
    print("  one maps back outside the picture.")
    print()

    genuinely_outside = 0
    back = warp.invert(rot30)
    for oy in range(height):
        for ox in range(width):
            if inverse30[oy, ox] == pattern.FILL:
                sx, sy = warp.apply_point(back, (ox + 0.5, oy + 0.5))
                if not (0 <= math.floor(sx) < width and 0 <= math.floor(sy) < height):
                    genuinely_outside += 1
    print(f"  fill-valued output pixels: {unset}")
    print(f"  of those, sourced from outside the input: {genuinely_outside}")
    assert unset == genuinely_outside, (unset, genuinely_outside)
    assert genuinely_outside > 0
    print("  All of them. Inverse mapping has no holes at all.")
    print()

    assert np.array_equal(inverse30, warp.warp_nearest(img, rot30, fill=pattern.FILL))

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


if __name__ == "__main__":
    main()

"""03 — Inverse mapping and nearest-neighbour, with exact answers asserted.

Run from the examples directory:

    ../.venv/bin/python3 03_inverse_mapping.py

The loop is turned inside out. Instead of asking "where does this input pixel
go", ask "where did this output pixel come from". Every output pixel is visited
exactly once, so holes are impossible by construction rather than by care.
"""

import math

import numpy as np

import pattern
import warp

SCRIPT = "03_inverse_mapping.py"


def main():
    img = pattern.make_pattern()
    height, width = img.shape

    print("The loop, inside out")
    print("=" * 60)
    print("  for each OUTPUT pixel:")
    print("      take its centre, (x + 0.5, y + 0.5)")
    print("      send it BACKWARD through the inverse matrix")
    print("      take the value of the input pixel containing that point")
    print()
    print("  Every output pixel is assigned. That is not a claim about the")
    print("  arithmetic; it is a property of iterating over the array you are")
    print("  filling. This is why every real implementation works this way.")
    print()

    # ----------------------------------------------------------------------
    print("A quarter turn, with the answer known in advance")
    print("=" * 60)
    print("  90 degrees is the angle to test first, because the right answer is")
    print("  exact: every pixel lands on a pixel, so this can be asserted with")
    print("  == rather than with a tolerance.")
    print()

    quarter = warp.about_centre(warp.rotation_quarter_turns(1), width, height)
    turned = warp.warp_nearest(img, quarter, fill=pattern.FILL)

    print("  before                     after")
    for a, b in zip(pattern.as_text(img).split("\n"),
                    pattern.as_text(turned).split("\n")):
        print(f"    {a}                {b}")
    print()

    print("  The F turned CLOCKWISE, from a counter-clockwise rotation matrix.")
    print("  That is the y-down coordinate system, not a sign error. Day 102's")
    print("  graphs had y growing upward; here row 0 is the top.")
    print()

    expected = np.rot90(img, -1)
    same = np.array_equal(turned, expected)
    print(f"  identical to numpy.rot90(img, -1)?  {same}")
    print(f"  differing pixels: {int((turned != expected).sum())}")
    print()
    assert same
    assert int((turned != expected).sum()) == pattern.PIXEL_TOL

    print("  No pixel took the fill value, because a quarter turn of a square")
    print("  image lands entirely inside the frame:")
    print(f"    pixels at the fill value: {int((turned == pattern.FILL).sum())}")
    print(f"    ink pixels before: {len(pattern.ink_cells(img))}, "
          f"after: {len(pattern.ink_cells(turned))}")
    print()
    assert int((turned == pattern.FILL).sum()) == 0
    assert len(pattern.ink_cells(turned)) == len(pattern.ink_cells(img)) == 24

    # ----------------------------------------------------------------------
    print("Where individual pixels went, checked one at a time")
    print("=" * 60)
    print("  The corner mark is the easiest thing to follow, because there is")
    print("  exactly one pixel of that value in the whole image.")
    before_mark = tuple(int(v) for v in np.argwhere(img == pattern.MARK)[0])
    after_mark = tuple(int(v) for v in np.argwhere(turned == pattern.MARK)[0])
    print(f"    before: row {before_mark[0]}, column {before_mark[1]}")
    print(f"    after:  row {after_mark[0]}, column {after_mark[1]}")
    print("  Bottom-right to bottom-left, which is what a clockwise quarter turn")
    print("  does to a bottom-right corner.")
    print()
    assert before_mark == (8, 8)
    assert after_mark == (8, 0)

    print("  And the whole top bar of the F, which was row 0, columns 1 to 6:")
    top_bar_after = sorted(
        (int(r), int(c)) for r, c in np.argwhere(turned == pattern.INK)
        if c == 8
    )
    print(f"    it is now column 8, rows {[r for r, _ in top_bar_after]}")
    print()
    assert [r for r, _ in top_bar_after] == list(range(1, 7))

    # ----------------------------------------------------------------------
    print("All four quarter turns")
    print("=" * 60)
    for turns in (1, 2, 3, 4):
        matrix = warp.about_centre(warp.rotation_quarter_turns(turns), width, height)
        out = warp.warp_nearest(img, matrix, fill=pattern.FILL)
        reference = np.rot90(img, -turns)
        ok = np.array_equal(out, reference)
        print(f"  {90 * turns:>3} degrees: matches numpy.rot90(img, -{turns})  ->  {ok}")
        assert ok, turns
    print()

    # ----------------------------------------------------------------------
    print("Does the float version of a quarter turn agree?")
    print("=" * 60)
    print(f"  math.cos(math.pi / 2) = {math.cos(math.pi / 2)!r}")
    print("  Not 0.0 -- the Day 102 result. So the trigonometric rotation matrix")
    print("  is a hair off the exact one. Does it change any pixel?")
    print()
    trig = warp.about_centre(warp.rotation(math.pi / 2), width, height)
    trig_out = warp.warp_nearest(img, trig, fill=pattern.FILL)
    print(f"    matrices identical?      "
          f"{trig == warp.about_centre(warp.rotation_quarter_turns(1), width, height)}")
    print(f"    matrices within {warp.TOL:g}?   "
          f"{warp.matrices_close(trig, quarter)}")
    print(f"    output images identical? {np.array_equal(trig_out, turned)}")
    print()
    assert not warp.matrices_close(trig, quarter, tol=0.0)
    assert warp.matrices_close(trig, quarter)
    assert np.array_equal(trig_out, turned)
    print("  The matrices differ; the images do not. Nearest-neighbour rounds to")
    print("  a whole pixel, and an error of 1e-16 never reaches a rounding")
    print("  boundary. The float noise is real and it is absorbed. That will not")
    print("  be true of every angle, which is why the exact-integer matrix exists")
    print("  and is used wherever the answer is meant to be checkable.")
    print()

    # ----------------------------------------------------------------------
    print("A full turn: exact, and only because it is ONE matrix")
    print("=" * 60)
    full = warp.about_centre(warp.rotation(2.0 * math.pi), width, height)
    full_out = warp.warp_nearest(img, full, fill=pattern.FILL)
    print(f"  rotation(2*pi) applied once, differing pixels: "
          f"{int((full_out != img).sum())}")
    assert np.array_equal(full_out, img)

    twelve = img
    for _ in range(12):
        twelve = warp.warp_nearest(
            twelve,
            warp.about_centre(warp.rotation(math.radians(30)), width, height),
            fill=pattern.FILL,
        )
    n_lost = int((twelve != img).sum())
    print(f"  twelve separate 30 degree rotations, differing pixels: {n_lost}")
    print()
    print(pattern.as_text(twelve))
    print()
    assert n_lost == 16, n_lost
    assert not np.array_equal(twelve, img)

    print("  Both routes are 360 degrees. One is exact and one loses 16 of 81")
    print("  pixels. The difference is not the angle -- it is the number of times")
    print("  the image was RESAMPLED. Each nearest-neighbour pass throws away the")
    print("  sub-pixel position, and twelve passes cannot recover what the first")
    print("  one discarded. Compose the matrices, resample once.")
    print()

    composed = warp.identity()
    for _ in range(12):
        composed = warp.compose(
            warp.about_centre(warp.rotation(math.radians(30)), width, height),
            composed,
        )
    composed_out = warp.warp_nearest(img, composed, fill=pattern.FILL)
    print(f"  the same twelve rotations COMPOSED into one matrix and applied")
    print(f"  once, differing pixels: {int((composed_out != img).sum())}")
    print()
    assert np.array_equal(composed_out, img)
    assert warp.matrices_close(composed, warp.identity(), tol=1e-12)

    # ----------------------------------------------------------------------
    print("Colour is not a new problem")
    print("=" * 60)
    colour = pattern.make_colour_pattern()
    turned_colour = warp.warp_colour(colour, quarter, fill=pattern.FILL)
    print(f"  input  shape {colour.shape}")
    print(f"  output shape {turned_colour.shape}")
    for c, name in enumerate(("red", "green", "blue")):
        ok = np.array_equal(turned_colour[:, :, c], np.rot90(colour[:, :, c], -1))
        print(f"  {name:<6} plane rotated correctly: {ok}")
        assert ok
    print("  Same matrix, three planes. The transformation acts on coordinates,")
    print("  and the three planes share their coordinates.")
    print()
    assert turned_colour.shape == (9, 9, 3)

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


if __name__ == "__main__":
    main()

"""04 — Scale, shear and flip, each checked against a known-exact answer.

Run from the examples directory:

    ../.venv/bin/python3 04_scale_shear_flip.py

Every transformation here is chosen so that the correct output is something
NumPy can produce a different way -- `numpy.fliplr`, `numpy.kron`, a strided
slice. Agreement between two independent routes to the same array is worth more
than any single implementation's say-so.
"""

import numpy as np

import pattern
import warp

SCRIPT = "04_scale_shear_flip.py"


def side_by_side(left, right, gap="        "):
    a = pattern.as_text(left).split("\n")
    b = pattern.as_text(right).split("\n")
    width = max(len(line) for line in a)
    for i in range(max(len(a), len(b))):
        la = a[i] if i < len(a) else ""
        lb = b[i] if i < len(b) else ""
        print(f"    {la:<{width}}{gap}{lb}")


def main():
    img = pattern.make_pattern()
    height, width = img.shape

    # ----------------------------------------------------------------------
    print("Flip: a reflection plus a translation, in one matrix")
    print("=" * 60)
    print("  A bare reflection sends x to -x, which puts the whole picture off")
    print("  the left edge. What is wanted is a reflection about the image's")
    print("  centre line: x becomes width - x. That is a reflection FOLLOWED BY")
    print("  a translation -- and translation is not linear, which is exactly why")
    print("  the matrix is 3 by 3 and not 2 by 2.")
    print()
    print("  flip_horizontal(9):")
    for row in warp.flip_horizontal(width):
        print("      [" + "  ".join(f"{v:5.1f}" for v in row) + " ]")
    print()

    flipped = warp.warp_nearest(img, warp.flip_horizontal(width), fill=pattern.FILL)
    print("  before                  after")
    side_by_side(img, flipped)
    print()
    print(f"  identical to numpy.fliplr(img)?  {np.array_equal(flipped, np.fliplr(img))}")
    print(f"  pixels at the fill value:        {int((flipped == pattern.FILL).sum())}")
    print()
    assert np.array_equal(flipped, np.fliplr(img))
    assert int((flipped == pattern.FILL).sum()) == 0

    flipped_v = warp.warp_nearest(img, warp.flip_vertical(height), fill=pattern.FILL)
    print(f"  flip_vertical matches numpy.flipud(img)? "
          f"{np.array_equal(flipped_v, np.flipud(img))}")
    assert np.array_equal(flipped_v, np.flipud(img))

    print("  Individual pixels, checked by hand:")
    mark_after = tuple(int(v) for v in np.argwhere(flipped == pattern.MARK)[0])
    print(f"    the corner mark was at (row 8, column 8); after the horizontal")
    print(f"    flip it is at (row {mark_after[0]}, column {mark_after[1]}).")
    print("    Column 8 became column 9 - 1 - 8 = 0, and the row did not move.")
    print()
    assert mark_after == (8, 0)

    print("  Two flips are the identity, and the matrices say so before any")
    print("  pixel is touched:")
    twice = warp.compose(warp.flip_horizontal(width), warp.flip_horizontal(width))
    print(f"    flip . flip == identity?  {warp.matrices_close(twice, warp.identity())}")
    print(f"    determinant of one flip:  {warp.determinant(warp.flip_horizontal(width))}")
    print("    A negative determinant is a reflection -- Day 102's signed area,")
    print("    unchanged in size and reversed in orientation.")
    print()
    assert warp.matrices_close(twice, warp.identity())
    assert warp.determinant(warp.flip_horizontal(width)) == -1.0
    assert np.array_equal(
        warp.warp_nearest(flipped, warp.flip_horizontal(width), fill=pattern.FILL), img
    )

    # ----------------------------------------------------------------------
    print("Scale up: exact pixel replication, and no new information")
    print("=" * 60)
    doubled = warp.warp_nearest(
        img, warp.scaling(2.0, 2.0), out_shape=(height * 2, width * 2),
        fill=pattern.FILL,
    )
    print(f"  output shape {doubled.shape}, "
          f"pixels at the fill value: {int((doubled == pattern.FILL).sum())}")
    print()
    print(pattern.as_text(doubled))
    print()
    kron = np.kron(img, np.ones((2, 2), dtype=np.uint8))
    print(f"  identical to numpy.kron(img, ones((2, 2)))?  "
          f"{np.array_equal(doubled, kron)}")
    print()
    assert np.array_equal(doubled, kron)
    assert int((doubled == pattern.FILL).sum()) == 0
    assert doubled.shape == (18, 18)

    print("  Zero holes -- compare script 02, where forward mapping left 243 of")
    print("  these 324 pixels unwritten. And notice what the enlargement did NOT")
    print("  do: it made every pixel into a 2 by 2 block. There are four times as")
    print("  many pixels and exactly as much information. Nearest-neighbour")
    print("  cannot invent detail, and nothing else can either.")
    print(f"    distinct values before: {len(np.unique(img))}, "
          f"after: {len(np.unique(doubled))}")
    print()
    assert len(np.unique(doubled)) == len(np.unique(img)) == 3

    # ----------------------------------------------------------------------
    print("Scale down: information is thrown away, and you can name which")
    print("=" * 60)
    halved = warp.warp_nearest(
        img, warp.scaling(0.5, 0.5), out_shape=(4, 4), fill=pattern.FILL
    )
    print(f"  output shape {halved.shape}")
    print()
    print(pattern.as_text(halved))
    print()
    strided = img[1::2, 1::2]
    print(f"  identical to img[1::2, 1::2]?  {np.array_equal(halved, strided)}")
    print("  Which is to say: nearest-neighbour downscaling by 2 keeps every")
    print("  second pixel starting at index 1, and discards the rest. The odd")
    print("  starting index is the half-pixel sampling offset showing itself --")
    print("  output pixel 0 has its centre at 0.5, which doubles to 1.0.")
    print()
    assert np.array_equal(halved, strided)

    print("  The corner mark did not survive, and that is correct behaviour:")
    print(f"    mark pixels in the input:  {int((img == pattern.MARK).sum())}")
    print(f"    mark pixels in the output: {int((halved == pattern.MARK).sum())}")
    print("    Row 8 and column 8 are both odd-one-out under this sampling, so")
    print("    the pixel at (8, 8) is one of the ones dropped. A single-pixel")
    print("    feature disappearing under a downscale is not a bug; it is what")
    print("    downscaling is. Averaging instead of sampling would have kept a")
    print("    trace of it, which is the argument for the next script's blending.")
    print()
    assert int((img == pattern.MARK).sum()) == 1
    assert int((halved == pattern.MARK).sum()) == 0

    # ----------------------------------------------------------------------
    print("Shear: each row slid sideways in proportion to its own y")
    print("=" * 60)
    print("  shear_x(k) sends (x, y) to (x + k*y, y). Row 0 has y = 0, so the")
    print("  textbook says row 0 does not move. Watch what the half-pixel")
    print("  sampling offset does to that claim.")
    print()

    for k, out_w in ((0.5, 14), (1.0, 18), (2.0, 27)):
        sheared = warp.warp_nearest(
            img, warp.shear_x(k), out_shape=(height, out_w), fill=pattern.FILL
        )
        # Where did the left edge of the stem (input column 1) end up in row 0?
        row0 = sheared[0]
        first_ink = int(np.argmax(row0 == pattern.INK))
        shift = first_ink - 1
        print(f"  shear_x({k}): row 0's first ink pixel is at column {first_ink}, "
              f"a shift of {shift}")
        assert warp.determinant(warp.shear_x(k)) == 1.0

    print()
    print("  k = 0.5 and k = 1.0 leave row 0 alone; k = 2.0 moves it by one whole")
    print("  pixel. The reason is arithmetic, not magic. The output pixel in row 0")
    print("  is sampled at its CENTRE, y = 0.5, not at y = 0:")
    print("    x_source = floor((x + 0.5) - k * 0.5)")
    for k in (0.5, 1.0, 2.0):
        inside = 0.5 - k * 0.5
        shift = -int(np.floor(inside))
        print(f"    k = {k}:  floor(x {inside:+.2f}) = x {-shift:+d}"
              f"   ->  the content moves right by {shift}")
    print()
    print("  Day 102 flagged this as the open question and deferred it to today.")
    print("  It is settled: the shear term is multiplied by the pixel CENTRE's y,")
    print("  and row 0's centre is at y = 0.5, so a large enough k moves row 0.")
    print()

    sheared2 = warp.warp_nearest(
        img, warp.shear_x(2.0), out_shape=(height, 27), fill=pattern.FILL
    )
    assert int(np.argmax(sheared2[0] == pattern.INK)) == 2, \
        int(np.argmax(sheared2[0] == pattern.INK))
    sheared1 = warp.warp_nearest(
        img, warp.shear_x(1.0), out_shape=(height, 18), fill=pattern.FILL
    )
    assert int(np.argmax(sheared1[0] == pattern.INK)) == 1

    print("  A shear with k = 1.0, drawn:")
    print()
    print(pattern.as_text(sheared1))
    print()
    print("  Straight lines are still straight and the two vertical edges of the")
    print("  stem are still parallel. Every transformation in this lab is AFFINE,")
    print("  and that is exactly what affine guarantees.")
    print()

    print("  Area is unchanged: a shear's determinant is 1, so it cannot lose or")
    print("  gain ink. Counting the surviving ink pixels checks that claim, once")
    print("  the output is wide enough to hold the sheared glyph:")
    ink_before = int((img == pattern.INK).sum())
    ink_after = int((sheared1 == pattern.INK).sum())
    print(f"    ink before: {ink_before}, ink after: {ink_after}")
    print()
    assert ink_before == ink_after == 24

    print("  Shear the other way and it comes back -- with one caveat that is")
    print("  worth more than the rule it breaks.")
    print()
    round_trip = warp.compose(warp.shear_x(-1.0), warp.shear_x(1.0))
    print(f"    as MATRICES, shear_x(-1) . shear_x(1) == identity?  "
          f"{warp.matrices_close(round_trip, warp.identity())}")
    one_pass = warp.warp_nearest(img, round_trip, fill=pattern.FILL)
    print(f"    applied as ONE matrix, differing pixels: "
          f"{int((one_pass != img).sum())}")
    assert warp.matrices_close(round_trip, warp.identity())
    assert np.array_equal(one_pass, img)
    print()

    print("    applied as TWO separate resampling passes:")
    for k, wide in ((0.5, 14), (1.0, 18)):
        out = warp.warp_nearest(
            img, warp.shear_x(k), out_shape=(height, wide), fill=pattern.FILL
        )
        back = warp.warp_nearest(
            out, warp.shear_x(-k), out_shape=(height, width), fill=pattern.FILL
        )
        print(f"      k = {k}:  differing pixels {int((back != img).sum()):>3}")
    print()
    print("    k = 0.5 round-trips exactly. k = 1.0 does not, and loses 28 of 81")
    print("    pixels -- the whole image slides one column left. That is not a")
    print("    bug that was found and left in; it is a boundary case that is")
    print("    worth naming, because it will bite you in real code.")
    print()
    print("    Why: with k = 1.0 the sampled position is (x + 0.5) - 1.0 * 0.5,")
    print("    which is x EXACTLY -- a pixel boundary, where floor has to make an")
    print("    arbitrary choice between two neighbours. With k = 0.5 the position")
    print("    is x + 0.25, safely inside one pixel, and floor is unambiguous.")
    print("    Every integer-valued shear coefficient puts every sample on a")
    print("    boundary at once, so the arbitrary choice is made 81 times in the")
    print("    same direction and the error accumulates into a visible shift.")
    print()
    half_out = warp.warp_nearest(
        img, warp.shear_x(0.5), out_shape=(height, 14), fill=pattern.FILL
    )
    half_back = warp.warp_nearest(
        half_out, warp.shear_x(-0.5), out_shape=(height, width), fill=pattern.FILL
    )
    assert np.array_equal(half_back, img)
    unit_out = warp.warp_nearest(
        img, warp.shear_x(1.0), out_shape=(height, 18), fill=pattern.FILL
    )
    unit_back = warp.warp_nearest(
        unit_out, warp.shear_x(-1.0), out_shape=(height, width), fill=pattern.FILL
    )
    assert int((unit_back != img).sum()) == 28, int((unit_back != img).sum())
    print("    The lesson is the same one script 03 drew from twelve rotations:")
    print("    compose the matrices and resample ONCE. Done that way, this round")
    print("    trip is exact for every k, including 1.0.")
    print()

    # ----------------------------------------------------------------------
    print("What none of this can do")
    print("=" * 60)
    print("  Every matrix in this script is affine: straight lines stay straight,")
    print("  parallel lines stay parallel, and the ratio of lengths along any one")
    print("  line is preserved. Six numbers, and that is the whole family.")
    print()
    print("  What that rules out:")
    print("    * perspective -- railway tracks converging toward a horizon needs")
    print("      a PROJECTIVE transform, whose bottom row is not (0, 0, 1), so")
    print("      the third coordinate stops being 1 and has to be divided out;")
    print("    * lens distortion -- a barrel or pincushion bend is not linear in")
    print("      the coordinates at all, and no matrix of any size expresses it;")
    print("    * warping one face into another -- that is a dense displacement")
    print("      field, a different vector for every pixel.")
    print()
    print("  A quick proof that affine cannot do perspective: an affine map sends")
    print("  parallel lines to parallel lines, because it sends the direction")
    print("  vector of a line through the LINEAR part only, and two lines with")
    print("  the same direction keep the same direction.")
    top_edge = (1.0, 0.0)
    for name, matrix in (
        ("rotation(0.7)", warp.rotation(0.7)),
        ("shear_x(2)", warp.shear_x(2.0)),
        ("scaling(3, 0.5)", warp.scaling(3.0, 0.5)),
    ):
        d0 = warp.apply_point(matrix, (0.0, 0.0))
        d1 = warp.apply_point(matrix, top_edge)
        d2 = warp.apply_point(matrix, (0.0, 5.0))
        d3 = warp.apply_point(matrix, (1.0, 5.0))
        v1 = (d1[0] - d0[0], d1[1] - d0[1])
        v2 = (d3[0] - d2[0], d3[1] - d2[1])
        cross = v1[0] * v2[1] - v1[1] * v2[0]
        print(f"    {name:<16} two parallel edges stay parallel "
              f"(cross product {cross:.1e})")
        assert abs(cross) <= warp.TOL
    print()

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


if __name__ == "__main__":
    main()

"""05 — Homogeneous coordinates, and several transformations folded into one.

Run from the examples directory:

    ../.venv/bin/python3 05_homogeneous_and_composition.py

Day 102 proved that a linear map cannot move the origin, and translation moves
the origin, so translation is not linear and no 2 by 2 matrix performs it. This
script shows the fix -- a third coordinate, fixed at 1 -- and then uses it to
fold a rotation about the image centre into a single matrix.
"""

import math

import numpy as np

import pattern
import warp

SCRIPT = "05_homogeneous_and_composition.py"


def show_matrix(label, matrix):
    print(f"  {label}")
    for row in matrix:
        print("      [" + "  ".join(f"{v:8.4f}" for v in row) + " ]")
    print()


def main():
    img = pattern.make_pattern()
    height, width = img.shape

    # ----------------------------------------------------------------------
    print("The problem: translation is not linear")
    print("=" * 60)
    print("  Day 102's test for linearity was whether the map preserves addition")
    print("  and scalar multiplication, and one consequence is that a linear map")
    print("  must send the origin to the origin. Translation does not.")
    print()

    def move_by_three(point):
        return (point[0] + 3.0, point[1] + 0.0)

    origin = (0.0, 0.0)
    print(f"    moving by (3, 0) sends the origin to {move_by_three(origin)}")
    print("    -- so no 2 by 2 matrix can do it, because every 2 by 2 matrix")
    print("    sends (0, 0) to (0, 0) by construction: the sum of the columns")
    print("    weighted by zero and zero.")
    print()
    for name, matrix in (
        ("rotation(1.1)", warp.rotation(1.1)),
        ("scaling(4, 0.2)", warp.scaling(4.0, 0.2)),
        ("shear_x(9)", warp.shear_x(9.0)),
    ):
        landed = warp.apply_point(matrix, origin)
        print(f"    {name:<16} sends the origin to {landed}")
        assert landed == (0.0, 0.0)
    print()
    assert move_by_three(origin) != origin

    # ----------------------------------------------------------------------
    print("The fix: add a coordinate that is always 1")
    print("=" * 60)
    print("  Write the point (x, y) as the triple (x, y, 1). Now a 3 by 3 matrix")
    print("  can add a constant, because the constant is multiplied by that 1:")
    print()
    print("      [ 1  0  tx ] [ x ]   [ x + tx ]")
    print("      [ 0  1  ty ] [ y ] = [ y + ty ]")
    print("      [ 0  0   1 ] [ 1 ]   [   1    ]")
    print()
    print("  The third coordinate is not a z axis and the picture is not 3-D.")
    print("  It is a bookkeeping device: one extra slot whose only job is to")
    print("  give the translation something to multiply.")
    print()

    shift = warp.translation(3.0, -2.0)
    show_matrix("translation(3, -2):", shift)
    for point in ((0.0, 0.0), (1.0, 1.0), (8.0, 8.0)):
        print(f"    {point} -> {warp.apply_point(shift, point)}")
    print()
    assert warp.apply_point(shift, (0.0, 0.0)) == (3.0, -2.0)
    assert warp.apply_point(shift, (8.0, 8.0)) == (11.0, 6.0)

    print("  And now translation composes with everything else, because it is")
    print("  the same kind of object. A translation's determinant is 1 -- it")
    print("  moves the picture without changing its area -- and its inverse is")
    print("  the opposite translation:")
    print(f"    determinant(translation(3, -2)) = {warp.determinant(shift)}")
    back = warp.invert(shift)
    print(f"    its inverse is translation(-3, 2)?  "
          f"{warp.matrices_close(back, warp.translation(-3.0, 2.0))}")
    print()
    assert warp.determinant(shift) == 1.0
    assert warp.matrices_close(back, warp.translation(-3.0, 2.0))

    print("  On the actual image, translation by (2, 1):")
    moved = warp.warp_nearest(img, warp.translation(2.0, 1.0), fill=pattern.FILL)
    print()
    print("  before                  after")
    a = pattern.as_text(img).split("\n")
    b = pattern.as_text(moved).split("\n")
    for la, lb in zip(a, b):
        print(f"    {la}        {lb}")
    print()
    print("  Two columns right and one row down, with the vacated edges taking")
    print("  the fill value. Check it against a plain NumPy slice, which is what")
    print("  an integer translation ought to be:")
    reference = np.full_like(img, pattern.FILL)
    reference[1:, 2:] = img[:-1, :-2]
    print(f"    matches the slice-and-pad reference?  "
          f"{np.array_equal(moved, reference)}")
    print()
    assert np.array_equal(moved, reference)
    # Two vacated columns (2 * 9 pixels) plus one vacated row (9 pixels), minus
    # the 2 pixels counted twice where the vacated row and columns overlap.
    assert int((moved == pattern.FILL).sum()) == 2 * 9 + 9 - 2

    # ----------------------------------------------------------------------
    print("Composition: three matrices become one")
    print("=" * 60)
    print("  Rotating about the image's centre rather than its top-left corner")
    print("  is three steps: move the centre to the origin, rotate, move it back.")
    print("  It is also ONE matrix, and this is where homogeneous coordinates")
    print("  earn their place -- without them the middle step is a matrix and the")
    print("  two outer steps are not, so they cannot be multiplied together.")
    print()

    cx, cy = width / 2.0, height / 2.0
    to_origin = warp.translation(-cx, -cy)
    turn = warp.rotation_quarter_turns(1)
    back_again = warp.translation(cx, cy)

    show_matrix(f"1. translation({-cx}, {-cy})", to_origin)
    show_matrix("2. rotation, a quarter turn", turn)
    show_matrix(f"3. translation({cx}, {cy})", back_again)

    combined = warp.compose(back_again, turn, to_origin)
    show_matrix("combined = T(+c) . R . T(-c)", combined)

    print("  Read the product RIGHT to LEFT: the rightmost matrix acts first.")
    print("  That is the Day 101 convention and it has not changed.")
    print()
    print(f"  same as about_centre(rotation_quarter_turns(1), 9, 9)?  "
          f"{warp.matrices_close(combined, warp.about_centre(turn, width, height))}")
    print()
    assert warp.matrices_close(combined, warp.about_centre(turn, width, height))

    print("  One matrix must reproduce the three separate steps on every point.")
    print("  Checked on all 81 pixel centres:")
    worst = 0.0
    for y in range(height):
        for x in range(width):
            point = (x + 0.5, y + 0.5)
            stepwise = warp.apply_point(
                back_again, warp.apply_point(turn, warp.apply_point(to_origin, point))
            )
            at_once = warp.apply_point(combined, point)
            worst = max(worst, max(abs(p - q) for p, q in zip(stepwise, at_once)))
    print(f"    largest disagreement over 81 points: {worst:.3e}  "
          f"(tolerance {warp.TOL:g})")
    print()
    assert worst <= warp.TOL

    # ----------------------------------------------------------------------
    print("Order matters, and the images prove it")
    print("=" * 60)
    rotate = warp.about_centre(warp.rotation_quarter_turns(1), width, height)
    stretch = warp.scaling(1.0, 2.0)

    rotate_then_stretch = warp.compose(stretch, rotate)
    stretch_then_rotate = warp.compose(rotate, stretch)

    print(f"  the two products are different matrices?  "
          f"{not warp.matrices_close(rotate_then_stretch, stretch_then_rotate)}")
    print()
    a_img = warp.warp_nearest(
        img, rotate_then_stretch, out_shape=(18, 9), fill=pattern.FILL
    )
    b_img = warp.warp_nearest(
        img, stretch_then_rotate, out_shape=(18, 9), fill=pattern.FILL
    )
    print("  rotate then stretch      stretch then rotate")
    for la, lb in zip(pattern.as_text(a_img).split("\n"),
                      pattern.as_text(b_img).split("\n")):
        print(f"    {la}                {lb}")
    print()
    print(f"  identical images?  {np.array_equal(a_img, b_img)}")
    print("  Two different pictures from the same two operations. Matrix")
    print("  multiplication does not commute, and neither does the darkroom.")
    print()
    assert not warp.matrices_close(rotate_then_stretch, stretch_then_rotate)
    assert not np.array_equal(a_img, b_img)

    # ----------------------------------------------------------------------
    print("One matrix, one resample: the whole argument for composing")
    print("=" * 60)
    print("  A rotation, then a shear, then a scale. Two ways to get there.")
    print()
    steps = [
        ("rotate 30 degrees about the centre",
         warp.about_centre(warp.rotation(math.radians(30)), width, height)),
        ("shear x by 0.5", warp.shear_x(0.5)),
        ("scale by 1.5", warp.about_centre(warp.scaling(1.5, 1.5), width, height)),
    ]
    out_shape = (18, 18)

    sequential = img
    for _, matrix in steps:
        sequential = warp.warp_nearest(
            sequential, matrix, out_shape=out_shape, fill=pattern.FILL
        )

    single = warp.identity()
    for _, matrix in steps:
        single = warp.compose(matrix, single)
    at_once = warp.warp_nearest(img, single, out_shape=out_shape, fill=pattern.FILL)

    differing = int((sequential != at_once).sum())
    print(f"  three resampling passes vs one: {differing} of {sequential.size} "
          f"pixels differ ({100.0 * differing / sequential.size:.1f}%)")
    print()
    print("  three passes            one composed matrix")
    for la, lb in zip(pattern.as_text(sequential).split("\n"),
                      pattern.as_text(at_once).split("\n")):
        print(f"    {la}      {lb}")
    print()
    print("  The composed version is the correct one. Each intermediate resample")
    print("  in the three-pass version quantised the picture to whole pixels and")
    print("  threw the remainder away, and the next pass had no way to know. The")
    print("  matrices cost nine multiplications each to combine; the pixels cost")
    print("  a full pass over the image. Composing is both more accurate and")
    print("  cheaper, which is a rare combination and worth taking.")
    print()
    assert differing > 0
    assert single[2] == [0.0, 0.0, 1.0]

    print("  Composition also keeps the determinant honest -- the area factor of")
    print("  the whole is the product of the parts, exactly as on Day 102:")
    product = 1.0
    for name, matrix in steps:
        d = warp.determinant(matrix)
        product *= d
        print(f"    {name:<38} det {d:.6f}")
    print(f"    {'composed':<38} det {warp.determinant(single):.6f}")
    print(f"    {'product of the three':<38}     {product:.6f}")
    print()
    assert abs(warp.determinant(single) - product) <= 1e-12

    # ----------------------------------------------------------------------
    print("The inverse, and when there is not one")
    print("=" * 60)
    combo = warp.compose(warp.shear_x(0.5), warp.about_centre(
        warp.rotation(math.radians(37)), width, height))
    inverse = warp.invert(combo)
    identity_check = warp.compose(inverse, combo)
    print(f"  M . M^-1 is the identity within {warp.TOL:g}?  "
          f"{warp.matrices_close(identity_check, warp.identity())}")
    print(f"  det(M) = {warp.determinant(combo):.12f}, "
          f"det(M^-1) = {warp.determinant(inverse):.12f}, "
          f"product = {warp.determinant(combo) * warp.determinant(inverse):.12f}")
    print()
    assert warp.matrices_close(identity_check, warp.identity())

    print("  A transformation that flattens the image onto a line has")
    print("  determinant 0 and no inverse -- and because inverse mapping needs")
    print("  the inverse, such a transformation cannot be applied at all:")
    collapse = [[1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
    print(f"    determinant: {warp.determinant(collapse)}")
    try:
        warp.warp_nearest(img, collapse, fill=pattern.FILL)
    except warp.SingularTransform as exc:
        print(f"    warp_nearest raised {type(exc).__name__}")
        print(f"    message: {exc}")
        raised = type(exc).__name__
    else:
        raised = "NOTHING"
    print()
    assert raised == "SingularTransform"
    assert issubclass(warp.SingularTransform, ValueError)
    print("  SingularTransform is a ValueError, the same relationship")
    print("  numpy.linalg.LinAlgError has, so an existing `except ValueError`")
    print("  keeps working. Day 102 established that; nothing here changes it.")
    print()

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


if __name__ == "__main__":
    main()

"""06 — Twenty lines of ours against a mature library, on the same input.

Run from the examples directory:

    ../.venv/bin/python3 06_against_pillow.py

This is the day's strongest artifact and the reason the from-scratch code was
written with plain lists. If `warp.py` had been built out of NumPy helpers,
agreeing with a library would prove nothing. It was not, so agreement is
evidence.

It also settles the question Day 102 deliberately left open: what exactly is
Pillow's sampling convention, and why does a shear coefficient appear to move
row 0 when the mathematics says row 0 has y = 0?
"""

import math
import os
import random
import tempfile

import numpy as np
from PIL import Image

import pattern
import warp

SCRIPT = "06_against_pillow.py"

FILL = pattern.FILL


def pillow_affine(array, coefficients, out_shape=None, resample=None):
    """Run Pillow's own affine transform on the same array and coefficients."""
    height, width = array.shape
    out_h, out_w = out_shape or (height, width)
    image = Image.fromarray(array, mode="L")
    result = image.transform(
        (out_w, out_h),
        Image.Transform.AFFINE,
        coefficients,
        resample=resample or Image.Resampling.NEAREST,
        fillcolor=FILL,
    )
    return np.asarray(result)


def main():
    import PIL

    img = pattern.make_pattern()
    height, width = img.shape

    print(f"Pillow {PIL.__version__}, NumPy {np.__version__}")
    print()

    # ----------------------------------------------------------------------
    print("1. The convention: Pillow's coefficients run OUTPUT to INPUT")
    print("=" * 66)
    print("  Image.transform(size, AFFINE, (a, b, c, d, e, f)) means")
    print()
    print("      input_x  =  a * output_x  +  b * output_y  +  c")
    print("      input_y  =  d * output_x  +  e * output_y  +  f")
    print()
    print("  -- the inverse of the effect you see. Day 102 confirmed the")
    print("  DIRECTION by experiment. Here it is again, in one line, because a")
    print("  convention you have not checked today is a convention you are")
    print("  guessing at.")
    print()

    probe = np.zeros((1, 8), dtype=np.uint8)
    probe[0, 3] = 255
    shifted = pillow_affine(probe, (1, 0, 1, 0, 1, 0))
    before = int(np.argmax(probe[0]))
    after = int(np.argmax(shifted[0]))
    print(f"    a single bright pixel at input x = {before}")
    print(f"    coefficients (1, 0, 1, 0, 1, 0), so c = +1")
    print(f"    the bright pixel comes out at x = {after}")
    print(f"    the content moved LEFT by {before - after} when c said +1.")
    print()
    print("  That is the output-to-input direction, confirmed. If you want the")
    print("  picture to move right, you pass a NEGATIVE c -- or, better, you")
    print("  build the matrix you mean and let `to_pillow_coefficients` invert")
    print("  it for you, which is what this lab does everywhere below.")
    print()
    assert before == 3 and after == 2

    coeffs = warp.to_pillow_coefficients(warp.translation(1.0, 0.0))
    print(f"    to_pillow_coefficients(translation(1, 0)) = "
          f"{tuple(round(v, 12) for v in coeffs)}")
    print("    -- the c is -1, because the coefficients are read off the")
    print("    INVERSE of the matrix you asked for.")
    print()
    assert tuple(round(v, 12) for v in coeffs) == (1.0, 0.0, -1.0, 0.0, 1.0, 0.0)

    # ----------------------------------------------------------------------
    print("2. The open question from Day 102: where is the sample taken?")
    print("=" * 66)
    print("  Two candidate rules. Both agree on integer translations, which is")
    print("  why Day 102 could not tell them apart and said so rather than")
    print("  guessing:")
    print()
    print("    A (pixel centres):  source = floor(a*(x + 0.5) + b*(y + 0.5) + c)")
    print("    B (integer corners): source = floor(a*x + b*y + c + 0.5)")
    print()
    print("  A scale factor separates them in one measurement. Take the row")
    print("  0, 10, 20, ..., 70 and halve the image with a = 2:")
    print()

    row = (np.arange(8, dtype=np.uint8) * 10).reshape(1, 8)
    observed = [int(v) for v in pillow_affine(row, (2, 0, 0, 0, 1, 0))[0]]

    def predict(rule):
        out = []
        for x in range(8):
            i = rule(x)
            out.append(int(row[0, i]) if 0 <= i < 8 else FILL)
        return out

    model_a = predict(lambda x: math.floor(2 * (x + 0.5)))
    model_b = predict(lambda x: math.floor(2 * x + 0.5))
    print(f"    input             {[int(v) for v in row[0]]}")
    print(f"    Pillow observed   {observed}")
    print(f"    rule A predicts   {model_a}")
    print(f"    rule B predicts   {model_b}")
    print()
    print(f"    matches rule A?  {observed == model_a}")
    print(f"    matches rule B?  {observed == model_b}")
    print()
    assert observed == model_a
    assert observed != model_b

    print("  Rule A. Pillow evaluates the affine at the output pixel's CENTRE,")
    print("  (x + 0.5, y + 0.5), and takes the input pixel whose unit square")
    print("  contains the result. That is the answer Day 102 deferred, and it")
    print("  explains the shear puzzle exactly:")
    print()

    strip = np.zeros((3, 9), dtype=np.uint8)
    strip[:, 4] = 255
    sheared = pillow_affine(strip, (1, 2, 0, 0, 1, 0))
    print("    a vertical line at x = 4, with b = 2 (a shear in the")
    print("    output-to-input direction):")
    for y in range(3):
        found = np.flatnonzero(sheared[y] == 255)
        where = int(found[0]) if found.size else None
        predicted = math.floor(0.5 + 2 * (y + 0.5))
        print(f"      row {y}: line now at x = {where!s:<4} "
              f"shift predicted by rule A: {predicted}")
    print()
    print("    Row 0 MOVED, by one pixel, even though the shear term is")
    print("    multiplied by y and row 0 is 'at y = 0'. It is not at y = 0.")
    print("    Its centre is at y = 0.5, and 2 * 0.5 = 1. There is nothing")
    print("    mysterious left in it.")
    print()
    assert int(np.flatnonzero(sheared[0] == 255)[0]) == 3
    assert int(np.flatnonzero(sheared[1] == 255)[0]) == 1

    print("  `warp.py` uses the same rule -- see SAMPLE_OFFSET = 0.5 -- which is")
    print("  why the comparison below can be exact rather than approximate.")
    print()
    assert warp.SAMPLE_OFFSET == 0.5

    # ----------------------------------------------------------------------
    print("3. Ours against theirs, on 510 affine transformations")
    print("=" * 66)
    print("  500 random rotate-scale-shear-translate combinations plus 10")
    print("  deliberate edge cases, each handed to both implementations as the")
    print("  identical six numbers. Nearest-neighbour, same fill colour, same")
    print("  output size.")
    print()

    rng = random.Random(105)
    cases = []
    for _ in range(500):
        theta = rng.uniform(-math.pi, math.pi)
        scale = rng.uniform(0.4, 2.5)
        skew = rng.uniform(-2.5, 2.5)
        cos_t, sin_t = math.cos(theta), math.sin(theta)
        cases.append((
            scale * cos_t,
            scale * (cos_t * skew - sin_t),
            rng.uniform(-6, 6),
            scale * sin_t,
            scale * (sin_t * skew + cos_t),
            rng.uniform(-6, 6),
        ))
    edge_cases = [
        (1, 0, 0, 0, 1, 0),        # identity
        (1, 0, 1, 0, 1, 0),        # whole-pixel translation
        (1, 0, 0.5, 0, 1, 0),      # half-pixel translation: floor on a boundary
        (1, 2, 0, 0, 1, 0),        # the shear that moves row 0
        (2, 0, 0, 0, 2, 0),        # exact halving
        (0.5, 0, 0, 0, 0.5, 0),    # exact doubling
        (0, -1, 9, 1, 0, 0),       # a quarter turn
        (-1, 0, 9, 0, -1, 9),      # a half turn
        (1, 0, -3, 0, 1, -3),      # translation clean off the edge
        (1, 0.5, 0, 0, 1, 0),      # a gentle shear
    ]
    cases.extend(edge_cases)

    mismatches = 0
    worst_pixels = 0
    for coefficients in cases:
        inverse = warp.coefficients_to_matrix(coefficients)
        mine = warp.warp_nearest_with_inverse(img, inverse, fill=FILL)
        theirs = pillow_affine(img, coefficients)
        differing = int((mine != theirs).sum())
        worst_pixels = max(worst_pixels, differing)
        if differing:
            mismatches += 1

    print(f"    transformations compared:            {len(cases)}")
    print(f"    transformations matching EXACTLY:    {len(cases) - mismatches}")
    print(f"    worst case, pixels differing:        {worst_pixels}")
    print(f"    stated tolerance for this comparison: {pattern.PIXEL_TOL} "
          f"(exact equality)")
    print()
    assert mismatches == 0, mismatches
    assert worst_pixels == pattern.PIXEL_TOL

    print("  Every pixel of every one of them. Not 'close enough' -- identical.")
    print("  Twenty-odd lines of arithmetic in `warp_nearest_with_inverse` and a")
    print("  library maintained since 2010 produce byte-for-byte the same array.")
    print()

    # ----------------------------------------------------------------------
    print("4. Where the agreement DOES break, and why")
    print("=" * 66)
    print("  It would be easy to stop at the line above. It would also be")
    print("  misleading. Sweep every whole-degree rotation about the centre --")
    print("  360 transformations chosen to be nothing like random -- and the")
    print("  picture changes:")
    print()

    disagreeing = []
    worst_rot = 0
    furthest_from_boundary = 0.0
    for degrees in range(360):
        matrix = warp.about_centre(warp.rotation(math.radians(degrees)), width, height)
        coefficients = warp.to_pillow_coefficients(matrix)
        mine = warp.warp_nearest(img, matrix, fill=FILL)
        theirs = pillow_affine(img, coefficients)
        wrong = np.argwhere(mine != theirs)
        if len(wrong):
            disagreeing.append((degrees, len(wrong)))
            worst_rot = max(worst_rot, len(wrong))
        a, b, c, d, e, f = coefficients
        for oy, ox in wrong:
            xs = a * (ox + 0.5) + b * (oy + 0.5) + c
            ys = d * (ox + 0.5) + e * (oy + 0.5) + f
            gap = min(abs(xs - round(xs)), abs(ys - round(ys)))
            furthest_from_boundary = max(furthest_from_boundary, gap)

    print(f"    rotations compared:                 360")
    print(f"    identical, pixel for pixel:         {360 - len(disagreeing)}")
    print(f"    disagreeing in at least one pixel:  {len(disagreeing)}")
    print(f"    worst case, pixels differing:       {worst_rot} of 81")
    print(f"    the angles: {[deg for deg, _ in disagreeing]}")
    print()
    assert len(disagreeing) == 8, disagreeing
    assert worst_rot == 2, worst_rot

    print("  Eight angles out of 360, never more than 2 pixels out of 81. Now the")
    print("  useful part -- every single disagreeing sample landed within")
    print(f"  {furthest_from_boundary:.3e} of a pixel BOUNDARY:")
    print()
    assert furthest_from_boundary < 1e-9, furthest_from_boundary

    matrix = warp.about_centre(warp.rotation(math.radians(30)), width, height)
    coefficients = warp.to_pillow_coefficients(matrix)
    mine = warp.warp_nearest(img, matrix, fill=FILL)
    theirs = pillow_affine(img, coefficients)
    wrong = np.argwhere(mine != theirs)
    print("  30 degrees about the centre, the smallest failing case.")
    print(f"  coefficients passed to Pillow: "
          f"({', '.join(f'{v!r}' for v in coefficients)})")
    print()
    print("  ours                    Pillow")
    for a_line, b_line in zip(pattern.as_text(mine).split("\n"),
                              pattern.as_text(theirs).split("\n")):
        print(f"    {a_line}                {b_line}")
    print()
    a, b, c, d, e, f = coefficients
    for oy, ox in wrong:
        xs = a * (ox + 0.5) + b * (oy + 0.5) + c
        ys = d * (ox + 0.5) + e * (oy + 0.5) + f
        print(f"    output pixel (row {oy}, column {ox}): ours {int(mine[oy, ox])}, "
              f"Pillow {int(theirs[oy, ox])}")
        print(f"      our source y = {float(ys)!r}")
        print(f"      floor of that = {math.floor(ys)}; Pillow took row "
              f"{math.floor(ys) + 1}")
    print()
    assert len(wrong) == 1, len(wrong)
    assert int(wrong[0][0]) == 4 and int(wrong[0][1]) == 3

    print("  The source row is 4.999999999999999. The exact answer is 5. Ours")
    print("  floors to 4, Pillow's arithmetic reaches 5.0 or a hair above and")
    print("  floors to 5. Neither is wrong: the true sample sits exactly on the")
    print("  boundary between two pixels, and which one you get is decided by the")
    print("  ORDER the floating-point additions happen in. Pillow's C loop walks")
    print("  along each output row accumulating the source coordinate step by")
    print("  step; our Python evaluates the whole expression per pixel. Same")
    print("  formula, different rounding in the last bit.")
    print()
    print("  Look at WHICH angles failed: 30, 60, 120, 150, 210, 240, 300, 330.")
    print("  Every one of them is a 'nice' angle whose sine or cosine is exactly")
    print("  0.5 or exactly half the square root of 3. Nice angles are precisely")
    print("  the ones that put samples on boundaries. The angles nobody would")
    print("  choose for a test -- 37 degrees, 113 degrees -- all agreed. That is")
    print("  the opposite of the usual intuition and it is worth remembering:")
    print("  round numbers are where floating-point ties live.")
    print()
    print("  This is the real shape of the result, and it is more useful than")
    print("  'they always agree' would have been:")
    print("    * for a transformation whose numbers are not special, the two")
    print("      implementations are byte-for-byte identical -- 510 of 510;")
    print("    * for transformations that place samples exactly on pixel")
    print("      boundaries, they can differ by a pixel, and 8 of the 360 whole-")
    print("      degree rotations do;")
    print("    * the disagreement is never larger than the rounding step, and it")
    print("      is a property of floating point, not of either implementation.")
    print()
    print("  If you need bit-identical output across libraries, do not rely on")
    print("  ties breaking the same way. Use angles and offsets that keep samples")
    print("  away from boundaries, or accept a one-pixel tolerance and say so.")
    print()

    print("  And a non-square output, to check that the two agree about which")
    print("  way round a size tuple goes -- Pillow takes (width, height) and")
    print("  NumPy reports (height, width), which is one more place the two")
    print("  orderings can be swapped without any error being raised:")
    wide = warp.warp_nearest(img, warp.shear_x(1.0), out_shape=(9, 18), fill=FILL)
    wide_theirs = pillow_affine(
        img, warp.to_pillow_coefficients(warp.shear_x(1.0)), out_shape=(9, 18)
    )
    print(f"    ours   shape {wide.shape}")
    print(f"    Pillow shape {wide_theirs.shape}")
    print(f"    identical:   {np.array_equal(wide, wide_theirs)}")
    print()
    assert wide.shape == wide_theirs.shape == (9, 18)
    assert np.array_equal(wide, wide_theirs)

    # ----------------------------------------------------------------------
    print("5. A full turn: exact, and Pillow agrees it is exact")
    print("=" * 66)
    full = warp.about_centre(warp.rotation(2.0 * math.pi), width, height)
    ours_full = warp.warp_nearest(img, full, fill=FILL)
    theirs_full = pillow_affine(img, warp.to_pillow_coefficients(full))
    print(f"  rotation(2*pi) as ONE matrix:")
    print(f"    ours   differs from the original in {int((ours_full != img).sum())} pixels")
    print(f"    Pillow differs from the original in {int((theirs_full != img).sum())} pixels")
    print()
    assert np.array_equal(ours_full, img)
    assert np.array_equal(theirs_full, img)
    print("  Exactly zero, in both. Not 'within a tolerance' -- exact, and the")
    print("  reason is worth being precise about. The matrix is not exactly the")
    print("  identity: cos(2*pi) is 1.0 but sin(2*pi) is "
          f"{math.sin(2 * math.pi):.6e},")
    print("  not 0. The residual displacement is around 1e-15 of a pixel, and")
    print("  nearest-neighbour rounds to a whole pixel, so an error fifteen")
    print("  orders of magnitude below the rounding step cannot change the")
    print("  answer. The float error is real and it is absorbed.")
    print()
    assert math.sin(2 * math.pi) != 0.0
    assert not warp.matrices_close(full, warp.identity(), tol=0.0)
    assert warp.matrices_close(full, warp.identity(), tol=1e-12)

    print("  Twelve 30-degree passes is a different story, and script 03 measured")
    print("  it: 16 of 81 pixels lost. Same 360 degrees, resampled twelve times")
    print("  instead of once. Both implementations lose about the same amount,")
    print("  because the loss is in the METHOD and not in either of them:")
    ours_twelve, theirs_twelve = img, img
    step = warp.about_centre(warp.rotation(math.radians(30)), width, height)
    step_coeffs = warp.to_pillow_coefficients(step)
    for _ in range(12):
        ours_twelve = warp.warp_nearest(ours_twelve, step, fill=FILL)
        theirs_twelve = pillow_affine(theirs_twelve, step_coeffs)
    print(f"    ours   differs from the original in "
          f"{int((ours_twelve != img).sum())} pixels")
    print(f"    Pillow differs from the original in "
          f"{int((theirs_twelve != img).sum())} pixels")
    print(f"    ours and Pillow differ from each other in "
          f"{int((ours_twelve != theirs_twelve).sum())} pixels")
    print()
    assert int((ours_twelve != img).sum()) == 16
    assert int((theirs_twelve != img).sum()) == 17
    assert int((ours_twelve != theirs_twelve).sum()) == 3
    print("  16 against 17, and the two results differ from EACH OTHER in 3")
    print("  pixels. That is section 4 compounding: 30 degrees is one of the")
    print("  eight tie-prone angles, so each of the twelve passes can take a")
    print("  different branch, and twelve passes of a one-pixel difference is")
    print("  three pixels apart at the end rather than one. A single pass agreed")
    print("  exactly. Repeated resampling does not just lose information -- it")
    print("  amplifies the disagreements too. One more reason to compose.")
    print()

    # ----------------------------------------------------------------------
    print("6. Bilinear: where the agreement stops, stated plainly")
    print("=" * 66)
    print("  Nearest-neighbour picks the closest pixel. Bilinear averages the")
    print("  four surrounding pixels, weighted by distance -- which is what you")
    print("  want when the inverse-mapped position lands between pixels, because")
    print("  it usually does.")
    print()
    print("  The visible difference: nearest-neighbour gives hard, stair-stepped")
    print("  edges and keeps every value exactly as it was; bilinear gives smooth")
    print("  edges and INVENTS intermediate values that were not in the input.")
    print()

    small_shift = warp.translation(0.5, 0.0)
    nn = warp.warp_nearest(img, small_shift, fill=FILL)
    bl = warp.warp_bilinear_with_inverse(
        img, warp.invert(small_shift), fill=float(FILL)
    )
    print(f"    distinct values, input:              {sorted(np.unique(img).tolist())}")
    print(f"    distinct values, nearest-neighbour:  {sorted(np.unique(nn).tolist())}")
    print(f"    distinct values, bilinear:           "
          f"{len(np.unique(np.round(bl, 6)))} different levels")
    print()
    assert len(np.unique(nn)) <= 4
    assert len(np.unique(np.round(bl, 6))) > len(np.unique(nn))

    print("  Now the honest part. Our bilinear does NOT reproduce Pillow's")
    print("  bilinear pixel-for-pixel, and the lab says so rather than quietly")
    print("  loosening a tolerance until it passes.")
    print()

    print("  The split turns out to be clean, and it is worth stating exactly")
    print("  rather than as 'roughly agrees'. Separate the output pixels into")
    print("  those whose four contributing input pixels are ALL inside the")
    print("  image, and those where at least one contributor lies outside it.")
    print()

    bilinear_cases = [
        ("translate (0.25, 0.25)", warp.translation(0.25, 0.25)),
        ("rotate 30 about the centre",
         warp.about_centre(warp.rotation(math.radians(30)), width, height)),
        ("rotate 17 about the centre",
         warp.about_centre(warp.rotation(math.radians(17)), width, height)),
        ("scale 1.5 about the centre",
         warp.about_centre(warp.scaling(1.5, 1.5), width, height)),
        ("shear x by 0.4", warp.shear_x(0.4)),
    ]

    worst_inside = 0.0
    worst_anywhere = 0.0
    print(f"    {'transformation':<28}{'all 4 inside':>14}{'anywhere':>12}")
    for name, matrix in bilinear_cases:
        inverse = warp.invert(matrix)
        ours_bl = warp.warp_bilinear_with_inverse(img, inverse, fill=0.0)
        theirs_bl = pillow_affine(
            img,
            warp.to_pillow_coefficients(matrix),
            resample=Image.Resampling.BILINEAR,
        ).astype(float)
        difference = np.abs(ours_bl - theirs_bl)

        inside = np.zeros((height, width), dtype=bool)
        for oy in range(height):
            for ox in range(width):
                sx, sy = warp.apply_point(inverse, (ox + 0.5, oy + 0.5))
                x0 = math.floor(sx - warp.SAMPLE_OFFSET)
                y0 = math.floor(sy - warp.SAMPLE_OFFSET)
                inside[oy, ox] = (
                    0 <= x0 and x0 + 1 < width and 0 <= y0 and y0 + 1 < height
                )

        in_max = float(difference[inside].max()) if inside.any() else 0.0
        any_max = float(difference.max())
        worst_inside = max(worst_inside, in_max)
        worst_anywhere = max(worst_anywhere, any_max)
        print(f"    {name:<28}{in_max:>14.3f}{any_max:>12.3f}")
    print()
    print(f"    worst difference where all four contributors are inside: "
          f"{worst_inside:.3f}")
    print(f"    worst difference anywhere:                               "
          f"{worst_anywhere:.3f}")
    print()
    assert worst_inside <= 1.0, worst_inside
    assert worst_anywhere > 100.0, worst_anywhere

    print("  So the claim this lab makes is precise, and it is a better claim")
    print("  than 'they agree' would have been:")
    print()
    print("    * NEAREST-NEIGHBOUR: identical to Pillow, 510 of 510 random and")
    print("      edge cases, zero differing pixels; and identical on 352 of the")
    print("      360 whole-degree rotations, the other 8 differing by at most 2")
    print("      pixels at floating-point ties.")
    print()
    print("    * BILINEAR: wherever all four contributing pixels are inside the")
    print(f"      image, the two agree to within {worst_inside:.0f} grey level --")
    print("      which is exactly the rounding of a float average back into a")
    print("      byte, and cannot be improved on. Wherever a contributor lies")
    print("      OUTSIDE the image, they diverge by up to "
          f"{worst_anywhere:.0f} levels,")
    print("      because they extrapolate differently: ours averages the fill")
    print("      value in, Pillow does not.")
    print()
    print("  The border behaviour was measured, not assumed, and it was not")
    print("  chased further. Naming the boundary of what agrees is more useful")
    print("  than widening a tolerance until a test goes green.")
    print()

    # ----------------------------------------------------------------------
    print("7. Through a real file, and cleaned up afterwards")
    print("=" * 66)
    print("  Everything above happened in memory. One round trip through an")
    print("  actual PNG, to show that the file format is not where information")
    print("  is lost -- PNG is lossless, so the array survives byte for byte.")
    print()
    with tempfile.TemporaryDirectory() as tmp:
        path = os.path.join(tmp, "pattern.png")
        Image.fromarray(img, mode="L").save(path)
        size = os.path.getsize(path)
        reloaded = np.asarray(Image.open(path).convert("L"))
        print(f"    written, {size} bytes on disk for {img.size} pixels")
        print(f"    reloaded array identical to the original?  "
              f"{np.array_equal(reloaded, img)}")
        assert np.array_equal(reloaded, img)
        assert os.path.exists(path)
    print(f"    temporary directory removed?  {not os.path.exists(path)}")
    print()
    assert not os.path.exists(path)
    print("  The file lived in the operating system's temporary directory and is")
    print("  gone. This lab writes no image into its own tree, which is why there")
    print("  is nothing to commit and nothing to clean up by hand.")
    print()

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


if __name__ == "__main__":
    main()

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

Both `examples/` and `starter/` contain modules called `warp` and `pattern`,
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 `warp` 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
`warp` or `pattern` 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 ("warp", "pattern", "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 test image, generated in code rather than downloaded.

Nothing here loads a file. The pattern is built from arithmetic so that every
pixel value in this lab has a reason, and so that the lab needs no network and
ships no photograph.

The glyph is a capital F on a 9 by 9 grid. An F is the standard test shape in
image processing for one reason: it is asymmetric under every operation this
lab performs. A square survives a horizontal flip unchanged, so a square would
let a broken flip pass. An F does not survive anything -- flip it, rotate it,
transpose it, and you can see immediately which one happened.

Coordinates. This is the single most common source of confusion in the whole
subject, so it is stated once, here, and every function in `warp.py` obeys it:

    image[row, column]  ==  image[y, x]

A NumPy array is indexed rows first. The maths of Week 15 wrote points as
(x, y). Those are the SAME two numbers in the OPPOSITE order. Row is y, column
is x. And y grows DOWNWARD, because row 0 is the top of the picture, which is
the reverse of the graphs on Day 102 where y grew upward.
"""

import numpy as np

# Grey levels. Distinct, well separated, and none of them 0 or 255 by accident.
INK = 255  # the glyph itself
PAPER = 0  # the background
MARK = 96  # a single corner pixel, so a 180 degree turn is distinguishable
FILL = 32  # what lands in an output pixel whose source is off the image

SIZE = 9  # the pattern is SIZE by SIZE

# The exact cells of the F, written out rather than computed, so the values in
# the tests can be read off this list by eye.
TOP_BAR = [(0, c) for c in range(1, 7)]  # row 0, columns 1..6
MIDDLE_BAR = [(4, c) for c in range(1, 5)]  # row 4, columns 1..4
STEM = [(r, 1) for r in range(1, 9)] + [(r, 2) for r in range(1, 9)]

# Where the single corner mark goes: bottom-right, the corner the F never
# reaches, so it can never be confused with part of the glyph.
MARK_CELL = (8, 8)


def make_pattern():
    """Return the 9 by 9 greyscale test image as a uint8 array of shape (9, 9).

    Shape is (height, width) -- rows then columns -- which is (y, x).
    """
    img = np.full((SIZE, SIZE), PAPER, dtype=np.uint8)
    for row, col in TOP_BAR + MIDDLE_BAR + STEM:
        img[row, col] = INK
    img[MARK_CELL] = MARK
    return img


def make_colour_pattern():
    """Return the same glyph as a colour image of shape (9, 9, 3).

    Colour is three greyscale planes stacked on the last axis: red, green,
    blue. Each plane is a matrix in its own right, and every transformation in
    this lab acts on the coordinates, which the three planes share. That is
    why transforming a colour image is the same work as transforming a
    greyscale one, done three times.
    """
    grey = make_pattern()
    img = np.zeros((SIZE, SIZE, 3), dtype=np.uint8)
    img[:, :, 0] = grey  # red plane: the glyph
    img[:, :, 1] = np.fliplr(grey)  # green plane: the glyph mirrored
    img[:, :, 2] = MARK  # blue plane: flat, so a channel mix-up is obvious
    return img


def as_text(img, ink_char="#", paper_char=".", mark_char="o", fill_char="~"):
    """Render a greyscale array as ASCII so a transformation can be SEEN.

    Anything that is not one of the four known levels prints as `?`, which is
    how an interpolation artefact announces itself.
    """
    table = {INK: ink_char, PAPER: paper_char, MARK: mark_char, FILL: fill_char}
    return "\n".join(
        "".join(table.get(int(v), "?") for v in row) for row in np.asarray(img)
    )


def ink_cells(img):
    """Return the sorted (row, column) pairs whose value is INK.

    Comparing two images by their ink cells is exact -- these are integers, not
    floats -- which is why the rotation tests in this lab can assert equality
    rather than a tolerance.
    """
    rows, cols = np.nonzero(np.asarray(img) == INK)
    return sorted(zip(rows.tolist(), cols.tolist()))


# The values the tests assert against, written here once so that a change to
# the pattern cannot quietly change what "correct" means somewhere else.
EXPECTED_SHAPE = (9, 9)
EXPECTED_COLOUR_SHAPE = (9, 9, 3)
EXPECTED_INK_COUNT = len(set(TOP_BAR + MIDDLE_BAR + STEM))
EXPECTED_TEXT = "\n".join(
    [
        ".######..",
        ".##......",
        ".##......",
        ".##......",
        ".####....",
        ".##......",
        ".##......",
        ".##......",
        ".##.....o",
    ]
)

# Tolerances. Every float comparison in this lab names one of these.
TOL = 1e-12  # for matrix arithmetic done in floating point
PIXEL_TOL = 0  # for nearest-neighbour pixel values: they must match EXACTLY

"""The reference test suite: real pixels, real values, real disagreements.

Run from the LAB DIRECTORY:

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

Every float comparison states a tolerance. Every pixel comparison states
whether it is exact or within a stated number of grey levels, and the exact
ones say `==` on purpose, because that is the stronger claim.
"""

import math
import os
import random
import tempfile

import numpy as np
import pytest
from PIL import Image

import pattern
import warp

TOL = pattern.TOL
FILL = pattern.FILL


def pillow_affine(array, coefficients, out_shape=None, resample=None):
    height, width = array.shape
    out_h, out_w = out_shape or (height, width)
    image = Image.fromarray(array, mode="L")
    return np.asarray(
        image.transform(
            (out_w, out_h),
            Image.Transform.AFFINE,
            coefficients,
            resample=resample or Image.Resampling.NEAREST,
            fillcolor=FILL,
        )
    )


@pytest.fixture
def img():
    return pattern.make_pattern()


# -- The image itself --------------------------------------------------------


def test_pattern_has_the_documented_shape_and_type(img):
    assert img.shape == pattern.EXPECTED_SHAPE == (9, 9)
    assert img.dtype == np.uint8
    assert img.nbytes == 81


def test_pattern_pixel_values_are_exactly_as_written(img):
    assert pattern.as_text(img) == pattern.EXPECTED_TEXT
    assert int(img[0, 1]) == pattern.INK
    assert int(img[0, 0]) == pattern.PAPER
    assert int(img[8, 8]) == pattern.MARK
    assert int((img == pattern.INK).sum()) == pattern.EXPECTED_INK_COUNT == 24
    assert int((img == pattern.MARK).sum()) == 1


def test_row_and_column_are_not_interchangeable(img):
    # The ordering trap, asserted rather than described.
    assert int(img[4, 3]) == pattern.INK
    assert int(img[3, 4]) == pattern.PAPER
    assert img[4, 3] != img[3, 4]


def test_the_pattern_is_asymmetric_under_every_operation_tested(img):
    assert not np.array_equal(img, np.fliplr(img))
    assert not np.array_equal(img, np.flipud(img))
    assert not np.array_equal(img, img.T)
    assert not np.array_equal(img, np.rot90(img))


def test_colour_is_three_stacked_planes(img):
    colour = pattern.make_colour_pattern()
    assert colour.shape == (9, 9, 3)
    assert np.array_equal(colour[:, :, 0], img)
    assert np.array_equal(colour[:, :, 1], np.fliplr(img))
    assert int(colour[:, :, 2].min()) == int(colour[:, :, 2].max()) == pattern.MARK


def test_the_pattern_is_generated_not_loaded():
    # Two independent calls must agree, and neither may touch the filesystem.
    assert np.array_equal(pattern.make_pattern(), pattern.make_pattern())


# -- Matrices ---------------------------------------------------------------


def test_every_affine_matrix_has_the_bottom_row_0_0_1():
    for matrix in (
        warp.identity(),
        warp.translation(3.0, -2.0),
        warp.scaling(2.0, 0.5),
        warp.rotation(1.1),
        warp.shear_x(2.0),
        warp.shear_y(-0.5),
        warp.flip_horizontal(9),
        warp.flip_vertical(9),
        warp.about_centre(warp.rotation(0.3), 9, 9),
    ):
        assert matrix[2] == [0.0, 0.0, 1.0]


def test_translation_is_not_linear_but_every_2x2_part_is():
    # A linear map must fix the origin; translation does not.
    assert warp.apply_point(warp.translation(3.0, -2.0), (0.0, 0.0)) == (3.0, -2.0)
    for matrix in (warp.rotation(1.1), warp.scaling(4.0, 0.2), warp.shear_x(9.0)):
        assert warp.apply_point(matrix, (0.0, 0.0)) == (0.0, 0.0)


def test_translation_composes_and_inverts():
    shift = warp.translation(3.0, -2.0)
    assert warp.determinant(shift) == 1.0
    assert warp.matrices_close(warp.invert(shift), warp.translation(-3.0, 2.0), TOL)
    assert warp.matrices_close(
        warp.compose(warp.invert(shift), shift), warp.identity(), TOL
    )


def test_quarter_turn_matrix_is_exact_where_trigonometry_is_not():
    exact = warp.rotation_quarter_turns(1)
    trig = warp.rotation(math.pi / 2)
    assert exact == [[0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
    assert math.cos(math.pi / 2) != 0.0
    assert trig != exact
    assert warp.matrices_close(trig, exact, TOL)


def test_four_quarter_turns_return_to_the_identity():
    turn = warp.rotation_quarter_turns(1)
    assert warp.matrices_close(
        warp.compose(turn, turn, turn, turn), warp.identity(), TOL
    )


def test_composition_order_is_right_to_left():
    first = warp.translation(1.0, 0.0)
    second = warp.scaling(2.0, 2.0)
    combined = warp.compose(second, first)
    stepwise = warp.apply_point(second, warp.apply_point(first, (3.0, 5.0)))
    assert combined is not first
    for a, b in zip(warp.apply_point(combined, (3.0, 5.0)), stepwise):
        assert abs(a - b) <= TOL
    # And the other order is genuinely different.
    assert not warp.matrices_close(warp.compose(first, second), combined, TOL)


def test_determinant_of_a_composition_is_the_product():
    parts = [warp.scaling(1.5, 1.5), warp.shear_x(0.5), warp.rotation(0.7)]
    combined = warp.compose(*parts)
    product = 1.0
    for matrix in parts:
        product *= warp.determinant(matrix)
    assert abs(warp.determinant(combined) - product) <= TOL


def test_determinants_of_the_named_transformations():
    assert warp.determinant(warp.identity()) == 1.0
    assert warp.determinant(warp.translation(7.0, -3.0)) == 1.0
    assert warp.determinant(warp.shear_x(9.0)) == 1.0
    assert warp.determinant(warp.flip_horizontal(9)) == -1.0
    assert warp.determinant(warp.scaling(2.0, 3.0)) == 6.0


def test_inverse_round_trips_within_tolerance():
    matrix = warp.compose(
        warp.shear_x(0.5), warp.about_centre(warp.rotation(math.radians(37)), 9, 9)
    )
    assert warp.matrices_close(
        warp.compose(warp.invert(matrix), matrix), warp.identity(), TOL
    )
    assert warp.matrices_close(
        warp.compose(matrix, warp.invert(matrix)), warp.identity(), TOL
    )


def test_a_collapsing_transformation_cannot_be_inverted(img):
    collapse = [[1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
    assert warp.determinant(collapse) == 0.0
    with pytest.raises(warp.SingularTransform):
        warp.invert(collapse)
    # And therefore it cannot be applied at all, because inverse mapping needs
    # the inverse.
    with pytest.raises(warp.SingularTransform):
        warp.warp_nearest(img, collapse, fill=FILL)


def test_singular_transform_is_catchable_as_a_value_error():
    assert issubclass(warp.SingularTransform, ValueError)


def test_affine_maps_send_parallel_lines_to_parallel_lines():
    for matrix in (warp.rotation(0.7), warp.shear_x(2.0), warp.scaling(3.0, 0.5)):
        p0 = warp.apply_point(matrix, (0.0, 0.0))
        p1 = warp.apply_point(matrix, (1.0, 0.0))
        p2 = warp.apply_point(matrix, (0.0, 5.0))
        p3 = warp.apply_point(matrix, (1.0, 5.0))
        v1 = (p1[0] - p0[0], p1[1] - p0[1])
        v2 = (p3[0] - p2[0], p3[1] - p2[1])
        assert abs(v1[0] * v2[1] - v1[1] * v2[0]) <= TOL


# -- Forward mapping: the failure, asserted ---------------------------------


def test_forward_mapping_leaves_holes_on_a_rotation(img):
    matrix = warp.about_centre(warp.rotation(math.radians(30)), 9, 9)
    _, holes = warp.warp_forward(img, matrix, fill=FILL)
    assert int(holes.sum()) == 22
    assert holes.shape == (9, 9)


def test_forward_mapping_holes_are_arithmetically_unavoidable_when_enlarging(img):
    _, holes = warp.warp_forward(
        img, warp.scaling(2.0, 2.0), out_shape=(18, 18), fill=FILL
    )
    # 81 input pixels cannot fill 324 output pixels, whatever the rounding does.
    assert int(holes.sum()) == 243
    assert int(holes.sum()) >= holes.size - img.size


def test_forward_mapping_overwrites_when_shrinking(img):
    _, holes = warp.warp_forward(
        img, warp.scaling(0.5, 0.5), out_shape=(5, 5), fill=FILL
    )
    assert int(holes.sum()) == 0  # no holes ...
    landings = {}
    for y in range(9):
        for x in range(9):
            fx, fy = warp.apply_point(warp.scaling(0.5, 0.5), (x + 0.5, y + 0.5))
            landings.setdefault((math.floor(fy), math.floor(fx)), []).append((y, x))
    # ... but 24 of the 25 output pixels were written more than once.
    assert sum(1 for v in landings.values() if len(v) > 1) == 24
    assert max(len(v) for v in landings.values()) == 4


def test_inverse_mapping_has_no_holes_only_clipping(img):
    matrix = warp.about_centre(warp.rotation(math.radians(30)), 9, 9)
    out = warp.warp_nearest(img, matrix, fill=FILL)
    back = warp.invert(matrix)
    filled = np.argwhere(out == FILL)
    assert len(filled) == 12
    for oy, ox in filled:
        sx, sy = warp.apply_point(back, (ox + 0.5, oy + 0.5))
        inside = 0 <= math.floor(sx) < 9 and 0 <= math.floor(sy) < 9
        assert not inside, (oy, ox)


# -- Inverse mapping: exact answers ------------------------------------------


@pytest.mark.parametrize("turns", [1, 2, 3, 4])
def test_quarter_turns_match_numpy_rot90_exactly(img, turns):
    matrix = warp.about_centre(warp.rotation_quarter_turns(turns), 9, 9)
    out = warp.warp_nearest(img, matrix, fill=FILL)
    assert np.array_equal(out, np.rot90(img, -turns))
    assert int((out == FILL).sum()) == 0


def test_a_quarter_turn_moves_named_pixels_to_named_places(img):
    matrix = warp.about_centre(warp.rotation_quarter_turns(1), 9, 9)
    out = warp.warp_nearest(img, matrix, fill=FILL)
    # The corner mark: bottom-right to bottom-left.
    assert tuple(int(v) for v in np.argwhere(img == pattern.MARK)[0]) == (8, 8)
    assert tuple(int(v) for v in np.argwhere(out == pattern.MARK)[0]) == (8, 0)
    # The top bar, row 0 columns 1..6, becomes column 8 rows 1..6.
    bar = sorted(int(r) for r, c in np.argwhere(out == pattern.INK) if c == 8)
    assert bar == list(range(1, 7))
    assert len(pattern.ink_cells(out)) == len(pattern.ink_cells(img)) == 24


def test_the_trigonometric_quarter_turn_gives_identical_pixels(img):
    exact = warp.warp_nearest(
        img, warp.about_centre(warp.rotation_quarter_turns(1), 9, 9), fill=FILL
    )
    trig = warp.warp_nearest(
        img, warp.about_centre(warp.rotation(math.pi / 2), 9, 9), fill=FILL
    )
    assert np.array_equal(exact, trig)


def test_flips_match_numpy_exactly(img):
    assert np.array_equal(
        warp.warp_nearest(img, warp.flip_horizontal(9), fill=FILL), np.fliplr(img)
    )
    assert np.array_equal(
        warp.warp_nearest(img, warp.flip_vertical(9), fill=FILL), np.flipud(img)
    )


def test_flipping_twice_restores_the_original(img):
    once = warp.warp_nearest(img, warp.flip_horizontal(9), fill=FILL)
    twice = warp.warp_nearest(once, warp.flip_horizontal(9), fill=FILL)
    assert np.array_equal(twice, img)


def test_doubling_is_exact_pixel_replication(img):
    out = warp.warp_nearest(img, warp.scaling(2.0, 2.0), out_shape=(18, 18), fill=FILL)
    assert np.array_equal(out, np.kron(img, np.ones((2, 2), dtype=np.uint8)))
    assert int((out == FILL).sum()) == 0
    # Four times the pixels, exactly the same information.
    assert len(np.unique(out)) == len(np.unique(img)) == 3


def test_halving_is_exactly_a_strided_slice(img):
    out = warp.warp_nearest(img, warp.scaling(0.5, 0.5), out_shape=(4, 4), fill=FILL)
    assert np.array_equal(out, img[1::2, 1::2])
    # The single-pixel corner mark is one of the pixels thrown away.
    assert int((out == pattern.MARK).sum()) == 0


def test_translation_matches_a_slice_and_pad(img):
    out = warp.warp_nearest(img, warp.translation(2.0, 1.0), fill=FILL)
    reference = np.full_like(img, FILL)
    reference[1:, 2:] = img[:-1, :-2]
    assert np.array_equal(out, reference)
    assert int((out == FILL).sum()) == 2 * 9 + 9 - 2


@pytest.mark.parametrize(
    "k,out_width,expected_row0_shift", [(0.5, 14, 0), (1.0, 18, 0), (2.0, 27, 1)]
)
def test_the_half_pixel_offset_decides_whether_row_0_moves(
    img, k, out_width, expected_row0_shift
):
    """The question Day 102 deferred, asserted as a number.

    Row 0's output pixels are sampled at y = 0.5, not y = 0, so the shear term
    contributes k * 0.5 even in the top row.
    """
    out = warp.warp_nearest(
        img, warp.shear_x(k), out_shape=(9, out_width), fill=FILL
    )
    first_ink = int(np.argmax(out[0] == pattern.INK))
    assert first_ink - 1 == expected_row0_shift
    assert first_ink - 1 == -math.floor(warp.SAMPLE_OFFSET - k * warp.SAMPLE_OFFSET)


def test_a_shear_preserves_the_ink_count(img):
    out = warp.warp_nearest(img, warp.shear_x(1.0), out_shape=(9, 18), fill=FILL)
    assert warp.determinant(warp.shear_x(1.0)) == 1.0
    assert int((out == pattern.INK).sum()) == int((img == pattern.INK).sum()) == 24


def test_colour_is_transformed_plane_by_plane_with_the_same_matrix():
    colour = pattern.make_colour_pattern()
    matrix = warp.about_centre(warp.rotation_quarter_turns(1), 9, 9)
    out = warp.warp_colour(colour, matrix, fill=FILL)
    assert out.shape == (9, 9, 3)
    for channel in range(3):
        assert np.array_equal(out[:, :, channel], np.rot90(colour[:, :, channel], -1))


def test_warp_rejects_a_colour_array_with_a_helpful_message():
    colour = pattern.make_colour_pattern()
    with pytest.raises(ValueError) as caught:
        warp.warp_nearest(colour, warp.identity(), fill=FILL)
    assert "(height, width)" in str(caught.value)


# -- Composing versus repeating ---------------------------------------------


def test_a_full_turn_as_one_matrix_is_pixel_exact(img):
    out = warp.warp_nearest(
        img, warp.about_centre(warp.rotation(2.0 * math.pi), 9, 9), fill=FILL
    )
    assert np.array_equal(out, img)


def test_twelve_separate_thirty_degree_turns_are_not(img):
    step = warp.about_centre(warp.rotation(math.radians(30)), 9, 9)
    out = img
    for _ in range(12):
        out = warp.warp_nearest(out, step, fill=FILL)
    assert int((out != img).sum()) == 16
    assert not np.array_equal(out, img)


def test_the_same_twelve_turns_composed_into_one_matrix_are_exact(img):
    step = warp.about_centre(warp.rotation(math.radians(30)), 9, 9)
    combined = warp.identity()
    for _ in range(12):
        combined = warp.compose(step, combined)
    assert warp.matrices_close(combined, warp.identity(), TOL)
    assert np.array_equal(warp.warp_nearest(img, combined, fill=FILL), img)


def test_one_composed_matrix_agrees_with_three_separate_steps_on_every_point():
    steps = [
        warp.about_centre(warp.rotation(math.radians(30)), 9, 9),
        warp.shear_x(0.5),
        warp.about_centre(warp.scaling(1.5, 1.5), 9, 9),
    ]
    combined = warp.identity()
    for matrix in steps:
        combined = warp.compose(matrix, combined)
    for y in range(9):
        for x in range(9):
            point = (x + 0.5, y + 0.5)
            stepwise = point
            for matrix in steps:
                stepwise = warp.apply_point(matrix, stepwise)
            at_once = warp.apply_point(combined, point)
            for a, b in zip(stepwise, at_once):
                assert abs(a - b) <= TOL


def test_an_integer_shear_round_trip_is_lossy_across_two_passes(img):
    """A boundary case that is asserted rather than hidden.

    k = 1.0 puts every sample exactly on a pixel boundary, where floor must
    make an arbitrary choice; k = 0.5 keeps samples well inside a pixel.
    """
    half = warp.warp_nearest(img, warp.shear_x(0.5), out_shape=(9, 14), fill=FILL)
    half_back = warp.warp_nearest(half, warp.shear_x(-0.5), out_shape=(9, 9), fill=FILL)
    assert np.array_equal(half_back, img)

    unit = warp.warp_nearest(img, warp.shear_x(1.0), out_shape=(9, 18), fill=FILL)
    unit_back = warp.warp_nearest(unit, warp.shear_x(-1.0), out_shape=(9, 9), fill=FILL)
    assert int((unit_back != img).sum()) == 28

    # As ONE composed matrix it is exact for both.
    for k in (0.5, 1.0):
        combined = warp.compose(warp.shear_x(-k), warp.shear_x(k))
        assert warp.matrices_close(combined, warp.identity(), TOL)
        assert np.array_equal(warp.warp_nearest(img, combined, fill=FILL), img)


# -- Against Pillow ----------------------------------------------------------


def test_pillow_coefficients_run_output_to_input():
    probe = np.zeros((1, 8), dtype=np.uint8)
    probe[0, 3] = 255
    out = pillow_affine(probe, (1, 0, 1, 0, 1, 0))
    # A positive c moved the content LEFT: the coefficients are the inverse.
    assert int(np.argmax(out[0])) == 2


def test_to_pillow_coefficients_inverts_the_matrix():
    coefficients = warp.to_pillow_coefficients(warp.translation(1.0, 0.0))
    assert tuple(round(v, 12) for v in coefficients) == (1.0, 0.0, -1.0, 0.0, 1.0, 0.0)
    assert warp.coefficients_to_matrix(coefficients)[2] == [0.0, 0.0, 1.0]


def test_pillow_samples_at_pixel_centres_not_at_integer_corners():
    """Settles the sampling question Day 102 deferred, by measurement."""
    row = (np.arange(8, dtype=np.uint8) * 10).reshape(1, 8)
    observed = [int(v) for v in pillow_affine(row, (2, 0, 0, 0, 1, 0))[0]]

    def predict(rule):
        return [
            int(row[0, rule(x)]) if 0 <= rule(x) < 8 else FILL for x in range(8)
        ]

    centres = predict(lambda x: math.floor(2 * (x + 0.5)))
    corners = predict(lambda x: math.floor(2 * x + 0.5))
    assert observed == centres
    assert observed != corners
    assert warp.SAMPLE_OFFSET == 0.5


def test_pillows_shear_moves_row_zero_and_the_offset_explains_it():
    strip = np.zeros((3, 9), dtype=np.uint8)
    strip[:, 4] = 255
    out = pillow_affine(strip, (1, 2, 0, 0, 1, 0))
    for y, expected_shift in ((0, 1), (1, 3)):
        found = np.flatnonzero(out[y] == 255)
        assert found.size == 1
        assert int(found[0]) == 4 - expected_shift
        assert expected_shift == math.floor(0.5 + 2 * (y + 0.5))


def test_ours_and_pillow_agree_exactly_on_510_affine_transformations(img):
    rng = random.Random(105)
    cases = []
    for _ in range(500):
        theta = rng.uniform(-math.pi, math.pi)
        scale = rng.uniform(0.4, 2.5)
        skew = rng.uniform(-2.5, 2.5)
        cos_t, sin_t = math.cos(theta), math.sin(theta)
        cases.append((
            scale * cos_t,
            scale * (cos_t * skew - sin_t),
            rng.uniform(-6, 6),
            scale * sin_t,
            scale * (sin_t * skew + cos_t),
            rng.uniform(-6, 6),
        ))
    cases += [
        (1, 0, 0, 0, 1, 0),
        (1, 0, 1, 0, 1, 0),
        (1, 0, 0.5, 0, 1, 0),
        (1, 2, 0, 0, 1, 0),
        (2, 0, 0, 0, 2, 0),
        (0.5, 0, 0, 0, 0.5, 0),
        (0, -1, 9, 1, 0, 0),
        (-1, 0, 9, 0, -1, 9),
        (1, 0, -3, 0, 1, -3),
        (1, 0.5, 0, 0, 1, 0),
    ]
    assert len(cases) == 510
    for coefficients in cases:
        mine = warp.warp_nearest_with_inverse(
            img, warp.coefficients_to_matrix(coefficients), fill=FILL
        )
        theirs = pillow_affine(img, coefficients)
        assert int((mine != theirs).sum()) == pattern.PIXEL_TOL == 0


def test_where_ours_and_pillow_disagree_the_sample_is_on_a_pixel_boundary(img):
    """The honest half of the comparison.

    Eight of the 360 whole-degree rotations differ, by at most 2 pixels of 81,
    and every disagreeing sample lands within one ulp of a pixel boundary.
    """
    disagreeing = []
    furthest = 0.0
    for degrees in range(360):
        matrix = warp.about_centre(warp.rotation(math.radians(degrees)), 9, 9)
        coefficients = warp.to_pillow_coefficients(matrix)
        mine = warp.warp_nearest(img, matrix, fill=FILL)
        theirs = pillow_affine(img, coefficients)
        wrong = np.argwhere(mine != theirs)
        if len(wrong):
            disagreeing.append((degrees, len(wrong)))
        a, b, c, d, e, f = coefficients
        for oy, ox in wrong:
            xs = a * (ox + 0.5) + b * (oy + 0.5) + c
            ys = d * (ox + 0.5) + e * (oy + 0.5) + f
            furthest = max(
                furthest, min(abs(xs - round(xs)), abs(ys - round(ys)))
            )

    assert [deg for deg, _ in disagreeing] == [30, 60, 120, 150, 210, 240, 300, 330]
    assert max(count for _, count in disagreeing) == 2
    assert furthest < 1e-9


def test_the_thirty_degree_disagreement_is_one_pixel_at_a_tie(img):
    matrix = warp.about_centre(warp.rotation(math.radians(30)), 9, 9)
    coefficients = warp.to_pillow_coefficients(matrix)
    mine = warp.warp_nearest(img, matrix, fill=FILL)
    theirs = pillow_affine(img, coefficients)
    wrong = np.argwhere(mine != theirs)
    assert len(wrong) == 1
    oy, ox = int(wrong[0][0]), int(wrong[0][1])
    assert (oy, ox) == (4, 3)
    a, b, c, d, e, f = coefficients
    ys = d * (ox + 0.5) + e * (oy + 0.5) + f
    assert abs(ys - 5.0) < 1e-14
    assert ys != 5.0
    assert math.floor(ys) == 4


def test_pillow_agrees_a_full_turn_is_exact(img):
    matrix = warp.about_centre(warp.rotation(2.0 * math.pi), 9, 9)
    assert math.sin(2.0 * math.pi) != 0.0
    assert warp.matrices_close(matrix, warp.identity(), TOL)
    assert np.array_equal(warp.warp_nearest(img, matrix, fill=FILL), img)
    assert np.array_equal(
        pillow_affine(img, warp.to_pillow_coefficients(matrix)), img
    )


def test_a_non_square_output_agrees_about_which_way_round_the_size_goes(img):
    matrix = warp.shear_x(1.0)
    mine = warp.warp_nearest(img, matrix, out_shape=(9, 18), fill=FILL)
    theirs = pillow_affine(
        img, warp.to_pillow_coefficients(matrix), out_shape=(9, 18)
    )
    assert mine.shape == theirs.shape == (9, 18)
    assert np.array_equal(mine, theirs)


# -- Interpolation -----------------------------------------------------------


def test_bilinear_invents_values_that_nearest_neighbour_cannot(img):
    matrix = warp.translation(0.5, 0.0)
    nearest = warp.warp_nearest(img, matrix, fill=FILL)
    blended = warp.warp_bilinear_with_inverse(
        img, warp.invert(matrix), fill=float(FILL)
    )
    assert set(np.unique(nearest).tolist()) <= set(np.unique(img).tolist()) | {FILL}
    assert len(np.unique(np.round(blended, 6))) > len(np.unique(nearest))


def test_bilinear_at_a_whole_pixel_offset_reduces_to_nearest_neighbour(img):
    matrix = warp.translation(2.0, 1.0)
    nearest = warp.warp_nearest(img, matrix, fill=0)
    blended = warp.warp_bilinear_with_inverse(img, warp.invert(matrix), fill=0.0)
    # No fractional part, so every weight is 0 or 1 and no blending happens.
    assert np.abs(blended - nearest.astype(float)).max() <= 1e-9


@pytest.mark.parametrize(
    "name,matrix",
    [
        ("translate", warp.translation(0.25, 0.25)),
        ("rotate30", warp.about_centre(warp.rotation(math.radians(30)), 9, 9)),
        ("rotate17", warp.about_centre(warp.rotation(math.radians(17)), 9, 9)),
        ("scale", warp.about_centre(warp.scaling(1.5, 1.5), 9, 9)),
        ("shear", warp.shear_x(0.4)),
    ],
)
def test_bilinear_matches_pillow_wherever_all_four_neighbours_are_inside(
    img, name, matrix
):
    """The precise boundary of the agreement, stated as a tolerance.

    Where all four contributing pixels are inside the image, ours and Pillow
    agree to within 1 grey level -- the rounding of a float average into a
    byte. Where a contributor lies outside, they extrapolate differently.
    """
    inverse = warp.invert(matrix)
    mine = warp.warp_bilinear_with_inverse(img, inverse, fill=0.0)
    theirs = pillow_affine(
        img,
        warp.to_pillow_coefficients(matrix),
        resample=Image.Resampling.BILINEAR,
    ).astype(float)

    inside = np.zeros((9, 9), dtype=bool)
    for oy in range(9):
        for ox in range(9):
            sx, sy = warp.apply_point(inverse, (ox + 0.5, oy + 0.5))
            x0 = math.floor(sx - warp.SAMPLE_OFFSET)
            y0 = math.floor(sy - warp.SAMPLE_OFFSET)
            inside[oy, ox] = 0 <= x0 and x0 + 1 < 9 and 0 <= y0 and y0 + 1 < 9

    assert inside.any()
    assert np.abs(mine - theirs)[inside].max() <= 1.0


def test_bilinear_and_pillow_do_diverge_at_the_border(img):
    matrix = warp.about_centre(warp.rotation(math.radians(30)), 9, 9)
    mine = warp.warp_bilinear_with_inverse(img, warp.invert(matrix), fill=0.0)
    theirs = pillow_affine(
        img,
        warp.to_pillow_coefficients(matrix),
        resample=Image.Resampling.BILINEAR,
    ).astype(float)
    # Asserted rather than glossed over: the border difference is large.
    assert np.abs(mine - theirs).max() > 100.0


# -- Files and hygiene -------------------------------------------------------


def test_a_png_round_trip_is_lossless_and_leaves_nothing_behind(img):
    with tempfile.TemporaryDirectory() as tmp:
        path = os.path.join(tmp, "pattern.png")
        Image.fromarray(img, mode="L").save(path)
        assert os.path.getsize(path) > 0
        reloaded = np.asarray(Image.open(path).convert("L"))
        assert np.array_equal(reloaded, img)
    assert not os.path.exists(path)


def test_the_lab_writes_no_image_into_its_own_directory():
    """No image file in the lab's own tree -- the pattern is generated.

    `.venv` is skipped deliberately. It is the documented setup from the
    README, not litter, and Pillow ships a large collection of its own test
    images inside site-packages. Walking into it would fail this test for
    following the installation instructions.
    """
    here = os.path.dirname(os.path.abspath(__file__))
    lab = os.path.dirname(here)
    found = []
    for root, dirs, files in os.walk(lab):
        dirs[:] = [d for d in dirs if d != ".venv"]
        for name in files:
            if name.lower().endswith((".png", ".jpg", ".jpeg", ".bmp", ".gif")):
                found.append(os.path.join(root, name))
    assert found == []


def test_the_installed_versions_are_the_ones_this_lab_was_written_against():
    from importlib.metadata import version

    assert version("numpy").split(".")[0] == "2"
    assert int(version("pillow").split(".")[0]) >= 12

"""Image transformation from first principles: no NumPy in the matrix algebra.

The 3 by 3 matrices here are plain nested lists and the arithmetic is written
out by hand, for the same reason Day 102 gave: if `rotation` returned a NumPy
array built by a NumPy helper, then checking it against NumPy would be checking
NumPy against itself. The arrays that hold PIXELS are NumPy arrays, because an
image is exactly the kind of thing NumPy exists for. The MATHS is ours.

Two conventions, fixed here and obeyed everywhere:

1. A point is written (x, y): x is the COLUMN, y is the ROW. An image array is
   indexed `img[y, x]`. Row is y, column is x, and y grows downward.

2. A transformation matrix is a 3 by 3 homogeneous matrix that maps an INPUT
   point to an OUTPUT point -- the direction you can see. Translation is not a
   linear map (Day 102 proved a linear map cannot move the origin), so a third
   coordinate fixed at 1 is added and translation becomes an ordinary matrix
   multiply:

       [ a  b  tx ] [ x ]   [ a*x + b*y + tx ]
       [ c  d  ty ] [ y ] = [ c*x + d*y + ty ]
       [ 0  0   1 ] [ 1 ]   [        1       ]

   The top-left 2 by 2 block is exactly the linear part from Day 102 -- its
   columns are still where the basis vectors land. The third column is the
   translation, and the bottom row is always (0, 0, 1) for an affine map.
"""

import math

# Every output pixel is a little square. Its CENTRE, not its corner, is the
# point the transformation is evaluated at. Pixel (x, y) covers the square
# from (x, y) to (x + 1, y + 1), so its centre is at (x + 0.5, y + 0.5).
#
# This half is not a detail. It is what makes a shear coefficient move row 0
# even though the shear term is multiplied by y, and it is exactly what Pillow
# does -- verified in `06_against_pillow.py`, not assumed.
SAMPLE_OFFSET = 0.5

TOL = 1e-12


class SingularTransform(ValueError):
    """Raised when a transformation cannot be inverted, so it cannot be applied.

    Inheriting from ValueError mirrors numpy.linalg.LinAlgError, which Day 102
    showed is also a ValueError -- so existing `except ValueError` handlers
    keep working.
    """


# --------------------------------------------------------------------------
# Building the matrices
# --------------------------------------------------------------------------


def identity():
    """The transformation that changes nothing."""
    return [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]


def translation(tx, ty):
    """Move every point by (tx, ty). NOT a linear map -- it moves the origin.

    This is the whole reason homogeneous coordinates exist. In 2 by 2 there is
    no matrix that adds a constant; in 3 by 3 there is, and it is this one.
    """
    return [[1.0, 0.0, float(tx)], [0.0, 1.0, float(ty)], [0.0, 0.0, 1.0]]


def scaling(sx, sy):
    """Stretch x by sx and y by sy about the origin (the top-left corner)."""
    return [[float(sx), 0.0, 0.0], [0.0, float(sy), 0.0], [0.0, 0.0, 1.0]]


def rotation(theta):
    """Rotate by theta radians about the origin.

    In a y-UP coordinate system this turns counter-clockwise, which is what
    Day 102 drew. On an image y grows DOWNWARD, so the same matrix turns
    CLOCKWISE on screen. Nothing about the matrix changed; the picture is
    upside down relative to the graph paper.
    """
    cos_t, sin_t = math.cos(theta), math.sin(theta)
    return [[cos_t, -sin_t, 0.0], [sin_t, cos_t, 0.0], [0.0, 0.0, 1.0]]


def rotation_quarter_turns(turns):
    """Rotate by an exact multiple of 90 degrees, with integer entries.

    `rotation(math.pi / 2)` is correct but its cosine is 6.123233995736766e-17
    rather than 0.0 -- the Day 102 result. For the cases where the answer
    should be checkable to the exact pixel, build the matrix from integers
    instead of from trigonometry and the float noise never enters.
    """
    cos_t, sin_t = [(1, 0), (0, 1), (-1, 0), (0, -1)][turns % 4]
    return [
        [float(cos_t), float(-sin_t), 0.0],
        [float(sin_t), float(cos_t), 0.0],
        [0.0, 0.0, 1.0],
    ]


def shear_x(k):
    """Slide each row sideways in proportion to its y: x becomes x + k*y."""
    return [[1.0, float(k), 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]


def shear_y(k):
    """Slide each column vertically in proportion to its x: y becomes y + k*x."""
    return [[1.0, 0.0, 0.0], [float(k), 1.0, 0.0], [0.0, 0.0, 1.0]]


def flip_horizontal(width):
    """Mirror left-to-right inside an image `width` pixels wide.

    A bare reflection `x -> -x` sends the picture off the left edge. What is
    wanted is a reflection about the image's vertical centre line, which is a
    reflection FOLLOWED BY a translation of `width` -- and in homogeneous
    coordinates that is one matrix, not two steps.
    """
    return [[-1.0, 0.0, float(width)], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]


def flip_vertical(height):
    """Mirror top-to-bottom inside an image `height` pixels tall."""
    return [[1.0, 0.0, 0.0], [0.0, -1.0, float(height)], [0.0, 0.0, 1.0]]


def about_centre(matrix, width, height):
    """Do `matrix` about the image's centre instead of about its top-left corner.

    Move the centre to the origin, transform, move it back. Three matrices,
    composed into one:  T(+c) . M . T(-c).
    """
    cx, cy = width / 2.0, height / 2.0
    return compose(translation(cx, cy), matrix, translation(-cx, -cy))


# --------------------------------------------------------------------------
# Matrix arithmetic, written out
# --------------------------------------------------------------------------


def matmul(a, b):
    """The 3 by 3 product a @ b, computed by hand."""
    return [
        [sum(a[i][k] * b[k][j] for k in range(3)) for j in range(3)] for i in range(3)
    ]


def compose(*matrices):
    """Combine transformations into ONE matrix, applied RIGHT to LEFT.

    `compose(B, A)` means "do A first, then B" -- the same order as Day 101's
    matrix product and the same order as reading `B(A(x))` inside out.
    """
    if not matrices:
        return identity()
    result = matrices[0]
    for m in matrices[1:]:
        result = matmul(result, m)
    return result


def apply_point(matrix, point):
    """Send one (x, y) point through the matrix and return the new (x, y)."""
    x, y = point
    return (
        matrix[0][0] * x + matrix[0][1] * y + matrix[0][2],
        matrix[1][0] * x + matrix[1][1] * y + matrix[1][2],
    )


def determinant(matrix):
    """The determinant of the LINEAR part -- the area factor, as on Day 102.

    The third row of an affine matrix is (0, 0, 1), so the full 3 by 3
    determinant equals the 2 by 2 determinant of the top-left block.
    """
    return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]


def invert(matrix):
    """Invert an affine 3 by 3 matrix, using its structure rather than brute force.

    Write the matrix as a linear part A and a translation t. Then the inverse
    is the inverse of A, followed by minus that inverse applied to t:

        M = [ A  t ]      M^-1 = [ A^-1   -A^-1 t ]
            [ 0  1 ]             [  0         1   ]

    Raises SingularTransform when the determinant is 0 -- a transformation that
    flattens the image onto a line throws information away, and no amount of
    arithmetic can put it back.
    """
    det = determinant(matrix)
    if abs(det) <= TOL:
        raise SingularTransform(
            f"determinant is {det!r}: this transformation collapses the image "
            "and cannot be undone"
        )
    a, b, tx = matrix[0]
    c, d, ty = matrix[1]
    ia, ib = d / det, -b / det
    ic, idd = -c / det, a / det
    # `+ 0.0` normalises negative zero. Negating a 0.0 entry -- which every
    # translation and every axis-aligned scale has -- produces -0.0, which
    # compares equal to 0.0 but PRINTS as "-0.0". That is noise in every
    # coefficient tuple this function feeds, so it is cleaned up once, here.
    return [
        [ia + 0.0, ib + 0.0, -(ia * tx + ib * ty) + 0.0],
        [ic + 0.0, idd + 0.0, -(ic * tx + idd * ty) + 0.0],
        [0.0, 0.0, 1.0],
    ]


def matrices_close(a, b, tol=TOL):
    """True when two 3 by 3 matrices agree entry by entry within `tol`."""
    return all(abs(a[i][j] - b[i][j]) <= tol for i in range(3) for j in range(3))


# --------------------------------------------------------------------------
# The two ways to move pixels
# --------------------------------------------------------------------------


def warp_forward(image, matrix, out_shape=None, fill=0):
    """Forward mapping: push every INPUT pixel to where it lands. This is wrong.

    It is written out in full because seeing it fail is the argument for the
    method that follows. Walk the input, send each pixel through the matrix,
    round to the nearest output pixel, and write the value there.

    The output is full of holes. Nothing in the arithmetic guarantees that the
    input pixels land on output pixels one-to-one: a rotation spreads them out
    so some output pixels are never written at all, and a shrink piles several
    input pixels onto the same output pixel so others are written repeatedly.

    Returns (output_image, hole_mask) where hole_mask is a boolean array that
    is True at every output pixel no input pixel ever reached.
    """
    src = _as_2d(image)
    height, width = src.shape
    out_h, out_w = out_shape or (height, width)

    out = _new_like(src, (out_h, out_w), fill)
    written = _zeros_bool(out_h, out_w)

    for y in range(height):
        for x in range(width):
            # The input pixel's centre, sent through the matrix.
            fx, fy = apply_point(
                matrix, (x + SAMPLE_OFFSET, y + SAMPLE_OFFSET)
            )
            # Which output pixel square does that point fall in?
            ox, oy = math.floor(fx), math.floor(fy)
            if 0 <= ox < out_w and 0 <= oy < out_h:
                out[oy, ox] = src[y, x]
                written[oy, ox] = True
    return out, ~written


def warp_nearest(image, matrix, out_shape=None, fill=0):
    """Inverse mapping with nearest-neighbour sampling. This is the right way.

    Walk the OUTPUT, not the input. For each output pixel, take its centre,
    send it BACKWARD through the inverse matrix to find the place in the input
    it came from, and take the value of whichever input pixel contains that
    place. Every output pixel is visited exactly once, so there are no holes --
    not because the holes were patched, but because the loop is over the array
    being filled.

    `matrix` is the transformation you can SEE: input to output. The inverse is
    taken here, once, rather than being demanded of the caller.
    """
    inverse = invert(matrix)
    return warp_nearest_with_inverse(image, inverse, out_shape=out_shape, fill=fill)


def warp_nearest_with_inverse(image, inverse, out_shape=None, fill=0):
    """The same as `warp_nearest`, but given the OUTPUT-to-INPUT matrix directly.

    This is the form Pillow's `Image.transform` takes its coefficients in, and
    having it separately is what lets `06_against_pillow.py` hand the two
    implementations the identical six numbers.
    """
    src = _as_2d(image)
    height, width = src.shape
    out_h, out_w = out_shape or (height, width)
    out = _new_like(src, (out_h, out_w), fill)

    for oy in range(out_h):
        for ox in range(out_w):
            sx, sy = apply_point(
                inverse, (ox + SAMPLE_OFFSET, oy + SAMPLE_OFFSET)
            )
            # floor, not round: the input pixel whose SQUARE contains the point.
            ix, iy = math.floor(sx), math.floor(sy)
            if 0 <= ix < width and 0 <= iy < height:
                out[oy, ox] = src[iy, ix]
            # Otherwise leave the fill value: the source lies outside the
            # picture. This is the clipping the corners of a rotated image run
            # into, and it is a decision, not an error.
    return out


def warp_bilinear_with_inverse(image, inverse, out_shape=None, fill=0.0):
    """Inverse mapping again, but blending the four pixels around the landing point.

    The inverse-mapped position almost never lands on a pixel exactly.
    Nearest-neighbour answers "which pixel is closest"; bilinear answers "what
    would the value be here", by taking a weighted average of the four
    surrounding pixels, weighted by how close the point is to each.

    The half-pixel bookkeeping is the fiddly part. Pixel (i, j) has its VALUE
    at its centre, (i + 0.5, j + 0.5). So to interpolate between pixel centres,
    subtract the half back off before splitting into whole and fractional
    parts. Getting this wrong shifts the whole image by half a pixel, which
    looks like a mysterious blur rather than like an offset.

    Returns a float array, because an average of integers is not an integer.
    Out-of-range contributions count as `fill`.
    """
    src = _as_2d(image).astype(float)
    height, width = src.shape
    out_h, out_w = out_shape or (height, width)
    out = _zeros_float(out_h, out_w)

    for oy in range(out_h):
        for ox in range(out_w):
            sx, sy = apply_point(
                inverse, (ox + SAMPLE_OFFSET, oy + SAMPLE_OFFSET)
            )
            gx, gy = sx - SAMPLE_OFFSET, sy - SAMPLE_OFFSET
            x0, y0 = math.floor(gx), math.floor(gy)
            tx, ty = gx - x0, gy - y0
            total = 0.0
            for dy in (0, 1):
                for dx in (0, 1):
                    weight = (tx if dx else 1.0 - tx) * (ty if dy else 1.0 - ty)
                    px, py = x0 + dx, y0 + dy
                    inside = 0 <= px < width and 0 <= py < height
                    total += weight * (src[py, px] if inside else float(fill))
            out[oy, ox] = total
    return out


def to_pillow_coefficients(matrix):
    """Turn a visible input-to-output matrix into Pillow's six coefficients.

    Pillow's `Image.transform(..., Image.Transform.AFFINE, coeffs)` takes
    `(a, b, c, d, e, f)` meaning

        input_x = a * output_x + b * output_y + c
        input_y = d * output_x + e * output_y + f

    -- the OUTPUT-to-INPUT direction, which is the inverse of the effect you
    see. Day 102 verified that direction by experiment. So the coefficients are
    read off the INVERSE of the matrix, and forgetting to invert is the single
    most common way to get a Pillow transform backwards.
    """
    inverse = invert(matrix)
    return (
        inverse[0][0],
        inverse[0][1],
        inverse[0][2],
        inverse[1][0],
        inverse[1][1],
        inverse[1][2],
    )


def coefficients_to_matrix(coeffs):
    """The reverse of `to_pillow_coefficients`: six numbers back to a 3 by 3."""
    a, b, c, d, e, f = coeffs
    return [[a, b, c], [d, e, f], [0.0, 0.0, 1.0]]


# --------------------------------------------------------------------------
# Small array helpers, kept apart so the maths above reads cleanly
# --------------------------------------------------------------------------


def _as_2d(image):
    import numpy as np

    arr = np.asarray(image)
    if arr.ndim != 2:
        raise ValueError(
            f"expected a 2-D greyscale array of shape (height, width), got "
            f"shape {arr.shape}. For a colour image, transform each of the "
            f"three planes: img[:, :, 0], img[:, :, 1], img[:, :, 2]."
        )
    return arr


def _new_like(src, shape, fill):
    import numpy as np

    return np.full(shape, fill, dtype=src.dtype)


def _zeros_bool(h, w):
    import numpy as np

    return np.zeros((h, w), dtype=bool)


def _zeros_float(h, w):
    import numpy as np

    return np.zeros((h, w), dtype=float)


def warp_colour(image, matrix, out_shape=None, fill=0):
    """Transform a (height, width, 3) colour image plane by plane.

    There is no new mathematics here, and that is the point worth noticing: the
    transformation acts on COORDINATES, and the three planes share their
    coordinates, so the same matrix does all three.
    """
    import numpy as np

    arr = np.asarray(image)
    if arr.ndim != 3 or arr.shape[2] != 3:
        raise ValueError(
            f"expected a colour array of shape (height, width, 3), got {arr.shape}"
        )
    planes = [
        warp_nearest(arr[:, :, c], matrix, out_shape=out_shape, fill=fill)
        for c in range(3)
    ]
    return np.stack(planes, axis=2)

lesson_id: D105
day: 105
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-105-transforming-images-with-matrices
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - .venv/bin/python3 -c "import numpy, PIL; print(numpy.__version__, PIL.__version__)"
run_commands:
  - 'cd examples && ../.venv/bin/python3 01_an_image_is_a_matrix.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_forward_mapping_leaves_holes.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_inverse_mapping.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_scale_shear_flip.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_homogeneous_and_composition.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_against_pillow.py && cd ..'
  - .venv/bin/pytest examples -q -p no:cacheprovider
  - .venv/bin/pytest starter -q -p no:cacheprovider
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - "find . -type d -name '__pycache__' -prune -exec rm -rf -- {} +"
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 35
last_executed: '2026-08-17'
executed_on: 'macOS 26.5.2 (Apple Silicon, arm64), Python 3.14.0, numpy 2.5.2, Pillow 12.3.0, pytest 9.1.1, bash 3.2.57 — bash tests/run_tests.sh -> 79 checks, 0 failure(s), exit 0; pytest examples -> 64 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, Pillow and pytest; nothing else in the lab opens a socket, and the test image is GENERATED in code rather than downloaded — section 7 of the harness greps the sources for network calls and also asserts that no image file exists anywhere under the lab. Section 6 re-runs the harness with one expectation deliberately swapped for the naive belief that Pillow samples at integer pixel corners, 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. Three measured results are asserted rather than smoothed over. (1) Pillow evaluates an affine transform at each output pixel CENTRE, (x + 0.5, y + 0.5), and floors — this settles the half-pixel sampling question Day 102 raised and deliberately deferred, and it is why a shear coefficient of 2.0 moves row 0 by one whole pixel despite the shear term being multiplied by y. (2) This lab implementation and Pillow produce byte-for-byte identical output on 510 of 510 affine transformations, but disagree on 8 of the 360 whole-degree rotations (30, 60, 120, 150, 210, 240, 300, 330) by at most 2 pixels of 81; every disagreeing sample lands within 2.220e-15 of a pixel boundary, where the order of floating-point additions decides which side of the tie you get. The claim asserted is the boundary property, not the specific angle list. (3) The from-scratch bilinear agrees with Pillow within 1.0 grey level wherever all four contributing pixels are inside the image, and diverges by up to 118 levels at the border because the two extrapolate differently — both halves asserted rather than the tolerance being widened.'

# Dependencies for the Day 105 lab

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

| Package | Pinned version | Why this lab needs it |
| --- | --- | --- |
| `numpy` | `2.5.2` | Holds the pixels. An image is a 2-D array of bytes and this is what NumPy exists for. Also supplies the independent answers the lab checks against: `numpy.rot90`, `numpy.fliplr`, `numpy.flipud`, `numpy.kron` and ordinary strided slicing. |
| `pillow` | `12.3.0` | The mature library your from-scratch code is compared against. `Image.transform` with `Image.Transform.AFFINE` is the whole of the comparison, plus one PNG round trip. |
| `pytest` | `9.1.1` | The test runner from Days 071–074. Nothing new here except what it is pointed at. |

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

`examples/warp.py` and `starter/warp.py` build and multiply their 3 by 3
matrices with plain lists and hand-written arithmetic. If `rotation` returned a
NumPy array built by a NumPy helper, then checking it against NumPy would be
checking NumPy against itself, and agreeing with Pillow would prove nothing
either — both sit on the same numerical machinery.

NumPy holds the **pixels**, because a 9 by 9 array of bytes is exactly what it
is for and writing that by hand would teach nothing. The **mathematics** is
ours. That split is what makes section 3 of `06_against_pillow.py` mean
something: 510 transformations, byte-for-byte identical output, from two
implementations that share no code.

## Why the versions are pinned

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

Two places the version could genuinely matter, both handled honestly rather
than pinned to a last digit:

1. **Pillow's sampling rule.** This lab establishes by measurement that Pillow
   evaluates an affine transform at each output pixel's *centre* and takes the
   input pixel whose square contains the result. That was measured on Pillow
   12.3.0 on the authoring machine. It is a long-standing behaviour rather than
   a documented guarantee, so the lab measures it every run instead of
   asserting it from memory — `test_pillow_samples_at_pixel_centres_not_at_integer_corners`
   will fail loudly if a future version changes it, which is the correct
   outcome.

2. **Floating-point tie-breaking.** Eight of the 360 whole-degree rotations
   produce output that differs between this lab's implementation and Pillow's,
   by at most 2 pixels out of 81, and every one of those disagreements is a
   sample landing within one ulp of a pixel boundary. That set of eight angles
   is specific to this machine's floating-point behaviour and this version of
   Pillow. The lab asserts the *count* and asserts that every disagreement is
   at a boundary; it does not claim the two implementations agree everywhere,
   because they do not.

The versions were read from the installed packages rather than guessed:

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

On the authoring machine, on 17 August 2026, that printed
`2.5.2 12.3.0 9.1.1`.

## Licences

NumPy is distributed under the BSD 3-Clause licence, Pillow under the MIT-CMU
licence, and pytest under the MIT licence, each stated on that project's own
documentation site. All three 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, PIL; print(numpy.__version__, PIL.__version__)"
```

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

## Network

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

In particular, the test image is *generated in code*, not downloaded. That was
a deliberate choice: a lab that fetches a photograph is a lab that breaks on a
train, ships a file whose licence someone has to check, and hides its own test
data behind a URL. `pattern.py` builds a 9 by 9 capital F from arithmetic, and
every pixel value in it is asserted. Section 7 of `tests/run_tests.sh` greps
every file under `examples/` and `starter/` for the patterns that would
indicate a socket being opened.

## Running without a lab-local environment

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

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

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

## What you would give up without Pillow

Less than you might fear, and the part you would lose is the best part.

Exercise 1 — writing all twelve functions — needs `math` and NumPy only, and
every test of your rotation, scaling, shear, flip, composition and inverse
mapping still runs, because those are checked against `numpy.rot90`,
`numpy.fliplr`, `numpy.kron` and strided slices rather than against Pillow.

What you lose is the comparison: the 510-transformation agreement, the
measurement that settles where Pillow takes its sample, the shear that moves
row 0, and the bilinear border finding. Those are the day's strongest
artifacts and they cannot be faked. The lab does not pretend otherwise.

numpy==2.5.2
pillow==12.3.0
pytest==9.1.1

# Day 105 lab — Rotate It Yourself

One idea holds this whole lab together:

> **An image is a matrix, so you transform its COORDINATES, not its pixels.**

Everything below is a consequence of that sentence. You will rotate, scale,
shear and flip a picture with matrices you wrote yourself, and then hand the
identical six numbers to Pillow and watch a library maintained since 2010 agree
with you pixel for pixel.

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.

---

## The picture

Nothing is downloaded. `pattern.py` builds the test image from arithmetic — a
capital F on a 9 by 9 greyscale grid:

```
.######..
.##......
.##......
.##......
.####....
.##......
.##......
.##......
.##.....o
```

`.` is 0, `#` is 255, `o` is 96. It is an F for one reason: an F is asymmetric
under every operation in this lab. A square survives a horizontal flip
unchanged, so a square would let a broken flip pass its own test. An F does not
survive anything.

Print your own results the same way at any time:

```python
import pattern
print(pattern.as_text(my_image))
```

Anything that is not one of the four known grey levels prints as `?`, which is
how an interpolation artefact announces itself.

---

## The two conventions, and the one that will trip you

**Rows are y, columns are x.** A NumPy array is indexed `img[row, column]`,
which is `img[y, x]`. All week a point has been written `(x, y)`. Those are the
same two numbers in the opposite order, and swapping them does not raise — it
silently reads the wrong pixel. This is the single most common source of
confusion in the subject.

**y grows downward.** Row 0 is the top of the picture. Day 102's graphs had y
growing upward. Nothing about the matrices changes, but a counter-clockwise
rotation matrix turns an image *clockwise* on screen. That is not a sign error
and you should not "fix" it.

---

## Exercise 1 — `warp.py` (twelve functions)

Write the twelve functions marked `raise NotImplementedError`. Each docstring
gives the derivation and a worked example you can check on paper. Use only
`math` for the matrix arithmetic; NumPy is fine for holding pixels, which is
what it is for.

| Step | Function | What it must do |
| --- | --- | --- |
| 1.1 | `translation` | The reason homogeneous coordinates exist. Add a constant by multiplying the third coordinate. |
| 1.2 | `scaling` | Derive it: where does one step right go when the picture stretches? |
| 1.3 | `rotation` | Day 102's matrix, dropped into the top-left of a 3 by 3. |
| 1.4 | `shear_x` | Derive it: the sideways push is proportional to the **height**. |
| 1.5 | `flip_horizontal` | A mirror **plus** a translation, in one matrix. Getting this wrong sends the picture off the left edge. |
| 1.6 | `matmul` | The 3 by 3 product, by hand. `compose` is written for you on top of it. |
| 1.7 | `apply_point` | The third coordinate is always 1, so you never build the triple. |
| 1.8 | `determinant` | The area factor. Compute it directly so whole numbers stay exact. |
| 1.9 | `invert` | Use the affine structure, not brute force. Raise `SingularTransform` at a determinant of 0. |
| 1.10 | `warp_forward` | Forward mapping, **done wrong on purpose**. Return the holes; do not patch them. |
| 1.11 | `warp_nearest_with_inverse` | The real one. Loop over the OUTPUT. |
| 1.12 | `to_pillow_coefficients` | Read the six numbers off the **inverse**. |

Seven helpers are written for you at the bottom of the file — `identity`,
`rotation_quarter_turns`, `shear_y`, `flip_vertical`, `about_centre`,
`matrices_close` and `coefficients_to_matrix`. Read them; the tests use them,
and `about_centre` in particular is worth understanding because it is three of
your matrices folded into one.

### The three places people lose an afternoon

**1.10 is meant to fail.** `warp_forward` leaves 22 of 81 output pixels
unwritten on a 30 degree rotation, including pixels punched through the middle
of solid ink. That is the correct answer. Do not add a second pass to fill
them; the whole point of 1.11 is that turning the loop inside out makes the
problem disappear rather than needing a patch.

**1.11 needs the half.** Every output pixel is a little square, and you sample
its **centre**, `(x + 0.5, y + 0.5)`, not its corner. `SAMPLE_OFFSET` is
already defined at the top of `warp.py`. Leave it out and every result is half
a pixel adrift — which looks like a mysterious blur rather than like an offset,
and which is exactly why `test_1_11_halving_is_exactly_a_strided_slice` exists.

**1.11 needs `math.floor`.** Not `round`, not `int`. You want the input pixel
whose *square contains* the point. `int` truncates toward zero, which is wrong
for negative coordinates; `round` is a different rule that will disagree with
Pillow on half the cases.

---

## Exercises 2 to 6 — `answers.py` (predictions)

Twenty-six predictions, every one of which can be reasoned out on paper before
you run anything. Replace each `None` with your answer. Anything still `None`
is skipped rather than failed, so your score only ever counts work you actually
attempted.

| Exercise | About |
| --- | --- |
| 2 | An image is a matrix, and the `(row, column)` versus `(x, y)` trap. |
| 3 | Why forward mapping leaves holes, and why the fill-valued pixels after inverse mapping are *not* holes. |
| 4 | Inverse mapping, and four results that are exactly checkable against NumPy. |
| 5 | Why no 2 by 2 matrix can translate, and what composition buys you. |
| 6 | Pillow's coefficient direction, and where exactly it takes its sample. |

Exercise 6 is the one to slow down on. Question 6.4 asks by how many pixels row
0 moves under a shear whose coefficient is 2.0. The mathematics says row 0 has
`y = 0` and therefore cannot move. Predict the number *before* you run
anything, then check it. Day 102 raised this question and deliberately left it
open; today you settle it.

---

## Finishing

```bash
bash tests/run_tests.sh
```

The harness runs the six reference scripts, both pytest suites, and a set of
checks that read real values rather than reading source. It prints
`N checks, 0 failure(s)` and exits 0 when everything holds.

If you want to see the finished versions, `examples/` has them — `warp.py`
there is the complete implementation and the six numbered scripts walk through
every result in this brief with the numbers printed. Read them **after** you
have attempted the exercises; reading them first turns a lab into a
transcription task.

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

Every one of these can be reasoned out on paper. That is the point: a lab about
image 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
"""

# =============================================================================
# Exercise 2 -- an image is a matrix, and the ordering trap
# =============================================================================
#
# The test pattern is a capital F drawn on a 9 by 9 greyscale grid:
#
#     .######..
#     .##......
#     .##......
#     .##......
#     .####....
#     .##......
#     .##......
#     .##......
#     .##.....o
#
# where `.` is 0, `#` is 255 and `o` is 96. Rows are numbered 0 to 8 from the
# TOP; columns 0 to 8 from the LEFT.

# 2.1 What does `img.shape` report? A tuple of two ints.
#     Careful: NumPy reports (height, width), not (width, height).
SHAPE = None

# 2.2 The colour version stacks three of those planes. What is ITS shape?
COLOUR_SHAPE = None

# 2.3 The `o` pixel sits at row 8, column 8. Written as a POINT in the
#     language of Week 15 -- (x, y) -- what is it? A tuple of two ints.
MARK_AS_POINT = None

# 2.4 One ink pixel of the F sits at row 4, column 3 (in the middle bar).
#     What is the VALUE of img[3, 4] -- the same two numbers, swapped?
#     An int. This is the trap: swapping them does not raise, it returns a
#     wrong answer silently.
VALUE_AT_SWAPPED_INDEX = None

# 2.5 How many pixels of the image are ink (value 255)? An int.
#     Count them off the picture above: the top bar, the two-wide stem down
#     all nine rows, and the middle bar. Mind the overlaps.
INK_PIXEL_COUNT = None


# =============================================================================
# Exercise 3 -- forward mapping leaves holes
# =============================================================================
#
# Forward mapping walks the INPUT and writes each pixel where it lands.

# 3.1 Scale the 9 by 9 image up by 2 into an 18 by 18 output, by forward
#     mapping. AT LEAST how many output pixels must be holes? An int.
#     This is a counting argument, not a rounding one: how many output pixels
#     are there, and how many input pixels are available to fill them?
MINIMUM_HOLES_WHEN_DOUBLING = None

# 3.2 Now shrink instead: forward-map into a 5 by 5 output. How many holes?
#     An int. Think about whether there are enough input pixels this time.
HOLES_WHEN_SHRINKING = None

# 3.3 Shrinking leaves no holes but loses information a different way. In one
#     word, what happens when two input pixels land on the same output pixel?
#     Answer with the string "overwriting" or "blending" -- only one of them
#     is what the algorithm in warp_forward actually does.
WHAT_SHRINKING_DOES = None

# 3.4 Inverse mapping leaves SOME output pixels at the fill value after a 30
#     degree rotation. Is that the same failure as a hole?
#     Answer True if those pixels are holes, or False if they are something
#     else. If False, exercise 3.5 asks you to name it.
FILL_PIXELS_ARE_HOLES = None

# 3.5 One word for what those fill-valued pixels actually are: the source
#     position fell outside the input image. The string is one of
#     "clipping", "aliasing", "quantisation".
NAME_FOR_FILL_PIXELS = None


# =============================================================================
# Exercise 4 -- inverse mapping, with exact answers
# =============================================================================

# 4.1 Rotate the F a quarter turn using rotation_quarter_turns(1) about the
#     image centre. The corner mark starts at (row 8, column 8). Where does it
#     end up? A tuple (row, column) of two ints.
#     Remember: y grows downward, so the counter-clockwise matrix turns the
#     PICTURE clockwise.
MARK_AFTER_QUARTER_TURN = None

# 4.2 That same quarter turn equals one of NumPy's own rotations exactly.
#     Which value of k makes numpy.rot90(img, k) identical to your result?
#     An int in the range -3 to 3.
NUMPY_ROT90_K = None

# 4.3 How many pixels take the fill value after a quarter turn of a SQUARE
#     image? An int.
FILL_COUNT_AFTER_QUARTER_TURN = None

# 4.4 Scale up by exactly 2 with nearest-neighbour into an 18 by 18 output.
#     How many DISTINCT pixel values does the result contain? An int.
#     The input contains three. Does enlarging invent any new ones?
DISTINCT_VALUES_AFTER_DOUBLING = None

# 4.5 Scale DOWN by exactly a half into a 4 by 4 output. The result turns out
#     to be a plain NumPy strided slice of the input. Which one?
#     Answer with the string "img[0::2, 0::2]" or "img[1::2, 1::2]".
#     Hint: output pixel 0 is sampled at its centre, 0.5, which doubles to 1.0.
DOWNSCALE_IS_THE_SLICE = None


# =============================================================================
# Exercise 5 -- homogeneous coordinates and composition
# =============================================================================

# 5.1 Why can no 2 by 2 matrix perform a translation? Answer with the string
#     "it cannot move the origin", "it cannot change area", or
#     "it cannot rotate".
WHY_2X2_CANNOT_TRANSLATE = None

# 5.2 What is the determinant of translation(7, -3)? A float.
DETERMINANT_OF_A_TRANSLATION = None

# 5.3 `compose(B, A)` applies which one first? Answer with the string "A"
#     or "B".
COMPOSE_APPLIES_FIRST = None

# 5.4 Rotate the image 30 degrees twelve times in a row, resampling each time,
#     then compare with the original. Will it come back exactly?
#     True or False.
TWELVE_SEPARATE_ROTATIONS_ARE_EXACT = None

# 5.5 Compose those same twelve rotations into ONE matrix and apply it once.
#     Will THAT come back exactly? True or False.
#     If your two answers differ, you have understood the day's main practical
#     lesson.
TWELVE_COMPOSED_ROTATIONS_ARE_EXACT = None


# =============================================================================
# Exercise 6 -- against Pillow, and the half-pixel question
# =============================================================================

# 6.1 Pillow's affine coefficients (a, b, c, d, e, f) express the map in which
#     direction? Answer with the string "input to output" or
#     "output to input".
PILLOW_COEFFICIENT_DIRECTION = None

# 6.2 So if you pass c = +1 with everything else at identity, which way does
#     the picture appear to move? The string "left" or "right".
PICTURE_MOVES_WHEN_C_IS_POSITIVE = None

# 6.3 Pillow evaluates the transformation at which point of each output pixel?
#     The string "its top-left corner" or "its centre".
#     Exercise 6 in the reference scripts settles this by measurement; predict
#     it first.
PILLOW_SAMPLES_AT = None

# 6.4 A shear whose coefficient b is 2.0. The mathematics says row 0 has y = 0
#     and therefore does not move. By how many whole pixels does row 0
#     ACTUALLY move? An int.
#     Work it out from your answer to 6.3.
ROW_ZERO_SHIFT_WITH_B_EQUALS_2 = None

# 6.5 Rotate by 2*pi -- a full turn -- as a single matrix, with
#     nearest-neighbour. How many pixels of the 81 differ from the original?
#     An int.
#     math.sin(2*math.pi) is about -2.4e-16, not 0. Does an error that size
#     survive rounding to a whole pixel?
PIXELS_CHANGED_BY_A_FULL_TURN = None

# 6.6 Ours and Pillow disagree on 8 of the 360 whole-degree rotations. Those
#     eight are 30, 60, 120, 150, 210, 240, 300 and 330 degrees. What do those
#     angles have in common? Answer with the string
#     "their sines and cosines land samples exactly on pixel boundaries",
#     "they are all multiples of 30", or
#     "they are randomly distributed".
#     The second is a true statement about the list and explains nothing --
#     90, 180 and 270 are multiples of 30 too, and they agreed. Pick the one
#     that says WHY.
WHY_THOSE_ANGLES_DISAGREE = None

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

Both `examples/` and `starter/` contain modules called `warp` and `pattern`,
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 `warp` 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
`warp` or `pattern` 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 ("warp", "pattern", "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 test image, generated in code rather than downloaded.

Nothing here loads a file. The pattern is built from arithmetic so that every
pixel value in this lab has a reason, and so that the lab needs no network and
ships no photograph.

The glyph is a capital F on a 9 by 9 grid. An F is the standard test shape in
image processing for one reason: it is asymmetric under every operation this
lab performs. A square survives a horizontal flip unchanged, so a square would
let a broken flip pass. An F does not survive anything -- flip it, rotate it,
transpose it, and you can see immediately which one happened.

Coordinates. This is the single most common source of confusion in the whole
subject, so it is stated once, here, and every function in `warp.py` obeys it:

    image[row, column]  ==  image[y, x]

A NumPy array is indexed rows first. The maths of Week 15 wrote points as
(x, y). Those are the SAME two numbers in the OPPOSITE order. Row is y, column
is x. And y grows DOWNWARD, because row 0 is the top of the picture, which is
the reverse of the graphs on Day 102 where y grew upward.
"""

import numpy as np

# Grey levels. Distinct, well separated, and none of them 0 or 255 by accident.
INK = 255  # the glyph itself
PAPER = 0  # the background
MARK = 96  # a single corner pixel, so a 180 degree turn is distinguishable
FILL = 32  # what lands in an output pixel whose source is off the image

SIZE = 9  # the pattern is SIZE by SIZE

# The exact cells of the F, written out rather than computed, so the values in
# the tests can be read off this list by eye.
TOP_BAR = [(0, c) for c in range(1, 7)]  # row 0, columns 1..6
MIDDLE_BAR = [(4, c) for c in range(1, 5)]  # row 4, columns 1..4
STEM = [(r, 1) for r in range(1, 9)] + [(r, 2) for r in range(1, 9)]

# Where the single corner mark goes: bottom-right, the corner the F never
# reaches, so it can never be confused with part of the glyph.
MARK_CELL = (8, 8)


def make_pattern():
    """Return the 9 by 9 greyscale test image as a uint8 array of shape (9, 9).

    Shape is (height, width) -- rows then columns -- which is (y, x).
    """
    img = np.full((SIZE, SIZE), PAPER, dtype=np.uint8)
    for row, col in TOP_BAR + MIDDLE_BAR + STEM:
        img[row, col] = INK
    img[MARK_CELL] = MARK
    return img


def make_colour_pattern():
    """Return the same glyph as a colour image of shape (9, 9, 3).

    Colour is three greyscale planes stacked on the last axis: red, green,
    blue. Each plane is a matrix in its own right, and every transformation in
    this lab acts on the coordinates, which the three planes share. That is
    why transforming a colour image is the same work as transforming a
    greyscale one, done three times.
    """
    grey = make_pattern()
    img = np.zeros((SIZE, SIZE, 3), dtype=np.uint8)
    img[:, :, 0] = grey  # red plane: the glyph
    img[:, :, 1] = np.fliplr(grey)  # green plane: the glyph mirrored
    img[:, :, 2] = MARK  # blue plane: flat, so a channel mix-up is obvious
    return img


def as_text(img, ink_char="#", paper_char=".", mark_char="o", fill_char="~"):
    """Render a greyscale array as ASCII so a transformation can be SEEN.

    Anything that is not one of the four known levels prints as `?`, which is
    how an interpolation artefact announces itself.
    """
    table = {INK: ink_char, PAPER: paper_char, MARK: mark_char, FILL: fill_char}
    return "\n".join(
        "".join(table.get(int(v), "?") for v in row) for row in np.asarray(img)
    )


def ink_cells(img):
    """Return the sorted (row, column) pairs whose value is INK.

    Comparing two images by their ink cells is exact -- these are integers, not
    floats -- which is why the rotation tests in this lab can assert equality
    rather than a tolerance.
    """
    rows, cols = np.nonzero(np.asarray(img) == INK)
    return sorted(zip(rows.tolist(), cols.tolist()))


# The values the tests assert against, written here once so that a change to
# the pattern cannot quietly change what "correct" means somewhere else.
EXPECTED_SHAPE = (9, 9)
EXPECTED_COLOUR_SHAPE = (9, 9, 3)
EXPECTED_INK_COUNT = len(set(TOP_BAR + MIDDLE_BAR + STEM))
EXPECTED_TEXT = "\n".join(
    [
        ".######..",
        ".##......",
        ".##......",
        ".##......",
        ".####....",
        ".##......",
        ".##......",
        ".##......",
        ".##.....o",
    ]
)

# Tolerances. Every float comparison in this lab names one of these.
TOL = 1e-12  # for matrix arithmetic done in floating point
PIXEL_TOL = 0  # for nearest-neighbour pixel values: they must match EXACTLY

"""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 the tolerance TOL stated in warp.py. Pixel comparisons
are EXACT, because nearest-neighbour produces whole pixel values and there is
nothing to round.
"""

import math

import numpy as np
import pytest
from PIL import Image

import answers
import pattern
import warp

TOL = warp.TOL
FILL = pattern.FILL


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 centred(matrix, width=9, height=9):
    """`about_centre`, routed through `written` so an unwritten helper skips.

    `about_centre` is written for you, but it calls YOUR `compose`, `matmul`
    and `translation`. Calling it directly at an argument position would let a
    NotImplementedError escape and be reported as a failure rather than as
    "not attempted", so every test goes through here instead.
    """
    return written(warp.about_centre, matrix, width, height)


def close(a, b, tol=TOL):
    """Elementwise closeness for points and for 3 by 3 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))


def pillow_affine(array, coefficients, out_shape=None, resample=None):
    height, width = array.shape
    out_h, out_w = out_shape or (height, width)
    return np.asarray(
        Image.fromarray(array, mode="L").transform(
            (out_w, out_h),
            Image.Transform.AFFINE,
            coefficients,
            resample=resample or Image.Resampling.NEAREST,
            fillcolor=FILL,
        )
    )


@pytest.fixture
def img():
    return pattern.make_pattern()


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


def test_0_the_environment_is_ready():
    """Always passes once the install worked. Everything below is your work."""
    import PIL

    assert np.__version__, "numpy is importable"
    assert PIL.__version__, "Pillow is importable"
    assert warp.identity() == [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
    assert warp.SAMPLE_OFFSET == 0.5
    assert pattern.make_pattern().shape == (9, 9)
    assert pattern.as_text(pattern.make_pattern()) == pattern.EXPECTED_TEXT


# -- Exercise 1.1: translation ------------------------------------------------


def test_1_1_translation_has_the_right_shape():
    M = written(warp.translation, 3.0, -2.0)
    assert len(M) == 3 and all(len(row) == 3 for row in M), "3 by 3, rows of 3"
    assert M[2] == [0.0, 0.0, 1.0], "the bottom row of an affine matrix"


def test_1_1_translation_is_derived_correctly():
    M = written(warp.translation, 3.0, -2.0)
    assert close(M, [[1.0, 0.0, 3.0], [0.0, 1.0, -2.0], [0.0, 0.0, 1.0]])


# -- Exercise 1.2: scaling ----------------------------------------------------


def test_1_2_scaling_is_derived_correctly():
    M = written(warp.scaling, 2.0, 3.0)
    assert close(M, [[2.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 1.0]])


# -- Exercise 1.3: rotation ---------------------------------------------------


def test_1_3_rotation_is_derived_correctly():
    M = written(warp.rotation, math.pi / 2)
    assert close(M, [[0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]])


def test_1_3_rotation_by_zero_is_the_identity():
    M = written(warp.rotation, 0.0)
    assert close(M, warp.identity())


# -- Exercise 1.4: shear ------------------------------------------------------


def test_1_4_shear_is_derived_correctly():
    M = written(warp.shear_x, 2.0)
    assert close(M, [[1.0, 2.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])


# -- Exercise 1.5: flip -------------------------------------------------------


def test_1_5_flip_mirrors_about_the_centre_line_not_about_zero():
    M = written(warp.flip_horizontal, 9.0)
    assert close(M, [[-1.0, 0.0, 9.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])


# -- Exercise 1.6: matrix product --------------------------------------------


def test_1_6_matmul_matches_numpy():
    a = [[1.0, 2.0, 3.0], [0.0, 1.0, -1.0], [0.0, 0.0, 1.0]]
    b = [[2.0, 0.0, 1.0], [1.0, 3.0, 0.0], [0.0, 0.0, 1.0]]
    mine = written(warp.matmul, a, b)
    theirs = (np.array(a) @ np.array(b)).tolist()
    assert close(mine, theirs)


def test_1_6_compose_applies_right_to_left():
    first = written(warp.translation, 1.0, 0.0)
    second = written(warp.scaling, 2.0, 2.0)
    combined = written(warp.compose, second, first)
    # Do the translation FIRST, then the scale: (3, 5) -> (4, 5) -> (8, 10)
    assert close(combined, [[2.0, 0.0, 2.0], [0.0, 2.0, 0.0], [0.0, 0.0, 1.0]])


# -- Exercise 1.7: applying to a point ---------------------------------------


def test_1_7_apply_point_matches_the_worked_example():
    M = written(warp.translation, 3.0, -2.0)
    assert close(written(warp.apply_point, M, (8.0, 8.0)), (11.0, 6.0))


def test_1_7_a_linear_part_fixes_the_origin_and_a_translation_does_not():
    assert close(
        written(warp.apply_point, written(warp.scaling, 4.0, 0.2), (0.0, 0.0)),
        (0.0, 0.0),
    )
    assert close(
        written(warp.apply_point, written(warp.translation, 3.0, 0.0), (0.0, 0.0)),
        (3.0, 0.0),
    )


# -- Exercise 1.8: determinant ------------------------------------------------


def test_1_8_determinants_are_exact_for_whole_numbers():
    assert written(warp.determinant, written(warp.shear_x, 9.0)) == 1.0
    assert written(warp.determinant, written(warp.flip_horizontal, 9.0)) == -1.0
    assert written(warp.determinant, written(warp.scaling, 2.0, 3.0)) == 6.0
    assert written(warp.determinant, written(warp.translation, 7.0, -3.0)) == 1.0


# -- Exercise 1.9: inverse ----------------------------------------------------


def test_1_9_inverse_undoes_a_translation_and_a_shear():
    assert close(
        written(warp.invert, written(warp.translation, 3.0, -2.0)),
        written(warp.translation, -3.0, 2.0),
    )
    assert close(
        written(warp.invert, written(warp.shear_x, 2.0)),
        written(warp.shear_x, -2.0),
    )


def test_1_9_inverse_round_trips_to_the_identity():
    M = written(
        warp.compose,
        written(warp.shear_x, 0.5),
        centred(written(warp.rotation, math.radians(37))),
    )
    assert close(written(warp.compose, written(warp.invert, M), M), warp.identity())


def test_1_9_inverse_refuses_a_collapsing_transformation():
    collapse = [[1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
    try:
        warp.invert(collapse)
    except NotImplementedError as exc:
        pytest.skip(f"not written yet: {exc}")
    except warp.SingularTransform:
        return
    pytest.fail("invert should raise SingularTransform on a determinant of 0")


# -- Exercise 1.10: forward mapping, and its holes ---------------------------


def test_1_10_forward_mapping_leaves_22_holes_on_a_30_degree_rotation(img):
    matrix = centred(written(warp.rotation, math.radians(30)))
    out, holes = written(warp.warp_forward, img, matrix, fill=FILL)
    assert out.shape == (9, 9)
    assert holes.shape == (9, 9)
    assert holes.dtype == bool
    assert int(holes.sum()) == 22


def test_1_10_the_holes_are_not_all_at_the_edges(img):
    """The point of the exercise: holes appear INSIDE the glyph."""
    matrix = centred(written(warp.rotation, math.radians(30)))
    _, holes = written(warp.warp_forward, img, matrix, fill=FILL)
    interior = holes[2:7, 2:7]
    assert int(interior.sum()) > 0


def test_1_10_enlarging_forward_cannot_fill_the_output(img):
    _, holes = written(
        warp.warp_forward, img, written(warp.scaling, 2.0, 2.0),
        out_shape=(18, 18), fill=FILL,
    )
    assert int(holes.sum()) == 243


# -- Exercise 1.11: inverse mapping ------------------------------------------


def test_1_11_a_quarter_turn_is_exactly_numpy_rot90(img):
    matrix = centred(warp.rotation_quarter_turns(1))
    inverse = written(warp.invert, matrix)
    out = written(warp.warp_nearest_with_inverse, img, inverse, fill=FILL)
    assert np.array_equal(out, np.rot90(img, -1))


def test_1_11_no_pixel_is_left_unassigned(img):
    """Inverse mapping has no holes. Not fewer holes -- none."""
    matrix = centred(warp.rotation_quarter_turns(1))
    out = written(warp.warp_nearest, img, matrix, fill=FILL)
    assert int((out == FILL).sum()) == 0


def test_1_11_flips_match_numpy_exactly(img):
    assert np.array_equal(
        written(warp.warp_nearest, img, written(warp.flip_horizontal, 9.0), fill=FILL),
        np.fliplr(img),
    )
    assert np.array_equal(
        written(warp.warp_nearest, img, warp.flip_vertical(9), fill=FILL),
        np.flipud(img),
    )


def test_1_11_doubling_is_exact_pixel_replication(img):
    out = written(
        warp.warp_nearest, img, written(warp.scaling, 2.0, 2.0),
        out_shape=(18, 18), fill=FILL,
    )
    assert np.array_equal(out, np.kron(img, np.ones((2, 2), dtype=np.uint8)))


def test_1_11_halving_is_exactly_a_strided_slice(img):
    """If this fails by half a pixel, you left SAMPLE_OFFSET out."""
    out = written(
        warp.warp_nearest, img, written(warp.scaling, 0.5, 0.5),
        out_shape=(4, 4), fill=FILL,
    )
    assert np.array_equal(out, img[1::2, 1::2])


def test_1_11_translation_matches_a_slice_and_pad(img):
    out = written(warp.warp_nearest, img, written(warp.translation, 2.0, 1.0),
                  fill=FILL)
    reference = np.full_like(img, FILL)
    reference[1:, 2:] = img[:-1, :-2]
    assert np.array_equal(out, reference)


def test_1_11_a_full_turn_as_one_matrix_is_pixel_exact(img):
    matrix = centred(written(warp.rotation, 2.0 * math.pi))
    out = written(warp.warp_nearest, img, matrix, fill=FILL)
    assert np.array_equal(out, img)


def test_1_11_out_of_range_sources_take_the_fill_value_and_nothing_else(img):
    matrix = centred(written(warp.rotation, math.radians(30)))
    out = written(warp.warp_nearest, img, matrix, fill=FILL)
    back = written(warp.invert, matrix)
    for oy, ox in np.argwhere(out == FILL):
        sx, sy = warp.apply_point(back, (ox + 0.5, oy + 0.5))
        assert not (0 <= math.floor(sx) < 9 and 0 <= math.floor(sy) < 9)


# -- Exercise 1.12: Pillow's coefficients, and the comparison ----------------


def test_1_12_coefficients_are_read_off_the_inverse():
    coeffs = written(warp.to_pillow_coefficients, written(warp.translation, 1.0, 0.0))
    assert tuple(round(v, 12) for v in coeffs) == (1.0, 0.0, -1.0, 0.0, 1.0, 0.0)


def test_1_12_yours_and_pillow_agree_pixel_for_pixel(img):
    """The day's strongest claim, checked on your own code."""
    cases = [
        (1, 0, 0, 0, 1, 0),
        (1, 0, 1, 0, 1, 0),
        (1, 0, 0.5, 0, 1, 0),
        (1, 2, 0, 0, 1, 0),
        (2, 0, 0, 0, 2, 0),
        (0.5, 0, 0, 0, 0.5, 0),
        (0, -1, 9, 1, 0, 0),
        (-1, 0, 9, 0, -1, 9),
        (1, 0, -3, 0, 1, -3),
        (0.9231, 0.3129, -1.2044, -0.3129, 0.9231, 2.1177),
    ]
    for coefficients in cases:
        mine = written(
            warp.warp_nearest_with_inverse,
            img,
            warp.coefficients_to_matrix(coefficients),
            fill=FILL,
        )
        theirs = pillow_affine(img, coefficients)
        assert int((mine != theirs).sum()) == 0, f"coefficients {coefficients}"


def test_1_12_a_rotation_built_by_you_matches_pillow(img):
    matrix = centred(written(warp.rotation, math.radians(17)))
    coefficients = written(warp.to_pillow_coefficients, matrix)
    mine = written(warp.warp_nearest, img, matrix, fill=FILL)
    theirs = pillow_affine(img, coefficients)
    assert np.array_equal(mine, theirs)


# -- Exercise 2: an image is a matrix ----------------------------------------


def test_2_1_shape(img):
    assert tuple(predicted("SHAPE")) == img.shape == (9, 9)


def test_2_2_colour_shape():
    assert tuple(predicted("COLOUR_SHAPE")) == pattern.make_colour_pattern().shape


def test_2_3_the_mark_as_a_point(img):
    guess = tuple(predicted("MARK_AS_POINT"))
    assert guess == (8, 8)
    # x is the column and y is the row, so reading it back needs img[y, x].
    assert int(img[guess[1], guess[0]]) == pattern.MARK


def test_2_4_swapping_row_and_column_reads_a_different_pixel(img):
    assert predicted("VALUE_AT_SWAPPED_INDEX") == int(img[3, 4]) == pattern.PAPER
    assert int(img[4, 3]) == pattern.INK


def test_2_5_ink_pixel_count(img):
    assert predicted("INK_PIXEL_COUNT") == int((img == pattern.INK).sum()) == 24


# -- Exercise 3: forward mapping ---------------------------------------------


def test_3_1_minimum_holes_when_doubling(img):
    assert predicted("MINIMUM_HOLES_WHEN_DOUBLING") == 18 * 18 - img.size == 243


def test_3_2_holes_when_shrinking(img):
    _, holes = written(
        warp.warp_forward, img, written(warp.scaling, 0.5, 0.5),
        out_shape=(5, 5), fill=FILL,
    )
    assert predicted("HOLES_WHEN_SHRINKING") == int(holes.sum()) == 0


def test_3_3_what_shrinking_does():
    assert predicted("WHAT_SHRINKING_DOES") == "overwriting"


def test_3_4_and_3_5_fill_pixels_are_clipping_not_holes(img):
    assert predicted("FILL_PIXELS_ARE_HOLES") is False
    assert predicted("NAME_FOR_FILL_PIXELS") == "clipping"


# -- Exercise 4: inverse mapping ---------------------------------------------


def test_4_1_where_the_mark_goes(img):
    matrix = centred(warp.rotation_quarter_turns(1))
    out = written(warp.warp_nearest, img, matrix, fill=FILL)
    where = tuple(int(v) for v in np.argwhere(out == pattern.MARK)[0])
    assert tuple(predicted("MARK_AFTER_QUARTER_TURN")) == where == (8, 0)


def test_4_2_which_numpy_rotation_it_equals(img):
    k = predicted("NUMPY_ROT90_K")
    matrix = centred(warp.rotation_quarter_turns(1))
    out = written(warp.warp_nearest, img, matrix, fill=FILL)
    assert k == -1
    assert np.array_equal(out, np.rot90(img, k))


def test_4_3_no_fill_pixels_after_a_square_quarter_turn(img):
    matrix = centred(warp.rotation_quarter_turns(1))
    out = written(warp.warp_nearest, img, matrix, fill=FILL)
    assert predicted("FILL_COUNT_AFTER_QUARTER_TURN") == int((out == FILL).sum()) == 0


def test_4_4_enlarging_invents_no_new_values(img):
    out = written(
        warp.warp_nearest, img, written(warp.scaling, 2.0, 2.0),
        out_shape=(18, 18), fill=FILL,
    )
    assert predicted("DISTINCT_VALUES_AFTER_DOUBLING") == len(np.unique(out)) == 3


def test_4_5_which_slice_the_downscale_is(img):
    out = written(
        warp.warp_nearest, img, written(warp.scaling, 0.5, 0.5),
        out_shape=(4, 4), fill=FILL,
    )
    assert predicted("DOWNSCALE_IS_THE_SLICE") == "img[1::2, 1::2]"
    assert np.array_equal(out, img[1::2, 1::2])


# -- Exercise 5: homogeneous coordinates -------------------------------------


def test_5_1_why_2x2_cannot_translate():
    assert predicted("WHY_2X2_CANNOT_TRANSLATE") == "it cannot move the origin"
    for matrix in (warp.rotation_quarter_turns(1), warp.shear_y(3.0)):
        assert warp.apply_point(matrix, (0.0, 0.0)) == (0.0, 0.0)


def test_5_2_determinant_of_a_translation():
    M = written(warp.translation, 7.0, -3.0)
    assert predicted("DETERMINANT_OF_A_TRANSLATION") == written(warp.determinant, M)


def test_5_3_compose_order():
    assert predicted("COMPOSE_APPLIES_FIRST") == "A"


def test_5_4_twelve_separate_rotations_lose_pixels(img):
    assert predicted("TWELVE_SEPARATE_ROTATIONS_ARE_EXACT") is False
    step = centred(written(warp.rotation, math.radians(30)))
    out = img
    for _ in range(12):
        out = written(warp.warp_nearest, out, step, fill=FILL)
    assert int((out != img).sum()) == 16


def test_5_5_the_same_twelve_composed_are_exact(img):
    assert predicted("TWELVE_COMPOSED_ROTATIONS_ARE_EXACT") is True
    step = centred(written(warp.rotation, math.radians(30)))
    combined = warp.identity()
    for _ in range(12):
        combined = written(warp.compose, step, combined)
    assert np.array_equal(written(warp.warp_nearest, img, combined, fill=FILL), img)


# -- Exercise 6: Pillow and the half pixel -----------------------------------


def test_6_1_and_6_2_the_coefficient_direction():
    assert predicted("PILLOW_COEFFICIENT_DIRECTION") == "output to input"
    assert predicted("PICTURE_MOVES_WHEN_C_IS_POSITIVE") == "left"
    probe = np.zeros((1, 8), dtype=np.uint8)
    probe[0, 3] = 255
    assert int(np.argmax(pillow_affine(probe, (1, 0, 1, 0, 1, 0))[0])) == 2


def test_6_3_pillow_samples_at_pixel_centres():
    assert predicted("PILLOW_SAMPLES_AT") == "its centre"
    row = (np.arange(8, dtype=np.uint8) * 10).reshape(1, 8)
    observed = [int(v) for v in pillow_affine(row, (2, 0, 0, 0, 1, 0))[0]]
    centres = [
        int(row[0, math.floor(2 * (x + 0.5))])
        if math.floor(2 * (x + 0.5)) < 8 else FILL
        for x in range(8)
    ]
    assert observed == centres


def test_6_4_row_zero_moves_under_a_shear_of_two():
    """The question Day 102 deferred, answered as a number."""
    assert predicted("ROW_ZERO_SHIFT_WITH_B_EQUALS_2") == 1
    strip = np.zeros((3, 9), dtype=np.uint8)
    strip[:, 4] = 255
    out = pillow_affine(strip, (1, 2, 0, 0, 1, 0))
    assert int(np.flatnonzero(out[0] == 255)[0]) == 3  # was 4, moved by 1


def test_6_5_a_full_turn_changes_nothing(img):
    assert predicted("PIXELS_CHANGED_BY_A_FULL_TURN") == 0
    assert math.sin(2.0 * math.pi) != 0.0
    matrix = centred(written(warp.rotation, 2.0 * math.pi))
    assert int((written(warp.warp_nearest, img, matrix, fill=FILL) != img).sum()) == 0


def test_6_6_why_those_angles_disagree():
    assert (
        predicted("WHY_THOSE_ANGLES_DISAGREE")
        == "their sines and cosines land samples exactly on pixel boundaries"
    )

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

Twelve functions to write. Each 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 (the one above this file):

    .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 for the MATRIX arithmetic. NumPy is
allowed for holding pixels -- that is what it is for -- and it appears in the
tests, where it checks your work. Writing the maths yourself and then having
Pillow agree with you pixel for pixel is the whole point of the day, and it
only means something if you did not build your answer out of someone else's.

Two conventions, fixed and used everywhere:

1. A point is (x, y): x is the COLUMN, y is the ROW, and y grows DOWNWARD.
   An image array is indexed `img[y, x]`.

2. A transformation matrix is 3 by 3 and maps INPUT to OUTPUT -- the direction
   you can see:

       [ a  b  tx ] [ x ]   [ a*x + b*y + tx ]
       [ c  d  ty ] [ y ] = [ c*x + d*y + ty ]
       [ 0  0   1 ] [ 1 ]   [        1       ]

   The top-left 2 by 2 block is Day 102's linear part; its columns are still
   where the basis vectors land. The third column is the translation.
"""

from __future__ import annotations

import math

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

# Written for you, and it matters more than it looks. Every output pixel is a
# little square; the point the transformation is evaluated at is its CENTRE,
# not its corner. Pixel (x, y) covers (x, y) to (x + 1, y + 1), so its centre
# is (x + 0.5, y + 0.5).
#
# Exercise 6 measures why this half is not optional.
SAMPLE_OFFSET = 0.5

TOL = 1e-12


class SingularTransform(ValueError):
    """Raised when a transformation has no inverse, so it cannot be applied.

    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 translation(tx: float, ty: float) -> Matrix:
    """Move every point by (tx, ty).

    This is the function homogeneous coordinates exist for. Day 102 proved a
    linear map cannot move the origin, and this moves the origin, so no 2 by 2
    matrix can do it. In 3 by 3 it is easy: the constant gets multiplied by the
    third coordinate, which is always 1.

    Worked example. translation(3, -2) must be

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

    and it must send (0, 0) to (3, -2) and (8, 8) to (11, 6).

    Return floats, not ints, so that later arithmetic does not surprise you.
    """
    raise NotImplementedError("Exercise 1.1: build the translation matrix")


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


def scaling(sx: float, sy: float) -> Matrix:
    """Stretch x by sx and y by sy, about the origin -- the top-left corner.

    Same derivation as Day 102: where does one step right go? To (sx, 0). Where
    does one step down go? To (0, sy). Those are the two columns. The
    translation column is zero, because scaling about the origin leaves the
    origin alone.

    Worked example. scaling(2, 3) sends (1, 1) to (2, 3) and (0, 0) to (0, 0).
    """
    raise NotImplementedError("Exercise 1.2: build the scaling matrix")


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


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

    The linear part is exactly Day 102's:

        [[cos, -sin],
         [sin,  cos]]

    Put it in the top-left of a 3 by 3 with a zero translation column.

    One thing to expect and NOT to treat as a bug: on an image, y grows
    DOWNWARD, so this turns the picture CLOCKWISE on screen even though it is
    the counter-clockwise matrix from Day 102. Nothing about the matrix
    changed; the picture is flipped relative to the graph paper.

    Worked example. rotation(math.pi / 2) sends (1, 0) to approximately
    (0, 1) -- approximately, because math.cos(math.pi / 2) is
    6.123233995736766e-17 and not 0.0. That is why every test here states a
    tolerance.
    """
    raise NotImplementedError("Exercise 1.3: build the rotation matrix")


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


def shear_x(k: float) -> Matrix:
    """Slide each row sideways in proportion to its y: x becomes x + k*y.

    Derive it the Day 102 way. Where does one step right, (1, 0), go? Its y is
    0, so it does not move: to (1, 0). Where does one step down, (0, 1), go?
    Its y is 1, so it slides k to the right: to (k, 1). Those two landings are
    the two columns.

    Worked example. shear_x(2) must be

        [[1, 2, 0],
         [0, 1, 0],
         [0, 0, 1]]

    and it must leave (5, 0) exactly where it was.
    """
    raise NotImplementedError("Exercise 1.4: build the shear matrix")


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


def flip_horizontal(width: float) -> Matrix:
    """Mirror left-to-right inside an image `width` pixels wide.

    The trap: a bare reflection sends x to -x, which puts the whole picture off
    the left-hand edge. What you want is a mirror about the image's own centre
    line, which sends x to width - x. That is a reflection FOLLOWED BY a
    translation of `width` -- and because you now have homogeneous
    coordinates, it is ONE matrix rather than two steps.

    Worked example. flip_horizontal(9) sends (0, 4) to (9, 4) and (9, 4) to
    (0, 4), and leaves every y alone.

    Hint: the linear part is [[-1, 0], [0, 1]] and the translation column is
    (width, 0).
    """
    raise NotImplementedError("Exercise 1.5: build the horizontal-flip matrix")


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


def matmul(a: Matrix, b: Matrix) -> Matrix:
    """The 3 by 3 matrix product a @ b, computed by hand.

    Entry (i, j) of the result is row i of `a` dotted with column j of `b`:

        result[i][j] = sum over k of a[i][k] * b[k][j]

    Three nested loops, or one comprehension. Do not import NumPy for this --
    checking your matrix product against NumPy's is one of the tests, and it
    proves nothing if you used NumPy to compute it.
    """
    raise NotImplementedError("Exercise 1.6: multiply two 3 by 3 matrices")


def compose(*matrices: Matrix) -> Matrix:
    """Combine transformations into ONE matrix, applied RIGHT to LEFT.

    `compose(B, A)` means "do A first, then B" -- the Day 101 convention, and
    the same order as reading B(A(x)) from the inside out.

    Written for you, on top of your `matmul`, because the order convention is
    the part worth getting right and the loop is not. Read it.
    """
    if not matrices:
        return identity()
    result = matrices[0]
    for m in matrices[1:]:
        result = matmul(result, m)
    return result


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


def apply_point(matrix: Matrix, point: Point) -> Point:
    """Send one (x, y) point through the matrix and return the new (x, y).

    The third coordinate is always 1 going in, and for an affine matrix it is
    always 1 coming out, so you never have to build the triple explicitly:

        new_x = matrix[0][0]*x + matrix[0][1]*y + matrix[0][2]
        new_y = matrix[1][0]*x + matrix[1][1]*y + matrix[1][2]

    Worked example. With translation(3, -2), the point (8, 8) becomes
    (11.0, 6.0).
    """
    raise NotImplementedError("Exercise 1.7: apply a matrix to a point")


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


def determinant(matrix: Matrix) -> float:
    """The determinant of the LINEAR part -- the area factor, as on Day 102.

    The third row of an affine matrix is (0, 0, 1), so the full 3 by 3
    determinant equals the 2 by 2 determinant of the top-left block:

        a*d - b*c   where the block is [[a, b], [c, d]]

    Compute it directly rather than via NumPy, so that whole numbers stay
    exact. Day 102 measured this: numpy.linalg.det can return 7.000000000000001
    where the direct formula returns exactly 7.

    Worked example. determinant(shear_x(9)) is exactly 1.0.
    determinant(flip_horizontal(9)) is exactly -1.0.
    """
    raise NotImplementedError("Exercise 1.8: compute the determinant")


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


def invert(matrix: Matrix) -> Matrix:
    """Invert an affine 3 by 3 matrix, using its structure rather than brute force.

    Write the matrix as a linear part A and a translation column t:

        M = [ A  t ]      M^-1 = [ A^-1   -A^-1 t ]
            [ 0  1 ]             [  0         1   ]

    So: invert the 2 by 2 block the Day 102 way (swap a and d, negate b and c,
    divide everything by the determinant), then apply that inverted block to
    the translation column and negate the result.

        A^-1 = (1/det) * [[ d, -b],
                          [-c,  a]]

    Raise `SingularTransform` when abs(determinant) <= TOL. A transformation
    that flattens the picture onto a line has thrown information away, and no
    arithmetic puts it back.

    One cosmetic detail worth copying: add 0.0 to each entry before returning
    it. Negating a 0.0 -- which every translation and every axis-aligned scale
    has -- produces -0.0, which compares equal to 0.0 but PRINTS as "-0.0" and
    turns up as noise in every coefficient tuple. Adding 0.0 normalises it.

    Worked example. invert(translation(3, -2)) is translation(-3, 2).
    invert(shear_x(2)) is shear_x(-2).
    """
    raise NotImplementedError("Exercise 1.9: invert the matrix")


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


def warp_forward(image, matrix: Matrix, out_shape=None, fill=0):
    """Forward mapping: push every INPUT pixel to where it lands. This is WRONG.

    You are writing it anyway, because seeing it fail is the argument for the
    method in 1.11. Do not fix it. The holes are the result.

    The algorithm:

        make an output array of `out_shape`, filled with `fill`
        make a boolean array the same shape, all False, called `written`
        for each input pixel (y, x):
            send its CENTRE, (x + SAMPLE_OFFSET, y + SAMPLE_OFFSET), through
                the matrix with apply_point
            floor both results to get the output pixel (ox, oy)
            if that pixel is inside the output:
                out[oy, ox] = image[y, x]
                written[oy, ox] = True

    Return the pair `(out, holes)` where `holes` is the boolean array that is
    True wherever nothing was ever written -- that is, `~written`.

    NumPy is fine here: `numpy.full(shape, fill, dtype=src.dtype)` and
    `numpy.zeros(shape, dtype=bool)`. It is holding pixels, not doing your
    maths.

    Worked example. A 30 degree rotation of the 9 by 9 test pattern about its
    centre leaves 22 of the 81 output pixels unwritten. If you get 22, you have
    it right -- including the ones punched through the middle of the glyph.
    """
    raise NotImplementedError("Exercise 1.10: forward mapping, holes and all")


# -- Exercise 1.11 ------------------------------------------------------------


def warp_nearest_with_inverse(image, inverse: Matrix, out_shape=None, fill=0):
    """Inverse mapping with nearest-neighbour sampling. This is the RIGHT way.

    Turn the loop inside out. Walk the OUTPUT, not the input:

        make an output array of `out_shape`, filled with `fill`
        for each output pixel (oy, ox):
            take its CENTRE, (ox + SAMPLE_OFFSET, oy + SAMPLE_OFFSET)
            send it through `inverse` with apply_point
            floor both results to get the input pixel (ix, iy)
            if that pixel is inside the input:
                out[oy, ox] = image[iy, ix]
            otherwise leave the fill value -- the source is off the picture,
                which is clipping, and clipping is a decision, not an error

    Every output pixel is visited exactly once, so holes are impossible. Not
    because you were careful: because of which array the loop is over.

    Note the argument. This takes the OUTPUT-to-INPUT matrix, already inverted,
    because that is the form Pillow's coefficients come in and exercise 6
    hands both implementations the identical six numbers.

    Two details that decide whether you match Pillow exactly:
      * SAMPLE_OFFSET, the half. Leave it out and everything is half a pixel
        adrift, which looks like a mysterious blur rather than like an offset.
      * `math.floor`, not `round` and not `int`. You want the input pixel whose
        SQUARE contains the point. `int` truncates toward zero, which is wrong
        for negatives; `round` is a different rule entirely.

    Worked example. A quarter turn of the test pattern about its centre must
    equal numpy.rot90(img, -1) EXACTLY -- every one of the 81 pixels.
    """
    raise NotImplementedError("Exercise 1.11: inverse mapping, nearest neighbour")


def warp_nearest(image, matrix: Matrix, out_shape=None, fill=0):
    """`warp_nearest_with_inverse`, but taking the transformation you can SEE.

    Written for you, on top of your 1.9 and 1.11. It inverts here, once, rather
    than making every caller remember to.
    """
    return warp_nearest_with_inverse(
        image, invert(matrix), out_shape=out_shape, fill=fill
    )


# -- Exercise 1.12 ------------------------------------------------------------


def to_pillow_coefficients(matrix: Matrix):
    """Turn a visible input-to-output matrix into Pillow's six coefficients.

    Pillow's `Image.transform(size, Image.Transform.AFFINE, coeffs)` takes
    `(a, b, c, d, e, f)` meaning

        input_x = a * output_x + b * output_y + c
        input_y = d * output_x + e * output_y + f

    -- the OUTPUT-to-INPUT direction, which is the INVERSE of the effect you
    see. Day 102 confirmed that direction by experiment and exercise 6 confirms
    it again. Forgetting to invert is the single most common way to get a
    Pillow transform backwards, and it does not raise: it just moves the
    picture the wrong way.

    So: invert the matrix, then read the six numbers off the first two rows,
    left to right, top row first.

    Worked example. to_pillow_coefficients(translation(1, 0)) must be
    (1.0, 0.0, -1.0, 0.0, 1.0, 0.0). Note the minus.
    """
    raise NotImplementedError("Exercise 1.12: read off Pillow's coefficients")


# =============================================================================
# Written for you. Read these -- the tests use them, and two of them are the
# answers to questions the exercises above ask you to think about.
# =============================================================================


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


def rotation_quarter_turns(turns: int) -> Matrix:
    """Rotate by an exact multiple of 90 degrees, with integer entries.

    `rotation(math.pi / 2)` is correct but its cosine is 6.123233995736766e-17
    rather than 0.0. Where the answer is meant to be checkable to the exact
    pixel, build the matrix from integers and the float noise never enters.
    """
    cos_t, sin_t = [(1, 0), (0, 1), (-1, 0), (0, -1)][turns % 4]
    return [
        [float(cos_t), float(-sin_t), 0.0],
        [float(sin_t), float(cos_t), 0.0],
        [0.0, 0.0, 1.0],
    ]


def shear_y(k: float) -> Matrix:
    """The other shear: y becomes y + k*x."""
    return [[1.0, 0.0, 0.0], [float(k), 1.0, 0.0], [0.0, 0.0, 1.0]]


def flip_vertical(height: float) -> Matrix:
    """Mirror top-to-bottom inside an image `height` pixels tall."""
    return [[1.0, 0.0, 0.0], [0.0, -1.0, float(height)], [0.0, 0.0, 1.0]]


def about_centre(matrix: Matrix, width: float, height: float) -> Matrix:
    """Do `matrix` about the image's centre instead of its top-left corner.

    Move the centre to the origin, transform, move it back: T(+c) . M . T(-c).
    Three matrices folded into one, which is only possible because translation
    became a matrix. Uses your `compose` and your `translation`.
    """
    cx, cy = width / 2.0, height / 2.0
    return compose(translation(cx, cy), matrix, translation(-cx, -cy))


def matrices_close(a: Matrix, b: Matrix, tol: float = TOL) -> bool:
    """True when two 3 by 3 matrices agree entry by entry within `tol`."""
    return all(abs(a[i][j] - b[i][j]) <= tol for i in range(3) for j in range(3))


def coefficients_to_matrix(coeffs) -> Matrix:
    """Six Pillow coefficients back into a 3 by 3 matrix."""
    a, b, c, d, e, f = coeffs
    return [[a, b, c], [d, e, f], [0.0, 0.0, 1.0]]

#!/usr/bin/env bash
# Tests for the Day 105 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:
#
#   * an image is a matrix, its shape is (height, width), and img[4, 3] is a
#     different pixel from img[3, 4] -- the ordering trap, measured;
#   * forward mapping leaves 22 of 81 output pixels unwritten on a 30 degree
#     rotation and 243 of 324 when doubling, and inverse mapping leaves none;
#   * a quarter turn is EXACTLY numpy.rot90(img, -1), a flip is exactly
#     numpy.fliplr, doubling is exactly numpy.kron, and halving is exactly the
#     strided slice img[1::2, 1::2];
#   * translation is not linear, needs a third coordinate, and then composes
#     with everything else -- checked as a matrix and as pixels;
#   * Pillow's affine coefficients run OUTPUT to INPUT, and Pillow samples at
#     each output pixel's CENTRE -- both settled by measurement, which is the
#     question Day 102 deferred to today;
#   * a shear coefficient of 2.0 moves row 0 by one whole pixel, and the half-
#     pixel offset explains exactly why;
#   * this implementation and Pillow produce byte-for-byte identical output on
#     510 affine transformations, and differ on 8 of the 360 whole-degree
#     rotations by at most 2 pixels -- every disagreement a floating-point tie
#     at a pixel boundary, which is asserted rather than hidden;
#   * a 360 degree rotation as ONE matrix is pixel-exact; as twelve separate
#     resampling passes it loses 16 of 81 pixels;
#   * nothing is downloaded, no image is written into the lab, and nothing is
#     left behind on disk.
#
# Everything runs offline. Nothing binds a port, nothing writes outside the
# lab or a temporary directory, nothing needs a key. Deterministic,
# non-interactive, exits 0 only if every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1

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

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

failures=0
checks=0

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

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

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

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

# The Python that owns that pytest is the one with numpy and Pillow 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

for module in numpy PIL; do
  if ! "${python_bin}" -c "import ${module}" >/dev/null 2>&1; then
    echo "FAIL: ${module} 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
done

echo "Day 105 — Rotate It Yourself"
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", "pillow", "pytest"):
    print(f"{name:<8} {version(name)}")
print(f"platform {platform.platform()}")
print(f"exe      {sys.executable.rsplit('/', 3)[-1]}")
PY
)"
echo "${versions}" | sed 's/^/  /'

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

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

pil_major="$("${python_bin}" -c "import PIL; print(PIL.__version__.split('.')[0])")"
check_eq "Pillow is version 12 or later" "12" "${pil_major}"

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

for script in 01_an_image_is_a_matrix 02_forward_mapping_leaves_holes \
              03_inverse_mapping 04_scale_shear_flip \
              05_homogeneous_and_composition 06_against_pillow; 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 pixels, real values"
# --------------------------------------------------------------------------

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 60 ]; then
  check "the reference suite ran at least 60 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 60 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 `warp` and
# `pattern`, 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 os
import random
import tempfile

import numpy as np
from PIL import Image

import pattern
import warp

FILL = pattern.FILL
img = pattern.make_pattern()
H, W = img.shape


def pil(array, coefficients, out_shape=None, resample=None):
    h, w = array.shape
    oh, ow = out_shape or (h, w)
    return np.asarray(
        Image.fromarray(array, mode="L").transform(
            (ow, oh), Image.Transform.AFFINE, coefficients,
            resample=resample or Image.Resampling.NEAREST, fillcolor=FILL,
        )
    )


# -- the image itself
print("shape", img.shape)
print("colour_shape", pattern.make_colour_pattern().shape)
print("dtype", img.dtype)
print("ink_count", int((img == pattern.INK).sum()))
print("row_col_swap_differs", int(img[4, 3]), int(img[3, 4]))
print("asymmetric", (not np.array_equal(img, np.fliplr(img)),
                     not np.array_equal(img, np.flipud(img)),
                     not np.array_equal(img, img.T)))
print("generated_twice_identical",
      np.array_equal(pattern.make_pattern(), pattern.make_pattern()))

# -- forward mapping
rot30 = warp.about_centre(warp.rotation(math.radians(30)), W, H)
_, holes = warp.warp_forward(img, rot30, fill=FILL)
print("forward_holes_rot30", int(holes.sum()))
print("forward_holes_inside_glyph", int(holes[2:7, 2:7].sum()) > 0)
_, big_holes = warp.warp_forward(img, warp.scaling(2.0, 2.0),
                                 out_shape=(18, 18), fill=FILL)
print("forward_holes_double", int(big_holes.sum()), 18 * 18 - img.size)

inv30 = warp.warp_nearest(img, rot30, fill=FILL)
back = warp.invert(rot30)
filled = np.argwhere(inv30 == FILL)
outside = sum(
    1 for oy, ox in filled
    if not (0 <= math.floor(warp.apply_point(back, (ox + .5, oy + .5))[0]) < W
            and 0 <= math.floor(warp.apply_point(back, (ox + .5, oy + .5))[1]) < H)
)
print("inverse_fill_count", len(filled))
print("inverse_fill_all_outside", len(filled) == outside and len(filled) > 0)

# -- exact answers
q = warp.about_centre(warp.rotation_quarter_turns(1), W, H)
turned = warp.warp_nearest(img, q, fill=FILL)
print("quarter_turn_is_rot90", np.array_equal(turned, np.rot90(img, -1)))
print("quarter_turn_no_fill", int((turned == FILL).sum()))
print("mark_before", tuple(int(v) for v in np.argwhere(img == pattern.MARK)[0]))
print("mark_after", tuple(int(v) for v in np.argwhere(turned == pattern.MARK)[0]))
print("flip_h_is_fliplr", np.array_equal(
    warp.warp_nearest(img, warp.flip_horizontal(W), fill=FILL), np.fliplr(img)))
print("flip_v_is_flipud", np.array_equal(
    warp.warp_nearest(img, warp.flip_vertical(H), fill=FILL), np.flipud(img)))
print("double_is_kron", np.array_equal(
    warp.warp_nearest(img, warp.scaling(2.0, 2.0), out_shape=(18, 18), fill=FILL),
    np.kron(img, np.ones((2, 2), dtype=np.uint8))))
print("halve_is_strided_slice", np.array_equal(
    warp.warp_nearest(img, warp.scaling(0.5, 0.5), out_shape=(4, 4), fill=FILL),
    img[1::2, 1::2]))

# -- translation and homogeneous coordinates
print("translation_moves_origin",
      warp.apply_point(warp.translation(3.0, -2.0), (0.0, 0.0)))
print("linear_fixes_origin", all(
    warp.apply_point(m, (0.0, 0.0)) == (0.0, 0.0)
    for m in (warp.rotation(1.1), warp.scaling(4.0, .2), warp.shear_x(9.0))))
print("translation_det", warp.determinant(warp.translation(7.0, -3.0)))
print("bottom_row_always", all(
    m[2] == [0.0, 0.0, 1.0] for m in (
        warp.translation(3, -2), warp.scaling(2, .5), warp.rotation(1.1),
        warp.shear_x(2), warp.flip_horizontal(9), q)))
det_parts = [warp.scaling(1.5, 1.5), warp.shear_x(0.5), warp.rotation(0.7)]
prod = 1.0
for m in det_parts:
    prod *= warp.determinant(m)
print("det_of_composition_is_product",
      abs(warp.determinant(warp.compose(*det_parts)) - prod) <= warp.TOL)

# -- resampling once versus repeatedly
full = warp.about_centre(warp.rotation(2.0 * math.pi), W, H)
print("full_turn_one_matrix_diff",
      int((warp.warp_nearest(img, full, fill=FILL) != img).sum()))
cur = img
for _ in range(12):
    cur = warp.warp_nearest(cur, rot30, fill=FILL)
print("full_turn_twelve_passes_diff", int((cur != img).sum()))
comb = warp.identity()
for _ in range(12):
    comb = warp.compose(rot30, comb)
print("full_turn_twelve_composed_diff",
      int((warp.warp_nearest(img, comb, fill=FILL) != img).sum()))

# -- the Pillow conventions, measured
probe = np.zeros((1, 8), dtype=np.uint8)
probe[0, 3] = 255
print("pillow_positive_c_moves_left",
      int(np.argmax(probe[0])), int(np.argmax(pil(probe, (1, 0, 1, 0, 1, 0))[0])))
print("to_pillow_coefficients_inverts",
      tuple(round(v, 12) for v in warp.to_pillow_coefficients(warp.translation(1.0, 0.0))))

row = (np.arange(8, dtype=np.uint8) * 10).reshape(1, 8)
obs = [int(v) for v in pil(row, (2, 0, 0, 0, 1, 0))[0]]
centres = [int(row[0, math.floor(2 * (x + .5))])
           if math.floor(2 * (x + .5)) < 8 else FILL for x in range(8)]
corners = [int(row[0, math.floor(2 * x + .5)])
           if math.floor(2 * x + .5) < 8 else FILL for x in range(8)]
print("pillow_matches_centres", obs == centres)
print("pillow_matches_corners", obs == corners)
print("sample_offset", warp.SAMPLE_OFFSET)

strip = np.zeros((3, 9), dtype=np.uint8)
strip[:, 4] = 255
sh = pil(strip, (1, 2, 0, 0, 1, 0))
print("shear_row0_line_at", int(np.flatnonzero(sh[0] == 255)[0]))
print("shear_row1_line_at", int(np.flatnonzero(sh[1] == 255)[0]))
print("shear_row0_shift_predicted", math.floor(0.5 + 2 * 0.5))

# -- ours against theirs
rng = random.Random(105)
cases = []
for _ in range(500):
    th = rng.uniform(-math.pi, math.pi)
    sc = rng.uniform(0.4, 2.5)
    sk = rng.uniform(-2.5, 2.5)
    ca, sa = math.cos(th), math.sin(th)
    cases.append((sc * ca, sc * (ca * sk - sa), rng.uniform(-6, 6),
                  sc * sa, sc * (sa * sk + ca), rng.uniform(-6, 6)))
cases += [(1, 0, 0, 0, 1, 0), (1, 0, 1, 0, 1, 0), (1, 0, .5, 0, 1, 0),
          (1, 2, 0, 0, 1, 0), (2, 0, 0, 0, 2, 0), (.5, 0, 0, 0, .5, 0),
          (0, -1, 9, 1, 0, 0), (-1, 0, 9, 0, -1, 9), (1, 0, -3, 0, 1, -3),
          (1, .5, 0, 0, 1, 0)]
worst = 0
for co in cases:
    mine = warp.warp_nearest_with_inverse(
        img, warp.coefficients_to_matrix(co), fill=FILL)
    worst = max(worst, int((mine != pil(img, co)).sum()))
print("random_cases", len(cases))
print("random_worst_differing_pixels", worst)

bad_angles, worst_rot, furthest = [], 0, 0.0
for deg in range(360):
    M = warp.about_centre(warp.rotation(math.radians(deg)), W, H)
    co = warp.to_pillow_coefficients(M)
    wrong = np.argwhere(warp.warp_nearest(img, M, fill=FILL) != pil(img, co))
    if len(wrong):
        bad_angles.append(deg)
        worst_rot = max(worst_rot, len(wrong))
    a, b, c, d, e, f = co
    for oy, ox in wrong:
        xs = a * (ox + .5) + b * (oy + .5) + c
        ys = d * (ox + .5) + e * (oy + .5) + f
        furthest = max(furthest, min(abs(xs - round(xs)), abs(ys - round(ys))))
print("rotation_sweep_disagreeing", bad_angles)
print("rotation_sweep_worst_pixels", worst_rot)
print("rotation_sweep_all_at_boundaries", furthest < 1e-9)

M30 = warp.about_centre(warp.rotation(math.radians(30)), W, H)
co30 = warp.to_pillow_coefficients(M30)
w30 = np.argwhere(warp.warp_nearest(img, M30, fill=FILL) != pil(img, co30))
oy, ox = int(w30[0][0]), int(w30[0][1])
ys30 = co30[3] * (ox + .5) + co30[4] * (oy + .5) + co30[5]
print("thirty_degree_disagreements", len(w30))
print("thirty_degree_pixel", (oy, ox))
print("thirty_degree_source_y_is_just_under_5",
      abs(ys30 - 5.0) < 1e-14 and ys30 != 5.0)

# -- bilinear
worst_inside, worst_any = 0.0, 0.0
for M in (warp.translation(.25, .25),
          warp.about_centre(warp.rotation(math.radians(30)), W, H),
          warp.about_centre(warp.rotation(math.radians(17)), W, H),
          warp.about_centre(warp.scaling(1.5, 1.5), W, H),
          warp.shear_x(0.4)):
    inv = warp.invert(M)
    mine = warp.warp_bilinear_with_inverse(img, inv, fill=0.0)
    theirs = pil(img, warp.to_pillow_coefficients(M),
                 resample=Image.Resampling.BILINEAR).astype(float)
    d = np.abs(mine - theirs)
    inside = np.zeros((H, W), dtype=bool)
    for y in range(H):
        for x in range(W):
            sx, sy = warp.apply_point(inv, (x + .5, y + .5))
            x0, y0 = math.floor(sx - .5), math.floor(sy - .5)
            inside[y, x] = 0 <= x0 and x0 + 1 < W and 0 <= y0 and y0 + 1 < H
    worst_inside = max(worst_inside, float(d[inside].max()))
    worst_any = max(worst_any, float(d.max()))
print("bilinear_worst_all_inside", round(worst_inside, 6))
print("bilinear_worst_anywhere_exceeds_100", worst_any > 100.0)

# -- files
with tempfile.TemporaryDirectory() as tmp:
    p = os.path.join(tmp, "pattern.png")
    Image.fromarray(img, mode="L").save(p)
    reloaded = np.asarray(Image.open(p).convert("L"))
    print("png_round_trip_lossless", np.array_equal(reloaded, img))
print("temp_file_removed", not os.path.exists(p))
PY
)"

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

check_eq "the image is a (height, width) array of 9 by 9" "(9, 9)" "$(get shape)"
check_eq "the colour image is three stacked planes" "(9, 9, 3)" "$(get colour_shape)"
check_eq "one byte per pixel" "uint8" "$(get dtype)"
check_eq "the pattern has 24 ink pixels" "24" "$(get ink_count)"
check_eq "img[4, 3] and img[3, 4] are DIFFERENT pixels" \
  "255 0" "$(get row_col_swap_differs)"
check_eq "the pattern is asymmetric under mirror, flip and transpose" \
  "(True, True, True)" "$(get asymmetric)"
check_eq "the pattern is generated, not loaded: two calls agree" \
  "True" "$(get generated_twice_identical)"

check_eq "forward mapping leaves 22 holes on a 30 degree rotation" \
  "22" "$(get forward_holes_rot30)"
check_eq "and some of those holes are INSIDE the glyph, not at the edge" \
  "True" "$(get forward_holes_inside_glyph)"
check_eq "forward mapping cannot fill a doubled output: 243 holes of 324" \
  "243 243" "$(get forward_holes_double)"
check_eq "inverse mapping leaves 12 pixels at the fill value" \
  "12" "$(get inverse_fill_count)"
check_eq "and every one of them is clipping, not a hole" \
  "True" "$(get inverse_fill_all_outside)"

check_eq "a quarter turn is EXACTLY numpy.rot90(img, -1)" \
  "True" "$(get quarter_turn_is_rot90)"
check_eq "a quarter turn of a square image clips nothing" \
  "0" "$(get quarter_turn_no_fill)"
check_eq "the corner mark starts at row 8, column 8" "(8, 8)" "$(get mark_before)"
check_eq "and a clockwise quarter turn puts it at row 8, column 0" \
  "(8, 0)" "$(get mark_after)"
check_eq "a horizontal flip is EXACTLY numpy.fliplr" "True" "$(get flip_h_is_fliplr)"
check_eq "a vertical flip is EXACTLY numpy.flipud" "True" "$(get flip_v_is_flipud)"
check_eq "doubling is EXACTLY numpy.kron with a 2 by 2 block of ones" \
  "True" "$(get double_is_kron)"
check_eq "halving is EXACTLY the strided slice img[1::2, 1::2]" \
  "True" "$(get halve_is_strided_slice)"

check_eq "translation moves the origin, so it is not linear" \
  "(3.0, -2.0)" "$(get translation_moves_origin)"
check_eq "every purely linear part leaves the origin alone" \
  "True" "$(get linear_fixes_origin)"
check_eq "a translation's determinant is exactly 1" \
  "1.0" "$(get translation_det)"
check_eq "every affine matrix has the bottom row (0, 0, 1)" \
  "True" "$(get bottom_row_always)"
check_eq "the determinant of a composition is the product of the parts" \
  "True" "$(get det_of_composition_is_product)"

check_eq "a full turn as ONE matrix changes no pixel" \
  "0" "$(get full_turn_one_matrix_diff)"
check_eq "the same full turn as twelve passes loses 16 pixels" \
  "16" "$(get full_turn_twelve_passes_diff)"
check_eq "and composing those twelve into one matrix is exact again" \
  "0" "$(get full_turn_twelve_composed_diff)"

check_eq "a POSITIVE Pillow c coefficient moves the picture LEFT" \
  "3 2" "$(get pillow_positive_c_moves_left)"
check_eq "to_pillow_coefficients reads off the INVERSE, so c is negative" \
  "(1.0, 0.0, -1.0, 0.0, 1.0, 0.0)" "$(get to_pillow_coefficients_inverts)"

# Section 6 re-runs this script with D105_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_centres="True"
if [ -n "${D105_SELF_TEST:-}" ]; then
  expected_centres="False"   # the naive belief, deliberately wrong here
fi
check_eq "Pillow samples at each output pixel's CENTRE" \
  "${expected_centres}" "$(get pillow_matches_centres)"
check_eq "Pillow does NOT sample at integer corners" \
  "False" "$(get pillow_matches_corners)"
check_eq "this lab uses the same half-pixel offset" "0.5" "$(get sample_offset)"

check_eq "a shear coefficient of 2 moves row 0 by one whole pixel" \
  "3" "$(get shear_row0_line_at)"
check_eq "and moves row 1 by three" "1" "$(get shear_row1_line_at)"
check_eq "the half-pixel offset predicts that row 0 shift exactly" \
  "1" "$(get shear_row0_shift_predicted)"

check_eq "510 affine transformations were compared with Pillow" \
  "510" "$(get random_cases)"
check_eq "and every one of them agreed byte for byte" \
  "0" "$(get random_worst_differing_pixels)"
check_eq "8 of the 360 whole-degree rotations DO disagree, and they are named" \
  "[30, 60, 120, 150, 210, 240, 300, 330]" "$(get rotation_sweep_disagreeing)"
check_eq "no disagreement is larger than 2 pixels of 81" \
  "2" "$(get rotation_sweep_worst_pixels)"
check_eq "and every disagreeing sample sits on a pixel boundary" \
  "True" "$(get rotation_sweep_all_at_boundaries)"
check_eq "the 30 degree case differs in exactly one pixel" \
  "1" "$(get thirty_degree_disagreements)"
check_eq "that pixel is row 4, column 3" "(4, 3)" "$(get thirty_degree_pixel)"
check_eq "its source row is 4.999999999999999 rather than 5" \
  "True" "$(get thirty_degree_source_y_is_just_under_5)"

check_eq "bilinear agrees with Pillow within 1 grey level away from the border" \
  "1.0" "$(get bilinear_worst_all_inside)"
check_eq "and diverges by more than 100 levels AT the border, which is stated" \
  "True" "$(get bilinear_worst_anywhere_exceeds_100)"

check_eq "a PNG round trip is lossless" "True" "$(get png_round_trip_lossless)"
check_eq "and the temporary file is gone afterwards" "True" "$(get temp_file_removed)"

# --------------------------------------------------------------------------
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 Pillow samples at integer pixel corners, and
# asserts that the re-run reports the failure and exits non-zero. If this
# section passes, section 5 is not decorative.
if [ -z "${D105_SELF_TEST:-}" ]; then
  self_out="$(D105_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: Pillow samples at each output pixel's CENTRE"*)
      check "the failing check is named in the output with both values" "yes" ;;
    *) check "the failing check is named in the output with both values" "no" ;;
  esac
  case "${self_out}" in
    *", 1 failure(s)."*)
      check "the summary line counts exactly one failure" "yes" ;;
    *) check "the summary line counts exactly one failure" "no" ;;
  esac
else
  echo "  (self-test run: section 6 does not recurse)"
fi

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

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

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

# The test image is GENERATED, not downloaded and not committed. If an image
# file ever appears in the lab's own tree, either something was committed by
# mistake or a script wrote one and failed to clean up. Pillow ships a pile of
# its own test images inside site-packages, so .venv is pruned here too.
image_files="$(find "${lab_dir}" -name '.venv' -prune -o -type f \
  \( -name '*.png' -o -name '*.jpg' -o -name '*.jpeg' -o -name '*.bmp' \
     -o -name '*.gif' -o -name '*.tif' -o -name '*.tiff' \) -print 2>/dev/null \
  | wc -l | tr -d ' ')"
check_eq "no image file in the lab's own tree: the pattern is generated" \
  "0" "${image_files}"

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

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